High-Order Numerical Methods for 2D Parabolic Problems in Single and Composite Domains

Abstract

In this work, we discuss and compare three methods for the numerical approximation of constant- and variable-coefficient diffusion equations in both single and composite domains with possible discontinuity in the solution/flux at interfaces, considering (i) the Cut Finite Element Method; (ii) the Difference Potentials Method; and (iii) the summation-by-parts Finite Difference Method. First we give a brief introduction for each of the three methods. Next, we propose benchmark problems, and consider numerical tests—with respect to accuracy and convergence—for linear parabolic problems on a single domain, and continue with similar tests for linear parabolic problems on a composite domain (with the interface defined either explicitly or implicitly). Lastly, a comparative discussion of the methods and numerical results will be given.

1 Introduction

Designing methods for the high-order accurate numerical approximation of partial differential equations (PDE) posed on composite domains with interfaces, or on irregular and geometrically complex domains, is crucial in the modeling and analysis of problems from science and engineering. Such problems may arise, for example, in materials science (models for the evolution of grain boundaries in polycrystalline materials), fluid dynamics (the simulation of homogeneous or multi-phase fluids), engineering (wave propagation in an irregular medium or a composite medium with different material properties), biology (models of blood flow or the cardiac action potential), etc. The analytic solutions of the underlying PDE may have non-smooth or even discontinuous features, particularly at material interfaces or at interfaces within a composite medium. Standard numerical techniques involving finite-difference approximations, finite-element approximation, etc., may fail to produce an accurate approximation near the interface, leading one to consider and develop new techniques.

There is extensive existing work addressing numerical approximation of PDE posed on composite domains with interfaces or irregular domains, for example, the boundary integral method [11, 56], difference potentials method [3, 6, 26, 27, 58, 65], immersed boundary method [30, 42, 61, 74], immersed interface method [2, 46, 47, 49, 69], ghost fluid method [31, 32, 50, 51], the matched interface and boundary method [82, 84,85,86], Cartesian grid embedded boundary method [19, 41, 57, 83], multigrid method for elliptic problems with discontinuous coefficients on an arbitrary interface [18], virtual node method [9, 39], Voronoi interface method [35, 36], the finite difference method [8, 10, 24, 75, 78, 79] and finite volume method [22, 34] based on mapped grids, or cut finite element method [13,14,15, 37, 38, 71, 76]. Indeed, there have been great advances in numerical methods for the approximation of PDE posed on composite domains with interfaces, or on irregular domains. However, it is still a challenge to design high-order accurate and computationally-efficient methods for PDE posed in these complicated geometries, especially for time-dependent problems, problems with variable coefficients, or problems with general boundary/interface conditions.

The aim of this work is to establish benchmark (test) problems for the numerical approximation of parabolic PDE defined in irregular or composite domains. The considered models (Sect. 2) arise in the study of mass or heat diffusion in single or composite materials, or as simplified models in other areas (e.g., biology, materials science, etc.). The formulated test problems (Sect. 4) are intended (a) to be suitable for comparison of high-order accurate numerical methods—and will be used as such in this study—and (b) to be useful in further research. Moreover, the proposed problems include a wide variety of possibilities relevant in applications, which any robust numerical method should resolve accurately, including constant diffusion; time-varying diffusion; high frequency oscillations in the analytical solution; large jumps in diffusion coefficients, solution, and/or flux; etc. For now, we will consider a simplified geometrical setting, with the intent of setting a “baseline” from which further research, or more involved comparisons, might be conducted. Therefore, in Sect. 2 we will introduce two circular geometries, which are defined either explicitly, or implicitly via a level set function.

In Sect. 3, we briefly introduce the numerical methods we will consider in this work, i.e., second- and fourth-order versions of (i) the Cut Finite Element Method (cut-FEM); (ii) the Difference Potentials Method (DPM), with Finite Difference approximation as the underlying discretization in the current work; and (iii) the summation-by-parts Finite Difference Method combined with the simultaneous approximation term technique (SBP–SAT–FD). These three methods are all modern numerical methods which may be designed for problems in irregular or composite domains, allowing for high-order accurate numerical approximation, even at points close to irregular interfaces or boundaries. We will apply each method to the formulated benchmark problems, and compare results. From the comparisons, we expect to learn what further developments of the methods at hand would be most important.

To resolve geometrical features of irregular domains, both cut-FEM and DPM use a Cartesian grid on top of the domain, which need not conform with boundaries or interfaces. These types of methods are often characterized as “immersed” or “embedded”. In the finite difference framework, embedded methods for parabolic problems are developed in [1, 23]. For comparison with cut-FEM and DPM, however, in this paper we use a finite difference method based on a conforming approach. The finite difference operators we use satisfy a summation-by-parts principle. Then, in combination with the SAT method to weakly impose boundary and interface conditions, an energy estimate of the semi-discretization can be derived to ensure stability. In addition, we use curvilinear grids and transfinite interpolation to resolve complex geometries.

For recent work on SBP–SAT–FD for wave equations in composite domains, see [10, 17, 75, 78], and the two review papers [21, 73]; for recent work in DPM for elliptic/parabolic problems in composite domains with interface defined explicitly, see [3,4,5,6, 26,27,28, 58, 59, 65]; and for recent work in cut-FEM, see [13,14,15, 37, 53, 71].

The paper is outlined as follows. In Sect. 2, we give brief overview of the continuous formulation of the parabolic problems in a single domain or a composite domain. In Sect. 3, we give introductions to the basics of the three proposed methods: cut-FEM, DPM, and SBP–SAT–FD. In Sect. 4, we formulate the numerical test problems. In Sect. 5, we present extensive numerical comparisons of errors and convergence rates, between the second- and fourth-order versions of each method. The comparisons include single domain problems with constant or time-dependent diffusivity; and interface problems with interface defined explicitly, or implicitly by a level set function. In Sect. 6, we give a comparative discussion of the three methods and the numerical results, together with a discussion on future research directions. Lastly, in Sect. 7, we give our concluding remarks.

Fig. 1

The (a) single domain $\varOmega$ and (b) composite domain $\varOmega = \varOmega _1 \cup \varOmega _2$. In (b), $\partial \varOmega _1$ has two connected components: the boundary $\partial \varOmega$ and interface $\varGamma = \partial \varOmega _2$

2 Statement of Problem

In this section, we describe two diffusion problems, which will be the setting for our proposed benchmark (test) problems in Sect. 4. (Recall from Sect. 1 that these models arise, for example, in the study of mass or heat diffusion.) For brevity, in the following discussion, we denote $u := u(x, y, t)$ and $u_s := u_s(x, y, t)$, with $s = 1, 2$.

2.1 The Single Domain Problem

First, we consider the linear parabolic PDE on a single domain $\varOmega$ (e.g., Fig. 1a), with variable diffusion $\lambda (t)$:

$$\begin{aligned} \dfrac{\partial u}{\partial t} = \nabla \cdot (\lambda (t) \nabla u) + f(x, y, t), \ \ (x, y, t) \in \varOmega \times (0, T], \end{aligned}$$

(1)

subject to initial and Dirichlet boundary conditions:

$$\begin{aligned} u(x, y, 0) = u^0(x, y), \ \ (x, y) \in \varOmega \quad \text{ and } \quad u = \psi (x, y, t), \ \ (x, y, t) \in \partial \varOmega \times (0, T]. \end{aligned}$$

(2)

Here, the initial and boundary data $u^0(x, y)$ and $\psi (x, y, t)$, the diffusion coefficient $\lambda (t)$, the forcing function f(x, y, t), and the final time T are known (given) data.

2.2 The Composite Domain Problem

Next, we consider the linear parabolic PDE on a composite domain $\varOmega := \varOmega _1 \cup \varOmega _2$ (e.g., Fig. 1b), with constant diffusion coefficients $(\lambda _1, \lambda _2)$:

$$\begin{aligned} \dfrac{\partial u_1}{\partial t}&= \nabla \cdot (\lambda _1 \nabla u_1) + f_1(x, y, t), \ \ (x, y, t) \in \varOmega _1 \times (0, T], \end{aligned}$$

(3)

$$\begin{aligned} \dfrac{\partial u_2}{\partial t}&= \nabla \cdot (\lambda _2 \nabla u_2) + f_2(x, y, t), \ \ (x, y, t) \in \varOmega _2 \times (0, T], \end{aligned}$$

(4)

subject to initial conditions:

$$\begin{aligned} u_1(x, y, 0)&= u_1^0(x, y), \ \ (x, y) \in \varOmega _1, \end{aligned}$$

(5)

$$\begin{aligned} u_2(x, y, 0)&= u_2^0(x, y), \ \ (x, y) \in \varOmega _2, \end{aligned}$$

(6)

Dirichlet boundary conditions:

$$\begin{aligned} u_1 = \psi (x, y, t), \ \ (x, y, t) \in \partial \varOmega \times (0, T], \end{aligned}$$

(7)

and interface/matching conditions:

$$\begin{aligned} u_1 - u_2&= \mu _1(x, y, t), \ \ (x, y, t) \in \varGamma \times (0, T], \end{aligned}$$

(8)

$$\begin{aligned} \lambda _1 \frac{\partial u_1}{\partial n} - \lambda _2 \frac{\partial u_2}{\partial n}&= \mu _2(x, y, t), \ \ (x, y, t) \in \varGamma \times (0, T]. \end{aligned}$$

(9)

In formula (9), $\frac{\partial u_s}{\partial n}$, $s = 1, 2$ denotes the normal derivative at the interface $\varGamma$, i.e., $\frac{\partial u_s}{\partial n} = \nabla u_s \cdot \mathbf {n}$, where $\mathbf {n}$ is the outward unit normal vector at the interface $\varGamma$. The initial, boundary, and interface data $u_1^0(x, y)$, $u_2^0(x, y)$, $\psi (x, y, t)$, $\mu _1(x, y, t)$, and $\mu _2(x, y, t)$; the diffusion coefficients $(\lambda _1, \lambda _2)$; the forcing functions $f_1(x,y,t)$ and $f_2(x,y,t)$; and the final time T are some known (given) data.

Remark 1

We consider the circular geometries depicted in Fig. 1 as the geometrical setting for our proposed benchmark problems in this work. In applications (Sect. 1), other geometries will likely be considered, some much more complicated than Fig. 1. While our methods can handle more complicated geometry, this is (to the best of our knowledge) the first work looking to establish benchmarks—and compare numerical methods—for parabolic interface problems (3–9). As such, we think that the geometries in Fig. 1 are a good “baseline”—without all the added complexities that more complicated geometries might produce—from which further research, or more involved comparisons, might be done.

To be more specific, we aim to define a simple set of test problems that can be easily implemented and tested for any numerical scheme of interest. With circular domains, it suffices for us to compare/contrast performance of the numerical methods on a simple geometry with smooth boundary versus on a composite domain with fixed interface (explicit or implicit). The approximation of the solution to such composite-domain problems are already challenging for any numerical methods, since (i) the solution may fail to be smooth (or may be discontinuous) at the interface, and (ii) there may be discontinuous material coefficients ($\lambda _1 \not = \lambda _2$).

Remark 2

For both the single and composite domain problems, we could also consider other boundary conditions, e.g., a Neumann boundary condition as in [6, 13], etc.

Fig. 2

The (a) subdomain $\varOmega _1$ immersed in a mesh $\mathcal {T}_{1}$, (b) subdomain $\varOmega _2$ immersed in a mesh $\mathcal {T}_{2}$ and domain $\varOmega$ immersed in $\mathcal {T}$, and (c) intersected elements $\mathcal {T}_\varGamma$

3 Overview of Numerical Methods

3.1 Cut-FEM

In this section, we give a brief presentation of the cut-FEM method. For a more detailed presentation of cut-FEM, see, for example, [13, 14, 53].

Let $\varOmega _s$ be covered by a structured triangulation, $\mathcal {T}_{s}$, so that each element $T\in \mathcal {T}_{s}$ has some part inside of $\varOmega _s$; see Fig. 2a and b. Here, $s = 1,2$ is an index for the composite domain problem (3–9), which will be omitted when referring to the single domain problem (1, 2). (For the latter, note that $\mathcal {T}$ covers $\varOmega$.) Typically $\mathcal {T}_1$ and $\mathcal {T}_2$ would be created from a larger mesh by removing some of the cells. Further, let $\mathcal {T}_\varGamma =\{T\in \mathcal {T}: T \cap \varGamma \ne \emptyset \}$ be the set of intersected elements; see Fig. 2c. In the following, we shall use $\varGamma$ both for the immersed boundary of the single domain problem and for the immersed interface of the composite domain problem, in order to make the connection to the set $\mathcal {T}_\varGamma$ clearer.

