A Correction Function Method for Poisson Problems with Interface Jump Conditions.
Abstract
In this paper we present a method to treat interface jump conditions for constant coefficients Poisson problems that allows the use of standard “black box” solvers, without compromising accuracy. The basic idea of the new approach is similar to the Ghost Fluid Method (GFM). The GFM relies on corrections applied on nodes located across the interface for discretization stencils that straddle the interface. If the corrections are solutionindependent, they can be moved to the righthandside (RHS) of the equations, producing a problem with the same linear system as if there were no jumps, only with a different RHS. However, achieving high accuracy is very hard (if not impossible) with the “standard” approaches used to compute the GFM correction terms.
In this paper we generalize the GFM correction terms to a correction function, defined on a band around the interface. This function is then shown to be characterized as the solution to a PDE, with appropriate boundary conditions. This PDE can, in principle, be solved to any desired order of accuracy. As an example, we apply this new method to devise a 4 order accurate scheme for the constant coefficients Poisson equation with discontinuities in 2D. This scheme is based on (i) the standard 9point stencil discretization of the Poisson equation, (ii) a representation of the correction function in terms of bicubics, and (iii) a solution of the correction function PDE by a least squares minimization. Several applications of the method are presented to illustrate its robustness dealing with a variety of interface geometries, its capability to capture sharp discontinuities, and its high convergence rate.
keywords:
Poisson equation, interface jump condition, Ghost Fluid method, Gradient augmented level set method, High accuracy, Hermite cubic splinePacs:
47.11j, 47.11.BcMsc:
[2010] 76M20, 35N06sort&compress \newdefinitionrmkRemark
1 Introduction.
1.1 Motivation and background information.
In this paper we present a new and efficient method to solve the constant coefficients Poisson equation in the presence of discontinuities across an interface, with a high order of accuracy. Solutions of the Poisson equation with discontinuities are of fundamental importance in the description of fluid flows separated by interfaces (e.g. the contact surfaces for immiscible multiphase fluids, or fluids separated by a membrane) and other multiphase diffusion phenomena. Over the last three decades, several methods have been developed to solve problems of this type numerically peskin:77; sussman:94; leveque:94; leveque:97; johansen:98; fedkiw:99; fedkiwetal:99; kang:00; liu:00; lai:00; li:01; nguyen:01; lee:03; gibou:07; gong:08; dolbow:09; bedrossian:10. However, obtaining a high order of accuracy still poses great challenges in terms of complexity and computational efficiency.
When the solution is known to be smooth, it is easy to obtain highly accurate finitedifference discretizations of the Poisson equation on a regular grid. Furthermore, these discretizations commonly yield symmetric and banded linear systems, which can be inverted efficiently trefthen:97. On the other hand, when singularities occur (e.g. discontinuities) across internal interfaces, some of the regular discretization stencils will straddle the interface, which renders the whole procedure invalid.
Several strategies have been proposed to tackle this issue. Peskin peskin:77 introduced the Immersed Boundary Method (IBM) peskin:77; lai:00, in which the discontinuities are reinterpreted as additional (singular) source terms concentrated on the interface. These singular terms are then “regularized” and appropriately spread out over the regular grid — in a “thin” band enclosing the interface. The result is a first order scheme that smears discontinuities. In order to avoid this smearing of the interface information, LeVeque and Li leveque:94 developed the Immersed Interface Method (IIM) leveque:94; leveque:97; li:01; lee:03, which is a methodology to modify the discretization stencils, taking into consideration the discontinuities at their actual locations. The IIM guarantees second order accuracy and sharp discontinuities, but at the cost of added discretization complexity and loss of symmetry.
The new method advanced in this paper builds on the ideas introduced by the Ghost Fluid Method (GFM) mayo:84; fedkiw:99; fedkiwetal:99; liu:00; kang:00; nguyen:01; gibou:07. The GFM is based on defining both actual and “ghost” fluid variables at every node on a narrow band enclosing the interface. The ghost variables work as extensions of the actual variables across the interface — the solution on each side of the interface is assumed to have a smooth extension into the other side. This approach allows the use of standard discretizations everywhere in the domain. In most GFM versions, the ghost values are written as the actual values, plus corrections that are independent of the underlying solution to the Poisson problem. Hence, the corrections can be precomputed, and moved into the source term for the equation. In this fashion the GFM yields the same linear system as the one produced by the problem without an interface, except for changes in the righthandside (sources) only, which can then be inverted just as efficiently.
The key difficulty in the GFM is the calculation of the correction terms, since the overall accuracy of the scheme depends heavily on the quality of the assigned ghost values. In fedkiw:99; fedkiwetal:99; liu:00; kang:00; nguyen:01 the authors develop first order accurate approaches to deal with discontinuities. In the present work, we show that for the constant coefficients Poisson equation we can generalize the GFM correction term (at each ghost point) concept to that of a correction function defined on a narrow band enclosing the interface. Hence we call this new approach the Correction Function Method (CFM). This correction function is then shown to be characterized as the solution to a PDE, with appropriate boundary conditions on the interface — see § 4. Thus, at least in principle, one can calculate the correction function to any order of accuracy, by designing algorithms to solve the PDE that defines it. In this paper we present examples of 2 and 4 order accurate schemes (to solve the constant coefficients Poisson equation, with discontinuities across interfaces, in 2D) developed using this general framework.
A key point (see § 5) in the scheme developed here is the way we solve the PDE defining the correction function. This PDE is solved in a weak fashion using a least squares minimization procedure. This provides a flexible approach that allows the development of a robust scheme that can deal with the geometrical complications of the placement of the regular grid stencils relative to the interface. Furthermore, this approach is easy to generalize to 3D, or to higher orders of accuracy.
1.2 Other related work.
It is relevant to note other developments to solve the Poisson equation under similar circumstances — multiple phases separated by interfaces — but with different interface conditions. The Poisson problem with Dirichlet boundary conditions, on an irregular boundary embedded in a regular grid, has been solved to second order of accuracy using fast Poisson solvers mayo:84, finite volume johansen:98, and finite differences udaykumar:99; gibou:02; jomaa:05; gibou:05 approaches. In particular, Gibou et al. gibou:02 and Jomaa and Macaskill jomaa:05 have shown that it is possible to obtain symmetric discretizations of the embedded Dirichlet problem, up to second order of accuracy. Gibou and Fedkiw gibou:05 have developed a fourth order accurate discretization of the problem, at the cost of giving up symmetry. More recently, the same problem has also been solved, to second order of accuracy, in nongraded adaptive Cartesian grids by Chen et al. chen:07. Furthermore, the embedded Dirichlet problem is closely related to the Stefan problem modeling dendritic growth, as described in sethian:92; chen:97.
The finiteelement community has also made significant progress in incorporating the IIM and similar techniques to solve the Poisson equation using embedded grids. In particular the works by Gong et al. gong:08, Dolbow and Harari dolbow:09, and Bedrossian et al. bedrossian:10 describe second order accurate finiteelement discretizations that result in symmetric linear systems. Moreover, in these works the interface (or boundary) conditions are imposed in a weak fashion, which bears some conceptual similarities with the CFM presented here, although the execution is rather different.
1.3 Interface representation.
Another issue of primary importance to multiphase problems is the representation of the interface (and its tracking in unsteady cases). Some authors (see peskin:77; mayo:84; udaykumar:99) choose to represent the interface explicitly, by tracking interface particles. The location of the neighboring particles is then used to produce local interpolations (e.g. splines), which are then applied to compute geometric information — such as curvature and normal directions. Although this approach can be quite accurate, it requires special treatment when the interface undergoes either large deformations or topological changes — such as mergers or splits. Even though we are not concerned with these issues in this paper, we elected to adopt an implicit representation, to avoid complications in future applications. In an implicit representation, the interface is given as the zero level of a function that is defined everywhere in the regular grid — the level set function osher:88. In particular, we adopted the GradientAugmented Level Set (GALS) method nave:10. With this extension of the level set method, we can obtain highly accurate representations of the interface, and other geometric information, with the additional advantage that this method uses only local grid information. We discuss the question of interface representation in more detail in § 5.6.
1.4 Organization of the paper.
The remainder of the paper is organized as follows. In § 2 we introduce the Poisson problem that we seek to solve. In § 3 the basic idea behind the solution method, and its relationship to the GFM, are explored. Next, in § 4, we introduce the concept of the correction function and show how it is defined by a PDE problem. In § 5 we apply this new framework to build a 4 order accurate scheme in 2D. Since the emphasis of this paper is on highorder schemes, we describe the 2 order accurate scheme in appendix C. Next, in § 6 we demonstrate the robustness and accuracy of the 2D scheme by applying it to several example problems. The conclusions are in § 7. In appendix A we review some background material, and notation, on bicubic interpolation. Finally, in appendix B we discuss some technical issues regarding the construction of the sets where the correction function is solved for.
2 Definition of the problem.
Our objective is to solve the constant coefficients Poisson’s equation in a domain in which the solution is discontinuous across a codimension 1 interface , which divides the domain into the subdomains and , as illustrated in figure 1. We use the notation and to denote the solution in each of the subdomains. Let the discontinuities across be given in terms of two functions defined on the interface: for the jump in the function values, and for the jump in the normal derivatives. Furthermore, assume Dirichlet boundary conditions on the “outer” boundary (see figure 1). Thus the problem to be solved is
(1a)  
(1b)  
(1c)  
(1d) 
where
(2a)  
(2b) 
Throughout this paper, is the spatial vector (where , or ), and is the Laplacian operator defined by
(3) 
Furthermore,
(4) 
denotes the derivative of in the direction of , the unit vector normal to the interface pointing towards (see figure 1).
It is important to note that our method focuses on the discretization of the problem in the vicinity of the interface only. Thus, the method is compatible with any set of boundary conditions on , not just Dirichlet.
3 Solution method – the basic idea.
To achieve the goal of efficient high order discretizations of Poisson equation in the presence of discontinuities, we build on the idea of the Ghost Fluid Method (GFM). In essence, we use a standard discretization of the Laplace operator (on a domain without an interface ) and modify the righthandside (RHS) to incorporate the jump conditions across . Thus, the resulting linear system can be inverted as efficiently as in the case of a solution without discontinuities.
Let us first illustrate the key concept in the GFM with a simple example, involving a regular grid and the standard second order discretization of the 1D analog of the problem we are interested in. Thus, consider the problem of discretizing the equation in some interval , where is discontinuous across some point — hence (respectively ) is the domain (respectively ). Then, see figure 2, when trying to approximate at a grid point such that , we would like to write
(5) 
where is the grid spacing. However, we do not have information on , but rather on . Thus, the idea is to estimate a correction for , to recover , such that Eq. (5) can be applied:
(6) 
where is the correction term. Now we note that, if can be written as a correction that is independent on the solution , then it can be moved to the RHS of the equation, and absorbed into . That is
(7) 
This allows the solution of the problem with prescribed discontinuities using the same discretization as the one employed to solve the simpler problem without an interface — which leads to a great efficiency gain.
The error in estimating is crucial in determining the accuracy of the final discretization. Liu, Fedkiw, and Kang liu:00 introduced a dimensionbydimension linear extrapolation of the interface jump conditions, to get a first order approximation for . Our new method is based on generalizing the idea of a correction term to that of a correction function, for which we can write an equation. One can then obtain high accuracy representations for by solving this equation, without the complications into which dimensionbydimension (with Taylor expansions) approaches run into.
An additional advantage of the correction function approach is that can be calculated at any point near the interface . Hence it can be used with any finite differences discretization of the Poisson equation, without regard to the particulars of its stencil (as would be the case with any approach based on Taylor expansions).
4 The correction function and the equation defining it.
As mentioned earlier, the aim here is to generalize the correction term concept to that of a correction function, and then to find an equation (a PDE, with appropriate boundary conditions) that uniquely characterizes the correction function. Then, at least in principle, one can design algorithms to solve the PDE in order to obtain solutions to the correction function of any desired order of accuracy.
Let us begin by considering a small region enclosing the interface , defined as the set of all the points within some distance of , where is of the order of the grid size . As we will see below, we would like to have as small as possible. On the other hand, has to include all the points where the CFM requires corrections to be computed, which means^{1}^{1}1For the particular discretization of the Laplace operator that we use in this paper. that cannot be smaller than . In addition, algorithmic considerations (to be seen later) force to be slightly larger than this last value.
Next, we assume we that can extrapolate both and , so that they are valid everywhere within , in such a way that they satisfy the Poisson equations
(8a)  
(8b) 
where and are smooth enough (see remark 4 below) extensions of the source term to . In particular, notice that the introduction of and allows the possibility of the source term changing (i.e. a discontinuous source term) across . The correction function is then defined by .
Taking the difference between the equations in (8), and using the jump conditions (1b1c), yields
(9a)  
(9b)  
(9c) 
This achieves the aim of having the correction function defined by a set of equations, with some provisos — see remark 4 below. Note that:

If , for , then , for .

Equation (9c) imposes the true jump condition in the normal direction, whereas some versions of the GFM rely on a dimensionbydimension approximation of this condition (see Ref. liu:00).
The smoothness requirement on and is tied up to how accurate an approximation to the correction term is needed. For example, if a 4 order algorithm is used to find , this will (generally) necessitate to be at least for the errors to actually be 4 order. Hence, in this case, must be .
Equation (9) is an elliptic Cauchy problem for in . In general, such problems are illposed. However, we are seeking for solutions within the context of a numerical approximation where

There is a frequency cutoff in both the data and , and the description of the curve .

We are interested in the solution only a small distance away from the interface , where this distance vanishes simultaneously with the inverse of the cutoff frequency in point (a).
What (a) and (b) mean is that the arbitrarily large growth rate for arbitrarily small perturbations, which is responsible for the illposedness of the Cauchy problem in Eq. (9), does not occur within the special context where we need to solve the problem. This large growth rate does not occur because, for the solutions of the Poisson equation, the growth rate for a perturbation of wave number along some straight line, is given by — where is the distance from the line. However, by construction, in the case of interest to us is bounded.
Let us be more precise, and define a number characterizing how well posed the discretized version of Eq. (9) is, by
where growth is defined relative to the size of a perturbation to the solution on the interface. This number is determined by (the “radius” of ) as the following calculation shows: First of all, there is no loss of generality in assuming that the interface is flat, provided that the numerical grid is fine enough to resolve . In this case, let us introduce an orthogonal coordinate system on , and let be the signed distance to (say, in ). Expanding the perturbations in Fourier modes along the interface, the typical mode has the form
where is the Fourier wave vector, and . The shortest wavelength that can be represented on a grid with mesh size corresponds to . Hence, we obtain the estimate
for the maximum growth rate.
Clearly, is intimately related to the condition number for the discretized problem — see § 5. In fact, at leading order, the two numbers should be (roughly) proportional to each other — with a proportionality constant that depends on the details of the discretization. For the discretization used in this paper (described further below), , which leads to the rough estimate . On the other hand, the observed condition numbers vary between 5,000 and 10,000. Hence, the actual condition numbers are only slightly higher than for the ranges of grid sizes that we used (we did not explore the asymptotic limit ).
Equation (9) depends on the known inputs for the problem only. Namely: , , , and . Consequently does not depend on the solution . Hence, after solving for , we can use a discretization for that does not involve the interface: Whenever is discontinuous, we evaluate where the correction is needed, and transfer these values to the RHS.
When developing an algorithm for a linear Cauchy problem, such as the one in Eq. (9), the two key requirements are consistency and stability. In particular, when the solution depends on the “initial conditions” globally, stability (typically) imposes stringent constraints on the “time” step for any local (explicit) scheme. This would seem to suggest that, in order to solve Eq. (9), a “global” (involving the whole domain ) method will be needed. This, however, is not true: because we need to solve Eq. (9) for one “time” step only — i.e. within an distance from , stability is not relevant. Hence, consistency is enough, and a fully local scheme is possible. In the algorithm described in § 5 we found that, for (local) quadrangular patches, the Cauchy problem leads to a well behaved algorithm when the length of the interface contained in each patch is of the same order as the diagonal length of the patch. This result is in line with the calculation in remark 4: we want to keep the “wavelength” (along ) of the perturbations introduced by the discretization as long as possible. In particular, this should then minimize the condition number for the local problems — see remark 4.
5 A 4 Order Accurate Scheme in 2D.
5.1 Overview.
In this section we use the general ideas presented earlier to develop a specific example of a 4 order accurate scheme in 2D. Before proceeding with an indepth description of the scheme, we highlight a few key points:

We approximate using bicubic interpolations (bicubics), each valid in a small neighborhood of the interface. This guarantees local 4 order accuracy with only 12 interpolation parameters — see nave:10. Each corresponds to a point in the grid at which the standard discretization of Poisson’s equation involves a stencil that straddles the interface .

The domains are rectangular regions, each enclosing a portion of , and all the nodes where is needed to complete the discretization of the Poisson equation at the ()th stencil. Each is a subdomain of .

Starting from (b) and (c), we design a local solver that provides an approximation to inside each domain .

The interface is represented using the GradientAugmented Level Set approach — see nave:10. This guarantees a local 4 order representation of the interface, as required to keep the overall accuracy of the scheme.

In each , we solve the PDE in (9) in a least squares sense. Namely: First we define an appropriate positive quadratic integral quantity — Eq. (17) — for which the solution is a minimum (actually, zero). Next we substitute the bicubic approximation for the solution into , and discretize the integrals using Gaussian quadrature. Finally, we find the bicubic parameters by minimizing the discretized .
Solving the PDE in a least squares sense is crucial, since an algorithm is needed that can deal with the myriad ways in which the interface can be placed relative to the fixed rectangular grid used to discretize Poisson’s equation. This approach provides a scheme that (i) is robust with respect to the details of the interface geometry, (ii) has a formulation that is (essentially) dimension independent — there are no fundamental changes from 2D to 3D, and (iii) has a clear theoretical underpinning that allows extensions to higher orders, or to other discretizations of the Poisson equation.
5.2 Standard Stencil.
We use the standard 4 order accurate 9point discretization of Poisson’s equation^{2}^{2}2Notice that here we allow for the possibility of different grid spacings in each direction.:
(10) 
where is the second order 5point discretization of the Laplace operator:
(11) 
and
(12)  
(13) 
The terms and may be given analytically (if known), or computed using appropriate second order discretizations.
In the absence of discontinuities, Eq. (10) provides a compact 4 order accurate representation of Poisson’s equation. In the vicinity of the discontinuities at the interface , we define an appropriate domain , and compute the correction terms necessary to Eq. (10) — as described in detail next.
To understand how the correction terms affect the discretization, let us consider the situation depicted in figure 3. In this case, the node lies in while the nodes , , and are in . Hence, to be able to use Eq. (10), we need to compute , , and .
After having solved for where necessary (see § 5.3 and § 5.4), we modify Eq. (10) and write
(14) 
which differs from Eq. (10) by the terms on the RHS only. Here the are the CFM correction terms needed to complete the stencil across the discontinuity at . In the particular case illustrated in figure 3, we have
(15) 
Similar formulas apply for the other possible arrangements of the Poisson’s equation stencil relative to the interface .
5.3 Definition of .
There is some freedom on how to define . The basic requirements are

