Solving Elliptic Interface Problems with Jump Conditions on Cartesian Grids

TL;DR

This paper introduces a straightforward numerical method for elliptic equations with discontinuities across irregular interfaces, achieving high accuracy and stable condition numbers on Cartesian grids in 2D and 3D.

Contribution

The paper presents a simple, second-order accurate algorithm for elliptic interface problems with jump conditions on Cartesian grids, maintaining bounded condition numbers regardless of diffusion coefficient ratios.

Findings

Achieves second-order accuracy in solutions

Attains first-order accuracy in gradients in the $L^ abla$-norm

Condition number remains bounded regardless of diffusion coefficient ratio

Abstract

We present a simple numerical algorithm for solving elliptic equations where the diffusion coefficient, the source term, the solution and its flux are discontinuous across an irregular interface. The algorithm produces second-order accurate solutions and first-order accurate gradients in the $L^\infty$-norm on Cartesian grids. The condition number is bounded, regardless of the ratio of the diffusion constant and scales like that of the standard 5-point stencil approximation on a rectangular grid with no interface. Numerical examples are given in two and three spatial dimensions.

Figures

Tables

Particle number	$x_c$	$y_c$	$r_0$	$\delta$	$m$	${\theta }_0$	Charge
1	0.294320	-0.731980	0.000948	0.177428	5	3.779031	+
2	0.292603	-0.722570	0.000196	0.120466	2	0.575194	+
3	0.296621	-0.730048	0.000553	0.195375	3	5.983689	+
4	0.301392	-0.723932	0.000374	0.161709	3	1.178485	+
5	0.293268	-0.730192	0.000798	0.193548	3	6.112317	-
6	0.289541	-0.732616	0.000755	0.088156	4	2.266495	-
7	0.290527	-0.736414	0.000341	0.140005	3	6.140219	-
8	0.295213	-0.726113	0.000792	0.121836	5	0.884964	-
9	0.300884	-0.733236	0.000857	0.078646	4	1.588758	-
10	0.935243	0.156877	0.000462	0.026599	4	1.152778	+
11	0.936367	0.165333	0.000946	0.129853	3	2.585702	+
12	0.937713	0.161907	0.000114	0.111576	3	3.069164	+
13	0.927721	0.160662	0.000653	0.059185	4	4.986538	-
14	0.926407	0.164860	0.000243	0.009826	2	3.613796	-
15	0.939303	0.170991	0.000443	0.141560	3	2.988958	+
16	0.926755	0.162948	0.000482	0.025398	5	3.380732	+
17	0.936908	0.160965	0.000713	0.179751	3	0.252419	+
18	0.934454	0.162162	0.000448	0.011111	6	2.899990	+
19	0.160228	-0.743358	0.000329	0.052804	6	4.507911	+
20	0.164293	-0.747939	0.000792	0.015578	3	3.753404	+
21	0.173756	-0.744021	0.000264	0.186352	4	6.101247	+
22	0.168733	-0.761421	0.000999	0.001466	5	4.162976	-
23	0.170695	-0.747981	0.000389	0.025428	4	0.952773	-
24	0.160737	-0.760843	0.000209	0.199829	6	0.814897	-
25	0.725447	0.889409	0.000966	0.064583	3	5.719127	-
26	0.740474	0.893693	0.000453	0.004282	3	5.779404	-
27	0.736547	0.884784	0.000708	0.036598	4	0.871440	-
28	0.723810	0.875899	0.000367	0.044844	3	0.309913	-
29	0.742820	0.874970	0.000306	0.191154	4	1.616747	-
30	0.728138	0.878376	0.000857	0.026387	4	5.913726	-
31	0.728254	0.880038	0.000451	0.069948	6	2.499319	+
32	0.052907	-0.344586	0.000660	0.184134	4	2.412480	-
33	0.034801	-0.356601	0.000834	0.079140	4	4.794653	-
34	0.048034	-0.350902	0.000961	0.175580	4	1.000997	-
35	0.038434	-0.350375	0.000546	0.119999	5	5.330543	-
36	0.037001	-0.361473	0.000859	0.104759	6	2.104655	-
37	0.048676	-0.344655	0.000590	0.018488	3	3.979605	-
38	-0.182697	0.433937	0.000642	0.191264	5	0.196232	-
39	-0.180690	0.449481	0.000499	0.107538	4	3.958288	+
40	-0.171622	0.442325	0.000714	0.162207	3	4.889946	-
41	-0.177232	0.435290	0.000706	0.108419	6	2.765306	+
42	-0.182553	0.446933	0.000927	0.139933	5	0.682710	-
43	-0.169332	0.440621	0.000652	0.197255	6	0.524252	-
44	0.082096	0.665721	0.000309	0.018026	6	3.925875	+
45	0.070903	0.669204	0.000229	0.189876	3	1.973639	+
46	0.086667	0.652859	0.000256	0.009732	4	5.337312	-
47	0.089332	0.668903	0.000778	0.189929	4	0.506358	+
48	-0.416502	-0.873554	0.000168	0.038622	2	4.772787	+
49	-0.435866	-0.873537	0.000167	0.040327	6	4.357852	+
50	-0.428412	-0.875426	0.000329	0.067637	4	4.422504	-
51	0.367410	-0.528637	0.000239	0.004636	6	1.018504	-
52	0.364188	-0.533808	0.000271	0.072883	3	2.058938	-
53	0.357943	-0.522965	0.000682	0.097692	3	4.725618	+
54	0.354170	-0.530541	0.000809	0.035313	5	0.726723	-
55	0.355478	-0.527060	0.000107	0.166186	4	5.450653	+
56	0.356304	-0.525395	0.000985	0.004988	2	3.471888	+
57	0.368282	-0.530008	0.000107	0.133467	5	5.475823	-
58	0.370185	-0.524812	0.000830	0.145549	5	6.126808	-
59	0.352045	-0.519527	0.000198	0.011507	4	3.886958	+
60	-0.187080	0.207027	0.000928	0.121431	2	1.730282	+
61	-0.177875	0.213721	0.000126	0.001767	3	0.076936	+
62	-0.177947	0.216082	0.000647	0.130215	4	6.002824	-
63	-0.193546	0.211842	0.000749	0.088286	6	4.366932	-
64	-0.190933	0.205428	0.000170	0.122096	3	0.951230	-
65	-0.184008	0.210547	0.000752	0.030592	3	3.986605	-
66	-0.192313	0.214133	0.000925	0.029170	4	1.284015	-
67	-0.190761	0.223399	0.000728	0.195173	4	5.881409	-

