A fourth-order accurate compact difference scheme for solving the   three-dimensional Poisson equation with arbitrary boundaries

Shirzad Hosseinverdi; Hermann F. Fasel

arXiv:1901.08269·physics.comp-ph·December 20, 2024

A fourth-order accurate compact difference scheme for solving the three-dimensional Poisson equation with arbitrary boundaries

Shirzad Hosseinverdi, Hermann F. Fasel

PDF

TL;DR

This paper introduces a high-order compact difference scheme combined with a preconditioned BiCGSTAB method for efficiently solving the 3D Poisson equation with irregular boundaries, maintaining fourth-order accuracy without jump corrections.

Contribution

The work develops a fourth-order accurate 3D Poisson solver that remains precise and efficient near immersed boundaries without requiring jump corrections, unlike previous methods.

Findings

01

Achieves fourth-order accuracy on irregular domains.

02

Maintains efficiency regardless of boundary complexity.

03

Demonstrates superior convergence and computational performance.

Abstract

This paper presents an efficient high-order sharp-interface method for solving the three-dimensional (3D) Poisson equation with Dirichlet boundary conditions on a nonuniform Cartesian grid with irregular domain boundaries. The new approach is based on the combination of the fourth-order compact finite difference scheme and the preconditioned stabilized biconjugate-gradient (BiCGSTAB) method. Contrary to the original immersed interface method by LeVeque and Li [1], the new method does not require jump corrections, instead, the (regular) compact finite difference stencil is adjusted at the irregular grid points (in the vicinity of the interfaces of the immersed bodies) to obtain a solution that is sharp across the interface while keeping the fourth-order global accuracy. The contribution of the present work is the design of a fourth-order-accurate 3D Poisson solver whose accuracy and…

Equations190

\displaystyle\bigg{(}\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+\frac{\partial^{2}}{\partial z^{2}}\bigg{)}~{}u(x,y,z)

\displaystyle\bigg{(}\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+\frac{\partial^{2}}{\partial z^{2}}\bigg{)}~{}u(x,y,z)

u (x, y, z)

(L_{xx} u_{xx})_{j, k}

(L_{xx} u_{xx})_{j, k}

(L_{y y} u_{y y})_{i, k}

(L_{z z} u_{z z})_{i, j}

(L_{xx} u_{xx})_{j, k}

(L_{xx} u_{xx})_{j, k}

(L_{y y} u_{y y})_{i, k}

(L_{z z} u_{z z})_{i, j}

(R_{xx} u)_{j, k}

(R_{xx} u)_{j, k}

(R_{y y} u)_{i, k}

(R_{z z} u)_{i, j}

a_{i}^{x}

a_{i}^{x}

b_{i}^{x}

c_{i}^{x}

a_{j}^{y} a_{k}^{z} (L_{xx} u_{xx} - R_{xx} u)_{j - 1, k - 1}

a_{j}^{y} a_{k}^{z} (L_{xx} u_{xx} - R_{xx} u)_{j - 1, k - 1}

\displaystyle a^{y}_{j}~{}b^{z}_{k}~{}(L_{xx}u_{xx}-R_{xx}u)_{j-1,k{\color[rgb]{1,1,1}+1}}

a_{j}^{y} c_{k}^{z} (L_{xx} u_{xx} - R_{xx} u)_{j - 1, k + 1}

a_{i}^{x} a_{k}^{z} (L_{y y} u_{y y} - R_{y y} u)_{i - 1, k - 1}

\displaystyle a^{x}_{i}~{}b^{z}_{k}~{}(L_{yy}u_{yy}-R_{yy}u)_{i-1,k{\color[rgb]{1,1,1}+1}}

a_{i}^{x} c_{k}^{z} (L_{y y} u_{y y} - R_{y y} u)_{i - 1, k + 1}

a_{i}^{x} a_{j}^{y} (L_{z z} u_{z z} - R_{z z} u)_{i - 1, j - 1}

a_{i}^{x} b_{j}^{y} (L_{z z} u_{z z} - R_{z z} u)_{i - 1, j}

a_{i}^{x} c_{j}^{y} (L_{z z} u_{z z} - R_{z z} u)_{i - 1, j + 1}

\displaystyle A_{i,j,k}^{1{\color[rgb]{1,1,1}1}}~{}u_{i-1,j+1,k-1}

\displaystyle A_{i,j,k}^{1{\color[rgb]{1,1,1}1}}~{}u_{i-1,j+1,k-1}

\displaystyle A_{i,j,k}^{4{\color[rgb]{1,1,1}1}}~{}u_{i-1,j{\color[rgb]{1,1,1}+1},k-1}

\displaystyle A_{i,j,k}^{7{\color[rgb]{1,1,1}1}}~{}u_{i-1,j-1,k-1}

\displaystyle A_{i,j,k}^{10}~{}u_{i-1,j+1,k{\color[rgb]{1,1,1}+1}}

\displaystyle A_{i,j,k}^{13}~{}u_{i-1,j{\color[rgb]{1,1,1}+1},k{\color[rgb]{1,1,1}+1}}

\displaystyle A_{i,j,k}^{16}~{}u_{i-1,j-1,k{\color[rgb]{1,1,1}+1}}

A_{i, j, k}^{19} u_{i - 1, j + 1, k + 1}

\displaystyle A_{i,j,k}^{22}~{}u_{i-1,j{\color[rgb]{1,1,1}+1},k+1}

A_{i, j, k}^{25} u_{i - 1, j - 1, k + 1}

Q_{i, j, k}

Q_{i, j, k}

\displaystyle~{}+q_{i,j,k}^{4{\color[rgb]{1,1,1}1}}~{}f_{i-1,j{\color[rgb]{1,1,1}+1},k-1}

\displaystyle~{}+q_{i,j,k}^{7{\color[rgb]{1,1,1}1}}~{}f_{i-1,j-1,k-1}

\displaystyle~{}+q_{i,j,k}^{10}~{}f_{i-1,j+1,k{\color[rgb]{1,1,1}+1}}

\displaystyle~{}+q_{i,j,k}^{13}~{}f_{i-1,j{\color[rgb]{1,1,1}+1},k{\color[rgb]{1,1,1}+1}}

\displaystyle~{}+q_{i,j,k}^{16}~{}f_{i-1,j-1,k{\color[rgb]{1,1,1}+1}}

+ q_{i, j, k}^{19} f_{i - 1, j + 1, k + 1}

\displaystyle~{}+q_{i,j,k}^{22}~{}f_{i-1,j{\color[rgb]{1,1,1}+1},k+1}

+ q_{i, j, k}^{25} f_{i - 1, j - 1, k + 1}

