SOME NEW FINITE DIFFERENCE METHODS FOR HELMHOLTZ EQUATIONS ON IRREGULAR DOMAINS OR WITH INTERFACES

Solving a Helmholtz equation Δ u + λ u = f efficiently is a challenge for many applications. For example, the core part of many efficient solvers for the incompressible Navier-Stokes equations is to solve one or several Helmholtz equations. In this paper, two new finite difference methods are proposed for solving Helmholtz equations on irregular domains, or with interfaces. For Helmholtz equations on irregular domains, the accuracy of the numerical solution obtained using the existing augmented immersed interface method (AIIM) may deteriorate when the magnitude of λ is large. In our new method, we use a level set function to extend the source term and the PDE to a larger domain before we apply the AIIM. For Helmholtz equations with interfaces, a new maximum principle preserving finite difference method is developed. The new method still uses the standard five-point stencil with modifications of the finite difference scheme at irregular grid points. The resulting coefficient matrix of the linear system of finite difference equations satisfies the sign property of the discrete maximum principle and can be solved efficiently using a multigrid solver. The finite difference method is also extended to handle temporal discretized equations where the solution coefficient λ is inversely proportional to the mesh size.


Introduction
The Helmholtz equation (1) is used to describe many mathematical and physical applications including electromagnetics, wave propagation, potential theory, and many others. It has also been used in time stepping numerical schemes for diffusion-reaction systems and for Navier-Stokes equations. For example, the projection method [2,3,8] and the stream-function vorticity method [4,18,20] for the incompressible Navier-Stokes equations require solving two or more such Helmholtz equations at each time step in which λ ~ O(1/Δt), where Δt ~ h is the time step size, and h is mesh size in the space. While there are many existing fast algorithms in the literature for solving Helmholtz equations, it is still challenging to solve the problem on irregular domains or with an interface across which the coefficient λ is discontinuous.
In this paper, we propose two new finite difference methods for solving the Helmholtz equation. One is for irregular domains, either interior or exterior in a bounded domain, see Fig. 1(a) for an illustration; the other one is for interface problems in which the coefficient λ can have a finite jump across an arbitrary interface, see Fig. 1(b). Both methods use Cartesian grids and can handle the case when the magnitude of λ is large, say O(1/Δt). This is significant for Navier-Stokes solvers using Helmholtz equations. We assume that the boundary or the interface is smooth, that is, Γ ∈ C 2 .
Note that, in [6,7,10,17], the authors have developed augmented immersed interface method (AIIM) for Poisson/Helmholtz equation on irregular domains. The method works well with second order accuracy in the L ∞ norm for fixed and modest λ. But the accuracy deteriorates as the magnitude of λ gets bigger. From analysis, we know that this is due to a term O(hλ) which is O(1) when λ ~ 1/h in the local truncation error. Our new idea in this paper is to extend the entire equation, or more precisely the source term f to a wider domain via a level set function to offset the order-one error. An augmented variable, the jump in the normal derivative of the solution is introduced along the boundary Γ so that a fast Poisson solver based on FFT can be utilized. Numerical results show that our method is indeed second order accurate in the L ∞ norm, and it requires only 4 to 5 calls of fast Poisson solvers independent of the mesh size for the test problem. This is because the condition number of the coefficient matrix for the augmented variable is order one. This has been proved in [23] via the operator theory.
For the interface problem with a jump in λ, our strategy is to modify the finite difference scheme dimension by dimension when the interface cuts through the stencils. Thus we still have five point stencil. The finite difference scheme is constructed in such a way that the discrete maximum principle is preserved. The new method is much simplified compared with the one developed in [13] that is designed for more general problems. The new finite difference scheme provides an alternative explicit approach to discretize the Helmholtz equation with interfaces. This paper is organized as follows. In Section 2, we introduce the level set representation of the boundary or interface. We also explain how to carry out the local coordinate transformations and how to evaluate the normal and tangential directions, and the curvature information from the level set representation. In Section 3, we describe the fast solver for Helmholtz equation on irregular domain. In Section 4, we derive the maximum principle preserving scheme for the Helmholtz equation where both λ and f can be piecewise smooth functions across the interface. In Section 4.1, we develop an important extension of the new scheme to handle cases where λ ~ O (1/h). Discussion and conclusions are in Section 5.