Equations

$$k^{\pm }{u}^{\pm }-\nabla \cdot \left({\mu }^{\pm }\nabla u^{\pm }\right)$$

$$\left[u\right]$$

$$\left[\mu {\partial }_{\mathbf{n}}u\right]$$

$$\Delta x=\frac{x_{\max }-x_{\min }}{n_x-1},\Delta y=\frac{y_{\max }-y_{\min }}{n_y-1}$$

$$\underset{\text{Linear term}}{\underbrace{\sum _{s=+,-}{\int }_{{\Omega }^s\cap {\mathcal{V}}_{i,j}}{k}^s{u}^{s}d\Omega }}-\underset{\text{Flux between finite volumes}}{\underbrace{\sum _{s=+,-}{\int }_{{\Omega }^s\cap \partial {\mathcal{V}}_{i,j}}{\mu }^s{\partial }_{{\mathbf{n}}^s}{u}^{s}d\Gamma }}=\underset{\text{Volumetric generation}}{\underbrace{\sum _{s=+,-}{\int }_{{\Omega }^s\cap {\mathcal{V}}_{i,j}}{f}^{s}d\Omega }}+\underset{\text{Surface generation}}{\underbrace{{\int }_{\Gamma \cap {\mathcal{V}}_{i,j}}\left[\mu {\partial }_{\mathbf{n}}u\right]d\Gamma ,}}$$

$$\displaystyle\begin{aligned} \sum_{s=+,-}k_{i,j}^{s}u^{s}_{i,j}|\mathcal{V}_{i,j}^{s}|-\sum_{s=+,-}\Bigg{(}\mu^{s}_{i-{\frac{1}{2}},j}A_{i-{\frac{1}{2}},j}^{s}\frac{u^{s}_{i-1,j}-u^{s}_{i,j}}{\Delta x}&+\mu^{s}_{i+{\frac{1}{2}},j}A_{i+{\frac{1}{2}},j}^{s}\frac{u^{s}_{i+1,j}-u^{s}_{i,j}}{\Delta x}+\\ \mu^{s}_{i,j-{\frac{1}{2}}}A_{i,j-{\frac{1}{2}}}^{s}\frac{u^{s}_{i,j-1}-u^{s}_{i,j}}{\Delta y}&+\mu^{s}_{i,j+{\frac{1}{2}}}A_{i,j+{\frac{1}{2}}}^{s}\frac{u^{s}_{i,j+1}-u^{s}_{i,j}}{\Delta y}\Bigg{)}\\ =\sum_{s=+,-}f_{i,j}^{s}|\mathcal{V}_{i,j}^{s}|&+\int_{\Gamma\cap\mathcal{V}_{i,j}}\beta\,d\Gamma+\mathcal{O}\left(h^{\mathcal{D}}\right),\end{aligned}$$