A_{i, j, k}^{9 α + 1} A_{i, j, k}^{9 α + 2} A_{i, j, k}^{9 α + 3} A_{i, j, k}^{9 α + 4} A_{i, j, k}^{9 α + 5} A_{i, j, k}^{9 α + 6} A_{i, j, k}^{9 α + 7} A_{i, j, k}^{9 α + 8} A_{i, j, k}^{9 α + 9} = a r_{i}^{x} c_{j}^{y} β_{k}^{z} + a_{i}^{x} c r_{j}^{y} β_{k}^{z} + a_{i}^{x} c_{j}^{y} β r_{k}^{z}, = b r_{i}^{x} c_{j}^{y} β_{k}^{z} + b_{i}^{x} c r_{j}^{y} β_{k}^{z} + b_{i}^{x} c_{j}^{y} β r_{k}^{z}, = c r_{i}^{x} c_{j}^{y} β_{k}^{z} + c_{i}^{x} c r_{j}^{y} β_{k}^{z} + c_{i}^{x} c_{j}^{y} β r_{k}^{z}, = a r_{i}^{x} b_{j}^{y} β_{k}^{z} + a_{i}^{x} b r_{j}^{y} β_{k}^{z} + a_{i}^{x} b_{j}^{y} β r_{k}^{z}, = b r_{i}^{x} b_{j}^{y} β_{k}^{z} + b_{i}^{x} b r_{j}^{y} β_{k}^{z} + b_{i}^{x} b_{j}^{y} β r_{k}^{z}, = c r_{i}^{x} b_{j}^{y} β_{k}^{z} + c_{i}^{x} b r_{j}^{y} β_{k}^{z} + c_{i}^{x} b_{j}^{y} β r_{k}^{z}, = a r_{i}^{x} a_{j}^{y} β_{k}^{z} + a_{i}^{x} a r_{j}^{y} β_{k}^{z} + a_{i}^{x} a_{j}^{y} β r_{k}^{z}, = b r_{i}^{x} a_{j}^{y} β_{k}^{z} + b_{i}^{x} a r_{j}^{y} β_{k}^{z} + b_{i}^{x} a_{j}^{y} β r_{k}^{z}, = c r_{i}^{x} a_{j}^{y} β_{k}^{z} + c_{i}^{x} a r_{j}^{y} β_{k}^{z} + c_{i}^{x} a_{j}^{y} β r_{k}^{z} . q_{i, j, k}^{9 α + 1} q_{i, j, k}^{9 α + 2} q_{i, j, k}^{9 α + 3} q_{i, j, k}^{9 α + 4} q_{i, j, k}^{9 α + 5} q_{i, j, k}^{9 α + 6} q_{i, j, k}^{9 α + 7} q_{i, j, k}^{9 α + 8} q_{i, j, k}^{9 α + 9} = a_{i}^{x} c_{j}^{y} β_{k}^{z}, = b_{i}^{x} c_{j}^{y} β_{k}^{z}, = c_{i}^{x} c_{j}^{y} β_{k}^{z}, = a_{i}^{x} b_{j}^{y} β_{k}^{z}, = b_{i}^{x} b_{j}^{y} β_{k}^{z}, = c_{i}^{x} b_{j}^{y} β_{k}^{z}, = a_{i}^{x} a_{j}^{y} β_{k}^{z}, = b_{i}^{x} a_{j}^{y} β_{k}^{z}, = c_{i}^{x} a_{j}^{y} β_{k}^{z} .

A_{i, j, k}^{9 α + 1} A_{i, j, k}^{9 α + 2} A_{i, j, k}^{9 α + 3} A_{i, j, k}^{9 α + 4} A_{i, j, k}^{9 α + 5} A_{i, j, k}^{9 α + 6} A_{i, j, k}^{9 α + 7} A_{i, j, k}^{9 α + 8} A_{i, j, k}^{9 α + 9} = a r_{i}^{x} c_{j}^{y} β_{k}^{z} + a_{i}^{x} c r_{j}^{y} β_{k}^{z} + a_{i}^{x} c_{j}^{y} β r_{k}^{z}, = b r_{i}^{x} c_{j}^{y} β_{k}^{z} + b_{i}^{x} c r_{j}^{y} β_{k}^{z} + b_{i}^{x} c_{j}^{y} β r_{k}^{z}, = c r_{i}^{x} c_{j}^{y} β_{k}^{z} + c_{i}^{x} c r_{j}^{y} β_{k}^{z} + c_{i}^{x} c_{j}^{y} β r_{k}^{z}, = a r_{i}^{x} b_{j}^{y} β_{k}^{z} + a_{i}^{x} b r_{j}^{y} β_{k}^{z} + a_{i}^{x} b_{j}^{y} β r_{k}^{z}, = b r_{i}^{x} b_{j}^{y} β_{k}^{z} + b_{i}^{x} b r_{j}^{y} β_{k}^{z} + b_{i}^{x} b_{j}^{y} β r_{k}^{z}, = c r_{i}^{x} b_{j}^{y} β_{k}^{z} + c_{i}^{x} b r_{j}^{y} β_{k}^{z} + c_{i}^{x} b_{j}^{y} β r_{k}^{z}, = a r_{i}^{x} a_{j}^{y} β_{k}^{z} + a_{i}^{x} a r_{j}^{y} β_{k}^{z} + a_{i}^{x} a_{j}^{y} β r_{k}^{z}, = b r_{i}^{x} a_{j}^{y} β_{k}^{z} + b_{i}^{x} a r_{j}^{y} β_{k}^{z} + b_{i}^{x} a_{j}^{y} β r_{k}^{z}, = c r_{i}^{x} a_{j}^{y} β_{k}^{z} + c_{i}^{x} a r_{j}^{y} β_{k}^{z} + c_{i}^{x} a_{j}^{y} β r_{k}^{z} . q_{i, j, k}^{9 α + 1} q_{i, j, k}^{9 α + 2} q_{i, j, k}^{9 α + 3} q_{i, j, k}^{9 α + 4} q_{i, j, k}^{9 α + 5} q_{i, j, k}^{9 α + 6} q_{i, j, k}^{9 α + 7} q_{i, j, k}^{9 α + 8} q_{i, j, k}^{9 α + 9} = a_{i}^{x} c_{j}^{y} β_{k}^{z}, = b_{i}^{x} c_{j}^{y} β_{k}^{z}, = c_{i}^{x} c_{j}^{y} β_{k}^{z}, = a_{i}^{x} b_{j}^{y} β_{k}^{z}, = b_{i}^{x} b_{j}^{y} β_{k}^{z}, = c_{i}^{x} b_{j}^{y} β_{k}^{z}, = a_{i}^{x} a_{j}^{y} β_{k}^{z}, = b_{i}^{x} a_{j}^{y} β_{k}^{z}, = c_{i}^{x} a_{j}^{y} β_{k}^{z} .

a_{j}^{y} a_{k}^{z} (\overline{L_{xx} u_{xx}} - \overline{R_{xx} u})_{j - 1, k - 1}

a_{j}^{y} a_{k}^{z} (\overline{L_{xx} u_{xx}} - \overline{R_{xx} u})_{j - 1, k - 1}

\displaystyle a^{y}_{j}~{}b^{z}_{k}~{}(L_{xx}u_{xx}-R_{xx}u)_{j-1,k{\color[rgb]{1,1,1}+1}}

a_{j}^{y} c_{k}^{z} (L_{xx} u_{xx} - R_{xx} u)_{j - 1, k + 1}

a_{i}^{x} a_{k}^{z} (L_{y y} u_{y y} - R_{y y} u)_{i - 1, k - 1}

\displaystyle a^{x}_{i}~{}b^{z}_{k}~{}(L_{yy}u_{yy}-R_{yy}u)_{i-1,k{\color[rgb]{1,1,1}+1}}

a_{i}^{x} c_{k}^{z} (L_{y y} u_{y y} - R_{y y} u)_{i - 1, k + 1}

a_{i}^{x} a_{j}^{y} (L_{z z} u_{z z} - R_{z z} u)_{i - 1, j - 1}

a_{i}^{x} b_{j}^{y} (L_{z z} u_{z z} - R_{z z} u)_{i - 1, j}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

A fourth-order accurate compact difference scheme for solving the three-dimensional Poisson equation with arbitrary boundaries

Shirzad Hosseinverdi

[email protected]

Hermann F. Fasel

Department of Aerospace and Mechanical Engineering, University of Arizona, Tucson, AZ 85721, USA

Abstract