The level set representation and related preparations
We use the zero level set function of a Lipschitz continuous function φ(x, y) to represent the irregular boundary or interface [19,21]. That is, (2) Let x * = (x * , y * ) be a point on Γ. We can build a local Cartesian coordinate system (x, ỹ) at x * by using the unit normal direction n(x * ) pointing from Ω − to Ω + and the unit tangential direction τ(x * ) counter-clockwise of n(x * ). See Fig. 2 for an illustration of the local Cartesian coordinate system. The coordinate transformation matrix A is (3) Note that A is a special orthogonal matrix with detA = 1, and A −1 = A T . The coordinate transformation rule is (4) In the neighborhood of x * (or 0̃ in the local coordinates), the interface curve Γ can be parameterized as (5) where χ is a smooth function of ỹ. It is important to know that by definition of the local coordinates, (6) The domain Ω is assumed to be a rectangle and can be discretized as (7) with a uniform mesh size h.
For a given grid point x = (x i , y j ), we define (8) We call (x i , y j ) an irregular grid point in reference to the central five point stencil if . For an irregular grid point, there exists an orthogonal projection point x * = (x * , y * ) on the interface satisfying (10) where α is to be determined and (11) is the unit normal direction at x. The parameter α can be determined approximately by solving the following quadratic equation explicitly (12) where H is the Hessian matrix (13) All the quantities related to φ and its derivatives are evaluated at the grid point (x i , y j ).
The unit tangential direction τ at x can be calculated as (14) To calculate the surface derivative , we differentiate the level set function φ(χ(ỹ), ỹ) = 0 twice with respect to ỹ and rearrange terms, then we can get [15,22] (15) Note that from the coordinate transformation rule, we have (16) and (17) It is also noted that (18) where is the curvature and used in the immersed interface method. The above formula offers another way to calculate the surface derivatives without using the local coordinate system.

A fast solver of the Helmholtz equation on an irregular domain
We consider the Helmholtz equation (19) (20) where q(u, u n ) is a linear functional of u and that includes the Dirichlet and Neumann boundary conditions. For an exterior bounded irregular domain, we also need to prescribe the boundary condition on the outer boundary ∂R, see Fig. 1(a) for an illustration.
Note that, in [6,7,10,17], the authors have developed the augmented immersed interface method for Poisson/Helmholtz equations on irregular domains. The method works well with second order accuracy in the L ∞ norm for fixed and modest |λ|. But the accuracy deteriorates as the magnitude of λ gets bigger even for λ < 0 since the local truncation errors have an order of O(hλ) at irregular grid points. The main contribution of our work here is to develop a new method that can also work for large |λ|.
In order to use a fast Poisson/Helmholtz solver such as the one from FISHPACK [1], we embed the domain Ω into a rectangle domain R ⊃ Ω. The PDE and the source term are extended to the entire rectangle domain R: (21) ( 22) where f e is an extension of f(x, y) to the entire domain. In other words, we change the original problem to an interface problem1. The solution u to the interface problem above is a functional of g, and g is one dimension lower than u. We determine g such that the solution u(g) satisfies the boundary condition q(u, u n ) = 0. Note that, given g, we can solve the above interface problem using the immersed interface method with a single call to a fast Poisson solver.

Extension of the source term
If |λ| is modest, we can simply take f e = 0 as in the second order accurate augmented immersed interface method [14]. However, the error constant depends on |λ| and the accuracy will deteriorate as |λ| gets large. This can be seen in the interface relation of the second order partial derivative, see (3.19) on page 40 of [14]. Our strategy to overcome this difficulty is to extend the PDE and the source term beyond the boundary Γ. There are several ways to extend the source term off the front. One is the fast marching method, see [5,21], the other one is the PDE approach using the following evolution equation (23) 1 We have also tested using [u] = 0, [u n ] = g. The approach presented in this paper gives better results especially in terms of the where t is an artificial variable, see for example [9,10,16,17]. The sign is determined from the problem whether it is exterior or interior problem. Note that is the unit normal direction. In our numerical test, we simply use the upwinding scheme to solve (23). The cost is about O(δN) operations, where δ is the width of the extension.
After the extension, there is no such O(hλ) term in the local truncation error anymore.
Nevertheless, we still need to interpolate the boundary condition in the augmented approach as described in the following subsection.

