High-order numerical methods for 2D parabolic problems in single and composite domains
Gustav Ludvigsson Kyle R. Steffen Simon Sticko Siyang Wang Qing Xia Yekaterina Epshteyn Gunilla Kreiss
math.NA 65M22 (Primary) 65M06, 65M12, 65M55, 65M60, 65M70, 35K20 (Secondary)







Wigner tomography of multispin quantum states
    Steffen J. Glaser
    January 15, 2018
=============================================




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. 

Keywords: parabolic problems; interface models; level set; complex geometry; discontinuous solutions; SBP–SAT finite difference; difference potentials; spectral approach; finite element method; cut elements; immersed boundary; stabilization; higher order accuracy and convergence;

AMS Subject Classification: 65M06, 65M12, 65M22, 65M55, 65M60, 65M70, 35K20








§ 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 <cit.>, 
difference potentials method <cit.>, 
immersed boundary method <cit.>, 
immersed interface method <cit.>, 
ghost fluid method <cit.>,  
the matched interface and boundary method <cit.>, 
Cartesian grid embedded boundary method <cit.>, 
multigrid method for elliptic problems with discontinuous coefficients on an arbitrary interface <cit.>, 
virtual node method <cit.>, 
Voronoi interface method <cit.>, 
summation–by–partsthe finite difference method <cit.>  and finite volume method <cit.> based on mapped grids,Comment 11 of reviewer 1 
or cut finite element method <cit.>. 
In spite ofIndeed, 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 efficientcomputationally-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.* = Comment 2 of Reviewer 1

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 <ref>) 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.).Comment 3 Reviewer 1 
In particular, the formulatedThe formulated testComment 3 Reviewer 1 problems (Section <ref>) 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 <ref> we will introduce two circular geometries, which are defined either explicitly, or implicitly via a level set function. 

In Section <ref>, 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 <cit.>.   For comparison with cut–FEM and DPM, however, in this paper we use a conforming approach based on the finite difference method – the SBP–SAT–FD method, which resolves geometrical features by curvilinear mapping.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.Comment 10 of reviewer 1

For recent work on SBP–SAT–FD for wave equations in composite domains, see <cit.>, and the two review papers <cit.>; for recent work in DPM for elliptic/parabolic problems in composite domains with interface defined explicitly, see <cit.>; and for recent work in cut–FEM see <cit.>. 

The paper is outlined as follows. In Section <ref>, we give brief overview of the continuous formulation of the parabolic problems in a single domain or a composite domain.  In Section <ref>, we give introductions to the basics of the three proposed methods: cut–FEM, DPM, and SBP–SAT–FD. In Section <ref>, we formulate the numerical test problems. In Section <ref>, 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 <ref>, we give a comparative discussion of the three methods and the numerical results, together with a discussion on future research directions. Lastly, in Section <ref>, we give our concluding remarks. 





















§ STATEMENT OF PROBLEM
 





In this section, we describe two diffusion problems, which will be the setting for our proposed benchmark (test) problems in Section <ref>.  (Recall from Section <ref> that these models arise, for example, in the study of mass or heat diffusion.)Comment 3 Reviewer 1  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. 



 §.§ The single domain problem


  
  First, we consider the linear parabolic PDE on a single domain Ω (e.g., Figure <ref>), with variable diffusion λ(t): 
  
  
    ∂ u∂ t   = ∇· (λ(t) ∇ u) + f(x, y, t),    (x, y, t) ∈Ω× (0, T], 
    subject to initial and Dirichlet boundary conditions:
        u(x, y, 0)    = u^0(x, y),    (x, y) ∈Ω     u = ψ(x, y, t),    (x, y, t) ∈∂Ω× (0, T].

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



 §.§ The composite domain problem

  
  Next, we consider the linear parabolic PDE on a composite domain Ω := Ω_1 ∪Ω_2 (e.g., Figure <ref>), with constant diffusion coefficients (λ_1, λ_2): 
  
  
    ∂ u_1∂ t   = ∇· (λ_1 ∇ u_1) + f_1(x, y, t),    (x, y, t) ∈Ω_1 × (0, T], 
    ∂ u_2∂ t   = ∇· (λ_2 ∇ u_2) + f_2(x, y, t),    (x, y, t) ∈Ω_2 × (0, T], 
    subject to initial conditions:
        u_1(x, y, 0)    = u_1^0(x, y),    (x, y) ∈Ω_1, 
    
        u_2(x, y, 0)    = u_2^0(x, y),    (x, y) ∈Ω_2, 
    Dirichlet boundary conditions:
        u_1    = ψ(x, y, t),    (x, y, t) ∈∂Ω× (0, T], 
    and interface/matching conditions:
        u_1 - u_2    = μ_1(x, y, t),    (x, y, t) ∈Γ× (0, T], 
    λ_1 ∂ u_1/∂ n - λ_2 ∂ u_2/∂ n   = μ_2(x, y, t),    (x, y, t) ∈Γ× (0, T].

  
  In formula (<ref>), ∂ u_s/∂ n, s = 1, 2 denotes the normal derivative at the interface Γ, i.e., ∂ u_s/∂ n = ∇ u_s ·𝐧, where 𝐧 is the outward unit normal vector at the interface Γ.  
  The initial, boundary, and interface data u_1^0(x, y), u_2^0(x, y), ψ(x, y, t), μ_1(x, y, t), and μ_2(x, y, t); the diffusion coefficients (λ_1, λ_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.
  
  
  



 
  
  We consider the circular geometries depicted in Figure <ref> as the geometrical setting for our proposed benchmark problems in this work.  In applications (Section <ref>), other geometries will likely be considered, some much more complicated than Figure <ref>.  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 (<ref>–<ref>).  As such, we think that the geometries in Figure <ref> 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).Comment 1 Reviewer 1



  For both the single and composite domain problems, we could also consider other boundary conditions, e.g., a Neumann boundary condition as in <cit.>, etc. 



  
  
  
  
  
  
  
  
  
  
  
  
  




§ OVERVIEW OF NUMERICAL METHODS
 

  

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