To construct the finite element spaces we use Lagrange elements with Gauss–Lobatto nodes of order p ($Q_p$-elements). Let $V_h^s$ denote a continuous finite element space on $\varOmega _s$, consisting of $Q_p$-elements on the mesh $\mathcal {T}_s$:

$$\begin{aligned} V_h^s=\left\{ v\in C^0(\varOmega _s): \left. v\right| _{T}\in Q_p(T), \, T\in \mathcal {T}_s\right\} . \end{aligned}$$

(10)

For the single domain problem (1, 2) we solve for the solution $u \in V_h$; while for the composite domain problem (3–9), we solve for the pair $\{ {u}_{1} , {u}_{2} \}\in V_h^1 \times V_h^2$. For the latter problem, this means that the degrees of freedom are doubled over elements belonging to $\mathcal {T}_\varGamma$.

We begin by stating the weak formulation for the single domain problem (1, 2). Let $(\cdot ,\cdot )_X$ and $\left\langle \cdot ,\cdot \right\rangle _Y$ be the $L_2$ scalar products taken over the two- and one-dimensional domains $X\subset \mathbb {R}^2$ and $Y\subset \mathbb {R}^{1}$, respectively. The present method is based on modifying the weak formulation by using Nitsche’s method [60] to enforce the boundary condition (2). By multiplying (1) with a test function $v\in V_h$, and integrating by parts, we obtain:

$$\begin{aligned} (\dot{u},v)_{\varOmega }+(\lambda \nabla u, \nabla v)_{\varOmega } -\left\langle \lambda \frac{\partial u}{\partial n},v\right\rangle _{\varGamma } =(f,v)_{\varOmega }, \quad \forall v\in V_h. \end{aligned}$$

(11)

Note that (2) is consistent with the following terms:

$$\begin{aligned} \frac{\gamma _{D}}{h_T}\left\langle \lambda u,v\right\rangle _{\varGamma }&= \frac{\gamma _{D}}{h_T}\left\langle \lambda \psi ,v\right\rangle _{\varGamma }, \end{aligned}$$

(12)

$$\begin{aligned} -\left\langle u,\lambda \frac{\partial v}{\partial n}\right\rangle _{\varGamma }&= -\left\langle \psi ,\lambda \frac{\partial v}{\partial n}\right\rangle _{\varGamma }, \end{aligned}$$

(13)

where $\gamma _D$ is a constant, and $h_T$ is the side length of the quadrilaterals in the triangulation. Now, adding (12, 13) to (11) gives the following weak form: Find $u \in V_h$ such that

$$\begin{aligned} (\dot{u},v)_{\varOmega }+a(u,v)=L(v), \quad \forall v\in V_h, \end{aligned}$$

(14)

where

$$\begin{aligned} a(u,v)&=(\lambda \nabla u, \nabla v)_{\varOmega } -\left\langle \lambda \frac{\partial u}{\partial n},v\right\rangle _{\varGamma } -\left\langle u,\lambda \frac{\partial v}{\partial n}\right\rangle _{\varGamma }+\frac{\gamma _{D}}{h_T} \left\langle \lambda u,v\right\rangle _{\varGamma }, \end{aligned}$$

(15)

$$\begin{aligned} L(v)&=(f,v)_{\varOmega } +\left\langle \lambda \psi ,\frac{\gamma _{D}}{h_T} v- \frac{\partial v}{\partial n}\right\rangle _{\varGamma }. \end{aligned}$$

(16)

For $\mathcal {T}_\varGamma$ (the elements intersected by $\varGamma$), note that one must integrate only over the part of the element that lies inside $\varOmega$. A problem with this is that one cannot control how the intersections (cuts) between $\varOmega$ and $\mathcal {T}$ are made. Depending on how $\varOmega$ is located with respect to the triangulation, some elements can have an arbitrarily small intersection with the domain—see, for example, Fig. 3a. If $\varOmega$ is moved with respect to $\mathcal {T}$ to make the cut arbitrarily small, then the condition numbers of the mass and stiffness matrices can become arbitrarily large.

To mitigate this issue, in this work we add a stabilizing term j—defined shortly in (19)—to the mass and stiffness matrices, so that their condition numbers are bounded, independently of how the domain $\varOmega$ is located with respect to the triangulation $\mathcal {T}$ [14, 53]. Adding stabilization to (14) results in the following weak form: Find $u \in V_h$ such that

$$\begin{aligned} (\dot{u},v)_{\varOmega }+ \gamma _M j(\dot{u},v)+a(u,v) + \gamma _A h_T^{-2}\lambda j(u,v)=L(v), \quad \forall v\in V_h, \end{aligned}$$

(17)

where $\gamma _M$ and $\gamma _A$ are scalar constants.

Fig. 3

(a) An element having a small intersection, shown in gray, with the domain; (b) faces belonging to $\mathcal {F}_{2}$; and (c) faces belonging to $\mathcal {F}_{1}$

In order to state the definition of stabilization (19), denote by $\mathcal {F}_{s}$ the set of faces, as seen in Fig. 3b and c. That is, $\mathcal {F}_{s}$ is the set of all faces of the elements in $\mathcal {T}_\varGamma$, excluding the boundary faces of $\mathcal {T}_s$:

$$\begin{aligned} \mathcal {F}_{s}=\{F=T_A \cap T_B \; : \; T_A \in \mathcal {T}_\varGamma \;\text { or }\; T_B \in \mathcal {T}_\varGamma , \quad T_A,T_B\in \mathcal {T}_{s} \}. \end{aligned}$$

(18)

Then, the stabilization term is defined as:

$$\begin{aligned} j_s(u,v)=\sum _{F\in \mathcal {F}_{s}} \sum _{k=1}^p \frac{h_T^{2k+1}}{(2k+1)(k!)^2}\left\langle [\partial _n^k u],[\partial _n^k v]\right\rangle _F, \end{aligned}$$

(19)

where $[u] = u|_{F_+}-u|_{F_-}$ is the jump over a face, F; n refers to a normal of F; and $\partial _n^k u$ denotes the k-th order normal derivative. The scaling with respect to k of the terms in (19) is based on how the stabilization was derived. In particular, the k!-factors come from the Taylor-expansion and the factor $2k+1$ comes from integrating each term once.

We now consider the composite domain problem (3–9). To derive the weak formulation, one follows essentially the same steps as for the single domain problem, namely:

• 1. For both (3) and (4), multiply the equation for ${u}_{s}$ with a test function ${v}_{s} \in V_h^s$, and then integrate by parts;

• 2. Add terms consistent with the interface and boundary conditions; and

• 3. Add stabilization terms $j_1$ and $j_2$ over $\mathcal {F}_1$ and $\mathcal {F}_2$, respectively.

This results in the following weak formulation for (3–9). Find $u=\{ {u}_{1} , {u}_{2} \} \in V_h^1 \times V_h^2$ such that:

$$\begin{aligned} M(\dot{u},v)&+A(u,v)+a_\varGamma (u,v)+a_{\partial \varOmega }(u,v) \nonumber \\&=L_\varOmega (v)+L_\varGamma (v)+L_{\partial \varOmega }, \quad \forall v=\{ {v}_{1} , {v}_{2} \}\in V_h^1 \times V_h^2, \end{aligned}$$

(20)

where the bilinear forms M and A correspond to the stabilized mass and stiffness matrices:

$$\begin{aligned} M(\dot{u},v)&=\sum _{s=1}^2({\dot{u}}_{s},{v}_{s})_{\varOmega _s} + \gamma _M j_s({\dot{u}}_{s},{v}_{s}), \end{aligned}$$

(21)

$$\begin{aligned} A(u,v)&=\sum _{s=1}^2(\lambda \nabla {u}_{s}, \nabla {v}_{s})_{\varOmega _s} +\gamma _A h_T^{-2} \lambda j_s({u}_{s},{v}_{s}); \end{aligned}$$

(22)

$L_\varOmega$ corresponds to the forcing function:

$$\begin{aligned} L_\varOmega (v)&=\sum _{s=1}^2(f_s,{v}_{s})_{\varOmega _s}; \end{aligned}$$

(23)

$a_\varGamma$ and $L_\varGamma$ consistently enforce the interface conditions (8, 9):

$$\begin{aligned} a_\varGamma (u,v)&=-\left\langle [u], \{\lambda \frac{\partial v}{\partial n}\}\right\rangle _{\varGamma }-\left\langle \{\lambda \frac{\partial u}{\partial n}\},[v]\right\rangle _{\varGamma } + \left\langle \frac{\gamma _\varGamma }{h_T} [u],[v]\right\rangle _{\varGamma }, \end{aligned}$$

(24)

$$\begin{aligned} L_\varGamma (v)&=\left\langle \frac{\gamma _\varGamma }{h_T} \mu _1,[v]\right\rangle _\varGamma + \left\langle \kappa _1\mu _2,{v}_{2}\right\rangle _\varGamma + \left\langle \kappa _2\mu _2,{v}_{1}\right\rangle _\varGamma -\left\langle \mu _1,\{\lambda \frac{\partial v}{\partial n}\}\right\rangle _\varGamma ; \end{aligned}$$

(25)

and the terms $a_{\partial \varOmega }$ and $L_{\partial \varOmega }$ enforce the boundary condition (7) along the outer boundary, $\partial \varOmega$:

$$\begin{aligned} a_{\partial \varOmega }(u,v)&=-\left\langle \lambda \frac{\partial {u}_{1}}{\partial n},{v}_{1}\right\rangle _{\partial \varOmega } -\left\langle {u}_{1},\lambda \frac{\partial {v}_{1}}{\partial n}\right\rangle _{\partial \varOmega }+\frac{\gamma _{D}}{h_T} \left\langle \lambda {u}_{1},{v}_{1}\right\rangle _{\partial \varOmega }, \end{aligned}$$

(26)

$$\begin{aligned} L_{\partial \varOmega }(v)&=\left\langle \lambda \psi ,\frac{\gamma _D}{h_T} {v}_{1} - \frac{\partial {v}_{1}}{\partial n}\right\rangle _{\partial \varOmega }. \end{aligned}$$

(27)

In (24–27), n denotes the outward pointing normal at either $\varGamma$ or $\partial \varOmega$ (depending on the domain of integration); $\kappa _1 + \kappa _2 = 1$, so that $\{ v \} = \kappa _1 v_1+ \kappa _2 v_2$ is a convex combination; and $\gamma _\varGamma$, $\kappa _1$, $\kappa _2$ are chosen as in [13]:

$$\begin{aligned} \kappa _1 = \frac{{\lambda }_{2}}{{\lambda }_{1} + {\lambda }_{2}}, \quad \kappa _2 = \frac{{\lambda }_{1}}{{\lambda }_{1} + {\lambda }_{2}}, \quad \gamma _{\varGamma } = \gamma _D \frac{{\lambda }_{1}{\lambda }_{2}}{{\lambda }_{1} + {\lambda }_{2}}. \end{aligned}$$

(28)

The remaining parameters (appearing in Eqs. 21, 22, 26–28) are given by:

$$\begin{aligned} \gamma _M=0.75,\quad \gamma _A=1.5, \quad \gamma _D=5p^2. \end{aligned}$$

(29)

The scaling of $\gamma _D$ with respect to p follows from an inverse inequality. When $p=1$ these reduce to the same parameters as the ones used in [71], where $\gamma _M$ was chosen based on numerical experiments on the condition number of the mass matrix. This also agrees with the choice of $\gamma _A$ and $\gamma _D$ in [14], where $\gamma _A$ was investigated numerically.

In order to use cut-FEM, one needs a way to perform integration over the intersected elements $\mathcal {T}_\varGamma$. For example, with the interface problem, on each element $K \in \mathcal {T}_\varGamma$, we need a quadrature rule for the $K\cap \varOmega _1$, $K\cap \varOmega _2$ and $K\cap \varGamma$. For the numerical tests in this work (Sect. 4), we represent the geometry by a level set function, and compute high-order accurate quadrature rules with the algorithm from [68].