$$u_{i,j}=\{\begin{array}{l}{u}_{i,j}^{+},{\mathbf{r}}_{i,j}\in {\Omega }^{+},\\ u_{i,j}^{-},{\mathbf{r}}_{i,j}\in {\Omega }^{-}.\end{array}$$

$$u_{i,j}^{\pm }=u^{\pm }({\mathbf{r}}_{i,j}^{\text{pr}})+{\delta }_{i,j}{\partial }_{\mathbf{n}}{u}^{\pm }({\mathbf{r}}_{i,j}^{\text{pr}})+\mathcal{O}\left(h^2\right),$$

$$\mathbf{n}({\mathbf{r}}_{i,j}^{\text{pr}})$$

$${\mathbf{r}}_{i,j}^{\text{pr}}$$

$${\delta }_{i,j}$$

$$u_{i,j}^{+}-u_{i,j}^{-}=\alpha ({\mathbf{r}}_{i,j}^{\text{pr}})+{\delta }_{i,j}\left({\partial }_{\mathbf{n}}{u}^{+}({\mathbf{r}}_{i,j}^{\text{pr}})-{\partial }_{\mathbf{n}}{u}^{-}({\mathbf{r}}_{i,j}^{\text{pr}})\right)+\mathcal{O}\left(h^2\right).$$

$$u_{i,j}^{+}-u_{i,j}^{-}=\{\begin{array}{l}\alpha ({\mathbf{r}}_{i,j}^{\text{pr}})+{\delta }_{i,j}\frac{\beta ({\mathbf{r}}_{i,j}^{\text{pr}})}{{\mu }^{+}({\mathbf{r}}_{i,j}^{\text{pr}})}-{\delta }_{i,j}\frac{{\mu }^{+}({\mathbf{r}}_{i,j}^{\text{pr}})-{\mu }^{-}({\mathbf{r}}_{i,j}^{\text{pr}})}{{\mu }^{+}({\mathbf{r}}_{i,j}^{\text{pr}})}{\partial }_{\mathbf{n}}{u}^{-}({\mathbf{r}}_{i,j}^{\text{pr}})\\ \alpha ({\mathbf{r}}_{i,j}^{\text{pr}})+{\delta }_{i,j}\frac{\beta ({\mathbf{r}}_{i,j}^{\text{pr}})}{{\mu }^{-}({\mathbf{r}}_{i,j}^{\text{pr}})}-{\delta }_{i,j}\frac{{\mu }^{+}({\mathbf{r}}_{i,j}^{\text{pr}})-{\mu }^{-}({\mathbf{r}}_{i,j}^{\text{pr}})}{{\mu }^{-}({\mathbf{r}}_{i,j}^{\text{pr}})}{\partial }_{\mathbf{n}}{u}^{+}({\mathbf{r}}_{i,j}^{\text{pr}})\end{array}+\mathcal{O}\left(h^2\right).$$

$${\partial }_{\mathbf{n}}{u}^{\pm }({\mathbf{r}}_{i,j}^{\text{pr}})=\mathbf{n}({\mathbf{r}}_{i,j}^{\text{pr}})\nabla u_I^{\pm }({\mathbf{r}}_{i,j}^{\text{pr}}).$$

$$u_I^{\pm }(\mathbf{r})=u_{i,j}^{\pm }+{\left(\mathbf{r}-{\mathbf{r}}_{i,j}\right)}^T{\left(\nabla u^{\pm }\right)}_{i,j}+\mathcal{O}\left(h^2\right),$$

$$u_{i+p,j+q}=u_{i,j}^{\pm }+{\left({\mathbf{r}}_{i+p,j+q}-{\mathbf{r}}_{i,j}\right)}^T{\left(\nabla u^{\pm }\right)}_{i,j},(p,q)\in N_{i,j}^{\pm }.$$

$$N_{i,j}^{\pm }$$

$${\underline{\underline{\mathbf{X}}}}_{i,j}{\underline{\underline{\mathbf{W}}}}_{i,j}^{\pm }{\left(\nabla u^{\pm }\right)}_{i,j}={\underline{\underline{\mathbf{W}}}}_{i,j}^{\pm }\left(\begin{array}{c}{u}_{i-1,j-1}-u_{i,j}^{\pm }\\ u_{i,j-1}-u_{i,j}^{\pm }\\ \dots \\ u_{i+1,j+1}-u_{i,j}^{\pm }\end{array}\right),$$

$${\left(\nabla u^{\pm }\right)}_{i,j}={\underline{\underline{\mathbf{D}}}}_{i,j}^{\pm }\left(\begin{array}{c}{u}_{i-1,j-1}-u_{i,j}^{\pm }\\ u_{i,j-1}-u_{i,j}^{\pm }\\ \dots \\ u_{i+1,j+1}-u_{i,j}^{\pm }\end{array}\right),{\underline{\underline{\mathbf{D}}}}_{i,j}^{\pm }={\left({\underline{\underline{\mathbf{X}}}}_{i,j}^T{\underline{\underline{\mathbf{W}}}}_{i,j}^{\pm }{\underline{\underline{\mathbf{X}}}}_{i,j}\right)}^{-1}{\left({\underline{\underline{\mathbf{W}}}}_{i,j}^{\pm }{\underline{\underline{\mathbf{X}}}}_{i,j}\right)}^T,$$

$${\underline{\underline{\mathbf{X}}}}_{i,j}=\left(\begin{array}{c}{\left({\mathbf{r}}_{i-1,j-1}-{\mathbf{r}}_{i,j}\right)}^T\\ {\left({\mathbf{r}}_{i,j-1}-{\mathbf{r}}_{i,j}\right)}^T\\ \dots \\ {\left({\mathbf{r}}_{i+1,j+1}-{\mathbf{r}}_{i,j}\right)}^T\end{array}\right)\text{ and }{\underline{\underline{\mathbf{W}}}}_{i,j}^{\pm }=\left(\begin{array}{cccc}{\omega }_{i,j}^{\pm }(-1,-1)& & & \\ & {\omega }_{i,j}^{\pm }(0,-1)& & \\ & & \ddots & \\ & & & {\omega }_{i,j}^{\pm }(1,1)\end{array}\right),$$

$${\omega }_{i,j}^{\pm }(p,q)=\{\begin{array}{l}1,(p,q)\in N_{i,j}^{\pm }\\ 0,(p,q)\in N_{i,j}^{\mp }\end{array}.$$

$${\underline{\underline{\mathbf{D}}}}_{i,j}^{\pm }=\left(\begin{array}{cccc}{\mathbf{d}}_{i,j,-1,-1}^{\pm }& {\mathbf{d}}_{i,j,0,-1}^{\pm }& \cdots & {\mathbf{d}}_{i,j,1,1}^{\pm }\end{array}\right),$$

$${\partial }_{\mathbf{n}}{u}^{\pm }({\mathbf{r}}_{i,j}^{\text{pr}})$$

$$c_{i,j,p,q}^{\pm }={\mathbf{n}}_{i,j}^T{\mathbf{d}}_{i,j,p,q}^{\pm },(p,q)\in N_{i,j}^{\pm },\text{and}{c}_{i,j}^{\pm }=-\sum _{(p,q)\in N_{i,j}^{\pm }}{c}_{i,j,p,q}^{\pm }.$$

$$u_{i,j}^{-}$$

$$u_{i,j}^{+}$$

$$u_{i,j}^{-}$$

$$u_{i,j}^{+}$$

$$\varphi (x,y)=\sqrt{x^2+y^2}-r_0\left(1+\sum _{k=1}^3{\beta }_k\cos \left(n_k\left(\arctan \left(\frac{y}{x}\right)-{\theta }_k\right)\right)\right),$$

$$r_0=0.483,\left(\begin{array}{c}{n}_1\\ {\beta }_1\\ {\theta }_1\end{array}\right)=\left(\begin{array}{c}3\\ 0.1\\ 0.5\end{array}\right),\left(\begin{array}{c}{n}_2\\ {\beta }_2\\ {\theta }_2\end{array}\right)=\left(\begin{array}{c}4\\ -0.1\\ 1.8\end{array}\right)\text{and}\left(\begin{array}{c}{n}_3\\ {\beta }_3\\ {\theta }_3\end{array}\right)=\left(\begin{array}{c}7\\ 0.15\\ 0\end{array}\right).$$

$$\varphi (x,y,z)=\sqrt{x^2+y^2+z^2}-r_0\left(1+{\left(\frac{x^2+y^2}{x^2+y^2+z^2}\right)}^2\sum _{k=1}^3{\beta }_k\cos \left(n_k\left(\arctan \left(\frac{y}{x}\right)-{\theta }_k\right)\right)\right),$$

$$\varphi (r,\theta )=\sqrt{(r\cos (\theta -{\theta }_0)-x_c{)}^2+(r\sin (\theta -{\theta }_0)-y_c{)}^2}-r_0\left(1+\delta \cos (m(\theta -{\theta }_0))\right).$$

Full text

1 Introduction

It is crucial, for simulating important processes in the physical and life sciences, to find the numerical solution of elliptic equations with discontinuities in the diffusion coefficient, the source term, the solution and its flux. In the case of interfacial flows for example, jump conditions describe the discontinuity in stress that is balanced by forces that exist between phases [11]. In the simulation of protein folding, it is the electrostatic potential that has a jump across the protein’s Solvent-Excluded Surface [18, 45, 44, 64]. Other examples include solidification of multicomponent alloys [61, 37, 9] or any diffusion dominated processes with different materials properties. At the macroscale, changes across the surface can only be represented by sharp jumps, hence the need to numerically represent them as such. Failure to do so introduces errors that change the characteristics of the problem.

Numerical approximations to solve such problems have been proposed and fall into two categories, depending on whether the interface is represented explicitly or implicitly. For example, finite element discretizations approximate the space in which the solution is defined and rely on a mesh that explicitly describes the surface [4]. It is straightforward to impose boundary conditions in that framework, which is ideally suited for cases where deformations are small. For large deformations, difficulties associated with the mesh generation process are severe. Consequently, in this case, implicit representations of the interface have proved to be a better choice; imposing jump conditions, however, is a difficulty task in that framework. One of the first attempts is the Immersed Interface Method (IIM), where the jump conditions are combined with Taylor expansions of the solution on each side of the interface in order to modify the stencils of grid points adjacent to the interface. The main difficulties are the need to evaluate high-order jump conditions and surface derivatives. Several authors have further developed numerical methods within the IIM framework, e.g. [13, 39, 40, 63, 8, 1, 2, 3, 61]. Another approach is the Ghost Fluid Method (GFM) [20], first developed to treat shocks and contact discontinuities in compressible flows. The idea is to define a ghost fluid in the regions across the discontinuities by adding the interface jump to the true fluid. This simple treatment avoids the large error incurred by differentiating discontinuous solutions, and thus gives an elegant framework to manage jump conditions. The idea of the GFM was used for solving the Poisson equation with jump conditions in [41]. In this case, the jump in the normal derivative of the solution is projected onto the Cartesian directions in order to use a dimension-by-dimension approach. The authors showed that the normal jump is accurately captured, while the tangential jump is smeared; this, in turn, leads to a lack of convergence in the flux. The Voronoi Interface Method [30] solved that problem by first constructing a local Voronoi mesh adjacent to the interface and by then considering a GFM treatment. In that case, the solution is second-order accurate in the $L^{\infty}$-norm with first-order accurate fluxes in the same norm. This method has been applied to electroporation problems [31, 46], where the unknown is the electric potential at each grid points. While this method produces symmetric positive definite linear systems and only requires the right-hand side of the linear system to be modified, it requires the generation of a local Voronoi mesh and interpolation of numerical solutions from such unstructured meshes back onto Cartesian grids, which may add some challenges, especially in three spatial dimensions. The literature on solving elliptic problems with jump conditions is quite vast and we refer the interested reader to the review [28] and to other approaches, such as cut-cell approaches [16, 52], discontinuous Galerkin and the eXtended Finite Element Method (XFEM) [38, 32, 48, 17, 7, 47, 34, 21, 29, 62], the Virtual Node Method [49, 6, 49, 60, 56, 36] or other fictitious domain approaches [15, 14, 23].

In this work, we propose a finite volume discretization for elliptic interface problems in a similar vein as in [51, 55, 10] for the treatment of Neuman and Robin boundary conditions. To take into account the jump conditions, we adopt the ideas of relating the values of discontinuous functions using Taylor expansions in the normal direction and employing one-sided local least-square interpolations. In these aspects, the present method is similar to the variant of the augmented immersed interface method used in [33, 65]. The differences are, however, substantial: first, a finite volume approach is used instead of finite differences; second, the resulting linear system contains no augmented variable, which makes it straightforward to invert with “black-box” linear solvers (BiCGSTAB is used in this work; in [33, 65] the system is solved iteratively using GMRES in conjunction with a fast Poisson solver); third, the method does not involve quadratic terms in Taylor expansions and local interpolations are linear (compared to cubic ones in [33, 65]), which keeps the discretization stencil quite compact while still resulting in second-order accurate solutions (thanks to the finite volume approach). The small stencil size and the simplicity of inverting the resulting linear system make the method a good candidate for parallelization and application in the context of adaptive grids, which will be demonstrated in this work as well. We consider a level-set representation of the interface so that the method can be used in free boundary problems [54, 53, 58, 27].

2 Numerical Discretization

Consider a rectangular domain $\Omega=[x_{\min};x_{\max}]\times[y_{\min};y_{\max}]$ with an immersed irregular interface $\Gamma$ that splits $\Omega$ into two sets $\Omega^{-}$ and $\Omega^{+}$ as illustrated in Fig. 1(a). We seek a numerical solution $u=u(\boldsymbol{r})$, with $\boldsymbol{r}=(x,y)$, to the following problem:

[TABLE]

where the functions $k^{\pm}=k^{\pm}(\boldsymbol{r})$, $\mu^{\pm}=\mu^{\pm}(\boldsymbol{r})\geq\epsilon>0$, $f^{\pm}=f^{\pm}(\boldsymbol{r})$, $\boldsymbol{r}\in\Omega^{\pm}$, and $\alpha=\alpha(\boldsymbol{r})$, $\beta=\beta(\boldsymbol{r})$, $\boldsymbol{r}\in\Gamma$, are given. We denote by $[q]$ the jump of a scalar quantity $q$ across $\Gamma$, i.e. $[q]=q^{+}-q^{-}$. For simplicity, we impose Dirichlet boundary conditions on the boundary of the computation domain, i.e. $u^{+}=g$ on $\partial\Omega$, where $g=g(\boldsymbol{r})$ is given.

We discretize the domain $\Omega$ into a uniform rectangular grid of $N_{x}\times N_{y}$ points with spatial steps

[TABLE]

and associate with each point $\boldsymbol{r}_{i,j}=\left(x_{i},\,y_{j}\right)=\left(x_{\min}+(i-1)\Delta x,\,y_{\min}+(j-1)\Delta y\right)$ a finite volume $\mathcal{V}_{i,j}=\left[x_{i}-{\frac{1}{2}}\Delta x;x_{i}+{\frac{1}{2}}\Delta x\right]\times\left[y_{j}-{\frac{1}{2}}\Delta y;y_{j}+{\frac{1}{2}}\Delta y\right]$, $i\in\left[2;N_{x}-1\right]$, $j\in\left[2;N_{y}-1\right]$ (see Fig. 1(b)). The Level-Set Method [54] is used to describe the irregular interface $\Gamma$. That is, we use a Lipschitz-continuous function $\phi(\boldsymbol{r})$ such that $\Omega^{+}=\left\{\boldsymbol{r}:\,\phi(\boldsymbol{r})>0\right\}$, $\Omega^{-}=\left\{\boldsymbol{r}:\,\phi(\boldsymbol{r})<0\right\}$ and $\Gamma=\left\{\boldsymbol{r}:\,\phi(\boldsymbol{r})=0\right\}$.

At the grid points for which the interface $\Gamma$ does not cross the finite volumes, equation (1) is discretized using the standard five-point stencil. Let us consider a point $\boldsymbol{r}_{i,j}$ with its finite volume $\mathcal{V}_{i,j}$ crossed by $\Gamma$. Integrating equations (1) over $\mathcal{V}_{i,j}$ and applying the divergence theorem, one gets the following expression:

[TABLE]

where the superscript s refers to the sign ±.

Following [51, 55], that is, approximating the domain integrals by the integrand value multiplied by the corresponding volumes, and estimating the fluxes between cells using values at nearest-neighbor grid points and central difference formulas, one obtains:

[TABLE]

where $\mathcal{D}$ is the problem dimensionality, $h=\max(\Delta x,\Delta y)$, $|\mathcal{V}_{i,j}^{\pm}|$ denotes the volume of $\mathcal{V}_{i,j}\cap\Omega^{\pm}$, $u_{i,j}=u(\boldsymbol{r}_{i,j})$, $A^{\pm}_{i\pm{\frac{1}{2}},j}$ and $A^{\pm}_{i,j\pm{\frac{1}{2}}}$ are face areas of $\mathcal{V}_{i,j}^{\pm}$ in the $x$- and $y$-directions, respectively. To compute the boundary and domain integrals required by the proposed discretization, we use the geometric reconstruction approach from [42]. In case when an immersed interface is only piece-wise smooth the method from [10] can be used.

The discretization given by equation (4) requires that both values of $u^{-}$ and $u^{+}$ be available at grid points with a control volume crossed by $\Gamma$, thus, one more equation is required at such grid points for the system of equations to be uniquely invertible. We derive the additional equation based on the jump conditions (2)-(3) and Taylor expansions of $u^{\pm}$ in the normal to the interface direction. Moreover, equations derived in this way can be used to express $u^{\pm}_{i,j}$ near the interface as a function of $u^{\mp}_{i,j}$ and the jump conditions (2)-(3). As a result, it enables us to eliminate the additional degrees of freedom, i.e., reduce the system’s size back to $N_{x}\times N_{y}$, and make its structure more homogeneous (all of the equations in the linear system are of the same type). This is expected to make the system of equations less difficult to invert by black-box linear solvers. We select the $N_{x}\times N_{y}$ unknowns to solve for as:

[TABLE]

We then develop formulas to express $u^{+}_{i,j}$ for $\boldsymbol{r}_{i,j}\in\Omega^{-}$ and for $u^{-}_{i,j}$ for $\boldsymbol{r}_{i,j}\in\Omega^{+}$ as a function of the unknowns $u_{i,j}$. This is described next.

Consider a grid point $\boldsymbol{r}_{i,j}$ near the interface $\Gamma$ and its projection, $\boldsymbol{r}_{i,j}^{\text{pr}}$, onto the interface (see Fig. 1(c)). Taylor expansion relates the values of $u^{\pm}$ at $\boldsymbol{r}_{i,j}$ and $\boldsymbol{r}_{i,j}^{\text{pr}}$ as:

[TABLE]

where $\delta_{i,j}$ is the signed distance from $\boldsymbol{r}_{i,j}$ to $\boldsymbol{r}_{i,j}^{\text{pr}}$ ($\pm\delta_{i,j}>0$ if $\boldsymbol{r}_{i,j}\in\Omega^{\pm}$). The geometrical quantities $\boldsymbol{n}(\boldsymbol{r}_{i,j}^{\text{pr}})$, $\boldsymbol{r}_{i,j}^{\text{pr}}$ and $\delta_{i,j}$ are estimated from the level-set function as:

[TABLE]

Subtracting $u^{-}_{i,j}$ from $u^{+}_{i,j}$ given in (6) and taking into account the jump condition (2) one obtains:

[TABLE]

Furthermore, eliminating either $\partial_{\boldsymbol{n}}u^{+}(\boldsymbol{r}_{i,j}^{\text{pr}})$ or $\partial_{\boldsymbol{n}}u^{-}(\boldsymbol{r}_{i,j}^{\text{pr}})$ in the above expression using the jump condition (3) results in the following two equations:

[TABLE]

If one approximates either $\partial_{\boldsymbol{n}}u^{+}(\boldsymbol{r}_{i,j}^{\text{pr}})$ or $\partial_{\boldsymbol{n}}u^{-}(\boldsymbol{r}_{i,j}^{\text{pr}})$ using $u_{i,j}^{\pm}$ and $\left\{u_{p,q},p\in[1,N_{x}],q\in[1,N_{y}]\right\}$, then these formulas can be used to eliminate additional degrees of freedom. A straightforward way to do that is to use a suitable local interpolants $u^{\pm}_{I}=u^{\pm}_{I}(\boldsymbol{r})$ of unknown functions $u^{\pm}$ as:

[TABLE]

Specifically, in this work we use the linear interpolation:

[TABLE]

where the gradient $\left(\nabla u^{\pm}\right)_{i,j}$ is found as the least-square solution satisfying the constraints:

[TABLE]

$N^{\pm}_{i,j}$ denotes the set of neighboring grid points of $\boldsymbol{r}_{i,j}$, lying in the region $\Omega^{\pm}$, that is:

[TABLE]

In other words, the local linear interpolants are constructed using available values at the nearest-neighbor grid points of $\boldsymbol{r}_{i,j}$ (in Cartesian and diagonal directions). Note also that $u_{i+p,j+q}=u^{\pm}_{i+p,j+q}$ if $(p,q)\in N^{\pm}_{i,j}$.

Thus, the gradient $\left(\nabla u^{\pm}\right)_{i,j}$ is the least-squares solution of the following linear system:

[TABLE]

that is:

[TABLE]

where the $3^{\mathcal{D}}\times\mathcal{D}$ and $3^{\mathcal{D}}\times 3^{\mathcal{D}}$ matrices $\underline{\underline{\boldsymbol{X}}}_{i,j}$ and $\underline{\underline{\boldsymbol{W}}}_{i,j}$ are given by:

[TABLE]

with

[TABLE]

Substitution of $u^{\pm}_{I}(\boldsymbol{r})$ into $\eqref{eq::norm_deriv}$ yields approximations of $\partial_{\boldsymbol{n}}u^{\pm}(\boldsymbol{r}_{i,j}^{\text{pr}})$ as linear combinations of $\left\{u_{i+p,j+q},p=-1,0,1,q=-1,0,1\right\}$. Specifically, let us write the $\mathcal{D}\times 3^{\mathcal{D}}$ matrix $\underline{\underline{\boldsymbol{D}}}^{\pm}_{i,j}$ as:

[TABLE]

where vectors $\boldsymbol{d}_{i,j,-1,-1}^{\pm}$, $\ldots$, $\boldsymbol{d}_{i,j,1,1}^{\pm}$ represent the columns of the matrix $\underline{\underline{\boldsymbol{D}}}^{\pm}_{i,j}$. Then the normal derivative can be expressed as:

[TABLE]

where the coefficients are given by:

[TABLE]

Substitution of (10) into (7) produces formulas expressing $u^{+}_{i,j}$ and $u^{-}_{i,j}$ in terms of the selected $N_{x}\times N_{y}$ unknowns $\left\{u_{p,q},p\in[1,N_{x}],q\in[1,N_{y}]\right\}$. Combining them with the definition (5), we get the following two sets of rules (which are $\mathcal{O}\left(h^{2}\right)$ accurate in the value of $u$): one is based on approximating $\partial_{\boldsymbol{n}}u^{+}(\boldsymbol{r}_{i,j}^{\text{pr}})$:

[TABLE]

while the other one is based on approximating $\partial_{\boldsymbol{n}}u^{-}(\boldsymbol{r}_{i,j}^{\text{pr}})$:

[TABLE]

Thus, one has a certain flexibility in constructing the final discretization. For example, one could choose, for each $\boldsymbol{r}_{i,j}$, the formula based on approximating $\partial_{\boldsymbol{n}}u^{-}(\boldsymbol{r}_{i,j}^{\text{pr}})$ or $\partial_{\boldsymbol{n}}u^{+}(\boldsymbol{r}_{i,j}^{\text{pr}})$ depending on the largest number of neighboring points of $\boldsymbol{r}_{i,j}$ that are in $\Omega^{-}$ or in $\Omega^{+}$ (let us denote this scheme as Random). However, this choice would ignore the magnitude of the diffusion constants $\mu^{-}$ and $\mu^{+}$ and their influence on the condition number of the linear system. To investigate this issue, we consider two additional schemes: the first one (referred to as Bias Fast) uses interpolation in the fast-diffusion region (i.e., if $\mu^{-}>\mu^{+}$ then the formula based on $\partial_{\boldsymbol{n}}u^{-}(\boldsymbol{r}_{i,j}^{\text{pr}})$ is used); the second scheme (referred to as Bias Slow) uses interpolation in the slow-diffusion region (i.e., if $\mu^{-}>\mu^{+}$ then we use the formula based on $\partial_{\boldsymbol{n}}u^{+}(\boldsymbol{r}_{i,j}^{\text{pr}})$).

Remarks:

In the limiting cases $\frac{\mu^{-}}{\mu^{+}}\rightarrow\infty$ or $\frac{\mu^{-}}{\mu^{+}}\rightarrow 0$, only the scheme Bias Slow remains well defined, thus, we expect it to perform the best and be well-conditioned for any ratio of diffusion coefficients. We will illustrate in section 3 that only the scheme Bias Slow produces a condition number that is bounded. This is consistent with the results reported in [66] describing a variant of the augmented immersed interface method.

• 2.

In the limiting case $\frac{\mu^{-}}{\mu^{+}}\equiv 1$, the three schemes coincide. Moreover, the matrix associated with the resulting linear system is the same as for the case when no interface is present (that is, as for the standard five-point stencil in 2D) and only the right-hand is changed to account for jump conditions.

• 3.

The truncation error is the same for all three schemes. Therefore, we expect them to have similar accuracies. Specifically, the truncation error111After scaling the resulting discretization by the cell volume $\mathcal{O}\left(h^{\mathcal{D}}\right)$ to account for the integration of the PDEs over a finite volume. is $\mathcal{O}\left(h^{2}\right)$ for grid points away from the immersed interface and $\mathcal{O}\left(1\right)$ for cells crossed by the interface. Following the results of [35, 26, 57, 12, 50, 55, 22, 24, 10], we expect the schemes to produce second-order accurate numerical solutions with first-order accurate gradients.

• 4.

The truncation error can be improved to be $\mathcal{O}\left(h\right)$ for cells crossed by the interface by making the following changes in the discretization scheme: 1) Estimate fluxes between cells at the centroids of cell faces using linear interpolation as done, for example, in [35, 10]; 2) Retain the quadratic term in Taylor expansion (6); 3) Use a quadratic interpolant (instead of linear (9)) to approximate $\partial_{\boldsymbol{n}}u^{\pm}$ at projection points. This is expected to increase the accuracy of solution gradients to second order. We leave this investigation to future work.

