The level set equation (2.15) is easy to handle and use. The tricky part of the level set modeling is the reinitialization process. There exists other approaches for the reinitialization:
• It is possible to apply ENO1-schemes in order to minimize effects of oscilla-tions and Runge phenomena at discontinuities. Examples are [63, 66]. ENO and WENO schemes are fairly easily applied on structured grids. On general unstructured grids the procedure becomes much more complicated. The use of WENO methods has been proposed [77], but remains a challenge.
• Explicit calculation of the distance to the interface. An example is [67], where an explicit 1D surface function is created, on the formy=h(x) or x=h(y) depending on the rotation of the surface. Then the distance from a point in space (x, y) to the surface function is calculated, using a Newton iteration proces.
A note about the contour plots produced in this thesis: We do not have a contour plotter which can produce high sub-element accuracy of the contours. The contour plotter is based on a triangulation of the nodes within each element. Within each element, linear interpolation between the nodes is used, giving a number of straight lines as contours, as in Figure 5.28, where the +-marker marks the line end points. Hence, it is onlyO((∆x/n)2) order accurate,nbeing the degree of the
−0.1 0 0.1
−0.1 0 0.1
Figure 5.28: Contour plot accuracy polynomial spacePn.
1Essentially Non-Oscillatory
C H A P T E R 6 Numerical Tests
I can show you that the art of calculation has to do with odd and even numbers in their numerical relations to themselves and to each other.
— Plato 427–347 BC In the following we shall present a few examples validating the performance of the proposed schemes.
The examples are intended for validation of the schemes, not for measuring of computational performance and accuracy. Thus we have not made any attempt to optimize the meshes and elements, which in all tests are of same size in the entire domain. However, in all tests there are regions which could be resolved using much larger elements, effectively decreasing the total number of elements and the computing time, and still maintaining accuracy.
All test have been carried out in Matlab, either on a SUN ultra-SPARCIII 1000Mhz system or a Intel PentiumIIIm 1200Mhz laptop. During implementation and testing, focus has been on flexibility to add new features and components, and not computational aspects as optimal memory usage and time efficiency. Although many operations in Matlab are based on optimized BLAS routines, there are still parts of the implementation that runs very inefficiently, taking much more time than it would in an optimized implementation. We will therefore not report any timings, and we will not compare time-accuracy-efficiency of the methods here with other existing methods.
The Matlab implementation is based on the USEMe package, [74]. The USEMe package handles the complicated logic in loading the mesh and setting up the standard operators.
6.1 Zalesaks Problem
Zalesak [76] proposed a test to evaluate how well a method transports an interface.
The Zalesak disk is designed to mimic many of the interface configurations which
may appear in a simulation: It has smooth round areas, sharp corners, and inter-faces located close to each other. Testing on the Zalesak disk should therefore give a good idea of the performance of the method. The Zalesak disk, including the mesh which we will use, is depicted in Fig. 6.1. The disc is rotated one revolution
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
Figure 6.1: Zalesaks problem, including underlying mesh. Entire domain is [−4; 4]× [−4,4] hence only a part of the domain is shown
around the center of the domain using
φt−u· ∇φ= 0, u= (−πy, πx). (6.1)
We shall use the fourth order explicit Runge Kutta scheme of [10] for time step-ping Equation (6.1). The rotation is a shape preserving transformation, hence reinitialization is not necessary. In each Runge Kutta stage we will use a small amount of filtering to control oscillations coming from the discontinuities. Figure
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
Figure 6.2: Zalesaks problem, solution after one revolution for 3rd (left) and 5th (right) order elements.
6.2 shows the initial disk and the disk after one full revolution using 3rd and 5th order elements respectively. The mesh has 894 elements. For comparison with
6.1. Zalesaks Problem 101 other methods, counting the number of unknown, the 3rd order element has 10 and the 5th has 21 unknowns, in total 8940 and 18774 unknowns respectively, corresponding approximately to rectangular meshes of size 95x95 and 137x137.
Results using the 5th order elements show that the DG method can evolve the Zalesak disk well and solve it accurately without the need for any special techniques or tricks. Note that the width of the gap through the center of the circle is less than the side length of an element, and the disk itself is only resolved with about 30 elements, leading to a relatively coarse mesh.
Results using 3rd order elements show a typical degradation of the solution, when the DG method have difficulties representing the smallest scales of the level set. This may, however, be improved by using techniques presented in [66], designed to improve on the area conservation properties of the reinitialization.
The DG method works very well where the level set is smooth, on the circular part of the disk. Where the surface has sharp corners, the DG method smoothes the corner. In general we cannot expect high order of accuracy at such corners for the DG method, since the corners involve a discontinuity in the gradient of the level set, and high order cannot handle discontinuities well. On the other hand will surface tension quickly smoothen out such corners, and make sure that they will not reappear, at least for well resolved problems. Finally, the DG method is less accurate where two surfaces are very close, in the gap of the disc. This is again due to a discontinuity in the gradient of the level set, coming from the construction of the level set as a signed distance function.
Compared to lower order methods using a 137x137 mesh, the results are excel-lent, [31, 59, 66, 71, 76]