should be a rectangle.

the edges of should be parallel to the grid lines.

should contain all the nodes where is needed. For the example, in figure 3 we need to know , , and . Hence, in this case, should include the nodes , , and .

should contain a segment of , with a length that is as large as possible — i.e. comparable to the length of the diagonal of . This follows from the calculation in remark 4, which indicates that the wavelength of the perturbations (along ) introduced by the discretization should be as long as possible. This should then minimize the condition number for the local problem — see remark 4.
Requirements (i) and (ii) are needed for algorithmic convenience only, and do not arise from any particular argument in § 4. Thus, in principle, this convenience could be traded for improvements in other areas — for example, for better condition numbers for the local problems, or for additional flexibility in dealing with complex geometries. However, for simplicity, in this paper we enforce (i) and (ii). As explained earlier (see remark 5.1), we solve Eq. (9) in a least squares sense. Hence integrations over are required. It is thus useful to keep as simple as possible.
A discussion of various aspects regarding the proper definition of can be found in appendix B. For instance, the requirement in item (ii) is convenient only when an implicit representation of the interface is used. Furthermore, although the definition of presented here proved robust for all the applications of the 4 order accurate scheme (see § 6), there are specific geometrical arrangements of the interface for which (ii) results in extremely elongated . These elongated geometries can have negative effects on the accuracy of the scheme. We noticed this effect in the 2 order accurate version of the method described in appendix C. These issues are addressed by the various algorithms (of increasing complexity) presented in appendix B.
With the points above in mind, here (for simplicity) we define as the smallest rectangle that satisfies the requirements in (i), (ii), (iv), and (v) — then (iii) follows automatically. Hence can be constructed using the following three easy steps:

Find the coordinates () and () of the smallest rectangle satisfying condition (ii), which completely encloses the section of the interface contained by the region covered by the 9point stencil.

Find the coordinates () and () of the smallest rectangle satisfying condition (ii), which completely encloses all the nodes at which needs to be known.

Then is the smallest rectangle that encloses the two previous rectangles. Its edges are given by
(16a) (16b) (16c) (16d)
Figure 4 shows an example of defined using these specifications.
Notice that for each node next to the interface we construct a domain . When doing so, we allow the domains to overlap. For example, the domain shown in figure 4 is used to determine . It should be clear that (used to determine ), and (used to determine ), each will overlap with .
The consequence of these overlaps is that different computed values for at the same node can (in fact, will) happen — depending on which domain is used to solve the local Cauchy problem. However, because we solve for — within each — to 4 order accuracy, any differences that arise from this multiple definition of lie within the order of accuracy of the scheme. Since it is convenient to keep the computations local, the values of resulting from the domain , are used to evaluate the correction term .
While rare, cases where a single interface crosses the same stencil multiple times can occur. In § 6.3 we present such an example. A simple approach to deal with situations like this is as follows: First associate each node where the correction function is needed to a particular piece of interface crossing the stencil (say, the closest one). Then define one for each of the individual pieces of interface crossing the stencil.
For example, figure 5(a) depicts a situation where the stencil is crossed by two pieces of the same interface ( and ), with needed at the nodes , , , and . Then, first associate: (i) , , and to , and (ii) to . Second, define