• 5.

In general the case where $\mu^{+}\neq\mu^{-}$, the resulting linear system is nonsymmetric. In the worst case scenario the computational stencil involves nearest neighbors (both in Cartesian and diagonal directions) of the standard five-point stencil as illustrated in Fig. 2(a). An example of the matrix associated with the resulting linear system structure is shown in Fig. 2(b).

3 Numerical tests

To numerically illustrate the properties of the proposed schemes, we study three characteristics: the order of accuracy of the numerical solution in the $L^{\infty}$-norm, the order of accuracy of the numerical gradients in the $L^{\infty}$-norm, and the condition number of the linear system, estimated by the MATLAB condest function. We consider two tests: the first one, the convergence test, studies the dependence of those three characteristics on the grid resolution. The second one, the conditioning test, focuses on the dependence of the three characteristics on the ratio, $\frac{\mu^{-}}{\mu^{+}}$, of the diffusion coefficients. We perform both tests in two and three spatial dimensions. In all the examples, we use the implementation of the BiCGStab algorithm provided by PETSc [5] with the Hypre preconditioner [19].

3.1 Two-dimensional case

Consider an annular region222We enclose an immersed interface inside another region in order to be able to obtain results for different placements of the immersed interface and the computational grid without changing the problem statement. On the boundaries of the enclosing region, Dirichlet boundary conditions can be imposed with any of the methods [26, 25, 59] with inner and outer radii $r_{i}=0.151$ and $r_{e}=0.911$, and an immersed star-shaped interface (see Fig. 3(a)), described by the following level-set function:

[TABLE]

with parameters:

[TABLE]

Using the method of manufactured solutions, we take the exact solution to be $u^{-}=\sin(2x)\cos(2y)$ and $u^{+}=\left(16\left(\frac{y-x}{3}\right)^{5}-20\left(\frac{y-x}{3}\right)^{3}+5\left(\frac{y-x}{3}\right)\right)\log\left(x+y+3\right)$ (see Fig. 3(c)). For the convergence test, we set the diffusion coefficients to $\mu^{-}=10\left(1+\frac{1}{5}\cos(2\pi(x+y))\sin(2\pi(x-y))\right)$ and $\mu^{+}=1$ (see Fig. 3(b)), and we vary the grid resolution from $2^{-4}$ to $2^{-9}$. For the conditioning test, we fix the grid resolution at $2^{-6}$ and $\mu^{+}=1$ and vary $\mu^{-}$ from $10^{-4}$ to $10^{4}$. The results are presented in Fig. 4 and 5, where each data point represents the maximum value among $10\times 10=100$ different relative placements of the immersed interface on the computational grid (as done in [10]). The different placements thus account for cases where the interface defines a control volume that is arbitrarily small or large, relative to an elementary grid cell. Section 3.3 will draw some conclusions from these results.

3.2 Three-dimensional case