This paper presents an efficient high-order sharp-interface method for solving the three-dimensional (3D) Poisson equation with Dirichlet boundary conditions on a nonuniform Cartesian grid with irregular domain boundaries. The new approach is based on the combination of the fourth-order compact finite difference scheme and the preconditioned stabilized biconjugate-gradient (BiCGSTAB) method. Contrary to the original immersed interface method by LeVeque & Li [1], the new method does not require jump corrections, instead, the (regular) compact finite difference stencil is adjusted at the irregular grid points (in the vicinity of the interfaces of the immersed bodies) to obtain a solution that is sharp across the interface while keeping the fourth-order global accuracy. The contribution of the present work is the design of a fourth-order-accurate 3D Poisson solver whose accuracy and efficiency does not deteriorate in the presence of an immersed boundary. This is attributed to (i) the modification of the discrete operators near immersed boundaries does not lead to a wide grid stencil due to the compact nature of the discretization and (ii) a preconditioning technique whose efficiency and cost are independent of the complexity of the geometry and the presence or not of an immersed boundary. The accuracy and computational efficiency of the proposed algorithm is demonstrated and validated over a range of problems including smooth and irregular boundaries. The test cases show that the new method is fourth-order accurate in the maximum norm whether an immersed boundary is present or not, on uniform or nonuniform grids. Furthermore, the efficiency of the preconditioned BiCGSTAB is demonstrated with regard to convergence rate and “extra” floating-point operation ( $FLOP_{extra}$ ) which is due to the presence of immersed boundaries. It is shown that the solution method is equally efficient for domains with and without irregular boundaries, with a negligible $FLOP_{extra}$ in the presence of immersed boundaries.

keywords:

3D Poisson equation , immersed interface method , compact finite difference , high-order

††journal: Journal of Computational Physics

1 Introduction

In this paper, a new numerical method is presented for solving the three-dimensional (3D) Poisson equation with Dirichlet boundary conditions in the 3D rectangular box $\Omega$ which can contain an arbitrary immersed body (IB) with boundary $\Gamma$ ,

[TABLE]

where $u(x,y,z)$ is an unknown function defined in the domain $\Omega$ with prescribed boundary values $g(x,y,z)$ on its boundary $\partial\Omega$ and $f(x,y,z)$ is a known source function. There has been a great deal of research work on the development of high-order finite difference numerical solution of 3D Poisson equations in the past two decades [2, 3, 4, 5, 6], with particular focus on the fourth-order accurate compact difference scheme. However, the majority of these schemes are limited to simple domains with uniform mesh distributions. For many practical applications, however, the numerical computation of the 3D Poisson equation in a domain with irregular boundaries is required. In the presence of an IB, the domain $\Omega$ is divided by the surface $\Gamma$ into two subdomains $\Omega^{+}$ and $\Omega^{-}$ corresponding to the outside and inside of the immersed body, respectively (see Fig. 1). One would typically solve Eq. (1) defined on the open region $\Omega^{+}$ with boundary conditions on $\partial\Omega$ (the outer boundary which conforms to the computational domain), and $\Gamma$ , the immersed boundary. The solution in the region $\Omega^{-}$ may or may not be of interest. In either case, the immersed boundary $\Gamma$ represents a singularity, thus field variables and/or their derivatives will be discontinuous across the immersed boundary. In the present work, we refer to $\Omega^{-}$ and $\Omega^{+}$ solid and fluid regions, respectively, and the solution inside the immersed body, $\Omega^{-}$ , is trivial and set to zero.

Most of the numerical algorithms capable of handling complex geometries use body-fitted structured or unstructured grids. However, generating high-quality structured grids is generally cumbersome and it becomes very laborious as the complexity of the geometry increases. The grid generation process becomes a very difficult task, even for the simplest geometries, when more than one body is located within the domain. Grids with poor qualities (smoothness, orthogonality, aspect ratio) could negatively impact the accuracy and convergence properties of the numerical method for solving the Poisson equation. On the other hand, unstructured grids which are more suited for complex geometries, suffer from slow rate of convergence which could lead to a substantial increase in the computation time [7]. Compared to structured grids, unstructured grids require a larger amount of memory and the extension of the numerical schemes to higher-order is not straightforward.

An alternative approach is to discretize Eq. (1) on a fixed Cartesian grid which is allowed to intersect with an immersed boundary $\Gamma$ . As a result, the grid does not conform to the solid boundaries. The main advantage in this case lies in the use of Cartesian grids for which there is almost no grid generation cost and high-order numerical methods developed for Cartesian grids, can be employed. A main challenge here is the imposition of boundary conditions on the immersed boundaries such that the accuracy and efficiency of the Cartesian solvers is maintained. Several methods have been proposed in the past to handle the singularity associated with the immersed boundary. Generally, they can be classified into two categories: the so-called immersed boundary methods (IBM), and immersed interface methods (IIM).

The original IBM was pioneered by Peskin [8, 9] to handle elastic boundaries for simulating blood flow in the heart. In Peskin‘s approach, the boundary conditions are enforced through a smooth forcing term added to the Navier-Stokes equations. This type of IBM is classified as continuous forcing (diffuse) approach. One disadvantage of diffuse methods is that the effect of the boundary is distributed over a band of several grid points which causes a smearing effect near the boundary. This smearing has a detrimental effect on the accuracy of the numerical scheme. The other approach in IBM is discrete forcing (sharp) methods [10]. In the discrete IBM, the numerical discretization near the immersed boundary is modified such as to account directly for the presence of the boundary, so that the interface remains “sharp“. The “ghost cell“ in the finite-difference method and the “cut-cell“ within the finite volume framework fall into this category.

In the ghost cell method, boundary conditions at the immersed boundary are enforced through ghost cells. The ghost cell refers to a fictitious cell that lies inside the solid region and whose value is extrapolated from the boundary condition at the immersed boundary and the surrounding fluid points (in the fluid region). While the approach is well suited for achieving second-order accuracy, extension to higher-order formulations is problematic. Higher-order formulations require large interpolation stencils which could lead to robustness issues. For example, for a fourth-order accurate scheme in three dimensions, at least 35 points are required to evaluate the ghost cell [11]. Hence, most of the existing immersed-boundary ghost cell methods are of second-order accuracy [12, 13, 14]. On the other hand, being based on the finite volume approach, the cut-cell method is designed to provide better conservation properties, especially for cells at the immersed boundary [15, 16]. Grid cells cut by the immersed boundary, whose cell centers are inside the fluid region, are reshaped by discarding the portion of these cells that lies in the solid. The disadvantages of this method are related with the process of cutting the cells. The reshaping may in some cases result in very small grid cells with an adverse impact on numerical stability. Hence, the cut cells near boundaries must be adjusted, modified and/or merged. In case of very complex geometries, this process may even fail to preserve the geometrical representation of bodies. The extension of this approach to 3D problems becomes very complicated and a non-trivial task. Similar to the ghost cell method, most of the cut-cell based IMB methods are second-order accurate such as the one adopted by Seo and Mittal [17] to solve the pressure Poisson equation as part of the solution of Navier-Stokes equations.

Standard finite-difference schemes fail when applied to non-smooth functions because the underlying Taylor expansions upon which they are based are invalid. To avoid this problem, jump correction terms need to be added to the finite difference schemes at jump discontinuities for the function value and its derivatives. Based on this key idea, LeVeque & Li [1] developed a second-order accurate sharp immersed interface method to solve elliptic problems with discontinuous and non-smooth solutions. The advantages of this technique are that boundary conditions are imposed directly at the location of the boundary and high-order local accuracy can be achieved around the immersed boundary. An important application of sharp IIM is the treatment of problems defined in an irregular domain where the solution inside or the outside the interface is not of interest and is trivial [18, 19]. Apart from the IBM and IIM, other approaches have been proposed in the literature to solve 3D elliptic equations with irregular boundaries such as the second-order accurate Shortley–Weller embedded finite-difference method for the solution of the 3D Poisson equation [20] and high-order matched interface and boundary method for solving elliptic equations with discontinuous coefficients and non-smooth interfaces [21].