is the smallest rectangle, parallel to the grid lines, that includes and the nodes , , and .

is the smallest rectangle, parallel to the grid lines, that includes and the node .
After the multiple are defined within a given stencil, the local Cauchy problem is solved for each separately. For example, in the case shown in figure 5(a), the solution for inside is done completely independent of the solution for inside . The decoupling between multiple crossings renders the CFM flexible and robust enough to handle complex geometries without any special algorithmic considerations.
When multiple distinct interfaces are involved, a single stencil can be crossed by different interfaces – e.g.: see § 6.5 and § 6.6. This situation is similar to the one described in remark 5.3, but with an additional complication: there may occur distinct domain regions that are not separated by an interface, but rather by a third (or more) regions between them. An example is shown in figure 5(b), where and are not part of the same interface. Here is the interface between and , while is the interface between and . There is no interface separating from , hence no jump conditions between these regions are provided. Nonetheless, is needed at .
Situations such as these can be easily handled by noticing that we can distinguish between primary (e.g. and ) and secondary correction functions, which can be written in terms of the primary functions (e.g. ) and need not be computed directly. Hence we can proceed exactly as in remark 5.3, except that we have to make sure that the intersections of the regions where the primary correction functions are computed include the nodes where the secondary correction functions are needed. For example, in the particular case in figure 5(b), we define

is the smallest rectangle, parallel to the grid lines, that includes and the nodes , , and .

