Immersed Boundary Smooth Extension: A high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods

# Immersed Boundary Smooth Extension: A high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods

David B. Stein Robert D. Guy Becca Thomases Department of Mathematics, University of California, Davis, Davis, CA 95616-5270, USA
###### Abstract

The Immersed Boundary method is a simple, efficient, and robust numerical scheme for solving PDE in general domains, yet it only achieves first-order spatial accuracy near embedded boundaries. In this paper, we introduce a new high-order numerical method which we call the Immersed Boundary Smooth Extension (IBSE) method. The IBSE method achieves high-order accuracy by smoothly extending the unknown solution of the PDE from a given smooth domain to a larger computational domain, enabling the use of simple Cartesian-grid discretizations (e.g. Fourier spectral methods). The method preserves much of the flexibility and robustness of the original IB method. In particular, it requires minimal geometric information to describe the boundary and relies only on convolution with regularized delta-functions to communicate information between the computational grid and the boundary. We present a fast algorithm for solving elliptic equations, which forms the basis for simple, high-order implicit-time methods for parabolic PDE and implicit-explicit methods for related nonlinear PDE. We apply the IBSE method to solve the Poisson, heat, Burgers’, and Fitzhugh-Nagumo equations, and demonstrate fourth-order pointwise convergence for Dirichlet problems and third-order pointwise convergence for Neumann problems.

###### keywords:
Embedded boundary, Immersed Boundary, Fourier spectral method, Complex geometry, Partial Differential Equations, High-order
\newdefinition

remarkRemark

## 1 Introduction

The Immersed Boundary (IB) method was originally developed for the study of moving, deformable structures immersed in a fluid, and it has been widely applied to such problems since its introduction [1, 2, 3]. Recently, the method has been adapted to more general fluid-structure problems, including the motion of rigid bodies immersed in a fluid [4] and fluid flow through a domain with either stationary boundaries or boundaries with prescribed motion [5, 6]. In this broadened context, we use the term Immersed Boundary method to refer only to methods in which (i) the boundary is treated as a Lagrangian structure embedded in a geometrically simple domain, (ii) the background PDE (e.g. the Navier-Stokes equations) are solved on a Cartesian grid everywhere in that domain, and (iii) all communication between the Lagrangian structure and the underlying PDE is mediated only by convolutions with regularized -functions. These methods have many desirable properties: they make use of robust and efficient Cartesian-grid methods for solving the underlying PDE, are flexible to a wide range of problems, and are simple to implement, requiring minimal geometric information and processing to describe the boundary.

The IB method belongs to the broad category of methods known as embedded boundary (EB) methods, including the Immersed Interface [7], Ghost Fluid [8], and Volume Penalty methods [9]. These methods share a common feature: they enable solutions to PDE on nontrivial domains to be computed using efficient and robust structured-grid discretizations; yet these methods differ largely in how boundary conditions are enforced and whether or not the solution is produced in the entirety of a simple domain. Methods which compute the solution everywhere in a -dimensional rectangle admit the simplest discretizations and enable the use of high-order discretizations such as Fourier spectral methods. Unfortunately, this simplicity comes coupled with a fundamental difficulty: the analytic solution to these problems is rarely globally smooth on the entire domain. Consider the one-dimensional Poisson problem on the periodic interval with Dirichlet boundary conditions for . Even if , the solution will typically display jumps in its derivative at the values and (see Figure 2). The lack of regularity in the analytic problem leads to low-order convergence in many numerical schemes, including the Immersed Boundary method.