Eq. (1) is an elliptic partial differential equation with a broad range of applications in electromagnetism, geophysics, astrophysics and fluid mechanics. Therefore, there is great interest to develop highly-accurate and computationally efficient numerical methods for the numerical solution of Eq. (1) for simple and complex geometries. Building on our previous research [22], the main goal of this paper is to present a solver for the 3D Poisson equation in a domain with immersed boundaries on a nonuniform grid which combines high accuracy and high efficiency. The objective of the present work is twofold: (1) Developing a uniformly fourth-order-accurate finite-difference algorithm on an irregularly shaped boundary, and (2) designing an efficient and cost-effective iterative solver which easily and efficiently accommodates the irregular immersed boundaries. The paper is organized as follows: The discretization of Eq. (1) on a nonuniform grid is presented in Section 2. It contains the numerical procedure to construct a fourth-order accurate compact difference stencils for regular and irregular grid points, as well as a formal proof for the order of accuracy. A spectral analysis of the resulting coefficient matrices for simple and irregular domains and the solution strategy of the discretized equations is explained in Section 3. Then, in Section 4, the proposed method is validated for several test cases, which demonstrate the high efficiency and confirm the fourth-order convergence. A summary and conclusions are provided in Section 5.

2 Discretization

Eq. (1) is solved in a cubic domain $\Omega$ defined on $[a_{1},a_{2}]\times[b_{1},b_{2}]\times[c_{1},c_{2}]$ . The domain is divided into $nx\times ny\times nz$ uniform/nonuniform cells by the points $a_{1}=x_{1}<x_{2}<\cdots<x_{nx}=a_{2}$ , $b_{1}=y_{1}<y_{2}<\cdots<y_{ny}=b_{2}$ and $c_{1}=z_{1}<z_{2}<\cdots<z_{nz}=c_{2}$ . The discretization of Eq. (1) for regular and irregular grid points are discussed in detail in this section. A grid point is said to be regular if all the 26 neighboring grid points are outside the immersed body, otherwise, it is defined as an irregular point. A grid point is called a solid point if it lies inside an immersed body. Furthermore, the locations where the immersed/irregular boundary intersects with the grid are called immersed/irregular boundary intersection (IBI) points. The IBI points are the locations where the boundary conditions can be enforced. Regular, irregular and IBI points are illustrated in Fig. 1.

2.1 Fourth-order compact difference scheme on nonuniform grids

In this section, the discretization for regular grid points away from the immersed boundary is discussed, so that the compact difference scheme is well-defined and valid. The discretization of Eq. (1) is based on three one-dimensional, fourth-order compact finite-difference schemes for second derivatives in $x$ , $y$ and $z$ :

[TABLE]

where the finite difference (FD) operators are given by

[TABLE]

and

[TABLE]

Here $u_{xx}$ , $u_{yy}$ and $u_{zz}$ represent numerical approximations to the second partial derivatives in $x$ , $y$ and $z$ directions, respectively. Coefficients for compact FD operators in the $x-$ direction in Eqs. (3) and (4), obtained by matching the coefficients in the Taylor expansion about $u_{i,j,j}$ in the $x-$ direction, are given by

[TABLE]

where $dx_{f}=x_{i+1}-x_{i}$ and $dx_{b}=x_{i}-x_{i-1}$ . The coefficients for compact FD operators in the $y-$ and $z-$ directions in Eqs. (3) and (4) are the same except that $x$ is replaced with $y$ and $z$ , respectively. Combining Eq. (2) at three consecutive $x-$ , $y-$ and $z-$ locations centered at point ( $i,j,k$ ) leads to

[TABLE]

Applying the FD operators in Eq. (6) and using Eq. (1) leads to a 27-points, fourth-order compact scheme at the regular grid point ( $i,j,k$ ) inside the computational domain as follows:

[TABLE]

where the RHS is obtained using the following relation

[TABLE]

The coefficients for LHS in Eq. (7) and the coefficents for RHS in Eq. (8) are given by

[TABLE]

In Eq. (9), $\alpha$ is 0, 1 and 2 with corresponding $\beta=a$ , $b$ and $c$ , respectively.

2.2 Treatment of irregular grid points: Fourth-order sharp immersed interface method

In this section, the method to determine the coefficients of the compact scheme stencil at an irregular grid point is presented. Our method falls into the sharp interface category. However, it distinguishes itself from other IIMs

in that jump corrections are no longer required. The key aspect of the new method is to modify and adjust the compact finite-difference operators, Eqs. (3)-(4), when they intersect an immersed boundary to obtain a solution that is sharp across the interface while keeping the fourth-order global accuracy. In particular, for an irregular grid point located at ( $i,j,k$ ), modified FD operators need to be employed along each grid line in the range $x_{i-1}\leq x\leq x_{i+1}$ , $y_{j-1}\leq y\leq y_{j+1}$ and $z_{k-1}\leq z\leq z_{k+1}$ , if there is an intersection with the immersed boundary. Otherwise, the standard FD operators can be used. For example, in the $y-$ direction, in the range $y_{j-1}\leq y\leq y_{j+1}$ , check if the grid lines passing through $x=x_{i-1}$ , $x=x_{i}$ and $x=x_{i+1}$ in the planes $z=z_{k-1}$ , $z=z_{k}$ and $z=z_{k+1}$ , have any intersection with the immersed boundary. If yes, modified FD operators have to be used along that line, otherwise, the standard FD operators in the $y-$ direction given by Eqs. (3)-(4) are valid and can therefore be used. Similar approach is employed in the $y-$ and $z-$ directions. The modified compact FD stencil will be built based on the combination of the standard and modified FD operators.

For illustration, we consider the $27-$ point 3D stencil centered at irregular grid point ( $i,j,k$ ) as shown in Fig. 2(a). In the $x-$ direction, the lines $y=y_{j-1}$ , $y=y_{j}$ , and $y=y_{j+1}$ in the plane $z=z_{k-1}$ , and the lines $y=y_{j}$ and $y=y_{j+1}$ in the plane $z=z_{k}$ do intersect the immersed boundary. In the $y-$ direction, the line passing through $x=x_{i+1}$ in the plane $z=z_{k}$ crosses the immersed boundary. Finally, the lines $y=y_{j-1}$ , $y=y_{j}$ , and $y=y_{j+1}$ in the plane $x=x_{i+1}$ do intersect the immersed boundary in the $z-$ direction. The FD operators along these lines need to be adjusted to take into account the immersed boundary while maintaining the formal fourth-order accuracy. Eq. (6) is rewritten as

[TABLE]

In the above equation, the FD operators with the overbar are modified, while the other operators are the same as those given in Eqs. (3)-(4). A key point is that the coefficients used in the modified LHS FD operators ( $\overline{L_{xx}},\overline{L_{yy}},\overline{L_{zz}}$ ) in Eq. (10) must be kept the same as those used in Eq. (3). However, the coefficients corresponding to the grid points inside the immersed body have to be dropped. Therefore, to maintain the formal fourth order accuracy, additional grid points are needed to determine the coefficients for the modified RHS FD operators ( $\overline{R_{xx}},\overline{R_{yy}},\overline{R_{zz}}$ ) in Eq. (10). Furthermore, the line passing through $z=z_{k-1}$ in the plane $x=x_{i+1}$ in the $y-$ direction falls into the subdomain $\Omega^{-}$ , hence, the corresponding FD operators are dropped (struck through in Eq. 10).

In the $x-$ direction along the planes $z=z_{k-1}$ and $z=z_{k}$ (see Fig. 2d), $(\overline{L_{xx}u_{xx}})_{J,K}=(\overline{R_{xx}u})_{J,K}$ where

[TABLE]

In Eqs. (11)-(12), $J$ is $j-1$ , $j$ , and $j+1$ for $K=k-1$ with corresponding $\alpha=1$ , $2$ and $3$ , respectively. For $K=k$ , $J$ is $j$ , and $j+1$ with $\alpha=4$ and $5$ . Furthermoe, $u_{{IBI}_{\alpha}}$ is the known function value at the body intercept location $IBI_{\alpha}$ . It can be seen that the term $c^{x}_{i}~{}u_{{xx}_{i+1,J,K}}$ is not included in Eq. (11) since it is inside the immersed body. For the right-hand side operator, Eq. (12), we use three additional grid points to keep the fourth-order accuracy (see Fig. 2d). Taylor series expansions about $u_{i,J,K}$ are used to find the coefficients by solving the following system of equations:

[TABLE]

In the above equation, the coefficient matrix $\mathcal{M}$ is defined as

[TABLE]

where $h_{0}=x_{i-1}-x_{i}$ , $h_{1}=x_{i-2}-x_{i}$ , $h_{2}=x_{i-3}-x_{i}$ , $h_{3}=x_{i-4}-x_{i}$ and $h_{b}=x_{{IBI}_{\alpha}}-x_{i}$ . Note that the coefficients for the RHS in Eq. (13) are given by Eq. (5). For the stencil passing through the line $x=x_{i+1}$ along $z=z_{k}$ as shown in Fig. 2(e), $(\overline{L_{yy}u_{yy}})_{i+1,k}=(\overline{R_{yy}u})_{i+1,k}$ where the modified FD operators are given by

[TABLE]

$u_{{IBI}_{6}}$ denoting the known boundary value at IBI location 6, $IBI_{6}$ . One should note that the Taylor expansion about $u_{i+1,j-1,k}$ is used to find the coefficients in Eq. (16) as follows

[TABLE]

where in the matrix $\mathcal{M}$ , given by Eq. (14), $h_{0}=y_{j-2}-y_{j-1}$ , $h_{1}=y_{j-3}-y_{j-1}$ , $h_{2}=y_{j-4}-y_{j-1}$ , $h_{3}=y_{j-5}-y_{j-1}$ and $h_{b}=y_{{IBI}_{6}}-y_{j-1}$ . Finally, for the stencil along the plane $x=x_{i+1}$ , the FD scheme takes the form $(\overline{L_{zz}u_{zz}})_{i+1,J}=(\overline{R_{zz}u})_{i+1,J}$

[TABLE]

The above equations hold for $J=j-1$ (see Fig. 2c). For $J=j$ and $J=j+1$ , the modified FD operators are given by

[TABLE]

with $\alpha=2$ and $\alpha=3$ for $J=j$ and $J=j+1$ , respectively (see Fig. 2b). Matching the Taylor series coefficients about $u_{i+1,j-1,k}$ in the $z-$ direction, the coefficients in Eq. (19) can be found from

[TABLE]

where $h_{0}=z_{k+1}-z_{k}$ , $h_{1}=z_{k+2}-z_{k}$ , $h_{2}=z_{k+3}-z_{k}$ , $h_{3}=z_{k+4}-z_{k}$ and $h_{b}=z_{{IBI}_{7}}-z_{k}$ are used in the matrix $\mathcal{M}$ . In Eq. (21), the coefficients are obtained by matching the Taylor coefficients about $u_{i+1,J,k+1}$

[TABLE]

where $h_{0}=z_{k+2}-z_{k+1}$ , $h_{1}=z_{k+3}-z_{k+1}$ , $h_{2}=z_{k+4}-z_{k+1}$ , $h_{3}=z_{k+5}-z_{k+1}$ and $h_{b}=z_{{IBI}_{\alpha}}-z_{k+1}$ are used in the coefficient matrix $\mathcal{M}$ in Eq. (14). Applying the above equations to Eq. (10), we get the modified compact scheme stencil at the irregular grid point ( $i,j,k$ )

[TABLE]

with the modified coefficients

[TABLE]

The modified RHS term $\overline{Q_{i,j,k}}$ is the same as the one in Eq. (8) except that all the coefficients corresponding to the grid points that lie inside the immersed body are set to zero. In Eq. (24), $\mathcal{B}_{i,j,k}$ includes the known function values at the immersed boundary intercept points,

[TABLE]

and $\mathcal{C}_{i,j,k}$ contains the additional points used to keep the fourth-order formal accuracy

[TABLE]

It is worth noting that the coefficients without the overbar in Eq. (24) are the same as the coefficients obtained for the regular grid points in Eq. (9). It should be noted that the coefficients for the solid grid points are not used in Eq. (24) since the function value is zero at those grid points (struck through in Eq. 24).

2.3 Special case: multiple intersections with immersed boundary

In the previous section, the immersed boundary has only one intersection with each grid line in the range of $x_{i-1}\leq x\leq x_{i+1}$ , $y_{j-1}\leq y\leq y_{j+1}$ and $z_{k-1}\leq z\leq z_{k+1}$ . In many practical applications, however, the immersed boundary $\Gamma$ can cross the grid lines more than once as shown in Fig. 3, for example. In this case, the compact finite-difference stencil centered at the irregular grid point ( $i,j,k$ ) takes the following form:

[TABLE]

According to Fig. 3(a), the line passing though $z=z_{k}$ in the plane $x=x_{i+1}$ in the $y-$ direction falls into the subdomain $\Omega^{-}$ , hence, the corresponding FD operators are dropped (struck through in Eq. 28). In the $x-$ direction, the plane passing through $z=z_{k}$ crosses the immersed boundary at the $IBI_{7}$ , $IBI_{8}$ and $IBI_{9}$ . Along this plane, $(\overline{L_{xx}u_{xx}}-\overline{R_{xx}u})_{J,k}=0$ where

[TABLE]

In Eqs. (29)-(30), $J$ is $j-1$ , $j$ and $j+1$ with corresponding $\alpha=1$ , $2$ and $3$ , respectively. The coefficients in Eq. (30) can be found from matching the Taylor series coefficients about $u_{i,J,k}$ in the $x-$ direction (see Eq. 13). In the $z-$ direction, however, the lines $y=y_{j-1}$ , $y=y_{j}$ and $y=y_{j+1}$ in the plane $x=x_{i+1}$ intersect the immersed boundary at more than one location as shown in Fig. 3(b), namely $IBI_{1}$ & $IBI_{4}$ along $y=y_{j-1}$ , $IBI_{2}$ & $IBI_{5}$ along $y=y_{j}$ and $IBI_{3}$ & $IBI_{6}$ along $y=y_{j+1}$ . Therefore, special care is needed here to modify FD operators along these lines.

For the stencil passing through the plane $x=x_{i+1}$ in the $z-$ direction, $(\widehat{L_{zz}u_{zz}})_{i+1,J}=(\widehat{R_{zz}u})_{i+1,J}$ where $(\widehat{L_{zz}u_{zz}})_{i+1,J}=(\overline{L_{zz}u_{zz}})_{i+1,J}+(\overline{\overline{L_{zz}u_{zz}}})_{i+1,J}$ and $(\widehat{R_{zz}u})_{i+1,J}=(\overline{R_{zz}u})_{i+1,J}+(\overline{\overline{R_{zz}u}})_{i+1,J}$ , with $J=j-1$ , $J=j$ and $J=j+1$ . The modified FD operators are given by

[TABLE]

and

[TABLE]

with $\alpha=1$ , $2$ and $3$ for $J=j-1$ , $j$ and $j+1$ , respectively. Here, the Taylor series expansions in the $z-$ direction about $u_{i+1,J,k+1}$ and $u_{i+1,J,k-1}$ are used to find the coefficients in Eqs. (32) and (34), respectively. Using the original FD operators provided by Eqs. (3)-(4) for the unmodified operators in Eq. (28) and applying the modified FD operators given by Eqs. (29)-(34) to Eq. (28), the modified compact scheme stencil at the irregular grid point ( $i,j,k$ ) can be obtained,

[TABLE]

The modified coefficients are

[TABLE]

The correction term, $\mathcal{C}_{i,j,k}$ , is given by

[TABLE]

The known function values at the points of intersection with the immersed boundary are included in $\mathcal{B}_{i,j,k}$ ,

[TABLE]

2.4 Order of original and modified compact finite-difference stencils

