An inverse problem formulation of the immersed boundary methodSubmitted 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.

# 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.

Jianfeng Yan Department of Mechanical, Aerospace, and Nuclear Engineering, Rensselaer Polytechnic Institute, Troy, New York, United States (, ).    Jason E. Hicken1
2footnotemark: 2
###### Abstract

We formulate the immersed-boundary 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 target-domain boundary. We begin by investigating a naïve problem formulation that we show is ill-posed: 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 ill-posedness 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 parameter-free 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 advection-diffusion equations using a high-order 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.

i
\newsiamremark

remarkRemark \newsiamremarkhypothesisHypothesis \newsiamthmclaimClaim \headersinverse immersed boundary methodJ. Yan and J.E. Hicken

mmersed-boundary method, inverse problem, PDE-constrained optimization

{AMS}

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 curved-element, anisotropic meshes around complex geometries, which has limited the adoption of high-order methods in industry.

One way to address the meshing bottleneck is to develop a high-order 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 immersed-interface, or cut-cell, methods. IBMs [29] impose the boundary conditions indirectly through a body force or modified boundary flux. Immersed-interface 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]. Immersed-interface 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 non-trivial task whose complexity is arguably on par with geometry-conforming mesh generation. Furthermore, the cutting process invariably creates relatively small and non-standard 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 finite-difference/volume methods see [30] and Sections 3 and 4.1 of [28]. A review of penalty-based IBMs used in finite-element 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, advection-dominated problems common in the aerospace industry. We believe this is primarily due to their limited accuracy; most IBMs are first-order accurate in practice, because the underlying solution is non-smooth at the interface that represents the boundary [30].

