## Coupling to Elliptic Solvers

Very often, the speed of the interface is determined by solving an associated elliptic equation, e.g. the pressure equation for incompressible fluid flow. This leads to an elliptic equation which must be solved on an irregularly shaped domain or where there is an internal boundary with jump conditions across the boundary. There are several strategies to handle this problem. When using finite elements to solve this elliptic equation, a mesh is dynamically generated so that it conforms to this irregular boundary. When using finite differences, special delta functions can be added at nodes near the interface to enforce the jump conditions, see e.g. .

In the context of the level set method, there are three strategies for setting up and solving the associated elliptic equation. They vary in generality, Figure 4.17: Plot of the characteristic curves along which F is constant.

complexity, and accuracy, and provide different advantages. All three strategies are presented here.

### The Extended Finite Element Method

The extended finite element method (X-FEM) [29,81,121] is a numerical method to model internal (or external) boundaries without the need for the mesh to conform to these boundaries. The X-FEM is based on a standard Galerkin procedure and uses the concept of partition of unity  to accommodate the internal boundaries in the discrete model. The partition of unity method  generalized finite element approximations by presenting a means to embed local solutions of boundary-value problems into the finite element approximation.

For a standard finite element approximation, consider a point x of R that lies inside a finite element e. Denote the nodal set N = {n1; n2,..., nm j, where m is the number of nodes of element e. The approximation for a vector-valued function u(x) : Rd ^ Rd assumes the form uh(x) 4>I(x)u/, (ui e Rd), (4.53)

ni eN

where the functions (x) are the finite element basis functions and uI are the weights.

The extended finite element method uses enrichment functions, extra basis functions which are sensitive to prescribed boundaries, to capture the boundary conditions and improve the solution in the neighborhood of regions which would otherwise require greater spatial resolution. Consider again a point x that lies inside a finite element e. The enriched approximation for the function u(x) becomes uh(x) = £ fa(x)uI + £ \$j(x)f (x)aj, (4.54)

n eN nj eNg

classical enriched where the nodal set Ng consists of nodes which are on elements cut by the boundary, for example, see Fig. 4.18. In general, the choice of the enrichment function ^ (x) that appears in Eq. 4.54 depends on the geometry, the boundary condition, and the elliptic equation being solved.

To illustrate the effectiveness of this approach, consider the following simple example. Suppose we wish to solve the radial heat equation on an annulus given Figure 4.18: Example of choosing enriched nodes. Enriched nodes are indicated by gray dots.

The exact solution is given by u(r) = — 10e ln(r) + 10e ln(L). (4.57)

If we solve this equation for € = 0.01, L = 9 using a standard finite element method with linear elements and with nodes at r = 0,..., 9, the solution for € < r < 1 is very unsatisfactory, as shown in Fig. 4.19. However, by using a simple enrichment function ^1(r) = ln(r), and using this enrichment function on the first two nodes (located at r = 0, 1), dramatically better results are achieved (Fig. 4.19). Of course, refining the finite element mesh would also improve the results, but this requires remeshing as the interface (in this example the left boundary) moves. The X-FEM achieves this accuracy without remeshing.

The merits of coupling level sets to the extended finite element method were first explored in , and subsequently its advantages further realized in [53, 61,82,115,117,120]. The two methods make a natural pair of methods where:

1. Level sets provide greater ease and simplification in the representation of geometric interfaces.

2. The X-FEM, given the right enrichment functions, can accurately compute solutions of elliptic equations which are often required for computing the interface velocity.

3. Geometric computations required for evaluating the enrichment functions (such as the normal or the distance to the interface) are readily computed from the level set function .

4. The nodes to be enriched are easily identified using the signed distance construction of the level set function [115,117,118,120].

Compared to the other methods to follow, this algorithm is more complex, but it is also much more general. Through the use of enrichment functions, this method provides a much better solution near the interface, providing subgrid resolution in that region without requiring additional mesh refinement. This is

Finite Element Method Exact Solution

Extended Finite Element Method

01 23456789

Finite Element Method Exact Solution

Extended Finite Element Method

Figure 4.19: Solutions of the radial heat equation: (a) whole domain e < r < L and (b) across first three nodes.

important when having to interpolate the data to determine the front speed on the boundary contour. Of the three methods, this is the only one that has this capability.

### The Immersed Interface Method

The immersed interface method, introduced by LeVeque and Li , has also been coupled to the level set method [76,78]. Like the X-FEM described above, the immersed interface method is designed to solve elliptic equations which arise in a variety of physical applications. The advantage of the immersed interface method is that it is second-order accurate, even near the interface where jump conditions may appear.

The immersed interface method is designed to solve equations of the form where the coefficient functions 3, k, and f may have discontinuities across an interface r. The function f may also have a delta function singularity, which often arises, for example, from surface tension in multiphase flow.

The keyideainthe immersed interface methodisto modifythe discretization of Eq. 4.58 in such a way that the jump discontinuities and singularities are accounted for, leading to a fully second-order method. At points away from the interface, where the coefficient functions and the solution are smooth, the standard central difference approximation is used. However, for grid points which are near the interface, an additional grid point is added to the usual central difference stencil to account for a second-order Taylor approximation around a point on the interface.

To illustrate how this method works, consider the one-dimensional problem

(8 Ux) x u — f, x € [0, 1] \ a, u+ — u- — a, at x — a, u+ — u— — b, at x — a,

where u— is the value of u on the interval [0, a], and u+ is the value of u on the interval [a, 1]. Suppose that the point a is located between the uniformly spaced grid points xi and xi+1. The idea is to calculate coefficients yi—i, yi, yi+i, and an additional constant, Ci, so that the approximation

Yi-iUi-i + YiUi + Yi+i Ui+i + KiUi = fi + Ci (4.62)

is second-order accurate, with jump conditions Eqs. 4.60 and 4.61.

To determine the Yi's, Taylor expansions are taken about the point x = a to get u(Xi-1) = u- + (Xi-1 - a)u- + 1(xi-1 - a)u—x + O (Ax3), (4.63)

u(xi) = u + (xi - a)ux + ^(xi - a)2uxx + O (Ax3), (4.64)

u(Xi+1) = u+ + (Xi+1 - a)u+ + ^(Xi+1 - a)2u+x + O (Ax3). (4.65)

These expansions are inserted into Eq. 4.62, and the u+ terms are eliminated from the equation by using the jump conditions Eqs. 4.60 and 4.61, combined with the equation

which comes from the continuity of f in Eq. 4.59. The function f on the right side of Eq. 4.62 is replaced with the approximation from the left side, f = (Pu- )x + ku-. This results in the following equation:

( _ 1 2 / bpx - k a + Yi+1 ( u + a + (xi+1 - a)(ux + b) + ^(xi+1 - a)21 uxx--p-

The coefficients Yi-1, Yi, Yi+1, and Ci are now chosen so that Eq. 4.67 holds up to second order. This leads to the following equations:

Yi-1 (xi-1 - a) + Yi(x - a) + Yi+1 (xi+1 - a) + k(xi - a) = px, (4.69)

Yi-1 (xi-1 - a)2 + Yi(xi - a)2 + Yi+1 (xi+1 - a)2 + k(xi - a)2 = 2p, (4.70)

( , \ 1 (bPx - ka)(xi+1 - a)2\ Yi+1 I a + b(x+1 - a) - ^-p-) = Ci. (4.71)  Figure 4.20: Choice of stencil for (a) points not crossed by the interface and (b) points where the interface crosses the stencil. Dashed lines indicate the points used in the stencil.

These equations are solved for yi_1, Yi, Yi+h and Q, thus determining the numerical approximation corresponding to the point xi using Eq. 4.62. A similar process is followed for the approximation centered at xi+1. This results in a specialized discretization at these two points and standard central difference approximations everywhere else.

For higher dimensional problems, a similar approach is taken. At grid points not crossed by the interface, the standard central difference stencil is used (see Fig. 4.20(a)) to approximate Eq. 4.58. At grid points where the interface crosses through the stencil, an additional grid point is chosen across the interface from the center of the stencil (see Fig. 4.20(b)).

When building the specialized discretization for the stencil at grid points as in Fig. 4.20(b), a point (x*, y*) is chosen for the point around which the approximation will be computed, and around which all Taylor expansions will be taken. Usually, the point (x*, y* ) is the point on the interface closest to the center of the stencil (in this example, point 2). Once (x*, y*) is chosen, a coordinate transformation is taken so that the interface normal maps onto the x-axis. Once this coordinate transformation is completed, the computation of the stencil is similar to the one-dimensional case described above.

As noted earlier, the advantage of this method is that it is truly second-order accurate, even in the neighborhood of the interface. However, the stencil that is produced is irregular, and it sometimes can be difficult to solve the resulting linear system. Also, the choice of the points (x*, y*) is somewhat arbitrary, and it is not clear what the best choices should be. Nonetheless, the method has been used successfully in a number of applications, e.g. see the review in .

### The Ghost Point Method

The ghost point method  is another method designed to solve elliptic equations with irregular and moving boundaries represented by the level set method. The idea behind this method is similar to the use of what are often called ghost points for discretizing boundary conditions in finite difference methods. In this context, ghost points are grid points located outside the computational domain, and are used to enforce boundary conditions.

The method presented in  is designed to solve equations of the form

in an irregularly shaped domain where S and f are smooth functions defined on and g is defined on 3^, the boundary of This is a more restrictive class of problems than can be handled by the previous two methods described, but it is a class of problems that often arises. By focusing on this simpler class, a second-order method with a simple discretization can be employed, which uses a stencil that has properties which make it easier to solve numerically than the system created by the previous methods.

To illustrate this method, consider first the one-dimensional problem

with 3^ = xj, and u(xi) = uj. Assume xj lies between the two grid points xi and xi+1. For points xj in the interior of the domain, the central difference discretization, similar to the one used in the immersed interface method, is used:

At the boundary, the discretization Eq. 4.74 is again employed, but the value of ui+l is not defined because xi+l is outside of Instead, a ghost value for ui+l is computed from the boundary condition using a linear extrapolation:

0 Ax

For stability reasons, if 0 < Ax, then Eq. 4.75 is replaced with ui+1 = uI. Using Eq. 4.75 in Eq. 4.74 produces the following discretization for the point near the boundary:

In multiple dimensions, this same extrapolation technique is carried out along each coordinate direction.

The resulting discretization is only first-order accurate near the boundary, but is second-order accurate overall. This is due to the confinement of the firstorder error to the nodes adjacent to the boundary. On the other hand, the linear system that comes from this discretization can be solved using faster conjugate gradient-type algorithms. Increasing the order of the extrapolation to compute ui+1 can result in a linear system that is more difficult to solve numerically, because of the non-symmetric stencil, and hence is not preferred.

This method is used primarily for its simplicity, while still yielding second-order convergence overall. For problems where the accuracy at the boundary is critical, this is probably not the preferred method, especially if the solution is difficult to resolve near the boundary. The method has been used in a handful of applications, for example, see .

### Comparison of the Elliptic Equation Solvers

The algorithms presented here, for solving elliptic equations in conjunction with the level set method, vary significantly in sophistication, complexity, and capability. The X-FEM approach is by far the most difficult to construct, but is also the most general, and has the greatest potential to solve challenging problems. In particular, the X-FEM approach provides a much more accurate representation of the solution near the boundary, a property that is of critical importance when the velocity of the interface depends on this very value.

The immersed interface method and ghost point method, on the other hand, are built much more easily, and still produce accurate solutions. The immersed interface method handles a larger range of equations than does the ghost point method, which is the most restrictive in this regard. Between these two methods, the immersed interface method is more accurate at the boundary, but at the expense of a more difficult system of equations to solve numerically.

The ghost point method is probably the fastest, due to its use of faster linear solvers, but an actual direct comparison has not been done. Both the immersed interface method and ghost point method will be faster than the X-FEM approach on the same mesh. However, to obtain the same accuracy near the interface, the X-FEM will not require as fine a mesh as the others, and hence can make up the difference in time by using a coarser mesh to obtain comparable results. A direct comparison of these three methods is the subject of current research.

### 4.3.4 Particle Level Set Method

Another modification of the level set method, called the particle level set method, was proposed by Enright et al. in . In the particle level set method, the level set function is compared with the motion of particles which move along the characteristics of the same velocity field. For an interface which is passively advected using the same velocity field, the particles, in theory, should not cross the interface. By comparing the motion of the particles with the moving interface, problems with the location of the interface can be identified and corrected.

Suppose the interface velocity is determined by a velocity field v(x, t). Given this velocity, the interface speed function, F, in Eq. 4.5 is given by

Substituting this expression for F into Eq. 4.77 gives the passive interface ad-vection equation d\$

At the same time, the particles themselves are moving with this same velocity, v. These two evolutions are coupled together when the particles are checked to see if any has crossed the interface, which in this case indicates that a particle has moved from a point where \$ > 0 to a point where \$ < 0, or vice versa. At that point, the level set function is "corrected."

In , a large number of particles are randomly distributed uniformly in the neighborhood of the interface \$ = 0. Each particle, p, is assigned a sign, sp, to indicate whether it is starting where \$ > 0 or \$ < 0, and is also assigned its distance, rp, to the interface. As the evolution of the interface and the particles proceeds, the particle locations are periodically checked to determine whether they have strayed across the level set function interface.

When a particle is determined to have strayed sufficiently far across the level set interface, the interface is reconstructed using the particle information. To do this, each particle, p, located at the point xp, is assigned a local signed distance function dp(x) = sp(rp - ||x - xp||). (4.79)

The level set function is now reconstructed in two steps. First, the functions 0+ and 0- are computed where

and where P + and P- are the sets of points which were assigned positive and negative sp respectively. The final 0 function is now recovered from 0+ and 0-by the equation

where

There is no guarantee that the resulting reconstructed level set function will be a signed distance function, so if this is desired, a reinitialization step will be applied to reform 0 into a signed distance function.

What is novel about this approach is the use of the Lagrangian and Eulerian methods to play against each other to ensure proper interface motion. However, one must carefully determine when the particle solution is correct, versus the level set evolution. This is determined by checking the local characteristics to see if they are colliding or expanding. The level set evolution tends to be better when characteristics are colliding, whereas the particle method will be more reliable when the interface is moving tangentially or stretching. Nonetheless, this combination tries to extract the positive capabilities of both the Lagrangian and Eulerian types of approaches to interface motion, while discounting the negatives.