In this section, the order of accuracy of the compact FD stencils for regular and irregular grid points developed in Sections 2.1 $-$ 2.3, is formally verified. It is worth noting that all the FD operators employed to construct the original and modified compact FD stencils, are based on the desired fourth-order compact finite-difference approximations. Using the fact that the eigenfunctions of the exact Laplace operator are known to be $\phi_{k,l,m}=cos(kx)~{}cos(ly)~{}cos(mz)$ with corresponding eigenvalues $\lambda=-(k^{2}+l^{2}+m^{2})$ , one can estimate the order of discrete FD operators where eigenvalues change according to the order of the accuracy of finite-difference approximations [23, 24]. In the derivation of the compact FD stencils, one operator acts on the field $u$ , and one operator operates on the source term $f$ , i.e. $A_{h}~{}\boldsymbol{u}=Q_{h}~{}\boldsymbol{f}$ , where $A_{h}$ and $Q_{h}$ are discrete operators, see Eqs. (7)-(8) for example. Therefore, a generalized eigenvalue problem, $A_{h}~{}\boldsymbol{u}=\lambda~{}Q_{h}~{}\boldsymbol{u}$ , has to be solved. As pointed out by Sutmann and Steffen [23], both $A_{h}$ and $Q_{h}$ have the desired property that eigenfunctions are sampled continuous eigenvectors. As a result, only the eigenvalues need to be considered. For practical purposes one may insert the eigenfunctions of the exact operator and estimate the order of FD approximations in $h=max(dx,dy,dz)$ , i.e.,

[TABLE]

where $n$ is the order of the difference approximation. Without loss of generality, it is assumed that the grid is uniform in each direction, however, the grid spacing in each direction could be different, i.e. $dx\neq dy\neq dz$ . Furthermore, the locations of the points of intersection with the immersed boundary are estimated according to Figs. (2)-(3). Applying the FD coefficients developed in Sections 2.1 $-$ 2.3, provides the order of the compact FD stencils for a regular grid point, Eq. (7), and irregular grid points, Eqs. (24) and (35), as follows:

[TABLE]

where H.O.T refers to higher order terms. For all the FD stencils developed for the regular and irregular grid points, the expected fourth-order accuracy is recovered.

3 Solution strategy

The discretization of the 3D Poisson equation on a uniform/nonuniform grid with $nx\times ny\times nz$ points leads to a set of $(nx-2)\times(ny-2)\times(nz-2)$ linear equations for all unknown function $u_{i,j,k}$ inside the computational domain, $2\leq i\leq nx-1$ , $2\leq j\leq ny-1$ and $2\leq k\leq nz-1$ . The system of linear equations can be represented in the matrix-vector form as

[TABLE]

where $\mathcal{A}$ is the coefficient matrix (see Eq. 7), $\boldsymbol{q}$ is the source term (see Eq. 8) and $\boldsymbol{u}$ is the unknown vector containing all $u_{i,j,k}$ in the interior of the computational domain. The coefficient matrix $\mathcal{A}$ has a 27-diagonal structure for a simple domain without immersed boundary. It is worth noting that the matrix $\mathcal{A}$ is not symmetric in general. The only situation where the matrix $\mathcal{A}$ is symmetric is for a domain with uniform grid spacing.

The discretization of Eq. (1) for all regular and irregular grid points inside the computation domain with an immersed boundary leads to a linear system

[TABLE]

where $\boldsymbol{b}$ includes the known function values at the points of intersection with the immersed boundary (see Eqs. 26 and 38). The coefficient matrix $\mathcal{B}$ has the same dimension as matrix $\mathcal{A}$ , i.e. $N\times N$ where $N=(nx-2)\times(ny-2)\times(nz-2)$ . However, as discussed in the previous sections, the finite-difference schemes are modified when the 27-point stencil intersects with the immersed boundary. As a result, this introduces additional points to the 27-point compact discretization that are solution dependent, so that the coefficient matrix $\mathcal{B}$ is no longer 27-diagonal. The matrix $\mathcal{B}$ is non-symmetric even for a domain with uniform grid spacing.

3.1 Spectral analysis: eigenvalue spectra of the coefficient matrices

The eigenvalue spectra of the resulting coefficient matrices for a simple domain and a domain with an immersed boundary is investigated in this section. The numerical stability involving a finite difference approximation for Eq. (1) requires that all eigenvalues ( $\lambda$ ) of $\mathcal{A}$ (or $\mathcal{B}$ in the case with immersed boundaries) satisfies $real(\lambda)<0$ [25]. The eigenvalues of the discrete finite-difference operator will depend upon the geometry and grid resolution of the problem. Therefore, it is necessary to evaluate the effect of the extended stencil due to the presence of an immersed boundary on the eigenvalue spectra of matrix $\mathcal{B}$ . The eigenvalue spectra of the fourth-order-accurate compact difference scheme for two different cases are investigated on a cuboid of dimensions $[0,1]^{3}$ : (i) A simple domain and (ii) a domain with an immersed sphere located at the center of the domain with radius 0.2.

The eigenvalue spectra of the resulting discretization matrix for the two cases when computed for a uniform $21^{3}$ grid are presented in Fig. 4. The eigenvalues are normalized with $dx^{2}$ . For the simple domain, all eigenvalues are real and negative as the coefficient matrix is symmetric. As can be seen the immersed boundary introduces complex eigenvalues, however, all the eigenvalues have negative real components and therefore satisfy the stability condition.

3.2 Preconditioned BiCGSTAB

The fourth-order compact finite-difference discretization of Eq. (1) at all grid points forms a large sparse linear system. In this work, a stabilized biconjugate gradient (BiCGSTAB) method is implemented for solving the linear system, which is efficiently adjusted to account for the sharp immersed boundaries [26]. Algorithm 1 shows a schematic of the BiCGSTAB iterative method for a simple domain (without immersed boundary) with starting guess $\boldsymbol{u}_{0}$ .

As discussed in the previous section, in the presence of immersed boundaries, the coefficient matrix $\mathcal{A}$ is replaced with the matrix $\mathcal{B}$ . The matrix $\mathcal{B}$ can be decomposed into

[TABLE]

where the modified matrix $\overline{\mathcal{A}}$ has the same structure as matrix $\mathcal{A}$ (27-diagonal matrix). In fact, all the rows corresponding to the regular grid points are the same as the ones in matrix $\mathcal{A}$ , however, some of the coefficients in the rows corresponding to the irregular grid points are modified. In Eq. (43), the matrix $\mathcal{C}$ represents irregular entries caused by the extended stencils at the irregular grid points, see Eqs. (27) and (37). The following changes are incorporated into algorithm 1. We have three matrix vector multiplications in algorithm 1 , lines 1, 5 and 9. Thus, they are modified accordingly $~{}~{}1:~{}~{}\boldsymbol{r}_{0}=\boldsymbol{q}-\overline{\mathcal{A}}\boldsymbol{u}_{0}+\big{<}\boldsymbol{b}-\mathcal{C}\boldsymbol{u}_{0}\big{>}$ $~{}~{}5:~{}~{}~{}\hskip 1.87772pt\boldsymbol{v}=\overline{\mathcal{A}}\boldsymbol{\hat{p}}+\big{<}\mathcal{C}\boldsymbol{\hat{p}}\big{>}$ $~{}~{}9:~{}~{}~{}\hskip 2.98741pt\boldsymbol{t}=\overline{\mathcal{A}}\boldsymbol{\hat{s}}+\big{<}\mathcal{C}\boldsymbol{\hat{s}}\big{>}$ In the modified lines 1, 5 and 9 of algorithm 1, the $<~{}>$ corresponds to the additional operation due to the presence of immersed boundaries. These operations are carried out only for irregular grid points.

3.3 Preconditioner based on second-order finite-difference approximation

