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:
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
and use the notation for the jump of across the interface
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 solute-solvent 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 Implicit-Solvent Model dzubiella_coupling_2006, the solute-solvent interface is the minimizer of a free-energy functional which includes the electric potential . When we use the level-set method to numerically minimize such a VISM free-energy 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_second-order-accurate_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 dimension-splitting approach, in each dimension, usually the principal second-order 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 second-order 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 12-14 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 second-order 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 first-order derivatives and principal second-order 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.
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 multi-index 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
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 first-order approximation for the term and .
2.1 Dimension-by-dimension discretization
Along coordinate direction , if the interface does not intersect the grid segment , then by Taylor expansion
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
Subtract the above two equations and write the right hand side in terms of jumps and quantities from :
We can approximate components of and by
Thus, (7) can be written as
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 first-order derivatives and principle second-order derivatives for . This leads to a system of linear equations of the following form:
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
In the coordinate direction, this gives
We use the following trick repeatedly in our derivation to decouple the jump
Taking , together with with , the jump condition can be rewritten as
By (8), we get
Notice that now we approximate the jump in the first-order derivative at the interface in terms of the first-order and second-order derivatives at grid point.
To remove the jumps in the principal second-order 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
By differentiating the jump in flux in tangential direction, we get another equations for ,
The final equation comes from the PDE
We arrive at a system of linear equations of where the variables are the jump in second-order derivatives
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
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
If we could approximate the mixed derivatives , in terms of -values and the jump in second-order derivatives, then by inverting , the jump in second-order derivatives can be approximated by -values, first-order derivatives, and principal second-order derivatives, which are the terms used in the coupling equation. By back substitution, the mixed derivatives and thus the jump in the first-order 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,
We can also use biased differencing as in Fig. 3 case 2 and case 3
In case 4, we can make use of the first-order derivatives
In case 5, we can make use of the first-order and second-order derivatives
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
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 pointsor , We would like this point to be as close to
as 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 numberhager_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.
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
Popcorn: , where
The exact solution and the coefficient are given by
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 4-point 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 “black-box” linear solvers. Fig. 8 shows the log-log 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 sub-linearly 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:
is a smoothed characteristic function
The molecule 1D63 has 486 atoms and has a double-helix 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
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
We use the forward Euler method for first-order 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
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