Consider a spherical shell333As in the two-dimensional case, Dirichlet boundary conditions are enforced on the boundaries of the enclosing region with inner and outer radii $r_{i}=0.151$ and $r_{e}=0.911$, and an immersed star-shaped interface described by the level-set function:

[TABLE]

with the same parameters (11) as for the two-dimensional case. The problem geometry is illustrated in Fig. 6. The exact solutions are taken to be $u^{-}=\sin(2x)\cos(2y)\exp(z)$ and $u^{+}=\left(16\left(\frac{y-x}{3}\right)^{5}-20\left(\frac{y-x}{3}\right)^{3}+5\left(\frac{y-x}{3}\right)\right)\log\left(x+y+3\right)\cos(z)$. In the convergence test, the diffusion coefficients are set to $\mu^{-}=10\left(1+\frac{1}{5}\cos(2\pi(x+y))\sin(2\pi(x-y))\cos(z)\right)$ and $\mu^{+}=1$. In the conditioning test, the grid resolution is fixed at $2^{-4}$, $\mu^{+}=1$ and $\mu^{-}$ is varied from $10^{-4}$ to $10^{4}$. The test results are presented in Fig. 7 and 8 ,where each data point is obtained as the maximum (worse) value among $5\times 5\times 5=125$ different relative placements of the immersed interface on the computational grid. As for the two dimensional case, section 3.3 will draw some conclusions from these results.