Let Ω_s be covered by a structured triangulation, _s, so that each element T∈_s has some part inside of Ω_s; see Figures <ref> and <ref>.
Here, s = 1,2 is an index for the composite domain problem (<ref>–<ref>), which will be omitted when referring to the single domain problem (<ref>, <ref>).
(For the latter, note that  covers Ω.)
Typically _1 and _2 would be created from a larger mesh by removing some of the cells.Comment 4, Reviewer 1
Further, let _Γ={T∈ : T ∩Γ≠∅} be the set of intersected elements; see Figure <ref>.
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 _Γ 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 Ω_s, consisting of Q_p-elements on the mesh _s:

    V_h^s={ v∈ C^0(Ω_s): vT∈ Q_p(T),   T∈_s}.

For the single domain problem (<ref>, <ref>) we solve for the solution u ∈ V_h; while for the composite domain problem (<ref>–<ref>), we solve for the pair u∈ V_h^1 × V_h^2.
For the latter problem, this means that the degrees of freedom are doubled over elements belonging to _Γ.




We begin by stating the weak formulation for the single domain problem (<ref>, <ref>).

Let (·,·)_X and ··_Y be the L_2 scalar products taken over the two- and one-dimensional domains X⊂^2 and Y⊂^1, respectively.
The present method is based on modifying the weak formulation by using Nitsche's method <cit.> to enforce the boundary condition (<ref>). 
By multiplying (<ref>) with a test function v∈ V_h, and integrating by parts, we obtain: 

    (u̇,v)_Ω+(λ∇ u, ∇ v)_Ω -λunv_Γ
    =(f,v)_Ω,   ∀ v∈ V_h.

Note that (<ref>) is consistent with the following terms:

    γ_D/λ uv_Γ   = γ_D/λψv_Γ,
     
    -uλ∂ v/∂ n_Γ   = -ψλ∂ v/∂ n_Γ,

where γ_D is a constant, and  is the side length of the quadrilaterals in the triangulation. Now, adding (<ref>, <ref>) to (<ref>) gives the following weak form: Find u ∈ V_h such that 

    (u̇,v)_Ω+a(u,v)=L(v),   ∀ v∈ V_h,

where

    a(u,v)     =(λ∇ u, ∇ v)_Ω -λunv_Γ
    -uλvn_Γ+γ_D/λ uv_Γ, 
    
    
    
    L(v)     =(f,v)_Ω
    +λψγ_D/ v- vn_Γ.


