Immersed Boundary Smooth Extension (IBSE): A highorder method for solving incompressible flows in arbitrary smooth domains
Abstract
The Immersed Boundary method is a simple, efficient, and robust numerical scheme for solving PDE in general domains, yet for fluid problems it only achieves firstorder spatial accuracy near embedded boundaries for the velocity field and fails to converge pointwise for elements of the stress tensor. In a previous work we introduced the Immersed Boundary Smooth Extension (IBSE) method, a variation of the IB method that achieves highorder accuracy for elliptic PDE by smoothly extending the unknown solution of the PDE from a given smooth domain to a larger computational domain, enabling the use of simple Cartesiangrid discretizations. In this work, we extend the IBSE method to allow for the imposition of a divergence constraint, and demonstrate highorder convergence for the Stokes and incompressible NavierStokes equations: up to thirdorder pointwise convergence for the velocity field, and secondorder pointwise convergence for all elements of the stress tensor. The method is flexible to the underlying discretization: we demonstrate solutions produced using both a Fourier spectral discretization and a standard secondorder finitedifference discretization.
keywords:
Embedded boundary, Immersed Boundary, Incompressible NavierStokes, Fourier spectral method, Complex geometry, HighorderremarkRemark
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 fluidstructure 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 NavierStokes 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 by convolutions with regularized functions. These methods have many desirable properties: they make use of robust and efficient Cartesiangrid 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 structuredgrid 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 highorder 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 onedimensional 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 . The lack of regularity in the analytic problem leads to loworder convergence in many numerical schemes, including the Immersed Boundary method: the addition of (regularized) singular forces supported at the boundary causes the solution to be , and solutions are accurate only to firstorder in the grid spacing [4, 5].
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 [27], we introduce another approach: the Immersed Boundary Smooth Extension (IBSE) method for the solution of the Poisson and related problems (e.g. the heat equation and Burgers equation). The IBSE method builds off of the basic framework of the direct forcing IB method. Drawing on ideas from the AP and FC methods, the IBSE method remedies the slow convergence of the direct forcing IB method by solving an auxiliary problem to extend the unknown solution from the physical domain into the entire computational domain. This extension is used to define the body forcing in the nonphysical portion of the domain, forcing the solution to be globally . Highorder accuracy is then achieved naturally, up to or the maximum accuracy supported by the underlying discretization for a smooth problem.
In this paper, we extend the work from [27] to apply to PDE requiring the imposition of a divergence constraint, including the Stokes and NavierStokes equations. This is a nontrivial difficulty for a highorder embedded boundary scheme based on smooth extensions. Prior work has either been restricted to the compressible NavierStokes equations [23], or has dealt with the pressure in an adhoc manner, limiting the overall accuracy of the numerical scheme to secondorder [26]. These methods either extend the forcing function or solution from a prior timestep, and are thus unable to provide solutions to the Stokes problem. For timedependent problems, they require either explicit timestepping or the use of AlternatingDirectionImplicit (ADI) schemes to take implicit timesteps, complicating the ability to achieve highorder accuracy in time.
In contrast, we work directly with the Stokes equations, rather than employing dimensional splitting or decoupling the problem into a set of Poisson equations. The smooth extensions to the solution of the Stokes equation define forcing functions in the nonphysical domain in a coupled manner, and the Lagrange multipliers supported on the boundary that define these extensions are fully coupled through the Stokes problem, eliminating any need to impose artificial boundary conditions for a pressure Poisson equation. The IBSE method smoothly extends the unknown solution, allowing for efficient direct solutions of the steady incompressible Stokes equation and simple highorder in time implicitexplicit time discretization of the NavierStokes equations. Remarkably, the fully coupled problem for the solution to the incompressible Stokes equations (combined with the equations that define the smooth extensions to the velocity and pressure fields) can be reduced to the solution of a relatively small^{1}^{1}1With size a small multiple of the number of points used to discretize the boundary dense system of equations, two Stokes solves on a simple domain ignoring the internal boundary, and a handful of fast Fourier transforms. The dense system depends only on the boundary and the discretization, and so it can be formed and prefactored to allow for efficient timestepping.
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. As with the IB method, the IBSE method is flexible to the choice of the underlying PDE discretization, so long as the mesh is uniform in the vicinity of the immersed boundary, facilitating simple coupling to existing code. In this paper, we demonstrate the IBSE solver with an underlying Fourier spectral discretization as well as a secondorder staggeredgrid finite difference discretization.
This paper is organized as follows. In Section 2, we introduce the methodology, ignoring the details of the numerical implementation. This section includes a review of the relevant methodology introduced in [27]. In Section 3, we discuss the numerical implementation of the method, not including the choice of the underlying discretization of the PDE. In Section 4, we study a steady Stokes problem with an analytic solution, using a Fourier spectral discretization. We analyze the convergence of the velocity and pressure fields, as well as derivatives of the velocity fields. We demonstrate up to thirdorder pointwise convergence of the velocity fields, as well as secondorder pointwise convergence of all elements of the stress tensor. In Section 5, we analyze the flow of a viscous incompressible fluid around a cylinder in a confined channel at zero Reynolds number. We solve the problem using both a Fourier spectral discretization and a secondorder finite difference discretization, comparing the accuracy and convergence of the solution, as well as the convergence of the scalar drag coefficient to known values from the literature. Finally, in Section 6, we solve the incompressible NavierStokes equations in both steady and unsteady cases, for the flow around a cylinder in a confined channel, demonstrating rapid convergence and quantitative agreement with a number of benchmarks from the literature.
2 Methods
We begin by considering the time dependent, incompressible NavierStokes equations: \cref@addtoresetequationparentequation
(1a)  
(1b)  
(1c) 
where the boundary of the domain is denoted by . One way to discretize these equations in time is to treat the nonlinearity explicitly, and to treat the Stokes operator implicitly. The simplest such time discretization, using forward Euler for the nonlinear terms and backward Euler for the remaining terms, is \cref@addtoresetequationparentequation
(2a)  
(2b)  
(2c) 
Although the NavierStokes equations are often discretized in time using projection methods [28], these discretizations produce large splitting errors near the boundary at low Reynolds number. ImplicitExplicit timediscretizations of the form of Equation 2 eliminate these errors, and are simpler to generalize to higherorder in time schemes (see Section 6). The unsplit Stokes system may be inverted efficiently using GMRES with a projection preconditioner [29]. We find this inversion strategy to be very efficient (convergence to an algebraic tolerance of achieved in 34 iterations) for moderate to high Reynolds numbers (), and moderately efficient (convergence achieved in approximately 35 iterations) for the Stokes problem where projection methods cannot be used.
Motivated by the form of Equation 2, in this section we focus primarily on developing the methodology required to solve the general Stokes problem: \cref@addtoresetequationparentequation
(3a)  
(3b)  
(3c) 
When and , Equation 3 reduces to the steady incompressible Stokes equations. The case where arises in the study of porous media [30] and from discretizing the NavierStokes equations in time using an implicitexplicit discretization, as in Equation 2. The nonzero divergence constraint () adds no additional complexity and is included for completeness. Before discussing the IBSE method as applied to Equation 3, we first review relevant material from [27], in order to introduce the central ideas in the simplified context of the Poisson problem.
2.1 Review of [27]: the IBSE method for the Poisson problem
Suppose that we wish to solve the Poisson problem \cref@addtoresetequationparentequation
(4a)  
(4b) 
in an arbitrary smooth domain . The IBSE method works by smoothly extending the unknown solution of a PDE from to a simplier domain . The domain will be referred to as the physical domain and will be referred to as the computational domain. It is assumed that . The interior boundary will be assumed to be smooth, not selfintersecting, and must separate into the two disjoint regions and , referred to as the extension domain. Neither nor need be connected. An additional simple computational domain is required to solve the auxiliary equation that defines the extension of the unknown solution. Although may coincide with , it need only contain . The boundary of the computational domain will be denoted by . Two typical domains are shown in Figure 1. In Figure 1, is periodic, and no boundary conditions need to be specified on . In Figure 1, Dirichlet conditions are imposed on . To prevent confusion with the interior boundary , the boundary and its respective boundary condition will only be explicitly included in equations when required to prevent confusion, e.g. in Sections 5.2 and 3.4.
We will denote the extension of the unknown solution by . This extension of the solution is then used to define a volumetric forcing in the region . With this forcing, an extended problem in all of the simple domain may be solved: \cref@addtoresetequationparentequation
(5a)  
(5b) 
The solution gives the desired solution in and is equal to in . Because was chosen to be a smooth extension to , the function is globally smooth in .
The extension to the unknown function is defined as a solution to a highorder PDE which takes for its boundary conditions matching criteria of the form . This allows the extension to be defined by a small number of unknowns (proportional to the number of points used to discretize the boundary). The extension PDE for is solved efficiently in the simple domain using an Immersed Boundary type method.
In order to succinctly describe the methodology we will require some notation. We define the spread operator:
(6) 
and the interpolation operator:
(7) 
for the normal derivative, where is a parametrization for with in the parameter interval . The integral on the right hand side of Equation 7 is notation for the action of the distribution on the smooth function . The traditional Immersed Boundary spread and interpolation operators, and are denoted by and respectively. The interpolation operator maps function values (or the normal derivatives of the function) from to , while the spread operator maps singular (and hypersingular) forces supported on to . To simplify equations presented in the forthcoming material, we also define the operators , , and by: \cref@addtoresetequationparentequation
(8a)  
(8b)  
(8c) 
The operator provides an interpolation of a function and its first normal derivatives to the boundary; provides an interpolation of the first normal derivatives to the boundary, but excludes the values; the spread operator represents a set of singular forces (like) and hypersingular forces (like the first normal derivatives of the function) on the boundary.
The central challenge of the IBSE method is to compute the smooth extension to an unknown solution. We first discuss how to compute an extension to a given function. Let be given. To compute a extension to , we solve the following highorder PDE in the region : \cref@addtoresetequationparentequation
(9a)  
(9b) 
Here is an appropriate differential operator such as the polyharmonic operator ; we discuss the choice of this operator in more detail, from a numerical perspective, in Section 3.2. This problem may be solved on the simpler domain using methodology directly analogous to the direct forcing Immersed Boundary method. The boundary conditions given in Equation 9b that force to share its first normal derivatives with along are enforced by the addition of unknown singular and hypersingular forces supported on the boundary: \cref@addtoresetequationparentequation
(10a)  
(10b) 
Future equations similar to Equation 10b will be shortened simply to , and assumed to hold for all . Notice that is not actually an extension to : that is, in . We will only be interested in the function in , and so need not form its literal extension (which is .
To solve the Poisson problem given by Equation 4 using the IBSE method, we instead solve the extended problem given in Equation 5. The forcing function that is specified in must be chosen so that it forces the extended solution to be . Let smoothly extend , that is, we ask that is globally smooth in and that it satisfies the constraints
(11) 
at the interface . These constraints require that the first normal derivatives of agree with the first normal derivatives of on the boundary. Notice that it is not enforced that the values of agree with the values of on ; we will see momentarily that this is unnecessary. The forcing function is then defined as:
(12) 
Coupling these equations together, we obtain the IBSE formulation for the Poisson problem given by Equation 4: \cref@addtoresetequationparentequation
(13a)  
(13b)  
(13c)  
(13d) 
These equations are (13a) the extended Poisson problem, (13b) the extension PDE, (13c) the interface constraints that force and to share their first normal derivatives, and (13d) the physical boundary condition to the Poisson problem. Both the interface constraints and the physical boundary condition are enforced by the unknown forces in the extension PDE. We will subsequently drop the notation and refer to the solution of the IBSE system of equations simply as .
In , satisfies the physical Poisson equation; while in , satisfies: \cref@addtoresetequationparentequation
(14a)  
(14b) 
and thus in up to a constant. Since implies that , then for , the forcing , and so , and is thus at least (for spatial dimensions )^{2}^{2}2The notation denotes the space of measurable functions with the property that the function itself, as well as its first two weak derivatives, are square integrable. The Sobolev embedding theorem implies that any function in is equal almost everywhere to a continuous function, at least in spatial dimensions and [31].. This continuity, along with the smoothness constraints guaranteed by the constraint that imply that is globally .
In [27], we verify that the IBSE formulation of the problem produces solutions that converge at a rate of in the gridspacing , for the Poisson problem, as well as the heat equation, Burgers equation, and the FitzhughNagumo equations.
In [27], both a volumetric force and a singular force distribution are added to the Poisson equation, giving a slightly different formulation for the IBSE method:
(15)  
(16)  
(17)  
(18) 
The singular force distribution acts to enforce that the boundary values of equal those of , this is an unnecessary requirement for the IBSE method. We will proceed with the simpler formulation given in Equation 13 when developing the methodology for the Stokes equation.
2.2 Solution of the general Stokes problem
We now focus on the modifications to the IBSE method required to solve the general Stokes problem: \cref@addtoresetequationparentequation
(19a)  
(19b)  
(19c) 
We assume that and are sufficiently smooth functions defined on (for ) and (for ). Here denotes the Helmholtz operator and denotes a boundary operator. We discuss the boundary operator associated with the Dirichlet problem (); imposition of other boundary conditions is similar and discussed in [27]. Note that this equation reduces to the incompressible Stokes problem when and are zero.
In the direct forcing Immersed Boundary (IB) method, this problem is solved in a simple computational domain by adding singular forces supported on the boundary that act as Lagrange multipliers which enforce the boundary condition. For Dirichlet problems this can be represented as \cref@addtoresetequationparentequation
(20a)  
(20b)  
(20c) 
The singular forces supported on induce jumps in the normal derivatives of at the boundary; the velocities produced by the IB method are generically while the pressure field typically has jump discontinuities. For Eulerian discretizations of the underlying PDE that are ignorant of the boundary, this leads to the slow convergence rate for the velocity field in . The pressure field fails to converge near to the boundary, but does converge at the rate in .
In order to improve this convergence, the Immersed Boundary Smooth Extension (IBSE) method extends the unknown solutions and to be and , respectively. Suppose that the solution vector to Equation 19 is known in . Define by \cref@addtoresetequationparentequation
(21a)  
(21b) 
and by \cref@addtoresetequationparentequation
(22a)  
(22b) 
Recall that is a highorder differential operator such as , which will be defined precisely in Section 3.2. We make two remarks concerning the extension of the pressure. {remark} The pressure function is typically thought of as a Lagrange multiplier: it equals whatever is required to enforce the divergence constraint that . One may suspect therefore that it is unnecessary to construct an explicit extension of the pressure function as we do in Equation 22. Unfortunately, this is not true. Consider the problem: \cref@addtoresetequationparentequation
(23a)  
(23b)  
(23c) 
along with the modified problem where is replaced by for any scalar function . The pressure will be changed: . But the solution will be left unchanged, and thus its extension functions will be unchanged. If the extension of the pressure function depends only locally on , it cannot extend both and smoothly. {remark} The choice to extend the pressure function to be rather than yields improved numerical stability. When , the structure of the Stokes equations implies that . Applying the operator to the function with that global regularity yields an inconsistent estimate of , i.e. , and thus enforcing that does not actually enforce that .
From these extensions, we may define a set of forces : \cref@addtoresetequationparentequation
(24a)  
(24b) 
Analogous to Equation 5, to smoothly extend the Stokes problem given by Equation 19 from to , we add to Equation 19a and to Equation 19b in the extension region to obtain \cref@addtoresetequationparentequation
(25a)  
(25b)  
(25c) 
In the physical domain , we have that \cref@addtoresetequationparentequation
(26a)  
(26b)  
(26c) 
and thus solves the original Stokes problem given in Equation 19. In the extension domain , we have that \cref@addtoresetequationparentequation
(27a)  
(27b)  
(27c) 
Thus and , both up to a constant (in ). Since is enforced to match at the boundary by Equation 22, they are equal at the boundary, and thus in , implying that is globally smooth in , while the global regularity for is implied by the same argument given in Section 2.2. Numerical approximations of these solutions may converge rapidly, even when the PDEs are discretized using an underlying Cartesian grid discretization that is ignorant of the boundary. The difficulty in this formulation comes from the fact that the extensions and depend on the unknown solution .
Combining Equations 19, 22 and 21, we obtain a coupled system for the solution together with the smooth extensions to that solution:
\cref@addtoreset
equationparentequation
(28a)
(28b)
(28c)
(28d)
(28e)
(28f)
(28g)
We will refer to this system of equations as the IBSE equations. The remainder of this paper is dedicated to demonstrating that these equations provide accurate solutions to the Stokes and NavierStokes equations, that they allow a flexible choice of discretization, and that they admit an efficient inversion strategy. We organize the presentation of information as follows:

In Section 3, we discuss the details of numerical implementation that do not depend upon the underlying choice of discretization.

In Section 4, we solve a Stokes problem with an analytic solution that allows us to carefully verify the convergence rates of the solution (up to thirdorder for and secondorder for , in ), as well as the convergence of elements of the fluid stress tensor (up to secondorder in ).

In Section 5, we describe both a Fourier spectral and a finitedifference discretization to the confined channel flow around a cylinder problem for incompressible Stokes flow. We validate the convergence rates for solutions produced by the IBSE method and compare to known benchmarks.

In Section 6, we discuss timestepping in the IBSE framework, solve both steady and unsteady NavierStokes problems, and demonstrate rapid convergence of the solutions and agreement with known benchmarks.
3 Numerical implementation
In this section, we discretize the coupled IBSE system given in Equation 28. We assume that the functions are discretized using an underlying Cartesian mesh with uniform^{3}^{3}3For the discretization described here, the mesh needs to be uniform only in a small neighborhood of the boundary, with the same width as the regularized function used to discretize the operators and . In fact, the mesh need not even be uniform, but this would add some complexity to the discretization [32]. gridspacing but otherwise make no assumptions regarding the underlying discretization of the functions and differential operators. These final details will be treated in Section 4 for a Fourier spectral discretization, and in Section 5.2 for a standard secondorder finite difference scheme. Little varies in the details.

In Section 3.1, we discretize the spread and interpolation operators for the normal derivative that were introduced in Equations 7 and 6. Note that the discretization of these operators automatically induces a discretization for the composite operators , , and , defined in Equation 8.

In Section 3.2, we discuss how we choose the extension operator introduced in Equation 10.

In Section 3.3, we describe an efficient inversion strategy for the IBSE system given by Equation 28.

In Section 3.4, we provide a brief summary of the algorithm.
3.1 Discretization of singular integrals
Let the boundary be parametrized by the function . In all examples in this manuscript, we work with a twodimensional fluid, and the boundary is onedimensional and closed, with the single parameter defined on the periodic interval . The discrete version of the spread operator for the normal derivative, defined in Equation 6 requires a regularized function and a discretization of the integral over . The regularized function that we use is analytically three times differentiable, satisfies four discrete moment conditions, and has a support width of . The function and its properties are introduced in [27, 33]; this discrete function will be represented by . Multivariate functions are computed as Cartesian products of the univariate and also denoted by . Normal derivatives of are computed by the formula [34]
(29) 
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
(30) 
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
(31) 
We do not adopt explicit notation to distinguish between the analytic and discrete operators. The number of points in the quadrature is chosen so that . This choice of nodespacing 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 directforcing IB methods [4]. We define the interpolation operator by the adjoint property , but note that the discrete interpolation operator produces a discrete function
(32) 
Discrete integrals over are straightforward sums computed over the underlying uniform Cartesian mesh, and may be efficiently evaluated due to the finite support of . When acting on vector functions, all operators defined in this section are assumed to operate elementwise.
3.2 Solution of the extension equations and the choice of extension operator
Due to the difficulty in accurately and efficiently inverting the highorder differential operator needed to define the smooth extensions, we solve these equations utilizing a Fourier spectral discretization regardless of the discretization used for the Stokes equations. Let be the largest wavenumber associated with the discrete Fourier transform (DFT) on the discretized domain . We choose the extension operator introduced in Equation 9 to be:
(33) 
Here is a positive scalar function that depends on the smoothness of the extension () and the largest wavenumber (). The function is chosen to mitigate the numerical condition number of the operator :
(34) 
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. We set the length scale as a function of , in effect allowing the extension to decay to zero over some number of grid cells. We thus choose
(35) 
where is proportional to the number of gridpoints over which the extension decays. The choice of can impact the accuracy of the scheme: if is too large, the condition number of may be high with respect to the precision of the computer being used, and numerical instability causes a degradation of the quality of the solution. If is too small, short length scales that are not resolvable by the discretization are introduced to the problem. For all examples presented in this paper, is chosen to be . In Section 4, we investigate the effect of the choice of on a specific problem.
3.3 Inversion of the IBSE system, Equation 28
In this section, we describe an efficient inversion strategy for Equation 28. We will assume that the discrete Stokes operator is invertible^{4}^{4}4The nullspace in the pressure function should be fixed by the discrete operator; e.g. by enforcing that .. The details for when it is not invertible (e.g. when solving on a periodic domain with ) are analogous to the Poisson problem, and the solution is presented in [27]. In matrix form, Equation 28 is given by:
(36) 
This system is large, but the number of constraints (the last three equations) is relatively small. By block Gaussian elimination, we can find a Schurcomplement for this matrix, which we label :
(37) 
along with an associated system of equations for the Lagrange multipliers and ,
(38) 
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, and may be found by solving Equations 28d and 28c. With and known, and may be computed simply by solving the Stokes equation in , ignoring the boundary (that is, by solving Equations 28b and 28a).
Rapid inversion of the system of equations for and given in Equation 38 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 Schurcomplement is then factored by the LU algorithm provided by LAPACK [35]. Once this factorization is computed Equation 38 can be solved rapidly. This Schurcomplement depends only on the domain and the discretization, so the LUdecomposition can be reused to solve multiple problems on the same domain or in each timestep when solving timedependent problems.
3.4 Outline of the methodology
We provide a brief outline of the steps required for solving the general Stokes problem given by Equation 19 using the IBSE method. The algorithm proceeds as follows:

Form the right hand side of the Schurcomplement equation, as defined in Equation 38. This may be done by solving \cref@addtoresetequationparentequation
(39a) (39b) (39c) for and then computing , and . Note that the internal boundary is ignored when computing the solution to Equation 39.

Compute and by solving the Schurcomplement system given by Equation 38, given the right hand side computed above. We assume that the Schur complement has been formed and its LU decomposition found as a precomputation.

Compute and from and by solving Equations 28d and 28c. We take to always have periodic boundary conditions, and these equations may be simply inverted using the fast Fourier transform.

With and known, find and by solving Equations 28b and 28a. To be precise, we solve \cref@addtoresetequationparentequation
(40a) (40b) (40c) Note again that the internal boundary is ignored when computing this solution.
Step two of this algorithm requires that the Schurcomplement be formed and factored as a preprocessing step. The Schur complement operator maps a set of singular (and hypersingular) forces to the constraints given by equations Equations 28g, 28f and 28e. We form the matrix by applying to basis vectors. To apply to any set of singular forces, first compute and from and by solving Equations 28d and 28c. Then and may be found by solving \cref@addtoresetequationparentequation
(41a)  
(41b)  
(41c) 
and finally with , , and known, then , , and can be evaluated. Note that , and do not enter into this computation. The Schur complement depends only on the geometry and the discrete operators , , .
4 Results: Stokes equation with analytic solution
To demonstrate the improved convergence properties of the IBSE method as compared to the IB method, we construct a simple example with an analytic solution and compare the accuracy of the solutions with those produced by the Immersed Boundary method. For the numerical comparisons with the IB method provided in this section, we show two sets of values: those produced using a regularized function commonly used in IB simulations [10]:
(42) 
as well as those produced using the smoother and more accurate function introduced in [27] and used in all IBSE simulations in this paper. These two methods will be denoted by and , respectively. We solve the general Stokes problem with and Dirichlet boundary conditions: \cref@addtoresetequationparentequation
(43a)  
(43b)  
(43c) 
where the solution is given by \cref@addtoresetequationparentequation
(44a)  
(44b)  
(44c) 
the forcing function is defined as
(45) 
and the boundary condition is found by evaluating on . The physical domain is taken to be , where denotes the 2torus identified with the periodic rectangle , and where denotes the ball of radius centered at . The computational domain is taken to be . Plots of the analytic solution in the physical domain are shown in Figure 2.
For this problem, it is convenient to use a Fourier spectral discretization for and to choose . We discretize the domain using a regular Cartesian mesh with points discretizing each dimension: and for . Differential operators are discretized in the usual way.
A refinement study for the and errors in the solutions for and is shown in Figure 3. The refinement path for is omitted; it is similar to what is shown for . In , the velocity field converges at first, second, and thirdorder accuracy in for the IB, IBSE and IBSE methods, respectively, while the pressure field converges at first and secondorder accuracy for IBSE and IBSE, but at the slow rate of for the IB method. In , the rates of convergence for the velocity remain unchanged in all methods, while for the pressure the difference is starker: the IBSE and IBSE methods maintain first and secondorder convergence, while the IB method fails to converge pointwise.
It may be surprising that the velocity fields produced by the Immersed Boundary method converge at only firstorder in . Indeed, for a given and , the velocity fields produced by solving \cref@addtoresetequationparentequation
(46a)  
(46b)  
(46c) 
converge at , , and in the , , and norms, respectively [36]. However, is not given, it is computed by inverting a Schurcomplement equation whose elements include terms of the form where , and thus contain errors.
In the refinement study shown in Figure 3, the pressure function converges more slowly than the velocity function; in the case of the IB method, it fails to converge pointwise. This will be true in general for elements of the stress tensor, which are one derivative less smooth than the velocity field. In Figure 4, we show a refinement study for the and errors in the derivatives and (the convergence paths for and are omitted, but essentially the same). A similar pattern to the convergence of the pressure function is observed: for all of the derivatives, the IB method converges slowly (at ) in , but fails to converge pointwise; the IBSE and IBSE method provide first and secondorder convergence in both and . We emphasize that with pointwise convergence of the pressure function as well as pointwise convergence of these derivatives, the IBSE method is able to capture the stress tensor pointwise up to the boundary .
Finally, we use this test problem to numerically investigate the choice of the free parameter in the definition of the extension operator . This parameter controls the rate at which the extension functions and decay to away from the boundary. We compute solutions with , using extensions across a wide range of choices of (from to ). The error for , and the condition number of the Schur complement, as a function of , are shown in Figures 5 and 4(a). For very small values of , the error is dominated by the inability of the discretization to resolve the small length scale introduced to the problem. As is increased, the error produced by the IBSE method improves, but the condition number of the Schurcomplement worsens. As grows large, the error again begins to degrade, dominated by the error produced when inverting the illconditioned Schurcomplement. This presents a balance between resolution and conditioning. For this test, the minimum error is somewhere between and , a result that has been fairly robust across the range of test problems we have examined. Other than for this test, all numerical results in this paper are produced using with to construct . In Figures 5 and 5, we show the extension function for the solution to Equation 43 for and .