3.3 Analysis

From the results presented in sections 3.1 and 3.2, it is clear that the numerical schemes have the same behavior in two and three spatial dimensions. The convergence test results (see Fig. 4 and 7) indicate that, for a moderate diffusion coefficient ratio, all three schemes have comparable convergence properties: the numerical solutions are second-order accurate with first-order accurate gradients in the $L^{\infty}$-norm. The condition number scales with the grid resolution as $h^{-2}$, which is similar to the scaling of the condition number for the standard five-point (nine-point) stencil in 2D (3D). The only difference between the three schemes is the magnitude of the errors and the magnitude of the condition numbers, with the scheme Bias Slow giving the best results.

On the other hand, the conditioning test in two and three spatial dimensions demonstrate that the three schemes behaviors are drastically different when the ratio $\frac{\mu^{-}}{\mu^{+}}$ varies (Fig. 5 and 8). In particular, the condition numbers for the schemes Random and Bias Fast grow unboundedly as the ratio of the diffusion coefficients either decreases or increases away from 1. As a result, the magnitude of the errors in the solution and its gradient grow significantly. We also note that, for approximately $\frac{\mu^{-}}{\mu^{+}}>10$ and $\frac{\mu^{-}}{\mu^{+}}<10^{-1}$, the linear solver is not able to invert the resulting linear system in a given number of iterations (we set that number to 50 in those numerical examples). In contrast, the condition number for the scheme Bias Slow converges to finite values as $\frac{\mu^{-}}{\mu^{+}}\rightarrow 0$ or $\frac{\mu^{-}}{\mu^{+}}\rightarrow\infty$. As a result, the linear solver is able to invert the resulting linear sytem for any values of $\frac{\mu^{-}}{\mu^{+}}$ (the number of iterations depends only on the grid resolution). Moreover, the errors of the numerical solutions and their gradients are only moderately affected by small or large ratios $\frac{\mu^{-}}{\mu^{+}}$.