For _Γ (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  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 <ref>.
If Ω is moved with respect to  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 (<ref>) – 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  <cit.>. 
Adding stabilization to (<ref>) results in the following weak form: Find u ∈ V_h such that 

    (u̇,v)_Ω+ γ_M j(u̇,v)+a(u,v) + γ_A ^-2λ j(u,v)=L(v),   ∀ v∈ V_h,

where γ_M  and γ_A are scalar constants.








In order to state the definition of stabilization (<ref>), denote by _s the set of faces, as seen in Figures <ref> and <ref>. 
That is, _s is the set of all faces of the elements in _Γ, excluding the boundary faces of _s: 


    _s={F=T_A ∩ T_B   :   T_A ∈_Γ  or   T_B ∈_Γ,    T_A,T_B∈_s}.

Then, the stabilization term is defined as: 

    j_s(u,v)=∑_F∈_s∑_k=1^p ^2k+1/(2k+1)(k!)^2[∂_n^k u][∂_n^k v]_F,

where [u] = u|_F_+-u|_F_-  is the jump over a face, F; n refers to a normal of F;  and ∂_n^k u denotes the k-th order normal derivative. 

The scaling with respect to k of the terms in (<ref>) 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.
Addition to Comment 5, Reviewer 1






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


  
  * For both (<ref>) and (<ref>), multiply the equation for us with a test function  vs∈ V_h^s, and then integrate by parts; 
  
  * Add terms consistent with the interface and boundary conditions; and 
  
  * Add stabilization terms j_1 and j_2 over _1 and _2, respectively.

This results in the following weak formulation for (<ref>–<ref>). Find u=u∈ V_h^1 × V_h^2 such that:

    M(u̇,v)   +A(u,v)+a_Γ(u,v)+a_∂Ω(u,v) 
       =L_Ω(v)+L_Γ(v)+L_∂Ω,   ∀ v=v∈ V_h^1 × V_h^2,
    where the bilinear forms M and A correspond to the stabilized mass and stiffness matrices:
    
    M(u̇,v)   =∑_s=1^2(u̇s,vs)_Ω_s + γ_M j_s(u̇s,vs), 
    
    A(u,v)   =∑ _s=1^2(λ∇us, ∇vs)_Ω_s +γ_A ^-2λ j_s(us,vs); L_Ω corresponds to the forcing function:
    
    L_Ω(v)     =∑ _s=1^2(f_s,vs)_Ω_s; 
    a_Γ and L_Γ consistently enforce the interface conditions (<ref>, <ref>):
    
    a_Γ (u,v)     =-[u]{λvn}_Γ-{λun}[v]_Γ + γ_Γ/ [u][v]_Γ, 
    
    L_Γ(v)     =γ_Γ/μ_1[v]_Γ + κ_1μ_2v2_Γ + κ_2μ_2v1_Γ -μ_1{λvn}_Γ; and the terms a_∂Ω and L_∂Ω enforce the boundary condition (<ref>) along the outer boundary, ∂Ω:
    
    a_∂Ω(u,v)     =-λu1nv1_∂Ω -u1λv1n_∂Ω+γ_D/λu1v1_∂Ω, 
    
    L_∂Ω(v)     =λψγ_D/v1 - v1n_∂Ω.

In (<ref>–<ref>), 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 <cit.>:

    κ_1 = λ2/λ1 + λ2,   κ_2 = λ1/λ1 + λ2,   γ_Γ = γ_D λ1λ2/λ1 + λ2.

The remaining parameters (appearing in Equations <ref>, <ref>, <ref>–<ref>) are given by:

    γ_M=0.75,  γ_A=1.5,   γ_D=5p^2.

The scaling of γ_D with respect to p follows from an inverse inequality. 
and the constants are based on numerical experiments,
When p=1 these reduce to the same parameters as the ones used in <cit.>, 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 <cit.>, where γ_A was investigated numerically.
Comment 5, Reviewer 1

In order to use cut–FEM, one needs a way to perform integration over the intersected elements _Γ.

For example, with the interface problem, on each element K ∈_Γ, we need a quadrature rule for the K∩Ω_1, K∩Ω_2 and K∩Γ.
For the numerical tests in this work (Section <ref>), we represent the geometry by a level set function, and compute high-order accurate quadrature rules with the algorithm from <cit.>.


Optimal (second-order) convergence was rigorously proven for cut–FEM applied to the Poisson problem in <cit.>.
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
 









matrix[1][]
  #1
  -ifnextchar@ifnextchar
  *@̧MaxMatrixCols c






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 <cit.>, and see <cit.> 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 (<ref>) or (<ref>, <ref>) 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 <cit.> and <cit.>.  We also introduce new developments here for handling implicitly-defined geometries.  
(The reader can consult <cit.> 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 Ω 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 <cit.>). 



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 𝒪(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 (<ref>–<ref>).  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 (<ref>, <ref>), after omitting the index s and replacing interface conditions with boundary conditions; see <cit.>.) 

Let Ω_s (s=1,2) be embedded in a rectangular auxiliary domain Ω_s^0. 
Introduce a uniform, Cartesian grid denoted M_s^0 on Ω_s^0, with grid-spacing h_s. 
Let M_s^+ = M_s^0 ∩Ω_s denote the grid points inside each subdomain Ω_s, and 
M_s^- = M_s^0 ∖ M_s^+ the grid points outside each subdomain Ω_s. 
Note that the auxiliary domains Ω_1^0, Ω_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, α, with α = 5, 9, to be the stencil of the standard five-point or a wide nine-point Laplacian, i.e.,  



    N_j, k^s, 5   = { (x_j,y_k), (x_j±1,y_k), (x_j, y_k±1) }  or
    
        N_j, k^s, 9   = { (x_j,y_k), (x_j±1,y_k), (x_j, y_k±1), (x_j±2,y_k), (x_j, y_k±2) }.


Next, with α fixed, 
define the point-sets 


    N_s^0 = ⋃_(x_j, y_k) ∈ M_s^0 N_j, k^s, α,    N_s^+ = ⋃_(x_j, y_k) ∈ M_s^+ N_j, k^s, α,   and   N_s^- = ⋃_(x_j, y_k) ∈ M_s^+ N_j, k^s, α.


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 Ω_s, so that N_s^+ ⊄Ω_s, even though M_s^+ ⊂Ω_s.  (Likewise, N_s^- ⊄Ω∖Ω_s, even though M_s^- ⊂Ω∖Ω_s). 



Lastly, we now define the important point-set 


    γ_s = N_s^+ ⋂ N_s^-,


which we call the discrete grid boundary.  In words, γ_s is the set of grid points that straddle the continuous interface Γ.  (See Figure <ref> 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 (<ref>, <ref>), and then define the Auxiliary Problem. 

Indeed, the discretization we consider is 


    L_Δ t, h^s u_s^i+1 = F_s^i+1,    (x_j, y_k) ∈ M_s^+,


where (i) L_Δ t, h^s u_s^i+1 := λ_s(t^i+1) Δ_h u_s^i+1 - σ u_s^i+1, (ii) Δ_h is either a five- or nine-point Laplacian 
in each subdomain, 
and (iii) σ 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 <ref>. 

The choices of discretization (<ref>) in each subdomain need not be the same.  As in <cit.>, 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. 


 
    At time t^i+1, given the right-hand side grid function q_s^i+1: M_s^0 →ℝ, the following difference equations (<ref>, <ref>) are defined as the discrete AP. 
    
    
    
    
    
        L_Δt, h^s u_s^i+1    = q_s^i+1,          (x_j,y_k)∈M^0_s  

        u_s^i+1    =0,          (x_j,y_k)∈N^0_s\M^0_s  
    





 
    
    
    For a given right-hand side q_s^i + 1, the solution of the discrete AP (<ref>, <ref>) defines a discrete Green's operator G_Δ t, h^s q_s^i + 1. 
    The choice of boundary conditions (<ref>) will affect the resulting grid function G_Δ t, h^s q_s^i + 1, and thus the Boundary Equations with Projection defined below. 
    However, the choice of boundary conditions 
    (<ref>) in the AP 
    will not affect the numerical approximation of (<ref>–<ref>), so long as the discrete AP is uniquely solvable and well-posed.



Let us denote by G_Δ t, h^s F_s^i + 1 the particular solution on N^+_s of the fully-discrete problem (<ref>), defined by solving the AP (<ref>, <ref>) with 


    q_s^i + 1 = 
            F_s^i + 1,     (x_j, y_k) ∈ M_s^+, 
    
            0,     (x_j, y_k) ∈ M_s^-,


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


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


    The Difference Potential of a given density v_γ_s^i + 1 is the grid function P_N_s^+^i + 1 v_γ_s^i + 1 on N_s^+, defined by solving the AP (<ref>, <ref>) with 
    
    
    q_s^i + 1 = 
                0,     (x_j, y_k) ∈ M_s^-, 
    
                L_Δ t, h^s [v_γ_s^i + 1],     (x_j, y_k) ∈ M_s^+,

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


Note that P_N_s^+^i + 1 : 𝐕_γ_s→ N_s^+ is a linear operator on the space 𝐕_γ_s of densities. 
Moreover, the coefficients of P_N_s^+^i + 1 can be computed by solving the AP (Definition <ref>) with the appropriate density v_γ_s^i+1 defined at the points (x_j, y_k) ∈γ_s. 

 
    Given a grid function v_s^i + 1, we denote by _γ_s [v_s^i + 1] the Trace (or Restriction) from N_s^+ to γ_s. 

    Moreover, for a given density v_γ_s^i + 1, denote the trace of the Difference Potential of v_γ_s^i + 1 by P_γ_s v_γ_s^i + 1.  In other words, P_γ_s v_γ_s^i + 1 = _γ_s [P_N_s^+^i + 1 v_γ_s^i + 1]. 








    



    
    



































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

 
At time-level t^i+1, the discrete density u^i+1_γ_s (s=1,2) is the trace of some solution u^i+1_s on domain Ω_s to the Difference Equations (<ref>), i.e., u^i+1_γ_s:=_γ_s[u^i+1_s], if and only if the following BEP holds

    u^i+1_γ_s-P_γ_s^i+1u^i+1_γ_s=_γ_s[G^i+1_Δ t, hF^i+1_s],    (x_j,y_k)∈γ_s,


with _γ_s[·] and P_γ_s defined in Definition <ref>. 



    See <cit.> for the general theory of DPM (including the proof for general elliptic PDE), or one of <cit.> for the proof in the case of parabolic interface problems. 



    A given density v_γ_s^i+1 is the trace of some solution of the fully-discrete finite-difference equations (<ref>) 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_γ_s^i+1.  As originally disucssed in <cit.>, 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 (<ref>–<ref>) at the discrete grid boundary γ_s, using the Cauchy data of (<ref>–<ref>) on the continuous interface Γ, up to the desired second- or fourth-order accuracy.  (By Cauchy data, 
    we mean the trace of the solution of (<ref>–<ref>), together with the trace of its normal derivative, on Γ.) Below, we will define an Extension Operator 
    
    which will extend the Cauchy data of (<ref>–<ref>) 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 (<ref>–<ref>) on Γ.  
    Then, we will use the Extension Operator, together with the BEP (<ref>) and the interface conditions (<ref>, <ref>), 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 (<ref>–<ref>) 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 <ref> 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 


    v_s^i+1(d, ϑ) = v_s^i+1(0, ϑ) + ∑_l = 1^p 1/l!∂^l v_s^i+1(0, ϑ)/∂ n^l d^l,


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 <ref>). 

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


    v_s^i+1(0, ϑ) = u_s^i + 1|_Γ,   ∂ v_s^i+1(0, ϑ)/∂ n = . ∂ u_s^i + 1/∂ n|_Γ,   and  ∂^2 v_s^i+1(0, ϑ)/∂ n^2 = . ∂^2 u_s^i + 1/∂ n^2|_Γ,


where u_s^i + 1 := u_s(x, y, t^i + 1), ∂ u_s^i + 1/∂ n := ∂ u_s(x, y, t^i + 1)/∂ n, etc. 
As a last step, a straightforward sequence of calculations (see Appendix <ref>) shows that 


    ∂^2 u_s^i + 1/∂ n^2≈1/λ_s( 3 u_s^i + 1 - 4 u_s^i + u_s^i - 1/2 Δ t - f^i + 1) - ∂^2 u_s^i + 1/∂ϑ^2 + κ∂ u_s^i+1/∂ n,


where κ denotes the curvature of Γ. 
Therefore, with v_s^i+1(d, ϑ) defined by (<ref>–<ref>), 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 ∂ u_s^i + 1/∂ n. 
The Extension Operator will incorporate the interface conditions (<ref>, <ref>) when it is combined with the BEP (<ref>), so that the only independent unknowns at each time step t^i + 1 will be . (u_1^i + 1, ∂ u_1^i + 1/∂ n) |_Γ or . (u_2^i + 1, ∂ u_2^i + 1/∂ n) |_Γ. 

(This is also true for the fourth-order numerical method – see Appendices <ref> and <ref>.) 

Now we are ready to define the Extension Operator that extends the Cauchy data of (<ref>–<ref>) from Γ to γ_s. 




 
    Let v_s^i + 1(d, ϑ) be defined by (<ref>–<ref>).  Let 𝔲_s, Γ^i + 1 denote the Cauchy data of (<ref>–<ref>) at t^i + 1 on Γ, i.e., 𝔲_s, Γ^i + 1 = . (u_s^i + 1, ∂ u_s^i + 1/∂ n) |_Γ.  The extension operator Ex_s that extends 𝔲_s, Γ^i + 1 from Γ to γ_s is 
    
    
    Ex_s 𝔲_s, Γ^i + 1 := v_s^i + 1(d, ϑ)|_γ_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 𝔲_s, Γ^i + 1.  Indeed, we wish to choose a basis ϕ_ν(ϑ) on Γ (ν = 1, 2, 3, …) in order to accurately approximate the two components of the Cauchy data 𝔲_s, Γ^i + 1.  To be specific, whichever basis we choose, we require that 


    ε_𝒩^0, 𝒩^1(𝔲_s, Γ^i + 1)    =*min_c_1, ν^s, i+1, c_2, ν^s, i+1∫_Γ(| u_s^i + 1 - ∑_ν=1^𝒩^0 c_1, ν^s, i+1ϕ_ν(ϑ)|^2 
                       ⋯ + | ∂ u_s^i + 1/∂ n -∑_ν=1^𝒩^1 c_2, ν^s, i+1ϕ_ν(ϑ)|^2)  dϑ


tends to zero as 𝒩^0, 𝒩^1 →∞, for some sequence of real numbers (c_1, ν^s, i+1)_ν=1^𝒩^0 and (c_2, ν^s, i+1)_ν=1^𝒩^1.  In other words, we require 


    lim_𝒩^0, 𝒩^1 →∞ε_𝒩^0, 𝒩^1(𝔲_s, Γ^i + 1) = 0.



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 <ref>, each function u considered in the test problems on a composite domain (<ref>, <ref>, <ref>) 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 𝔲_1, Γ^i + 1 and 𝔲_2, Γ^i + 1 are smooth, periodic functions of arclength ϑ.  (Note that 𝔲_1, Γ^i + 1 and 𝔲_2, Γ^i + 1 need not agree, and indeed do not – neither μ_1(x, y, t) nor μ_2(x, y, t) in (<ref>, <ref>) 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 


    ϕ_1(ϑ) = 1,   ϕ_2k(ϑ) = cos(2 π k ϑ/|Γ|),   ϕ_2k + 1(ϑ) = sin(2 π k ϑ/|Γ|),


and k > 1. 
Moreover, at every time step t^i + 1, we will discretize the Cauchy data 𝔲_s, Γ^i + 1 = . (u_s^i + 1, ∂ u_s^i + 1/∂ n) |_Γ using this basis.  Therefore, we let 


    𝔲̃_s, Γ^i+1 = ∑_ν=1^𝒩^0 c^0, i+1_1,νΦ^0_ν(ϑ) + ∑_ν=1^𝒩^1 c^1, i+1_2,νΦ^1_ν(ϑ)   and  ũ^i+1_s,Γ≈𝔲^i+1_s,Γ,


where Φ^0_ν = (ϕ_ν,0) and Φ^1_ν = (0,
ϕ_ν) are the set of basis functions used to represent the Cauchy
data on the interface Γ. 


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, <cit.> use this approach to design a high-order accurate numerical method for the Helmholtz equation, in a geometry with a reentrant corner.  Furthermore, <cit.> 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 <ref>, we derive a linear system for the coefficients (c_1, ν^s, i+1)_ν=1^𝒩^0 and (c_2, ν^s, i+1)_ν=1^𝒩^1, by 
combining the interface conditions (<ref>, <ref>), the BEP (<ref>), the Extension Operator (<ref>–<ref>), and the spectral discretization (<ref>). 






Then, the numerical approximation u_s^i+1≈ u_s(x_j, y_k, t^i + 1) of (<ref>–<ref>) at all grid-points (x_j, y_k) ∈ N_s^+ follows directly from the discrete, generalized Green's formula, which we state now. 


    At each time step t^i + 1, the numerical approximation u_s^i+1≈ u_s(x_j, y_k, t^i + 1)|_(x_j, y_k) ∈ N_s^+ of (<ref>–<ref>) is given by 
    
    
    u_s^i + 1 := P_N_s^+^i + 1 u_γ_s^i + 1 + G_Δ t, h^s F_s^i + 1.

    
    Here, u_γ_s^i + 1 = Ex_s 𝔲̃_s, Γ^i + 1, and 𝔲̃_s, Γ^i + 1 is constructed from (i) the solution of 
    the BEP (see Appendix <ref>) and (ii) the spectral 
    discretization (<ref>).  (Recall that P_N_s^+^i + 1 u_γ_s^i + 1 
    is the Difference Potential of the density u_γ_s^i + 1, 
    while G_Δ 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 <cit.> and <cit.> to the composite domain problem (<ref>–<ref>) 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 interpolant F̃(x, y) of F(x,y) on a subset of M^0
 
near the continuous interface Γ. Next, we parameterize Γ by arc-length using numerical quadrature. With this parameterization, we (i) compute the Fourier series expansion from initial conditions for the Cauchy data 
𝔲_s, Γ^i + 1
on the implicit interface Γ, and (ii) construct the extension operators (Definition <ref>) with p=2 or p=4. 

 
    Due to the second or fourth-order accuracy (in both space and time) of the underlying discretization (<ref>), the extension operator (<ref>) 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 
    <cit.> and <cit.>), we expect second- and fourth-order accuracy in the maximum norm for the error in the computed solution (<ref> or <ref>) for both the single and composite domain parabolic problems.  



    Indeed, in the numerical results (Section <ref>) we see that the computed solution (<ref>) at every time level t^i+1 has accuracy 𝒪(h^2 + Δ t^2) for the second-order method, and 𝒪(h^4 + Δ t^4) for the fourth-order method, for both the single and composite domain problems, with explicit or implicit geometry.  See <cit.> 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 Ω_s^0 (s = 1,2) and formulate the Auxiliary Problem (AP; Definition <ref>). 
    
    
  * Step 2: At each time step t^i+1, compute the Particular Solution u_s^i+1 = G^i+1_Δ t, h F^i+1_s, (x_j,y_k)∈ N_s^+, using the AP with the right-hand side (<ref>). 

    
  * Step 3: Construct the matrix in the boundary equations (<ref>) (discussed in Appendix <ref>), derived from the Boundary Equation with Projection (BEP) (<ref>), via several solutions of the AP. (When the diffusion coefficients λ_s are constant, this is done once, as a pre-processing step before the first time step.) 
    
    
  * Step 4: Compute the approximation of the density u^i+1_γ_s, by applying the Extension Operator (<ref>) to the solution of (<ref>). 
    
    
  * Step 5: Construct the Difference Potentials P_N_s^+γ_s u^i+1_γ_s of the density u_γ_s^i + 1, using the AP with the right-hand side (<ref>). 

    
  * Step 6: Compute the numerical approximation u_s^i+1≈ u_s(x_j, y_k, t^i + 1) of the PDE (<ref>–<ref>) using the discrete, generalized Green's formula (<ref>). 





  

 §.§ SBP–SAT–FD
 



We continue in this section with a brief presentation of SBP–SAT–FD, for solving the parabolic problems presented in Section <ref>. For more detailed discussions of the SBP–SAT–FD method, we refer the reader to two review papers <cit.>. 



In SBP–SAT–FD, to resolve geometrical features, the physical domain is mapped to a reference domain with a simple geometry, e.g., the unit square.The SBP–SAT–FD method was originally used on Cartesian grids. To resolve complex geometries, we consider a grid mapping approach by transfinite interpolation <cit.>. Comment 10 of reviewer 1 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 (<ref>, <ref>), shown in Figure <ref>, 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  (<ref>–<ref>) is divided into nine subdomains, as shown in Figure <ref>.  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 <cit.>. Comment 12 of reviewer 1











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 

    u_t = u_xx+u_yy,   (x,y)∈Ω

in a physical domain will be transformed to 

    Ju_t = (α u_ξ)_ξ+(β u_η)_ξ+(β u_ξ)_η+(γ u_η)_η,    (ξ,η)∈ [0,1]^2,

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 (<ref>) and the derivation of the grid transformation are presented in Section 3.2 of <cit.>.Comment 3 from reviewer 2 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 (<ref>). 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 𝐮,𝐯, respectively.  Integration by parts states:

    ∫_0^1 u_xv dx= uv|_0^1 - ∫_0^1 uv_x dx.

The SBP operator D_1 mimics integration by parts:
 
    (D_1 𝐮)^T H 𝐯= 𝐮^T B𝐯 - 𝐮^T H D_1𝐯,

 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 <cit.>. It is easy to verify that (<ref>) is equivalent to 
  
    D_1^TH+HD_1=B,

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 (<ref>). 

The SBP operators D_1 were first constructed in <cit.> and later revisited in <cit.>. 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 <cit.>. 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 <cit.>.




 
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 <cit.>, 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 <cit.> and the discontinuous Galerkin method <cit.>. 

We note that in <cit.>, SBP–SAT–FD methods were developed for the wave equation 
 
    Jv_tt = (av_ξ)_ξ+(bv_η)_ξ+(bv_ξ)_η+(cv_η)_η,   (ξ,η)∈ [0,1]^2,

with Dirichlet boundary conditions, Neumann boundary conditions, and interface conditions.  Comparing equation (<ref>) with (<ref>), 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 (<ref>) and (<ref>) are the same. 
 
 Assuming homogeneous boundary data for simplified notation, we write the SBP–SAT–FD discretization of (<ref>) as 
 
    𝐯_tt = Q𝐯,

 where Q is the spatial discretization operator including the boundary implementation. For the scheme developed in <cit.>, stability is proved by the energy method by multiplying (<ref>) by 𝐯_t^TH_2 from the left, 
  
    𝐯^T_tH_2𝐯_tt = 𝐯^T_tH_2Q𝐯,

  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  <cit.> that H_2Q is symmetric and negative semi-definite. Therefore, we can write (<ref>) as 
    
    d/d t (𝐯^T_tH_2𝐯_t - 𝐯^TH_2Q𝐯)=0,

 where the discrete energy, 𝐯^T_tH_2𝐯_t - 𝐯^TH_2Q𝐯,  for (<ref>) is conserved. 

 
If we use the same operator Q to discretize the heat equation (<ref>) with the same boundary condition as the wave equation (<ref>), then the scheme
 
    𝐯_t = Q𝐯,

 is also stable. To see this, we multiply (<ref>) by 𝐯^TH_2 from the left, and obtain

    d/d t (𝐯^TH_2𝐯) = 𝐯^TH_2Q𝐯≤ 0,

 where 𝐯^TH_2𝐯 is the discrete energy for (<ref>). In this paper, we use the spatial discretization operators developed in <cit.> to solve both the single (<ref>, <ref>) and composite domain problems (<ref>–<ref>).
 
In <cit.>, 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 (<ref>), 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 <ref>), and then briefly motivate and discuss these choices (in Section <ref>).  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 (<ref>, <ref>) (for the single domain problem) or (<ref>–<ref>) (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. 



 §.§ List of test problems
 



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

  
    
  * Constant diffusion (Test Problem 1A; TP–1A):  Consider the PDE (<ref>, <ref>), with λ(t) ≡ 1, Ω = { (x, y) ∈^2 : x^2 + y^2 ≤ 1 }, and the final time T=1.0.  Then, TP–1A (adapted from <cit.>), is given by 
    
    
    u(x, y, t) = x^9 y^8 e^-t. TP–1A

    

    
  * Time-varying diffusion (Test Problem 3A; TP–3A):  Same as TP–1A, but with diffusion coefficient 
    
    λ(t) = 11/10 + sin (10 π t). TP–3A
 

  



  
  * 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 (<ref>–<ref>), with Ω = [-2, 2] × [-2, 2], Ω_2 = { (x, y) ∈^2 : x^2 + y^2 ≤ 1 }, Ω_1 = Ω∖Ω_2, Γ = { (x, y) ∈^2 : x^2 + y^2 = 1 }, and the final time T = 1.0. 
  
  
    
  * (Test Problem 2A; TP–2A):  A modified version of the test adapted from <cit.>.  Let (λ_1, λ_2) = (10, 1), and 
    
    
    u(x, y, t) = 
            e^-tsin x cos y,     (x, y) ∈Ω_1, 
    
            e^-t (x^2 - y^2),     (x, y) ∈Ω_2. 
          TP–2A


    
  * High-frequency oscillations (Test Problem 2B; TP–2B):  A modified version of the test adapted from <cit.>.  Let (λ_1, λ_2) = (10, 1), and 
    
    
    u(x, y, t) = 
            e^-tsin (3 π x) cos (7 π y),     (x, y) ∈Ω_1, 
    
            e^-t (x^2 - y^2),     (x, y) ∈Ω_2. 
          TP–2B


    
  * Large contrast in diffusion coefficients, and large jumps in both solution and flux at interface (Test Problem 2C; TP–2C):  Let (λ_1, λ_2) = (1000, 1), and 
    
    
    u(x, y, t) = 
            0,     (x, y) ∈Ω_1, 
    
            1000 sin(10 t) x^4 y^5,     (x, y) ∈Ω_2. 
          TP–2C

    

  
  
  




 §.§ Motivation of the chosen test problems
 







Test Problem 1A (<ref>) 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, (<ref>) adds on (incrementally) the complication of time-varying diffusion. 

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

Then, (<ref>) adds additional challenges onto (<ref>) in the form of much higher-frequency oscillations; while (<ref>) adds onto (<ref>) in the form of both (i) large contrast in diffusion, and (ii) large jumps in the analytical solution and its flux.








§ NUMERICAL RESULTS
 

  

 §.§ Time discretization
 


The spatial discretization for each method is discussed in Section <ref>. 
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 


    Δ t=0.5h.


However, note that h in (<ref>) 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 (<ref>)); and for SBP–SAT–FD, h is the minimum grid spacing in the reference domain.

  

 §.§ Measure for comparison
 




Let u_j, k^i denote the computed numerical approximation of u(x, y, t) at the grid-point (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 (<ref>, <ref>), the maximum error is computed as:


    E     := max_t^i ∈ (0, T]max_(x_j, y_k) ∈Ωu(x_j, y_k, t^i) - u_j, k^i, 
    and for the composite domain problem (<ref>–<ref>) as:
    
    
    E     := max_t^i ∈ (0, T]max_(x_j, y_k) ∈Ω_1 ∪Ω_2u(x_j, y_k, t^i) - u_j, k^i.



    

    

    

  

 §.§ Convergence results
 

    
    




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



In Tables <ref>–<ref>, 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 


    ρ_n = log(E_n - 1 / E_n)/log(DOF_n - 1 / DOF_n),


which is the estimated rate of convergence, denoted in Tables <ref>–<ref> by “Rate”. 
 

In Figures <ref>, <ref>, <ref>–<ref>, 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 <ref>–<ref> 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 (<ref>, <ref>) are similar, up to a constant factor; while for the composite-domain problems (<ref>, <ref>, <ref>) we do see differences of one or two orders of magnitude, with the DPM having the smallest errors. 





In Table <ref> and Figure <ref>, we observe that the measured rates of convergence for the numerical approximation of Test Problem 1A (<ref>) 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. 






Table <ref> and Figure <ref> show the numerical results for (<ref>).  This test problem has the same manufactured solution as (<ref>), 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 (<ref>).  This similarity in the numerical results demonstrates that the three methods can robustly handle time-varying diffusion coefficients.Comment 16 Reviewer 1





The plots of spatial error at the final time T = 1.0, shown in Figure <ref>, 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 (<ref>) are shown in Figure <ref>.  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 <ref>, 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 Ω, where the grid is most stretched.  








Regarding the max-norm error in presented in Table <ref> and Figure <ref>, 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 <ref>, 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. Comment 15 of reviewer 1








In Table <ref> and Figure <ref>, we see the numerical results for (<ref>).  The analytical solution is similar to (<ref>), though much more oscillatory – this additional challenge is manifested by an increase in error by several orders of magnitude. 







In Table <ref> and Figure <ref>, we see the numerical results for (<ref>), 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 <ref>–<ref> are almost identical, which demonstrates the robustness and flexibility of the DPM.




§ DISCUSSION
 
  


There are many possible methods (Section <ref>) 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 (<ref>, <ref> or <ref>–<ref>).  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 Figures <ref>, <ref>. 
(We experimented with different stabilization parameters in the numerical tests in Section <ref> and have seen that the error can be reduced up to a factor of 15, so large gains in accuracy are possible.)(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.)Comment 6, Reviewer 1
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 equations (<ref>) and the Boundary Equations with Projection (Thm. <ref>).  The formulation outlined in Section <ref> allows for an immersed (non-conforming) grid; fast 𝒪(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 <ref> 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 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 (<ref>) 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 <ref>), which incorporates higher-order normal derivatives at the boundary/interface Γ.

Returning to Section <ref>, we see (in Tables <ref>–<ref> and Figures <ref>–<ref>) 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 Figure <ref>, 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 <ref>). 

As seen in the error plots (Figures <ref>, <ref>), 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 <ref> 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 <cit.>,
(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 <ref>, 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 <ref>.
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 <ref> 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 (Section <ref>). 




 §.§ Fully-discrete formulation of (<ref>, <ref>).



The fully-discrete, finite-difference discretization introduced in (<ref>) is

    L^s_Δ t,hu_s^i+1=F_s^i+1,   (x_j,y_k)∈ M^+_s.

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

    Δ_h u_j,k=1/12h^2(-u_j-2,k   +16u_j-1,k+16u_j+1,k-u_j+2,k-u_j,k-2
       +16u_j,k-1+16u_j,k+1-u_j,k+2-60u_j,k)

for points sufficiently far away from the boundary of the auxiliary domain Ω_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 

    Δ_h u_1,1=1/12h^2(   10u_0,1-4u_2,1+14u_3,1-6u_4,1+u_5,1
    
             +   10u_1,0-4u_1,2+14u_1,3-6u_1,4+u_1,5-30u_1,1),

where u_0,1 and u_1,0 will be from the boundary condition (<ref>).

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

    F^i+1_s=-f^i+1_s-σ/3(4u_s^i-u^i-1_s),

and for BDF4–DPM4 by

    F^i+1_s=-f^i+1_s-σ/25(48u_s^i-36u^i-1_s+16u^i-2_s-3u^i-3_s).

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.




 §.§ Equation-based extension



Let us now expand the discussion surrounding (<ref>–<ref>) leading up to Definition <ref> of the Extension Operator (<ref>). 

An important step in this discussion is to recast the original PDE (<ref>, <ref>) into a curvilinear form, for points (x, y) in the vicinity of Γ.  Following the notation <cit.>, let us first introduce the coordinate system (d, ϑ) for points in the vicinity of Γ.  Recall from Definition <ref> 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 (<ref>, <ref>) becomes 




    ∂ u_s/∂ t - λ_s(1/H_ϑ[∂/∂ n(H_ϑ∂ u_s/∂ n)+∂/∂ϑ(1/H_ϑ∂ u_s/∂ϑ)]) =f_s,


where where H_ϑ=1 - d κ is the Lamé coefficient, and κ is the signed curvature along the interface Γ. 


From (<ref>), a straightforward calculation gives the second-order normal derivative ∂^2 u_s/∂ n^2 (used in the calculation of (<ref>)), which is  

    ∂^2 u_s/∂ n^2 = 1/λ_s(∂ u_s/∂ t-f_s)-∂^2u_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 (<ref>) with respect to n, we see that 

    ∂^3 u_s/∂ n^3=1/λ_s(∂^2u_s/∂ t∂ n-∂ f_s/∂ n)-∂^3u_s/∂ n∂ϑ^2+κ∂^2 u_s/∂ n^2


and

    ∂^4u_s/∂ n^4   =1/λ_s^2(∂^2u_s/∂ t^2-∂ f_s/∂ t) +1/λ_s(-2∂^3u_s/∂ t∂ϑ^2 +κ∂^2u_s/∂ n∂ t. 
               . ⋯ -∂^2f_s/∂ n^2+∂^2f_s/∂ϑ^2+∂^4u_s/∂ϑ^4-κ∂^3u_s/∂ n∂ϑ^2)+κ∂^3 u_s/∂ n^3.


Next, let us follow-up on comments made in the text following (<ref>).  There, it was pointed out that the unknown Dirichlet and Neumann data (u_s,∂ u_s/∂ n) are the only data required for the Extension Operator (<ref>) 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 ∂ u_s/∂ t, ∂^2u_s/∂ t^2, ∂ f_s/∂ t, ∂^2 u_s/∂ n∂ t, and ∂^3u_s/∂ t∂ϑ^2 can be approximated by the backward difference formula. In BDF2–DPM2,

    ∂ u_s^i+1/∂ t   ≈3u_s^i+1-4u_s^i+u_s^i-1/2Δ t   and 
    ∂^2 u_s^i+1/∂ t^2   ≈2u_s^i+1-5u_s^i+4u_s^i-1-u_s^i-2/Δ t^2,

while in BDF4–DPM4,

    ∂ u_s^i+1/∂ t   ≈25u_s^i+1-48u_s^i+36u_s^i-1-16u_s^i-2+3u_s^i-3/12Δ t   and 
    ∂^2 u_s^i+1/∂ t^2   ≈1/Δ t^2(15/4u_s^i+1-77/6u_s^i + 107/6u_s^i-1. 
                   . ⋯-13u_s^i-2+61/12u_s^i-3-5/6u_s^i-4).



  * The derivatives in terms of arclength ϑ can be computed from u_s or ∂ u_s/∂ n. For example, denoting u_s=∑_ν=1^𝒩^0 c_1,ν^i+1ϕ_ν(ϑ) (using notation following from (<ref>)), then it comes handy that

    ∂ u_s/∂ϑ   =∑_ν=1^𝒩^0 c_1,ν^s, i+1ϕ_ν'(ϑ), 
    ∂^2 u_s/∂ϑ^2   =∑_ν=1^𝒩^0 c_1,ν^s, i+1ϕ_ν”(ϑ),   and
    ∂^4 u_s/∂ϑ^4   =∑_ν=1^𝒩^0 c_1,ν^s, i+1ϕ_ν^(4)(ϑ).


 







 §.§ The system of equations at each time step.
  


With the Cauchy data 𝔲_s, Γ^i + 1 and Extension Operator Ex_s 𝔲_s, Γ^i + 1 from Γ to γ_s introduced in Definition <ref>, and the spectral representation introduced in (<ref>), we now give a sketch of the linear system for the coefficients (c_1, ν^s, i+1)_ν=1^𝒩^0 and (c_2, ν^s, i+1)_ν=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 𝔲̃_s, Γ^i + 1 (<ref>, <ref>) into the BEP (<ref>), the resulting linear systems are 


    ∑_k = 0^1 ∑_ν=1^𝒩^k( c_k, j^s, i + 1Ex_s Φ^k_ν - c_k, j^s, i + 1 P_γ_s^i + 1Ex_s Φ^k_ν) = _γ_s [G_Δ t, h^i+1 F_s^i + 1].


This can be further elucidated by introducing the vector of unknowns 


    𝐜_s^i + 1 = [
            c_1, 1^s, i + 1   c_1, 2^s, i + 1  ⋯   c_1, 𝒩^0^s, i + 1_𝐜_s, 1^i + 1  c_2, 1^s, i + 1   c_2, 2^s, i + 1  ⋯   c_2, 𝒩^1^s, i + 1_𝐜_s, 2^i + 1]^⊤


(so that 𝐜_s^i + 1 = [ 𝐜_s,1^i + 1, 𝐜_s,2^i + 1 ]^⊤), and the matrix 


    A_s = [    (I - P_γ_s^i + 1) Ex_s Φ^0_1,    (I - P_γ_s^i + 1) Ex_s Φ^0_2,   ⋯   (I - P_γ_s^i + 1) Ex_s Φ^0_𝒩^0_A_s, 1, 
       (I - P_γ_s^i + 1) Ex_s Φ^1_1,    (I - P_γ_s^i + 1) Ex_s Φ^1_2,   ⋯   (I - P_γ_s^i + 1) Ex_s Φ^1_𝒩^1]_A_s, 2.


Then, the full system of equations (<ref>) is 


    A [ [1.5]
            𝐜_1^i + 1;               𝐜_2^i + 1 ] = [ [1.5]
            _γ_1 [G_Δ t, h^i + 1 F^i + 1_1];               _γ_2 [G_Δ t, h^i + 1 F^i + 1_2] ],   with  A = [ [1.5]
            A_1                 0;                 0               A_2 ].


However, note that 𝐜_1^i + 1 and 𝐜_2^i + 1 are related by the interface conditions (<ref>, <ref>), so that the number of unknowns in (<ref>) is equal the dimension of either 𝐜_1^i + 1 or 𝐜_2^i + 1, depending on which one is considered the independent unknown. 











Therefore, the dimension of A is (|γ_1| + |γ_2|) × (𝒩^0 + 𝒩^1), where 𝒩^0 + 𝒩^1 is the dimension of 𝐜_1^i + 1 or 𝐜_2^i + 1 (whichever is the independent unknown). 


    The independent unknown (𝐜_1^i + 1 or 𝐜_2^i + 1) is chosen so that the finite-dimensional, spectral representation (<ref>) of the Cauchy data 𝔲_s,Γ^i+1 accurately resolves the Cauchy data with a small number of basis functions, in the consideration of both accuracy and computational efficiency. For (<ref>) and (<ref>), we choose 𝐜_2^i + 1 as the independent unknown, while for (<ref>) we choose 𝐜_1^i + 1.  With these choices for the independent unknown, we have 𝒩^0=𝒩^1=1 for the three considered test problems.


Since each column involves the Difference Potentials operator P_γ_s^i+1 applied to a vector Ex_s Φ_ν^k, each column is therefore constructed via one solution of the Auxiliary Problem (Definition <ref>).  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. 





The authors are grateful to the anonymous referees for their valuable remarks and questions, which led to significant improvements to the manuscript.  We 
gratefully acknowledge the support of 

the Swedish Research Council (Grant No. 2014-6088); the Swedish Foundation for International Cooperation in Research and Higher Education (Grant No. STINT-IB2016-6512); Uppsala University, Department of Information Technology; and the University of Utah, Department of Mathematics. 

Y. Epshteyn, K. R. Steffen, and Q. Xia also acknowledge partial support of Simons Foundation Grant No. 415673.


plainurl      
   