Although BiCGSTAB is one of the most effective Krylov subspace methods, it can suffer from slow rate of convergence or even total breakdown. Preconditioning is a key ingredient for the improvement of both the efficiency and the robustness of Krylov subspace methods. In this section, the development of a preconditioning technique whose efficiency and cost are independent of the presence or not of an immersed boundary is presented. In algorithm 1, $\mathcal{K}$ is a preconditioner matrix which approximates $\mathcal{A}$ in some sense (or $\overline{\mathcal{A}}$ in the case with immersed boundaries). From a practical point of view, an efficient preconditioner should be low-cost to construct and the preconditioned system $\mathcal{K}\boldsymbol{\phi}=\boldsymbol{\rho}$ should be inexpensive to solve (see lines 4 and 8 in algorithm 1).

The matrix $\mathcal{K}$ is built based on the discretization of $\nabla^{2}\phi=\rho$ with zero Dirichlet boundary condition on the boundaries. The second partial derivatives are approximated using a second-order standard (non-compact) finite-difference schemes in $x$ , $y$ and $z$ :

[TABLE]

where the coefficients in the $x-$ direction are given by

[TABLE]

where $dx_{f}=x_{i+1}-x_{i}$ and $dx_{b}=x_{i}-x_{i-1}$ . The coefficients in $y-$ and $z-$ directions in Eq. (44) are the same except that we replace $x$ with $y$ and $z$ , respectively. Incorporating Eq. (44) into $\nabla^{2}\phi=\rho$ leads to a seven-point, second-order scheme at the regular grid point ( $i,j,k$ ) as shown in Fig. 5(a)

[TABLE]

Now consider the seven-point stencil centered at the irregular grid point located at ( $i,j,k$ ) which intersects with an immersed boundary as illustrated in Fig. 5(b). We use the same irregular grid point and immersed boundary as shown in Fig. 2. The immersed boundary intercepted point, $IBI_{1}$ is located at $x_{i}<x_{{IBI}_{1}}<x_{i+1}$ . As a result, the second derivative in the $x-$ direction needs to be adjusted to take into account the immersed boundary. Therefore

[TABLE]

with the modified coefficients

[TABLE]

where $\overline{dx_{f}}=x_{{IBI}_{1}}-x_{i}$ . It should be noted that $\phi_{{IBI}_{1}}=0$ . This leads to the modified stencil at the irregular grid point centered at ( $i,j,k$ )

[TABLE]

Note that the grid point located at ( $i+1,j,k$ ) falls into $\Omega^{-}$ , hence, it is dropped in the above equation (struck through in Eq. 49). It is important to note that the preconditioner matrix $\mathcal{K}$ has a seven-diagonal structure

with or without the presence of an immersed boundary. It is thus a convenient algorithm that easily accommodates immersed boundaries, yet the cost to build and solve the preconditioned system is independent of the complexity of the geometry and the presence or not of an immersed boundary. In the present work, this preconditioned system is solved iteratively using the modified strongly implicit (MSI) procedure [27]. The MSI method has very good convergence properties was shown to outperform Stone’s strongly implicit procedure [stone], incomplete lower-upper decomposition (ILU) and successive over-relaxation methods.

3.4 Floating-point operation counts

Another important aspect of the solution of Poisson equation with a high-order discretization method is the extra floating-point operation (FLOP) of the algorithm components due to the presence of an immersed body. The preconditioned BiCGSTAB algorithm is built from several basic components: Sparse matrix vector products (SPMV) such as $\mathcal{A}\boldsymbol{\hat{v}}$ , inner products of two vectors (DOT-PROD) like $(\boldsymbol{r}_{0},\boldsymbol{r}_{i})$ and combined scalar vector multiplication and vector addition (SVPV) for example $\boldsymbol{s}-\omega\boldsymbol{t}$ . Moreover, algorithm 1 contains two preconditioned systems (PCS) that need to be solved. Table 1 lists the FLOP counts for the preconditioned BiCGSTAB algorithm per outer iteration for a simple domain and a domain with immersed boundary.

In Table 1, $iter_{prec}$ refers to the number of inner iteration for solving the preconditioned system, $N$ is the total number of grid points and $N_{irr}$ is the number of irregular grid points. Therefore, the extra FLOP counts due to the presence of an immersed body as a percentage of FLOPs for a simple domain is

[TABLE]

It should be noted that the FLOP counts for $<~{}>$ , which corresponds to the operations carried out for irregular grid points come from the fact that there are nine FD operators for each direction and for each FD operator, there is a maximum of four additional grid points, resulting in 8 FLOPs (see Eq. 27 for example). Therefore, the maximum FLOP count for each irregular grid point will be $3\times 9\times 8=216$ FLOPs. For $iter_{prec}=4$ which is typical for

a preconditioned BiCGSTAB solver, we could have $FLOP_{extra}(\%)\approx N_{irr}/N\times 100$ which in most cases is negligible. It is worth mentioning that a different iterative solver for the preconditioned system ( $\mathcal{K}\boldsymbol{\phi}=\boldsymbol{\rho}$ ) will result in a different FLOP count (PCS). However, the number of FLOPs for PCS will be the same regardless of the presence or absence of an immersed body.

4 Numerical experiments

In this section, numerical experiments are conducted to demonstrate the accuracy and efficiency of the present method for solving Eq. (1) on uniform and nonuniform grids with immersed boundaries. All problems are solved on the cuboid domain $\Omega$ with the same number of grid points in the $x$ , $y$ and $z$ directions. The order of accuracy is computed using the definition

[TABLE]

where $||E_{n_{1}}||_{\infty}$ and $||E_{n_{2}}||_{\infty}$ are the infinity norms of the errors between the analytical, $u^{(ex)}$ , and the numerical, $u$ , solutions for two different grids with $n_{1}$ and $n_{2}$ points in each direction. For the preconditioner solver, eight inner iterations are performed for all grid sizes. The initial guess is the zero vector.

4.1 Problem 1

To demonstrate the accuracy of the new method and the efficiency of the preconditioned BiCGSTAB solver, solving the following Poisson equation with Dirichlet boundary conditions is considered:

[TABLE]

where $\omega=2\pi$ . Eq. (52) is solved on a cuboid of dimensions $[0,1]^{3}$ with a uniform grid in the $x$ , $y$ and $z$ directions for two cases: (i) simple domain without immersed boundary and (ii) a domain with one immersed toroidal surface with an outer radius of 0.35 and an inner radius of 0.05 plus an immersed sphere with radius 0.08 (see Fig. 6, bottom plots). Computed solutions and the corresponding errors, $e=|u^{(ex)}-u|$ are shown in Fig. 6 on a $129^{3}$ grid. The reader’s attention is drawn to the sharp interface of the obtained solution, and the relative smoothness of the error distribution in the vicinity of the immersed boundary. The solution and the corresponding error is computed for several grids. The errors in the infinity-norm, $||E_{N}||_{\infty}$ , are plotted as a function of grid size in Fig. 7(a). The results confirm

the fourth-order accuracy of the new method for both the simple domain and, most importantly, also for the domain with immersed boundaries. No loss of accuracy for the domain with immersed boundary is observed.

In the following, the convergence characteristics and the $FLOP_{extra}$ of the new method for solving Eq. (52) is compared between the simple domain and domain with immersed bodies. Figure 7(b) displays the convergence rates for various grid sizes for the two cases described above. It shows the infinity-norm of the error after each outer iteration. The convergence behaviour is very similar for the two cases for the different grid sizes. With a proper design of the discretization operator and the preconditioning strategy, the new solution strategy is proven to be equally efficient for domains with immersed boundaries and for simple domains (without immersed boundaries). Table 2 presents the extra FLOP counts for the problem with immersed boundaries. For the problem sizes of interest, Table 2 shows that the $FLOP_{extra}$ due to the presence of the immersed boundary is negligible for the preconditioned BiCGSTAB solver. For the largest grid size investigated, the $FLOP_{extra}$ is about 0.11%.

4.2 Problem 2

Consider the following equation with Dirichlet boundary conditions:

[TABLE]

