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

Interface problems are common in many practical problems such as modeling flows in composite materials, oil reservoir simulations, nuclear waste disposal, and other flows in porous media [13]. For example, in groundwater or oil reservoir modelling the permeability of the porous medium can change drastically across the interfaces between various geological layers and this can significantly affect the transport process [30]. The coefficients of heterogeneous and anisotropic diffusion problems may also be highly oscillatory, and may contain a wide range of various spatial scales, so very fine meshes are required in any standard finite difference/element discretization in order to capture the small scale features. Thus, speed and storage are two important criteria in choosing suitable algorithms for solving such problems.

High jumps in the coefficient functions can cause severe singularities in the exact solutions of the equations [4, 17, 20, 2, 28, 29, 26, 19, 27, 11, 16, 18]. In general, the solutions of such problems have limited smoothness and the error analysis of their approximations by finite elements, [9], and finite differences, [31], for problems with weak solutions could be used. However, such error estimates show convergence rates that are lower than the observed in the computational practice for interface problems. Thus, an accurate tailored approximation and an error analysis which takes into account the specificity of such problem is an important and challenging task.

One possible approach to the resolution of this issue is provided by the so called immersed interface methods (IIM) proposed by LeVeque and Li (see [23, 21, 24, 25] and the references therein). The main idea behind this approach is to adjust the finite difference approximation of the differential operators in the vicinity of the interface using Taylor expansions, so that the approximation order remains similar to the order of the approximation in the regions where no singularities are present, thus avoiding the need of a local grid refinement. The explicit-jump immersed interface method, introduced in [37], is based on the same idea, however, instead of modification of the discrete operators, it modifies explicitly the right hand side of the problem, and derives a second order finite difference scheme for problems with discontinuous, piecewise constant coefficients. In fact this approach is quite similar to the famous immersed boundary method (IBM) of Peskin [40]. Exploiting the idea of the IIM, in [10] the authors construct a fourth order compact finite difference method for the Helmholtz equation with discontinuous coefficients across straight vertical line interfaces. [6] derives a compact finite difference method for piecewise smooth coefficients, so that the solution and its gradient can both achieve a second order of accuracy. [7] considers anisotropic elliptic interface problems whose coefficient matrix is symmetric semi-positive definite and derives a hybrid discretization involving finite elements away of the interfaces, and an immersed interface finite difference approximation near or at the interfaces. The error in the maximum norm is order $O(h^{2}\log\frac{1}{h})$. In addition to the treatment of interface problems, Taylor expansions can be used to derive high order compact finite difference schemes for regular elliptic problems. A family of fourth and sixth order compact finite difference methods for the three-dimensional Poisson equation are derived in [38]. [32] concludes that the highest order for a compact finite difference method for the two-dimensional Poisson’s equation on uniform grids is sixth.

Another source of singularity in the solution of elliptic problems is the presence of singularities in the source term. One possible regularization of Dirac delta functions is analyzed in [33], and [39] introduces the high order matched interface and boundary (MIB) method for elliptic equations with singular sources. Similarly, the finite difference discretization of such problems are considered by [34]. In [3] the idea of the IBM is combined with a Discontinuous Galerkin (DG) spatial discretization to derive a high-order method for elliptic problems with discontinuous coefficients and singular sources. In [14], the authors combine the idea of the IIM with a continuous finite element discretization to derive a high order finite element method for elliptic problems with jumps in the solution and its flux across smooth interfaces. Elliptic problems with point-located Dirac delta source terms are considered in [8]. A second order approximation to the singular source is combined with a second order finite difference approximation of the operator on Cartesian grids with hanging nodes, that allow for local refinements around the singular points. Another finite difference version of the immersed interface method is used to solve the heat diffusion with singular sources in [15].

Singular solutions, induced by discontinuous coefficients of singular sources can also be approximated using a continuous finite element approximation, by enriching the basis with singular functions located in the proper spatial locations, as considered in [4, 17, 20, 19, 11, 1, 12]. Alternatively, a posteriori error estimates can be used to devise grid refinement algorithms, as demonstrated for example in [29], where such estimates are provided in case of interface problems with discontinuous coefficients.