Remark 3

Optimal (second-order) convergence was rigorously proven for cut-FEM applied to the Poisson problem in [14]. As far as we know, there is no rigorous proof of higher-order convergence for cut-FEM, though such a proof would likely be similar to the second-order case.

3.2 DPM

We continue in this section with a brief introduction to the Difference Potentials Method (DPM), which was originally proposed by Ryaben’kii (see [64,65,66], and see [29, 33] for papers in his honor). Our aim is to consider the numerical approximation of PDEs on arbitrary, smooth geometries (defined either explicitly or implicitly) using the DPM together with standard, finite-difference discretizations of (1) or (3, 4) on uniform, Cartesian grids, which need not conform with boundaries or interfaces. To this end, we work with high-order methods for interface problems based on Difference Potentials, which were originally developed in [67] and [3,4,5,6, 26, 28]. We also introduce new developments here for handling implicitly-defined geometries. (The reader can consult [65] for the general theory of the Difference Potentials Method.)

Broadly, the main idea of the DPM is to reduce uniquely solvable and well-posed boundary value problems in a domain $\varOmega$ to pseudo-differential Boundary Equations with Projections (BEP) on the boundary of $\varOmega$. First, we introduce a computationally simple auxiliary domain as part of the method. The original domain is embedded into the auxiliary domain, which is then discretized using a uniform Cartesian grid. Next, we define a Difference Potentials operator via the solution of a simple Auxiliary Problem (defined on the auxilairy domain), and construct the discrete, pseudo-differential Boundary Equations with Projections (BEP) at grid points near the continuous boundary or interface $\varGamma$. (This set of grid points is called the discrete grid boundary.) Once constructed, the BEP are then solved together with the boundary/interface conditions to obtain the value of the solution at the discrete grid boundary. Lastly, using these reconstructed values of the solution at the discrete grid boundary, the approximation to the solution in the domain $\varOmega$ is obtained through the discrete, generalized Green’s formula.

Mathematically, the DPM is a discrete analog of the method of Calderón’s potentials in the theory of partial differential equations. The DPM, however, does not require explicit knowledge of Green’s functions. Although we use an Auxiliary Problem (AP) discretized by finite differences, the DPM is not limited to this choice of spatial discretization. Indeed, numerical methods based on the idea of Difference Potentials can be designed with whichever choice of spatial discretization is most natural for the problem at hand (e.g., see [25]).

Practically, the main computational complexity of the DPM reduces to the required solutions of the AP, which can be done very efficiently using fast, standard $\mathcal {O}(N \log N)$ solvers. Moreover, in general the DPM can be applied to problems with general boundary or interface conditions, with no change to the discretization of the PDE.

Let us now briefly introduce the DPM for the numerical approximation of parabolic interface models (3–9). First, we must introduce the point-sets that will be used throughout the DPM. (Note that the main construction of the method below applies to the single domain problem (1, 2), after omitting the index s and replacing interface conditions with boundary conditions; see [6].)

Let $\varOmega _s$ ($s=1,2$) be embedded in a rectangular auxiliary domain $\varOmega _{s}^0$. Introduce a uniform, Cartesian grid denoted $M_s^0$ on $\varOmega _s^0$, with grid-spacing $h_s$. Let $M_s^+ = M_s^0 \cap \varOmega _s$ denote the grid points inside each subdomain $\varOmega _s$, and $M_s^- = M_s^0 {\setminus } M_s^+$ the grid points outside each subdomain $\varOmega _s$. Note that the auxiliary domains $\varOmega _1^0$, $\varOmega _2^0$ and auxiliary grids $M_1^0$, $M_2^0$ need not agree, and indeed may be selected completely independently, given considerations regarding accuracy, adaptivity, or efficiency.

Define a finite-difference stencil $N_{j, k}^{s, \alpha }$, with $\alpha = 5, 9$, to be the stencil of the standard five-point or a wide nine-point Laplacian, i.e.,

$$\begin{aligned} N_{j, k}^{s, 5}&= \{ (x_j,y_k), (x_{j\pm 1},y_k), (x_j, y_{k\pm 1}) \} \quad \text {or}\nonumber \\ N_{j, k}^{s, 9}&= \{ (x_j,y_k), (x_{j\pm 1},y_k), (x_j, y_{k\pm 1}), (x_{j\pm 2},y_k), (x_j, y_{k\pm 2}) \}. \end{aligned}$$

(30)

Next, with $\alpha$ fixed, define the point-sets

$$\begin{aligned} N_s^0 = \mathop {\bigcup }\limits _{(x_j, y_k) \in M_s^0}\!\!\!\!\!\! N_{j, k}^{s, \alpha }, \quad N_s^+ = \mathop {\bigcup }\limits _{(x_j, y_k) \in M_s^+}\!\!\!\!\!\! N_{j, k}^{s, \alpha }, \quad \text {and} \quad N_s^- = \mathop {\bigcup }\limits _{(x_j, y_k) \in M_s^-} N_{j, k}^{s, \alpha }. \end{aligned}$$

(31)

The point-set $N_s^+$ ($N_s^-$, $N_s^0$) enlarges the point-set $M_s^+$ ($M_s^-$, $M_s^0$), by taking the union of finite-difference stencils at every point in $M_s^+$ ($M_s^-$, $M_s^0$). Therefore, $N_s^+$ contains a “thin row” of points belonging to the complement of $\varOmega _s$, so that $N_s^+ \not \subset \varOmega _s$, even though $M_s^+ \subset \varOmega _s$. (Likewise, $N_s^- \not \subset \varOmega {\setminus } \varOmega _s$, even though $M_s^- \subset \varOmega {\setminus } \varOmega _s$).

Lastly, we now define the important point-set

$$\begin{aligned} \gamma _s = N_s^+ \bigcap N_s^-, \end{aligned}$$

(32)

which we call the discrete grid boundary. In words, $\gamma _s$ is the set of grid points that straddle the continuous interface $\varGamma$. (See Fig. 4 for an example of these point-sets, given a single elliptical domain $\varOmega$.) Note that the point-sets $M_s^+$, $N_s^+$, and $\gamma _s$ will be used throughout the Difference Potentials Method.

Fig. 4

An example of the point-sets for the second-order Difference Potentials Method, applied to a single domain (s omitted for brevity), with a rotated ellipse $\varOmega$ and a rectangular auxiliary domain $\varOmega ^0$, showing (a) the physical grids $M^+$ (solid dots) and $N^+$ (open circles), with $M^+\subset N^+$, and (b) the discrete grid boundary $\gamma$ (open circles)

Here, we define the fully-discrete finite-difference discretization of (3, 4), and then define the Auxiliary Problem. Indeed, the discretization we consider is

$$\begin{aligned} L_{\Delta t, h}^s u_{s}^{i+1} = F_s^{i+1}, \quad (x_j, y_k) \in M_s^+, \end{aligned}$$

(33)

where (i) $L_{\Delta t, h}^s u_{s}^{i+1} := \lambda _s(t^{i+1}) \Delta _h u_{s}^{i+1} - \sigma u_{s}^{i+1}$, (ii) $\Delta _h$ is either a five- or nine-point Laplacian in each subdomain, and (iii) $\sigma$ and $F_s^{i+1}$ follow from the choice of time- and spatial-discretizations. (Here, we have simplified notation slightly by assuming that $h := h_1 = h_2$, which need not be the same in general.) For full details of the discretization, including the choice of BDF2 or BDF4 in the time discretization, we refer the reader to “Appendix 8.1”.

The choices of discretization (33) in each subdomain need not be the same. As in [3, 6], one could choose a second- and fourth-order discretization on $M_1^+$ and $M_2^+$, respectively, given considerations about accuracy, adaptivity, expected regularity of the analytical solution in each domain, etc.

Next, we define the discrete Auxiliary Problem, which plays a central role in the construction of the Difference Potentials operator, the resulting Boundary Equations with Projection at the discrete grid boundary, and in the numerical approximation of the solution via the discrete, generalized Green’s formula.

Definition 1

(Discrete Auxiliary Problem (AP)) At time $t^{i+1}$, given the right-hand side grid function $q_s^{i+1}$: $M_s^0 \rightarrow \mathbb {R}$, the following difference Eqs. (34, 35) are defined as the discrete AP.

$$\begin{aligned} L_{\Delta t, h}^s u_s^{i+1}= & {} q_s^{i+1}, \quad (x_j,y_k)\in M^0_s \end{aligned}$$

(34)

$$\begin{aligned} u_s^{i+1}= & {} 0, \quad (x_j,y_k)\in N^0_s\backslash M^0_s \end{aligned}$$

(35)

Remark 4

For a given right-hand side $q_s^{i + 1}$, the solution of the discrete AP (34, 35) defines a discrete Green’s operator $G_{\Delta t, h}^s q_s^{i + 1}$. The choice of boundary conditions (35) will affect the resulting grid function $G_{\Delta t, h}^s q_s^{i + 1}$, and thus the Boundary Equations with Projection defined below. However, the choice of boundary conditions (35) in the AP will not affect the numerical approximation of (3–9), so long as the discrete AP is uniquely solvable and well-posed.

Let us denote by $G_{\Delta t, h}^s F_s^{i + 1}$ the particular solution on $N^+_s$ of the fully-discrete problem (33), defined by solving the AP (34, 35) with

