Antoine Michael
email: antoine.jost@u-bordeaux.fr

Diego Jost

Stéphane Glockner
email: glockner@bordeaux-inp.fr

Direct forcing immersed boundary methods: Improvements to the Ghost-Cell Method

Keywords: Immersed boundary method, Direct forcing, Ghost-cell method, Discretization stencil, Poisson problem, Incompressible Navier-Stokes

. The quadratic ghost-cell method [2, 3] further improves the accuracy and convergence whilst maintaining a maximum stencil size of 3 while the currently proposed quadratic square shifting method makes it possible to increase the Lagrange polynomial interpolation order whilst maintaining a maximum stencil of 2 resulting in an improved order of convergence for Neumann immersed boundary conditions. The proposed methods are evaluated by considering their increase in performance (CPU time speed-up), accuracy, and convergence thanks to a comprehensive verification and validation process. Firstly, the canonical verification 2D and 3D Poisson test problems for various analytical solutions and immersed boundaries is considered and is followed by various verification and validation test cases for the Navier-Stokes governing equations, with and without heat transfer, in 2D and 3D.

Introduction

Immersed boundary methods (IBMs) were first developed by Peskin [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] and are widely used in numerical simulations involving complex geometries and complex fluid-structure interactions, such as heart valves [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF]. The term IBM refers to any method that solves the discretized governing equations, such as the Navier-Stokes equations, on grids that do not conform to the boundary of the immersed boundary [START_REF] Mittal | Immersed boundary methods[END_REF], i.e. non-body conformal grids. As IBMs can use Cartesian and rectilinear grids, regardless of the complexity of the geometry, some of the advantages associated with these methods are significantly simplified grid generation, sustained grid quality (and complexity), and stationary nondeforming grids for moving geometries (i.e. no-remeshing) [START_REF] Mittal | Immersed boundary methods[END_REF]. These qualities become increasingly more advantageous as the underlying physics of the problem increases.

Depending on the approach followed to implement the immersed boundary conditions (IBCs), the IBMs are divided into either the continuous or the discrete forcing method. The continuous forcing method implements the IBCs through a source term, denoted as forcing function, in the non-discretized governing equations while the discrete forcing method imposes the IBCs through modifications to the discretized governing equations [START_REF] Mittal | Immersed boundary methods[END_REF], similar to the implementation of domain boundary conditions. Picot and Glockner [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF] showed that the stencil sizes of numerical methods for direct ghost-cell IBMs are important as they greatly affect the linear solvers, computational memory requirements, and accuracy of the results. Although the linear shifted method proposed by [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF] significantly reduces the stencil size for rectilinear grids these proposed improvements are not applicable to Cartesian grids, for which the linear shifted method cannot be differentiated from the linear method proposed by Mittal et al. [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF]. This article fills the gap and focuses on improvements to the linear method for Cartesian grids that could potentially be extended to the shifted linear method [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF] for rectilinear grids.

Although the extensively used linear method of Mittal et al. [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] uses one probe point to implement the IBCs, various other IBMs that use two probe points have been proposed. Sandberg and Jones [START_REF] Sandberg | Direct numerical simulations of low reynolds number flow over airfoils with trailing-edge serrations[END_REF] and Luo et al. [START_REF] Luo | A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions[END_REF] used two image points to evaluate the wall-normal derivatives of the fluid variables while Das et al. [START_REF] Das | Direct numerical simulation for flow and heat transfer through random open-cell solid foams: Development of an ibm based cfd model[END_REF][START_REF] Das | A sharp-interface immersed boundary method to simulate convective and conjugate heat transfer through highly complex periodic porous structures[END_REF] used two fluid cell values and two image points to implement the Dirichlet and Neumann IBCs, respectively. Limited information of the density flow field at the immersed boundary led Dhamankar et al. [START_REF] Dhamankar | Implementation of a sharp immersed boundary method in a 3-d multi-block large eddy simulation tool for jet aeroacoustics[END_REF] to use a different number of image points based on the flow variables. Two image points were used for the linear extrapolation of density whilst only one image point was used for the velocity flow field. The use of two image points to evaluate the Dirichlet and Neumann IBCs was pursued by Chi et al. [START_REF] Chi | An improved ghost-cell immersed boundary method for compressible flow simulations[END_REF] and Auguste et al. [START_REF] Auguste | Implementation of an immersed boundary method in the meso-nh v5. 2 model: applications to an idealized urban environment[END_REF]. It should be noted that despite the use of quadratic methods for both Dirichlet and Neumann IBCs, to the knowledge of the authors no comprehensive analysis (stencil size, increase in accuracy, computational performance, etc.) of these methods has been pursued.

To address the shortcomings of the linear method for Cartesian grids and the analysis of the quadratic methods, in this article we propose two IBMs to reduce the stencil requirements for square/cubic cells, a quadratic interpolation method with a reduced stencil size, and we perform a comprehensive analysis of the effect of quadratic interpolations for Dirichlet and Neumann IBCs. This article is organized as follows. Section 2 presents the canonical Poisson verification test problem and numerical discretization methods and is followed by an overview of the linear method [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF]. In Section 3 the currently proposed linear square shift method that decreases the maximum stencil size obtained with the linear method for Cartesian grids is presented and is followed by a discussion on the proposed quadratic square shift method in Section 4, that further increases the accuracy of the linear and linear square shift methods. Numerical results for the canonical 2D and 3D verification test problems for the linear square shift method and the various proposed improvements are subsequently presented in Section 5. Section 6 focuses on the Navier-Stokes equations, with and without the energy governing equation, and presents various 3D validation cases and discusses the results obtained with various IBMs. Lastly, Section 7 presents the conclusions.

In the spirit of reproducibility, the currently proposed methods are implemented and the various verification and validation test cases considered are in the massively parallel incompressible opensource CFD code Notus [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF][START_REF] Lemoine | Moment-of-fluid analytic reconstruction on 2d cartesian grids[END_REF][START_REF] Coquerelle | A fourth-order accurate curvature computation in a level set framework for two-phase flows subjected to surface tension forces[END_REF][START_REF] Notus | Notus computational fluid dynamics[END_REF].

Immersed Boundary Method for the Poisson Problem

In this section the finite-difference direct forcing and numerical discretization schemes are presented. For the sake of clarity, the discussions are focused on the canonical 2D Poisson problem as the extensions to 3D and vector flow fields, such as velocity u, are straightforward.

Poisson Problem

Figure 1(a) shows an immersed boundary Γ that divides the 2D computational domain Ω ⊂ IR 2 into an inner domain Ω i and an outer domain Ω o . As a result, the constant coefficient Poisson problem for a scalar field φ with a source field f consists in

∆φ = f, in Ω i , φ = 0, in Ω o , φ = φ D , on ∂Ω D ∪ Γ D , ∂φ ∂n = φ N , on ∂Ω N ∪ Γ N , (1)
where ∂Ω and Γ are split into the Dirichlet and Neumann parts, and φ D and φ N are the known Dirichlet and Neumann boundary conditions, respectively.

Numerical Methods

As shown in Figure 1(b), the computational domain is discretized, with a finite volume method approach, into a Cartesian grid of size m × n with ∆h ι being the cell size in direction ι and is constant throughout the domain. The domain boundary conditions are implemented through c domain ghost cells at each domain boundary such that the grid size becomes (m + 2c) × (n + 2c).

Although the current article focuses on Cartesian grids, for sake of completeness the proposed methods are presented for the general case of rectilinear grids.

The discretized inner domain is solely composed of inner cells while, as seen from Figure 1(b), the discretized outer domain Ω o is further divided into ghost and outer cells. A ghost cell is a cell X i,j in Ω o such that at-least one of its non-diagonal neighboring cells is an inner cell while an outer cell is a cell in Ω o whose non-diagonal neighbors are in Ω o . For 2D computational domains, a cell X i,j in Ω o is considered a ghost-cell if {X i+1,j ,X i-1,j ,X i,j+1 ,X i,j-1 }∩Ω i =0.

The standard second-order central difference is used to discretize equation (1) over the inner computational domain yielding where Φ i,j =φ(X i,j ), F i,j =f (X i,j ), and ∆h=max(∆h ι ). Taking into account the domain and immersed boundary conditions, equation ( 2) is rewritten in matrix form as

Φ i+1,j -2Φ i,j + Φ i-1,j ∆h 2 x + Φ i,j+1 -2Φ i,j + Φ i,j-1 ∆h 2 y = F i,j + O(∆h 2 ), (2)
Ω i Ω o ∂Ω Γ n (a) (b)
(L + E + G + O)Φ i,j = F i,j + B i,j + R i,j , (3)
where L is the coefficient matrix for the matrix form of equation ( 2) without boundary conditions, E is the coefficient matrix that takes into account the boundary conditions, G and O are the coefficient matrices that take into account the IBCs (ghost and outer-cells, respectively), B i,j contains the boundary condition values, and R i,j is the truncation term for the IBCs (presented in subsequent sections) and the spatial discretization of equation [START_REF] Das | Direct numerical simulation for flow and heat transfer through random open-cell solid foams: Development of an ibm based cfd model[END_REF]. It should be noted that F i,j + B i,j + R i,j =0 and Φ i,j is equal to the trivial solution of zero ∀(i, j) ∈ C o , where C o is the set of outer cell indices.

Removing the truncation term from equation (3) results in the closed linear system evaluated for the numerical solution Φi,j ∀(i, j) ∈ C igo , where C igo is the set of inner, ghost, and outer cell indices.

Linear Method

For an accurate implementation of the linear IBM at a ghost-cell X i,j , the scalar field value at the normal reflection of X i,j about Γ, denoted as the image point x IP , the Euclidean distances ∆l = ||x IP -X i,j || 2 and ∆l h = ∆l/2 = ||x BI -X i,j || 2 , the boundary condition at the boundary intercept x BI (the intersection point of the line segment [X i,j , x IP ] and Γ), and the local surface unit normal n of Γ at x BI are needed. Figure 2 shows representative locations of these supplemental data points relative to each other. The linear method implements the IBCs along the linear segment [X i,j , x IP ] using second-order linear approximations and second-order central difference for Dirichlet and Neumann boundary conditions, respectively. As the image point is a normal reflection of the ghost-cell about the boundary point, the Dirichlet boundary condition follows

X i,j X k,l x IP x BI n ∆l h ∆l Figure 2:
Close up view of supplemental data points and values required for the linear method for a sample ghost-cell X i,j . BI: Boundary intercept; IP : Image point; n:

local surface unit normal; ∆l = ||x IP -X i,j || 2 ; ∆l h = ||X i,j -x BI || 2 = ||x IP -x BI || 2 . φ BI = Φ i,j + φ IP 2 + O(∆l 2 h ), (4)
where φ χ =φ(x χ ) is the exact scalar field value of φ at x χ . The Neumann boundary condition is given by

∂φ ∂n BI = φ IP -Φ i,j ∆l + O(∆l 2 h ), (5)
where (∂φ/∂n) BI = ∂(φ(x BI ))/∂n is the exact derivative of the scalar field at x BI . As shown in Figure 2, the location of the IP may not coincide with X i,j ∀(i, j) ∈ Ω i and therefore the value at x IP is obtained through a p-th order Lagrange interpolation using the scalar field values of the surrounding cells following

φ IP = (k,l)∈I P ρ k,l Φ k,l + O(∆h p ), (6)
where ρ k,l are the interpolation coefficients and I P is the set of interpolation indices that may include the ghost-cell value Φ i,j and excludes any outer-cells. Further discussions on the set of interpolation indices are found in [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF].

Combining equations ( 4) to [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF] results in the following discretized numerical approximations of the boundary conditions

φ BI = 1 2 (k,l)∈I IP ρ k,l Φ k,l + Φ i,j + O(∆l 2 h ) + O(∆h p ), ∂φ ∂n BI = 1 ∆l (k,l)∈I IP ρ k,l Φ k,l -Φ i,j + O(∆l 2 h ) + O(∆h p-1 ). ( 7
)
One of the major drawbacks of the linear IBM is the need for a large discretization stencil, even for cell size ratios a=∆h/∆µ=1, where ∆µ = min(∆h ι ). The stencil size measurement metrics presented in [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF] are also used in this article: Figure 3 shows the domains for the stencil sizes of 1 and 2. Picot and Glockner [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF] showed that for the linear method for a cell size ratio a=1 (Cartesian grids) and a second-order Lagrange interpolation (p=2) the maximum stencil size is 2, as shown in Figure 4, which is larger than the stencil size for second-order central difference of the Laplacian operator. As IBMs rely heavily on Cartesian grids with adaptive mesh refinement (AMR) to resolve flow features with large gradients, such as boundary layers, and geometrical features with large curvatures, such as the leading and trailing edges of wings, reducing the stencil size of the method to 1 may improve the precision and reduce the computational effort.

X i,j X i+2,j X i-2,j X i,j-2 X i,j+2
Figure 3: Interpolation quadrants (colored regions) and domains around X i,j of stencil size 1 and 2. Dashed blue line: stencil size of 1; Dashed orange line: stencil size of 2; Red region: Q1 interpolation quadrant; Blue region: Q2 interpolation quadrant; Orange region: Q3 interpolation quadrant; Green region: Q4 interpolation quadrant.

Linear Square Shifting Method

Despite the stencil size improvements of the shifted linear method for rectilinear grids [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF], this method reverts back to the linear method for Cartesian grids and therefore does not address the reduction of the stencil size for square/cubic cells. A stencil size larger than 1 does not yield band matrices with 9 or 27 points in 2D or 3D, respectively, and as a result the more efficient geometric multi-grid algorithms, such as the SMG [START_REF] Brown | Semicoarsening multigrid on distributed memory machines[END_REF] and PFMG [START_REF] Ashby | A parallel multigrid preconditioned conjugate gradient algorithm for groundwater flow simulations[END_REF] in the hypre library [START_REF] Falgout | hypre: A library of high performance preconditioners[END_REF][START_REF] Falgout | The design and implementation of hypre, a library of parallel high performance preconditioners[END_REF], cannot be used. Instead the more generic, and less efficient [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF][START_REF] Yang | Boomeramg: a parallel algebraic multigrid solver and preconditioner[END_REF][START_REF] Langer | Comparison of geometrical and algebraic multigrid preconditioners for data-sparse boundary element matrices[END_REF][START_REF] Sterck | Reducing complexity in parallel algebraic multigrid preconditioners[END_REF][START_REF] Baker | Scaling hypres multigrid solvers to 100,000 cores[END_REF][START_REF] Müller | Massively parallel solvers for elliptic partial differential equations in numerical weather and climate prediction[END_REF], preconditioned algebraic multi-grid algorithms, such as BoomerAMG [START_REF] Yang | Boomeramg: a parallel algebraic multigrid solver and preconditioner[END_REF] in the hypre library, must be used.

The currently proposed square/cubic shifting method, denoted as linear square shift, with its several variations presented in this section focus on decreasing the maximum stencil sizes for Cartesian grids through various methods. Modifications to the IBM result in a set of supplemental interpolated data points needed for the proper implementation of the IBCs. For sake of clarity the various modified supplemental data points and values required for the linear square shift method are shown in Figure 5 and are defined, in a staggered manner, in the subsequent sections.

Quadrant Intersection Point

As shown in Figure 5, for a ghost-cell with a stencil size of 2 the intersection of the line segment [X i,j , x IP ] and the domain of stencil size of 1 corresponds to the quadrant intersection point (QIP)

x QIP . In fact x QIP is the upper threshold for a maximum stencil size of 1 and therefore any point located at x QIP + α n, for α >0, is at a stencil size larger than 1. Consider the boundary intercept x BI , its local surface unit normal n, and the four quadrants used to determine the set of interpolation indices for the boundary intercept I BI [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF]. The four quadrants, shown in Figure 3, are determined by n=(n x , n y ) and are

• Quadrant 1 (Q1): n x ≥0 & n y ≥0, • Quadrant 2 (Q2): n x <0 & n y ≥0, • Quadrant 3 (Q3): n x <0 & n y <0, • Quadrant 4 (Q4): n x ≥0 & n y <0.
As the relative location of x QIP , and not its absolute location, is necessary to calculate ||x BIx QIP || 2 for the linear square shift method, Q2 to Q4 may be subsequently transformed into Q1. This change in coordinate direction simplifies the algorithm used to determine the x QIP . As seen from Figure 6, in Q1 and for a local coordinate axis centered at x BI , x QIP must lay on either the line segment Two different approaches are currently pursued to maintain a maximum stencil size of 1, such that

[X k,l , X k,l+1 ] or the line segment [X i,l+1 , X k,l+1 ]. X k,l x BI X k,l+1 x y X i,l+1 X i,j n x QIP x QIP
||X i,j -x IP || 2 ≤ ||X i,j -x QIP || 2 . ( 8
)
The first approach, denoted as ghost-cell shift (GS) shifts the ghost-cell such that the mirror image point satisfies equation [START_REF] Luo | A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions[END_REF], while the second approach satisfies equation (8) but solely shifting the image point and is called image point shift (IS).

Ghost-Cell Shift (GS)

Considering the ghost-cell X i,j , defined in Section 2.2, and the quadrant intersection point x QIP , a supplemental interpolated ghost point x GP is defined as

x GP = x BI -s n, (9)
where s=min{||x BI -X i,j || 2 , α||x BI -x QIP || 2 }, ||x BI -X i,j || 2 = ∆l h as x BI is the midpoint of the line segment [X i,j , x IP ], and α is a constant with a value between 0 and 1, inclusive. In this article we choose α=1.

As seen from Figure 5, equation ( 9) shifts the ghost point x GP towards its boundary intercept x BI such that its image point x IP is located inside Q1 with a maximum stencil size of 1. As the ghost-cell shift of the linear square shift method maintains x BI as the midpoint of the line segment [x GP , x IP ], the general form of the discretized numerical approximations of the IBCs in equation ( 7) remains the same. However, as x GP is not a cell value, an additional p-th order Lagrange interpolation is required and results in a new set of IBCs that follow

φ BI = 1 2 (k,l)∈I IP β k,l Φ k,l + (k,l)∈I GP γ k,l Φ k,l + O(∆l 2 h ) + O(∆h p ), ∂φ ∂n BI = 1 ∆l (k,l)∈I IP β k,l Φ k,l - (k,l)∈I GP γ k,l Φ k,l + O(∆l 2 h ) + O(∆h p-1 ), (10)
where ∆l h =∆l /2, γ k,l , and β k,l are the interpolation coefficients for the x GP and x IP , respectively, and I GP and I IP are the set of interpolation indices for x GP and x IP , respectively. It should be noted that based on the definition of x GP in equation ( 9), a maximum stencil size of 1 is guaranteed for p=2.

Remark It is interesting to notice that as ∆l → 0 (α → 0), the ghost-cell shift of the linear square shift method becomes the IBM proposed by Coco and Russo [START_REF] Coco | Finite-difference ghost-point multigrid methods on cartesian grids for elliptic problems in arbitrary domains[END_REF], hereby denoted as direct method. Indeed, as at ∆l =0, the interpolation coefficients of the image and ghost points are the same and the Dirichlet IBC becomes

φ BI = (k,l)∈I BI ρ k,l Φ k,l + O(∆h p ), (11)
which corresponds to the direct method.

For the Neumann IBC, as ∆l → 0 the error of the second-order central difference becomes increasingly negligible and the finite-difference approximation approaches the analytical partial differential equation ∂φ/∂n = φ n = ∇φ • n. Therefore, the Neumann IBC becomes

∂φ ∂n BI = (k,l)∈I BI ∂ρ k,l ∂x n x + ∂ρ k,l ∂y n y Φ k,l + O(∆h p-1 ), (12)
which corresponds to the direct method.

Image Point Shift (IS)

Reductions to the maximum stencil size can also be achieved by solely shifting the location of the image point such that it is located inside Q1 for a stencil size of 1. Considering the image point x IP , defined in Section 2.3, and the quadrant intersection point x QIP a modified image point x IP is defined as

x IP = x BI + s n, (13)
where s =min{||x BI -x IP || 2 , α||x BI -x QIP || 2 } and α is a constant with a value between 0 and 1, inclusive. In this article we choose α=1. As seen from Figure 5, equation ( 9) shifts the image point x IP towards its boundary intercept x BI such that it is located inside Q1 with a maximum stencil size of 1 for p=2. This method differs from the ghost-cell shift method in that the ghost-cell X i,j is used and no supplemental ghost point x GP is required, and instead only the image point x IP is shifted to the new image point x IP if it falls outside Q1.

In this method the boundary intercept x BI is no longer the midpoint of the line segment [X i,j , x IP ] and to take in account that ∆l 1 = ∆l 2 , where ∆l 1 =||X i,j -x BI || 2 and ∆l 2 =||x BI -x IP || 2 , linear interpolation is used to generalize equations ( 4) and ( 5) such that the generalized second-order Dirichlet boundary condition follows,