Augmented IIM and GMRES iterative method
In the discrete form, denote G as the approximation of g at the orthogonal projections of irregular grid points from the φ ≥ 0 side. The finite difference equations given G have the following form, (24) where C ij is the correction term that can be determined by the formula described in [14] at irregular grid points. The correction term is zero at regular grid points. In the matrix-vector form, it can be written as (25) where A is the discrete Helmholtz operator, BG corresponds to the correction terms given the jump in [u] at the discrete points. The system of equations can be solved by one step fast Poisson solver.
We should choose G such that the boundary condition is satisfied at the discrete points. This is done using an interpolation scheme that can be written as (26) where C and D are two matrices, and F 2 is a vector. The residual vector (27) is a measurement that describes how well the boundary condition is approximated by given G.
If we put those two matrix-vector equations (25) and (26) together, we get (28) In practice, we do not necessarily form the matrices A, B, C, and D. Note that the dimension of G is O(m) (assuming m ~ n), which is much smaller than that of U which is O(n 2 ).
Eliminating U from Eq. (28), we get the Schur complement system for G: Note that E = (D−CA −1 B) is not symmetric in general. We can either use a direct method, for example, the LU decomposition; or an iterative method, for example, the GMRES method, to solve the linear system (29) for G.
Using the GMRES method, first, we set G = 0 and then solve the Poisson equation. The residual of the linear system (29) (or the difference between the exact and the computed boundary condition), is actually the right hand side of the Schur complement with an opposite sign. This is because which gives the right hand side of the Schur complement system with an opposite sign.
Next, we explain how to find the matrix-vector multiplication given G. This again involves only two steps: (1) solving (25) by given G, to get U(G); (2) interpolating U(G) at Γ via the least squares interpolation. Once we know the matrix-vector multiplication, we can apply the GMRES or other iterative method directly. The procedure is illustrated in the following derivation: Thus the matrix-vector multiplication is the result of the difference of the residual of the interface condition that is to be interpolated.
The idea of the augmented method has been explained in [12,14] for elliptic interface problems with piecewise constant coefficient, and Poisson equations on irregular domains.

A numerical example for the Helmholtz equation on an irregular domain
We show an example of an exterior bounded irregular domain with a Neumann boundary condition, which is usually more difficult compared with problems with a Dirichlet boundary condition. We construct an exact solution: (40) The exterior domain is composed of two pieces: the boundary ∂Ω 1 of the rectangle: −1 ≤ x, y ≤ 1; and the boundary of ∂Ω 2 of the ellipse (41) The Dirichlet boundary condition on the rectangle −1 ≤ x, y ≤ 1, the Neumann boundary boundary on the ellipse, and the source term are determined from the exact solution.
In Table 1, we show the grid refinement analysis for the Helmholtz equation with large |λ|. We see average second order convergence in the infinity norm. The rate of convergence is defined as (42) And the number of GMRES iterations is denoted as k in the tables. For the modified Helmholtz equation where λ < 0, the condition number of the coefficient matrix is almost a constant (cond(A) ~ 10) independent of the mesh size. In this case, the coefficient matrix is positive definite and diagonally dominant. The number of GMRES iteration with tolerance 10 −6 is also a constant between 4 and 5. In [23], some theoretical analysis based on the operator theory explains why the condition number of coefficient matrix is a constant. In our new method, only the right hand side of the equations are modified, thus the condition number of the coefficient matrix remains the same. Thus the number of the GMRES iterations is independent of the mesh and λ if λ < 0. The number of the GMRES iterations does depends on the geometry and λ if λ > 0. Without extension, we can see that the errors are much larger with a rough first order convergence. For Helmholtz equation with λ > 0, it is well-known that a reasonable mesh is needed to resolve the solution. The solution may not exist if λ happens to be one of eigenvalues of the coefficient matrix.
In the example above, the right hand side f(x, y) is from the exact solution and the PDE, that is, f(x, y) = Δu e + λu e , where the u e is the true solution. Now we test another example in which the solution depends on λ. We choose f(x, y) = Δu, where u(x, y) is defined in (40).
The Dirichlet boundary condition along ∂R and the Neumann boundary condition along the ellipse are also from u(x, y) in (40). In this case, we do not have the true solution and the rate of convergence r cannot be assessed directly. Thus we compare the computed solution with the one obtained from a fine mesh 512 × 512 and present the error ratio instead. The error ratio is defined as the ratio of two consecutive errors. In Table 2, we show the grid refinement analysis result. We can see that the error ratio is close to number five, which indicates a true second order convergence, see page 81 of [11].