While uncommon, high-order 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 ghost-fluid method [11, 10, 21, 17] to high-order 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 right-hand side of the discretized equations. The correction is independent of the solution for linear problems with Dirichlet boundary conditions, so conditioning of the left-hand 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, high-order immersed boundary method (IBM). Our approach is similar in spirit to the fictitious domain method of Glowinski and He [13] and optimization-based 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 easy-to-generate boundary, are used to satisfy the desired boundary conditions on the true boundary. One strength of this inverse-problem 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 ill-posedness of the basic formulation. Section 3 presents an investigation of a model problem, which we use to better understand the nature of the ill-posedness. Regularization of the inverse IBM is discussed in Section 4, including a parameter-free 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 ill-posedness

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

 ∇2u =0, ∀x∈Ω, (1a) u =uΓ, ∀x∈Γ. (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 PDE-constrained optimization problem, which is given in variational form: find and that satisfy

 infc,u J(u)≡12∫Γ(u−uΓ)2dΓ, (2) s.t. b(u,v)≡∫~Ω∇u⋅∇vdΩ=0,∀v∈H10(~Ω).
###### 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 boundary-control 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 distributed-control 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 PDE-constrained 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

 b(u,v)=∫~Ω∇u⋅∇vdΩ=0,∀v∈H10(~Ω).

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:

 J(c)=J(u(c))=12a(u(c),u(c))−l(u(c))+J(0)

where . The forms and are given by

 a(c,d)≡∫Γu(c)u(d)dΓ, (3) andl(c)≡∫Γu(c)uΓdΓ, (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

Let be an open subset. The form defined by (3) is bilinear, symmetric, continuous, and -coercive. The form defined by (4) is linear and continuous.

With the groundwork laid, we can state and prove the main existence and uniqueness result.

###### Theorem 1

Under the assumptions of Lemmas 1 and 2, there is a unique solution to (2). That is, there is a unique element such that

 J(c)=infc′∈H12(~Γ)J(c′).

###### Proof

The theorem follows from, for example, [9, Theorem 6.1-1] 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 Ill-posedness 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 well-posedness, 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:

 uΓ(θ)=1nsin(nθ),∀θ∈[0,2π),

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

 u(r,θ)=rnnsin(nθ),∀(r,θ)∈~Ω.

Consequently, the corresponding (unique) boundary control that solves the inverse IBM problem (2) is the trace of on :

 c(θ)=Rnnsin(nθ),∀θ∈[0,2π).

This control can be made arbitrarily large by choosing sufficiently large; however, as the boundary-condition data on tends to zero. We conclude that the solution does not depend continuously on the data and the inverse IBM is not well-posed 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

 ∇⋅(λu) =0, ∀x∈Ω, (5a) u =uΓ, ∀x∈Γ−, (5b)

where is the vector-valued velocity field. Boundary conditions are imposed on the inflow boundary , where is the velocity component normal to the boundary and is the unit, outward-pointing normal vector on .

###### Remark 3

We will assume that (5) is well-posed 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

 infc,u J(u)≡12∫Γ−(u−uΓ)2dΓ, (6) s.t. b(u,v)≡∫~Ωv∇⋅(~λu)dΩ=0,∀v∈L2(~Ω),

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 ill-posedness 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 non-unique in the context of (6), i.e. the inverse IBM problem is ill-posed 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 no-penetration 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 no-penetration 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 ill-conditioning in the finite-dimensional case that ultimately lead to ill-posedness 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 finite-dimensional 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 semi-major and semi-minor 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 finite-dimensional 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 :

 u(c)=u(r,θ)=12π∫2π0(1−r2)c(~θ)1+r2−2rcos(~θ−θ)d~θ. (7)

To evaluate the bilinear form in a point-wise sense over , we use Dirac-delta control variations that are parameterized based on uniformly-spaced points on the unit circle:

 ci=c(θi)=δ(θ−θi),whereθi=2πi/K,i=1,2,…,K.

Substituting the distributions for and into the bilinear form , with and defined by (7), we obtain the following expression for :

 Hij=H(θi,θj)a(ci,cj)=1(2π)2∫Γ(1−r2)[1+r2−2rcos(θi−θ)](1−r2)[1+r2−2rcos(θj−θ)]dΓ. (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 semi-major axis of and progressively smaller semi-minor 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, semi-major axis at . In contrast, decreasing the semi-minor 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 control-mesh 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 semi-major and semi-minor 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 semi-minor and semi-major 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 major-minor 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 well-posedness in the limit , it may not prevent ill-conditioned discrete systems.

• The non-unique 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 ill-posed for both the Laplace and linear advection equations, albeit for different reasons. Furthermore, while bounding the ratio seems to resolve ill-posedness 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 well-posed and well-conditioned, 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

 Jreg(u,c)=J(u)+αS(u,c), (9)

where is a regularization parameter.

###### Remark 5

The regularized optimization problem based on corresponds to a bi-objective problem in which implicitly determines a particular Pareto optimal solution. Relatively small values of select solutions that emphasize the boundary-condition accuracy at the cost of ill-conditioning, 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

 STik(c)≡12∫~Γ(c−c0)2dΓ, (10)

where is a constant function that can be used to reduce the detrimental effect that the regularization has on solution accuracy. The TV-diminishing regularization takes the form

 STV(c)≡∫~Γ√∇~Γc⋅∇~Γc+ϵdΓ,

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 . TV-based 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 TV-diminishing 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 well-posedness, 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 finite-element discretization of the Laplace inverse IBM problem (2) with Tikhonov regularization: find and that satisfy

 minch,uh Jreg(uh,ch)≡12∫Γ(uh−uΓ)2dΓ+α2∫~Γ(ch−c0)2dΓ (11) s.t. bh(uh,vh)=∫~Ω∇uh⋅∇vhdΩ=0,∀vh∈V′h(~Ω),

where , and are appropriate finite-element 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

 α≤Mhp+1,∀h≤h⋆,

where is independent of .

###### Proof

Without loss of generality, consider Tikhonov regularization with . Then, the first-order optimality conditions for (11) in the reduced space are

 J′reg[ch](dh)=∫Γ(uh(ch)−uΓ)u′h[c]dhdΓ+α∫~ΓchdhdΓ=0,∀dh∈Vh(~Γ),

where the prime notation denotes Fréchet differentiation with respect to the term in brackets. In particular, for , we have

 ∫Γ(uh(ch)−uΓ)uh(ch)dΓ+α∫~Γc2hdΓ=0,

where follows from the linear dependence of on . Rearranging the above equality and applying Cauchy-Schwarz, we find

 α∥ch∥2~Γ≤∥uh(ch)−uΓ∥Γ∥uh(ch)∥Γ.

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

 α≤∥uh(ch)−uΓ∥Γ∥uh(ch)∥Γ∥ch∥2~Γ≤Λhp+1∥uΓ∥ΓB2,∀h

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 ill-conditioning, or ill-posedness in the case of the advection equation, will return as the mesh is refined. A similar conclusion can be drawn with TV-based regularization. This motivates our search for a regularization that is compatible with high-order 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 boundary-condition 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

 S(u,c)=12∫~Γ(u−c)2dΓ. (12)
###### Remark 7

The regularization (12) is analogous to a Tikhonov regularization (10) with .

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 ill-posedness 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 TV-diminishing 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 first-order optimality conditions in the reduced space imply (see the proof of Theorem 2 for further details)

 ∫Γ(uh(ch)−uΓ)uh(ch)dΓ+α∫~Γ(uh−ch)uh(ch)dΓ=0.

Previously, with Tikhonov regularization, it was necessary (but not sufficient) to have in order to satisfy the first-order 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, constant-coefficient, advection-diffusion equation:

 ∇⋅(λu−μ∇u) =f, ∀x∈Ω, (13) u =uΓ, ∀x∈Γ′,

where is the advection velocity and is the non-negative 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 Poisson-type diffusion problem;

• , for an advection-diffusion problem.

The boundary conditions are imposed on for pure diffusion and advection-diffusion, 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

 u(x,y)=ex+ysin(πx)sin(πy). (14)

### 5.1 Discretization of the inverse IBM formulation

#### 5.1.1 Discretization of the advection-diffusion equation

We use a standard discontinuous Galerkin (DG) finite-element method to discretize the advection-diffusion 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 finite-element 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 solution-convergence 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 finite-element 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 finite-element polynomial degree . We will present results of such a study in Section 5.3.

###### Remark 12

Recall the ill-posedness 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 ill-conditioning 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 finite-element state and control can be represented as

 uh(x)=n∑i=1(¯uh)iΦi(x),andch(x)=m∑i=1(¯ch)iϕi(x),

respectively, where and are the nodal coefficients. Substituting these expressions for and into the advection-diffusion trilinear form (18) and the discretized objective, the discretized inverse IBM problem can be expressed as the following finite-dimensional optimization problem:

 min¯ch,¯uh (15) s.t. bh(uh,ch,vh)=¯vTh(Au¯uh+Ac¯ch−¯fh)=0,∀¯vh∈Rn.

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 saddle-point 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 special-purpose 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 full-space 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 reduced-space optimization statement:

 min¯chJh(¯ch)=12¯cThH¯ch+gTc¯ch+J0,c

where is a constant, is the reduced gradient and

 H=ATcA−TuHuuA−1uAc−ATcA−TuHTcu−HcuA−1uAc+Hcc (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 Poisson-type 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 non-singular. 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 finite-element 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 model-problem 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 higher-order 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 finite-element 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 advection-diffusion 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 ill-posed 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 ill-posedness is addressed by ensuring the expanded domain tends to the target domain and by including a control-state penalty term in the objective. Applied to a discontinuous Galerkin finite-element discretization of the advection, advection-diffusion, and diffusion equations, the inverse IBM formulation achieves optimal solution convergence rates, and the reduced Hessian remains well-conditioned.

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 finite-element discretizations. Furthermore, the method is compatible with high-order 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 PDE-constrained optimization problem, which requires the solution of a linear saddle-point 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 PDE-constrained optimization problem associated with the inverse IBM. A related issue is the need for specialized preconditioners for the saddle-point 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

 b(u,v)=∫~Ω∇u⋅∇vdΩ=0,∀v∈H10(~Ω).

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

 b(u′,v)=∫~Ω∇u′⋅∇vdΩ=−∫~Ω∇ex(c)⋅∇vdΩ,∀v∈H10(~Ω).

Using the above equation, together with the Poincaré-Friedrichs and Cauchy-Schwarz inequalities, we find

 ∥u′∥21,~Ω≤M1|b(u′,u′)|≤M1∣∣∣∫~Ω∇ex(c)⋅∇u′dΩ∣∣∣≤M1∥ex(c)∥1,~Ω∥u′∥1,~Ω,

where is some constant. Therefore, for we have the following bound on :

 ∥u∥1,~Ω≤∥u′∥1,~Ω+∥ex(c)∥1,~Ω≤(1+M1)∥ex(c)∥1,~Ω≤(1+M1)M2∥c∥12,~Γ,

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

Let be an open subset. The form defined by (3) is bilinear, symmetric, continuous, and -coercive. The form defined by (4) is linear and continuous.

###### 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