An inverse problem formulation of the immersed boundary method^{†}^{†}thanks: Submitted to the editors DATE. Portions of this article were published in the conference proceedings [42]. \fundingThis material is based upon work supported by the National Science Foundation under Grant No. 1825991.
Abstract
We formulate the immersedboundary method (IBM) as an inverse problem. A control variable is introduced on the boundary of a larger domain that encompasses the target domain. The optimal control is the one that minimizes the mismatch between the state and the desired boundary value along the immersed targetdomain boundary. We begin by investigating a naïve problem formulation that we show is illposed: in the case of the Laplace equation, we prove that the solution is unique but it fails to depend continuously on the data; for the linear advection equation, even solution uniqueness fails to hold. The illposedness is rectified by two complimentary strategies. The first strategy is to ensure that the enclosing domain tends to the true domain as the mesh is refined. The second strategy is to include a specialized parameterfree regularization that is based on penalizing the difference between the control and the state on the boundary. The proposed inverse IBM is applied to the diffusion, advection, and advectiondiffusion equations using a highorder discontinuous Galerkin discretization. The numerical experiments demonstrate that the regularized scheme achieves optimal rates of convergence and that the reduced Hessian of the optimization problem has a bounded condition number as the mesh is refined.
remarkRemark \newsiamremarkhypothesisHypothesis \newsiamthmclaimClaim \headersinverse immersed boundary methodJ. Yan and J.E. Hicken
mmersedboundary method, inverse problem, PDEconstrained optimization
65M60, 65N12
1 Introduction
Mesh generation remains a significant bottleneck for many aerospace applications of computational fluid dynamics (CFD) [38] and the numerical solution of partial differential equations (PDEs) more generally. This bottleneck is particularly acute during the generation and adaptation of curvedelement, anisotropic meshes around complex geometries, which has limited the adoption of highorder methods in industry.
One way to address the meshing bottleneck is to develop a highorder discretization that does not require a conforming mesh. The immersed boundary method (IBM) offers a potential framework for constructing such a discretization and forms the basis for the method presented herein.
Before proceeding, we should distinguish the IBM from immersedinterface, or cutcell, methods. IBMs [29] impose the boundary conditions indirectly through a body force or modified boundary flux. Immersedinterface methods, on the other hand, modify the cells, elements, or stencil near the boundary such that the boundary condition can be applied directly [32, 5, 19, 1, 14, 12, 20, 8, 39, 24]. Immersedinterface methods offer an alternate approach to addressing the meshing bottleneck, but they are not without their difficulties. The process of “cutting” cells and elements is a nontrivial task whose complexity is arguably on par with geometryconforming mesh generation. Furthermore, the cutting process invariably creates relatively small and nonstandard element shapes or stencils that must be treated carefully to avoid poor conditioning and/or nonlinear solver convergence issues [39]. These problems are compounded when considering high Reynolds number flows that necessitate highly stretched grids to accurately resolve the boundary layer.
The present work is more closely related to classical IBMs that introduce a body force, or penalty, that imposes the boundary conditions indirectly. The IBM framework encompasses a wide range of approaches, so a comprehensive review is beyond the scope of this work. For reviews of the IBM in the context of finitedifference/volume methods see [30] and Sections 3 and 4.1 of [28]. A review of penaltybased IBMs used in finiteelement methods can be found in Section 2.2 of [20].
IBMs are popular in some CFD applications, but they have not found widespread use for the steady, advectiondominated problems common in the aerospace industry. We believe this is primarily due to their limited accuracy; most IBMs are firstorder accurate in practice, because the underlying solution is nonsmooth at the interface that represents the boundary [30].
While uncommon, highorder IBMs have been proposed and we mention a few here. Mayo [25] solved the Poisson’s and biharmonic equations on irregular domains by constructing a discontinuous extension of the solution onto a regular (e.g. square) domain. The discontinuities in the solution and its derivatives could then be determined by solving an integral equation, and these jumps could subsequently be introduced into the Poisson’s or biharmonic equations to solve the PDEs on the regular domain using fast solvers. Mayo later extended this technique to 4th order in [26].
Marques et al. [23] and Marques [22] generalized the ghostfluid method [11, 10, 21, 17] to highorder by defining a correction function that smoothly extends the solution on either side of an interface or boundary. Once the correction function is determined, it is included on the righthand side of the discretized equations. The correction is independent of the solution for linear problems with Dirichlet boundary conditions, so conditioning of the lefthand side is not affected in this case; however, more generally, the correction function depends on the solution for nonlinear problems and problems with more general boundary conditions.
In this paper we present the preliminary investigation of a novel, highorder immersed boundary method (IBM). Our approach is similar in spirit to the fictitious domain method of Glowinski and He [13] and optimizationbased interface treatment in [18]. The key idea is to formulate the IBM as an inverse problem in which unknown boundary fluxes, defined on an approximate easytogenerate boundary, are used to satisfy the desired boundary conditions on the true boundary. One strength of this inverseproblem formulation of the IBM is that it is straightforward both conceptually and in its implementation.
The remainder of the paper is organized as follows. Section 2 describes the basic inverse IBM formulation for both the Laplace and linear advection equations. This section also discusses the illposedness of the basic formulation. Section 3 presents an investigation of a model problem, which we use to better understand the nature of the illposedness. Regularization of the inverse IBM is discussed in Section 4, including a parameterfree regularization based on penalizing the control against the state. The proposed method is demonstrated in Section 5 using on a discontinuous Galerkin discretization. A summary and discussion are provided in Section 6.
2 Inverse IBM problem formulation and illposedness
This section presents and analyzes the basic formulation of the proposed IBM in the context of the Laplace and linear advection equations. We show that there is a unique solution when the inverse IBM is applied to the Laplace equation, but that the solution does not depend continuously on the data. In the case of the linear advection equation, we argue that not even solution uniqueness holds.
2.1 Application to the Laplace equation
Consider the Laplace equation on the open, bounded domain , with Dirichlet boundary conditions applied on the boundary : \cref@addtoresetequationparentequation
(1a)  
(1b) 
Suppose is a geometrically complex domain for which we do not want to generate a conforming mesh. Consequently, we introduce a geometrically simpler, bounded domain , with boundary . Figure 1 depicts an example of the domain , an encompassing domain , and their respective boundaries.
Having described this generic geometric setup, we now turn to the key question: how do we select the boundary conditions on such that the conditions (1b) are satisfied, or at least approximately satisfied, on ? Our answer to this question is to define the following PDEconstrained optimization problem, which is given in variational form: find and that satisfy
(2)  
s.t. 
Remark 1
Notice that the (variational) Laplace equation is posed on the larger domain . Furthermore, the boundary condition imposed on is defined by , rather than . The variable is the control variable used to minimize the cost functional in order to (indirectly) satisfy the boundary condition on .
Remark 2
Formulation (2) defines a boundarycontrol problem, and it is the prototype method that we will explore further in the rest of the paper; however, in some cases it is possible to define an inverse IBM using distributed control instead of boundary control. For example, if is defined such that the immersed geometry is bounded and coincides with , then a distributed control can be defined over the region interior to the immersed geometry. A distributedcontrol approach is the basis of the fictitious domain method in [13].
2.1.1 Existence and uniqueness of the inverse IBM Laplace solution
In this section we show that, under mild assumptions, the inverse problem (2) has a solution that is unique. This result follows from the Riesz representation theorem by establishing that the objective is associated with a symmetric, coercive bilinear form in the control space .
We will work in the reduced space of the control by eliminating the state from the PDEconstrained optimization problem 2. To this end, we first show that the state can be expressed as a linear function of the control.
Lemma 1
Let be a Lipschitz domain and define the operator as follows: for , find such that
Then is linear and continuous.
The proof is provided in Appendix A. In light of Lemma 1, we will use the notation to indicate the dependence of the state on the control.
Next, we express the objective, in terms of the control variable, as the sum of quadratic, linear, and constant terms:
where . The forms and are given by
(3)  
(4) 
respectively. The following lemma guarantees that the forms and have the properties we need in order to apply the Riesz representation theorem. The proof is somewhat lengthy, so it is relegated to Appendix B.
Lemma 2
With the groundwork laid, we can state and prove the main existence and uniqueness result.
Theorem 1
Proof
The theorem follows from, for example, [9, Theorem 6.11] and is an application of the Riesz representation theorem [33]. We need only verify that the assumptions of the Riesz theorem are met. First, is a Banach space. Second, the bilinear form is symmetric, continuous, and coercive by Lemma 2. Finally, we have that is a continuous linear form by the same lemma.
2.1.2 Illposedness in the context of the Laplace equation
We have shown that the inverse IBM formulation of the Laplace equation has a unique solution. To establish wellposedness, we need to additionally show that the solution depends continuously on the data that defines (2); unfortunately, as we now show, there exists inverse IBM problems for which this is not the case.
Consider an inverse IBM problem whose target domain is the unit disc, , and let the boundary condition on be given by the following sinusoid:
where is a positive integer. Now, suppose the larger, enclosing domain is also a disc, but one with radius ; thus, . It is easy to verify that the solution to the Laplace equation on that satisfies is given by
Consequently, the corresponding (unique) boundary control that solves the inverse IBM problem (2) is the trace of on :
This control can be made arbitrarily large by choosing sufficiently large; however, as the boundarycondition data on tends to zero. We conclude that the solution does not depend continuously on the data and the inverse IBM is not wellposed for the Laplace equation.
2.2 Application to the linear advection equation
We now turn our attention to the linear advection equation, which is an important model problem for inviscid fluid flows and hyperbolic PDEs more generally. Unlike the Laplace equation, we will see that the basic inverse IBM formulation does not even admit a unique solution for linear advection.
The linear advection equation on an open, bounded domain , with smooth boundary , is given by \cref@addtoresetequationparentequation
(5a)  
(5b) 
where is the vectorvalued velocity field. Boundary conditions are imposed on the inflow boundary , where is the velocity component normal to the boundary and is the unit, outwardpointing normal vector on .
Remark 3
We will assume that (5) is wellposed for the given data and geometry. Furthermore, in order to apply the inverse IBM formulation, we will assume that the velocity field can be extended to the larger computational domain . That is, we assume that there exists such that for all .
Next, we formulate an inverse IBM statement, analogous to (2), for the linear advection equation. Specifically, we replace the variational Laplace equation with the variational advection equation in (2) and obtain the following problem: find and that satisfy
(6)  
s.t. 
where and . Figure 2 provides an example of the domains and boundaries that appear in problem (6). We will also refer to this figure when explaining the illposedness that affects the inverse IBM in the context of the advection equation.
2.2.1 Nonuniqueness and other difficulties
Consider the point labeled in Figure 2. A compact control perturbation centered at determines the value of along the characteristics downstream of via the boundary condition on . However, if the perturbation is sufficiently local, then the corresponding characteristics will not pass through and the perturbation will have no impact on the objective. Consequently, such control perturbations are nonunique in the context of (6), i.e. the inverse IBM problem is illposed for linear advection.
To make matters worse, nonlinear hyperbolic PDEs present an additional difficulty for the proposed inverse IBM formulation. Namely, the flow field itself often depends on the solution, so the definition of and is not necessarily known a priori. For the discretized version of (6), this uncertainty in the set means that the dimension of the boundary control would need to change dynamically during the solution process.
Remark 4
This practical difficulty is exhibited by the Euler equations of gas dynamics. For example, if is a nopenetration wall in the Euler equations, it will have a known number of incoming and outgoing characteristics (one each); however, the surface may require additional incoming/outgoing characteristics in order to satisfy the nopenetration condition.
Our solution to this practical issue is to define the control on the entire boundary, as was done for the Laplace equation. On the one hand, this introduces another source of nonuniqueness: the value of the control on — for example, in Figure 2 — has no influence on the state, and, therefore, no influence either satisfying the constraint or minimizing the objective in (6). On the other hand, this form of nonuniqueness is easily eliminated with the regularization we propose in Section 4.
3 An investigation of conditioning of the discrete problem
It is clear that the basic inverse IBM formulation must be regularized to be useful in practice; however, before we discuss potential regularizations in Section 4, we want to better understand the sources of illconditioning in the finitedimensional case that ultimately lead to illposedness as the mesh size . To this end, we will model the discretization of the Laplace inverse IBM on a particular domain, and we use this model to investigate the conditioning of the finitedimensional version of (2).
3.1 Description of the model problem
We consider a model problem similar to the one described in Section 2.1.2. The computational domain is the unit disk, , and the immersed boundary is an ellipse with semimajor and semiminor lengths and , respectively. The model geometry is illustrated in Figure 4. We have chosen a circle for because it allows us to take advantage of Green’s function theory.
Our goal is to construct a model for the reduced Hessian , where , denotes the control evaluated at discrete locations along . We will then use to investigate the conditioning of the finitedimensional Hessian.
Thus, the first step is to get an expression for that we can use to evaluate the bilinear form . Using the appropriate Green’s function for the disk [31, §7.1.2], we obtain the following expression for , where is the control value on :
(7) 
To evaluate the bilinear form in a pointwise sense over , we use Diracdelta control variations that are parameterized based on uniformlyspaced points on the unit circle:
Substituting the distributions for and into the bilinear form , with and defined by (7), we obtain the following expression for :
(8) 
Note that, while both and appear in the integrand above, the point is constrained to lie on the ellipse .
Before we use (8) to explore the condition number of the discrete Hessian, it is worthwhile to plot versus for fixed , because these plots will help explain the trends in the condition number. Two such rows are shown in Figure 4: one corresponding to and one corresponding to . Each entry is evaluated for three different immersed ellipse boundaries with a fixed semimajor axis of and progressively smaller semiminor axes . We have approximated the integral in (8) using the trapezoid rule with intervals.
The trend is clear from the figure. There is little impact on the Hessian row near the fixed, semimajor axis at . In contrast, decreasing the semiminor axis stretches horizontally and compresses it vertically. The spreading of this row of the Hessian hints at conditioning issues on the horizon, since relatively “flat” rows are difficult to distinguish from one another, i.e. the rows are nearly linearly dependent.
3.2 Circular immersed boundary
Consider a circular immersed boundary , with . Figure 5(a) plots the condition number of the reduced Hessian versus the Hausdorff distance between and for three different controlmesh spacings . In the present case, where both and are circles, the Hausdorff distance is given by .
For fixed , the condition number grows exponentially with the distance between and . This growth is due to the diffusing of the analytical Hessian as the immersed boundary gets further from the domain boundary, as described earlier. Note that, for a fixed , we also observe a strong dependence on mesh spacing.
Figure 5(a) suggests that it is possible to keep the condition number of fixed as we refine the mesh, provided is also reduced. This relationship is clearly illustrated in Figure 5(b), which plots the condition number of versus the ratio . This figure is important because it shows that different Hausdorff distances essentially collapse onto the same curve, and it is the nondimensional ratio that determines the condition number. As long as is constant, the condition number of the reduced Hessian should also remain bounded.
3.3 Elliptical immersed boundary
We conclude our investigation of the model problem by considering the influence of differing semimajor and semiminor axis lengths. Figure 6(a) plots the reduced Hessian condition number versus for three different fixed values of . In all cases the mesh spacing is fixed at . In general, the condition number remains stable as increases until , i.e. the semiminor and semimajor axes switch, at which point the exponential growth is similar to that observed for circular . This behavior is consistent with the Hausdorff distance determining the condition number for fixed : for the Hausdorff distance is , which is constant for fixed .
Figure 6(b) plots the condition number of versus the nondimensional mesh spacing . We have included curves corresponding to several different majorminor offset ratios, . We observe that remains the key parameter that determines the condition number.
3.4 Discussion
The results of the above studies suggest that the discretized inverse IBM will be nonsingular, at least in the case of the Laplace equation, if the ratio remains bounded. This is not surprising, since this condition implies that , so the inverse IBM problem becomes essentially equivalent to the direct problem with .
However, there are good reasons why bounding alone may not be sufficient in practice.

