A Compact Coupling Interface Method with Accurate Gradient Approximation for Elliptic Interface Problems

by   Zirui Zhang, et al.
University of California, San Diego

Elliptic interface boundary value problems play a major role in numerous applications involving heat, fluids, materials, and proteins, to name a few. As an example, in implicit variational solvation, for the construction of biomolecular shapes, the electrostatic contributions satisfy the Poisson-Boltzmann equation with discontinuous dielectric constants across the interface. When interface motions are involved, one often needs not only accurate solution values, but accurate derivatives as well, such as the normal derivatives at the interface. We introduce here the Compact Coupling Interface Method (CCIM), a finite difference method for the elliptic interface problem with interfacial jump conditions. The CCIM can calculate solution values and their derivatives up to second-order accuracy in arbitrary ambient space dimensions. It combines elements of Chern and Shu's Coupling Interface Method and Mayo's approach for elliptic interface boundary value problems, leading to more compact finite difference stencils that are applicable to more general situations. Numerical results on a variety of geometric interfacial shapes and on complex protein molecules in three dimensions support the efficacy of our approach and reveal advantages in accuracy and robustness.



There are no comments yet.


page 10


A High Order Compact Finite Difference Scheme for Elliptic Interface Problems with Discontinuous and High-Contrast Coefficients

The elliptic interface problems with discontinuous and high-contrast coe...

Hybrid Finite Difference Schemes for Elliptic Interface Problems with Discontinuous and High-Contrast Variable Coefficients

For elliptic interface problems with discontinuous coefficients, the max...

A diffuse interface box method for elliptic problems

We introduce a diffuse interface box method (DIBM) for the numerical app...

On jump relations of anisotropic elliptic interface problems

Almost all materials are anisotropic. In this paper, interface relations...

Sixth Order Compact Finite Difference Scheme for Poisson Interface Problem with Singular Sources

Let Γ be a smooth curve inside a two-dimensional rectangular region Ω. I...

A Novel HOC-Immersed Interface Approach For Elliptic Problems

We present a new higher-order accurate finite difference explicit jump I...

Projection-based resolved interface mixed-dimension method for embedded tubular network systems

We present a flexible discretization technique for computational models ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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:


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.

Figure 1: schematic for elliptic interface problem

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.

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 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

Figure 2: and are approximated by Taylor expansion at interface

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


Expand every term, and with the help of (14) and (15), we get


By differentiating the jump in flux in tangential direction, we get another equations for ,


After expansion,


The final equation comes from the PDE


After expansion,


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

Figure 3: Approximation of mixed derivative at . The circles and disks are grid points in and . The squares are grid points that are used to approximate the mixed derivative. Case 1 is usual central differencing. Case 2 and 3 are biased differencing. Case 4 is coupled with the second derivative. Case 5 is coupled with the first and second derivative. Case 6 make use of the mixed derivative on the other side and the jump condition

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 points

or , 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 number


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 , .

1:  for  do
2:     for  do
3:        if the interface intersects at  then
4:           for ,  do
5:              Approximate the mixed derivative in terms of , and , as in section 2.3. If Fig.3 case 6 (using information from the other side) is used, then and the coefficient of is non-zero. Otherwise and the term vanishes.
6:           end for
7:           Differentiate the jump conditions and obtain a system of equations of the jump in second derivatives (23). substitute the approximation of the mixed derivatives. Invert the matrix , then the jumps in second derivatives , , are approximated in terms of , , , , .
8:           By back substitution, the mixed derivative (if case 6 exists) and thus (16) can be expressed in terms of , , , , .
9:           Substitute the expression for and into (10). After rearrangement, this gives one row of the coupling equations (11).
10:        else
11:           Taylor’s expansion (5) gives one row of the coupling equations (11).
12:        end if
13:     end for
14:  end for
Algorithm 1 Coupling equation at interface point

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




The source term and the jump conditions are calculated accordingly.

Figure 4: The six interfaces (a) eight balls; (b) ellipsoid; (c) peanut; (d) donut; (e) banana; (f) popcorn

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.

Figure 5: The log-log plot of the error versus for the six surfaces. In each figure, N ranges from 50 to 140 with the increment . Circles are the maximum errors of the solution . Diamonds are the maximum errors of the gradient at interface .

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.

(a) maximum condition number
(b) Convergence
Figure 6: Comparison of maximum condition numbers and convergence results with banana surface between scheme 1 and CCIM. In scheme 1, case 4 is preferred to case 3. In CCIM, case 3 and case 4 are chosen based on estimated condition numbers of the coupling matrices. (a) The maximum condition number of coupling matrices for different methods. (b) Log-log plot of the maximum errors in solution and gradient. The jump in error at for CCIM can be reduced by a surgical fix that choose the stencil with small local truncation error.

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.

Figure 7: Contour of at the grid point with maximum error in gradient with the banana interface at . The disks and the circles are grid points outside and inside the surface. Grid points marked with squares (Fig. 3 case 3) are used to approximate at due to its simplicity, while the three points stencil (Fig. 3 case 4)with , , has a smaller local truncation 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.

Figure 8: The log-log plot of the number of iterations versus for the six surfaces.

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.

Figure 9: Convergence result for 1D63 interface with and . (a) The smooth surface of 1D63. (b) log-log plot of error by CCIM. (b) log-log plot of error by ICIM. ranges from 100 to 340 with the increment
Figure 10: Convergence result for MDM2 interface with and . (a) The smooth surface of MDM2. (b) log-log plot of error by CCIM. (b) log-log plot of error by ICIM. ranges from 100 to 340 with the increment

As shown in Fig. 9 and  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.


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.

Figure 11: The log-log plot error with term versus for the six surfaces. In each figure, N ranges from 50 to 140 with the increment .

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