is the smallest rectangle, parallel to the grid lines, that includes and the node .
5.4 Solution of the Local Cauchy Problem.
Since we use a 4 order accurate discretization of the Poisson problem, we need to find with 4 order errors (or better) to keep the overall accuracy of the scheme — see § 5.7. Hence we approximate using cubic Hermite splines (bicubic interpolants in 2D), which guarantees 4 order accuracy — see nave:10. Note also that, even though the example scheme developed here is for 2D, this representation can be easily extended to any number of dimensions.
Given a choice of basis functions,^{3}^{3}3The basis functions that we use for the bicubic interpolation can be found in appendix A. we solve the local Cauchy problem defined in Eq. (9) in a least squares sense, using a minimization procedure. Since we do not have boundary conditions, but interface conditions, we must resort to a minimization functional that is different from the standard one associated with the Poisson equation. Thus we impose the Cauchy interface conditions by using a penalization method. The functional to be minimized is then
(17) 
where is the penalization coefficient used to enforce the interface conditions, and is a characteristic length associated with — we used the shortest side length. Clearly is a quadratic functional whose minimum (zero) occurs at the solution to Eq. (9).
In order to compute in the domain , its bicubic representation is substituted into the formula above for , with the integrals approximated by Gaussian quadratures — in this paper we used six quadrature points for the 1D line integrals, and 36 points for the 2D area integrals. The resulting discrete problem is then minimized. Because the bicubic representation for involves 12 basis polynomials, the minimization problem produces a (selfadjoint) linear system.
We explored the option of enforcing the interface conditions using Lagrange multipliers. While this second approach yields good results, our experience shows that the penalization method is better.
The scaling using in Eq. (17) is so that all the three terms in the definition of behave in the same fashion as the size of changes with (), or when the computational grid is refined.^{4}^{4}4The scaling also follows from dimensional consistency. This follows because we expect that
Hence each of the three terms in Eq. (17) should be .
Once all the terms in Eq. (17) are guaranteed to scale the same way with the size of , the penalization coefficient should be selected so that the three terms have (roughly) the same size for the numerical solution (they will, of course, not vanish). In principle, could be determined from knowledge of the fourth order derivatives of the solution, which control the error in the numerical solution. This approach does not appear to be practical. A simpler method is based on the observation that should not depend on the grid size (at least to leading order, and we do not need better than this). Hence it can be determined empirically from a low resolution calculation. In the examples in this paper we found that produced good results.
A more general version of in Eq. (17) would involve different penalization coefficients for the two line integrals, as well as the possibility of these coefficients having a dependence on the position along of . These modifications could be useful in cases where the solution to the Poisson problem has large variations — e.g. a very irregular interface , or a complicated forcing .
5.5 Computational Cost.
We can now infer something about the cost of the present scheme. To start with, let us denote the number of nodes in the and directions by
(18) 
assuming a 1 by 1 computational square. Hence, the total number of degrees of freedom is . Furthermore, the number of nodes adjacent to the interface is , since the interface is a 1D entity.
The discretization of Poisson’s equation results in a linear system. Furthermore, the present method produces changes only on the RHS of the equations. Thus, the basic cost of inverting the system is unchanged from that of inverting the system resulting from a problem without an interface . Namely: it varies from to , depending on the solution method.
Let us now consider the computational cost added by the modifications to the RHS. As presented above, for each node adjacent to the interface, we must construct , compute the integrals that define the local linear system, and invert it. The cost associated with these tasks is constant: it does not vary from node to node, and it does not change with the size of the mesh. Consequently the resulting additional cost is a constant times the number of nodes adjacent to the interface. Hence it scales as . Because of the (relatively large) coefficient of proportionality, for small this additional cost can be comparable to the cost of inverting the Poisson problem. Obviously, this extra cost becomes less significant as increases.
5.6 Interface Representation.
As far as the CFM is concerned, the framework needed to solve the local Cauchy problems is entirely described above. However, there is an important issue that deserves attention: the representation of the interface. This question is independent of the CFM. Many approaches are possible, and the optimal choice is geometry dependent. The discussion below is meant to shed some light on this issue, and motivate the solution we have adopted.
In the present work, generally we proceed assuming that the interface is not known exactly – since this is what frequently happens. The only exceptions to this are the examples in § 6.5 and § 6.6, which involve two distinct (circular) interfaces touching at a point. In the generic setting, in addition to a proper representation of the interfaces, one needs to be able to identify the distinct interfaces, regions in between, contact points, as well as distinguish between a single interface crossing the same stencil multiple times and multiple distinct interfaces crossing one stencil. While the CFM algorithm is capable of dealing with these situations once they have been identified (e.g. see remarks 5.3 and 5.3), the development of an algorithm with the capability to detect such generic geometries is beyond the scope of this paper, and a (hard) problem in interface representation. For these reasons, in the examples in § 6.5 and § 6.6 we use an explicit exact representation of the interface.
To guarantee the accuracy of the solution for , the interface conditions must be applied with the appropriate accuracy — see § 5.7. Since these conditions are imposed on the interface , the location of must be known with the same order of accuracy desired for . In the particular case of the 4 order implementation of the CFM algorithm in this paper, this means 4 order accuracy. For this reason, we adopted the gradientaugmented level set (GALS) approach, as introduced in nave:10. This method allows a simple and completely local 4 order accurate representation of the interface, using Hermite cubics defined everywhere in the domain. The approach also allows the computation of normal vectors in a straightforward and accurate fashion.
We point out that the GALS method is not the only option for an implicit 4 order representation of the interface. For example, a regular level set method osher:88, combined with a highorder interpolation scheme, could be used as well. Here we adopted the GALS approach because of the algorithmic coherence that results from representing both the level set, and the correction functions, using the same bicubic polynomial base.
5.7 Error analysis.
A naive reading of the discretized system in Eq. (14) suggests that, in order to obtain a fourth order accurate solution , we need to compute the CFM correction terms with fourth order accuracy. Thus, from Eq. (15), it would follow that we need to know the correction function with sixth order accuracy! This is, however, not correct, as explained below.
Since we need to compute the correction function only at gridpoints an distance away from , it should be clear that errors in the are equivalent to errors in and of the same order. But errors in and produce errors of the same order in — see Eq. (1) and Eq. (2). Hence, if we desire a fourth order accurate solution , we need to compute the correction terms with fourth order accuracy only. This argument is confirmed by the convergence plots in figures 7, 9, and 11.
5.8 Computation of gradients.
Some applications require not only the solution to the Poisson problem, but also its gradient. Hence, in § 6, we include plots characterizing the behavior of the errors in the gradients of the solutions. A key question is then: how are these gradients computed?
To compute the gradients near the interface, the correction function can be used to extend the solution across the interface, so that a standard stencil can be used. However, for this to work, it is important to discretize the gradient operator using the same nodes that are part of the 9point stencil — so that the same correction functions obtained while solving the Poisson equation can be used. Hence we discretize the gradient operator with a procedure similar to the one used to obtain the 9point stencil. Specifically, we use the following 4 order accurate discretization:
(19)  
(20) 
where
(21)  
(22) 
and and are defined by (12) and (13), respectively. The terms and may be given analytically (if known), or computed using appropriate second order accurate discretizations.
This discretization is 4 order accurate. However, since the error in the correction function is (generally) not smooth, the resulting gradient will be less than 4 order accurate (worse case scenario is 3 order accurate) next to the interface.
6 Results.
6.1 General Comments.
In this section we present five examples of computations in 2D using the algorithm introduced in § 5. We solve the Poisson problem in the unit square for five different configurations. Each example is defined below in terms of the problem parameters (source term , and jump conditions across ), the representation of the interface(s) — either exact or using a level set function, and the exact solution (needed to evaluate the errors in the convergence plots). Notice that