3.4 Application to adaptive quadtree and octree grids

Thanks to its fairly small stencil, it is simple to apply the proposed scheme in the context of adaptive grids. In this section, we demonstrate perhaps the easiest way of doing so: we consider adaptive Cartesian quadtree (octree) grids that are locally uniform around immersed interfaces. In the regions where an adaptive grid is non-uniform (and which are away from immersed interfaces), we use the second-order accurate superconvergent finite difference scheme of [43], while in the regions close to immersed interfaces (and where the grid is locally uniform), we use the proposed scheme for imposing interface jump conditions. Note that the discretization of [43] and the one described in this work reduce to the standard 5-point (9-point in three spatial dimensions) stencil on uniform grids and in the absence of immersed interfaces. This fact makes the combination of the two discretizations seamless.

In this example, we consider 10 clusters of small star-shaped uniformly charged dielectric particles in a $2$m$\times$2m vacuum domain and compute the electric field generated by this configuration. The number of particles in each cluster varies from 3 to 10, the average size of clusters is 1 cm, particles have between 2 and 6 bumps and sizes in the range $[0.1;1]$ mm. The exact problem geometry description is provided in A. The absolute permittivity of particles is $10\varepsilon_{0}$ and their charge densities are $\pm 10^{5}\varepsilon_{0}$ (the sign is assigned to each particle randomly), where $\varepsilon_{0}$ is the vacuum permittivity. The domain boundaries are assumed to be a good conductor. Thus, the electric potential $\varphi\equiv u$ satisfies the boundary value problem (1)-(3) with parameters $k^{\pm}=\alpha=\beta=0$, $\mu^{-}=10$, $\mu^{+}=1$, $f^{-}=\pm 10^{5}$ (depending which particle is considered), $f^{+}=0$, $g=0$, where we denote the particles as $\Omega^{-}$ and the vacuum as $\Omega^{+}$.

Figure 9 depicts the electric potential $\phi$ and the electric field $\boldsymbol{E}=-\nabla\phi$ computed on an adaptive grid with the coarsest and finest mesh sizes corresponding to uniform resolutions of $2^{10}\times 2^{10}$ and $2^{20}\times 2^{20}$ grid points, respectively. Such a computational grid contains $2754021$ points, which is approximately just $0.00025\%$ of the total number of points in a uniform $2^{20}\times 2^{20}$ grid.

Remark: While the strategy of enforcing adaptive grids to be locally uniform around immersed interfaces is quite common and justified from the point of view of accuracy in many situations, it may not be very efficient for certain applications; in particular those where parts of the interface has high curvature while other parts are rather flat. However, it seems straightforward to adapt the described method to fully-adaptive (non-uniform along immersed interfaces) and Voronoi grids (including Voronoi partitions generated based on adaptive quadtree and octree grids). This will be considered in future works.

4 Conclusions

We have presented a simple finite volume numerical method for solving Elliptic equations with jump conditions across irregular interfaces that are implicitly represented by a level-set function on Cartesian grids. Second-order accurate solutions and first-order accurate gradients are obtained in the $L^{\infty}$-norm. The linear system is non-symmetric but the condition number is bounded, regardless of the ratio of the diffusion coefficients, so that the linear system can be inverted in a constant number of iterations that depends only on the grid resolution: the condition number scales as $\mathcal{O}\left(h^{-2}\right)$, similarly to the linear system obtained from the standard five-point stencil.

Future work will be focused on the analysis of the numerical scheme and its properties, the improvement to a superconvergent scheme (i.e, with second-order accurate gradients) and the extension to fully-adaptive (non-uniform along immersed interfaces) grids.

Acknowledgement

This research was supported by ARO W911NF-16-1-0136 and ONR N00014-17-1-2676.