A maximum principle preserving solver for the modified Helmholtz equation with interfaces
In this section, we describe a maximum principle preserving solver for the 2D modified Helmholtz equation with interfaces, i.e. (43) with known jump conditions across an interface Γ, Both λ and f can have discontinuities across the interface Γ.
Grid points in this new finite difference scheme are categorized as either regular or irregular grid points. A grid point (x i , y j ) is called a regular grid point if the difference stencil centered at (x i , y j ) contains no grid points on the other side of the interface. For those majority regular grid points, we use the standard five-point stencil, i.e. (45) At those minority irregular grid points either close to or on the interface Γ, to derive the new scheme, we first locate the orthogonal projection point x * ∈ Γ of the irregular scheme point x = (x i , y j ) and define a local Cartesian coordinate system (x, ỹ) centered at x * with the coordinate transformation matrix A as defined in Section 2. Next we show the details of deriving finite difference equation for these irregular grid points.
For an irregular grid point (x i , y j ) ∈ Ω + , for example, without loss of generality, we assume that of the other four stencil grid points, only (46) in the five-point stencil. We also denote the local coordinates of (x i , y j−1 ) as (x̃i ,j−1 , ỹ i,j−1 ). If we still choose to use the standard five-point stencil, the truncation error reads (47) By applying the IIM with the jump conditions [u] and , we have [11,22] (48) where 0 means the evaluation at the origin in the local Cartesian coordinate system, xĩ ,j−1 is the x̃ coordinate of the grid point (x i , y j−1 ) in the local Cartesian coordinate system, and C + is the following correction term from the IIM [11,22]: Thus far, we have derived the local truncation error in the framework of the IIM. Now we will modify the five-point stencil at irregular grid points to satisfy both the sign property of the discrete maximum principle, and at the same time, to keep the local truncation error as O(h) with suitable correction terms.
The new scheme modifies the five-point stencil by taking the sign of [λ] into account. If [λ] (0) > 0, or λ + > λ − , it is easy to see that since (50) and the projection point (x * , y * ) should stay inside the box of [ Therefore, by plugging in Eq. (51, 52), the local truncation error can be written as (53) This means that we can modify the five-point scheme by changing the coefficient for and the local truncation error becomes T = C + + O(h). Note that the sign property for the discrete maximum principle is automatically satisfied since .
On the other hand, if [λ] (0) < 0, or λ + < λ − , after applying the jump condition [u] and following similar steps as in the case above, we have [11,22] (54) Therefore we can modify the five-point stencil by changing the coefficient for , with the correction term . Again the sign property is satisfied for the discrete maximum principle. This concludes the derivation of the new finite difference scheme for the irregular grid point (x i , y j ) ∈ Ω + . It is understood that corrections need to be made for every stencil point of (x i , y j ) which is not in Ω + . Since there are no interaction terms between those stencil points, all corrections are simply additive to each other.
Derivations of the modified five-point stencil at irregular grid points in Ω − are easily replicated from the case above where irregular grid points are in Ω + . For the sake of completeness, we still show major steps and results as follows. For an irregular grid point (x i , y j ) ∈ Ω − , we assume that only (55) in the five-point stencil. The local truncation error reads (56) Again by applying the IIM with the known jump conditions [u] and , we have [11,22] (57) Note that C − = −C + which reects symmetry of the interface problems. As we have shown so far, regardless of the sign of [λ](0) and the location of the irregular grid points ∈ Ω ± , the new finite difference scheme modified from the standard five-point stencil and based on the IIM always satisfies the sign property of the discrete maximum principle. The local truncation error for the minority irregular grid points is O(h) with certain correction terms. Since Γ is codimensional one object of the 2D domain Ω, it is known that for elliptic equations, the maximum norm error of discretized solution is still O(h 2 ) if the local truncation error at the irregular grid points is O(h).
There are various ways to invert the formed coefficient matrix, most part of which is the same as the one from the standard five-point stencil, with modifications only around the interface. We choose to solve the discretized equations by using the multigrid solver mgd9v [24], which is designed for standard nine-point finite difference stencil. Tolerance of the multigrid solver is set to be 10 −6 unless otherwise specified.

Extended scheme for temporal discretized equations
In Section 4, the maximum principle preserving scheme for the modified Helmholtz equation with interfaces is derived using a first-order grid point approximation of either u + (0) or u − (0). If we apply the scheme to solve temporal discretized equations arising from e.g. Navier-Stokes equations, note that λ is proportional to where μ is the kinematic viscosity of the fluid and Δt is the time step size [15,22]. It is expected that Δt is the same across the interface, therefore [λ] is proportional to . Usually the time step size Δt is in the order of grid width h. It is obvious that the first-order approximation of u + (0) or u − (0) is not sufficient to maintain the local truncation error as O(h). A straightforward extension of the new scheme is to leverage a second-order approximation of u + (0) or u − (0) for the temporal discretized equations where . We outline the extended scheme as follows.
For a given irregular grid point (x i , y j ), we use its five scheme points including (x i , y j ) itself to form the singular value decomposition (svd) approximation of u + (0) or u − (0) in the following form [15,22]: (58) where a m,n are the coefficients to be determined and b is a correction term resulting from the fact that interpolation grid points are from both sides of the interface Γ.
Using the local Cartesian coordinate system centered at x * ∈ Γ, the local coordinates for the interpolation grid points are denoted as (x̃m ,n , ỹ m,n ). The Taylor expansion of u m,n at the local origin is (59) The ± sign is determined according to whether (x m , y n ) ∈ Ω ± . The coefficients a m,n in the interpolation scheme are the solution of the following 3 × 5 linear system of equations: (60) The correction term to approximate u + (0) is (61) To calculate the correction term b, note that we have Since the system of equations for a m,n is under-determined, we use the svd method to solve the equations for the coefficients a m,n .
After the five coefficients a m,n is solved using the svd method, one can modify the finite difference scheme coefficients for the five scheme points in a similar fashion as shown in Section 4. And the correction term also needs to be updated using the value of b from the the svd approximations. It is not clear whether the maximum principle is still preserved for the extended scheme. However, it is noted that the svd solutions of a m,n only depend on the local geometry of the interfaces. Numerical examples in the following section show that the extended scheme can be solved efficiently using the same multigrid solver mgd9v [24].
The exact solution is constructed to be (66) where A, B, C, D, μ 1 , μ 2 are constants. And (67) Note that to calculate , it is obvious that (68) and (69) Then we have (70) In the following numerical examples, we first show the grid refinement and efficiency analyses of (1) the maximum principle preserving solver for the modified Helmholtz equation with interfaces, and (2) the extended solver for the modified Helmholtz equation with . Next we also demonstrate the applicability of the extended solver for Helmholtz equation with interfaces where λ is positive. Number of multigrid iterations (k) and solver execution time (T) in seconds are included for the three examples. Example 1. We choose the following sets of parameters for the exact solution in Ω + and Ω − : (71) or (72) The error analyzed is the ℒ ∞ error at all grid points. See Fig. 5 for the jump conditions of the solution and Fig. 6 for an example of computed solution u(x, y) with grid resolution of 81 × 81 and with parameters as in (71). The computed first and second order surface derivatives at orthogonal projection points are also shown in Fig. 7. Example 2. Same sets of parameters as in Example 1 are chosen for Example 2, except that λ is selected to be inversely proportional to the grid width h, i.e. (73) or (74) Note that in this example, the grid refinement analysis is different from usual case in that the differential equation actually depends on grid width. Example 3. The extended scheme is also applicable to the Helmholtz equation where λ is positive. We demonstrate the extended scheme by solving the Helmholtz equation with same sets of parameters as in Example 1, except that λ is selected to be positive numbers instead of negative numbers, i.e.

Conclusions and acknowledgments
In this paper, we have proposed some new methods for solving Helmholtz equations on irregular domains or with interface, particularly for the coefficient λ ≤ 0. The new methods are motivated for the case when the magnitude of λ is large and other methods may not perform well. It is noted that for the irregular domain problem, the number of calls to a fast Poisson solver of our new method is a small constant independent of the mesh size and λ. For interface problems, our new method preserve the discrete maximum principle without using an optimization solver.    Computed second-order surface derivatives at orthogonal projection points on the interface curve in the middle represented by the level set function.   Computed first and second order surface derivatives of jump condition w of the solution u(x, y) across the interface curve in the middle. Table 1 A grid refinement analysis for the irregular domain problem with a = 0.5, b = 0.4. The Dirichlet boundary condition is prescribed on ∂R, and the Neumann boundary condition is prescribed on Γ.  A grid refinement analysis for the irregular domain problem with a = 0.5, b = 0.4. The Dirichlet boundary condition is prescribed on ∂R, and the Neumann boundary condition is prescribed on Γ using u(x, y) in (40).  A grid refinement analysis for Example 1.   Table 4 A grid refinement analysis for Example 2.   A grid refinement analysis for Example 3.