For a desired mesh resolution , it may be challenging to generate a mesh for which the distance between and is sufficiently small to avoid large condition numbers; indeed, in the limit as , we are back to generating a conformal mesh. Thus, while bounding the ratio may ensure wellposedness in the limit , it may not prevent illconditioned discrete systems.

The nonunique solutions that the inverse IBM admits for linear advection are not excluded by bounding . Thus, without further intervention, the discretized IBM will produce a singular system when applied to the linear advection equation.
For the reasons listed above, we need to consider some form of regularization for the inverse IBM.
4 Regularization
We have seen that the basic formulation of the inverse IBM is illposed for both the Laplace and linear advection equations, albeit for different reasons. Furthermore, while bounding the ratio seems to resolve illposedness for the Laplace equation, this strategy does not suffice for the linear advection equation nor does it prevent high condition numbers for the Laplace equation. Therefore, the basic formulation must be modified if it is to be useful in practice.
Regularized formulations are commonly used to ensure inverse problems are both wellposed and wellconditioned, so this is an obvious strategy to examine for the inverse IBM. In this section, we briefly review two popular regularizations used for inverse problems and argue why they are poorly suited for the inverse IBM. We then introduce a novel regularization for the inverse IBM that is well suited for discretizations that use weakly imposed boundary conditions.
4.1 Review of regularization approaches for inverse problems
The popular regularizations that we review here are the Tikhonov and total variation (TV) diminishing regularizations. Both methods modify the objective function by adding a (scaled) convex term, which we denote by . Thus, the modified objective is given by
(9) 
where is a regularization parameter.
Remark 5
The regularized optimization problem based on corresponds to a biobjective problem in which implicitly determines a particular Pareto optimal solution. Relatively small values of select solutions that emphasize the boundarycondition accuracy at the cost of illconditioning, while large produce better conditioned problems whose solutions may have significant mismatch in the boundary condition.
In the case of Tikhonov regularization applied to the inverse IBM, we have
(10) 
where is a constant function that can be used to reduce the detrimental effect that the regularization has on solution accuracy. The TVdiminishing regularization takes the form
where and denotes the gradient of the control with respect to a parameterization of . The constant is a small parameter that is included to ensure the regularization is differentiable at points where . TVbased regularization was originally introduced in digital image processing for noise removal [35].
The Tikhonov term is quadratic, strongly convex, and it tends to produce smooth solutions. In contrast, the TVdiminishing term tends to produce piecewise constant solutions that may arise, for example, in the presence of discontinuous physical properties. Note that the TV regularization is not strongly convex — any constant minimizes — so it may not be sufficient to guarantee wellposedness, in general. Furthermore, the TV regularization is not quadratic, so iterative solution methods are required even for linear state equations.
For both regularizations, choosing an appropriate value for the parameter can be challenging: it should be sufficiently large to stabilize the problem, but not too large to significantly impact solution accuracy; see [40] for an error estimation of the solution with respect to . Some common techniques to select include the curve method [15] and the Morozov discrepancy principle [36].
In the case of the inverse IBM, selecting the parameter is further complicated by the requirements of solution convergence as . To explain this complication, we consider a conforming finiteelement discretization of the Laplace inverse IBM problem (2) with Tikhonov regularization: find and that satisfy
(11)  
s.t. 
where , and are appropriate finiteelement spaces.
Theorem 2
Let and define the unique solution to (2), and let and be the solution to (11). Let denote the nominal element size, and suppose the solution error on satisfies for some integer and sufficiently small . Furthermore, suppose the discrete control is bounded below, , for sufficiently small . Then there exists such that the regularization parameter satisfies
where is independent of .
Proof
Without loss of generality, consider Tikhonov regularization with . Then, the firstorder optimality conditions for (11) in the reduced space are
where the prime notation denotes Fréchet differentiation with respect to the term in brackets. In particular, for , we have
where follows from the linear dependence of on . Rearranging the above equality and applying CauchySchwarz, we find
The assumptions on the solution error and imply that there exists an such that both and , for some and all . Using these bounds in the above inequality we arrive at
The result follows with .
The implication of Theorem 2 is that, if we have optimal solution error, then the positive effect that Tikhonov regularization has on conditioning must necessarily vanish as . In other words, the illconditioning, or illposedness in the case of the advection equation, will return as the mesh is refined. A similar conclusion can be drawn with TVbased regularization. This motivates our search for a regularization that is compatible with highorder accuracy, but whose Hessian does not vanish with mesh refinement.
Remark 6
Theorem 2 gives a necessary but not sufficient condition on . Indeed, our numerical experiments (not reported here) indicate that in some cases optimal solution accuracy is achieved only if .
4.2 Regularization via boundarycondition penalization
Our proposed regularization for the inverse IBM is inspired by weakly imposed boundary conditions. To motivate the regularization, recall our discussion of the linear advection problem in Section 2.2. If the boundary control is defined on all of , we explained why a variation in on the downwind boundary does not influence the state. While the control is intended to define the boundary value, we can reverse the causation on and insist that determines . More generally, for weakly imposed boundary conditions, we can seek a solution for which the difference between state and control is small on the entire boundary . That is, we consider a regularization term of the form
(12) 
While the regularization (12) is compatible with strongly imposed boundary conditions, it is better suited for weakly imposed boundary conditions. For example, it can be used to eliminate illposedness due to defining on all of when the inverse IBM is applied to the linear advection equation; however, the regularization term has no further benefit when the boundary conditions are strongly imposed. In contrast, we have found that (12) improves conditioning significantly for both the Laplace and linear advection equations when weakly imposed boundary conditions are used; see Section 5.
Unlike Tikhonov and TVdiminishing regularizations, the discretized version of (12) tends to zero as . Consequently, any (constant) choice of the parameter in (9) is compatible with an optimal solution convergence rate. For example, consider the Laplace equation under the assumptions in Theorem 2 and suppose the optimal convergence rate is . Using the regularization (12), the firstorder optimality conditions in the reduced space imply (see the proof of Theorem 2 for further details)
Previously, with Tikhonov regularization, it was necessary (but not sufficient) to have in order to satisfy the firstorder optimality conditions as . In contrast, the proposed regularization can rely on the regularization itself to satisfy the above optimality condition, since both terms in the condition are asymptotically similar. That is, assuming optimal convergence of the state to the boundary value on , we have .
Remark 8
Based on the above discussion, any constant choice for is consistent with optimal solution convergence. We adopt for all subsequent numerical experiments.
5 A numerical investigation of accuracy and conditioning
The purpose of the following numerical experiments is to i) investigate the conditioning of the inverse IBM, and ii) verify optimal solution accuracy. For these studies we consider the steady, constantcoefficient, advectiondiffusion equation:
(13)  
where is the advection velocity and is the nonnegative diffusion coefficient. The problem domain is the unit disk, i.e., .
In the numerical experiments, three sets of coefficients are chosen to model different physics. These sets are as follows:

and for a pure advection problem;

, for a Poissontype diffusion problem;

, for an advectiondiffusion problem.
The boundary conditions are imposed on for pure diffusion and advectiondiffusion, and they are imposed on for pure advection.
We use the method of manufactured solutions [34], in which an exact solution is selected a priori and used to determine and in (13). For the following numerical experiments we use the exact solution
(14) 
5.1 Discretization of the inverse IBM formulation
5.1.1 Discretization of the advectiondiffusion equation
We use a standard discontinuous Galerkin (DG) finiteelement method to discretize the advectiondiffusion equation (13). The symmetric interior penalty Galerkin (SIPG) method [3, 37] is used to discretize the diffusion term while standard upwinding is used for the advection part; see, for example, [16]. For completeness, a detailed description of the discretization is provided in Appendix C.
Remark 9
While we have adopted a DG finiteelement method here, we emphasize that the proposed inverse IBM formulation is agnostic to the choice of discretization with the exception of the regularization (12), which requires weakly imposed boundary conditions.
5.1.2 Discretization of the objective function
Having described the discretization of the PDE constraint, we now turn our attention to the objective function .
Let us first consider the discretization of , and recall that the boundary is a circle in the following results. We uniformly divide the circle into line segments of length . On each line segment we introduce Gauss quadrature points, and we construct an interpolation to each quadrature point based on the finite element that contains the point. We can then use the interpolated values to evaluate on . The value of at the quadrature points can be determined from the manufactured solution (14).
For the following results, we employ a quadrature rule that is exact for polynomials of degree for each segment on , where is the degree of the Lagrange bases used in the discretization (18). In addition, we consider different values of when studying the condition number of the reduced Hessian; however, in the solutionconvergence study, we choose such that .
Remark 10
The size of the line segments, , as well as the number of quadrature points on each line segment, can be varied independently of the finiteelement discretization; however, as we shall see, there is a limit to how large the ratio can be made, for a given quadrature rule, before it impacts the conditioning of the discrete system.
Finally, the regularization term can be discretized directly using the trace of and the value of on the edges .
5.1.3 Definition of the expanded domain
The domain is determined using the following process. We begin by generating a uniform background triangular mesh of sufficient size to contain the unit circle. Then is defined as the set of all elements that are either interior to the disk or are intersected by . Figure 7 illustrates the coarsest two meshes and corresponding produced using this procedure. The nominal element size is defined as the square root of the triangle areas.
Remark 11
The proposed construction for ensures that the ratio remains bounded as , which is what we want in practice. On the other hand, the construction makes it impossible to vary the Hausdorff distance independently of , as we did for the model problem in Section 3; however, we can still investigate the impact of the ratio of to the spacing of the degrees of freedom by varying the finiteelement polynomial degree . We will present results of such a study in Section 5.3.
Remark 12
Recall the illposedness that affects the inverse IBM applied to the linear advection equation, e.g. point in Figure 2. We hypothesize that the proposed construction of helps avoid illconditioning produced by such points, since any element containing these points is necessarily coupled to elements interior to the domain ; however, further study is necessary to confirm this claim.
5.2 Solution of the discretized inverse IBM system
We now briefly discuss the solution of the discretized inverse IBM formulation, as well as the construction of the reduced Hessian matrix used to study the conditioning of the system.
Let be a triangularization of , let denote the basis for the discrete state space , and let be the basis for the discretized control space ; see Appendix C for definitions of and . The finiteelement state and control can be represented as
respectively, where and are the nodal coefficients. Substituting these expressions for and into the advectiondiffusion trilinear form (18) and the discretized objective, the discretized inverse IBM problem can be expressed as the following finitedimensional optimization problem:
(15)  
s.t. 
The optimization problem (15) is a quadratic program since it has a quadratic objective and a linear constraint. Furthermore, one can show that the objective is convex in and . Consequently, the solution to (15) is equivalent to the solution to the following saddlepoint problem:
(16) 
where are the Lagrange multipliers associated with the constraint; in the present context, the multipliers are equivalent to the nodal coefficients of the adjoint. We solve the linear system (16) using a sparse direct solver.
Remark 13
In practice, the efficient solution of (15) or (16) would require more sophisticated algorithms, such as sparse iterative solvers with specialpurpose preconditioners and Newton’s method for nonlinear PDEs. The application of such algorithms to (15) is the subject of ongoing research — see, for example, the preliminary results in [42] — but is beyond the scope of the current work.
The optimization statement (15) and linear system (16) are fullspace formulations, since the control, state, and multipliers are solved simultaneously; see, for example, [6, 7, 2]. Alternatively, we can eliminate the state and multipliers to arrive at the following reducedspace optimization statement:
where is a constant, is the reduced gradient and
(17) 
is the reduced Hessian. We will study the conditioning of the reduced Hessian in the following sections, since this provides an indication of the difficulty of solving (16) using iterative methods. Furthermore, we can use to verify our conclusions regarding the model problem from Section 3.
5.3 Conditioning of the reduced Hessian
5.3.1 Conditioning for the diffusion equation
We begin our investigation of the conditioning of , given in (17), with the Poissontype diffusion case, that is, with and in (13).
Recall that the circular immersed boundary is discretized into line segments of length , and that is independent of the volume mesh size . However, if there are control degrees of freedom, we expect that would be a necessary condition for the reduced Hessian to be nonsingular. This intuition is supported by the results in Figure 8(a), which plots the condition number of for the basic (unregularized) inverse IBM versus the ratio . Results for a range of element sizes are shown, with elements adopted in all cases. We see that the condition number is relatively constant until a threshold of is reached, at which point becomes singular to working precision.
Figure 8(b) shows the analogous results for the regularized objective. The figure highlights two benefits of the regularization for diffusive PDEs. First, rather than abruptly becoming singular, the condition number of gradually increases when passes a threshold of approximately one. Second, for values of below this threshold, the condition number of is approximately two orders of magnitude smaller than when no regularization is used.
Next, we would like to study the effect of the Hausdorff distance on the conditioning of like we did for the model problem in Section 3. Unfortunately, we cannot vary the Hausdorff distance independently from because of how is constructed; see the discussion in Section 5.1.3. We can, however, mimic the study from Section 3 by varying the degree of the polynomial basis for fixed and . That is, whereas was the control spacing in the earlier study, here we can use to define the nominal control spacing.
Figure 9 plots the condition number of the reduced Hessian versus the ratio for both the unregularized and regularized Poisson problems. We consider , , and degree finiteelement basis functions and a range of mesh sizes and a (dependent) range of Hausdorff distances. The results for the unregularized problem, plotted in Figure 9(a), behave similar to the modelproblem results shown in Figures 5(b) and 6(b): as the control spacing becomes small relative to the Hausdorff distance, the condition number “blows up.” The regularized problem — see Figure 9(b) — displays qualitatively the same pattern, but the condition number is two orders of magnitude smaller for the same value of . For both the unregularized and regularized cases, we observe no significant dependence on itself. Overall, the results in Figure 9 confirm the conclusions drawn earlier.
5.3.2 Conditioning for the linear advection equation
Here we investigate the conditioning of the reduced Hessian of the inverse IBM in the context of the linear advection equation. We can only consider the regularized variant of the inverse IBM, because, unlike the Poisson problem, the unregularized advection problem is singular to working precision.
Figure 10(a) plots the condition number of versus the size ratio for the discretization. We see that the regularized advection problem behaves more like the unregularized diffusion problem. In particular, there is a threshold, here, above which the reduced Hessian becomes singular. Consequently, it is critical to use a sufficiently fine discretization of when the inverse IBM is applied to hyperbolic PDEs.
Figure 10(b) shows how the condition number of the reduced Hessian responds to changes in the ratio . As with the diffusion problem, the condition number grows unbounded as the control spacing becomes small relative to the Hausdorff distance. Furthermore, the magnitude of the condition number is larger for linear advection for the same value of . This suggests that the inverse IBM, with the proposed regularization, may be limited to modest values of for hyperbolic problems.
5.4 Convergence study
We conclude the numerical experiments with a study of solution accuracy and convergence rate. To estimate the asymptotic convergence rate, we use a sequence of five uniformly refined triangular meshes. To obtain the next mesh in the sequence, each element is subdivided into four. Thus, the element sizes of the finer meshes are , , and . The coarsest two meshes, together with the immersed boundary , are shown in Figure 7. We can see that the physical boundary intersects mesh elements at different locations, including vertices; in our experience, the method is robust with respect to such coincident nodes.
The solution contours using and basis functions on the coarsest mesh are compared against the exact contours in Figure 11(a) and 11(b), respectively. We can see that the discrete solution agrees well with the exact solution. Furthermore, as expected for a smooth solution, the higherorder approximation produces qualitatively better results on the same mesh.
For a quantitative assessment of accuracy, it is common to use the error to measure the accuracy of a finiteelement solution ; however, evaluating the norm exactly on is not straightforward, because the standard quadrature rules do not apply on the elements cut by the boundary. The approach adopted in this paper is to set the solution error to zero at all quadrature points that lie in .
Figure 12 plots the solution error versus element size for the specific advection, diffusion, and advectiondiffusion problems defined earlier. For all problems, the inverse IBM achieves the optimal convergence rates of .
6 Summary and Discussion
First, a brief summary. The proposed inverse IBM introduces a control on the boundary of an expanded domain that encompasses the target domain . In the basic formulation, the value of the control is determined by minimizing the mismatch in the boundary condition on the true boundary subject to the state satisfying the desired PDE on the expanded domain . This basic formulation is illposed for the Laplace and linear advection equations: for the former, the solution does not depend continuously on the data, and for the latter the solution is not uniquely determined by the data. The illposedness is addressed by ensuring the expanded domain tends to the target domain and by including a controlstate penalty term in the objective. Applied to a discontinuous Galerkin finiteelement discretization of the advection, advectiondiffusion, and diffusion equations, the inverse IBM formulation achieves optimal solution convergence rates, and the reduced Hessian remains wellconditioned.
We believe the inverse IBM has a number of attractive features. The approach is agnostic to the underlying discretization, so it can be applied to finite difference, finite volume, and finiteelement discretizations. Furthermore, the method is compatible with highorder discretizations as our results demonstrate. Finally, once the discretization is chosen, the implementation is straightforward in the sense that there are no “corner cases” that require special treatment; all that is required is an interpolation/projection operator from to the boundary .
There is no “free lunch,” and the attractive features of the inverse IBM come at a price. Specifically, the method trades computational cost for ease of implementation and accuracy. If there are state variables and control variables, the inverse IBM has variables. This should be contrasted with other approaches that have only state variables. In addition, the inverse IBM is a PDEconstrained optimization problem, which requires the solution of a linear saddlepoint problem for linear PDEs and a nonlinear system in general. Nevertheless, we believe the increased computational cost is well worth the reduced human time required in mesh generation around complex configurations.
The proposed inverse IBM is promising, but there remain many challenges that must be addressed before it can be used on practical computational fluid dynamics problems. Foremost among these challenges is the robust and efficient solution of the nonlinear PDEconstrained optimization problem associated with the inverse IBM. A related issue is the need for specialized preconditioners for the saddlepoint systems that arise in the Newton iterations. Finally, there is the question of how the method can be extended to unsteady simulations.
Appendix A Proof of Lemma 1
NameIgnored (Lemma 1.)
Let be a Lipschitz domain and define the operator as follows: for , find such that
Then is linear and continuous.
Proof
The operator is linear because is bilinear and the boundary conditions are linear.
To show that is continuous, we introduce the extension operator , which is the right inverse of the trace operator, i.e. it extends a function on the trace space into the domain. Using this operator, we define , which satisfies the variational form of the Poisson equation
Using the above equation, together with the PoincaréFriedrichs and CauchySchwarz inequalities, we find
where is some constant. Therefore, for we have the following bound on :
where the last inequality follows from the boundedness of the extension operator [27, Theorem 3.37]. Note, if then , and we can use the same bound on . In conclusion, we have shown , so is continuous.
Appendix B Proof of Lemma 2
NameIgnored (Lemma 2.)
Proof
The form is clearly symmetric, so, to prove that it is bilinear, it is sufficient to show that is linear in its first argument. For and