1 Introduction
The elliptic interface problem arises from various physical and biological problems such as fluid dynamics, heat conduction, electrostatics, where material interfaces or phase boundaries are involved. The interface can be static or dynamic. In these problems, material properties and sources may be discontinuous across the interface, which lead to discontinuous solution or flux across the interface.
In this paper, we look at the following elliptic interface problem:
(1) 
Here is an interface that separates a cubical computational domain into an inside region and an outside region .
is the outward unit normal vector of the interface. See Fig.
1. g is the boundary condition. , , : are given functions that might be discontinuous across . For some function and , we denote the limiting value approaching from different side of the interface as(2) 
and use the notation for the jump of across the interface
(3) 
The term is sometime called the flux. , : are given interface jump conditions. We are interested in solving for and accurately on a Cartesian grid with finite difference method.
In some applications, the dynamics of the interface depends on gradient of the solution, therefore it is important to have accurate solutions and gradient. For a molecule in ionic solutions, a solutesolvent interface separates the solvent region from the solute region. In this context, is the electric potential, is the charge density and is the dielectric coefficient, which is about 2 in solute region and 80 in solvent region shu_accurate_2014. In Variational ImplicitSolvent Model dzubiella_coupling_2006, the solutesolvent interface is the minimizer of a freeenergy functional which includes the electric potential . When we use the levelset method to numerically minimize such a VISM freeenergy functional zhou_variational_2014, we evolve the surface under the effective dielectric boundary force, which depends on the jump of the normal derivative of the electrostatic potential at the interface. The Stefan problem gibou_secondorderaccurate_2002 models the dynamic interface where one material is converted into the other. In this context, is the temperature and the is thermal conductivity, and the interface velocity depends on the jump in normal derivative of the temperature .
There are many finite difference methods that achieve second order accuracy in the solution and first order accuracy in the gradient, such as the Immersed Interface Method (IIM) leveque_immersed_1994, Ghost Fluid Method (GFM) liu_boundary_2000, Coupling Interface Method (CIM) chern_coupling_2007. For a review on different finite difference methods, we refer the readers to the introduction in chern_coupling_2007.
In the above finite difference methods, the standard central difference is used at grid points away from the interface, and different approaches are taken at grid points next to the interface. For dimensionsplitting approach, in each dimension, usually the principal secondorder derivative is approximated by values from both side of the interface and the jump in coordinate direction. Since only the jump along the normal direction is known, a coupling equation is needed to connect information from each dimension.
The CIM chern_coupling_2007 includes CIM1 (first order) and CIM2 (second order). They are derived from linear and quadratic approximation on both side of the interface in each dimension. The coupling equations are obtained by decomposing the jump in flux along coordinate direction into normal and tangential components, and the unknowns are the principal secondorder derivatives. CIM2 is second order accurate in solution and gradient. However, CIM2 requires two grid points in the same region on each side of the interface, and on each plane the mixed derivative can be approximated by grid data. The resulting stencil for CIM2 includes 1214 grid points in three dimensions. If these requirements are not met, which is common for complex interfaces in three dimensions, then CIM1 is used. This hybrid approach is second order accurate in solution but only first order accurate in the gradient.
The Improved Coupling Interface Method (ICIM) shu_accurate_2014 achieves second order accurate gradient by proposing 2 recipes to handle the exceptional points in CIM2: Recipe 1 is “shifting”: at points where
approximation of second derivative is not available, the second derivative can be approximated by central difference at a neighboring point. Recipe 2 is “flipping”: at some points, the signature of domain (inside or outside) can be flipped, so that the usual second order CIM2 discretization can be applied. Then the solution is second order accurate everywhere for the “flipped” surface, but at the flipped grid point the solution should be understood as extrapolation of the solution from the other side. Then by interpolation from the neighbors that are not flipped, we can recover the secondorder accurate solution and gradient at those flipped points.
We proposed a Compact Coupling Interface Method (CCIM). Our method combines elements of the CIM and Mayo’s approach for elliptic interface boundary value problemmayo_fast_1984. The jump in each coordinate direction is obtained by differentiating the jump condition in tangential directions. And the coupling equation involve the firstorder derivatives and principal secondorder derivatives. Our method is second order accurate in the solution and the gradient at the interface. The stencil is more compact, requiring a minimum of 10 grid points in three dimensions. The compactness makes the method applicable to more general situations and leads to more stable convergence results.
This paper is organized as follows. Section 2 outlines the derivation and algorithm of our CCIM method. In section 3, we show the convergence tests in three dimension on geometric surfaces and two complex protein surfaces. We also test our method on a moving surface driven by the jump in gradient at the interface. Section 4 is the conclusion.
2 Method
In dimensions, let and discretize the domain uniformly with mesh size , where is the number of subintervals on one edge of the region . Let be the multiindex with for . The grid points are denoted as with . Let , be the unit coordinate vectors. We also write . Here we use for the Laplacian of and for the Hessian matrix of . We assume that is symmetric. We use to denote the grid segment between and , and assume that the interface intersects with any grid segment at most once.
Let be a grid point at which we try to discretize the PDE. For notational simplicity, we drop the argument and the dependency on is implicit. We rewrite the PDE (1) at as
(4) 
If , and are in the same region in each coordinate direction, then we call an interior point, otherwise is called an interface point. At interior points, standard central differencing gives a local truncation error of in direction. Our goal is to construct finite difference schemes with local truncation error at interface points. The overall accuracy will still be second order since the interface points belongs to a lower dimensional set. Next, we derive a firstorder approximation for the term and .
2.1 Dimensionbydimension discretization
Along coordinate direction , if the interface does not intersect the grid segment , then by Taylor expansion
(5) 
Suppose the interface intersects the grid segment at . Let and . Suppose is located in . Denote the limit of as approaches from by , and the limit from the other side by . See Fig 2. By Taylor expansion at the interface , we can express and as
(6)  
Subtract the above two equations and write the right hand side in terms of jumps and quantities from :
(7) 
We can approximate components of and by
(8) 
(9) 
Thus, (7) can be written as
(10) 
By the given jump condition, . We denote the set of neighboring grid points of as and call the radius of our finite difference stencil. In three dimensions, contains 27 grid points and forms a cube. Suppose we can approximate the jump and in terms of , and , with and for some stencil radius . Then in each coordinate direction, for , we can write down two equations, (10) or (5), by considering the two grid segments for . In dimensions we have equations and unknowns: the firstorder derivatives and principle secondorder derivatives for . This leads to a system of linear equations of the following form:
(11) 
where , , is some linear function of neighboring values and is some constant. We call (11) coupling equation and coupling matrix. By inverting , we attain the finite difference formula to approximate and in terms of values.
In the next section, we describe how to remove the jumps in (10).
2.2 Approximation of the jump condition
Let be the unit normal vector at the interface, and , …, be unit tangent vectors. The tangent vectors can be obtained by projecting the unit coordinate vector onto the tangent plane. We can write
(12) 
In the coordinate direction, this gives
(13) 
We use the following trick repeatedly in our derivation to decouple the jump
(14) 
Taking , together with with , the jump condition can be rewritten as
(15) 
By (8), we get
(16) 
Notice that now we approximate the jump in the firstorder derivative at the interface in terms of the firstorder and secondorder derivatives at grid point.
To remove the jumps in the principal secondorder derivatives , we solve a system of linear equations. The first set of equations are obtained by differentiating the interface boundary condition in tangential directions. For and , we get equations
(17) 
Expand every term, and with the help of (14) and (15), we get
(18) 
By differentiating the jump in flux in tangential direction, we get another equations for ,
(19) 
After expansion,
(20)  
The final equation comes from the PDE
(21) 
After expansion,
(22)  
We arrive at a system of linear equations of where the variables are the jump in secondorder derivatives
(23) 
where is a matrix that only depends on the normal and tangent vectors, stands for some general linear function and stands for some constant. In two and three dimensions, through a direct calculation, the absolute value of the determinant of is 1.
As an example, in two dimensions, let and , and assume that is a piecewise constant function, then the system of linear equations is given by
(24)  
By Taylor’s expansion as in (5) (8), and (9), and components of and can all be approximated by and components of and at the grid point. Therefore, after substitution, with different and , (23) becomes
(25) 
If we could approximate the mixed derivatives , in terms of values and the jump in secondorder derivatives, then by inverting , the jump in secondorder derivatives can be approximated by values, firstorder derivatives, and principal secondorder derivatives, which are the terms used in the coupling equation. By back substitution, the mixed derivatives and thus the jump in the firstorder derivative can also be approximated by these terms. Next, we describe how to approximate the mixed derivatives.
2.3 Approximation of the mixed derivative
Depending on the position of the interface with respect to the grid points, different schemes are needed in order to approximate , . In Fig. 3 case 1, we use the usual central difference formula,
(28)  
In case 4, we can make use of the firstorder derivatives
(29)  
In case 5, we can make use of the firstorder and secondorder derivatives
(30)  
In Fig. 3 case 6, when there are not enough grid points on the same side, we can make use of information from the other side of the interface and the jump in mixed derivative
(31) 
By approximating the mixed derivatives , in terms of , , and , where for some kernel radius , we can write the right hand side of (23) in terms of , and , , thus eliminating the mixed derivatives. Then by solving the linear system (23), the jump in second derivative can be approximated in terms of , and , .
Finally, substituting , into (7), we get an equation involving only , , and for . In each dimension, we have 2 equations (by looking forward and backward). Therefore, with equations and unknowns (first order derivatives and principal second order derivatives), we can solve for and in terms of values.
Multiple methods to approximate the mixed derivatives might be available at the same grid point. We would like the method to be simple and accurate. For simplicity, we prefer methods that only make use of
value, as in case 1, 2 and 3. For accuracy, we provide a heuristic criterion: for homogeneous polynomial of degree three, the approximation on the right hand side is exact for some nearby points
or , We would like this point to be as close toas possible. In this respect, we prefer case 4 to case 3. Another consideration for accuracy is the condition number of the coupling equation. Solving a linear system with large condition number is prone to large numerical errors. Therefore, in cases where both case 3 and case 4 are available, we choose the method with a smaller estimated condition number
hager_condition_1984.Though we can construct surfaces for a specific grid size such that none of the above schemes works, for smooth surfaces we can refine the grid such that the above schemes suffice. In addition, we note that case 5 and case 6 can be removed by refining the grid, while case 4 cannot shu_accurate_2014.
2.4 Algorithm
We describe our method to obtain the coupling equations at an interface point in algorithmic order in Algorithm 1. Once we have the coupling equations (11), by inverting the coupling matrix, and , can be approximated by linear combinations of , .
To get more stable convergence results, at grid points where case 1 and 2 not available, but case 3 and 4 are available, we use the algorithm to obtain two systems of coupling equations, and choose the system with a smaller estimated condition number of the coupling matrix. The effect of this criterion is demonstrated in Section 3.1.
3 Numerical results
We test our method in three dimensions with different surfaces. The first set of tests contains six geometric surfaces that are used in shu_accurate_2014. And the second set of tests consists of two complex biomolecular surfaces. These two sets are compared with ICIM shu_accurate_2014 with the same setup. As tests in shu_accurate_2014 do not include term, the third set of tests are the same six geometric surfaces with term. The last test is a sphere expanding under a normal velocity given by the derivative of the solution in normal direction. Let be the exact solution of (1), and be the numerical solution. For tests with a static interface, we look at the maximum error of the solution at all grid points, denoted as , and the maximum error of the gradient at all intersection of interface and grid lines, denoted as . For the expanding sphere, we look at the maximum error and Root Mean Square Error of the radius at all intersections of the interface and grid line. All the tests are performed on a 2017 iMac with 3.5 GHz Intel Core i5 and 16GB memory. We use the AMG method implemented in the HYPRE library falgout_hypre_2002 to solve the sparse linear systems to a tolerance of .
3.1 Example 1
We test several geometric interfaces as in shu_accurate_2014. The surfaces are shown in Fig.4. Their level set functions are given below.