This short literature review is far from being complete, however, one can conclude that there are not many available discretizations of elliptic problems with singularities, that can achieve higher-than-second order of accuracy, and not rely on a local grid refinement. It is particularly challenging to derive a compact finite difference scheme that can achieve the highest possible order of accuracy i.e. a compact scheme of order six. The derivation of such a scheme is the main goal of the present paper. In particular, we consider Poisson’s equations with discontinuous or singular source terms, such that the solution and the normal component of gradient can be discontinuous across an interface curve.

To fix the ideas, let $\Omega=(l_{1},l_{2})\times(l_{3},l_{4})$ be a two-dimensional rectangular region. Let also $\psi$ be a smooth two-dimensional function and consider a smooth curve $\Gamma:=\{(x,y)\in\Omega:\psi(x,y)=0\}$, which partitions $\Omega$ into two subregions: $\Omega^{+}:=\{(x,y)\in\Omega\;:\;\psi(x,y)>0\}$ and $\Omega^{-}:=\{(x,y)\in\Omega\;:\;\psi(x,y)<0\}$. We define $f_{\pm}:=f\chi_{\Omega_{\pm}}$ and $u_{\pm}:=u\chi_{\Omega_{\pm}}$. The particular examples for $\psi(x,y)=x^{2}+y^{2}-2$ and $\psi(x,y)=y-\cos(x)$ are illustrated in Fig. 1.

We now state the Poisson interface problem with singular sources as follows:

which, if $g_{1}=0$, can be equivalently rewritten as

Here $\vec{n}$ is the unit normal vector of $\Gamma$ pointing towards $\Omega^{+}$, and for a point $(x_{0},y_{0})\in\Gamma$,

The conditions in (1.3) and (1.4) are called jump conditions for interface problems.

The Poisson interface problem in (1.1) with singular sources arise in many applications. In chemical reaction-diffusion processes, the solution $u$ represents the chemical concentration [15, 5]. In case of catalytic reactions the catalyst is usually distributed in a very thin layer over an interface $\Gamma$, and therefore the reaction can be considered as occurring on a $d-1$-dimensional manifold in a $d$-dimensional space. Such reactions result in a continuous chemical concentration $u$, but a discontinuous gradient $\nabla u$ across the interface $\Gamma$, i.e., $g_{1}(x,y)=0$ and $g(x,y)\neq 0$ for all $(x,y)\in\Gamma$.

Problems with discontinuous solutions and/or discontinuous fluxes appear also in multicomponent incompressible flows with or without interfacial tension. As discussed by [22], if surface tension is present at the interface the incompressibility constraint, applied to the momentum equation yields a pressure Poisson equation with a dipole source (the gradient of the delta function representing the interfacial tension alongside the fluid-fluid interface). Since such a source function is difficult to approximate, its effect can be modeled via interface jump conditions for the pressure and its gradient. The solution for the velocity is always continuous across the interface, however, If the viscosities of the fluids on both sides of the interface differ, its flux has a jump there. So, both the velocity and the pressure can be subject to elliptic problems with jumps of the solution or its flux across fluid-fluid or fluid-elastic-structure interfaces.

Similar jumps appear in the solution for the pressure $u$ and velocity $\vec{v}$, of the set of the Darcy’s and continuity equations:

in case of a discontinuous mobility $\frac{k}{\mu}$ and/or singular source $f$:

In this paper we consider the Poisson interface problem in (1.1) under the following assumptions:

• •

  $f$ is smooth in each of $\Omega^{+}$ and $\Omega^{-}$, and $f$ may be discontinuous across the interface curve $\Gamma$.

• •

  All functions $\psi$, $g$, $g_{0}$ and $g_{1}$ are smooth.

• •

  The exact solution $u$ is piecewise smooth in the sense that $u$ has uniformly continuous partial derivatives of (total) order up to seven in each of the subregions $\Omega^{+}$ and $\Omega^{-}$.