$$\begin{aligned} q_s^{i + 1} = {\left\{ \begin{array}{ll} F_s^{i + 1}, &{} (x_j, y_k) \in M_s^+, \\ 0, &{} (x_j, y_k) \in M_s^-, \end{array}\right. } \end{aligned}$$

(36)

and restricting the solution from $N_s^0$ to $N_s^+$.

Let us also introduce a linear space $\mathbf {V}_{\gamma _s}$ of all grid functions denoted $v_{\gamma _s}^{i + 1}$, which are defined on $\gamma _s$ and extended by zero to the other points of $N_s^0$. These grid functions are referred to as discrete densities on $\gamma _s$.

Definition 2

(The Difference Potential of a density.) The Difference Potential of a given density $v_{\gamma _s}^{i + 1}$ is the grid function $P_{N_s^+}^{i + 1} v_{\gamma _s}^{i + 1}$ on $N_s^+$, defined by solving the AP (34, 35) with

$$\begin{aligned} q_s^{i + 1} = {\left\{ \begin{array}{ll} 0, &{} (x_j, y_k) \in M_s^-, \\ L_{\Delta t, h}^s [v_{\gamma _s}^{i + 1}], &{} (x_j, y_k) \in M_s^+, \end{array}\right. } \end{aligned}$$

(37)

and restricting the solution from $N_s^0$ to the point-set $N_s^+$.

Note that $P_{N_s^+}^{i + 1} : \mathbf {V}_{\gamma _s} \rightarrow N_s^+$ is a linear operator on the space $\mathbf {V}_{\gamma _s}$ of densities. Moreover, the coefficients of $P_{N_s^+}^{i + 1}$ can be computed by solving the AP (Definition 1) with the appropriate density $v_{\gamma _s}^{i+1}$ defined at the points $(x_{j}, y_{k}) \in \gamma _s$.

Definition 3

(The Trace operator.) Given a grid function $v_s^{i + 1}$, we denote by ${\text {Tr}}_{\gamma _s} [v_s^{i + 1}]$ the Trace (or Restriction) from $N_s^+$ to $\gamma _s$.

Moreover, for a given density $v_{\gamma _s}^{i + 1}$, denote the trace of the Difference Potential of $v_{\gamma _s}^{i + 1}$ by $P_{\gamma _s} v_{\gamma _s}^{i + 1}$. In other words, $P_{\gamma _s} v_{\gamma _s}^{i + 1} = {\text {Tr}}_{\gamma _s} [P_{N_s^+}^{i + 1} v_{\gamma _s}^{i + 1}]$.

Now we can state the central theorem of the Difference Potentials Method that will allow us to reformulate the finite-difference Eq. (33) on $M_s^+$ (without imposing any boundary or interface conditions yet) into the equivalent Boundary Equations with Projections on $\gamma _s$.

Theorem 1

(Boundary Equations with Projection (BEP)) At time-level $t^{i+1}$, the discrete density $u^{i+1}_{\gamma _s}$ ($s=1,2$) is the trace of some solution $u^{i+1}_s$ on domain $\varOmega _s$ to the Difference Eq. (33), i.e., $u^{i+1}_{\gamma _s}:={{\text {Tr}}_{\gamma _s}}[u^{i+1}_s]$, if and only if the following BEP holds

$$\begin{aligned} u^{i+1}_{\gamma _s}-P_{\gamma _s}^{i+1}u^{i+1}_{\gamma _s}={{\text {Tr}}_{\gamma _s}}[G^{i+1}_{\Delta t, h}F^{i+1}_s], \quad (x_j,y_k)\in \gamma _s, \end{aligned}$$

(38)

with ${\text {Tr}}_{\gamma _s}[\cdot ]$ and $P_{\gamma _s}$ defined in Definition 3.

Proof

See [65] for the general theory of DPM (including the proof for general elliptic PDE), or one of [3, 5, 6] for the proof in the case of parabolic interface problems. $\square$

Remark 5

A given density $v_{\gamma _s}^{i+1}$ is the trace of some solution of the fully-discrete finite-difference Eq. (33) if and only if it is a solution of the BEP.

However, since boundary or interface conditions have not yet been imposed, the BEP will have infinitely many solutions $u_{\gamma _s}^{i+1}$. As originally disucssed in [3,4,5,6, 26, 28], in this work we consider the following approach in order to find a unique solution of the BEP.

At each time level $t^{i + 1}$, one can approximate the solution of (3–9) at the discrete grid boundary $\gamma _s$, using the Cauchy data of (3–9) on the continuous interface $\varGamma$, up to the desired second- or fourth-order accuracy. (By Cauchy data, we mean the trace of the solution of (3–9), together with the trace of its normal derivative, on $\varGamma$.) Below, we will define an Extension Operator which will extend the Cauchy data of (3–9) from $\varGamma$ to $\gamma _s$.

As we will see, the Extension Operator in this work depends only on the given parabolic interface model. Moreover, we will use a finite-dimensional, spectral representation for the Cauchy data of (3–9) on $\varGamma$. Then, we will use the Extension Operator, together with the BEP (38) and the interface conditions (8, 9), to obtain a linear system of equations for the coefficients of the finite-dimensional, spectral representation. Hence, the derived BEP will be solved for the unknown coefficients of the Cauchy data. Using this obtained Cauchy data, we will construct the approximation of (3–9) using the Extension Operator, together with the discrete, generalized Green’s formula.

Let us now briefly discuss the Extension Operator for the second-order numerical method, and refer the reader to “Appendix 8.2” for details (including details for the fourth-order numerical method). For points in the vicinity of $\varGamma$, we define a coordinate system $(d, \vartheta )$, where $\vartheta$ is arclength from some reference point, and d is the signed distance in the normal direction from the point to $\varGamma$. Now, as a first step towards defining the Extension Operator, we define a new function

$$\begin{aligned} v_s^{i+1}(d, \vartheta ) = v_s^{i+1}(0, \vartheta ) + \sum _{l = 1}^p \frac{1}{l!} \frac{\partial ^l v_s^{i+1}(0, \vartheta )}{\partial n^l} d^l, \end{aligned}$$

(39)

where n is the unit outward normal vector at $\varGamma$. We choose $p = 2$ for the second-order method (which we will discuss now) and $p = 4$ for the fourth-order method (see “Appendix 8.2”).

As a next step for the second-order method (BDF2–DPM2), we define

$$\begin{aligned} v_s^{i+1}(0, \vartheta ) = u_s^{i + 1}|_\varGamma , \quad \frac{\partial v_s^{i+1}(0, \vartheta )}{\partial n} = \left. \frac{\partial u_s^{i + 1}}{\partial n} \right| _\varGamma , \quad \text {and} \quad \frac{\partial ^2 v_s^{i+1}(0, \vartheta )}{\partial n^2} = \left. \frac{\partial ^2 u_s^{i + 1}}{\partial n^2} \right| _\varGamma ,\nonumber \\ \end{aligned}$$

(40)

where $u_s^{i + 1} := u_s(x, y, t^{i + 1})$, $\frac{\partial u_s^{i + 1}}{\partial n} := \frac{\partial u_s(x, y, t^{i + 1})}{\partial n}$, etc. As a last step, a straightforward sequence of calculations (see “Appendix 8.2”) shows that

$$\begin{aligned} \frac{\partial ^2 u_s^{i + 1}}{\partial n^2} \approx \frac{1}{\lambda _s} \left( \frac{3 u_s^{i + 1} - 4 u_s^i + u_s^{i - 1}}{2 \Delta t} - f^{i + 1} \right) - \frac{\partial ^2 u_s^{i + 1}}{\partial \vartheta ^2} + \kappa \frac{\partial u_s^{i+1}}{\partial n}, \end{aligned}$$

(41)

where $\kappa$ denotes the curvature of $\varGamma$. Therefore, with $v_s^{i+1}(d, \vartheta )$ defined by (39–41), the only unknown data at each time step $t^{i + 1}$ are the unknown Dirichlet data $u_s^{i + 1}$ and the unknown Neumann data $\frac{\partial u_s^{i + 1}}{\partial n}$. The Extension Operator will incorporate the interface conditions (8, 9) when it is combined with the BEP (38), so that the only independent unknowns at each time step $t^{i + 1}$ will be $\left. \left( u_1^{i + 1}, \frac{\partial u_1^{i + 1}}{\partial n} \right) \right| _\varGamma$ or $\left. \left( u_2^{i + 1}, \frac{\partial u_2^{i + 1}}{\partial n} \right) \right| _\varGamma$. (This is also true for the fourth-order numerical method—see Appendices 8.2 and 8.3.)

Now we are ready to define the Extension Operator that extends the Cauchy data of (3–9) from $\varGamma$ to $\gamma _s$.

Definition 4

(The Extension Operator.) Let $v_s^{i + 1}(d, \vartheta )$ be defined by (39–41). Let $\mathfrak {u}_{s, \varGamma }^{i + 1}$ denote the Cauchy data of (3–9) at $t^{i + 1}$ on $\varGamma$, i.e., $\mathfrak {u}_{s, \varGamma }^{i + 1} = \left. \left( u_s^{i + 1}, \frac{\partial u_s^{i + 1}}{\partial n} \right) \right| _\varGamma$. The extension operator ${\text {Ex}}_s$ that extends $\mathfrak {u}_{s, \varGamma }^{i + 1}$ from $\varGamma$ to $\gamma _s$ is

$$\begin{aligned} {\text {Ex}}_s \mathfrak {u}_{s, \varGamma }^{i + 1} := v_s^{i + 1}(d, \vartheta )|_{\gamma _s}. \end{aligned}$$

(42)

For a given point $(x_j, y_k) \in \gamma _s$, note that d is the signed distance between $(x_j, y_k)$ and its orthogonal projection on $\varGamma$, while $\vartheta$ is the arclength along $\varGamma$ between a reference point and the orthogonal projection of $(x_j, y_k)$.

Next, we briefly discuss the finite-dimensional, spectral representation of Cauchy data $\mathfrak {u}_{s, \varGamma }^{i + 1}$. Indeed, we wish to choose a basis $\phi _{\nu }(\vartheta )$ on $\varGamma$ ($\nu = 1, 2, 3, \ldots$) in order to accurately approximate the two components of the Cauchy data $\mathfrak {u}_{s, \varGamma }^{i + 1}$. To be specific, whichever basis we choose, we require that

$$\begin{aligned}&\varepsilon _{\mathcal {N}^0, \mathcal {N}^1}(\mathfrak {u}_{s, \varGamma }^{i + 1})\nonumber \\&\quad ={\text {*}}{min}_{c_{1, \nu }^{s, i+1}, c_{2, \nu }^{s, i+1}} \int _{\varGamma } \Bigg (\Big | u_s^{i + 1} - \sum _{\nu =1}^{\mathcal {N}^0} c_{1, \nu }^{s, i+1} \phi _{\nu }(\vartheta )\Big |^2 + \Big | \frac{\partial u_s^{i + 1}}{\partial n} -\sum _{\nu =1}^{\mathcal {N}^1} c_{2, \nu }^{s, i+1} \phi _{\nu }(\vartheta )\Big |^2\Bigg ) \mathrm{d\vartheta } \quad \end{aligned}$$

(43)

tends to zero as $\mathcal {N}^0, \mathcal {N}^1 \rightarrow \infty$, for some sequence of real numbers $(c_{1, \nu }^{s, i+1})_{\nu =1}^{\mathcal {N}^0}$ and $(c_{2, \nu }^{s, i+1})_{\nu =1}^{\mathcal {N}^1}$. In other words, we require

$$\begin{aligned} \lim _{\mathcal {N}^0, \mathcal {N}^1 \rightarrow \infty } \varepsilon _{\mathcal {N}^0, \mathcal {N}^1}(\mathfrak {u}_{s, \varGamma }^{i + 1}) = 0. \end{aligned}$$

(44)

Now let us discuss a choice of basis. In this work, recall that we consider interfaces $\varGamma$ that are at least $C^2(\varGamma )$ (due to the choice of smooth, circular geometries). Also, as we will see in Sect. 4.1, each function u considered in the test problems on a composite domain (TP–2A, TP–2B, TP–2C) is locally smooth, in the sense that $u|_{\varOmega _1} = u_1$ and $u|_{\varOmega _2} = u_2$ are smooth in $\varOmega _1$ and $\varOmega _2$, respectively. Moreover, each component of the Cauchy data $\mathfrak {u}_{1, \varGamma }^{i + 1}$ and $\mathfrak {u}_{2, \varGamma }^{i + 1}$ are smooth, periodic functions of arclength $\vartheta$. (Note that $\mathfrak {u}_{1, \varGamma }^{i + 1}$ and $\mathfrak {u}_{2, \varGamma }^{i + 1}$ need not agree, and indeed do not—neither $\mu _1(x, y, t)$ nor $\mu _2(x, y, t)$ in (8, 9) are identically equal to zero, for any of our test problems on a composite domain.)

Therefore, in this work, we choose a standard trigonometric basis $\phi _\nu (\vartheta )$, with

$$\begin{aligned} \phi _1(\vartheta ) = 1, \quad \phi _{2K}(\vartheta ) = \cos \left( \frac{2 \pi K \vartheta }{|\varGamma |}\right) , \quad \phi _{2K + 1}(\vartheta ) = \sin \left( \frac{2 \pi K \vartheta }{|\varGamma |}\right) , \end{aligned}$$

(45)

and $K > 1$. Moreover, at every time step $t^{i + 1}$, we will discretize the Cauchy data $\mathfrak {u}_{s, \varGamma }^{i + 1} = \left. \left( u_s^{i + 1}, \frac{\partial u_s^{i + 1}}{\partial n} \right) \right| _\varGamma$ using this basis. Therefore, we let

$$\begin{aligned} \tilde{\mathfrak {u}}_{s, \varGamma }^{i+1} = \sum _{\nu =1}^{\mathcal {N}^0} c^{s, i+1}_{1,\nu } \Phi ^0_{\nu }(\vartheta ) + \sum _{\nu =1}^{\mathcal {N}^1} c^{s, i+1}_{2,\nu } \Phi ^1_{\nu }(\vartheta ) \quad \text {and} \quad \tilde{\mathrm{u}}^{i+1}_{s,\varGamma } \approx \mathfrak {u}^{i+1}_{s,\varGamma }, \end{aligned}$$

(46)

where $\Phi ^0_{\nu } = (\phi _{\nu },0)$ and $\Phi ^1_{\nu } = (0, \phi _{\nu })$ are the set of basis functions used to represent the Cauchy data on the interface $\varGamma$.

Remark 6

It should be also possible to relax regularity assumption on the domain under consideration. For example, one can consider piecewise-smooth, locally-supported basis functions (defined on $\varGamma$) as the part of the Extension Operator. For example, [52] use this approach to design a high-order accurate numerical method for the Helmholtz equation, in a geometry with a reentrant corner. Furthermore, [80, 81] combine the DPM together with the XFEM, and design a DPM for linear elasticity in a non-Lipschitz domain (with a cut).

Next, in “Appendix 8.3”, we derive a linear system for the coefficients $(c_{1, \nu }^{s, i+1})_{\nu =1}^{\mathcal {N}^0}$ and $(c_{2, \nu }^{s, i+1})_{\nu =1}^{\mathcal {N}^1}$, by combining the interface conditions (8, 9), the BEP (38), the Extension Operator (39–41), and the spectral discretization (46). Then, the numerical approximation $u_s^{i+1} \approx u_s(x_j, y_k, t^{i + 1})$ of (3–9) at all grid-points $(x_j, y_k) \in N_s^+$ follows directly from the discrete, generalized Green’s formula, which we state now.

Definition 5

(Discrete, generalized Green’s formula.) At each time step $t^{i + 1}$, the numerical approximation $u_s^{i+1} \approx u_s(x_j, y_k, t^{i + 1})|_{(x_j, y_k) \in N_s^+}$ of (3–9) is given by

$$\begin{aligned} u_s^{i + 1} := P_{N_s^+}^{i + 1} u_{\gamma _s}^{i + 1} + G_{\Delta t, h}^s F_s^{i + 1}. \end{aligned}$$

(47)

Here, $u_{\gamma _s}^{i + 1} = {\text {Ex}}_s \tilde{\mathfrak {u}}_{s, \varGamma }^{i + 1}$, and $\tilde{\mathfrak {u}}_{s, \varGamma }^{i + 1}$ is constructed from (i) the solution of the BEP (see “Appendix 8.3”) and (ii) the spectral discretization (46). (Recall that $P_{N_s^+}^{i + 1} u_{\gamma _s}^{i + 1}$ is the Difference Potential of the density $u_{\gamma _s}^{i + 1}$, while $G_{\Delta t, h}^s F_s^{i + 1}$ is the Particular Solution.)

In this work, we also propose a novel feature of DPM, extending the method originally developed in [67] and [3,4,5,6, 26, 28] to the composite domain problem (3–9) with implicitly-defined geometry. The primary difference between Difference Potentials Methods on explicitly-defined versus implicitly-defined composite domains is in the approximation of the interface $\varGamma$, which must be done accurately and efficiently, in order to maintain the desired second- or fourth-order accuracy.

The main idea of DPM-based methods for implicitly-defined geometry is to seek an accurate and efficient explicit parameterization of the implicit boundary/interface. First, we represent the geometry implicitly via a level set function F(x, y) on $M^0$. Then we construct a local interpolant $\tilde{F}(x, y)$ of F(x, y) on a subset of $M^0$ near the continuous interface $\varGamma$. Next, we parameterize $\varGamma$ by arc-length using numerical quadrature. With this parameterization, we (i) compute the Fourier series expansion from initial conditions for the Cauchy data $\mathfrak {u}_{s, \varGamma }^{i + 1}$ on the implicit interface $\varGamma$, and (ii) construct the extension operators (Definition 4) with $p=2$ or $p=4$.

Conjecture 1

(High-order accuracy of the DPM with implicit geometry) Due to the second- or fourth-order accuracy (in both space and time) of the underlying discretization (33), the extension operator (42) with $p=2$ or $p=4$, and the established error estimates and convergence results for the DPM for general linear elliptic boundary value problems on smooth domains (presented in [33, 62, 63, 65]), we expect second- and fourth-order accuracy in the maximum norm for the error in the computed solution (59 or 60) for both the single and composite domain parabolic problems.

Remark 7

Indeed, in the numerical results (Sect. 5) we see that the computed solution (47) at every time level $t^{i+1}$ has accuracy $\mathcal {O}(h^2 + \Delta t^2)$ for the second-order method, and $\mathcal {O}(h^4 + \Delta t^4)$ for the fourth-order method, for both the single and composite domain problems, with explicit or implicit geometry. See [3, 6, 27, 67] for more details and numerical tests involving explicit (circular and elliptical) geometries.

Main Steps of the algorithm Let us summarize the main steps for the Difference Potentials Method.

• Step 1 Introduce a computationally simple Auxiliary Domain $\varOmega _s^0$ ($s = 1,2$) and formulate the Auxiliary Problem (AP; Definition 1).

• Step 2 At each time step $t^{i+1}$, compute the Particular Solution $u_{s}^{i+1} = G^{i+1}_{\Delta t, h} F^{i+1}_{s}$, $(x_j,y_k)\in N_s^+$, using the AP with the right-hand side (36).

• Step 3 Construct the matrix in the boundary Eq. (79) (discussed in “Appendix 8.3”), derived from the Boundary Equation with Projection (BEP) (38), via several solutions of the AP. (When the diffusion coefficients $\lambda _s$ are constant, this is done once, as a pre-processing step before the first time step.)

• Step 4 Solve boundary Eq. (79) and compute the approximation of the density $u^{i+1}_{\gamma _s}$, by applying the Extension Operator (42) to the solution of (79).

• Step 5 Construct the Difference Potentials $P_{N_s^+\gamma _s} u^{i+1}_{\gamma _s}$ of the density $u_{\gamma _s}^{i + 1}$, using the AP with the right-hand side (37).

• Step 6 Compute the numerical approximation $u_s^{i+1} \approx u_s(x_j, y_k, t^{i + 1})$ of the PDE (3–9) using the discrete, generalized Green’s formula (47).

3.3 SBP–SAT–FD

We continue in this section with a brief presentation of SBP–SAT–FD, for solving the parabolic problems presented in Sect. 2. For more detailed discussions of the SBP–SAT–FD method, we refer the reader to two review papers [21, 73].

The SBP–SAT–FD method was originally used on Cartesian grids. To resolve complex geometries, we consider a grid mapping approach by transfinite interpolation [43]. A smooth mapping requires that the physical domain is a quadrilateral, possibly with smooth, curved sides. If the physical domain does not have the desired shape, we then partition the physical domain into subdomains, so that each subdomain can be mapped smoothly to the reference domain. As an example, the single domain of equation (1, 2), shown in Figure 5a, is divided into five subdomains. The five subdomains consist of one square subdomain, and four identical quadrilateral subdomains (modulo rotation by $\pi /2$) with curved sides. Similarly, the composite domain of equation (3–9) is divided into nine subdomains, as shown in Fig. 5b. Suitable interface conditions are imposed to patch the subdomains together.

Although the side-length of the centered square is arbitrary (as long as the square is strictly inside the circle), its size and position have a significant impact on the quality of the curvilinear grid. In a high-quality mesh, the elements should not be skewed too much, and the sizes of the elements should be nearly uniform. In practice, it is usually difficult to know a priori the optimal way of domain division.

A Cartesian grid in the reference domain is mapped to a curvilinear grid in each subdomain. The grids are aligned with boundaries and interfaces, thus avoiding small-cut difficulties sometimes associated with embedded methods. In this paper, we only consider conforming grid interfaces, i.e., the grid points from two adjacent blocks match on the interface. For numerical treatment of non-conforming grid interfaces in the SBP–SAT–FD framework, see [44, 55].

Fig. 5

A (a) circular domain (cf. Fig. 1), divided into five subdomains; and (b) composite domain, divided into nine subdomains

When a physical domain is mapped to a reference domain, the governing equation is transformed to the Cartesian coordinate in the reference domain. The transformed equation is usually in a more complicated form than the original equation. In general, a parabolic problem

$$\begin{aligned} u_t = u_{xx}+u_{yy},\quad (x,y)\in \varOmega \end{aligned}$$

(48)

in a physical domain will be transformed to

$$\begin{aligned} Ju_t = (\alpha u_\xi )_\xi +(\beta u_\eta )_\xi +(\beta u_\xi )_\eta +(\gamma u_\eta )_\eta ,\quad (\xi ,\eta )\in [0,1]^2, \end{aligned}$$

(49)

where $(\xi ,\eta )$ is the Cartesian coordinate in the unit square, and $J(\xi ,\eta )$, $\alpha (\xi ,\eta )$, $\beta (\xi ,\eta )$, $\gamma (\xi ,\eta )$ depend on the geometry of the physical domain and on the chosen mapping. In particular, we use transfinite interpolation for the grid mapping. In this case, the precise form of (49) and the derivation of the grid transformation are presented in Section 3.2 of [7]. Even though the original equation is in the simplest form with unit coefficients, the transformed equation has variable coefficients and mixed derivatives. Therefore, it is important to construct multi-block finite difference methods solving the transformed equation (49). Hence, we need two SBP operators, $D_1\approx \partial /\partial x$ to approximate a first derivative, and $D^{(b)}_2\approx \partial /\partial x(b(x)\partial /\partial x)$ to approximate a second derivative with variable coefficient, where $b(x)>0$ is a known function. Below we discuss SBP properties, and start with the first derivative.

Consider two smooth functions u(x), v(x) on $x\in [0,1]$. We discretize [0, 1] uniformly by N grid points, and denote the restriction of u(x), v(x) onto the grid by $\mathbf {u},\mathbf {v}$, respectively. Integration by parts states:

$$\begin{aligned} \int _0^{1} u_xv\ dx= uv\Big |_0^1 - \int _0^{1} uv_x\ dx. \end{aligned}$$

(50)

The SBP operator $D_1$ mimics integration by parts:

$$\begin{aligned} (D_1 \mathbf {u})^T H \mathbf {v}= \mathbf {u}^T B\mathbf {v} - \mathbf {u}^T H D_1\mathbf {v}, \end{aligned}$$

(51)

where H is symmetric positive definite—thus defining an inner product—and

$$\begin{aligned} B=\text {diag}(-1,0,\cdots ,0,1). \end{aligned}$$

In fact, H is also a quadrature [20]. It is easy to verify that (51) is equivalent to

$$\begin{aligned} D_1^TH+HD_1=B, \end{aligned}$$

(52)

which is the SBP property for the first derivative operator. At the grid points in the interior of the domain, standard, central, finite-difference stencils can be used in $D_1$, and the weights of the standard, discrete $L_2$-norm are used in H. At a few points close to boundaries, special stencils and weights must be constructed in $D_1$ and H, respectively, to satisfy (52).

The SBP operators $D_1$ were first constructed in [45] and later revisited in [72]. The SBP norm H can be diagonal or non-diagonal. While non-diagonal norm SBP operators have a better accuracy property than diagonal norm SBP operators, when terms with variable coefficients are present in the equation, a stability proof is only possible with diagonal norm SBP operators. Therefore, we use diagonal norm SBP operators in this paper.

For a second derivative with variable coefficients, the SBP operators $D^{(b)}_2$ were constructed in [54]. We remark that applying $D_1$ twice also approximates a second derivative, but is less accurate and more computationally expensive than $D^{(b)}_2$.

Due to the choice of centered difference stencils at interior grid points, the order of accuracy of the SBP operators is even at these points, and is often denoted by 2p. To fulfill the SBP property, at a few grid points near boundaries, the order of accuracy is reduced to p for diagonal norm operators. This detail notwithstanding, such a scheme is often referred to as $2p^{th}$-order accurate. In fact, for the second- and fourth-order SBP–SAT–FD schemes used in this paper to solve parabolic problems, we can expect a second- and fourth-order overall convergence rate, respectively [77].

An SBP operator only approximates a derivative. When imposing boundary and interface conditions, it is important that the SBP property is preserved and an energy estimate is obtained. For this reason, we consider the SAT method [16], where penalty terms are added to the semi-discretization, imposing the boundary and interface conditions weakly. This bears similarities with the Nitsche finite element method [60] and the discontinuous Galerkin method [40].

We note that in [75], SBP–SAT–FD methods were developed for the wave equation

$$\begin{aligned} Jv_{tt} = (av_\xi )_\xi +(bv_\eta )_\xi +(bv_\xi )_\eta +(cv_\eta )_\eta ,\quad (\xi ,\eta )\in [0,1]^2, \end{aligned}$$

(53)

with Dirichlet boundary conditions, Neumann boundary conditions, and interface conditions. Comparing Eq. (53) with (49), the only difference is that the wave equation has a second derivative in time, while the heat equation has a first derivative in time. The spatial derivatives of (53) and (49) are the same.

Assuming homogeneous boundary data for simplified notation, we write the SBP–SAT–FD discretization of (53) as

$$\begin{aligned} \mathbf {v}_{tt} = Q\mathbf {v}, \end{aligned}$$

(54)

where Q is the spatial discretization operator including the boundary implementation. For the scheme developed in [75], stability is proved by the energy method by multiplying (54) by $\mathbf {v}_t^TH_2$ from the left,

$$\begin{aligned} \mathbf {v}^T_tH_2\mathbf {v}_{tt} = \mathbf {v}^T_tH_2Q\mathbf {v}, \end{aligned}$$

(55)

where $H_2$ is a diagonal, positive-definite operator, obtained through a tensor product from the corresponding SBP norm, H, in one spatial dimension. It is shown in [75] that $H_2Q$ is symmetric and negative semi-definite. Therefore, we can write (55) as

$$\begin{aligned} \frac{d}{d t} (\mathbf {v}^T_tH_2\mathbf {v}_t - \mathbf {v}^TH_2Q\mathbf {v})=0, \end{aligned}$$

where the discrete energy, $\mathbf {v}^T_tH_2\mathbf {v}_t - \mathbf {v}^TH_2Q\mathbf {v}$, for (53) is conserved.

If we use the same operator Q to discretize the heat equation (49) with the same boundary condition as the wave equation (53), then the scheme

$$\begin{aligned} \mathbf {v}_t = Q\mathbf {v}, \end{aligned}$$

(56)

is also stable. To see this, we multiply (56) by $\mathbf {v}^TH_2$ from the left, and obtain

$$\begin{aligned} \frac{d}{d t} (\mathbf {v}^TH_2\mathbf {v}) = \mathbf {v}^TH_2Q\mathbf {v} \le 0, \end{aligned}$$

(57)

where $\mathbf {v}^TH_2\mathbf {v}$ is the discrete energy for (49). In this paper, we use the spatial discretization operators developed in [75] to solve both the single (1, 2) and composite domain problems (3–9).

In [10], SBP–SAT–FD methods are discussed for the one-dimensional heat equation with constant coefficients, both in a single domain and a composite domain. In theory, these schemes can also be generalized to solve Eq. (49), but are different from the ones used in this paper.

4 Test Problems

In this section, we first list the test problems that we will consider (in Sect. 4.1), and then briefly motivate and discuss these choices (in Sect. 4.2). The tests we propose are “manufactured solutions”, in the sense that we state an exact solution u(x, y, t) or $(u_1(x, y, t), \, u_2(x, y, t))$ and a diffusion coefficient $\lambda (t)$ or $(\lambda _1, \lambda _2)$. From (1, 2) (for the single domain problem) or (3–9) (for the composite domain problem) we compute the (i) right-hand side, (ii) initial conditions, (iii) boundary condition, and (iv) functions $(\mu _1(x, y, t), \, \mu _2(x, y, t))$ for the interface/matching conditions. Then, (i–iv), together with the diffusion coefficient, serve as the inputs for our numerical methods.

4.1 List of Test Problems

• 1. Single-domain, with an explicitly-defined boundary for DPM and SBP–SAT–FD, or an implicitly-defined boundary for cut-FEM.

  • (a) Constant diffusion (Test Problem 1A; TP–1A): Consider the PDE (1, 2), with $\lambda (t) \equiv 1$, $\varOmega = \{ (x, y) \in \mathbb {R}^2 : x^2 + y^2 \le 1 \}$, and the final time $T=1.0$. Then, TP–1A (adapted from [6]), is given by

    

  • (b) Time-varying diffusion (Test Problem 3A; TP–3A): Same as TP–1A, but with diffusion coefficient

    

• 2. Composite-domain, with an explicitly-defined interface (for DPM and SBP–SAT–FD) or implicitly-defined interface (for cut-FEM and DPM). Consider the PDE (3–9), with $\varOmega = [-2, 2] \times [-2, 2]$, $\varOmega _2 = \{ (x, y) \in \mathbb {R}^2 : x^2 + y^2 \le 1 \}$, $\varOmega _1 = \varOmega {\setminus } \varOmega _2$, $\varGamma = \{ (x, y) \in \mathbb {R}^2 : x^2 + y^2 = 1 \}$, and the final time $T = 1.0$.

  • (a) (Test Problem 2A; TP–2A): A modified version of the test adapted from [6, 48]. Let $(\lambda _1, \lambda _2) = (10, 1)$, and

    

  • (b) High-frequency oscillations (Test Problem 2B; TP–2B): A modified version of the test adapted from [6]. Let $(\lambda _1, \lambda _2) = (10, 1)$, and

    

  • (c) Large contrast in diffusion coefficients, and large jumps in both solution and flux at interface (Test Problem 2C; TP–2C): Let $(\lambda _1, \lambda _2) = (1000, 1)$, and

    

4.2 Motivation of the Chosen Test Problems

Test Problem 1A (TP–1A) involves a high-degree polynomial, with total degree of 17. This is a rather straightforward test problem, which allows us to establish a good “baseline” with which to compare each method. The choice of high degree ensures that there will be no cancellation of local truncation error, so that we should see—at most—second- or fourth-order convergence for the given methods, barring some type of superconvergence. Next, (TP–3A) adds on (incrementally) the complication of time-varying diffusion.

Likewise, (TP–2A) offers a straightforward “baseline” with which to consider the interface problem: The test problem is piecewise-smooth, and the geometry is simplified (see Remark 1). However, there is a jump in both the analytical solution and its flux, which requires a well-designed numerical method to accurately approximate. Moreover, (TP–2A) was first proposed in [48] (see also [6]), and is a good comparison with the immersed interface method therein.

Then, (TP–2B) adds additional challenges onto (TP–2A) in the form of much higher-frequency oscillations; while (TP–2C) adds onto (TP–2A) in the form of both (i) large contrast in diffusion, and (ii) large jumps in the analytical solution and its flux.

5 Numerical Results

5.1 Time Discretization

The spatial discretization for each method is discussed in Sect. 3. For the time discretization, the backward differentiation formulas of second- and fourth-order (BDF2 and BDF4) are used for the second- and fourth-order methods, respectively. In each case, the time-step is given by

$$\begin{aligned} \Delta t=0.5h. \end{aligned}$$

(58)

However, note that h in (58) bears different physical meanings for each method. Indeed, for cut-FEM, h is the average distance between the Gauss–Lobatto points; for DPM, h is the grid spacing in the uniform, Cartesian grid $M^0$ (see the text prior to (33)); and for SBP–SAT–FD, h is the minimum grid spacing in the reference domain.

5.2 Measure for Comparison

Let $\smash {u_{j, k}^i}$ denote the computed numerical approximation of u(x, y, t) at the grid-point $(x_j, y_k) \in \varOmega$ and time $t^i = i \Delta t \in (0, T]$. For the three methods, we will compare the size of the maximum error in u at the grid points, with respect to the number of degrees of freedom (DOF). For the single domain problem (1, 2), the maximum error is computed as:

$$\begin{aligned} E := \max _{t^i \in (0, T]} \max _{(x_j, y_k) \in \varOmega } |u(x_j, y_k, t^i) - u_{j, k}^i|, \end{aligned}$$

(59)

and for the composite domain problem (3–9) as:

$$\begin{aligned} E&:= \max _{t^i \in (0, T]} \max _{(x_j, y_k) \in \varOmega _1 \cup \varOmega _2} |u(x_j, y_k, t^i) - u_{j, k}^i|. \end{aligned}$$

(60)

5.3 Convergence Results

In the following tables and figures, we state the number of degrees of freedom in the grid, maximum error (59, 60 for the single- and composite-domain problems, respectively), and an estimate of the rate of convergence.

In Tables 1, 2, 3, 4 and 5, the estimate of rate of convergence is computed as follows. Let $(\text {DOF}_{n}, \, E_n)$ be given, with $n = 1, 2, 3$ referring to the first, second, and third grids (from coarsest to finest). Then, for $n = 2, 3$, compute the standard estimate

$$\begin{aligned} \rho _n = \frac{\log (E_{n - 1} / E_n)}{\log (\text {DOF}_{n - 1} / \text {DOF}_{n})}, \end{aligned}$$

(61)

which is the estimated rate of convergence, denoted in Tables 1, 2, 3, 4 and 5 by “Rate”.

In Figs. 6, 7, 10, 11 and 12, the estimate of rate of convergence is computed differently. Computing a least-square linear regression for the data $(\log _{10}(\sqrt{\text {DOF}_n}), \, \log _{10}(E_n))$ gives a line with slope m, where m is the estimate of rate of convergence, reported in the legend on the right side of each figure.

Table 1 Convergence in the maximum norm (59), for the second- and fourth-order versions of each method, applied to Test Problem 1A (TP–1A), with diffusion coefficient $\lambda = 1$, and time-step $\Delta t = 0.5 h$

Table 2 Convergence in the maximum norm (59), for the second- and fourth-order versions of each method, applied to Test Problem 3A (TP–3A), with diffusion coefficient $\lambda (t) = 1.1 + \sin (\pi t)$, and time-step $\Delta t = 0.5h$

Table 3 Convergence in the maximum norm (60), for the second- and fourth-order versions of each method, applied to Test Problem 2A (TP–2A), with diffusion coefficients $(\lambda _1, \lambda _2) = (10, 1)$, and time-step $\Delta t = 0.5h$. (DPM2-I/DPM4-I refers to the extension of the DPM method, to consider implicit geometry.)

Table 4 Convergence in the maximum norm (60), for the second- and fourth-order versions of each method, applied to Test Problem 2B (TP–2B), with diffusion coefficients $(\lambda _1, \lambda _2) = (10, 1)$, and time-step $\Delta t = 0.5h$. (DPM2-I/DPM4-I refers to the extension of the DPM method, to consider implicit geometry.)

Table 5 Convergence in the maximum norm (60), for the second- and fourth-order versions of each method, applied to Test Problem 2C (TP–2C), with diffusion coefficients $(\lambda _1, \lambda _2) = (1000, 1)$, and time-step $\Delta t = 0.5h$. (DPM2-I/DPM4-I refers to the extension of the DPM method, to consider implicit geometry.)

Overall, we see in Tables 1, 2, 3, 4 and 5 that the error for second-order methods (denoted, for brevity, as CUT2, DPM2, SBP2) on the finest mesh is similar, or sometimes larger, than the error for fourth-order methods (denoted CUT4, DPM4, SBP4) on the coarsest mesh—this illustrates the effectiveness of higher-order methods, when high accuracy is important. Additionally, comparing the three methods together, the size of the errors for the single-domain problems (TP–1A, TP–3A) are similar, up to a constant factor; while for the composite-domain problems (TP–2A, TP–2B, TP–2C) we do see differences of one or two orders of magnitude, with the DPM having the smallest errors.

Fig. 6

Log–log plot of absolute error (59) versus $\sqrt{\text {DOF}}$, and estimated rate of convergence, for the second- and fourth-order versions of each method, applied to Test Problem 1A (TP–1A). See Table 1 for more details

In Table 1 and Fig. 6, we observe that the measured rates of convergence for the numerical approximation of Test Problem 1A (TP–1A) are all $\approx 2$ (for the second-order versions) or $\approx 4$ (for the fourth-order versions), except for DPM4, which for this test problem is superconvergent, with fifth-order convergence. Such higher-than-expected convergence might occur due to several reasons—for example, (i) if the geometry is smooth; (ii) if the magnitude of the derivatives have fast decay (effectively reducing the local truncation error by a factor of h); or (iii) if there is cancellation of error due to symmetries in the geometry, or in the analytical solution.

Fig. 7

Log–log plot of absolute error (59) versus $\sqrt{\text {DOF}}$, and estimated rate of convergence, for the second- and fourth-order versions of each method, applied to Test Problem 3A (TP–3A). See Table 2 for more details

Table 2 and Fig. 7 show the numerical results for (TP–3A). This test problem has the same manufactured solution as (TP–1A), but with a time-varying diffusion coefficient. Despite this added complexity, the numerical results are the same order of accuracy, and in many cases the errors are the same up to seven digits, when compared with the results for (TP–1A). This similarity in the numerical results demonstrates that the three methods can robustly handle time-varying diffusion coefficients.

Fig. 8

Plot of error at the final time $T = 1.0$, for the fourth-order versions of (a) cut-FEM, (b) DPM, and (c) SBP–SAT–FD, respectively, applied to Test Problem 3A (TP–3A)

Fig. 9

Plot of error at the final time $T = 1.0$, for the fourth-order versions of (a) cut-FEM, (b) DPM, and (c) SBP–SAT–FD, respectively, applied to Test Problem 2A (TP–2A)

The plots of spatial error at the final time $T = 1.0$, shown in Fig. 8, are representative of other tests (not included in this text) on a single circular domain. The error in the cut-FEM solution presents largely at the boundary; the error in the DPM solution typically has smooth error, even for grid points very near $\varGamma$; while the error in the SBP–SAT–FD solution is not smooth at interfaces introduced by the domain partitioning.

The plots of spatial error at the final time $T =1.0$ for (TP–2A) are shown in Fig. 9. These plots are fairly representative of the other composite domain tests reported herein, and also of others test problems not included in this work. As in Fig. 8, the cut-FEM has its largest error at degrees of freedom on cut (intersected) elements; the DPM has piecewise smooth error, including even grid points at the boundary/interface; and the SBP–SAT–FD has its largest error at the interfaces between computational subdomains, with particularly pronounced error at the corners of $\varOmega$, where the grid is most stretched.

Fig. 10

Log–log plot of absolute error (60) versus $\sqrt{\text {DOF}}$, and estimated rate of convergence, for the second- and fourth-order versions of each method (cut-FEM, DPM with explicit geometry, SBP–SAT–FD), applied to Test Problem 2A (TP–2A). See Table 3 for more details

Fig. 11

Log–log plot of absolute error (60) versus $\sqrt{\text {DOF}}$, and estimated rate of convergence, for the second- and fourth-order versions of each method (cut-FEM, DPM with explicit geometry, SBP–SAT–FD), applied to Test Problem 2B (TP–2B). See Table 4 for more details

Regarding the max-norm error in presented in Table 3 and Fig. 10, we see that the DPM has smaller max-norm by more than an order of magnitude. We also observe that the convergence rate of the fourth-order SBP–SAT–FD is only three. This suboptimal convergence is inline with the error plot in Fig. 9c, which shows that the error at the corners of the domain is significantly larger than elsewhere. In addition, the error is only non-smooth along the interfaces on the two diagonal lines of the domain. We have also measured the $L_2$ error at the final time $T=1.0$ (not reported in this work), and fourth-order convergence is obtained.

In Table 4 and Fig. 11, we see the numerical results for (TP–2B). The analytical solution is similar to (TP–2A), though much more oscillatory—this additional challenge is manifested by an increase in error by several orders of magnitude.

Fig. 12

Log–log plot of absolute error (60) versus $\sqrt{\text {DOF}}$, and estimated rate of convergence, for the second- and fourth-order versions of each method (cut-FEM, DPM with explicit geometry, SBP–SAT–FD), applied to Test Problem 2C (TP–2C). See Table 5 for more details

In Table 5 and Fig. 12, we see the numerical results for (TP–2C), which shows that our numerical methods are robust to large jumps in diffusion coefficients, the analytical solution, and/or the flux of the true solution. Also, observe that the errors from DPM2/DPM4 (explicit geometry) and DPM2-I/DPM4-I (implicit geometry) in Tables 3, 4 and 5 are almost identical, which demonstrates the robustness and flexibility of the DPM.

6 Discussion

There are many possible methods (Sect. 1) for the numerical approximation of PDE posed on irregular domains, or on composite domains with interfaces. In this work, we consider three such methods, designed for the high-order accurate numerical approximation of parabolic PDEs (1, 2 or 3–9). Each implementation was written, tested, and optimized by the authors most experienced with the method—the cut-Finite Element Method (cut-FEM) by G. Ludvigsson, S. Sticko, G. Kreiss; the Difference Potentials Method (DPM) by K. R. Steffen, Q. Xia, Y. Epshteyn; and the Finite Difference Method satisfying Summation-By-Parts, with a Simultaneous Approximation Term (SBP–SAT–FD) by S. Wang, G. Kreiss. Although we consider only one type of boundary/interface (a circle), we hope that the benchmark problems considered will be a valuable resource, and the numerical results a valuable comparison, for researchers interested in numerical methods for such problems.

The primary differences between the cut-FEM and the standard finite element method are the stabilization terms for near-boundary degrees of freedom, and the quadrature over cut (intersected) elements. Tuning the free parameters in the stabilization terms could mitigate the errors observed in Figs. 8, 9. (We have done some preliminary experiments suggesting that the errors decrease when tuning these parameters, but further investigations are required in order to guarantee robustness.) Given a level-set description of the geometry, there are robust algorithms for constructing the quadrature over cut elements. Together, these differences allow for an immersed (non-conforming) grid to be used. The theoretical base for cut-FEM is well established.

The DPM is based on the equivalence between the discrete system of Eq. (33) and the Boundary Equations with Projection (Theorem 1). The formulation outlined in Sect. 3.2 allows for an immersed (non-conforming) grid; fast $\mathcal {O}(N \log N)$ algorithms, even for problems with general, smooth geometry; and reduces the size of the system to be solved at each time-step. The convergence theory is well-established for general, linear, elliptic boundary value problems, and we conjecture in Sect. 3.2 that this extends to the current setting. In this work, we have extended DPM to work with implicitly-defined geometries for the first time. This is a first step for solving problems where the interface moves with time.

In the finite difference framework (the SBP–SAT–FD method, in this work), the SBP property makes it possible to prove stability and convergence for high-order methods by an energy method. Combined with the SAT method to impose boundary and interface conditions, the SBP–SAT–FD method can be efficient to solve time-dependent PDE. Geometrical features are resolved by curvilinear mapping, which requires an explicit parameterization of boundaries and interfaces. High quality grid generation is important—our experiments, though not reported in this work, have shown that the error in the solution is sensitive to both the orthogonality of the grid and the grid stretching.

Similarities between the cut-FEM and the DPM (beyond the use of an immersed grid) include the thin layer of cut cells along the boundaries/interfaces (cut-FEM) and the discrete grid boundary $\gamma$ (DPM); and the use of higher-order normal derivatives in the stabilization term (cut-FEM) and extension operator (in the Boundary Equations with Projection; DPM). A similarity between the cut-FEM and SBP–SAT–FD is the weak imposition of boundary conditions, via Nitsche’s method (cut-FEM) or the SAT method (SBP–SAT–FD). In this work, the DPM and the SBP–SAT–FD method both use an underlying finite-difference discretization, but the DPM is not restricted to this type of discretization.

Although both the cut-FEM and the DPM use higher-order normal derivatives in their treatment of the boundary/interface, the precise usage differs. For cut-FEM, it is the normal of the element interfaces cut by $\varGamma$, while for DPM, it is the normal of the boundary/interface $\varGamma$. Moreover, in the cut-FEM, stabilization terms (19) involving higher-order normal derivatives at the boundaries of cut-elements are added to the weak form of the PDE, to control the condition number of the mass and stiffness matrices, with a priori estimation of parameters to guarantee positive-definiteness of these matrices; while in the DPM, the Boundary Equations with Projection is combined with the Extension Operator (Definition 4), which incorporates higher-order normal derivatives at the boundary/interface $\varGamma$.

Returning to Sect. 5.3, we see (in Tables 1, 2, 3, 4, 5 and Figs. 6, 7, 8, 9, 10, 11, 12) that the expected rate of convergence for the second- and fourth-order versions of DPM and cut-FEM is achieved, while the DPM has the smallest error constant across all tests. For the SBP–SAT–FD method, expected convergence rates are obtained in some experiments. A noticeable exception is Test Problem 2A, for which the fourth-order SBP–SAT–FD method only has a convergence rate of three.

From the error plot in Fig. 9c, we observe that the large error is localized at the four corners of the domain $\varOmega$, where the curvilinear grid is non-orthogonal and is stretched the most (see Fig. 5b).

As seen in the error plots (Figs. 8, 9), the error for the cut-FEM and the SBP–SAT–FD has “spikes”, while for the DPM the error is smooth. A surprising observation from Fig. 9 is that conforming grids (on which the SBP–SAT–FD method is designed) do not necessarily produce more accurate solutions than immersed grids (on which the cut-FEM and the DPM are designed). Indeed, it is challenging to construct a high-quality curvilinear grid for the considered composite domain problem.

Future directions we hope to consider (in the context of new developments and also further comparisons) include: (i) parabolic problems with moving boundaries/interfaces, (ii) comparison of numerical methods for interface problems involving wave equations [12, 70, 71, 75, 78], (iii) extending our methods to consider PDEs in 3D, (iv) design of fast algorithms, and (v) design of adaptive versions of our methods.

Indeed, for (i), difficulties for the cut-FEM might be the costly construction of quadrature, while for DPM difficulties might be the accurate construction of extension operators. Regarding (iii), this has already been done for the cut-FEM and SBP–SAT–FD; while for the DPM, this is current work, with the main steps extending from 2D to 3D in a straightforward manner.

7 Conclusion

In this work, we propose a set of benchmark problems to test numerical methods for parabolic partial differential equations in irregular or composite domains, in the simplified geometric setting of Sect. 2, with the interface defined either explicitly or implicitly. Next, we compare and contrast three methods for the numerical approximation of such problems: the (i) cut-FEM; (ii) DPM; and (iii) SBP–SAT–FD. Brief introductions of the three numerical methods are given in Sect. 3. It is noteworthy that the DPM has, for the first time, been extended to problems with an implicitly-defined interface.

For the three methods, the numerical results in Sect. 5.3 illustrate the high-order accuracy. Similar errors (different by a constant factor) are observed at grid points away from the boundary/interface, while the observed errors near the boundary/interface vary depending upon the given method. Although we consider only test problems with circular boundary/interface, the ideas underlying the three methods can readily be extended to more general geometries.

In general, all three methods require an accurate and efficient resolution of the explicitly- or implicitly-defined irregular geometry: cut-FEM relies on accurate quadrature rules for cut elements, and a good choice of stabilization parameters; DPM relies on an accurate and efficient representation of Cauchy data using a good choice of basis functions; and SBP–SAT–FD relies on the smooth parametrization to generate a high-quality curvilinear grid.

Appendix (DPM)

Let us now expand some details presented in the brief introduction to the Difference Potentials Method (Sect. 3.2).

1.1 Fully-Discrete Formulation of (3, 4)

The fully-discrete, finite-difference discretization introduced in (33) is

$$\begin{aligned} L^s_{\Delta t,h}u_s^{i+1}=F_s^{i+1},\quad (x_j,y_k)\in M^+_s. \end{aligned}$$

(62)

The general form of the operator is $L^s_{\Delta t,h}=\lambda _s(t^{i+1})\Delta _h-\sigma I$, where $\sigma =\frac{3}{2\Delta t}$ for second-order (BDF2–DPM2), $\sigma =\frac{25}{12\Delta t}$ for fourth-order (BDF4–DPM4), and $\Delta _h$ is either a standard five- or nine-point Laplacian. For the nine-point Laplacian, we have

$$\begin{aligned} \Delta _h u_{j,k}&=\frac{1}{12h^2}(-u_{j-2,k}+16u_{j-1,k}+16u_{j+1,k}-u_{j+2,k}-u_{j,k-2}\nonumber \\&\quad +\,16u_{j,k-1}+16u_{j,k+1}-u_{j,k+2}-60u_{j,k}) \end{aligned}$$

(63)

for points sufficiently far away from the boundary of the auxiliary domain $\varOmega _s^0$. For points that are close to the boundary, we use a modified, fourth-order stencil. For example, at the southwest corner, we take

$$\begin{aligned} \Delta _h u_{1,1}&=\frac{1}{12h^2}(10u_{0,1}-4u_{2,1}+14u_{3,1} -6u_{4,1}+u_{5,1}\nonumber \\&\quad +~10u_{1,0}-4u_{1,2}+14u_{1,3}-6u_{1,4}+u_{1,5}-30u_{1,1}), \end{aligned}$$

(64)

where $u_{0,1}$ and $u_{1,0}$ will be from the boundary condition (7).

Next, the right-hand side of (33) for BDF2–DPM2 is given by

$$\begin{aligned} F^{i+1}_s=-f^{i+1}_s-\frac{\sigma }{3}(4u_s^{i}-u^{i-1}_s), \end{aligned}$$

(65)

and for BDF4–DPM4 by

$$\begin{aligned} F^{i+1}_s=-f^{i+1}_s-\frac{\sigma }{25}(48u_s^{i}-36u^{i-1}_s+16u^{i-2}_s-3u^{i-3}_s). \end{aligned}$$

(66)

Lastly, the initialization at $t=0$ is done using the exact solutions for the terms $u^{0}_s$, $u^{-1}_s$, $u^{-2}_s$, and $u^{-3}_s$. Another possibility would be the use of lower-order BDF methods, then advancing to higher-order methods once enough stages are established. No significant differences were observed between the two approaches.

1.2 Equation-Based Extension

Let us now expand the discussion surrounding (39–41) leading up to Definition 4 of the Extension Operator (42).

An important step in this discussion is to recast the original PDE (3, 4) into a curvilinear form, for points (x, y) in the vicinity of $\varGamma$. Following the notation [58], let us first introduce the coordinate system $(d, \vartheta )$ for points in the vicinity of $\varGamma$. Recall from Definition 42 that d is the distance in the normal direction from a given point to its orthogonal projection on $\varGamma$, while $\vartheta$ is the arclength along $\varGamma$ from some reference point to the orthogonal projection. In this coordinate system, the PDE (3, 4) becomes

$$\begin{aligned} \frac{\partial u_s}{\partial t} - \lambda _s\left( \frac{1}{H_\vartheta }\left[ \frac{\partial }{\partial n}\left( H_\vartheta \frac{\partial u_s}{\partial n}\right) +\frac{\partial }{\partial \vartheta }\left( \frac{1}{H_\vartheta }\frac{\partial u_s}{\partial \vartheta }\right) \right] \right) =f_s, \end{aligned}$$

(67)

where where $H_\vartheta =1 - d \kappa$ is the Lamé coefficient, and $\kappa$ is the signed curvature along the interface $\varGamma$.

From (67), a straightforward calculation gives the second-order normal derivative $\frac{\partial ^2 u_s}{\partial n^2}$ (used in the calculation of (41)), which is

$$\begin{aligned} \frac{\partial ^2 u_s}{\partial n^2} = \frac{1}{\lambda _s}\left( \frac{\partial u_s}{\partial t}-f_s\right) -\frac{\partial ^2u_s}{\partial \vartheta ^2}+\kappa \frac{\partial u_s}{\partial n}. \end{aligned}$$

(68)

For the fourth-order numerical method, which uses an Extension Operator with $p=4$, we also need the third- and fourth-order normal derivatives, which we state now. Differentiating (68) with respect to n, we see that

$$\begin{aligned} \frac{\partial ^3 u_s}{\partial n^3}=\frac{1}{\lambda _s}\left( \frac{\partial ^2u_s}{\partial t\partial n}-\frac{\partial f_s}{\partial n}\right) -\frac{\partial ^3u_s}{\partial n\partial \vartheta ^2}+\kappa \frac{\partial ^2 u_s}{\partial n^2} \end{aligned}$$

(69)

and

$$\begin{aligned} \frac{\partial ^4u_s}{\partial n^4}&=\frac{1}{\lambda _s^2}\left( \frac{\partial ^2u_s}{\partial t^2}-\frac{\partial f_s}{\partial t}\right) \nonumber \\&\quad \cdots +\frac{1}{\lambda _s}\left( -2\frac{\partial ^3u_s}{\partial t\partial \vartheta ^2}+\kappa \frac{\partial ^2u_s}{\partial n\partial t}-\frac{\partial ^2f_s}{\partial n^2}+\frac{\partial ^2f_s}{\partial \vartheta ^2}+\frac{\partial ^4u_s}{\partial \vartheta ^4}-\kappa \frac{\partial ^3u_s}{\partial n\partial \vartheta ^2}\right) +\kappa \frac{\partial ^3 u_s}{\partial n^3}. \end{aligned}$$

(70)

Next, let us follow-up on comments made in the text following (41). There, it was pointed out that the unknown Dirichlet and Neumann data $\left( u_s,\frac{\partial u_s}{\partial n}\right)$ are the only data required for the Extension Operator (42) with $p = 2$. Moreover, it was pointed out that this is also true for the Extension Operator when $p = 4$. The reasoning is as follows.

• The time derivatives $\frac{\partial u_s}{\partial t}$, $\frac{\partial ^2u_s}{\partial t^2}$, $\frac{\partial f_s}{\partial t}$, $\frac{\partial ^2 u_s}{\partial n\partial t}$, and $\frac{\partial ^3u_s}{\partial t\partial \vartheta ^2}$ can be approximated by the backward difference formula. In BDF2–DPM2,

  $$\begin{aligned} \frac{\partial u_s^{i+1}}{\partial t}&\approx \frac{3u_s^{i+1}-4u_s^i+u_s^{i-1}}{2\Delta t} \quad \text { and }\end{aligned}$$

  (71)
  $$\begin{aligned} \frac{\partial ^2 u_s^{i+1}}{\partial t^2}&\approx \frac{2u_s^{i+1}-5u_s^i+4u_s^{i-1}-u_s^{i-2}}{\Delta t^2}, \end{aligned}$$

  (72)

  while in BDF4–DPM4,

  $$\begin{aligned} \frac{\partial u_s^{i+1}}{\partial t}&\approx \frac{25u_s^{i+1}-48u_s^i+36u_s^{i-1}-16u_s^{i-2}+3u_s^{i-3}}{12\Delta t} \quad \text { and }\end{aligned}$$

  (73)
  $$\begin{aligned} \frac{\partial ^2 u_s^{i+1}}{\partial t^2}&\approx \frac{1}{\Delta t^2}\left( \frac{15}{4}u_s^{i+1}-\frac{77}{6}u_s^i+\frac{107}{6}u_s^{i-1}-13u_s^{i-2}+\frac{61}{12}u_s^{i-3}-\frac{5}{6}u_s^{i-4}\right) .\nonumber \\ \end{aligned}$$

  (74)

• The derivatives in terms of arclength $\vartheta$ can be computed from $u_s$ or $\frac{\partial u_s}{\partial n}$. For example, denoting $u_s=\sum _{\nu =1}^{\mathcal {N}^0} c_{1,\nu }^{s,i+1}\phi _{\nu }(\vartheta )$ (using notation following from (46)), then it comes handy that

  $$\begin{aligned} \frac{\partial u_s}{\partial \vartheta }=\sum _{\nu =1}^{\mathcal {N}^0} c_{1,\nu }^{s, i+1}\phi _\nu '(\vartheta ),\quad \frac{\partial ^2 u_s}{\partial \vartheta ^2}=\sum _{\nu =1}^{\mathcal {N}^0} c_{1,\nu }^{s, i+1}\phi _\nu ''(\vartheta ),\quad \frac{\partial ^4 u_s}{\partial \vartheta ^4}=\sum _{\nu =1}^{\mathcal {N}^0} c_{1,\nu }^{s, i+1}\phi _\nu ^{(4)}(\vartheta ). \end{aligned}$$

  (75)

1.3 The System of Equations at Each Time Step

With the Cauchy data $\mathfrak {u}_{s, \varGamma }^{i + 1}$ and Extension Operator ${\text {Ex}}_s \mathfrak {u}_{s, \varGamma }^{i + 1}$ from $\varGamma$ to $\gamma _s$ introduced in Definition 4, and the spectral representation introduced in (46), we now give a sketch of the linear system for the coefficients $(c_{1, \nu }^{s, i+1})_{\nu =1}^{\mathcal {N}^0}$ and $(c_{2, \nu }^{s, i+1})_{\nu =1}^{\mathcal {N}^1}$, and moreover the approximation of the solution $u_s(x, y, t^{i+1})$ at $(x_j, y_k) \in N_s^+$.

Indeed, substituting ${\text {Ex}}_s \tilde{\mathfrak {u}}_{s, \varGamma }^{i + 1}$ (42, 46) into the BEP (38), the resulting linear systems are

$$\begin{aligned} \sum _{k = 0}^1 \sum _{\nu =1}^{\mathcal {N}^k} \left( c_{k+1, \nu }^{s, i + 1} {\text {Ex}}_s \Phi ^k_{\nu } - c_{k+1, \nu }^{s, i + 1} P_{\gamma _s}^{i + 1} {\text {Ex}}_s \Phi ^k_{\nu } \right) = {\text {Tr}}_{\gamma _s} [G_{\Delta t, h}^{i+1} F_s^{i + 1}]. \end{aligned}$$

(76)

This can be further elucidated by introducing the vector of unknowns

$$\begin{aligned} \mathbf {c}_s^{i + 1} = \big [ \underbrace{c_{1, 1}^{s, i + 1} \quad c_{1, 2}^{s, i + 1} \quad \cdots \quad c_{1, \mathcal {N}^0}^{s, i + 1}}_{\mathbf {c}_{s, 1}^{i + 1}} \quad \underbrace{c_{2, 1}^{s, i + 1} \quad c_{2, 2}^{s, i + 1} \quad \cdots \quad c_{2, \mathcal {N}^1}^{s, i + 1}}_{\mathbf {c}_{s, 2}^{i + 1}} \big ]^\top \end{aligned}$$

(77)

(so that $\mathbf {c}_s^{i + 1} = [ \mathbf {c}_{s,1}^{i + 1}, \mathbf {c}_{s,2}^{i + 1} ]^\top$), and the matrix

$$\begin{aligned} A_s = \Big [&\underbrace{(I - P_{\gamma _s}^{i + 1}) {\text {Ex}}_s \Phi ^0_1, \quad (I - P_{\gamma _s}^{i + 1}) {\text {Ex}}_s \Phi ^0_2, \quad \cdots \quad (I - P_{\gamma _s}^{i + 1}) {\text {Ex}}_s \Phi ^0_{\mathcal {N}^0}}_{A_{s, 1}}, \nonumber \\&\underbrace{(I - P_{\gamma _s}^{i + 1}) {\text {Ex}}_s \Phi ^1_1, \quad (I - P_{\gamma _s}^{i + 1}) {\text {Ex}}_s \Phi ^1_2, \quad \cdots \quad (I - P_{\gamma _s}^{i + 1}) {\text {Ex}}_s \Phi ^1_{\mathcal {N}^1} \Big ]}_{A_{s, 2}}. \end{aligned}$$

(78)

Then, the full system of Eq. (76) is

$$\begin{aligned} A \begin{bmatrix} \mathbf {c}_1^{i + 1} \\ \mathbf {c}_2^{i + 1} \end{bmatrix} = \begin{bmatrix} {\text {Tr}}_{\gamma _1} [G_{\Delta t, h}^{i + 1} F^{i + 1}_1] \\ {\text {Tr}}_{\gamma _2} [G_{\Delta t, h}^{i + 1} F^{i + 1}_2] \end{bmatrix}, \quad \text {with } A = \begin{bmatrix} A_1&0 \\ 0&A_2 \end{bmatrix}. \end{aligned}$$

(79)

However, note that $\mathbf {c}_1^{i + 1}$ and $\mathbf {c}_2^{i + 1}$ are related by the interface conditions (8, 9), so that the number of unknowns in (79) is equal the dimension of either $\mathbf {c}_1^{i + 1}$ or $\mathbf {c}_2^{i + 1}$, depending on which one is considered the independent unknown. Therefore, the dimension of A is $(|\gamma _1| + |\gamma _2|) \times (\mathcal {N}^0 + \mathcal {N}^1)$, where $\mathcal {N}^0 + \mathcal {N}^1$ is the dimension of $\mathbf {c}_1^{i + 1}$ or $\mathbf {c}_2^{i + 1}$ (whichever is the independent unknown).

Remark 8

The independent unknown ($\mathbf {c}_1^{i + 1}$ or $\mathbf {c}_2^{i + 1}$) is chosen so that the finite-dimensional, spectral representation (46) of the Cauchy data $\mathfrak {u}_{s,\varGamma }^{i+1}$ accurately resolves the Cauchy data with a small number of basis functions, in the consideration of both accuracy and computational efficiency. For (TP–2A) and (TP–2B), we choose $\mathbf {c}_2^{i + 1}$ as the independent unknown, while for (TP–2C) we choose $\mathbf {c}_1^{i + 1}$. With these choices for the independent unknown, we have $\mathcal {N}^0=\mathcal {N}^1=1$ for the three considered test problems.

Since each column involves the Difference Potentials operator $P_{\gamma _s}^{i+1}$ applied to a vector ${\text {Ex}}_s \Phi _\nu ^k$, each column is therefore constructed via one solution of the Auxiliary Problem (Definition 1). However, the Auxiliary Problems are posed on the computationally simple Auxiliary Domains, and can be computed using a fast FFT- or multigrid-based algorithm, which can significantly reduce the computational cost. Moreover, if $\lambda _s(t) \equiv \lambda _s$ is constant, then A can be computed and inverted once (as a pre-processing step), thus significantly reducing computational cost for long-time simulations.