φ BI = ζΦ i,j + ηφ IP + O(∆l 1 ∆l 2 ), ( 14
)
where ζ=∆l 2 /(∆l 1 +∆l 2 ) and η=∆l 1 /(∆l 1 +∆l 2 ). The generalized first-order Neumann boundary condition follows,

∂φ ∂n BI = φ IP -Φ i,j ∆l 1 + ∆l 2 + O(∆l 1 |(1 -r ∆l )|) + O(∆l 3 ), (15)
where r ∆l =∆l 2 /∆l 1 and has a values between 0 and 1, inclusive, and ∆l 3 = ∆l 1 2 |(r 2 ∆l -r ∆l + 1)|. Considering the generalized Dirichlet and Neumann boundary conditions, equations ( 14) and [START_REF] Brown | Semicoarsening multigrid on distributed memory machines[END_REF], respectively, the discretized numerical approximations of the IBCs for the image point shift of the linear square shift method follow

φ BI = ζΦ i,j + (k,l)∈I IP ηβ k,l Φ k,l + O(∆l 1 ∆l 2 ) + O(∆h p ), ∂φ ∂n BI = 1 ∆l 1 + ∆l 2 (k,l)∈I IP β k,l Φ k,l -Φ i,j + O(∆l 1 |(1 -r ∆l )|) + O(∆l 3 ) + O(∆h p-1
). ( 16) Similar to the ghost-cell shift of the linear square shift, based on the definition of the x IP in equation ( 13), a maximum stencil size of 1 is guaranteed for p=2.

Quadratic Square Shifting Method

In this section accuracy improvements to the image and ghost square shift methods through quadratic interpolation, denoted as the quadratic square shifting method, are presented. Although quadratic interpolation to increase the accuracy of the ghost-cell interpolation has been previously implemented for Dirichlet and Neumann IBCs [START_REF] Das | Direct numerical simulation for flow and heat transfer through random open-cell solid foams: Development of an ibm based cfd model[END_REF][START_REF] Das | A sharp-interface immersed boundary method to simulate convective and conjugate heat transfer through highly complex periodic porous structures[END_REF], as discussed in Section 1, a detailed analysis of these benefits and its use in conjunction with the square shift methods have not yet been pursued. Figure 7, shows representative locations for the supplemental data points, including the second image point x IP 2 , for the quadratic/quadratic square shifting method.

x BI Considering the boundary intercept x BI , its local surface unit normal n, and the image point x IP , the second image point x IP 2 required for the current quadratic interpolation follows

X i,j x IP x IP 2 X k,l n ∆l h
x IP 2 = x BI + 3∆l h n. ( 17
)
The generalized Dirichlet, equation [START_REF] Notus | Notus computational fluid dynamics[END_REF], is extended to the quadratic method and becomes

φ BI = λφ IP 2 + ηφ IP + ζΦ i,j + O(∆l 3 h ) + O(∆h p ), (18)
where ζ, η, and λ are the coefficients obtained from the 1D p-th order Lagrange interpolation.

The generalized Neumann boundary condition, equation [START_REF] Brown | Semicoarsening multigrid on distributed memory machines[END_REF], is extended to the quadratic method through the use of second-order one-sided finite-difference and 1D third-order Lagrange interpolation (p=3) [START_REF] Das | A sharp-interface immersed boundary method to simulate convective and conjugate heat transfer through highly complex periodic porous structures[END_REF]. As the Neumann IBC at x BI are known, the second-order forward finite-difference, which follows 19) is used to determine the value of the scalar field value at x BI , which is subsequently used to determine the value of the scalar field at the ghost-cell through equation [START_REF] Falgout | The design and implementation of hypre, a library of parallel high performance preconditioners[END_REF]. The discretized numerical approximations of the Dirichlet and Neumann IBCs therefore become

