High-order numerical methods for 2D parabolic problems in single and composite domains

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.


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, finiteelement approximation, etc., may fail to produce an accurate approximation near the interface, leading one to consider and develop new techniques.
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 (Section 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 (Section 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 Section 2 we will introduce two circular geometries, which are defined either explicitly, or implicitly via a level set function.
In Section 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.
The paper is outlined as follows. In Section 2, we give brief overview of the continuous formulation of the parabolic problems in a single domain or a composite domain. In Section 3, we give introductions to the basics of the three proposed methods: cut-FEM, DPM, and SBP-SAT-FD. In Section 4, we formulate the numerical test problems. In Section 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 Section 6, we give a comparative discussion of the three methods and the numerical results, together with a discussion on future research directions. Lastly, in Section 7, we give our concluding remarks.  Figure 1 The (a) single domain Ω and (b) composite domain Ω = Ω 1 ∪ Ω 2 . In (b), ∂Ω 1 has two connected components: the boundary ∂Ω and interface Γ = ∂Ω 2 .

Statement of problem
In this section, we describe two diffusion problems, which will be the setting for our proposed benchmark (test) problems in Section 4. (Recall from Section 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.
Remark 1. We consider the circular geometries depicted in Figure 1 as the geometrical setting for our proposed benchmark problems in this work. In applications (Section 1), other geometries will likely be considered, some much more complicated than Figure 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)(4)(5)(6)(7)(8)(9). As such, we think that the geometries in Figure 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 (λ 1 = λ 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.

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 Ω s be covered by a structured triangulation, T s , so that each element T ∈ T s has some part inside of Ω s ; see Figures 2a and 2b. 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 T covers Ω.) Typically T 1 and T 2 would be created from a larger mesh by removing some of the cells. Further, let T Γ = {T ∈ T : T ∩ Γ = ∅} be the set of intersected elements; see Figure 2c. In the following, we shall use Γ 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 T Γ clearer.  Figure 2 The (a) subdomain Ω 1 immersed in a mesh T 1 , (b) subdomain Ω 2 immersed in a mesh T 2 and domain Ω immersed in T , and (c) intersected elements T Γ .
To construct the finite element spaces we use Lagrange elements with Gauss-Lobatto nodes of order p (Q p -elements). Let V s h denote a continuous finite element space on Ω s , consisting of Q p -elements on the mesh T s : For the single domain problem (1, 2) we solve for the solution u ∈ V h ; while for the composite domain problem (3-9), we solve for the pair For the latter problem, this means that the degrees of freedom are doubled over elements belonging to T Γ .
We begin by stating the weak formulation for the single domain problem (1,2). Let (·, ·) X and ·, · Y be the L 2 scalar products taken over the two-and one-dimensional domains X ⊂ R 2 and Y ⊂ 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 ∈ V h , and integrating by parts, we obtain: Note that (2) is consistent with the following terms: where γ 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 where For T Γ (the elements intersected by Γ), note that one must integrate only over the part of the element that lies inside Ω. A problem with this is that one cannot control how the intersections (cuts) between Ω and T are made. Depending on how Ω is located with respect to the triangulation, some elements can have an arbitrarily small intersection with the domain -see, for example, Figure 3a. If Ω is moved with respect to 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 Ω is located with respect to the triangulation T [14,53]. Adding stabilization to (14) results in the following weak form: Find u ∈ V h such that where γ M and γ A are scalar constants.
In order to state the definition of stabilization (19), denote by F s the set of faces, as seen in Figures 3b and 3c. That is, F s is the set of all faces of the elements in T Γ , excluding the boundary faces of T s : Then, the stabilization term is defined as: where [u] = u| F+ − u| F− is the jump over a face, F ; n refers to a normal of F ; and ∂ k n 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)(4)(5)(6)(7)(8)(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 ∈ V s h , 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 F 1 and F 2 , respectively.
where the bilinear forms M and A correspond to the stabilized mass and stiffness matrices: L Ω corresponds to the forcing function: a Γ and L Γ consistently enforce the interface conditions (8,9): and the terms a ∂Ω and L ∂Ω enforce the boundary condition (7) along the outer boundary, ∂Ω: In (24)(25)(26)(27), n denotes the outward pointing normal at either Γ or ∂Ω (depending on the domain of integration); κ 1 + κ 2 = 1, so that {v} = κ 1 v 1 + κ 2 v 2 is a convex combination; and γ Γ , κ 1 , κ 2 are chosen as in [13]: The remaining parameters (appearing in Equations 21,22,[26][27][28] are given by: The scaling of γ 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 γ M was chosen based on numerical experiments on the condition number of the mass matrix. This also agrees with the choice of γ A and γ D in [14], where γ A was investigated numerically. In order to use cut-FEM, one needs a way to perform integration over the intersected elements T Γ . For example, with the interface problem, on each element K ∈ T Γ , we need a quadrature rule for the K ∩ Ω 1 , K ∩ Ω 2 and K ∩ Γ. For the numerical tests in this work (Section 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.

DPM
We continue in this section with a brief introduction to the Difference Potentials Method (DPM), which was originally proposed by V. S. Ryaben'kii (see [64,67,65], 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, finitedifference 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 [66] and [3,4,5,6,26,28]. We also introduce new developments here for handling implicitly-defined geometries. (The reader can consult [67] for the general theory of the Difference Potentials Method.) Broadly, the main idea of the DPM is to reduce uniquely solvable and wellposed boundary value problems in a domain Ω to pseudo-differential Boundary Equations with Projections (BEP) on the boundary of Ω. 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 Γ. (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 Ω 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 O(N log N ) solvers. Moreover, in general the DPM can be applied to problems with general boundary or interfaces 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)(4)(5)(6)(7)(8)(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 Ω s (s = 1, 2) be embedded in a rectangular auxiliary domain Ω 0 s . Introduce a uniform, Cartesian grid denoted M 0 s on Ω 0 s , with grid-spacing h s . Let M + s = M 0 s ∩ Ω s denote the grid points inside each subdomain Ω s , and M − s = M 0 s \ M + s the grid points outside each subdomain Ω s . Note that the auxiliary domains Ω 0 1 , Ω 0 2 and auxiliary grids M 0 1 , M 0 2 need not agree, and indeed may be selected completely independently, given considerations regarding accuracy, adaptivity, or efficiency.
Define a finite-difference stencil N s,α j,k , with α = 5, 9, to be the stencil of the standard five-point or a wide nine-point Laplacian, i.e., Next, with α fixed, define the point-sets Lastly, we now define the important point-set which we call the discrete grid boundary. In words, γ s is the set of grid points that straddle the continuous interface Γ. (See Figure 4 for an example of these points-sets, given a single elliptical domain Ω.) Note that the point-sets M + s , N + s , and γ s will be used throughout the Difference Potentials Method.
Here, we define the fully-discrete finite-difference discretization of (3,4), and then define the Auxiliary Problem. Indeed, the discretization we consider is 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 righthand side grid function q i+1 s : M 0 s → R, the following difference equations (34,35) are defined as the discrete AP.
Remark 4. For a given right-hand side q i+1 s , the solution of the discrete AP (34, 35) defines a discrete Green's operator G s ∆t,h q i+1 s . The choice of boundary conditions (35) will affect the resulting grid function G s ∆t,h q i+1 s , 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 s ∆t,h F i+1 s the particular solution on N + s of the fully-discrete problem (33), defined by solving the AP (34, 35) with and restricting the solution from N 0 s to N + s . Let us also introduce a linear space V γs of all grid functions denoted v i+1 γs , which are defined on γ s and extended by zero to the other points of N 0 s . These grid functions are referred to as discrete densities on γ s . Definition 2 (The Difference Potential of a density.). The Difference Potential of a given density v i+1 γs is the grid function P i+1 and restricting the solution from N 0 s to the point-set N + s .
Proof. See [67] 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.
Remark 5. A given density v i+1 γs is the trace of some solution of the fully-discrete finite-difference equations (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 i+1 γs . 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 γ s , using the Cauchy data of (3-9) on the continuous interface Γ, 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 Γ.) Below, we will define an Extension Operator which will extend the Cauchy data of (3-9) from Γ to γ 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 Γ. 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 finitedimensional, 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 Γ, we define a coordinate system (d, ϑ), where ϑ is arclength from some reference point, and d is the signed distance in the normal direction from the point to Γ. Now, as a first step towards defining the Extension Operator, we define a new function where n is the unit outward normal vector at Γ. 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 ∂n , etc. As a last step, a straightforward sequence of calculations (see Appendix 8.2) shows that where κ denotes the curvature of Γ. Therefore, with v i+1 s (d, ϑ) defined by (39)(40)(41), the only unknown data at each time step t i+1 are the unknown Dirichlet data u i+1 s and the unknown Neumann data ∂u i+1 s ∂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 u i+1 1 , . (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 Γ to γ s .
For a given point (x j , y k ) ∈ γ s , note that d is the signed distance between (x j , y k ) and its orthogonal projection on Γ, while ϑ is the arclength along Γ 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 u i+1 s,Γ . Indeed, we wish to choose a basis φ ν (ϑ) on Γ (ν = 1, 2, 3, . . .) in order to accurately approximate the two components of the Cauchy data u i+1 s,Γ . To be specific, whichever basis we choose, we require that tends to zero as N 0 , N 1 → ∞, for some sequence of real numbers (c s,i+1 and (c s,i+1 2,ν ) N 1 ν=1 . In other words, we require Now let us discuss a choice of basis. In this work, recall that we consider interfaces Γ that are at least C 2 (Γ) (due to the choice of smooth, circular geometries). Also, as we will see in Section 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| Ω1 = u 1 and u| Ω2 = u 2 are smooth in Ω 1 and Ω 2 , respectively. Moreover, each component of the Cauchy data u i+1 1,Γ and u i+1 2,Γ are smooth, periodic functions of arclength ϑ. (Note that u i+1 1,Γ and u i+1 2,Γ need not agree, and indeed do not -neither µ 1 (x, y, t) nor µ 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 φ ν (ϑ), with and k > 1. Moreover, at every time step t i+1 , we will discretize the Cauchy data u i+1 s,Γ = u i+1 s , ∂u i+1 s ∂n Γ using this basis. Therefore, we let where Φ 0 ν = (φ ν , 0) and Φ 1 ν = (0, φ ν ) are the set of basis functions used to represent the Cauchy data on the interface Γ.
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 Γ) 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).
Definition 5 (Discrete, generalized Green's formula.). At each time step t i+1 , the numerical approximation u i+1 s ≈ u s (x j , y k , t i+1 )| (xj ,y k )∈N + s of (3-9) is given by u i+1 Here, u i+1 γs = Ex sũ In this work, we also propose a novel feature of DPM, extending the method originally developed in [66] 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 Γ, 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 interpolantF (x, y) of F (x, y) on a subset of M 0 near the continuous interface Γ. Next, we parameterize Γ by arclength using numerical quadrature. With this parameterization, we (i) compute the Fourier series expansion from initial conditions for the Cauchy data u i+1 s,Γ on the implicit interface Γ, 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 [62,63,67] and [33]), 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 (Section 5) we see that the computed solution (47) at every time level t i+1 has accuracy O(h 2 + ∆t 2 ) for the secondorder method, and O(h 4 +∆t 4 ) for the fourth-order method, for both the single and composite domain problems, with explicit or implicit geometry. See [3,6,27,66] 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 2 : At each time step t i+1 , compute the Particular Solution u i+1 , (x j , y k ) ∈ N + s , using the AP with the right-hand side (36). • Step 4 : Compute the approximation of the density u i+1 γs , by applying the Extension Operator (42) to the solution of (79).
• Step 5 : Construct the Difference Potentials P N + s γs u i+1 γs of the density u i+1 γs , using the AP with the right-hand side (37).

SBP-SAT-FD
We continue in this section with a brief presentation of SBP-SAT-FD, for solving the parabolic problems presented in Section 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 π/2) with curved sides. Similarly, the composite domain of equation (3-9) is divided into nine subdomains, as shown in Figure 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].
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 in a physical domain will be transformed to where (ξ, η) is the Cartesian coordinate in the unit square, and J(ξ, η), α(ξ, η), β(ξ, η), γ(ξ, η) 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 ≈ ∂/∂x to approximate a first derivative, and D (b) 2 ≈ ∂/∂x(b(x)∂/∂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 ∈ [0, 1]. We discretize [0, 1] uniformly by N grid points, and denote the restriction of u(x), v(x) onto the grid by u, v, respectively. Integration by parts states: The SBP operator D 1 mimics integration by parts: where H is symmetric positive definite -thus defining an inner product -and B = diag(−1, 0, · · · , 0, 1).
In fact, H is also a quadrature [20]. It is easy to verify that (51) is equivalent to 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 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 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 with Dirichlet boundary conditions, Neumann boundary conditions, and interface conditions. Comparing equation (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 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) 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 2 Q is symmetric and negative semi-definite. Therefore, we can write (55) as d dt (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 is also stable. To see this, we multiply (56) by v T H 2 from the left, and obtain d dt where v T H 2 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)(4)(5)(6)(7)(8)(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 equation (49), but are different from the ones used in this paper.

Test Problems
In this section, we first list the test problems that we will consider (in Section 4.1), and then briefly motivate and discuss these choices (in Section 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 λ(t) or (λ 1 , λ 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 (µ 1 (x, y, t), µ 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. (TP-2C)

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.

Time discretization
The spatial discretization for each method is discussed in Section 3. For the time discretization, the backward differentiation formulas of second-and fourthorder (BDF2 and BDF4) are used for the second-and fourth-order methods, respectively. In each case, the time-step is given by 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.

Measure for comparison
Let u i j,k denote the computed numerical approximation of u(x, y, t) at the gridpoint (x j , y k ) ∈ Ω and time t i = i∆t ∈ (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: and for the composite domain problem (3-9) as:

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-5, the estimate of rate of convergence is computed as follows. Let (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 which is the estimated rate of convergence, denoted in Tables 1-5 by "Rate". In Figures 6, 7, 10-12, the estimate of rate of convergence is computed differently. Computing a least-square linear regression for the data (log 10 ( √ 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. Overall, we see in Tables 1-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.  Figure 6 Log-log plot of absolute error (59) versus √ 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 Figure 6, we observe that the measured rates of convergence for the numerical approximation of Test Problem 1A (TP-1A) are all ≈ 2 (for the second-order versions) or ≈ 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.  (59), for the second-and fourthorder versions of each method, applied to Test Problem 3A (TP-3A), with diffusion coefficient λ(t) = 1.1 + sin(πt), and time-step ∆t = 0.5h.  Table 2 and Figure 7 show the numerical results for (TP-3A). This test problem has the same manufactured solution as (TP-1A), but with a timevarying 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.
The plots of spatial error at the final time T = 1.0, shown in Figure 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 Γ; 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 Figure 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 Figure 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  Table 2 for more details.
the interfaces between computational subdomains, with particularly pronounced error at the corners of Ω, where the grid is most stretched.
Regarding the max-norm error in presented in Table 3 and Figure 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 Figure 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 Figure 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.
In Table 5 and Figure 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-5 are almost identical, which demonstrates the robustness and flexibility of the DPM.      (60), for the second-and fourthorder versions of each method, applied to Test Problem 2A (TP-2A), with diffusion coefficients (λ 1 , λ 2 ) = (10, 1), and time-step ∆t = 0.5h. (DPM2-I/DPM4-I refers to the extension of the DPM method, to consider implicit geometry. )

Discussion
There are many possible methods (Section 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- 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 Figures 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  Table 3 for more details.
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 equations (33) and the Boundary Equations with Projection (Thm. 1). The formulation outlined in Section 3.2 allows for an immersed (non-conforming) grid; fast 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 Section 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 γ (DPM); and the use of higher-order Table 4 Convergence in the maximum norm (60), for the second-and fourthorder versions of each method, applied to Test Problem 2B (TP-2B), with diffusion coefficients (λ 1 , λ 2 ) = (10, 1), and time-step ∆t = 0.5h. (DPM2-I/DPM4-I refers to the extension of the DPM method, to consider implicit geometry. ) 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 Γ, while for DPM, it is the normal of the boundary/interface Γ. 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 Γ.
Returning to Section 5.3, we see (in Tables 1-5 Figure 11 Log-log plot of absolute error (60) versus √ 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.
fourth-order SBP-SAT-FD method only has a convergence rate of three. From the error plot in Figure 9c, we observe that the large error is localized at the four corners of the domain Ω, where the curvilinear grid is non-orthogonal and is stretched the most (see Figure 5b).
As seen in the error plots (Figures 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 Figure 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.

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 Section 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 Section 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 Section 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.  Table 5 for more details.

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

Equation-based extension
Let us now expand the discussion surrounding (39)(40)(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 Γ. Following the notation [58], let us first introduce the coordinate system (d, ϑ) for points in the vicinity of Γ. Recall from Definition 42 that d is the distance in the normal direction from a given point to its orthogonal projection on Γ, while ϑ is the arclength along Γ from some reference point to the orthogonal projection. In this coordinate system, the PDE (3, 4) becomes where where H ϑ = 1 − dκ is the Lamé coefficient, and κ is the signed curvature along the interface Γ. From (67), a straightforward calculation gives the second-order normal derivative ∂ 2 us ∂n 2 (used in the calculation of (41)), which is ∂ 2 u s ∂n 2 = 1 λ s ∂u s ∂t − f s − ∂ 2 u s ∂ϑ 2 + κ ∂u s ∂n .
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 and 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 u s , ∂us ∂n 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 system of equations at each time step.
With the Cauchy data u i+1 s,Γ and Extension Operator Ex s u i+1 s,Γ from Γ to γ 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 s,i+1 1,ν ) N 0 ν=1 and (c s,i+1 2,ν ) N 1 ν=1 , and moreover the approximation of the solution u s (x, y, t i+1 ) at (x j , y k ) ∈ N + s . Indeed, substituting Ex sũ (so that c i+1 s = [c i+1 s,1 , c i+1 s,2 ] ), and the matrix A s = (I − P i+1 γs ) Ex s Φ 0 1 , (I − P i+1 γs ) Ex s Φ 0 2 , · · · (I − P i+1 γs ) Ex s Φ 0 N 0 As,1 , (I − P i+1 γs ) Ex s Φ 1 1 , (I − P i+1 γs ) Ex s Φ 1 2 , · · · (I − P i+1 γs ) Ex s Φ 1 Then, the full system of equations (76) is However, note that c i+1 1 and c i+1 2 are related by the interface conditions (8,9), so that the number of unknowns in (79) is equal the dimension of either c i+1 1 or c i+1 2 , depending on which one is considered the independent unknown. Therefore, the dimension of A is (|γ 1 | + |γ 2 |) × (N 0 + N 1 ), where N 0 + N 1 is the dimension of c i+1 1 or c i+1 2 (whichever is the independent unknown). 2 ) is chosen so that the finitedimensional, spectral representation (46) of the Cauchy data u i+1 s,Γ 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 c i+1 2 as the independent unknown, while for (TP-2C) we choose c i+1 1 . With these choices for the independent unknown, we have N 0 = N 1 = 1 for the three considered test problems.
Since each column involves the Difference Potentials operator P i+1 γs applied to a vector Ex s Φ 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 λ s (t) ≡ λ 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.