The advantages of EB methods are substantial enough that significant effort has been expended on improving their accuracy [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. Two different approaches are generally taken. The first approach involves locally altering the discretization of the PDE in the vicinity of the boundary to accommodate the lack of smoothness in the solution. One example of this approach is the Immersed Interface method [7]. Such approaches are particularly useful for interface problems where the solution is required on both sides of the embedded boundary. When the solution is only needed on one side of the interface, a second approach may be taken in which variables are redefined outside of the domain of interest to obtain higher global regularity. Improved convergence rates are achieved as a natural consequence of the properties of the discretization scheme when applied to smooth problems. Variations of this basic idea have been used by the Fourier Continuation (FC) [24, 25, 23] and the Active Penalty (AP) [26] methods to provide high order solutions to PDE on general domains.

In this paper, we introduce a new method termed the Immersed Boundary Smooth Extension (IBSE) method. This method is a first step in remedying two deficiencies of Immersed Boundary methods:

1. For generic problems, the IB method produces discretizations that are only first-order accurate in the vicinity of domain boundaries (or fluid-structure interfaces) [27].

2. The IB method is only able to specify Dirichlet boundary conditions (no-slip, for fluid problems) without specialized interpolation schemes subordinate to the geometry [28]. Although this is not of interest when solving traditional IB fluid-interface problems, it allows the IB method to be generalized for solving other PDE (i.e. reaction-diffusion equations).

In this paper we restrict our attention to PDE set on stationary domains, and consider only PDE without a global divergence constraint. Direct-forcing IB methods produce solutions to PDE that are , with jumps in the normal-derivative of the solution across the boundary, and converge at first order in the grid-spacing [4, 5]. Drawing on ideas from the AP and FC methods, we use Fourier spectral methods to obtain a highly accurate discretization, while adding a volumetric forcing to non-physical portions of the computational domain to force the solution to be globally . We will use the shorthand IBSE- to refer to our method when we need to explicitly denote the global regularity of the solution that is enforced. High-order accuracy is achieved naturally, as a simple consequence of the convergence properties of the Fourier transform.

The IBSE method retains the essential robustness and simplicity of the original IB method. All communication between the Lagrangian boundary and the underlying Cartesian grid is achieved by convolution with regularized -functions or normal derivatives of those -functions. This allows an absolute minimum of geometric information to be used. In the traditional IB method, only the position of the Lagrangian structure must be known; the IBSE method additionally requires normals to that structure and an indicator variable denoting whether Cartesian grid points lie inside or outside of the physical domain where the PDE is defined. Additionally, since normal derivatives can be accurately approximated by convolution with derivatives of regularized -functions, Neumann and Robin boundary conditions can be imposed in the same way that Dirichlet boundary conditions are in direct-forcing IB approaches.

In contrast to other approaches based on the smooth extension of the forcing function or the solution [24, 25, 23, 26, 19], we do not extend the forcing function or solution from a previous timestep. Instead, we smoothly extend the unknown solution to the entire computational domain. This approach directly enforces smoothness of the solution, and it allows for high-order implicit-time discretizations for parabolic equations and implicit-explicit discretizations for many nonlinear PDE. Remarkably, the coupled problem for the solution to an elliptic PDE and that solution’s smooth extension can be reduced to the solution of a relatively small dense system of equations (with size a small multiple of the number of points used to discretize the boundary), along with several FFTs. The dense system depends only on the boundary and the discretization, and so it can be formed and prefactored to allow for efficient time-stepping of parabolic equations.

This paper is organized as follows. In Section 2, we introduce the method, ignoring the details of the numerical implementation. The fundamental contributions of this paper are contained in the modification to the IB discretization that enforces smoothness of the solution, and in the system of equations (given in Section 2.2) that allows for the simultaneous solution of the PDE along with the smooth extension of that unknown solution. In Section 3, we discuss the particular numerical implementation that we choose, detailing a fast algorithm for solving elliptic equations. In Section 4, we compute solutions to the Poisson problem in one and two dimensions, demonstrating high-order pointwise convergence: up to fourth-order for Dirichlet problems and third-order for Neumann problems. In Section 5, we discuss discretization of the heat equation, detailing a fast algorithm for implicit timestepping, and demonstrate fourth-order convergence in space and time. These numerical tests include direct comparisons to the Fourier Continuation [24] and Active Penalty [26] methods. Finally, in Section 6, we solve two nonlinear problems: a 2D Burgers’ equation with homogeneous Dirichlet boundary conditions, and the Fitzhugh-Nagumo equations with homogeneous Neumann boundary conditions. We demonstrate fourth-order convergence in space and time for Burgers’ equation and third-order convergence for the Fitzhugh-Nagumo equations.

## 2 Methods

Let , , and . We will assume that is smooth and does not self-intersect. Two typical domains are shown in Figure 1. We refer to as the physical domain, as the extension domain, and as the computational domain. We first consider the Poisson problem with Dirichlet boundary conditions: \cref@addtoresetequationparentequation

 Δu =f in Ω, (1a) u =g on Γ. (1b)

For now, we assume that and are smooth () functions defined in and respectively. Let be a smooth extension of to , that is, and for all . The direct forcing formulation of the Immersed Boundary method provides a way to solve Equation 1: \cref@addtoresetequationparentequation

 Δu(x)+∫ΓG(s)δ(x−s)ds =fe(x) in C, (2a) ∫Cu(x)δ(x−s)dx =g(s) in Γ. (2b)

Equation 1a is solved in the entire computational domain , while the boundary condition is enforced by the addition of a singular force supported on the boundary that acts as a Lagrange multiplier. This singular force term leads to jumps in the normal derivative of the solution across the boundary . This lack of smoothness restricts Immersed Boundary formulations such as Equation 2 to first-order convergence in the vicinity of the boundary unless one-sided discrete -functions are used, which make the implementation of the method more difficult [10, 27].

We define an example one-dimensional problem to illustrate the limited regularity of solutions produced by Equation 2. Let the computational domain be , where the one-dimensional torus is identified with the periodic interval , and let the physical domain be given by the interval . The extension domain for this problem is . On this domain, we will solve the problem: \cref@addtoresetequationparentequation

 Δu =sinx in Ω=T∖[3,4], (3a) u =0 on Γ={3,4}. (3b)

As is chosen to be in , we may choose for all . The analytic solution to this problem is shown in Figure 2, along with the solution to the extended problem given by Equation 2. Note that the solution is only continuous despite the fact that . Choosing to be smooth leads to low regularity in the solution when the solution to the associated homogeneous problem is nontrivial, as is the case here. Without modification, direct discretizations of this problem will exhibit slow convergence due to the limited regularity of in the vicinity of the boundary.

Rather than choose a smooth extension to the forcing function , we choose an extension of the forcing function that gives a smooth solution on the entire computational domain. Let be any smooth extension of the unknown solution into . We can compute a forcing function associated with :

 Fe=Δue. (4)

The extended forcing function is then defined to be

 ~fe=χΩf+χEFe, (5)

where denotes the characteristic function of the domain . Figure 2 shows the extended forcing function , along with the associated smooth solution to Equation 3. Since the solution is smooth, we expect that standard discretizations of Equation 2 will converge more rapidly when computed using the extension than when using the extension . We emphasize that the extended forcing function depends on the unknown solution .

This motivates us to define a reformulated version of Equation 2, with the extended forcing function given by as defined in Equation 5: \cref@addtoresetequationparentequation

 Δu(x)−(χEΔEku)(x)+∫ΓG(s)δ(x−s)ds =χΩf(x) in C, (6a) ∫Cu(x)δ(x−s)dx =g(s) in Γ. (6b)

Here is any smooth extension operator that satisfies \cref@addtoresetequationparentequation

 Ek:C0(C)∩Ck(Ω) →Ck(C), (7a) (Eku)(x) =u(x)∀x∈Ω. (7b)

When restricted to the physical domain , Equation 6a reduces to , the problem of physical interest. When restricted to the extension domain , Equation 6a reduces to . As long as , then Equation 7b ensures that , and hence in . Again by Equation 7b, in , and so in . Since by Equation 7a, then . This additional regularity of the solution allows standard discretizations of Equation 6 to converge at a faster rate than discretizations of Equation 2. Using and a simple Fourier-spectral discretization, numerical solutions to Equation 6 should converge at [29].

In the remainder of this section, we lay out the remaining components needed to fully specify the IBSE method that do not depend on the particular choice of numerical implementation. In Section 2.1, we discuss how to smoothly extend a known function from to . In Section 2.2, we outline a system of equations that allows us to solve for and its smooth extension simultaneously. A numerical implementation of this method will be discussed in Section 3.

### 2.1 Smooth extension of a known function from Ω to C

One way to construct a smooth extension to a function is to solve a high-order PDE. Compared to the localized and effectively one-dimensional extension strategies taken in [24, 26], the decision to extend a function by solving a fully -dimensional PDE may appear to be needlessly complex. However, we will show in Sections 3 and 2.2 that this choice leads to a robust method that is straightforward to implement and requires relatively low computational cost.

In particular, to compute a extension to a function , we solve the following equation in the extension domain : \cref@addtoresetequationparentequation

 Hkξ =0 in E, (8a) ∂jξ∂nj =∂jv∂nj on Γ, ∀0≤j≤k. (8b)

Here denotes the normal derivative of on the boundary ; a computational formula for is given later in Equation 22. A simple choice of the operator is the polyharmonic operator . From an analytic perspective, this choice is sufficient, although we will show in Section 3.3 that other choices of will produce better results due to numerical issues. The smooth extension may then be constructed as

 ζ=χΩv+χEξ. (9)
{remark}

Our choice of smoothed forcing only depends on in the extension region (see Equation 5), so there is no need to compute the actual extension since the values of in are irrelevant. For our purposes, we can think of an ‘extension’ as simply a function in the computational domain that shares its first normal derivatives with on the boundary , i.e. the function .

As with Equation 1, we may solve Equation 8 in by extending the problem to all of . This can be accomplished by adding singular forces supported along the boundary: \cref@addtoresetequationparentequation

 Hkξ(x)+k∑j=0(−1)j∫ΓFj(s)∂jδ(x−s)∂njds =0 ∀x∈Td, (10a) (−1)j∫Cξ(x)∂jδ(x−s)∂njdx =∂jv∂nj(s) ∀s∈Γ, 0≤j≤k. (10b)

The integral on the left-hand side of Equation 10b is notation for the action of the distribution on the smooth function . The boundary condition given in Equation 10b forces the solution to be , so there is no need to alter this formulation to provide additional regularity.

For notational convenience, we define the operators and by: \cref@addtoresetequationparentequation

 TkF(x) =k∑j=0(−1)j∫ΓFj(s)∂jδ(x−s)∂njds, (11a) T∗kξ(s) =(∫Cξ(s)δ(x−s)ds∫Cξ(s)∂δ(x−s)∂n⋯∫Cξ(s)∂δk(x−s)∂nk)⊺. (11b)

In the language of the Immersed Boundary method, is an interpolation operator; is a spread operator. The notation is suggestive of the fact that these operators obey the adjoint property

 ⟨u,TkF⟩C=⟨T∗ku,F⟩Γ (12)

for all sufficiently smooth and , where the notation denotes the inner product on . Using the operators and , we may represent Equation 10 as

 (HkTkT∗k)(ξF)=(0T∗kv). (13)

### 2.2 A coupled system of equations for u and its extension

Let the spread () and interpolation () operators be defined as \cref@addtoresetequationparentequation

 SG(x) =∫ΓG(s)δ(x−s)ds, (14a) S∗ξ(s) =∫Cξ(x)δ(x−s)dx,∀s∈Γ. (14b)

Notice that and . We can now represent Equation 6 simply as \cref@addtoresetequationparentequation

 Δu−χEΔξ+SG =χΩf, (15a) S∗u =g, (15b)

where is defined by the extension equation

 Hkξ=0, (16)

along with the constraint that

 T∗kξ=T∗ku. (17)

Equations 17, 16 and 15 can be combined into one system of equations:

 ⎛⎜ ⎜ ⎜ ⎜⎝Δ−χEΔSHkTk−T∗kT∗kS∗⎞⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜⎝uξFG⎞⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝χΩf00g⎞⎟ ⎟ ⎟⎠. (18)

Equation 18 is the system of equations that allows the IBSE method to solve for and its extension simultaneously. The remainder of this paper will be concerned with the discretization and inversion of Equation 18 and numerical studies demonstrating the accuracy of solutions produced by those discretizations. We remark that Equation 18 is equally valid for operators of the Helmholtz type () that appear in the discretization of parabolic equations, which will be discussed in Section 5.

## 3 Numerical implementation

For concreteness, we will discuss the numerical implementation in the context of Fourier spectral methods, with the computational domain given by . We will discretize the domain using a regular Cartesian mesh with points discretizing each dimension: and . Differential operators are discretized in the usual way.

{remark}

Few of the details depend upon the choice to use Fourier spectral methods, and they are chosen because of their simple implementation, computational speed, and high order of accuracy. Inversion of the elliptic operator , as well as the extension operator is also greatly simplified when using Fourier spectral methods. However, any discretization based on a regular Cartesian mesh could be used with minimal modifications to the method.

In order to fully describe a discretization and solution strategy to Equation 18, we must describe several key elements.

1. In Section 3.1, we discretize the spread (, ) and interpolation (, ) operators.

2. In Section 3.2 we define a new regularization of the -function that is accurate and smooth.

3. In Section 3.3, we define our extension operator and show how this choice of controls the numerical conditioning of the extension problem.

4. In Section 3.4, we describe an efficient inversion strategy for the IBSE- system given by Equation 18.

5. Finally, in Section 3.5, we provide implementational details and briefly discuss the computational complexity of our inversion scheme.

### 3.1 Discretization of S, S∗, Tk, and T∗k

Let the boundary be parametrized by the function . In all examples in this manuscript, the boundary is one-dimensional and closed; the single parameter is defined on the periodic interval . The spread operator is defined as

 (SF)(x)=∫ΓF(s)δ(x−X(s))ds. (19)

The discrete version of this operator requires a regularized -function and a discretization of the integral over . Construction of a regularized -function with the properties necessary for the IBSE- method is non-trivial; we will delay discussion of this choice until Section 3.2 and simply denote the regularized -function as . Multivariate -functions are computed as Cartesian products of the univariate and also denoted by . Discretization of the integral over is made by choosing quadrature nodes , equally spaced in the parameter interval so that and . Quadrature weights are computed at each quadrature node to be ; this is a spectrally accurate quadrature rule for the integral of smooth periodic functions on . The discrete spread operator maps functions sampled at points in to by

 (SF)(x)=nbdy∑i=1F(si)~δ(x−Xi)Δsi. (20)

We do not adopt explicit notation to distinguish between the analytic and discretized operators. The number of points in the quadrature is chosen so that . This choice of node-spacing is wider than that recommended for the traditional IB method [1] but has been observed empirically to be the optimal choice in other studies of direct-forcing IB methods [4]. The interpolation operator may be defined by the adjoint property , but we note that the discrete interpolation operator produces a discrete function

 (S∗u)(sk)=∫Cu(x)δ(x−X(sk))dx. (21)

Discrete integrals over are straightforward sums computed over the underlying uniform Cartesian mesh.

Normal derivatives of are computed by the formula [30]

 ∂j~δ∂nj=ni1⋯ni2∂j~δ∂xi1⋯∂xij, (22)

where the Einstein summation convention has been used to indicate sums over repeated indices and is computed as Cartesian products of the appropriate derivatives of the one-dimensional . For example, in two dimensions is computed as

 ∂~δ∂n=nx~δ′⊗~δ+ny~δ⊗~δ′. (23)

Analogous to the definition of , we define the spread operator for the normal derivative, , as

 (T(j)F)(x)=(−1)jnbdy∑i=1F(si)∂j~δ∂nj(x−Xi)Δsi (24)

and again define the interpolation operator by the adjoint property . Analogous to the action of , the operator produces approximations of the normal derivative of a function defined on at the quadrature nodes of . Finally, the composite operator is defined by:

 (25)

while .

### 3.2 Construction of a smooth discretization of δ

The IBSE- method is capable of producing solutions that converge at . In order to achieve this accuracy, our choice of regularized -function must satisfy several conditions. For a general value of :

1. Its interpolation accuracy must be at least .

2. It must be at least , so that its normal derivative is continuous, allowing it to be used in the discretization of and .

In addition, we will require that the -function has compact (and small) support so that the spread and interpolation operators , , , and may be rapidly applied. The commonly used ‘four point’ -function is , has a support of four gridpoints, and produces approximations of , , and that are accurate to [31]. The use of this -function regularization is sufficient for the implementation of IBSE-, but it is not sufficient for higher order methods, due to both its limited accuracy and regularity.

In this manuscript, we discretize the IBSE- method for , , and , corresponding to second, third, and fourth order accuracy in for Dirichlet problems. We therefore need a regularization of that is and has an interpolation accuracy of for smooth functions. We are not aware of any functions with these properties currently defined in the literature, and so we construct such a function here. For simplicity, we do not attempt to impose other conditions that are often imposed on -functions, such as the even-odd condition or sum of squares condition [32]. The strategy that we follow is to choose a -function with sufficient accuracy but limited regularity and convolve it against itself to increase its smoothness. A similar approach was used in [33] to generate a smoother version of the traditional Peskin four point -function [1]. In order to provide an analytic formula, a base -function with a simple functional form must be used. We start with the function

 δIB4(r)=⎧⎪ ⎪⎨⎪ ⎪⎩1−12|r|−|r|2+12|r|30≤|r|≤1,1−116|r|+|r|2−16|r|31≤|r|≤2,02≤|r|, (26)

which has interpolation accuracy, regularity, and a support of four gridpoints [34, 35]. Define:

 ~δ=δIB4∗δIB4∗δIB4∗δIB4, (27)

where denotes convolution. It is clear that , its support is gridpoints, and it is a simple exercise to show that convolution preserves interpolation accuracy. The formula for is given in Appendix A.

{remark}

While our choice of -function regularization limits the spatial discretization in this paper to fourth-order accuracy, this is not a fundamental limitation of the method. One way to generalize to higher orders would be to construct smoother and more accurate regularizations of the -function similar to the one that we construct, which maintain finite support for computational efficiency. An alternative approach would be to use globally supported regularizations of the -function that are (i.e. sinc or Gaussian functions) and make use of fast sinc transforms or fast Gaussian transforms in order to rapidly apply the spread and interpolation operators [36, 37].

{remark}

Our construction of is likely not optimal, in that there probably exist -function regularizations accurate to with a support of less than 16 gridpoints (e.g. a -function with regularity and accuracy of with support of only six gridpoints is defined in [38]). The construction of such a -function would be worthwhile, allowing for coarser discretizations of problems with closely spaced boundaries and faster application of the operators , , , and .

### 3.3 Choice of extension operator Hk

Due to the nullspace of the periodic polyharmonic operator , we choose the operator used in the extension problem given by Equation 8 to be the invertible operator

 Hk=Δk+1+(−1)k+1Θ(k,n). (28)

Here is a positive scalar function that depends on the smoothness of the extension () and the number of Fourier modes () used in the discretization of the problem. The function is chosen to mitigate the numerical condition number of the operator :

 κ=1+(n/2)2(k+1)Θ. (29)

Minimizing the condition number (by taking to be large) must be balanced against the need to resolve the intrinsic length scale introduced to the problem. In practice, we find that taking to be

 Θ∗=max(1,αε(n2)2(k+1)), (30)

where is machine precision ( for double precision and for quadruple precision) and produces stable and accurate solutions. Numerical results are not particularly sensitive to the choice of (see Figure 3). This stability is demonstrated later (in Section 4.1) with highly accurate solutions shown for the IBSE-, IBSE-, and IBSE- methods up to . All numerical results in this paper are produced using to construct .

For large values of and , the value of can substantially impact the overall accuracy of the IBSE- algorithm. We demonstrate this effect in Figure 3 for , , and solutions to Equation 3. In Figure 3, we plot the error in the solution as a function of across 32 orders of magnitude. The minimum error observed is , at , which is within an order of magnitude of (denoted by a vertical grey line). Although there is some sensitivity to near the minimum, it is small: the error is bounded by over 6 orders of magnitude of . In Figure 3, we show a zoom of the extension near the domain boundary at for the solution produced using , which has an error of . The extension decays rapidly to , but is well resolved by the discretization: there are gridpoints between the boundary at and the peak at . In Figure 3, we show a zoom of the extension produced when , which gives a solution with an error of , four orders of magnitude worse than the error produced when . Although taking this large leads to a very well-conditioned operator , the extension decays very rapidly to , and the discretization is no longer able to fully resolve the length scales: there are only three gridpoints between the boundary at and the peak at . Finally, in Figure 3, we use . Although there are no fine length-scales to be resolved, ill-conditioning in the operator leads to an error of .

### 3.4 Inversion of the IBSE-k system (18)

In this section, we consider solving Equation 18 for an invertible elliptic operator in place of . The case where is the periodic Laplacian is complicated by the nullspace of : even though Equation 18 is invertible, the method for inversion that we outline in this section is not directly applicable; details for the resolution of this problem are provided in Appendix B. The system of equations is

 ⎛⎜ ⎜ ⎜ ⎜ ⎜⎝L−χELSHkTk−T∗kT∗kS∗⎞⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜⎝uξFG⎞⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝χΩf00g⎞⎟ ⎟ ⎟⎠. (31)

This system is large: in two spatial dimensions, the number of equations in the system is . By block Gaussian elimination, we can find a Schur-complement for this matrix, which we label :

 SC=(−T∗kT∗kS∗)(L−χELHk)−1(STk), (32)

along with an associated system of equations for the Lagrange multipliers and :

 (SC)(FG)=(S∗L−1χΩf−gT∗kL−1χΩf). (33)

The size of is comparatively small, only square. This equation is the key to the efficiency of the algorithm, as it allows the Lagrange multipliers and to be computed rapidly without first solving for or . Once and are determined, we may solve for :

 ξ=−(Hk)−1TkF, (34)

and once is known, may be computed as

 u=L−1(χDf+χELξ−SG). (35)

The computation of and requires only a series of fast operations: FFTs, multiplies, and applications of the spread operators.

Rapid inversion of the system of equations for and given in Equation 33 is not a trivial task. Because we have restricted to problems set on stationary domains and the size of is small, it is feasible to proceed using dense linear algebra. We form by repeatedly applying it to basis vectors. This operation is expensive: for two dimensions it is in the total number of unknowns . The Schur-complement is then factored by the LU algorithm provided by LAPACK [39]. Once this factorization is computed Equation 33 can be solved rapidly. This Schur-complement depends only on the domain and the discretization, so the LU-decomposition can be reused to solve multiple problems on the same domain or in each timestep when solving time-dependent problems. See Section 3.5 for a discussion of the computational complexity of the method and Table 2 for the actual numerical cost of the method applied to a test problem.

### 3.5 Summary of algorithm and numerical complexity

For a given geometry, choice of elliptic operator , and discretization size , the implementation of the IBSE- method proceeds as follows.

1. Setup

1. The boundary is discretized as described in Section 3.1, and stencils for the evaluation of and its first normal derivatives are computed at all nodes of the discrete boundary in to allow for the rapid application of the , , , and operators.

2. The Schur-complement matrix (Equation 32 for invertible , and Equation 70 for or other non-invertible ) is formed columnwise, by repeatedly applying it to basis vectors. For invertible , this task may be completed as follows (for non-invertible , the process is nearly identical):

1. One element of either or is set to , all other elements are set to .

2. These forces are spread to the Eulerian grid, by the application of and .

3. , and then , are computed by applying Equation 34 followed by Equation 35.

4. The quantities and are computed, and placed in the appropriate rows and column of the Schur-complement.

Note that item (ii) corresponds to application of the right-most matrix in Equation 32, item (iii) corresponds to inverting the middle matrix in Equation 32, while item (iv) corresponds to the application of the left-most matrix in Equation 32. We remark that although this operation is expensive (see Table 2), it is trivially parallelizable, as the task of applying the Schur-complement to each individual basis vector is independent.

3. The LU decomposition of the Schur-complement is computed.

4. Optionally, these items may be saved, bypassing the expensive computation and factorization of the Schur-complement in future uses of the IBSE- method when applied to the same geometry, , and discretization size .

2. Solve

1. The right-hand side of Equation 33 (or equivalently, Equation 70 when is not invertible) is computed.

2. Using the pre-computed LU decomposition of the Schur complement matrix, the Lagrange multipliers and are found by solving Equation 33 (or equivalently, Equation 70 when is not invertible).

3. The smooth extension , and finally , are computed via Equations 35 and 34 (or equivalently, Equation 72 when is not invertible)

In items 1(b), 2(a), and 2(c), application of the FFT and IFFT is used to move variables between frequency and physical space, while differential operators are applied in frequency space and multiplications with characteristic functions are carried out in physical space.

The computational work for this method can be broken into two main tasks: (i) the setup items, listed above, which are dominated by the formation and factorization of the Schur-complement defined in Equation 32 and (ii) the production of one solution given the Schur-complement, its factorization, the body forcing , and the boundary conditions . We summarize the scaling of the algorithm in Table 1; see Table 2 in Section 5.1 for the actual numerical cost of the method applied to a test problem. In two dimensions, the cost of a solve scales like the cost of an FFT; in three dimensions the cost is slightly worse: in the total number of unknowns . The higher cost of the algorithm in three-dimensions is due to the cost of factoring and inverting the Schur-complement, however, these matrices are highly structured. It is quite possible that the cost of inversion could be reduced with an appropriately preconditioned iterative method, by directly exploiting the structure of the matrix, e.g. with an inverse Multipole Method [40], or by indirectly exploiting that structure using an algorithm like HODLR [41] or Algebraic Multigrid [42].

## 4 Results: Poisson equation

### 4.1 One-dimensional test problem

To demonstrate both the power and the limitations of the IBSE method, we compute the solution to the one-dimensional example defined in Equation 3. Solutions are computed on grids ranging from to points, in both double- and quadruple-precision arithmetic, for the IB, IBSE-, IBSE-, and IBSE- methods. Figure 4 shows a refinement study demonstrating the expected first-, second-, third-, and fourth-order accuracy in for the IB, IBSE-, IBSE-, and IBSE- methods, respectively. The IB method achieves an error of with grid points. Imposing additional smoothness on the solution allows this error to be matched by the IBSE- method at , by the IBSE- method at , and by the IBSE- method at . For the IBSE- method, this is a factor of nearly 10000 less gridpoints; equivalent to a factor of 100 million less gridpoints for two-dimensional problems. In practice, obtaining solutions accurate to six digits with the traditional IB method is often impractical; with the IBSE method this kind of accuracy can be achieved on reasonably sized grids.

In double precision, all of the methods achieve the expected convergence rate up to some value of . The double-precision solutions begin to diverge from the quadruple-precision solutions at for IBSE-, at for IBSE-, and at for IBSE-. Using (see Section 3.3) to control the condition number of the extension operator provides sufficient numerical stability at all values of to enable the computation of solutions accurate to 10-11 digits. In quadruple precision, all methods exhibit the expected order of convergence across a wide range of values of , although the fourth order method begins to fall short of the expected path of convergence for extremely fine grids ().

{remark}

We note that this test problem, as well as the test problems in Sections 6.2, 6.1 and 5.1, are exterior problems set on periodic domains. Periodicity of solutions in the physical domain is enforced naturally in these problems through the use of the Fourier basis.

### 4.2 Two-dimensional test problem: solution inside a circle

Let be the circle of radius centered at , and identify with the square . We will solve the Dirichlet problem: \cref@addtoresetequationparentequation

 Δu =−4 in Ω, (36a) u =0 on Γ=∂Ω, (36b)

which has an analytic solution given by

 ua=4−(x−π)2−(y−π)2. (37)

We solve this problem in double-precision arithmetic with the IBSE-, IBSE-, and IBSE- methods for to . In Figure 5, we show error as a function of , indicating second-, third-, and fourth-order accuracy for IBSE-, IBSE-, and IBSE- respectively.

Using the solutions to this problem, we demonstrate the property that the IBSE method produces globally smooth solutions. In Figure 6, we show the solution to Equation 36, along with its first, second, and third derivatives in the -direction, produced using the IB, IBSE-, IBSE-, and IBSE- methods with . For simplicity, all functions are shown only along the slice , and the intersection of this slice with the boundary is shown as gray vertical lines. Derivatives are computed spectrally using the FFT. We expect solutions produced by the IBSE- method to exhibit discontinuities in the normal derivative across interfaces, and for all derivatives of lower order to be continuous. We see this expectation validated in the solutions presented in Figure 6: functions shown in plots with shading behind them in Figure 5(a) are at least continuous, all others are discontinuous. This property of global smoothness of the solutions enables the Fourier-spectral discretization to obtain high-order accuracy.

### 4.3 Two-dimensional test problem with Neumann boundary conditions

Traditional Immersed Boundary approaches are unable to provide solutions to Neumann problems since convolution-style estimation of the normal derivative is not convergent at the boundary due to the lack of regularity of the solution. The IBSE method does not have this limitation because of the additional smoothness of the solution. The only modification necessary to adapt Equation 18 to solve Neumann problems is to change the and operators that enforce and specify the boundary conditions. Consider the general Neumann problem: \cref@addtoresetequationparentequation

 Δu =f in Ω, (38a) ∂u∂n =g on Γ. (38b)

 T(1)F(x) =∫ΓFj(s)∂δ(x−s)∂nds, (39a) T∗(1)ξ(s) =−∫ξ(x)∂δ(x−s)∂ndx,∀s∈Γ. (39b)

The analogue to Equation 18 is

 ⎛⎜ ⎜ ⎜ ⎜ ⎜⎝Δ−χEΔT(1)Hk−T∗kT∗kT∗(1)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜⎝uξFG⎞⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝χΩf00g⎞⎟ ⎟ ⎟⎠. (40)
{remark}

It is just as simple to modify the IBSE method to allow the solution of Robin problems. Rather than replacing the upper-right and lower-left operators in Equation 18 with and , they should instead be replaced with linear combinations of , , , and .

Let . We will solve the problem \cref@addtoresetequationparentequation

 Δu =esinx(cos2x−sinx)−cosy in Ω, (41a) ∂u∂n =cos2xesinx−sin2y on Γ, (41b)

which has the analytic solution

 ua=esinx+cosy. (42)

We solve this problem with the IBSE-, IBSE-, and IBSE- methods, for which we expect first-, second-, and third-order convergence, respectively. The reason for the lower order of convergence than in Dirichlet problems is that our discretization of is one order less accurate than when acting on functions of the same regularity. Despite this, the solution produced by IBSE- still has global regularity. Solutions are computed in double precision, for to . Figure 7 shows error as a function of . Second- and third-order convergence is observed for the IBSE- and IBSE- methods. The IBSE- method converges at a rate that is slightly less than first-order.

## 5 Results: heat equation

We now consider solving the heat equation with Dirichlet boundary conditions: \cref@addtoresetequationparentequation

 ut−Δu =f(x,t) in Ω, (43a) u =g(x,t) on Γ. (43b)

Since the IBSE method can directly invert elliptic equations, it can easily be adapted to provide high-order implicit-time discretizations for the heat equation. Many implicit methods could be used; we choose BDF4 [43]:

 (I−1225ΔtΔ)un+1=125(12Δtf(t+Δt)+48un−36un−1+16un−2−3un−3). (44)

Applying this scheme to Equation 43b is equivalent to solving the problem: \cref@addtoresetequationparentequation

 Lun+1 =~f in Ω, (45a) un+1 =g(x,t+Δt) on Γ, (45b)

 L =I−1225ΔtΔ, (46a) ~f =125(12Δtf(t+Δt)+48un−36un−1+16un−2−3un−3). (46b)

This may be solved using the algorithm described in Section 3.4. There are several comments worth making regarding this discretization:

1. The use of an implicit scheme allows for large timesteps to be taken. For all test problems presented in this section, we choose to be:

 Δt=tfinal2⌈tfinalΔx⌉−1, (47)

where is the ceiling function, so that .

2. The fact that depends on implies that the Schur-complement depends on as well. This means that the Schur-complement must be reformed and refactored whenever is changed, complicating the use of a method that uses adaptive timestepping.

3. Modifying this formulation to solve the heat equation with Neumann and Robin boundary conditions may be done in the same way as shown for the Poisson equation in Section 4.3.

### 5.1 Solution in a periodic domain outside an obstacle

In this section we demonstrate high-order convergence to a heat-equation set on a periodic domain with a simple obstacle in it. To provide direct numerical comparison with results from the Active Penalty (AP) method, the following problem is from Section 6.2 of [26]: \cref@addtoresetequationparentequation

 ut−Δu =[cosy+esinx(sinx−cos2x)]cost−(esinx+cosy)sint in Ω, (48a) u(x,y,t=0) =esinx+cosy in Ω, (48b) u =(esinx+cosy)cost on Γ, (48c)

where the physical domain is . Equation 48 has the analytic solution

 ua=(esinx+cosy)cost. (49)

The initial condition is integrated from to by solving Equation 45 at each timestep. Startup values for BDF4 are computed by evaluating the analytic solution at , , and . In Figure 8, errors at for solutions generated with the IBSE-, IBSE-, and IBSE- methods are shown, demonstrating second-, third-, and fourth-order convergence in space and time, respectively. On this plot, we also show errors from the second- and third-order AP methods with , the finest resolution reported in [26]. Our method yields errors lower by several orders of magnitude.

In contrast to the AP method [26], we are able to use an implicit time discretization, enabling the use of large timesteps. The ability to take large timesteps efficiently comes at the cost of a significant precomputation for the IBSE method, a cost that the AP method does not incur. This precomputation prevents the IBSE method as described in this paper from being efficiently applied to moving boundary problems, a limitation not faced by the AP method.

In Table 2, we provide the computational time required by the IBSE-, IBSE-, and IBSE- methods to setup the problem (dominated by the formation and factorization of Equation 32), the time required to take one timestep, and the number of timesteps taken to advance the equation from to . Computational times are given as a multiple of the time required to take one FFT. In our code, 8 FFTs are used per timestep regardless of ; the time required to take a timestep does not vary significantly across the methods. The precomputation time increases approximately linearly in , but even for the finest discretization presented, the precomputation time is not prohibitive: for and the IBSE- method, the precomputation requires about 50 minutes using serial code that has not been carefully optimized.

### 5.2 Solution inside a parametrically defined region

In this section, we demonstrate high-order convergence to a heat equation set inside a parametrically defined domain. This problem is from Section 6.1 of [24], allowing direct comparison of the IBSE method with the Fourier Continuation method [24]. For convenience, we rescale the domain of the problem. We will solve \cref@addtoresetequationparentequation

 ut−Δu =π(2−Δϕ)cos(πϕ)+π2(ϕ2x+ϕ2y)sin(πϕ) in Ω, (50a) u(x,y,t=0) =sin(πϕ(x,y,t=0)) in Ω, (50b) u =sin(πϕ) on Γ, (50c)

where

 (51)

The domain is the region inside the boundary defined by the parametric equations: \cref@addtoresetequationparentequation

 x(θ) =(10sin2(2θ)+3cos3(2θ)+40)cosθ/20+π, (52a) y(θ) =(10sin2(2θ)+3cos3(2θ)+40)sinθ/20+π, (52b)

for . Equation 50 has the analytic solution . The initial condition is integrated from