∂φ ∂n BI = -φ IP 2 (x IP -x BI ) 2 + φ IP (x IP 2 -x BI ) 2 -φ BI [(x IP 2 -x BI ) 2 -(x IP -x BI ) 2 ] (x IP -x BI )(x IP 2 -x BI )(x IP 2 -x IP ) +O(∆l 2 h ), (
φ BI = (k,l)∈I IP 2 λρ k,l Φ k,l + (k,l)∈I IP ηβ k,l Φ k,l + ζΦ i,j + O(∆l 3 h ) + O(∆h p ), ∂φ ∂n BI = σ -[(x IP 2 -x BI ) 2 -(x IP -x BI ) 2 ]φ BI - (x IP -x BI ) 2 (k,l)∈I IP 2 ρ k,l Φ k,l + (x IP 2 -x BI ) 2 (k,l)∈I IP β k,l Φ k,l + O(∆l 2 h ) + O(∆h p-1 ), ( 20
) where σ=1/(x IP -x BI )(x IP 2 -x BI )(x IP 2 -x IP ).
As shown in Figure 8, a ghost-cell X i,j with a stencil size of 2 with the linear method may need a stencil size of 3 for the quadratic method. To maintain a maximum stencil size of 2, the quadratic square shifting method is proposed. This method uses the quadratic method in conjunction with the image point shift (solely for the second image point) such that

x IP 2 = x BI + s * n, (21)
where

s * =min{||x BI -x IP 2 || 2 , α||x BI -x QIP 2 || 2 }, x QIP 2
is the QIP for a stencil size of 2 (shown as the dashed red region in Figure 8), and α is a constant with a value between 0 and 1, inclusive.

In this article we choose α=1.

Remark The quadratic square shifting method does not always result in a stencil size larger than 1. It is possible that the Euclidean distance between the ghost-cell and boundary-intercept ∆l h is small enough such that ||x IP 2 -X i,j || 2 ≤ ||x QIP -X i,j || 2 resulting in identical interpolation indices for the first and second image points. The set of interpolation indices for the first and second image points may differ [START_REF] Das | Direct numerical simulation for flow and heat transfer through random open-cell solid foams: Development of an ibm based cfd model[END_REF][START_REF] Das | A sharp-interface immersed boundary method to simulate convective and conjugate heat transfer through highly complex periodic porous structures[END_REF][START_REF] Luo | A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions[END_REF], potentially increasing the accuracy of the solution [START_REF] Dhamankar | Implementation of a sharp immersed boundary method in a 3-d multi-block large eddy simulation tool for jet aeroacoustics[END_REF]. This is accomplished by again re-purposing the image point shift and using the following algorithm:

1. ||x IP -X i,j+1 || c ≤ 1 & 1 < ||x IP 2 -X i,j+1 || c ≤ 2; return 2. ||x IP -X i,j+1 || c ≤ 1 & ||x IP 2 -X i,j+1 || c ≤ 1 • x IP 2 = x BI + s * * n, s * * =(2 -α) × ||x BI -x QIP 1 || 2 ; α > 0 and α < 1 3. 1 < ||x IP -X i,j+1 || c ≤ 2 & 1 < ||x IP 2 -X i,j+1 || c ≤ 2 • x IP follows equation (13) • x IP 2 follows equation (21)
where the stencil norm

x c =max |x| ∆x1 , |y| ∆x2 is used to determine the stencil as c i,j = max{ X k,l - X i,j c | ∀(k, l) ∈ C ig , L k,l i,j = 0} [6].
As the definition of the second image point is based on the location of the boundary intercept, local surface unit normal, and first image point, the quadratic square shifting method is based on a straightforward integration of the quadratic and square shifting methods. To this end the appropriate boundary intercepts, first image points, and local surface unit normals associated with

X k,l X i,j x BI x IP x IP 2 n x QIP 2 x * IP 2 x QIP 1 X i,j+1 X k,l+1 x IP x * IP 2 x BI x IP 2 X k+1,l+1 X k+1,l+2
X k+1,l n Figure 8: Close up view of the supplemental data points and values required for the quadratic square shift method for the sample ghost-cells X i,j and X i,j+1 ; x BI : Boundary intercept; x IP : Image point; x IP 2 : Second image point; x * IP 2 : Shifted second image point; x QIP Λ : quadrant intersection point for stencil size of Λ; n: local surface unit normal; Dashed red region: stencil size of 2 for X i,j ; Dashed green region: stencil size of 2 for X i,j+1 ; Solid orange line: Γ. For plotting purposes for the ghost-cell

X i,j the distance ||x BI -x IP 2 || 2 < 3∆l h .
the ghost point/cell in the square shift method must be used in equations ( 17) to [START_REF] Sterck | Reducing complexity in parallel algebraic multigrid preconditioners[END_REF]. Henceforth, and without any loss in comprehension, the quadratic square shifting method will be denoted as the quadratic method.

The quadratic method can be used in conjunction with the square shift methods and results in the QGS, QIS, and QISS1 methods. The QGS and QIS are obtained by combining the GS and IS methods, respectively, with the quadratic method. For a second-order Lagrange interpolation (p=2) these two methods have a maximum stencil size of 2 with a stencil size of 1 for the first image point. The QISS1 method is based on the QIS method and also guarantees a maximum stencil size of 1 for the second image point. Table 1 summarizes the proposed quadratic ghost-cell methods.

Method Maximum Stencil Size Q 2 QGS 2 QGSO 2 QIS 2 QISS1 1
Table 1: Maximum stencil sizes for the proposed quadratic ghost-cell methods with a second-order Lagrange interpolation (p=2) for the image points.

Numerical Simulations for a Scalar Field

In this section the results of numerical simulations for the canonical 2D and 3D Poisson test problems with mixed domain boundary conditions and Dirichlet or Neumann IBCs for various immersed boundaries and IBMs are discussed. The nine IBMs used for the current verification and convergence studies are:

• L: linear method [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] (stencil size of 2),

• LGS: ghost-cell linear square shift method (stencil size of 1),

• LIS: image point linear square shift method (stencil size of 1),

• Q: quadratic ghost point(stencil size of 2),

• QX: quadratic ghost point with X ∈ {GS, IS} (stencil size of 2),

• QY O: quadratic ghost point with Y ∈ {GS} with different interpolation indices for x IP and x IP 2 (stencil size of 2),

• QISS1: quadratic ghost point with image point shift for a maximum stencil size of 1 (stencil size of 1).

For sake of clarity and compactness, the results for select representative methods are presented and discussed with appropriate comments on the results of the remaining methods. The set of representative methods depends on the analytical solution of the Poisson problem and the immersed boundary.

As the analytical solution to equation ( 3) with appropriate boundary conditions and source field can be determined, the L ∞ and L 2 norms of the error fields are evaluated and follow,

L ∞ (Φ e ) = max(|Φ e |) ∀(i, j) ∈ C i , L 2 (Φ e ) = i,j |Φ e | 2 V -i,j ∀(i, j) ∈ C i , (22)
where Φ e = Φi,j -Φ i,j is the error field and Vi,j is the inner volume of an inner cell, i.e. only the part of an inner cell that is located inside Ω i . The numerical simulations have a computational domain of Ω=[-1,1] ς and Ω=[0,1] ς , where ς is the number of spatial dimensions, for numerical simulations with parabolic and trigonometric analytical solutions, respectively. Dirichlet boundary conditions are used for the x-axis domain boundaries and Neumann boundary conditions are used for the remaining domain boundaries. Figure 9 shows the immersed boundaries Γs are the canonical circle and sphere with a radius r=0.325L x , where L x is the length of the x-axis domain edge, and the more complex, with convex and concave corners, flower-shaped geometries, denoted as flower, that follow

• 2D [6, 25]: (0.02 √ 5,0.02 √ 5)+(0.5+0.2sin(5θ))e r ∀θ ∈ [0, 2π],
• 3D:

(x 2 + y 2 + z 2 = r 2 ) - ((x ± r) 2 + y 2 + z 2 = r 2 f ) -(x 2 + (y ± r) 2 + z 2 = r 2 f ) -(x 2 + y 2 + (z ± r) 2 = r 2 f ),
where r f =r/2 and the centers of the immersed boundaries are at the centers of the computational domains. The Cartesian grid sizes considered are (16ν) ς , where ν=2 n-1 and n is the grid number in the grid convergence studies (i.e. first grid: n = 1, second grid: n = 2, etc.).

Remark

The second-order linear interpolations of the LGS and LIS methods for Dirichlet IBCs and first-order central difference of the LIS method for Neumann IBCs limit the use of third-order Lagrange interpolations (p=3) to the LGS method for Neumann IBCs. The stencil size requirements of the quadratic methods also limit the use of p=3 to the QISS1 method. As a result, unless stated otherwise second-order Lagrange interpolations (p=2) are used to determine the value at the various image and ghost points.

(-1, -1)

(1, 1) r (a) (-1, -1)
(1, 1)

(b) (c) (d)

Dirichlet IBCs

Dirichlet IBCs for parabolic and trigonometric analytical solutions are taken into account for the current series of verification and convergence studies. Calculating and keeping the appropriate values of R i,j in equation ( 3) to remove all non-IBM truncation errors, denoted as numerical correction, is necessary for a complete analysis of the accuracy and convergence of IBMs [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF], and will therefore be pursued. For sake of clarity, the results for select representative immersed boundaries for either a parabolic or a trigonometric analytic solution are presented and discussed in this section.

Parabolic Solution

One of the analytical solutions to the 3D Poisson problem follows

φ(x) = (1 + x) 2 , f (x) = 2, ∀(x) ∈ Ω i ∈ IR 3 , (23)
given that

φ(x) = (1 + x) 2 ∀(x) ∈ Γ, φ(-1, y, z) = 0, φ(1, y, z) = 4 ∀(y, z) ∈ [-1, 1], ∂φ ∂n (x, -1, z) = 0, ∂φ ∂n (x, 1, z) = 0 ∀(x, z) ∈ [-1, 1], ∂φ ∂n (x, y, -1) = 0, ∂φ ∂n (x, y, 1) = 0 ∀(x, y) ∈ [-1, 1], (24)
where numerical corrections are applied to the x-axis domain boundary conditions [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF]. The analytical solution for the corresponding 2D Poisson problem follows equation ( 23) with the appropriate 2D source term (e.g. f (x, y) = 2) and boundary conditions.

Trigonometric Solution

A second analytical solution to the Poisson problem is

φ(x) = sin(πx)cos(2πy), f (x) = -5π 2 sin(πx)cos(2πy) ∀(x) ∈ Ω i ∈ IR 2 , φ(x) = sin(πx)cos(2πy)cos(2πz), f (x) = -9π 2 sin(πx)cos(2πy)cos(2πz) ∀(x) ∈ Ω i ∈ IR 3 , (25) given that φ(x) = sin(πx)cos(2πy) ∀(x) ∈ Γ ∈ IR 2 , φ(0, y) = 0, φ(1, y) = 0 ∀(y) ∈ [0, 1], ∂φ ∂n (x, 0) = 0, ∂φ ∂n (x, 1) = 0 ∀(x) ∈ [0, 1], φ(x) = sin(πx)cos(2πy)cos(2πz) ∀(x) ∈ Γ ∈ IR 3 , φ(0, y, z) = 0, φ(1, y, z) = 0 ∀(y, z) ∈ [0, 1], ∂φ ∂n (x, 0, z) = 0, ∂φ ∂n (x, 1, z) = 0 ∀(x, z) ∈ [0, 1], ∂φ ∂n (x, y, 0) = 0, ∂φ ∂n (x, y, 1) = 0 ∀(x, y) ∈ [0, 1], (26)
where numerical corrections are not needed for any of the domain boundary conditions.

Reduce Stencil Size to 1

Figure 10 shows the L ∞ and L 2 norms of Φ e for various immersed boundaries for the Poisson problem with parabolic and trigonometric analytical solutions for the LGS and LIS methods and, as expected, they have a second-order limiting behavior. As seen from Figure 10(a), the LGS and LIS methods are more accurate than the L method, yield a more consistent second-order limiting behavior (especially for L ∞ ), and may be used interchangeably (i.e. negligible differences in the accuracy of the results). This similarity is attributed to the similar impact of the interpolation truncation errors of the ghost point in the LGS method and Dirichlet IBCs for the LIS method. For the LIS method X i,j = x GP and, therefore, uses one Lagrange interpolation less than the LGS method, but as ∆l 1 + ∆l 2 > ∆l the interpolation error associated with the Dirichlet IBCs is larger for the LIS method than the LGS method. As seen from Figure 10(b), for the sphere immersed boundaries the improvements in accuracy in the L ∞ and L 2 norms previously seen are increased and reduced, respectively. This behavior is not unique to the sphere immersed boundaries and is characteristic to all comparisons between 2D and 3D verification cases as it is also observed in the comparison of the error norms between the 2D and 3D flower immersed boundaries for a parabolic analytical solution (not shown). Despite this slight decrease in accuracy for the L 2 norm for all 3D immersed boundaries currently considered, the LGS and LIS methods still yield results with significantly smaller maximum errors and have a maximum stencil size of 1.

Figure 10(c) shows that for a 3D flower immersed boundary with a trigonometric analytical solution the LGS and LIS methods only slightly improve the accuracy of the L ∞ norm, yield a L ∞ norm with a more consistent second-order limiting behavior, and have little to no impact on the L 2 norm. This behavior is not unique to the 3D flower immersed boundary and is present for all comparisons between parabolic and trigonometric analytical solutions. It should be noted that despite the limited accuracy improvements of the LGS and LIS methods for the 3D flower immersed boundary for a trigonometric analytical solution, these methods have a maximum stencil size of 1.

Conclusion: Although for Dirichlet IBCs the LGS and LIS methods result in, at best, significantly smaller error norms, and, at worst, negligible improvements these methods tend to yield error norms with a more consistent second-order limiting behavior, compared to the L method, and have a maximum stencil size of 1, whereas the L method has a maximum stencil size of 2. downwards by factor of 10 for ease of reading.

Improved Accuracy with the Quadratic Method

In this section the improvements to the accuracy of the numerical solutions obtained with the proposed square shifting methods for various immersed boundaries for Poisson problems with parabolic and trigonometric analytical solutions are presented.

QX Methods Figure 11 shows the L ∞ and L 2 norms of Φ e for various immersed boundaries for a Poisson problem with parabolic and trigonometric analytical solutions for the Q, QGS, and QIS methods (denoted as the QX methods) and, as expected, have a second-order limiting behavior. As seen from Figure 11(a), for the same stencil size of 2 the QX methods are more accurate than the L method, and similar to the errors norms for the LGS and LIS methods, the QGS and QIS methods yield more consistent second-order limiting behaviors (especially for L ∞ ) and may be used interchangeably. As seen by comparing the QX methods in Figures 11(a) and 11(b), the Q method is least adept at smoothing out the less consistent second-order convergence of the L method, yields the largest errors norms, and cannot be used interchangeably with the QGS or QIS methods.

The trends observed in Figures 11(a) and 11(b) are also observed for all immersed boundaries and for the parabolic and trigonometric analytical solutions except for the 3D flower immersed boundary for a parabolic analytical solution, shown in Figure 11(c). For this verification case the QX methods, which have the same stencil size of 2 as the L method, are more accurate than the L method, yield a slightly more consistent second-order behavior, and may be used interchangeable. norms of Φ e for 2D flower and sphere immersed boundaries, respectively, for a Poisson problem with trigonometric analytical solution for the QIS and QISS1 methods. Although the QISS1 method with a maximum stencil size of 1 is more accurate and yields a more consistent second-order limiting behavior compared to the L method, it is marginally less accurate than the QIS method with a maximum stencil size of 2. It should be noted that the QISS1 method is a viable alternative to the QIS method as the slightly larger error norms of the former, compared to the latter, is accompanied by significant performance gains associated with its stencil size of 1.

Different Image Points Interpolation Sets

Third-Order Lagrange Interpolation

As third-order Lagrange (p=3) interpolations extend the stencil size of second-order Lagrange interpolations from 1 to 2 their use with IBMs with a stencil size of 2 is strongly discouraged, i.e. they extend the stencil size from 2 to 3. To prevent these additional performance penalties related to the stencil size of 3, p=3 interpolations should only be pursued with IBMs that have a maximum stencil size of 1, such as the linear square shift method and QISS1. Although the use of p=3 with these IBMs extends their maximum stencil size to 2 (the same as the L method), thereby rending their performance gains null, their increased accuracy still characterizes them as improvements of the L method. Figure 13 shows the L ∞ and L 2 norms of Φ e for a circle and sphere immersed boundary for a Poisson problem with a parabolic and trigonometric analytical solution for the LGS method with p=2 and p=3. As expected, the LGS method with third-order Lagrange interpolations also results in second-order convergence, as the truncation error on the order of max(O(∆l 2 h ), O(∆h 3 )), for both errors norms and yields more accurate results than the LGS method with p=2. Although the L ∞ and L 2 norms have a second-order convergence, the numerical error field for the parabolic analytical solution, shown in Figure 12, is not on the order of machine epsilon due to the linear interpolation used at the domain boundary conditions and the second-order Lagrange interpolations (p=2) used for the ghost and image points. However, for the Poisson problem with the parabolic analytical solution, the combined used of the QISS1 method, third-order Lagrange interpolations (p=3), and extending the quadratic method with second-order Lagrange interpolations (p=2) to the domain boundaries yields error fields with values on the order of machine epsilon. By also considering the Poisson problem with the trigonometric analytical solution, shown in Figure 14, it can be concluded the QISS1 method and p=3 for equation ( 20) not only significantly increase the accuracy regardless of the Poisson analytical solution, but also improves the order of convergence of the IBM to third-order (appropriate numerical corrections are used to isolate the IBM truncation errors). This behavior is expected as equation ( 20) now has a truncation error on the order of max(O(∆l 3 h ), O(∆h 3 )).

Conclusion

For Dirichlet IBCs and a maximum stencil size of 1, the proposed shifted linear methods (LGS and LIS) result in a numerical solution with the smallest L ∞ and L 2 norms for all immersed boundaries and analytical solutions currently considered, except for sphere immersed boundary with trigonometric analytical solution wherein the QISS1 method is slightly more accurate. These methods are also more accurate than the standard linear method and, as they have a smaller stencil size, yield more computationally efficient band matrices.

For Dirichlet IBCs and a maximum stencil size of 2, the proposed quadratic shifted methods (QGS and QIS) result in a numerical solution with the smallest L ∞ and L 2 norms for all current verification cases. Given that the standard linear and QX methods have the same maximum stencil size of 2, the latter methods increase the accuracy of the solution without any significant performance penalties. It should be noted that methods that use different interpolation sets (QY O) have negligible differences with the QX methods. As expected the QGS and QIS methods are more accurate than the LGS and LIS methods for all verification cases currently considered. By increasing the maximum stencil size of the QISS1 method from 1 to 2, through the use of third-order Lagrange interpolations (p=3), a more accurate numerical solution with an improved order of convergence of the L ∞ and L 2 norms of three is obtained without any significant performance penalty.

Table 2 summarizes the maximum stencil sizes and order of convergences of the various proposed methods for Dirichlet IBCs for parabolic and trigonometric analytical solutions for the 2D and 3D Poisson problems.

Method Maximum Stencil Size Order of Convergence

One Image Point L (p=2) 2 2 LGS (p=2) 1 2 LGS (p=3) 2 2 LIS (p=2) 1 2 Two Image Points Q (p=2) 2 2 QGS (p=2) 2 2 QGSO (p=2) 2 2 QIS (p=2) 2 2 QISS1 (p=2) 1 2 QISS1 (p=3) 2 3
Table 2: Maximum stencil sizes and order of convergences for the proposed methods for Dirichlet IBCs.

Neumann IBCs

Neumann IBCs for parabolic and trigonometric analytical solutions are also taken into account for the current series of verification and convergence studies. Similar to Dirichlet IBCs, numerical correction is taken into account to emphasise the accuracy of the IBMs currently considered.

Parabolic Solution

The analytical solution of the current 3D Poisson problem follows

φ(x) = (1 + x) 2 , f (x) = 2, ∀(x) ∈ Ω i , (27)
given that

∂φ ∂n (x) = 2x(1 + x) x 2 + y 2 + z 2 ∀(x) ∈ Γ, φ(-1, y, z) = 0, φ(1, y, z) = 4 ∀(y, z) ∈ [-1, 1], ∂φ ∂n (x, -1, z) = 0, ∂φ ∂n (x, 1, z) = 0 ∀(x, z) ∈ [-1, 1], ∂φ ∂n (x, y, -1) = 0, ∂φ ∂n (x, y, 1) = 0 ∀(x, y) ∈ [-1, 1], (28)
where numerical corrections only need to be applied to the x-axis domain boundary conditions [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF].

The analytical solution for the corresponding 2D Poisson problem follows equation ( 28) with the appropriate 2D source term and boundary conditions (e.g. ∂φ(x)/∂n = 2x(1 + x)/ x 2 + y 2 ).

Trigonometric Solution

A second analytical solution to the Poisson problem is

φ(x) = sin(πx)cos(2πy), f (x) = -5π 2 sin(πx)cos(2πy) ∀(x) ∈ Ω i ∈ IR 2 , φ(x) = sin(πx)cos(2πy)cos(2πz), f (x) = -9π 2 sin(πx)cos(2πy)cos(2πz) ∀(x) ∈ Ω i ∈ IR 3 , (29)
2πsin(πx)cos(2πy)sin(2πz)n z ∀(x) ∈ Γ ∈ IR 3 , φ(0, y, z) =0, φ(1, y, z) = 0 ∀(y, z) ∈ [0, 1], ∂φ ∂n (x, 0, z) =0, ∂φ ∂n (x, 1, z) = 0 ∀(x, z) ∈ [0, 1], ∂φ ∂n (x, y, 0) =0, ∂φ ∂n (x, y, 1) = 0 ∀(x, y) ∈ [0, 1], (30)
where n x , n y , and n z are the local surface unit normals in the x, y, and z-direction, respectively, and numerical corrections are not needed for any of the domain boundary conditions.

For an appropriate comparison with the results obtained with the first-order LIS method for Neumann IBCs, second-order Lagrange interpolations (p=2) for the ghost and image points are considered. As seen from equation [START_REF] Sandberg | Direct numerical simulations of low reynolds number flow over airfoils with trailing-edge serrations[END_REF], for Neumann IBCs a first-order limiting behavior is expected for second-order Lagrange interpolations.

Reduce Stencil Size to 1

Figure 15 shows the L ∞ and L 2 norms of Φ e for a circle and sphere immersed boundary for a Poisson problem with the parabolic analytical solution for the LGS and LIS methods, and, as expected, have a first-order limiting behavior. Similar to the results for Dirichlet IBCs, the LGS and LIS methods may be used interchangeably (i.e. negligible differences in the accuracy of the results). The LGS and LIS methods also yield the same L ∞ norms and slightly smaller L 2 norms compared to the L method. It should be noted that although the increase in accuracy obtained with the linear square shift methods for the parabolic analytical solution with Neumann IBCs is smaller than those obtained with their Dirichlet IBCs counterpart, the maximum stencil size of 1 of these methods lead to better performances compared to the maximum stencil size of 2 of the L method. norms are shifted downwards by factor of 10 for ease of reading.

Figure 16 shows that for a circle and a sphere immersed boundary for a Poisson problem with a trigonometric analytical solution the L ∞ and L 2 norms of Φ e for the LGS and LIS methods have between first and second-order limiting behavior. Unlike the behavior for a parabolic analytical solution with Neumann IBCs, for a trigonometric analytical solution the LGS and LIS methods are more accurate than the L method, increase the order of convergence from 1 to approximately 1.5 for a circle immersed boundary, and yield a more consistent first-order limiting behavior for a sphere immersed boundary (especially the LGS method). Although the LGS and LIS methods may only be used interchangeably for a circle immersed boundary, the differences in the errors norms between the two methods for a sphere immersed boundary are limited. The superconvergence seen in Figure 16(a) is not yet explained and it is therefore more prudent to retain an overall order of convergence of 1, as seen in the other test cases and in accordance with the truncation errors in equations ( 12) and ( 16). The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Improved Accuracy with the Quadratic Method

In this section the improvements to the accuracy of the numerical solutions obtained with the proposed square shifting methods for various immersed boundaries for Poisson problems with parabolic and trigonometric analytical solutions are presented.

QX Methods Similar to the linear square shift methods for Neumann IBCs, the QX methods yield the same L ∞ norm as the L method and only the QGS method yields slightly smaller values for the L 2 norm (an associated figure is not shown for conciseness purposes). As a result, the L ∞ and L 2 norms of Φ e for the circle and sphere immersed boundaries for a Poisson problem with the parabolic analytical solution for the QX methods with Neumann IBCs have a first-order limiting behavior. The limited improvement in the accuracy of the results despite using the quadratic method is attributed to the large first-order truncation error of the second-order Lagrange interpolations for Neumann IBCs which masks any accuracy gains associated with the quadratic method. It should be noted that as the grid size increases the QGS method yields slightly more accurate results than the QIS method (and by extension the Q method).

Figure 17 shows the L ∞ and L 2 norms of Φ e for a circle and a sphere immersed boundary for a Poisson problem with a trigonometric analytical solution for the Q, QGS, and QIS methods and have a first-order limiting behavior, except for the QGS method which increases the order of convergence from 1 to 1.6. Unlike the behavior observed for a parabolic analytical solution with Neumann IBCs, for a trigonometric analytical solution the QX methods may not be used interchangeably as the QGS method is more accurate than the Q and QIS methods. The superconvergence of the QGS method seen in Figure 17 is not yet explained. Similar to the results for Dirichlet IBCs, for Neumann IBCs with a parabolic and trigonometric analytical solution the QY O methods do not yield more accurate solutions and as they may be used interchangeably with the QX methods all comments on the latter methods are applicable to the former methods.

Quadratic Method with Stencil Size of 1 For Neumann IBCs with a parabolic analytical solution the QISS1 method do not results in a more accurate solution and may be used interchangeably with the QX and QY O methods.

Third-Order Lagrange Interpolation

As seen in Section 5.1.5, the use of third-order Lagrange (p=3) interpolations for IBMs with a maximum stencil size of 1 only increases the accuracy of results for Dirichlet IBCs as the leading term in the truncation error is related to second-order IBCs. However, as seen from equations ( 10), for Neumann IBCs the leading term in the truncation errors for the IBMs with a maximum stencil size of 1 is related to the order of interpolation used for the image points.

Figure 18 shows the L ∞ and L 2 norms of Φ e for a circle and sphere immersed boundary for a Poisson problem with a parabolic and trigonometric analytical solution for the LGS method with p=2 and p=3. As expected, the LGS method with third-order Lagrange interpolations with its maximum stencil size of 2 has a second-order convergence for both error norms and also yields significantly more accurate results than the LGS method with p=2.

Similar to the Dirichlet IBCs, the combined use of QISS1, p=3 for equation [START_REF] Langer | Comparison of geometrical and algebraic multigrid preconditioners for data-sparse boundary element matrices[END_REF], and extending the quadratic method with second-order Lagrange interpolations to the domain boundaries yields error fields with maximum values on the order of machine epsilon for the parabolic analytical solution.

Figure 19 shows the L ∞ and L 2 norms of Φ e for a circle and a sphere immersed boundary for a Poisson problem with a trigonometric analytical solution for the QIS and QISS1 methods. The QISS1 method with second-order Lagrange interpolations (p=2) is more accurate than the QIS method and increases the order of convergence from 1 to 1.6. It should be noted that a more accurate solution with the QISS1 is expected as the distances between the ghost-cell and image points decreases, thereby decreasing the error associated with the interpolations and discretizations. As seen from equation [START_REF] Langer | Comparison of geometrical and algebraic multigrid preconditioners for data-sparse boundary element matrices[END_REF] and shown in Figure 19, by extending the maximum stencil size of QISS1 from 1 to 2, through the use of p=3, the order of convergence of the error norms also increase from 1.6 to 2 for a trigonometric analytical solution.

Conclusion For Neumann IBCs and maximum stencil sizes of 1 and 2, the improvements in the L ∞ and L 2 norms of the error field Φ e obtained with the proposed shifted and quadratic methods varies with the analytical solution of the Poisson test problem.

For a maximum stencil size of 1 and a parabolic analytical solution, the proposed shifted methods (LGS, LIS, and QISS1) result in little to no improvements in the L ∞ and L 2 norms compared to the L method. However, for a trigonometric analytical solution the proposed shifted methods yield smaller L ∞ and L 2 norms and increase the order of convergence compared to the standard linear method. Regardless of the similarity in numerical accuracy obtained for a parabolic analytical solution, these square shift methods have a maximum stencil size of 1 and thereby yield computational efficient band matrices.

For a maximum stencil size of 2 and a parabolic analytical, the proposed quadratic methods (QX and QY O) also result in little to no improvements in the L ∞ and L 2 norms compared to the standard linear method and may be used interchangeably. However, for a trigonometric analytical solution the proposed quadratic methods result in smaller error norms and an increase in the order of convergence compared to the standard linear method. As expected, the proposed QISS1 method and third-order Lagrange interpolations result in significant decreases in the values of the L ∞ and L 2 norms and leads to second-order convergence. Table 3 summarizes the maximum stencil sizes and order of convergences of the various proposed methods for Neumann IBCs for parabolic and trigonometric analytical solutions for the 2D and 3D Poisson problems.

Method

Maximum Stencil Size Order of Convergence One Image Point L (p=2)

2 1 LGS (p=2) 1 1(1.5) LGS (p=3) 2 2 LIS (p=2) 1 1(1.5) Two Image Points Q (p=2) 2 1 QGS (p=2) 2 1(1.5) QGSO (p=2) 2 1 QIS (p=2) 2 1 QISS1 (p=2) 1 1 QISS1 (p=3) 2 2
Table 3: Maximum stencil sizes and order of convergences for the proposed methods for Neumann IBCs. The order of convergence shown in parenthesis corresponds to a superconvergence encountered for a specific test case.

Numerical Simulations for a Vector Flow Field

In this section the results of numerical simulations for the 3D incompressible Navier-Stokes equations with constant density and viscosity for various verification and validation test cases are presented. The governing equations follow

∂u ∂t + ∇ • (u ⊗ u) = -∇p + ν∆u + f, in Ω i , ∇ • u = 0, in Ω i , u = u D , on ∂Ω D ∪ Γ D , ∇u • n = u N , on ∂Ω N ∪ Γ N , (31)
where Ω i ⊂ IR 3 , u=(u, v, w) is the velocity flow field, p is the kinematic pressure, f is a source flow field, ν is the kinematic viscosity, u D is the Dirichlet domain boundary conditions, and u N is the Neumann domain boundary conditions. For sake of clarity, only the results for a subset of IBMs presented in previous sections are discussed. The pressure-correction method of Timmermans et al. [START_REF] Timmermans | An approximate projection scheme for incompressible flow using spectral elements[END_REF] on a staggered grid is used to solve the velocity/pressure coupling of the Navier-Stokes equations. Implicit second-order central difference is used for the non-linear advection and stress terms. Due to the current implicit advection scheme an additional immersed boundary extrapolation of the velocity field to the ghost-cells is required [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF]. Similar to the verification cases considered in Section 5.2, second-order Lagrange interpolations (p=2) for the ghost and image points are primarily considered. As will be shown, this yields pressure results with a first order convergence. Using third-order Lagrange interpolations (p=3) with the LGS or QISS1 methods yield error norms with increased order of convergence of 1.5 or 2, respectively.

Verification Test Cases

The verification test cases currently considered are the canonical 2D annular Couette flow and 3D pipe flow. As these flows have exact analytical solutions the following discussion focuses on both the L ∞ and L 2 error norms of u and pressure and on the CPU time speed-up, henceforth referred to as speed-up, attained with the currently proposed methods compared to the L method.

Annular Couette Flow

The first verification case is annular Couette flow for inner and outer cylinders with a radii of R i and R o of 0.02 and 0.1, respectively. The computational domain is [-9/80, 9/80] 2 with grid sizes of (16ν) 2 and with Dirichlet IBCs derived from their analytical solutions. The analytical solution of annular Couette flow follows

u = u θ θ , u θ (r) = Ar 2 + B r , p(r) = p 0 - (Ar) 2 8 + AB ln(r) - 1 2 B r 2 , ( 32
)
where r is the distance from the center of inner cylinder, p 0 =0, and the constants A and B follow,

A = 2 ω o R 2 o -ω i R 2 i R 2 o -R 2 i , B = (ω i -ω o ) (R i R o ) 2 R 2 o -R 2 i , (33)
where the ω i =0.5 rad and ω o =0 rad are the angular velocities for the inner and outer concentric cylinders, respectively. Figure 20 shows the L ∞ and L 2 error norms of u and pressure for various grid sizes and IBMs for annular Couette flow, wherein, as expected, all proposed methods yield more accurate results. The unexpected convergence (i.e. stair-like behavior) seen for the L ∞ error norm of u for the LGS, Q, and QGS methods is consistent with [START_REF] Picot | Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: The ghost node shifting method[END_REF] and is characteristic to many of the methods based on the linear method for the current test case. The LGS, Q, and QGS methods also locally increase the order of convergence of the L ∞ error norm of pressure between grids N 2 =64 and N 2 =256. The reduced stencil size of the LGS also results in a speed-up of 1.7 for a grid size N 2 = 128 2 .

The QISS1 method with third-order Lagrange interpolations both improves the convergence behavior (i.e. smooths out the previously seen stair-like behavior) and increase the order of convergence of for the L ∞ error norms of pressure and velocity to second order whilst maintaining the same stencil size of 2 as the L method. Extending the use of third-order Lagrange interpolations to the LGS method also yields improved convergence behaviors for the L ∞ error norm for both u and pressure (compared to the LGS method with p=2 -blue dashed line), whilst maintaining a stencil size of 2. LGS with third-order Lagrange interpolations (p=3); Green dotted line: QISS1 with p=3; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by a factor of 10 for ease of reading.

Pipe Flow

The second verification case is laminar pipe flow parallel to the x-axis with a Reynolds number Re=2Ru m /ν=20 and a radius R=3/4. The computational domain is [-1, 1] 3 with grid sizes of (16ν) 3 and with Dirichlet IBCs derived from their analytical solutions. The analytical solution of pipe Poiseuille flow follows u = u i i ,

u x (y, z) = u m 1 - r 2 R 2 , p x (x) = p 0 - ∂p ∂x x, ∂p ∂x = -4u m ν R 2 , (34)
where p 0 =0 is the initial pressure, u m =1/(2R) is the maximum velocity, ν=1/Re is the dynamic viscosity, and r 2 =(y -y c ) 2 +(z -z c ) 2 is the distance to the axis of the pipe (x-axis).

Figure 21 shows the L ∞ and L 2 error norms of u and pressure for various grid sizes and IBMs for laminar pipe flow and as expected the L method has a second and first-order limiting behavior, respectively. It should be noted that the first-order convergence for the L ∞ norm of pressure is expected as it is also encountered for the canonical validation case of 2D channel Poiseuille flow without immersed boundaries. Also as expected, the LIS and QIS methods are more accurate than the L method and the L ∞ error norm of pressure has a slightly better order of convergence.

Remark As the analytical solution is parabolic, the combined used of QISS1, third-order Lagrange interpolations (p=3), and extending the quadratic method with second-order Lagrange interpolations (p=2) to the domain boundaries yields error fields with values close to machine epsilon whilst maintaining a stencil size of 2 as the L method. Table 4 shows the total run times for the various IBMs considered for a grid size N 3 =64 3 and, as expected, the methods with a maximum stencil size of 1 (LGS, LIS, and QISS1) are faster than those methods with a maximum stencil size of 2 (L and QIS). The significance of a maximum stencil size of 1 for computational performance can be best seen by comparing the total run times of the various methods in Table 4 for the less efficient AMG preconditioner BoomerAMG of hypre library, wherein the IBMs with a maximum stencil size of 1 are at least 4 times faster than the L method. It should be noted that while QIS and L have the same stencil size of 2 and yield similar total wall clock times, the former method is more accurate than the latter.

IBM

Validation Test Cases

The validation test cases currently considered are the lid-driven cavity with an immersed cylinder, flow past a circular cylinder, flow past an isothermal and an isoflux heated circular cylinder, and flow past a sphere. As these flows do not have analytical solutions the following discussions primarily focus on the increased computational efficiency of the proposed methods for flows with more complex flow features.

Lid-Driven Cavity with an Immersed Cylinder

The first validation case is the lid-driven cavity with a stationary immersed cylinder case proposed by Cai et al. [START_REF] Cai | Moving immersed boundary method[END_REF]. A circular cylinder at Re=1000 with a diameter D=0.4L is placed in the middle of a square cavity with side length L=1 and a grid size N 2 =256 2 . The boundary conditions for the domain boundaries for velocity are uniform flow u=(1.0,0.0) for the top wall and no-slip boundary conditions for the remaining boundary, while pressure increment has zero normal gradient for all domain boundaries. Appropriate Dirichlet and Neumann IBCs for u and pressure increment, respectively, for no slip condition at the immersed boundaries are chosen. Figure 22 shows the good agreement between the current u and v-velocity profiles along the y and x-direction centerline, respectively, for the LIS method with those of Cai et al. [START_REF] Cai | Moving immersed boundary method[END_REF].

A lid-driven cavity with a rotating immersed boundary is also considered as part of the current validation test case. The numerical setup is identical to that for a stationary immersed boundary except for the angular velocity ω=1 rad of the circular cylinder and Re=5. Figure 23(a) shows the current u and v-velocity profiles along the y and x-direction centerline, respectively, and Figure 23(b) shows the distribution of the coefficient of pressure, as a function of θ, along the surface of the circular cylinder for the L and QISS1 with third-order Lagrange interpolation (p=3) methods. It should be noted that the reduced stencil size of 1 of LIS methods, compared to the stencil size of 2 of the L method, results in similar results as those presented in Figure 23 and a speed-up of 1.7 for a grid size N 2 =256 2 .

To further highlight the improvements in accuracy and convergence order of the proposed LGS and QISS1 with p=3 whilst maintaining the same stencil size of 2 as the L method, the L ∞ and L 2 error norms of pressure are considered. As the current validation test case does not have an exact solution, the results obtained on a highly resolved grid N 2 =512 2 are used as the baseline when computing the error norms. The errors norms of velocity u are not currently considered as the Dirichlet IBCs guarantee a second order convergence as seen in Section 6.1. Figure 24 shows L ∞ and L 2 error norms of pressure for the L method and, as expected, display first and in-between first and second order convergence, respectively. The results of the currently proposed LGS and QISS1 with p=3 display an improved order of convergence and are significantly more accurate than those obtained with the L method. The LGS method with p=3 has a 1.5 and 2nd order convergence for the L ∞ and L 2 error norms, respectively, while the QISS1 with p=3 displays a second order convergence for both error norms.

Laminar Flow past a Circular Cylinder

The second validation case is steady laminar flow past a circular cylinder at Re=40 with a diameter D=1.0 placed at (x,y)=(15,11)D of a computational domain Ω= [37.5D, 22D]. A non-uniform rectilinear grid of size 1320×720 with a uniform Cartesian grid of size 960×480 for the region around the circular cylinder enclosed by [-1.0, -1.0] and [3.0, 1.0] is used. Appropriate Dirichlet and Neumann IBCs for u and pressure increment, respectively, for no slip condition at the immersed boundaries are chosen. The boundary conditions for the domain boundaries for velocity are a uniform flow u=(1.0,0.0) at the inlet, zero-gradient at the outlet, and slip boundary condition for the remaining domain boundaries, while pressure increment has Dirichlet boundary conditions of zero for the outlet and zero normal gradient for the remaining domain boundaries. Table 5 shows the good agreement in the coefficient of drag and the non-dimensional wake bubble length L w /D, where L w is the dimensional wake bubble length, while Figure 25 shows the good agreement of the distribution of the coefficient of pressure, as a function of θ, along the surface of the circular cylinder for steady flow at Re=40. It should also be noted that the reduced stencil size of 1 of the LIS method yields a speed-up of 1.7 compared to the linear method.

Laminar Flow past a Heated Circular Cylinder

The third validation case is steady laminar flow past an isothermal and an isoflux heated circular cylinder at Re=20 and Prandtl number P r=0.70 with a diameter D=1.0 placed at (x,y)=( uniform Cartesian grid of sizes 640×320 for the region around the circular cylinder enclosed by [-1.0, -1.0] and [3.0, 1.0] are used. The boundary conditions for u and pressure increment are the same as those presented in Section 6.2.2. The boundary conditions for temperature at the immersed boundary is uniform temperature T =1 and uniform heat flux ∂T /∂n=1 for the isothermal and isoflux cases, respectively, while a uniform temperature T =0 is applied at the inlet and zerogradient ∂T /∂n is used for the remaining domain boundaries. Table 6 shows the good agreement in the average Nusselt number while Figure 27 shows the good agreement in the distribution of the local Nusselt number along the surface of the heated circular cylinder.

Isothermal Isoflux

Laminar Flow past a Sphere

The fourth validation case, shown in Figure 28, is steady laminar flow past a sphere at Re=50, 100, and 150 with a diameter D=1.0 placed at (x,y,z)=(12.5,8.5,8.5)D of a computational domain Ω=[31.25D, 17D, 17D]. A non-uniform rectilinear grid of size 275×150 2 with a uniform Cartesian grid of size 160×80 2 for the region around the sphere enclosed by [-1.0, -1.0, -1.0] and [3.0, 1.0, 1.0] is used for all Reynolds numbers considered. Appropriate Dirichlet and Neumann IBCs for u and pressure increment, respectively, for no slip condition at the immersed boundaries are chosen. The boundary conditions for the domain boundaries for velocity are a uniform flow u=(1.0,0.0,0.0) at the inlet, zero-gradient at the outlet and slip boundary condition for the remaining domain boundaries, while pressure increment has Dirichlet boundary conditions of zero for the outlet and zero normal gradient for the remaining domain boundaries. Table 7 shows the good agreement in the coefficient of drag C D and the non-dimensional wake bubble length L w /D, where L w is the wake bubble length. The need for a maximum stencil size of 1 is indeed important given that the LIS method is 3.1 times faster than the L method for steady flow at Re=50. Although the results in Table 7

Conclusion

Although IBMs with a stencil size of 2 for Cartesian grids, such as the linear method [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF], are widely used, they pose significant constraints in the pursuit of smaller total run times, smaller memory requirements, increased accuracy, and increase order of convergence. To address these shortcomings the linear square shifting and the quadratic square shifting ghost-cell methods, and their various combinations are proposed. The proposed linear square shifting methods consistently yield a maximum stencil size of 1 by shifting the ghost point and, consequently, the image point, or solely the image point, towards the immersed boundary, whilst the quadratic square shifting ghostcell methods derived from [START_REF] Das | A sharp-interface immersed boundary method to simulate convective and conjugate heat transfer through highly complex periodic porous structures[END_REF] increases the accuracy of the results through quadratic interpolations for the ghost point/cell values whilst maintaining a maximum stencil size of 2. The ghost-cell and image point linear square shift methods (LGS and LIS, respectively), with a maximum stencil size of 1, result in smaller total run times, smaller memory requirements, increased accuracy, and improved convergence, while the quadratic square shifting methods with ghost or image point shifts, with a maximum stencil size of 2, are proposed to further increase the accuracy and improve the convergence of the results. Lastly, the quadratic square shifting method with the first image point shift and with a maximum stencil size artificially increased from 1 to 2 through third-order Lagrange interpolations is also proposed to increase the order of convergence of Neumann immersed boundary conditions from first to second-order. Numerical simulations of the canonical Poisson test problem with parabolic and trigonometric analytical solutions for various 2D and 3D immersed boundaries with Dirichlet and Neumann IBCs for nine combinations of the square shift and quadratic square shifting methods further emphasises the need for a maximum stencil size of 1 and more accurate interpolation methods for the ghost point/cells. For Dirichlet IBCs the reduced stencil size of the linear square shift methods tends to yield more accurate results and improved second-order convergence when compared with the linear method. The accuracy of the results is further increased when the quadratic method is used in conjunction with the square shifting methods (i.e. QGS, QIS, etc.).

Although for Neumann IBCs the proposed methods yield smaller improvements in accuracy and negligible improvements in the first-order convergence of the linear method for the parabolic analytical solution to the Poisson test problem, the reduced stencil size results in smaller run times.

However, the combined use of the QISS1 method and third-order Lagrange interpolations for the image and ghost points yields more accurate results and error norms with second-order convergence while maintaining the same maximum stencil size of 2 as the linear method, which results in firstorder convergence.

The numerical accuracy and computational improvements of the square shifting methods, quadratic square shifting methods, and the combinations of these methods are also observed for the velocity and pressure variables for various verification and validation test cases for the Navier-Stokes governing equations. For the verification test cases, the proposed methods yield more accurate results and either improve the convergence behavior or increase the order of convergence of the error norms. For methods with a maximum stencil size of 1 (e.g. LIS) these improvements are accompanied by significant speed-ups. As seen from the various validation cases, the proposed methods are able to accurately capture the velocity, pressure, and temperature flow fields. In addition the combined use of the QISS1 method and third-order Lagrange interpolations for the image and ghost points results in second-order convergence for the pressure field with Neumann IBCs for a lid driven cavity with a rotating immersed cylinder. Similar to the verification test cases, the proposed methods with a maximum stencil size of 1 also yield speed-ups of at least 3 for 3D flows.

Figure 1 :

 1 Figure 1: (a) Inner and outer computational domains, Ω i and Ω o , respectively, for an immersed boundary Γ with an outward pointing local surface unit normal n and domain boundary conditions ∂Ω. (b) Finite volume discretization of (a) and classification of cells into inner (triangle), outer (circle), and ghost (diamond) cells.

Figure 4 :

 4 Figure 4: Close up view of an image point interpolation set and stencil size for a ghost-cell for Cartesian grids. Solid orange line: Γ; Dashed blue line: stencil size of 1; Dashed orange line: stencil size of 2; Dashed red line: interpolation set for enclosed image point; Black diamond: ghost-cell; Black circle: outer-cell; Black triangle: inner-cell; Black ×: image point.

Figure 6 :

 6 Figure 6: Close up view of possible quadrant intersection points x QIP for Q1 interpolation quadrant for a sample boundary intercept x BI . Solid orange line: Γ; Dotted red line: Local axis (x , y ) with origin at x BI ; Solid green vector: sample local surface unit normal to intersect line segments [X i,l+1 ,X k,l+1 ] and [X k,l+1 ,X k,l ]; Black diamond: ghost-cell; Black triangle: inner-cell.

Figure 7 :

 7 Figure 7: Close up view of the supplemental data points and values required for the quadratic methods for a sample ghost-cell X i,j . BI: Boundary intercept; IP : Image point; IP 2: Second image point; n: local surface unit normal.

Figure 9 :

 9 Figure 9: Computational domain for (a) circle, (b) 2D flower, (c) sphere, and (d) 3D flower immersed boundaries Γs for (a,b) 2D and (c,d) 3D verification test problems with a parabolic analytical solution. Solid orange line/surface: Γ; r: radius.

Figure 10 :

 10 Figure 10: L ∞ and L 2 norms of Φ e for various IBMs for a (a) 2D flower, (b) sphere, and (c) 3D flower immersed boundary with Dirichlet IBC for a (a,b) parabolic and (c) trigonometric analytical solution. Black solid line: L; Blue dashed line: LGS; Orange dotted line: LIS; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Figure 12

 12 (a) shows the L ∞ and L 2 norms of Φ e for a circle immersed boundary for a Poisson problem with a parabolic analytical solution for QX and QY O methods, where the image points have different sets of interpolation indices. As seen by comparing the Q and QO methods, and the QGS and QGSO methods, the QY O methods do not results in a more accurate solution and as they may be used interchangeably with the QX methods all comparison in Section 5.1.4 are applicable to the QY O methods. Quadratic Method with Stencil Size of 1 Figures 12(b) and 12(c) show the L ∞ and L 2

Figure 11 :

 11 Figure 11: L ∞ and L 2 norms of Φ e for various IBMs for a (a) circle, (b) sphere, and (c) 3D flower immersed boundary with Dirichlet IBC for a (a,b) trigonometric and (c) parabolic analytical solution. Black solid line: L; Red solid line: Q; Blue dotted line: QGS; Grey dashed line: QIS; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by factor of 10 for ease of reading. Blue dotted line and grey dashed line are superimposed on each other in (b) and (c).

Figure 12 :

 12 Figure 12: L ∞ and L 2 norms of Φ e for various IBMs for a (a) circle, (b) 2D flower, and (c) sphere immersed boundary with Dirichlet IBC for a (a) parabolic and (b,c) trigonometric analytical solution. Black solid line: L; Red solid line: Q; Grey dotted line: QO; Blue dashed line: QGS; Green dotted line: QGSO; Blue dotted line: QIS; Purple dashed line: QISS1; N : number of cells in the x-axis; Diamond symbols: L ∞ ; Circle symbols: L 2 norms. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Figure 13 :

 13 Figure 13: L ∞ and L 2 norms of Φ e for various IBMs for a (a,b) circle and (c) sphere immersed boundary with Dirichlet IBC for (a,c) parabolic and (b) trigonometric analytical solution. Black solid line: L; Orange solid line: LGS with p=2; Red solid line: LGS with p=3; N : number of cells in the x-axis; Diamond symbols: L ∞ ; Circle symbols: L 2 norms. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Figure 14 :

 14 Figure 14: L ∞ and L 2 norms of Φ e for a (a) circle and a (b) sphere immersed boundary with Dirichlet IBC and trigonometric analytical solution. Black line: p=2 ; Red line: p=3; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

  =πcos(πx)cos(2πy)n x -2πsin(πx)sin(2πy)n y ∀(x) ∈ Γ ∈ IR 2 , φ(0, y) =0, φ(1, y) =πcos(πx)cos(2πy)cos(2πz)n x -2πsin(πx)sin(2πy)cos(2πz)n y -

Figure 15 :

 15 Figure 15: L ∞ and L 2 norms of Φ e for various IBMs for a (a) circle and a (b) sphere immersed boundary with Neumann IBC. Black solid line: L; Blue dashed line: LGS; Orange dotted line: LIS; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Figure 16 :

 16 Figure 16: L ∞ and L 2 norms of Φ e for various IBMs for a (a) circle and a (b) sphere immersed boundary with Neumann IBC and trigonometric analytical solution. Black solid line: L; Blue dashed line: LGS; Orange dotted line: LIS; N : number of cells in the x-axis; Diamond symbols: L ∞ ; Circle symbols: L 2 norms. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Figure 17 :

 17 Figure 17: L ∞ and L 2 norms of Φ e for various IBMs for a (a) circle and a (b) sphere immersed boundary with Neumann IBC and trigonometric analytical solution. Black solid line: L; Red solid line: Q; Blue dotted line: QGS; Grey dashed line: QIS; N : number of cells in the x-axis; Diamond symbols: L ∞ ; Circle symbols: L 2 norms. The L 2 norms are shifted downwards by factor of 10 for ease of reading. The black and red solid line are superimposed on each other in (a) and (b).

Figure 18 :

 18 Figure 18: L ∞ and L 2 norms of Φ e for various IBMs for a (a,b) circle and (c) sphere immersed boundary with Dirichlet IBC for (a,c) parabolic and (b) trigonometric analytical solution. Black solid line: L; Orange solid line: LGS with p=2; Red solid line: LGS with p=3; N : number of cells in the x-axis; Diamond symbols: L ∞ ; Circle symbols: L 2 norms. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Figure 19 :

 19 Figure 19: L ∞ and L 2 norms of Φ e for a (a) circle and a (b) sphere immersed boundary with Neumann IBC and trigonometric analytical solution. Black line: L; Grey dotted line: QIS; Red line: QISS1 p=2 ; Orange line: QISS1 p=3; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by factor of 10 for ease of reading.

Figure 20 :

 20 Figure 20: L ∞ and L 2 norms of (a) u and (b) pressure for various IBMs for annular Couette flow. Black solid line: L; Red solid line: LGS; Blue dashed line: Q; Orange dotted line: QGS; Grey dotted line:LGS with third-order Lagrange interpolations (p=3); Green dotted line: QISS1 with p=3; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by a factor of 10 for ease of reading.

Figure 21 :

 21 Figure 21: L ∞ and L 2 norms of (a) u and (b) pressure for various IBMs for laminar pipe flow at Re=20. Black solid line: L; Blue dashed line: LIS; Orange dotted line: QIS; N : number of cells in the x-axis; Diamond symbols: L ∞ norm; Circle symbols: L 2 norm. The L 2 norms are shifted downwards by a factor of 10 for ease of reading.

Figure 22 :

 22 Figure 22: (a) Streamlines and (b) profiles of u and v-velocity along the y and x-direction centerline, respectively, for lid-driven cavity with a stationary immersed boundary at Re=1000. Blue line: uvelocity; Red line: v-velocity; Black +: Cai et al. [27]. y and x in (b) correspond to r/L, where r is the Euclidean distance from (-0.5,0.0) or (0.0,-0.5) and L=1 is the side length.

Figure 23 :

 23 Figure 23: Profiles of (a) u and v-velocity along the y and x-direction centerline and (b) the distribution of the coefficient of pressure along the surface of the immersed boundary for lid-driven cavity with a rotating immersed boundary with an angular velocity ω=1 rad and a Re=5. (a) Blue line: u-velocity; Red line: v-velocity. (b) Blue line: L method; Red line: QISS1 method with third-order Lagrange interpolation (p=3).

Figure 24 :

 24 Figure 24: (a)L ∞ and (b)L 2 norms of pressure for the L method, the LGS method with p=3, and the QISS1 method with p=3 for lid-driven cavity with a rotating immersed boundary with ω=1 rad at Re=5. Blue solid line: L method; Green solid solid: LGS method with p=3; Red solid line: QISS1 method with p=3; Blue dashed line: O(h); Red dashed line: O(h 2 );N : number of cells in the x-axis.

Figure 25 :

 25 Figure 25: (a) Pressure contour lines and (b) normalized coefficient of pressure C p distribution along the surface of a circular cylinder for steady flow at Re=40. Solid black line: Current results with the LIS method; Blue +: Dennis and Chang[START_REF] Dennis | Numerical solutions for steady flow past a circular cylinder at reynolds numbers up to 100[END_REF]; Orange square: Berthelsen and Faltinsen[START_REF] Berthelsen | A local directional ghost cell approach for incompressible viscous flow problems with irregular boundaries[END_REF]; Purple diamond: Grove et al.[START_REF] Grove | An experimental investigation of the steady separated flow past a circular cylinder[END_REF]; Green o: Tseng and Ferziger[START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF].

Figure 26 :

 26 Figure 26: Temperature contour lines for forced convection steady laminar flow past an (a) isothermal and (b) isoflux at Re=20.

Figure 27 :

 27 Figure 27: Local Nusselt number for forced convection steady laminar flow past an (a) isothermal and an (b) isoflux circle at Re=20 and P r=0.70. Blue +: Bharti et al. [38]; Orange line: Zhang et al. [37]; Black line: Current results.

Figure 28 :

 28 Figure 28: Streamlines coloured by velocity magnitude for steady laminar flow past a sphere immersed boundary at Re=150.

x GP x IP x IP x BI

  IP =x IP in this figure;x QIP : quadrant intersection point; ∆l:||X i,j -x IP || 2 (L); ∆l : ||x GP -x IP || 2 (ghost-cell shift); ∆l 1 : ||X i,j -x BI || 2 (image point shift); ∆l 2 : ||x BI -x IP || 2 (image point shift).

			X k,l+1
	∆l	∆l ∆l 1	n ∆l 2	x QIP
		X i,j		X k,l
	Figure 5: Close up view of the supplemental data points and values required for the linear (L) and
	linear square shift methods for a sample ghost-cell X i,j . x GP : ghost point (ghost-cell shift); x IP :
	image point (L); x IP : image point (ghost-cell shift); x IP (not shown): image point (image point
	shift) and x

Table 4 :

 4 Total wall clock time, stencil size, and speed-up with respect to the L method for various IBMs for laminar pipe flow at Re=20 with a grid size N 3 =64 3 .

		Stencil Size Total Wall Clock time (×10 2 )[s] Speed-up (vs. L method)
	L	2	5.02	-
	LGS	1	1.18	4.25
	LIS	1	1.21	4.15
	QIS	2	4.12	1.22
	QISS1	1	1.18	4.25

Table 5 :

 5 Coefficient of drag C D and non-dimensional wake bubble length L w /D for steady flow past a circle at Re=40.

	15,11)D

Table 6 :

 6 Average Nusselt number N u for forced convection steady laminar flow past an isothermal and an isoflux circle at Re=20 and P r=0.70.

	Luo et al. [8]	2.4336	2.7850
	Zhang et al. [37]	2.47	2.75
	Bharti et al. [38]	2.4653	2.7788
	Dennis et al. [39]	2.5216	-
	Lange et al. [40]	2.4087	-
	Soares et al. [41]	2.4300	-
	Ahmad and Qureshi [42]	-	2.6620
	Current(LIS)	2.4272	2.7916

Table 7 :

 7 are for the LIS method, similar results are obtained with the L method. L w /D C D L w /D C D L w /D Coefficient of drag C D and non-dimensional wake bubble length L w /D for steady flow past a sphere at Re=50, 100, and 150.

	Reynolds Number	50			100		150
	C D Mittal [43] 1.57	0.44	1.09	0.87	-	-
	Johnson and Patel et al. [44] 1.57	0.40	1.08	0.86	0.90	1.20
	Marella et al. [45]	1.56	0.39	1.06	0.90	0.85	1.19
	Current(LIS)	1.59	0.40	1.09	0.86	0.89	1.18

Acknowledgements

We acknowledge the French National Research Agency for its support (ANR-17-CE07-0029-SUPERFON) and the MCIA (Mésocentre de Calcul Intensif Aquitaine) for the HPC resources.