with $\omega=2\pi$ and $\sigma=-10$ . In this example, the computational domain is set to $[0,1]\times[−0.3,0.3]\times[0,1]$ and the immersed body is a finite tapered wing with a NACA 5514 airfoil section. Figure 8 illustrates the parameters that define the wing shape which are the wing span $b=0.7$ , the wing root chord $c_{r}=0.6$ , and the tip chord $c_{t}=0.35$ . Numerical investigations of the external flow over a finite wing represent one of the difficult challenges in the field of computational fluid dynamics. One of the main difficulties lies in the generation of a high-quality mesh around the wing tip. Sharp immerse interface method are therefore very useful for that type of problems.

Eq. (53) is solved on a uniform grid in the $x$ and $z$ directions while a non-uniform grid is used in the $y$ direction, for which the stretching function is given by

[TABLE]

The parameter $s$ is defined as

[TABLE]

where $L=0.3$ and $H=0.6$ . The parameter $\beta$ is the clustering parameter which controls the nonuniformity of the grid around $y=0$ . It is set to $\beta=3$ which puts more grid points around $y=0$ .

Figure 9 displays the numerical solution to Eq. (53) and the corresponding error for a grid with $129^{3}$ points. The solutions are sharp across the interfaces and the local error has a smooth distribution in the vicinity of the immersed boundary. Table 3 shows the numerical errors based on the infinity-norm of the current method for different sets of grids which demonstrates that the method is fourth-order-accurate. The extra FLOP as a result of the immersed

body in this problem is presented in Table 4, which varies from 0.11(%) to 1.68(%) depending on the grid size.

4.3 Problem 3

Next, we consider solving the following Poisson equation on a uniform grid

[TABLE]

where $r=(x-1)^{2}+(y-1)^{2}+(z-1)^{2}$ and $\omega=2$ . Eq. (56) is solved in the cuboid domain with dimensions $[0,2]^{3}$ in presence of four immersed bodies as illustrated in Fig. 10. In this example, the immersed bodies imitate tall buildings encountered in wind engineering, where wind-induced load effects on buildings with different heights and proximities are regularly investigated. Fig. 10 shows the numerical solution computed by the preconditioned BiCGSTAB method using $129^{3}$ grid points as well as the contour lines of the local errors.

The results of the error study are presented in Table 5 which demonstrates that the method is fourth-order-accurate based on the error measured in the infinity-norm.

4.4 Problem 4

For the final test case, the following equation with Dirichlet boundary condition on a nonuniform grid with irregular boundaries is solved:

[TABLE]

where $r=(x-0.5)^{2}+(y-0.5)^{2}+(z-0.5)^{2}$ and $\sigma=-10$ . The domain $\Omega^{+}$ is the region inside the boundary $\Gamma$ as shown in Fig. 11. In contrast to the other test problems, Eq. (57) is solved inside the immersed boundary. This test problem is designed to show the accuracy and efficiency of the new method for simulations of internal flows such as the flow through a channel or stenosed aorta.

Grid points are distributed uniformly in the $x$ direction for $0\leq x\leq 1$ while the stretching function for the nonuniform grids in the $y-$ and $z-$ directions is given by

[TABLE]

where $\beta=2$ and $\beta=4$ for Eqs. (58) and (59), respectively. The parameter $s$ given by Eq. (55) is computed with $L=0.14$ & $H=0.5$ for Eq. (58) and with $L=0.15$ & $H=0.3$ for Eq. (59). This test case could be a

model of stenosis pipe flow with 66% reduction of the area. Shown in Fig. 11 is a curved pipe centered at $x_{c}=0.5$ , $y_{c}=0.5$ and $z_{c}=0.5$ with an inner diameter of $D_{i}=0.24$ . The inner diameter of the reduced area is defined as

[TABLE]

where $\gamma=0.05$ is a coefficient to apply a 66% area reduction, $L_{x}(=0.35)$ denotes the length of the stenosis and its center is located at $x_{0}=0.75$ . The numerical solutions and error obtained for Eq. (57) are illustrated in Fig. 12 for a grid size of $129^{3}$ . Different grid sizes ranging from $33^{3}$ to $513^{3}$ were used to evaluate the order of accuracy summarized in Table 6. It can be observed that the new method preserves its fourth-order accuracy for this problem as well. Furthermore, the extra FLOP counts per outer iteration varies from 6.3% for the smallest grid size to 0.46% for the largest grid size.

5 Conclusion

In this paper, the development of a sharp-interface and high-order method to solve the three-dimensional Poisson equation with Dirichlet boundary conditions on nonuniform grids with immersed boundaries is presented. Fourth-order compact difference schemes are used to discretize the Poisson equation on the nonuniform grid resulting in 27-point FD regular stencils. The sharp representation of the interface of the immersed body is accomplished by modifying the finite difference stencils at irregular grid points (near the immersed boundary). In particular, additional grid points are added to the standard 27-point FD stencils to obtain a sharp solution across the interface while retaining the fourth-order formal accuracy.

A preconditioned BiCGSTAB method is designed to solve the large sparse algebraic system that arose from the discretization of the Poisson equation. A second-order standard (non-compact) finite-difference scheme is employed to discretize the preconditioned system leading to a seven-diagonal matrix. The extra floating-point operation ( $FLOP_{extra}$ ) count associated with the presence of immersed boundaries versus the case with simple domain ( $N$ grid points) showed that $FLOP_{extra}\approx O(N_{irr}/N)$ , where $N_{irr}$ is the number of irregular grid points. Two factors contribute to the cost-effectiveness and efficiency of the new method: (1) the preconditioner cost which is independent of the presence of an immersed boundary, and (2) the compact nature of the finite-difference discretization, which minimizes the additional operations per irregular grid point.

The accuracy and the efficiency of the new method was demonstrated through verification using problems with analytical solutions for domains with and without immersed boundaries on uniform and nonuniform grids. The numerical results demonstrate that the solutions are sharp across the interface and that the new method is fourth-order-accurate in the maximum norm, including the irregular grid points near the interface of the immersed boundary. For the test cases investigated, the solution technique was found to be very efficient in the sense that the accuracy and the convergence behavior do not deteriorate with the presence of irregular boundaries while the $FLOP_{extra}$ per outer iteration remains insignificant.

Acknowledgements

This work was supported in part by the Air Force Office of Scientific Research (AFOSR) under grant number FA9550-14-1-0184 and by National Science Foundation (NSF) under grant number 1805273, with Dr D. Smith and Dr R. Joslin serving as the program manager, respectively.

Bibliography27

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Leveque and Li [1994] R. J. Leveque, Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (1994) 1019–1044.
2Spotz and Carey [1996] W. F. Spotz, G. F. Carey, A high-order compact formulation for the 3D Poisson equation, Num. Meth. for Part. Dif. Eq. 12 (1996) 235–243.
3Wang et al. [2006] J. Wang, W. Zhong, J. Zhang, A general meshsize fourth-order compact difference discretization scheme for 3D Poisson equation, Applied Math. and Comp. 183 (2006) 804–812.
4Sutmann and Steffen [2006 a] G. Sutmann, B. Steffen, High-order compact solvers for the three-dimensional Poisson equation, J. Comput. Appl. Math. 187 (2006 a) 142–170.
5Ge [2010] Y. Ge, Multigrid method and fourth-order compact difference discretization scheme with unequal meshsizes for 3D Poisson equation, J. Comp. Phys. 229 (2010) 6381–6391.
6Ge et al. [2013] Y. Ge, F. Cao, J. Zhang, A transformation-free HOC scheme and multigrid method for solving the 3D Poisson equation on nonuniform grids, J. Comp. Phys. 234 (2013) 199–216.
7Nejat and Ollivier-Gooch [2003] A. Nejat, C. Ollivier-Gooch, A high-order accurate Unstructured GMRES solver for Poisson‘s equation, in: 11th annual conference of the computational fluid dynamics society of Canada, 2003.
8Peskin [1972] C. S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. (1972) 252–271.