The remainder of the paper is organized as follows. In Section 2, we derive the sixth order compact finite difference scheme at regular and irregular points, and discuss their consistency in Theorem 2.2. and Theorem 2.4, correspondingly. Here the center of a stencil is called a regular point if it, together with all other eight points in the stencil are completely inside $\Omega^{+}$ or are completely outside $\Omega^{+}$. Otherwise, it is called an irregular point. We also give a simple proof for the maximum order of compact schemes at regular points. In Theorem 2.3 we provide an expression for the jump of certain derivatives of the solution, due to the interface conditions. In 2D there are $72$ different configurations for the stencil, depending on how the interface curve partitions the nine points in it. Up to a symmetry and rigid motion, all configurations at an irregular point can be classified into nine typical cases, see Figs. 3, 4 and 5 for a graphical representation of these configurations. In Section 3, we provide numerical experiments to check the convergence rate measured in $l_{2}$ and $l_{\infty}$ norms. We test the numerical examples in four cases: (1) the exact solution is known and $\Gamma$ does not intersect $\partial\Omega$; (2) the exact solution is known and $\Gamma$ intersects $\partial\Omega$; (3) the exact solution is unknown and $\Gamma$ does not intersect $\partial\Omega$; and (4) the exact solution is unknown and $\Gamma$ intersects $\partial\Omega$. In Section 4, we summarize the main contributions of this paper. Finally, in Appendix A we shall provide the detailed proof for Theorem 2.3, which plays a key role in our development of the compact stencils at irregular points in Section 2.

Recall that the problem domain is $\Omega=(l_{1},l_{2})\times(l_{3},l_{4})$. Because we shall use uniform Cartesian meshes, we require that the longer side of $\Omega$ is a multiple of the shorter side of $\Omega$. Without loss of generality, we can assume $l_{4}-l_{3}=N_{0}(l_{2}-l_{1})$ for some positive integer $N_{0}$. For any positive integer $N_{1}\in\mathcal{N}$, we define $N_{2}:=N_{0}N_{1}$ and then the grid size is $h:=(l_{2}-l_{1})/N_{1}=(l_{4}-l_{3})/N_{2}$.

Let $x_{i}=l_{1}+ih$ and $y_{j}=l_{3}+jh$ for $i=1,\ldots,N_{1}-1$ and $j=1,\ldots,N_{2}-1$. Because in this paper we are only interested in compact finite difference schemes on uniform Cartesian grids, for a compact stencil centered at the center point $(x_{i},y_{j})$, the compact stencil involves nine points $(x_{i}+kh,y_{j}+lh)$ for $k,l\in\{-1,0,1\}$. It is convenient to use a level set function $\psi$, which is a two-dimensional smooth function, to describe a given smooth interface curve $\Gamma$ through

Then the interface curve $\Gamma$ splits the problem domain $\Omega$ into two subregions: $\Omega^{+}:=\{(x,y)\in\Omega\;:\;\psi(x,y)>0\}$ and $\Omega^{-}:=\{(x,y)\in\Omega\;:\;\psi(x,y)<0\}$. Now the interface curve $\Gamma$ splits these nine points into two groups depending on whether these points lie inside $\Omega^{+}$ or $\Omega^{-}$. If a grid point lies on the curve $\Gamma$, then the grid point lies on the boundaries of both $\Omega^{+}$ and $\Omega^{-}$. For simplicity we may assume that the grid point belongs to $\Omega^{+}$ and we can use the interface conditions to handle such a grid point. Therefore, we naturally define

and

That is, the interface curve $\Gamma$ splits the nine points in a compact stencil into two disjoint sets $\{(x_{i+k},y_{j+\ell})\;:\;(k,\ell)\in d_{i,j}^{+}\}\subseteq\Omega^{+}$ and $\{(x_{i+k},y_{j+\ell})\;:\;(k,\ell)\in d_{i,j}^{-}\}\subseteq\Omega^{-}$. We say that a grid/center point $(x_{i},y_{j})$ is a regular point if $d_{i,j}^{+}=\emptyset$ or $d_{i,j}^{-}=\emptyset$. That is, the center point $(x_{i},y_{j})$ of a stencil is regular if all its nine points are completely inside $\Omega^{+}$ (hence $d_{i,j}^{-}=\emptyset$) or inside $\Omega^{-}$ (i.e., $d_{i,j}^{+}=\emptyset$). See Fig. 2 for an example of regular points. Otherwise, the center point $(x_{i},y_{j})$ of a stencil is called an irregular point if $d_{i,j}^{+}\neq\emptyset$ and $d_{i,j}^{-}\neq\emptyset$. That is, the interface curve $\Gamma$ splits the nine points into two disjoint nonempty sets. As explained before, up to a symmetry and rigid motion, all the compact stencils at an irregular point can be classified into nine typical cases, see Figs. 3, 4 and 5 for these nine typical cases.

Let $M$ be a positive integer standing for a given desired accuracy order. Before we discuss the compact schemes at a regular or an irregular point $(x_{i},y_{j})$, let us introduce some notations and then outline the main ideas for finding the stencils at a general grid point $(x_{i},y_{j})$, achieving a given accuracy order $M$.