Eight balls: , where

Ellipsoid:

Peanut:

Donut:

Banana:

Popcorn: , where
The exact solution and the coefficient are given by
(32) 
and
(33) 
The source term and the jump conditions are calculated accordingly.
Fig 5 shows the convergence result of the six interfaces. The convergence of the solution at grid points is second order, and the convergence of the gradient at the interface is slightly below second order.
Next we demonstrate the effect of choosing the stencil for mixed derivatives based on estimated condition number of the coupling matrix. As mentioned in section 2.3, when both case 3 and case 4 are available to approximate the mixed derivatives, we choose the method that gives a smaller estimated condition number of the coupling matrix. We denote this scheme as CCIM. Alternatively, we can fix the order of preference for different methods. In scheme 1, we always prefer case 4 to case 3. Fig. 6 demonstrates the effect of this decision using the banana shape surface as an example. For different and for both schemes, Fig. 5(a) plots the maximum condition number of all the coupling matrices, and Fig. 5(a) plots the convergence results of these two schemes. For most of the values, CCIM and scheme 1 have roughly the same maximum error. We noticed that for , with scheme 1, at the interface point that produces the maximum error in the gradient, the coupling matrix has an exceptionally large condition number. By choosing the method with a smaller estimated condition number, we get a more stable convergence result in the gradient. If we prefer case 3 to case 4, then results are similar: at some grid points large condition number leads to large error, and CCIM gives more stable convergence results.
Though we can get a more stable convergence result by considering the condition number of the coupling matrices, there is a small jump for the banana interface at . A detailed analysis of the error reveals that it was caused by relatively large local truncation error when approximating . Fig. 7 shows the contour plot of the mixed derivative . Notice that change rapidly along the northeast direction. However, due to the alignment of the surface with the grid, at , our algorithm use the 4point stencil , , and to approximate and has a local truncation error 0.160723. If we use the three point stencil , , , the local truncation error would be 0.041853, and the coupling matrix does not have large condition number. With this surgical fix, the final error would be in line with the rest of the data points as shown in Fig 5(b)
at N = 110, marked as “Surgical fix”. This type of outliers happens rarely and does not affect the overall order of convergence. We apply this surgical fix only at this specific grid point to demonstrate a possible source of large error.
In summary, though the overall order of convergence is second order no matter which method is used to approximate the mixed derivatives, a relatively large error can be caused by large condition number of the coupling matrix, or a large local truncation error when approximating the mixed derivative. When different method to approximate the mixed derivative are available, ideally, we prefer the method that produces smaller local truncation error and smaller condition number of the coupling matrix. However these two goal might be incompatible sometimes. It’s time consuming to search through all available methods to approximate the mixed derivative and find the one that leads a small condition number of the coupling matrix. It’s also difficult to tell a priori which method gives smaller local truncation error. Therefore we try to find a middle ground by only considering the condition number when both case 3 and case 4 are available.
The resulting linear system for the PDE is sparse and asymmetric, and can be solved with any “blackbox” linear solvers. Fig. 8 shows the loglog plot for the number of iteration versus . We used BICGSTAB with ILU preconditioner and Algebraic Multigrid Method (AMG), both are implemented in the HYPRE library falgout_hypre_2002. The number of iterations grows linearly with for BICGSTAB and sublinearly for AMG. Though AMG has better scaling properties, for the range of in Fig 8, both solvers take approximately the same cpu time.
3.2 Example 2
Next we test our method on two complex molecular surfaces and compare CCIM with our implementation of ICIM shu_accurate_2014. The solvent accessible surface describes the interface between solute and solvent. Such interfaces are complex and important in applications. We construct the surfaces as in shu_accurate_2014: from the PDB file of 1D63 brown_crystal_1992 and MDM2 kussie_structure_1996, we use the PDB2PQR dolinsky_pdb2pqr_2004 software to assign charges and radii using the AMBER force field. The PQR files contain information of the positions and radii of the atoms. We scale the positions and radii such that the protein fit into the unit box. Then we construct the level set function of the interface as union of smoothed bumps:
(34) 
where
is a smoothed characteristic function
(35) 
The molecule 1D63 has 486 atoms and has a doublehelix shape, as shown in Fig. 9(a). MDM2 has 1448 atoms, and the surface has a deep pocket to which other proteins can bind, as shown in Fig 10(a). We also implement ICIM shu_accurate_2014 and compare the convergence results between CCIM and ICIM in Fig. 9 and Fig. 10.
, compared with our implementation of ICIM, the convergence results of CCIM is very robust even for complex interfaces. There is little fluctuation in the convergence results. In our ICIM implementation, the order of convergence exceeds second order because large errors at coarse grid points skew the fitting line to have a more negative slope. The results demonstrate the advantage of the compactness in our CCIM formulation when dealing with complex surfaces.
3.3 Example 3
We also test our problem with the same exact solution (32) and coefficients (33), but with an term, which is not handled in the ICIM formulation.
(36) 
As shown in Fig. 11, the convergence result is almost identical to those without the term. Therefore the convergence result mainly depends on the alignment of the surface with the grid line.
3.4 Example 4
In this example we look at the evolution of an interface driven by the jump in the normal derivative of the solution. Suppose the surface is evolved with normal velocity . The dynamics of the interface is given by the level set equation
(37) 
We use the forward Euler method for firstorder accurate time discretization, Godunov scheme for the Hamiltonian, and the Fast Marching Method osher_level_2003 to extend to the whole computational domain.
Consider the radially symmetric functions
(38) 
and
(39) 
The coefficient is the same as (33). The source term and the jump conditions are calculated accordingly.
If the surface is a sphere of radius , by symmetry, the normal velocity is uniform over the sphere and is given by
Comments
There are no comments yet.