As explained in § 5.6, in examples 1 through 3 we represent the interface(s) using the GALS method. Hence, the interface is defined by a level set function , with gradient — both of which are carried within the GALS framework nave:10.

Below the level set is described via an analytic formula. In examples 1 through 3 this formula is converted into the GALS representation for the level set, before it is fed into the code. Only this representation is used for the actual computations. This is done so as to test the code’s performance under generic conditions – where the interface would be known via a level set representation only.

Within the GALS framework we can, easily and accurately, compute the vectors normal to the interface — anywhere in the domain. Hence, it is convenient to write the jump in the normal derivative, , in terms of the jump in the gradient of dotted with the normal to the interface .
The last two examples involve touching circular interfaces and were devised to demonstrate the robustness of the CFM in the presence of interfaces that are very close together. In these last two examples, for the reasons discussed in § 5.6, we decided to use an exact representation of the circular interfaces.
6.2 Example 1.

Problem parameters:

Level set defining the interface: , where , , and .

Exact solution:
Figure 6 shows the numerical solution with a fine grid ( nodes). The discontinuity is captured very sharply, and it causes no oscillations in the solution. In addition, figure 7 shows the behavior of the error of the solution and its gradient in the and norms. As expected, the solution presents 4 order convergence as the grid is refined. Moreover, the gradient converges to 3 order in the norm and to 4 order in the norm, which is a reflection of the fact that the error in the solution is not smooth in a narrow region close to the interface only.
6.3 Example 2.

Problem parameters:

Level set defining the interface: , where , , , , , and .

Exact solution:
Figure 8 shows the numerical solution with a fine grid ( nodes). Once again, the overall quality of the solution is very satisfactory. Figure 9 shows the behavior of the error of the solution and its gradient in the and norms. Again, the solution converges to 4 order, while the gradient converges to 3 order in the norm and close to 4 order in the norm. However, unlike what happens in example 1, small wiggles are observed in the error plots. This behavior can be explained in terms of the construction of the sets — see § 5. The approach used to construct is highly dependent on the way in which the grid points are placed relative to the interface. Thus, as the grid is refined, the arrangement of the can vary quite a lot — specially for a “complicated” interface such as the one in this example. What this means is that, while one can guarantee that the correction function is obtained with 4 order precision, the proportionality coefficient is not constant — it may vary a little from grid to grid. This variation is responsible for the small oscillations observed in the convergence plot. Nevertheless, despite these oscillations, the overall convergence is clearly 4 order.
6.4 Example 3.

Problem parameters:

Level set defining the interface:
, where , , , and . 
Exact solution:
Figure 10 shows the numerical solution with a fine grid ( nodes). In this example, there are two circular interfaces in the solution domain. The two regions inside the circles make , while the remainder of the domain is . This example shows that the method is general enough to deal with multiple interfaces, keeping the same quality in the solution. Figure 11 shows that the solution converges to 4 order in both and norms, while the gradient converges to 3 order in the norm and close to 4 order in the norm.
6.5 Example 4.

Problem parameters:

Interface (exact representation):

Region 1: inside of the big circle.

Region 2: outer region.

Region 3: inside of the small circle.

Interface 1–2 (Big circle):

Interface 2–3 (Small circle):