We first pick up and fix a base point $(x_{i}^{*},y_{j}^{*})$ inside the open square $(x_{i}-h,x_{i}+h)\times(y_{j}-h,y_{j}+h)$. That is, we can write

For the sake of brevity, we shall use the following notions:

which are just their $(m,n)$th partial derivatives at the base point $(x_{i}^{*},y_{j}^{*})$. Define $\mathbb{N}_{0}:=\mathbb{N}\cup\{0\}$, the set of all nonnegative integers. For a nonnegative integer $K\in\mathbb{N}_{0}$, we define

For a smooth function $u$, its values $u(x+x_{i}^{*},y+y_{j}^{*})$ for small $x,y$ can be well approximated through its Taylor polynomial below:

In other words, in a neighborhood of the base point $(x_{i}^{*},y_{j}^{*})$, the function $u$ is well approximated and completely determined by the partial derivatives of $u$ of total degree less than $M+2$ at the base point $(x_{i}^{*},y_{j}^{*})$, i.e., by the unknown quantities $u^{(m,n)},(m,n)\in\Lambda_{M+1}$. We can similarly approximate $f(x+x_{i}^{*},y+y_{j}^{*})$ for small $x,y$. For $x\in\mathbb{R}$, the floor function $\lfloor x\rfloor$ is defined to be the largest integer less than or equal to $x$. For an integer $m$, we define

That is, $\operatorname{odd}(m)=m-2\lfloor m/2\rfloor$ and $\lfloor m/2\rfloor=\frac{m-\operatorname{odd}(m)}{2}$. Because the function $u$ is a solution to the partial differential equation in (1.1), we shall see that all the quantities $u^{(m,n)},(m,n)\in\Lambda_{M+1}$ are not independent of each other. In fact, we have:

For the convenience of the reader, see Fig. 6 for an illustration of the quantities $u^{(m,n)},(m,n)\in\Lambda_{M+1}^{1}$ and $(m,n)\in\Lambda_{M+1}^{2}$ in Lemma 2.1 with $M=6$.

For $M=6$, the identities in (2.5) of Lemma 2.1 for $u^{(m,n)},(m,n)\in\Lambda_{7}^{2}$ can be explicitly given by

Note that the cardinality of $\Lambda_{M+1}$ equals the sum of the cardinalities of $\Lambda_{M+1}^{1}$ and $\Lambda_{M-1}$. The identities in (2.5) of Lemma 2.1 show that every $u^{(m,n)},(m,n)\in\Lambda_{M+1}$ can be written as a linear combination of the quantities $u^{(m,n)},(m,n)\in\Lambda_{M+1}^{1}$ and $f^{(m,n)},(m,n)\in\Lambda_{M-1}$. Conversely, by (2.7), every $f^{(m,n)},(m,n)\in\Lambda_{M-1}$ and every $u^{(m,n)},(m,n)\in\Lambda_{M+1}^{1}$ can be trivially written as linear combinations of $u^{(m,n)}\in\Lambda_{M+1}$. Because the source term $f$ is known, this can reduce the number of constraints on $u^{(m,n)},(m,n)\in\Lambda_{M+1}$ for the function $u$ satisfying (1.1). Now using (2.5), we can rewrite the approximation of $u(x+x_{i}^{*},y+y_{j}^{*})$ in (2.4) as follows:

where $G_{m,n}$ and $H_{m,n}$ are polynomials uniquely determined by the identities in (2.5). Explicitly,

and

From (2.5) we observe that $G_{m,n}$ is a homogeneous polynomial of total degree $m+n$ for all $(m,n)\in\Lambda_{M+1}^{1}$, while $H_{m,n}$ is a homogeneous polynomial of total degree $m+n+2$ for all $(m,n)\in\Lambda_{M-1}$. For $M=6$, all the polynomials $G_{m,n},(m,n)\in\Lambda_{7}^{1}$ and $H_{m,n},(m,n)\in\Lambda_{5}$ are explicitly given by

and

Hence, by (2.4), the solution $u$ to (1.1) near the base point $(x_{i}^{*},y_{j}^{*})$ can be approximated by

for $x,y\in(-2h,2h)$. We shall use the above identity in (2.13) for finding compact stencils achieving a desired accuracy order $M$.