# Solving incompressible Navier–Stokes equations on irregular domains and
quadtrees by monolithic approach

Hyuntae Cho, Yesom Park, and Myungjoo Kang

###### Abstract

We present a second-order monolithic method for solving incompressible
Navier–Stokes equations on irregular domains with quadtree grids. A semi-
collocated grid layout is adopted, where velocity variables are located at
cell vertices, and pressure variables are located at cell centers. Compact
finite difference methods with ghost values are used to discretize the
advection and diffusion terms of the velocity. A pressure gradient and
divergence operator on the quadtree that use compact stencils are developed.
Furthermore, the proposed method is extended to cubical domains with octree
grids. Numerical results demonstrate that the method is second-order
convergent in $L^{\infty}$ norms and can handle irregular domains for various
Reynolds numbers.

Keywords Incompressible Navier–Stokes, Finite difference method, Quadtree,
Octree, Level set method

##  1 Introduction

The motion of an incompressible fluid flow is described by the Navier–Stokes
equations:

$\begin{gathered}\mathbf{u}_{t}+\mathbf{u}\cdot\nabla\mathbf{u}=\frac{1}{Re}\bigtriangleup\mathbf{u}-\nabla
p+\mathbf{f},\\\ \nabla\cdot\mathbf{u}=\mathbf{0}.\end{gathered}$ (1)

Here, $\mathbf{u}$ is the fluid velocity, $p$ is the pressure, $\mathbf{f}$ is
the external force such as gravity, and $Re$ is the Reynolds number. The
development of numerical methods for incompressible Navier–Stokes equations is
very important in both science and engineering. Applications of incompressible
flow, such as subsonic aerodynamics, ship hydrodynamics, and numerical weather
prediction, sometimes involve complex geometry and require that small-scale
turbulence and vortices be captured. An efficient way to solve the problem is
to generate a body-fitted mesh and adaptively refine the mesh where small-
scale phenomena are observed. Although various successful studies in recent
decades have used body-fitted meshes, we focus on Cartesian grids in this
study.

Peskin’s immersed boundary method (IBM) [35] is a popular and effective
approach to handling irregular domains on Cartesian grids. However, the
numerical delta function of the IBM lowers the order of convergence near the
boundary. To address this issue, Leveque and Li [25] developed the immersed
interface method (IIM) for elliptic interface problems, which shows second-
order convergence in the $L^{\infty}$ norm. The IBM also motivated ghost fluid
methods (GFMs) [16, 27, 23], which define ghost values at grid points outside
the domain. Although the IIM and GFM were developed to handle discontinuity
across the interface, the core concept of these sharp capturing methods has
been extended to treat various boundary conditions for elliptic equations on
irregular domains [18]. In addition to the IIM and GFM, various numerical
methods are used to treat irregular domains. For example, in the virtual node
method [3, 21], variational formulations on irregular domains embedded in a
Cartesian grid are discretized. Coco and Russo developed ghost point finite
difference methods [13, 14, 12], which adopt a multigrid method to impose
boundary conditions. In a finite volume framework, cut-cell methods are widely
used. Originally proposed by Purvis and Bulkhalter [37], these methods were
used to develop Hodge projections of vector fields [32]. Recent works on
finite volume cut-cell methods include the solution of Poisson’s equation [6]
and the convection–diffusion equation [2] with Robin boundary conditions.

Chorin introduced the projection method [11] to solve incompressible
Navier–Stokes equations on Cartesian grids. It uses Hodge decomposition of
vector fields. The intermediate velocities are first computed and then
projected to force a discrete divergence-free condition by solving Poisson’s
equation for the Hodge variable. Although the method developed by Chorin is
only first-order accurate, Brown et al. [7] analyzed a second-order projection
method. They showed that the boundary conditions of the intermediate
velocities are coupled with the Hodge variable, and thus appropriate boundary
conditions for the intermediate velocities are necessary to obtain second-
order-accurate solutions. For rectangular domains, the boundary conditions for
the intermediate velocities and the Hodge variable can be decoupled by using
the Hodge variable from previous time steps. However, for irregular domains,
it is challenging to impose an accurate boundary condition. Furthermore,
especially when the domain moves in time according to the flow, the
intermediate velocities and Hodge variable should be solved simultaneously in
a single linear system. Thus, the advantage of using the projection method is
lost.

An alternative approach to solving the incompressible Navier–Stokes equations
is to discretize the momentum equation while imposing a divergence-free
condition in a single linear system for the velocities and pressure. This
discretization is called the monolithic approach and results in a saddle-point
linear system. Although this method creates a larger linear system than
projection methods do, the boundary conditions for the velocities are no
longer coupled with other variables. Similarly, the jump conditions of two-
phase incompressible flow can be treated more easily than they can by
projection methods, and recent works [10, 8, 39] have reported success in
obtaining convergence in $L^{\infty}$ norms. For a single-phase flow problem
on irregular domains, a very recent study by Coco [12] solved the equations
using the monolithic approach. The author used a ghost-point finite difference
method to impose boundary conditions and obtained second-order convergence of
the velocity and divergence of the velocity in the $L^{1}$ and $L^{\infty}$
norms. These monolithic methods exhibit convergence with nontrivial analytical
solutions; note, however, that discretized linear systems are rank-deficient
(the pressure is unique up to a constant). In [10, 39], the right-hand side
was projected onto the range space of the linear system by formulating the
kernel of its transpose. By contrast, Coco [12] augmented the approach with an
additional scalar unknown. Although Yoon et al. [43] remarked that finding a
unique solution in the orthogonal complement of the kernel is equivalent to
the augmentation strategy, the range space of the resulting linear system was
not specified by Coco [12].

Staggered marker-and-cell (MAC) grids [20] have been widely used to solve the
incompressible Navier–Stokes equations (1) numerically. In a MAC grid, the
velocities are located on cell faces, whereas the pressures or Hodge variables
are located at cell centers. MAC grids do not exhibit the checkerboard
instability that occurs in a collocated grid, and they enable a simple
discretization for the incompressibility condition. However, there are
obstacles to developing numerical formulations on MAC grids for quad/octree
data structures. The extension of numerical methods designed for uniform grids
to the quad/octree grid is not trivial, and it requires additional effort.
Early works on solving incompressible flow on quad/octree grids focused on the
Euler equations, where the viscosity term in (1) was neglected. Popinet [36]
introduced a second-order projection method on a graded octree grid, that is,
a grid where the size ratio of adjacent cells must be 2:1. This work was later
extended by Losasso et al. [28] to relax the restrictions on adaptive
refinement. Although Popinet [36] discretized the advection term in the
Eulerian approach using an upwind method, Losasso et al. [28] used a semi-
Lagrangian method by extrapolating the velocities located on the cell faces to
the cell vertices. For cases where the viscosity terms are not neglected,
Olshanskii et al. [33] and Guittet et al. [19] developed second-order
projection methods on a graded and non-graded octree MAC grid, respectively.
The method of Olshanskii et al. [33] requires a relatively large stencil for
discretizing the viscous term. By contrast, Guittet et al. [19] developed a
finite volume method that requires a smaller stencil; however, a Voronoi mesh
was generated with respect to the cell faces.

Min and Gibou [29] developed a second-order projection method for
incompressible flow on quad/octree grids other than MAC grids, where the
velocities and Hodge variables are collocated at cell vertices. The main
advantage of the collocated grid on vertices is that second-order finite
differences with compact stencils can be used. As in [19], Hodge decomposition
that guarantees stability was proposed. However, the grid layout in [29]
cannot be used in the monolithic approach because there is no boundary
condition for the pressure, whereas homogeneous Neumann boundary conditions
are imposed for the Hodge variable in the projection method.

In this paper, we introduce a numerical method for solving the incompressible
Navier–Stokes equations on an irregular domain with non-graded quadtree
structure and on a cubical domain with non-graded octree structure. Unlike
many other methods on quad/octree grids, the proposed discretization is based
on the monolithic approach. The velocities are collocated at cell vertices,
and the pressure variable is stored at cell centers, which is also called the
Arakawa B grid layout [1]. Because the Arakawa B grid layout is adopted, the
effective and simple finite differences in [29] on the quadtree can be used.
Furthermore, we propose gradient and divergence operators that map from cell
centers to cell vertices and from cell vertices to cell centers, respectively.
These operators are extended to treat an irregular boundary, and numerical
results demonstrate that the proposed discretization is second-order accurate.
The rest of the paper is structured as follows. Section 2 describes a level
set method of capturing irregular domains and the proposed grid layout on a
quadtree data structure. Section 3 explains the numerical method and
discretization in detail. Numerical experiments and results are presented in
Section 4.

##  2  Computational domain and grid layout

###  2.1 Level set method

We are interested in solving the incompressible Navier–Stokes equations (1) in
an irregular domain $\Omega$ within a rectangular computational domain
$\mathbf{D}=\left[a,b\right]\times\left[c,d\right]$. The interface between the
fluid and the solid is captured by the level-set method [34]. The interface is
described by the zero level set of a level set function
$\phi:\mathbf{D}\to\mathbb{R}$ as
$\Gamma=\left\\{\mathbf{x}\in\mathbf{D}|\phi(\mathbf{x},t)=0\right\\}.$ In
particular, we define the fluid domain as
$\Omega=\left\\{\mathbf{x}\in\mathbf{D}|\phi(\mathbf{x},t)\leq 0\right\\},$
and thus the solid domain is $\mathbf{D}\setminus\Omega$.

Among infinitely many level sets that represent the interface as a zero level
set, a signed distance function has the desirable property that
$|\nabla\phi|=1$, and thus $\phi$ does not exhibit the instability caused by a
steep or shallow gradient. Therefore, the level set function is reinitialized
to become a signed distance function by solving the following eikonal
equation:

$\phi_{\tau}+\text{sgn}\left(\phi^{0}\right)\left(\left\|\nabla\phi\right\|-1\right)=0,$

where sgn denotes the signum function, and $\phi^{0}$ is the initial level set
function. In particular, we adopt the subcell resolution technique of Russo
and Smereka [38] to discretize the equation and the second-order
discretization of Min and Gibou [30]. An advantage of the level set method is
that $\phi$ can be used to compute geometric quantities such as the normal
vector and curvature according to the following formula:

$\mathbf{n}=\frac{\nabla\phi}{\left|\nabla\phi\right|},\quad\kappa=\nabla\cdot\frac{\nabla\phi}{\left|\nabla\phi\right|}.$

###  2.2 Quadtree data structure

For adaptive refinement, we define a root cell as the computational domain
$\mathbf{D}$ and recursively split each cell into quadrants until the desired
level of refinement is achieved; that is, a parent cell is divided into four
children cells of equal size. Here, the level of the cell is defined as the
number of refinements used for its construction. The notation $n/m$ describes
the level a quadtree grid, where the minimum level is $n$, and the maximum
level is $m$. Note that a quadtree grid of level $n/n$ is equivalent to a
$2^{n}\times 2^{n}$ uniform Cartesian grid. In this study, we consider non-
graded quadtree grids, and thus the size ratio of adjacent cells is not
restricted. When a partial differential equation on an irregular domain is
solved on the Cartesian grid, the order of the local truncation error near the
boundary is usually lower than that of interior grid points. Therefore, as
suggested in many works [30, 41, 19], we split a cell when its size exceeds
its distance to the boundary. More precisely, a cell $C$ is refined if the
following condition is satisfied:

$\underset{v\in\text{vertices}(C)}{\text{min}}\mid\phi(v)\mid\leq
L\sqrt{2}\Delta x,$

where $L$ refers to the Lipschitz constant of the level set function $\phi$,
and $\Delta x$ denotes the size of cell $C$ so that $\sqrt{2}\Delta x$
indicates the size of the cell diagonal. As we reinitialize the level set
function that allows $\phi$ to become a signed distance function, we take the
Lipschitz constant to be $L=1.2$. Figure 1 shows an example of a quadtree grid
generated in a circular domain. This locally refined grid structure enables a
storage reduction compared to the uniform Cartesian grid, whose size is the
same as that of the finest cell of the quadtree.

Figure 1:  An example of the quadtree mesh over a circular domain. The grid
has a maximum resolution of min level$=4$ and max level$=7$.

###  2.3 Grid layout

A staggered MAC grid has typically been used [19, 33] to solve the
incompressible Navier–Stokes equations (1) on quad/octree data structures. In
a two-dimensional domain, the horizontal and vertical velocities are located
on the horizontal and vertical cell faces, respectively, and the pressure
variables are located at cell centers. One major difficulty of using a MAC
grid on a quad/octree data structure is the discretization of the advection
and viscous terms. For example, consider the discretization at the cell face
marked with a circle in Figure 2(a). The cell face marked with a triangle
represents the $x$ component velocity $u$, which intersects the target cell
face or is a face of the same cell. Within this compact stencil, $u_{x}$ and
$u_{xx}$ cannot be discretized with errors of $o(\Delta x^{2})$ and $o(\Delta
x)$, respectively. To address this problem, Olshanskii et al. [33] used wider
stencils to discretize the advection and diffusion terms. By contrast, Guittet
et al. [19] addressed these difficulties by adopting a semi-Lagrangian scheme
for the advection term and a finite volume approach to treat the viscous term.
However, the quantities defined on the faces were interpolated to nodal values
to employ the semi-Lagrangian scheme, and a Voronoi mesh with respect to the
cell faces was constructed for the diffusion term. Furthermore, to ensure the
stability, Olshanskii et al. developed an explicit filter that modifies the
advection term near the coarse-to-fine grid interface, and Guittet et al.
derived a discrete divergence operator that uses faces from different cells.

Recognizing these difficulties, we use an Arakawa B grid [1] on the
quad/octree data structure for spatial discretization, as shown in Figure
2(b). That is, the velocity components are located at cell vertices, and the
pressure variables are located at cell centers. One advantage of this grid
layout is that the advection and diffusion terms can be discretized using
compact stencils with a consistent truncation error. The Arakawa B grid layout
is equivalent to the $Q1/P0$ finite element method and is known to exhibit
checkerboard instability. By adopting a stabilizing technique used in the
finite element community, we develop a second-order monolithic method with
reduced checkerboard instability for solving (1) in non-graded quad/octree
grids in the following section.

Figure 2: (a) Adjacent compact stencil in MAC grid. (b) General layout of the
Arakawa B-grid on the quadtree mesh. Velocity variables are defined at the
grid corners(red circles) whereas the pressure variables are stored at the
center of cells(blue squares).

##  3 Numerical method

In this section, we present a finite-difference-based second-order-accurate
monolithic scheme for solving the incompressible Navier–Stokes equations (1)
on the non-graded quadtree structure described in Section 2. Time is
discretized using the semi-Lagrangian method combined with the second-order
backward difference formula introduced by Xiu et al. [42]. By implicitly
discretizing the incompressibility condition, we obtain the following temporal
discretization for (1):

$\begin{cases}\frac{3\mathbf{u}^{n+1}-4\mathbf{u}_{d}^{n}+\mathbf{u}_{d}^{n-1}}{2\triangle
t}&=-\nabla p^{n+1}+\frac{1}{Re}\triangle\mathbf{u}^{n+1}+f^{n+1},\\\
\nabla\cdot\mathbf{u}^{n+1}&=S(p^{n+1}).\end{cases}$ (2)

Here, the subscript $d$ indicates departure points in the semi-Lagrangian
scheme, and $S$ is a stabilizing operator used to ease the checkerboard
instability of the pressure. In the next section, we describe the spatial
discretization and then discuss the induced linear system briefly. We point
out that the suggested discretization uses only the nodal values of adjacent
cells.

###  3.1 Advection term

The advection term in (1) is treated by a semi-Lagrangian method, which
discretizes the Lagrangian derivative of the velocity field in time instead of
the Eulerian derivative. The semi-Lagrangian scheme consists of tracking a
fluid particle backward by time integration of a characteristic equation and
recovering the values at the departure point using suitable high-order
interpolation schemes. The departure points are traced by a second-order
Runge–Kutta method as

$\displaystyle\mathbf{x}_{d}^{n}$
$\displaystyle=\mathbf{x}_{ij}^{n+1}-\triangle
t\mathbf{u}^{n+\frac{1}{2}}\left(\mathbf{x}_{ij}^{n+1}-\frac{\triangle
t}{2}\mathbf{u}^{n}\left(\mathbf{x}_{ij}^{n+1}\right)\right),$

and

$\displaystyle\mathbf{x}_{d}^{n-1}$
$\displaystyle=\mathbf{x}_{ij}^{n+1}-2\triangle
t\mathbf{u}^{n}\left(\mathbf{x}_{ij}^{n+1}-\triangle
t\mathbf{u}^{n}\left(\mathbf{x}_{ij}^{n+1}\right)\right).$

The velocity $\mathbf{u}^{n+\frac{1}{2}}$ at the intermediate time step is
approximated by the second-order extrapolation

$\mathbf{u}^{n+\frac{1}{2}}\left(\mathbf{x}\right)=\frac{3}{2}\mathbf{u}^{n}\left(\mathbf{x}\right)-\frac{1}{2}\mathbf{u}^{n-1}\left(\mathbf{x}\right).$

Because the points $\mathbf{x}_{ij}^{n+1}$ and $\mathbf{x}_{d}^{n}$ usually do
not coincide with grid nodes, appropriate spatial interpolations with nodal
values are required to obtain the velocity values at these points. As the
domain evolves with time, these points may lie outside the domain at preceding
time steps, and neighboring nodes may be external points, making interpolation
impossible. Here we design an interpolation scheme to compensate for these
problems using boundary ghost points. To describe the interpolation procedure
formally, we use the following notation. Suppose we would like to approximate
$\mathbf{u}^{n}$ at a point $\mathbf{x}$ and denote by $\mathbf{C}_{0}$ the
leaf cell in $\Omega^{n}$ containing $\mathbf{x}$.

Figure 3: Grid illustrations for discretizing advection term: (a) Iterative
procedure of finding cell for interpolation, (b) neighboring cells of a node
near the boundary.

When all the vertices of the cell $\mathbf{C}_{0}$ belong to $\Omega$,
$\mathbf{u}^{n}\left(\mathbf{x}\right)$ is approximated by the quadratic
essentially non-oscillatory (ENO) interpolation. Let
$\mathbf{x}_{00},\mathbf{x}_{10},\mathbf{x}_{01},\mathbf{x}_{11}$ denote the
lower left, lower right, upper left, and upper right vertices at the four
corners of the cell, and $u_{00},u_{10},u_{01},u_{11}$ represent the value of
$u$ at these four vertices, respectively. In addition, let
$\left(\theta_{x},\theta_{y}\right)=\frac{\mathbf{x}-\mathbf{x}_{00}}{\Delta
x}$ for the cell size $\Delta x$. Then the quadratic ENO interpolation is
defined as

$\displaystyle\mathbf{u}(x,y)$
$\displaystyle=(1-\theta_{x})(1-\theta_{y})\mathbf{u}_{00}+(1-\theta_{x})\theta_{y}\mathbf{u}_{01}+\theta_{x}(1-\theta_{y})\mathbf{u}_{10}+\theta_{x}\theta_{y}\mathbf{u}_{11}$
$\displaystyle-\frac{\theta_{x}(1-\theta_{x})\Delta
x^{2}}{2}D^{min}_{xx}\mathbf{u}-\frac{\theta_{y}(1-\theta_{y})\Delta
x^{2}}{2}D^{min}_{yy}\mathbf{u},$

where

$D^{min}_{xx}\mathbf{u}=\text{minmod}\left(D_{xx}\mathbf{u}_{ij}\right)$

for $i,j=0,1$. Here, $minmod$ denotes a function that returns the value that
has the minimum absolute value, and $D_{xx}u_{ij}$ is a discretization of
$u_{xx}$ at $\mathbf{x}_{ij}$ according to the second-order finite differences
in [29, 30]. For the detailed formula, see Section 3.2.

When one of the vertices does not belong to $\Omega$, interpolation is
performed as follows.

  1. 1.

Replace $\mathbf{C}_{0}$ with another cell that has at least one vertex
belonging to $\Omega$.

  2. 2.

Approximate the velocity using the moving least-squares (MLS) method.

When $\mathbf{C}_{0}$ has a vertex that belongs to $\Omega$, we skip step 1
above. By contrast, when there is no vertex that belongs to $\Omega$, we adopt
the bisection algorithm introduced in [12]. First, we set
$\mathbf{x}^{0}=\mathbf{x}$. Then we repeat the sequence

$\mathbf{x}^{k+1}=\mathbf{x}^{k}-\frac{\triangle
x}{2}\mathbf{n}(\mathbf{x}^{k})$

until the cell containing $\mathbf{x}^{k}$ has a vertex that belongs to
$\Omega$. After the iteration terminates, we set $\mathbf{C}_{0}$ to be the
cell containing $\mathbf{x}^{k}$. This process is schematically illustrated in
Figure 3(a).

After the cell $\mathbf{C}_{0}$ is selected, we apply the MLS method. A linear
polynomial $P_{u}$ with basis $\\{1,x,y\\}$ is constructed so as to minimize

$\sum_{i}w(\mathbf{x},\mathbf{x}_{i})\left(P_{u}(\mathbf{x_{i}})-u_{i}\right)^{2},$

where $w(\mathbf{x},\mathbf{y})=e^{-|\mathbf{x}-\mathbf{y}|^{2}}$ is a
Gaussian weight function. To construct a sample point
$\\{\left(\mathbf{x}_{i},u_{i}\right)\\}$, we first set
$\mathbf{x}_{0}=\left(x_{0},y_{0}\right)=argmax_{\mathbf{x}\in
vertices(\mathbf{C}_{0})}-\phi(\mathbf{x})$. Let $\mathbf{C}_{j}$ denote a
cell such that $\mathbf{x}_{0}\in\partial\mathbf{C}_{j}$. Then we define
$\\{\left(\mathbf{x}_{i}\right)\\}=\mathbf{X}_{in}\cup\mathbf{X}_{bd}$, where

$\mathbf{X}_{in}=\\{\mathbf{y}|\mathbf{y}\in
Vertices(\mathbf{C}_{j})\cap\Omega\\}$

and

$\mathbf{X}_{bd}=\\{\mathbf{y}|\mathbf{y}\in\mathbf{C}_{j}\cap\Gamma,\
\left(\mathbf{x_{0}}-\mathbf{y}\right)\mathbin{\\!/\mkern-5.0mu/\\!}\mathbf{e}_{i}\text{
for some }i=1,2\\}$

for the standard basis $e_{1}=(1,0),e_{2}=(0,1)$. Note that $\mathbf{X}_{in}$
consists of neighboring grid points that belong to $\Omega$, and
$\mathbf{X}_{bd}$ consists of the points that belong to $\Gamma$, where the
line segment with the point and $\mathbf{x}_{0}$ lies along a Cartesian axis.
In practice, we compute $\mathbf{X}_{bd}$ using the techniques in the GFM [27,
23]. For example, consider the situation illustrated in Figure 3(b). The grid
point $\mathbf{X}_{R}$ is not in the fluid domain $\Omega$. Let $\phi_{R}$ and
$\phi_{0}$ denote the values of the level set function $\phi$ at
$\mathbf{x}_{R}$ and $\mathbf{x}_{0}$, respectively. Then
$\hat{\mathbf{x}}_{R}=(x_{0}+\theta_{R}d_{R},y_{0})$ is in $\mathbf{X}_{bd}$,
where $\theta_{R}$ is computed as

$\theta_{R}=\frac{\mid\phi_{0}\mid}{\mid\phi_{0}\mid+\mid\phi_{R}\mid}.$

Sample points $\\{\left(\mathbf{x}_{i},u_{i}\right)\\}$ are then defined
naturally; $u_{i}$ denotes the value of $u$ on $\mathbf{x}_{i}$ when
$\mathbf{x}_{i}\in\Omega$, and $u_{i}$ is given by the boundary condition when
$\mathbf{x}_{i}\in\Gamma$. After finding the polynomial $P_{u}$, we simply
approximate

$\mathbf{u}(\mathbf{x})\approx P_{u}(\mathbf{x}).$

A complete description of the MLS method is given in [26, 24].

###  3.2 Diffusion term

The Laplacian velocity operator is discretized using the method introduced by
Min et al. [31]. We refer to a regular node as a node for which neighboring
nodes in all the Cartesian directions exist; see Figure 4(a). Let
$\mathbf{u}_{R},\ \mathbf{u}_{L},\ \mathbf{u}_{T},\mathbf{u}_{B}$ be the
values of $\mathbf{u}$ at the nodes immediately adjacent to $\mathbf{u}_{0}$,
and $d_{R},\ d_{L},\ d_{T},d_{B}$ denote the distances between
$\mathbf{x}_{0}$ and the neighboring nodes to the right, left, top, and
bottom, respectively. Then the Laplacian of $\mathbf{u}$ at the regular node
$\mathbf{x}_{0}$ is discretized as

$\triangle\mathbf{u}_{0}=\left[\frac{\mathbf{u}_{R}-\mathbf{u}_{0}}{d_{R}}-\frac{\mathbf{u}_{0}-\mathbf{u}_{L}}{d_{L}}\right]\frac{2}{d_{R}+d_{L}}+\left[\frac{\mathbf{u}_{T}-\mathbf{u}_{0}}{d_{T}}-\frac{\mathbf{u}_{0}-\mathbf{u}_{B}}{d_{B}}\right]\frac{2}{d_{T}+d_{B}}.$

Figure 4: Neighboring nodes of a regular node (a) and a T-junction node (b) in
quadtree. The dashed line is imaginary.

Nonuniform Cartesian grids produce hanging nodes, which appear only within the
region of finer resolution. Referring to figure 4(b), the Laplacian of
$\mathbf{u}$ at a T-junction node is discretized by

$\triangle\mathbf{u}_{0}=\left[\frac{\mathbf{u}_{R}-\mathbf{u}_{0}}{d_{R}}-\frac{\mathbf{u}_{0}-\mathbf{u}_{L}}{d_{L}}\right]\frac{2}{d_{R}+d_{L}}+\left(1-\frac{d_{LT}d_{LB}}{2d_{L}\left(d_{R}+d_{L}\right)}\right)\left[\frac{\mathbf{u}_{T}-\mathbf{u}_{0}}{d_{T}}-\frac{\mathbf{u}_{0}-\mathbf{u}_{B}}{d_{B}}\right]\frac{2}{d_{T}+d_{B}}$

for

$\mathbf{u}_{L}=\frac{d_{LB}}{d_{LB}+d_{LT}}\mathbf{u}_{LT}+\frac{d_{LT}}{d_{LB}+d_{LT}}\mathbf{u}_{LB}.$

When the node is close to the boundary $\Gamma$, some of the adjacent nodes
used in the discretization may lie outside the domain. For example, the right-
hand neighbor $\mathbf{x}_{R}$ may not belong to $\Omega$. In this case, as
explained in the previous section, we find the boundary point
$\hat{\mathbf{x}}_{R}$ by computing $\theta_{R}$. Then, by replacing $d_{R}$
with $\theta_{R}d_{R}$ and $u_{R}$ with the boundary condition, the Laplacian
can be discretized at the node close to the boundary.

###  3.3 Pressure gradient

We now consider the discretization of the pressure gradient operator. In the
Arakawa B grid structure, the pressure variables are located at cell centers.
Pressure nodes are usually not aligned horizontally or vertically with
velocity nodes. Therefore, when $p_{x}$ is discretized at the node
$\mathbf{x}_{0}$, the intermediate value $p_{L},p_{R}$ and corresponding
points $\mathbf{x}_{L},\mathbf{x}_{R}$, which have the same $y$ coordinate as
$\mathbf{x}_{0}$, are computed and used in the discretization.

First, consider the case when $\mathbf{x}_{0}$ is a regular node. Referring to
Figure 5(a), the pressure variables $p_{00},p_{01}$ on the left side of
$\mathbf{x}_{0}$ are linearly interpolated to define

$p_{L}\coloneqq\frac{s_{B}p_{01}+s_{T}p_{00}}{\theta_{T}+\theta_{B}}.$

The distance $d_{L}$ between $\mathbf{x}_{L}$ and $\mathbf{x_{0}}$ is then

$d_{L}\coloneqq 2\frac{s_{T}s_{B}}{s_{T}+s_{B}}.$

We can compute $p_{R}$ and $d_{R}$ in a similar manner. Then $p_{x}$ is
discretized as

$p_{x}=\frac{p_{R}-p_{L}}{d_{R}+d_{L}}.$

For the T-junction node $\mathbf{x}_{0}$ in the cell configuration shown in
Figure 5(b), the preceding discretization cannot be directly used to
approximate $p_{L}$ because there is only one adjacent cell to the left of
$\mathbf{x}_{0}$. In Figure 5(b), the $y$ coordinate of the node with $p_{00}$
is greater than that of $\mathbf{x_{0}}$. Then, the line segment connecting
the two pressure nodes with $p_{00}$ and $p_{10}$ intersects the line
$y=y_{0}$ on the left side of $\mathbf{x_{0}}$. Therefore, $p_{00},p_{10}$ and
its nodes are linearly interpolated to define $p_{L}$ and $d_{L}$ as

$p_{L}\coloneqq\frac{s_{B}p_{10}+s_{T}p_{00}}{\theta_{T}+\theta_{B}},$

and

$d_{L}\coloneqq\frac{\left(l-s_{T}\right)s_{B}}{s_{T}+s_{B}}.$

When the $y$ coordinate of the node with $p_{00}$ is less than that of
$\mathbf{x_{0}}$, $p_{11}$ and its node are used instead of $p_{10}$.

Figure 5: Neighboring nodes of a regular node (a) and a T-junction node (b)
for discretization of the pressure gradient. Dashed lines are imaginary.

The discretization of $\frac{\partial p}{\partial y}$ is defined similarly.
For nodes near the boundary, whenever one of the vertices of the cell belongs
to $\Omega$, ghost pressure nodes are defined at the cell. As a result, the
discretization described above can be directly applied to the nodes near the
boundary.

####  3.3.1 Alternative discretization in finite volume approach

There is a relatively simple approach to discretizing the pressure gradient
operator by the finite volume method. It does not require that an intermediate
value be computed, and the extension to octrees is natural and simple. As
shown in Figure 4, we define the control volume $V$ of the grid node
$\mathbf{x}_{0}=(x_{0},y_{0})$ as
$V=[x_{0}-\frac{d_{L}}{2},x_{0}+\frac{d_{R}}{2}]\times[y_{0}-\frac{d_{L}}{2},y_{0}+\frac{d_{T}}{2}]$.
Note that the control volume is consistent with the discretization of the
Laplacian discussed in Section 3.2. Then the divergence theorem yields

$\int_{V}p_{x}d\mathbf{x}=\int_{\partial
V}\left(p,0\right)\cdot\mathbf{n}_{V}ds,$

where $\mathbf{n}_{V}$ is the outward vector normal to $\partial V$. Let
$C_{i}$ denote the cell intersecting the control volume $V$ and $p_{i}$ denote
a pressure variable on the cell $C_{i}$. By interpreting the pressure variable
as a piecewise constant function on the cell, the integral can be approximated
in the finite volume sense as

$\int_{\partial
V}\left(p,0\right)\cdot\mathbf{n}_{V}ds=\sum_{i}\left(p_{i},0\right)\cdot\mathbf{n}_{V}|\partial
V\cap C_{i}|.$

Then, $p_{x}$ is discretized by dividing this approximation by the area of the
control volume $V$. For example, $\nabla p$ in Figure 4(a) is discretized by

$\displaystyle
p_{x}=\frac{4}{(d_{R}+d_{L})(d_{T}+d_{B})}\left((p_{11}-p_{01})\frac{d_{T}}{2}+(p_{10}-p_{00})\frac{d_{B}}{2}\right),$
$\displaystyle
p_{y}=\frac{4}{(d_{R}+d_{L})(d_{T}+d_{B})}\left((p_{11}-p_{10})\frac{d_{R}}{2}+(p_{01}-p_{00})\frac{d_{L}}{2}\right).$

On the hanging node shown in Figure 4(b), the discrete pressure gradient is
defined as

$\displaystyle
p_{x}=\frac{4}{(d_{R}+d_{L})(d_{T}+d_{B})}\left((p_{11}-p_{00})\frac{d_{T}}{2}+(p_{10}-p_{00})\frac{d_{B}}{2}\right),$
$\displaystyle
p_{y}=\frac{4}{(d_{R}+d_{L})(d_{T}+d_{B})}\left((p_{11}-p_{10})\frac{d_{R}}{2}+(p_{00}-p_{00})\frac{d_{L}}{2}\right).$

On the nodes near the boundary, the control volume is modified as described in
Section 3.2. If $\mathbf{x}_{R}\notin\Omega$, $\theta_{R}$ is computed such
that $\mathbf{x}_{0}+(\theta_{R}d_{R},0)\in\Gamma$, and $d_{R}$ is replaced
with $\theta_{R}d_{R}$. Note that the finite volume approach is motivated by
[37]. Although it is much simpler than the finite difference approach, it is
found numerically that the order of convergence is reduced to first order at
higher Reynolds number. The degradation of the order of accuracy is reported
in Table 1. For this reason, the finite difference approach is used in most of
the numerical experiments.

###  3.4 Divergence-free equation

The divergence of the velocities is discretized at each cell. Let
$\mathbb{\mathbf{C}}_{0}$ denote a cell, and define a set
$\mathbf{X}=\\{\mathbf{x}_{i}\\}_{i=1,\ldots,n}\subset\partial\mathbf{C}_{0}$
as

$\mathbf{x}\in\mathbf{X}\iff\mathbf{x}\in\partial\mathbf{C}_{0}\text{ is a
grid node }\quad\text{or}\quad\mathbf{x}\in\Gamma\cap\partial\mathbf{C}_{0},$

where $\mathbf{x}_{i}$ are ordered in the counterclockwise direction. Then
$\partial\left(\mathbf{C}_{0}\cap\Omega\right)$ is discretely approximated as
the union of line segments $\Gamma_{i}$:

$\partial\left(\mathbf{C}_{0}\cap\Omega\right)\approx\sum_{i=1}^{n}\Gamma_{i},$

where $\Gamma_{i}$ is a line segment connecting $\mathbf{x}_{i}$ and
$\mathbf{x}_{i+1}$. Here, $\mathbf{x}_{n+1}=\mathbf{x}_{0}$. The normal vector
$\mathbf{n}_{i}$ corresponding to $\Gamma_{i}$ is obtained by $90^{\circ}$
clockwise rotation of
$\left(\mathbf{x}_{i+1}-\mathbf{x}_{i}\right)/\|\mathbf{x}_{i+1}-\mathbf{x}_{i}\|$.
Then the divergence of the velocities is discretized in the finite volume
sense using the divergence theorem:

$\displaystyle\intop_{\mathbf{C}_{0}\cap\Omega}\nabla\cdot\mathbf{u}dxdy$
$\displaystyle=\intop_{\partial\left(\mathbf{C}_{0}\cap\Omega\right)}\mathbf{u}\cdot\mathbf{n}ds=\Sigma_{i=1}^{F}\intop_{\Gamma_{i}}\mathbf{u}\cdot\mathbf{n}_{i}ds=\Sigma_{i=1}^{n}\left(\frac{\mathbf{u}_{i}+\mathbf{u}_{i+1}}{2}\cdot\mathbf{n}_{i}\right)|\Gamma_{i}|.$
(3)

Precise discretization is demonstrated by considering the general cell
configuration shown in Figure 6. The set $\mathbf{X}$ consists of eight
points;
$\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3},\mathbf{x}_{6},\mathbf{x}_{7},\mathbf{x}_{8}$
are grid nodes that belongs to $\Omega$, and $\mathbf{x}_{4},\mathbf{x}_{5}$
are boundary points. As shown in the previous sections,
$\theta_{3},\theta_{5}$, which define the position of the boundary, are
computed by linear interpolation of $\phi$. Consequently, the divergence
operator is discretized as follows:

$\displaystyle\frac{1}{\left|\mathbf{C}_{0}\right|}\intop_{\Omega_{0}}\nabla\cdot\mathbf{u}dxdy$
$\displaystyle=\frac{1}{\left|\mathbf{C}_{0}\right|}\Sigma_{i=0}^{7}\left(\mathbf{u}_{i}\cdot\mathbf{n}_{i}\right)\mid\Gamma_{i}\mid$
$\displaystyle=\frac{1}{\left|\mathbf{C}_{0}\right|}\Sigma_{i=0}^{7}\left(\frac{\mathbf{u}_{i}+\mathbf{u}_{i+1}}{2}\cdot\mathbf{n}_{i}\right)d_{i}$
$\displaystyle=\left\\{\left(\frac{d_{2}}{2}u_{2}+\frac{d_{2}+\theta_{3}d_{3}}{d}u_{3}+\frac{\theta_{3}d_{3}}{2}\hat{u}_{4}\right)+d_{4}\frac{\hat{\mathbf{u}}_{4}+\hat{\mathbf{u}}_{5}}{2}\cdot\mathbf{n}_{4}\right.$
$\displaystyle+\left(\frac{\theta_{5}d_{5}}{2}\hat{v}_{5}+\frac{\theta_{5}d_{5}+d_{6}}{2}v_{6}+\frac{d_{6}+d_{7}}{2}v_{7}+\frac{d_{7}}{2}v_{8}\right)-\left(\frac{d_{8}}{2}u_{8}+\frac{d_{8}}{2}u_{1}\right)$
$\displaystyle\left.-\left(\frac{d_{1}}{2}v_{1}+\frac{d_{1}}{2}v_{2}\right)\right\\}\frac{1}{\left|\mathbf{C}_{0}\right|.}$

In the last expression, parentheses indicate boundary integrand values in the
counterclockwise direction from the right edge.

Figure 6: Fictitious cell configuration used in the discretization of
divergence operator.

Note that the discretization of the divergence-free equation (3) on a uniform
Cartesian grid is equivalent to the discretization of the bilinear-constant
velocity–pressure finite element method. It is also known as the $Q1-P0$
element and exhibits pressure checkerboard instability. To reduce the
instability, stabilizing techniques have been studied in the finite element
community. For example, the penalty method [22] reduces the instability by
decoupling the velocity and pressure variables. Other examples include local
and global stabilizers [5, 40], which add an artificial term to the
divergence-free equation, and the pressure gradient projection method [4, 15],
which was inspired by fractional step methods. Among many stabilizers, we
extend the global stabilizer in [40] to the quadtree grid as follows. Using
the same notation as in (3), let $\mathbf{C}_{i}$ be a neighboring cell to
$\mathbf{C}_{0}$ such that $\mathbf{C_{0}}\cap\mathbf{C_{i}}=\Gamma_{i}$, and
let $p_{i}$ denote the pressure variable in the cell $\mathbf{C}_{i}$. Then
the stabilizer is defined as

$S(p)=\epsilon\frac{1}{|\mathbf{C}_{0}|}\Sigma_{i=1}^{n}\left(\left(p_{i}-p_{0}\right)\cdot\mathbf{n}_{i}|\Gamma_{i}|\right)$

for $\epsilon=\Delta x_{s}^{2}$. In Section 4.1.2, the need for the stabilizer
is described.

##### Remarks

Note that the presented pressure gradient operator is rank-2 deficient on a
uniform Cartesian grid. For any real numbers $a,b\in\mathbb{R}$, if the
pressure variables are given as illustrated in Figure 7(a), $\nabla
p=\mathbf{0}$ at the interior vertices. However, for quadtree grids, this
problem is resolved. The discretization $\frac{\partial p}{\partial y}=0$ at
the vertex marked with a red circle in Figure 7(b) yields $p_{11}=p_{10}$;
consequently, the checkerboard pressure distribution shown in Figure 7(a) no
longer belongs to the kernel space unless $a=b$. Therefore, on nonuniform
Cartesian grids, the $p\to\nabla p$ operator is rank-1 deficient. To show that
the divergence operator $\nabla\cdot$ is rank-1 deficient, let us denote the
divergence operator by $\mathbf{D}$, and let $\mathbf{M}$ be a diagonal matrix
that maps a pressure variable to pressure variable and is defined as

$\mathbf{M}_{ii}=|\mathbf{C}_{i}|,$

where $|\mathbf{C}_{i}|$ is the area of the cell $\mathbf{C}_{i}$. It can be
shown that $\left(\mathbf{M}\circ\mathbf{D}\right)^{T}$, like the pressure
gradient operator, is rank-2 deficient on a uniform grid and is reduced to a
rank-1 deficient matrix on a quadtree grid. Note that the kernel of
$\left(\mathbf{M}\circ\mathbf{D}\right)^{T}$ is spanned by the vector of all
those at cell centers. Because the kernel space of the transpose matrix has
been determined, we project the right-hand side onto the range space before
solving the linear system. Let $r_{i}$ denote a variable located at the cell
center on the right-hand side. We define $\hat{r}_{i}$ as

$\hat{r}_{i}=r_{i}-\frac{1}{|\mathbf{C}_{i}|}\frac{\sum_{j}r_{j}|\mathbf{C}_{j}|}{\sum_{j}|\mathbf{C}_{j}|}.$

Then $r_{i}$ is replaced by $\hat{r}_{i}$.

Let $\mathbf{G}$ denote the linear operator
$p\to\left(\int_{V}p_{x},\int_{V}p_{y}\right)$ in the finite volume
discretization introduced in Section 3.3.1. Then it can be shown that the
discrete gradient operator is the negative adjoint of the discrete divergence,
that is, $-\mathbf{G}=\left(\mathbf{M}\circ\mathbf{D}\right)^{T}$. However,
this property does not hold for the proposed discretization of the pressure
gradient by the finite difference method. Despite the desirable matrix
property, we do not adopt the finite volume approach because of the lower
order of convergence at high Reynolds number.

(a) Checkerboard type kernel with rank 2 (b) Constant pressure kernel induced
by discretization at the hanging node

Figure 7: Illustrations of pressure variable that belongs to kernel space of
gradient operator.

###  3.5 Extension to three spatial dimensions

The proposed method can be extended to three dimensions on cubical domains
with octree data structure. The discretization of the advection and diffusion
terms can be naturally extended to three spatial dimensions as in [29, 31].
However, the discretization of the pressure gradient and incompressibility
condition require modifications.

We first describe the discretization of the pressure gradient. Consider the
discretization of the derivative of the pressure in the $x$ direction at a
node $\mathbf{x}_{0}$. As shown in Section 3.3, a standard approach is to
compute the intermediate values $p_{R}$ and $p_{L}$ with their locations
$\mathbf{x}_{R}$ and $\mathbf{x}_{L}$, respectively, which have the same $y$
and $z$ coordinates as $\mathbf{x}_{0}$. In contrast to the two-dimensional
case, three or four cells are used to approximate the intermediate values. If
four cells are chosen, bilinear interpolation on the $yz$ plane is applied to
approximate the intermediate values. If three cells are used, linear
interpolation is applied instead of bilinear interpolation. For example,
suppose we have chosen three cells whose center is
$\left(x_{i},y_{i},z_{i}\right)$ for $i=1,2,3$, and let us denote by $p_{i}$
the pressure at the corresponding cell. Linear polynomials $P_{p}$ and $P_{x}$
are constructed on the $yz$ plane so that

$\displaystyle P_{p}(y_{i},z_{i})=p_{i},\quad P_{x}(y_{i},z_{i})=x_{i}.$

Then $p_{R}$ and $x_{R}$ are approximated as $p_{R}=P_{p}(y_{0},z_{0})$ and
$x_{R}=P_{x}(y_{0},z_{0})$.

The remaining difficulty is the suitable choice of cells, which depends on the
number of adjacent cells in the positive $x$ direction with respect to
$\mathbf{x}_{0}$, that is, cells for which the maximum $x$ coordinate is
larger than $x_{0}$.

When there are three or four cells in the positive $x$ direction, the cells
are used to construct linear or bilinear polynomials, respectively. However,
when there are only two cells or one cell in the positive $x$ direction, we
use the cell information from the negative $x$ direction in the linear
interpolation. Specifically, we select cells so that the point
$\left(y_{0},z_{0}\right)$ lies within a triangle with vertices
$(y_{1},z_{1}),(y_{2},z_{2})$, and $(y_{3},z_{3})$. Figure 8 shows an example.

(a) Two cells are located in positive $x-$direction (b) One cell is located in
positive $x-$direction

Figure 8:  General cell configurations for three cases considered for the
discretization of the pressure gradient. These represent projected grids onto
$yz$-plane. Cells in negative $x-$directions are illustrated by shaded regions
and the red dot is the evaluation node $\mathbf{x}_{0}$.

Now, consider the discretization of the divergence operator. As in two spatial
dimensions, we use the divergence theorem to approximate the divergence
operator in a cell $\mathbf{C}$. In three dimensions, the line fraction is
replaced with a face fraction, and thus the line integration becomes an area
integration. The discrete form of the divergence operator is given by

$\nabla\cdot\mathbf{u}=\frac{1}{\mid\mathbf{C}\mid}\Sigma_{F\in\mathcal{F}\mathbf{(}C)}\mid
F\mid\left(\mathbf{u\cdot n}\right)_{F},$

where $\mathcal{F}\left(\mathbf{C}\right)$ denotes the set of all faces $F$ of
the cell $\mathbf{C}$, and $\mathbf{n}$ is the outward unit normal to each
face. $\mid\mathbf{C}\mid$ represents the volume of the cell $\mathbf{C}$,
whereas $\mid\mathbf{F}\mid$ is the area of the cell face. Moreover,
$\left(\mathbf{u\cdot n}\right)_{F}$ is discretized by the average of
$\mathbf{u}$ at the four corner vertices of the face. As shown in Figure 9,
the discretization of the area integration with the normal vector parallel to
the $x$ axis is calculated as

$\displaystyle\mid\mathbf{C}\mid\frac{\partial u}{\partial x}$
$\displaystyle\approx\left(\frac{\triangle
x}{2}\right)^{2}\frac{u_{5}+u_{6}+u_{8}+u_{9}}{4}+\left(\frac{\triangle
x}{2}\right)^{2}\frac{u_{6}+u_{7}+u_{9}+u_{10}}{4}$
$\displaystyle+\left(\frac{\triangle
x}{2}\right)^{2}\frac{u_{8}+u_{9}+u_{11}+u_{12}}{4}+\left(\frac{\triangle
x}{2}\right)^{2}\frac{u_{9}+u_{10}+u_{12}+u_{13}}{4}$ $\displaystyle-\triangle
x^{2}\frac{u_{1}+u_{2}+u_{3}+u_{4}}{4}.$

Figure 9: General configuration of an octree grid used for the discretization
of the divergence operator in the cell center of $\mathbf{C}$.

##  4 Numerical results

In this section, we present the results of several numerical experiments to
demonstrate the performance of the proposed method. We begin with numerical
tests that confirm the second-order accuracy of the proposed method. Next, we
provide numerical evidence that illustrates the performance of the method for
standard benchmark problems. We implemented all the numerical experiments in
C++ on a personal computer. To solve the saddle-point system (2) generated by
our discretization, we used the generalized minimal residual method with an
incomplete LU preconditioner.

###  4.1 Accuracy tests

To estimate the order of accuracy of the suggested method, we compare
approximated numerical solutions to an analytic solution. For a two-
dimensional example, consider a single-vortex problem whose exact solutions
are given by

$\displaystyle u(x,y,t)$ $\displaystyle=-\cos x\sin y\cos t,$ (4)
$\displaystyle v(x,y,t)$ $\displaystyle=\cos y\sin x\cos t,$ $\displaystyle
p(x,y,t)$ $\displaystyle=-\frac{1}{4}\cos^{2}t\left(\cos 2x+\cos 2y\right).$

The corresponding forcing term $f=\left(f_{1},f_{2}\right)$ is given by

$\displaystyle f_{1}(x,y,t)$ $\displaystyle=\cos x\sin y\left(\sin
t-2\frac{1}{Re}\cos t\right),$ $\displaystyle f_{2}(x,y,t)$
$\displaystyle=-\cos y\sin x\left(\sin t-2\frac{1}{Re}\cos t\right).$

The computational domain is taken to be
$\Omega=\left[-\frac{\pi}{2},\frac{\pi}{2}\right]^{2}$, and the terminal time
is set to $t=\pi$.

####  4.1.1 Random quadtree

To demonstrate the robustness of the proposed method, we first describe
numerical tests on randomly generated non-graded quadtrees of levels ranging
from $4/6$ to $7/9$ for a wide range of Reynolds numbers between $1$ and
$1000$. For each Reynolds number, a test was conducted in four randomly
generated grids with $\Delta t=\Delta x_{s}$. Figure 10 shows an example of
the randomly generated quadtree of level $4/6$ and a refined grid of level
$5/7$. The average numerical errors on four randomly generated grids and the
best-fit order of convergence are presented in Figure 11. The numerical
results indicate that the proposed method is second-order accurate for the
fluid velocity regardless of the grid formulation and Reynolds number. Second-
order convergence of the pressure of the velocity in $L^{\infty}$ norms was
not observed at relatively low Reynolds number; however, second-order
convergence in $L^{2}$ norms was observed.

Figure 10: A example of random quadtree of level (a) $4/6$ and refined (b)
$5/7$ used in example 4.1.1.

Figure 11: $L^{\infty}$ and $L^{2}$ errors of the velocity, pressure and
divergence on random qaudtree grid.

####  4.1.2 Irregular domains

To demonstrate the second-order accuracy in various irregular domains, we
conduct a numerical experiment with the same exact solutions as (4) with
$Re=1000$ on the four domains $\Omega=\\{(x,y)\in\mathbb{R}^{2}|\phi(x,y)\leq
0\\}$ in $\left[-\frac{\pi}{2},\frac{\pi}{2}\right]^{2}$:

Circle | : |  $\phi(x,y)=\sqrt{x^{2}+y^{2}}-1$ .
---|---|---
Ellipse | : |  $\phi(x,y)=\sqrt{\frac{\left(x\cos\frac{\pi}{6}-y\sin\frac{\pi}{6}\right)^{2}}{1.5^{2}}+\frac{\left(x\sin\frac{\pi}{6}+y\cos\frac{\pi}{6}\right)^{2}}{0.5^{2}}}-1$ .
Flower | : |  $\phi(r,\theta,t)=r-1-0.3\cos\left(5\theta\right)$ .
Moving Flower | : |  $\phi(r,\theta,t)=r-1-0.3\cos\left(5(\theta-t)\right)$ .

Note that the domains Circle, Ellipse, and Flower do not change with respect
to time, whereas the other domain, Moving Flower, rotates counterclockwise.
The domains and the generated quadtree grid of level $5/7$ at $t=\pi$ are
presented in Figure 12. Time step size is set to be $\Delta t=\Delta x_{s}$.
Second-order convergence of the velocity, pressure, and divergence in both the
$L^{\infty}$ and $L^{2}$ norms are observed in Figure 13. The results indicate
that the proposed method can handle irregular domains and moving domains.

In Table 1, we present the results obtained by discretizing $\nabla p$ in the
finite volume approach introduced in Section 3.3.1, which are compared with
those obtained by the finite difference approach. The approximation errors of
the velocity measured in the infinity norm and the order of accuracy at
Reynolds numbers of $Re=1,1000$ on the circular domain are reported. The
finite volume approach delivers second-order accuracy for $Re=1$; however, it
is reduced to first-order convergence at the higher Reynolds number,
$Re=1000$. Otherwise, the finite difference method shows satisfactory second-
order convergence independent of Reynolds number, and the magnitude of the
errors is lower at all resolutions.

We demonstrate the role of pressure stabilization by comparing the results
with those of an unstabilized system, that is, $\epsilon=0$. Table 2 reports
the order of convergence of the velocity and pressure on the Moving Flower
domain with Reynolds number $Re=1$. Although second-order accuracy of the
velocity was achieved in both cases, the convergence rate of the pressure was
degraded in the $L^{2}$ norm when the stabilizer was not applied. In addition,
Figure 15 shows the change in the $L^{\infty}$ errors of the pressure on a
quadtree grid of level $6/8$ with respect to time. Although both methods show
oscillations in the error, the magnitudes of the oscillations and errors are
smaller when the stabilizer is used.

(a) circle

(b) ellipse

(c) flower

(d) moving flower

Figure 12: Various test domains of example 4.1.2

Figure 13: $L^{\infty}$ and $L^{2}$ errors of the velocity, pressure and
divergence on various irregular domains.

Figure 14: $L^{\infty}$ and $L^{2}$ errors of the velocity, pressure and
divergence on random octree.

Figure 15: $L^{\infty}$ errors of pressure when (a) stabilizer is not applied and (b) stabilizer is applied for example 4.1.2. | $Re=1$ |  | $Re=1000$
---|---|---|---
| FV approach | order | FD approach | order |  | FV approach | order | FD approach | order
4/6 | 1.23.E-03 |  | 6.81.E-04 |  |  | 2.12.E-02 |  | 5.00.E-03 |
5/7 | 3.09.E-04 | 2.00 | 2.26.E-04 | 1.59 |  | 9.66.E-03 | 1.14 | 1.89.E-03 | 1.52
6/8 | 7.83.E-05 | 1.98 | 6.40.E-05 | 1.82 |  | 5.94.E-03 | 0.70 | 5.48.E-04 | 1.78
7/9 | 2.13.E-05 | 1.88 | 1.71.E-05 | 1.91 |  | 2.70.E-03 | 1.14 | 1.30.E-04 | 2.07
Table 1: Comparison of $L^{\infty}$ errors in velocities between finite volume and finite difference discretization for pressure gradient. | Not stabilized |  | Stabilized
---|---|---|---
|  $L^{\infty}$ errors of $\mathbf{u}$ | order |  $L^{2}$ errors of $p$ | order |  |  $L^{\infty}$ errors of $\mathbf{u}$ | order |  $L^{2}$ errors of $p$ | order
4/6 | 7.73.E-04 |  | 2.89.E-02 |  |  | 7.82.E-04 |  | 8.37.E-03 |
5/7 | 3.33.E-04 | 1.21 | 2.24.E-02 | 0.36 |  | 3.36.E-04 | 1.22 | 3.69.E-03 | 1.18
6/8 | 8.09.E-05 | 2.04 | 1.23.E-02 | 0.87 |  | 8.20.E-05 | 2.04 | 1.25.E-03 | 1.56
7/9 | 2.16.E-05 | 1.90 | 1.32.E-02 | -0.10 |  | 2.16.E-05 | 1.93 | 5.31.E-04 | 1.24
Table 2: Convergence comparison between results with/without the stabilizer
for example 4.1.2.

####  4.1.3 Random octree

As a three-dimensional example, we consider the following analytical solution
on the cubical domain $\Omega=\left[-\frac{\pi}{2},\frac{\pi}{2}\right]^{3}$:

$\displaystyle u(x,y,z,t)$ $\displaystyle=-2\cos x\sin y\sin z\cos t,$ (5)
$\displaystyle v(x,y,z,t)$ $\displaystyle=\sin x\cos y\sin z\cos t,$
$\displaystyle w(x,y,z,t)$ $\displaystyle=\sin x\sin y\cos z\cos t,$
$\displaystyle p(x,y,z,t)$ $\displaystyle=\frac{1}{4}\cos^{2}t\left(2\cos
2x+\cos 2y+\cos 2z\right).$

Simulations were conducted up to $t=\pi$ with $\Delta t=4\Delta x_{s}$ on a
random octree grid of levels ranging from $3/5$ to $6/8$. Two Reynolds
numbers, $Re=1,100$, are considered. Figure 14 presents the $L^{\infty}$ and
$L^{2}$ errors, which are compared to those of the analytical solutions.
Although the convergence rate of the pressure attained at $Re=1$ is lower than
1, we observe overall second-order convergence rates even for a random octree
grid.

###  4.2 Driven cavity flows

In this section, we consider the lid-driven cavity flow problem. Because there
is no exact solution of cavity flow, we compare the simulation results to
those reported in previous works. We first simulated cavity flow on a
rectangular domain and compared it with the reference results of Ghia et al.
[17]. The numerical simulation of driven cavity flow around a flower-shaped
object introduced by Coco [12] follows. In this example, we adaptively refine
the cell $\mathbf{C}$ when $\max_{\mathbf{x}\in
Vertices(\mathbf{C})}\frac{\|\nabla\mathbf{u}(\mathbf{X})\|_{2}}{\|\mathbf{u}\|_{\infty}}\Delta
x>\varepsilon$ for a tolerance $\varepsilon$ and cell size $\Delta x$. As a
two-dimensional example, we set the tolerance to $\varepsilon=0.05$.

####  4.2.1 Lid-driven cavity on a rectangular domain

We solve the lid-driven cavity flow problem on the computational domain
$\Omega=\left[0,1\right]^{2}$ with the boundary condition

$\mathbf{u}=\begin{cases}(1,0)\text{ if $y=1$, $0<x<1$}\\\ (0,0)\text{ if
$y=0$, $x=0,1$}\\\ \end{cases}.$

Numerical simulations were conducted for two Reynolds numbers, $Re=100$ and
$1000$, on quadtree grids of level $5/7$ and $6/8$, respectively.

In both cases, the simulations were run until $t=40$. Streamlines and the
quadtree grids are shown in Figure 16. In Figure 17, the results are
quantitatively compared with the reference solution of Ghia et al. [17]. The
results are in overall agreement with the reference solutions for all cases.

(a) $Re=100$

(b) $Re=1000$

Figure 16: Streamlines and quadtree grid at $t=40$ for lid driven cavity flow.

Figure 17: Comparison between velocity field in the driven cavity from example
4.2.1 and the reference results of Ghia et al. [17] of $(a),(b)$ $Re=100$, and
$(c),(d)$ $Re=1000$.

####  4.2.2 Lid-driven cavity flow around a flower-shaped object

Using the experimental setup designed by Coco [12], we tested the driven
cavity flow on the computational domain
$\Omega=\\{(x,y)\in\left[-1,1\right]^{2}\mid\phi(x,y)<0\\}$, where the level
set function is given in polar coordinates as

$\phi(r,\theta)=-r+0.5+0.15\sin(5\theta).$

As in 4.2.1, no-slip boundary conditions were used on all but the top wall,
where the wall velocity was specified by the unit horizontal velocity,
$\mathbf{u}=(1,0)$. A no-slip boundary condition was also imposed on the
boundary of the flower-shaped object. The cavity flow was studied for several
Reynolds numbers: at $Re=1,10,100$ with a quadtree grid of level $5/7$ and at
$Re=1000,5000$ with a quadtree grid of level $6/8$. Simulations were conducted
up to $t=10$ and $t=100$ for low ($Re=1,10,100$) and high ($Re=1000,5000$)
Reynolds numbers, respectively. Streamlines of the velocities at the terminal
time are plotted in Figure 18.

The results for the low Reynolds numbers show good agreement with those
reported in [12]. Two primary vortices are observed at the upper left and
upper right near the flower-shaped object for $Re=1,10,100$, and the vortex at
the upper left is not observed at $Re=1000$. According to [12], the flower-
shaped object prevents the formation of counterclockwise vortices in the lower
corners. However, we found that counterclockwise vortices appear at the lower
corners after a sufficiently long time at the high Reynolds numbers
($Re=1000,5000$). Furthermore, the primary vortex that appears in the standard
square cavity flow problem is split into five vortices around the flower-
shaped object.

(a) $Re=1,t=10$

(b) $Re=10,t=10$

(c) $Re=100,t=10$

(d) $Re=1000,t=10$

(e) $Re=1000,t=100$

(f) $Re=5000,t=100$

Figure 18: Streamlines of the results in example 4.2.2.

###  4.3 Rotating flower

As a third example, we consider the rotating flower-shaped object problem
proposed by Coco [12]. A flower-shaped object rotates counterclockwise about
the origin with angular velocity $\frac{2\pi}{5}$. Mathematically, the
computational domain at time $t$ can be expressed as
$\Omega=\\{(x,y)\in\left[-1,1\right]^{2}|\phi(x,y,t)<0\\},$ where the level
set function is given in polar coordinates by

$\phi(r,\theta,t)=-r+0.5+0.15\sin\left(5\theta-2\pi t\right).$

A no-slip boundary condition $\mathbf{u}=(0,0)$ is imposed on all four walls,
and the boundary condition $\mathbf{u}=\frac{2}{5}\pi\left(-y,x\right)$ is
imposed on the flower-shaped object, which agrees with the movement of the
interface.

Numerical experiments were conducted on a quadtree grid of level $6/8$ for two
Reynolds numbers, $Re=100$ and $10000$, whereas Coco performed experiments
using only $Re=100$. As in [12], we plotted streamlines of the resulting fluid
flow in the test with $Re=100$ every quarter rotation, that is, every 1.25 s,
as shown in Figure 19(f). The black dot is plotted on the same petal. Figure
19(f) shows that the positions of the vortices agree with those presented in
[12]. After a full rotation, we observe four clockwise vortices in the
corners. For the high Reynolds number, $Re=10000$, streamlines appear every 10
s in Figure 20(f). In contrast to the $Re=100$ case, vortices are formed in
each corner around $t=60$. The results from $t=60$ to $80$ shows that the
fluid converges to a steady-state solution.

(a) $t=0$

(b) $t=1.25$

(c) $t=2.5$

(d) $t=3.75$

(e) $t=5$

(f) $t=6.25$ Figure 19: Streamlines of the results in example 4.3 for
$Re=100$.

(g) $t=7.5$

(h) $t=8.75$

(i) $t=10$

(a) $t=0$

(b) $t=10$

(c) $t=20$

(d) $t=30$

(e) $t=40$

(f) $t=50$ Figure 20: Streamlines of the results of example 4.3 for
$Re=10000$.

(g) $t=60$

(h) $t=70$

(i) $t=80$

###  4.4 Application to two-dimensional rising bubble

Although the proposed method is applicable only to single-phase flow, it can
be easily extended to two-phase flow problems. The level set function $\phi$
divides $\Omega$ into $\Omega^{+}=\\{(x,y)\in\Omega\mid\phi>0\\}$ and
$\Omega^{-}=\\{(x,y)\in\Omega\mid\phi<0\\}$. The densities $\rho^{\pm}$ and
viscosities $\mu^{\pm}$ for $\Omega^{\pm}$ may differ. According to the delta
function equation, two-phase flow with surface tension can be described as

$\begin{gathered}\rho^{\pm}(\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u})=\mu^{\pm}\bigtriangleup\mathbf{u}-\nabla
p+\sigma\delta\kappa\mathbf{n}+\rho\mathbf{g},\\\
\nabla\cdot\mathbf{u}=0.\end{gathered}$ (6)

Here, $\mathbf{n}$ represents the outward normal pointing from $\Omega^{+}$ to
the interface, $\mathbf{g}=(0,-g)$ represents gravity, $\kappa$ denotes the
mean curvature of the interface, and $\sigma$ is the surface tension
coefficient. A time-step procedure similar to that for the single-phase flow
described in Section 3 is used to simulate two-phase flow, as follows:

  1. 1.

From time level $t^{n}$, compute the appropriate $\Delta t_{n}$, and define
$t^{n+1}=t^{n}+\Delta t_{n}$. Advect the level set to obtain $\phi^{n+1}$ at
time level $t^{n+1}$, and then reinitialize it.

  2. 2.

Generate a quadtree grid with respect to the level set $\phi^{n+1}$.

  3. 3.

Solve for $\mathbf{u}^{n+1}$ and $p^{n+1}$ at time level $t^{n+1}$.

For advection and the reinitialization of the level set function, we adopt the
second-order finite difference method on a quadtree grid [30]. By using the
numerical delta function, the Navier–Stokes equations can be discretized as

$\begin{gathered}\rho^{n+1}\left(\alpha_{0}\mathbf{u}^{n+1}+\alpha_{1}\mathbf{u}_{d}^{n}+\alpha_{2}\mathbf{u}_{d}^{n-1}\right)=\mu\bigtriangleup\mathbf{u}^{n+1}-\nabla
p^{n+1}+\sigma\delta(\phi^{n+1})\kappa+\rho^{n+1}\mathbf{g},\\\
\nabla\cdot\mathbf{u}^{n+1}=0,\end{gathered}$

where $\alpha_{0},\alpha_{1},\alpha_{2}$ are determined according to the
second-order backward differences as follows:

$\displaystyle\alpha_{0}$ $\displaystyle=\left(2\Delta t_{n}+\Delta
t_{n-1}\right)/(\Delta t_{n}^{2}+\Delta t_{n}\Delta t_{n-1}),\quad\alpha_{2}$
$\displaystyle=\Delta t_{n}/(\Delta t_{n}\Delta t_{n-1}+\Delta
t_{n-1}^{2}),\quad\alpha_{1}$ $\displaystyle=-\alpha_{0}-\alpha_{2}.$

The discretization of $\nabla p$ and $\nabla\cdot\mathbf{u}$ does not differ
from the discretization introduced in Section 3. Using notation similar to
that in Figure 4, the terms $\rho$, $\mu\bigtriangleup\mathbf{u}$, and
$\sigma\delta\kappa\mathbf{n}$ at $\mathbf{x}_{0}$ are discretized by

$\displaystyle\rho(\mathbf{x}_{0})$
$\displaystyle=\left(\rho^{+}-\rho^{-}\right)H(\phi_{0})+\rho^{-},$
$\displaystyle\mu\bigtriangleup\mathbf{u}(\mathbf{x}_{0})$
$\displaystyle=\left[\mu_{R}\frac{\mathbf{u}_{R}-\mathbf{u}_{0}}{d_{R}}-\mu_{L}\frac{\mathbf{u}_{0}-\mathbf{u}_{L}}{d_{L}}\right]\frac{2}{d_{R}+d_{L}}+\left[\mu_{T}\frac{\mathbf{u}_{T}-\mathbf{u}_{0}}{d_{T}}-\mu_{B}\frac{\mathbf{u}_{0}-\mathbf{u}_{B}}{d_{B}}\right]\frac{2}{d_{T}+d_{B}},$
$\displaystyle\sigma\delta\kappa\mathbf{n}(\mathbf{x}_{0})$
$\displaystyle=\sigma\delta(\phi_{0})\kappa\mathbf{n},$

where
$\mu_{R}=\left(\mu^{+}-\mu^{-}\right)H(\frac{\phi_{0}+\phi_{R}}{2})+\mu^{-}$
and $\mu_{L},\mu_{T},\mu_{B}$ are defined similarly. Here, $H$ and $\delta$
are numerical Heaviside and delta functions and are defined as

$H(\phi)=\left\\{\begin{array}[]{ll}0&\phi<-\epsilon\\\
\frac{1}{2}+\frac{\phi}{2\epsilon}+\frac{1}{2\pi}\sin\left(\frac{\pi\phi}{\epsilon}\right)&-\epsilon\leq\phi\leq\epsilon,\\\
1&\epsilon<\phi\end{array}\quad\delta(\phi)=\left\\{\begin{array}[]{ll}0&\phi<-\epsilon\\\
\frac{1}{2\epsilon}+\frac{1}{2\epsilon}\cos\left(\frac{\pi\phi}{\epsilon}\right)&-\epsilon\leq\phi\leq\epsilon\\\
0&\epsilon<\phi\end{array}\right.\right.$

with $\epsilon=1.5\Delta x$.

In this section, we consider the evolution of a bubble rising in a rectangular
domain $\left[-2,2\right]\times\left[-1,3\right]$. A bubble of radius $0.5$ is
initially located at the origin. We consider three-dimensionless quantities:
the Reynolds number $Re=7.77$, the Eötvös number $Eo=243$, and the Morton
number $Mo=266$. Then the corresponding density, viscosity, surface tension,
and gravity coefficient are determined:

$\rho^{+}=1,\quad\rho^{-}=0.001,\quad\mu^{+}=\frac{\rho^{+}}{Re},\quad\mu^{-}=0.01\mu^{+},\quad\beta=\frac{\mu^{+^{2}}}{\rho^{+}}\sqrt{\frac{Eo}{Mo}},\quad\mathrm{\leavevmode\nobreak\
g}=\left(0,-\frac{\rho^{+}\beta^{3}}{\mu^{+}}\right).$

Simulations were performed up to $t=2$, and $\Delta t$ was computed under the
following Courant–Friedrichs–Lewy condition:

$\Delta t=\min\left(\frac{\Delta
x}{\|\mathbf{U}\|_{\infty}},0.7\sqrt{\frac{\rho^{+}+\rho^{-}}{4\pi\sigma}\Delta
x^{3}}\right).$

In Figure 21, the terminal shapes of the bubble at time $t=2$ are shown on
uniform grids of $64\times 64,128\times 128$ and quadtree grids of level
$4/8,4/9$, and $5/9$. In addition, the relative area loss, rising velocity,
and center of mass are plotted in Figure 22. We observe that the results on
the quadtree grid are comparable to those on the uniform grid.

Figure 21: Visualization of the bubble at (a) uniform cartesian grid and (b)
adaptive quadtree grid for example 4.4.

(a) Relative area loss

(b) Rising velocity

(c) Center of mass

(d) Rising velocity (zoomed)

(e) Center of mass (Zoomed)

Figure 22: Relative area loss, rising velocity, and center of mass for rising
bubble in example 4.4.

##  5 Conclusion

We presented a second-order-accurate monolithic method for solving the two-
dimensional incompressible Navier–Stokes equations on irregular domains and
quadtree grids. In the proposed grid layout, velocities are located at cell
vertices, and pressure variables are located at cell centers. Simple and
effective second-order finite differences on the vertices of a quadtree grid
[29] were used to discretize the advection and diffusion terms of the
velocity. The pressure gradient was discretized using finite differences, and
ghost points were used for cells outside the domain. The divergence of the
velocity was discretized using the divergence theorem, whereas a cut-cell
method was used to treat the boundary conditions. Numerical experiments
included analytical and practical examples. For the analytical test, we
presented the order of convergence of the velocity and pressure in both the
$L^{\infty}$ and $L^{2}$ norms. The results indicate that the suggested method
is second-order accurate for velocity regardless of the Reynolds number and
domain shape. Second-order convergence of the pressure was observed only for
relatively high Reynolds numbers, and first-order convergence was obtained for
low Reynolds numbers. Classical lid-driven cavity flow and recent numerical
experiment of incompressible flow on irregular domains were conducted and
compared to the reference results. To demonstrate a more demanding
application, numerical simulation of a rising bubble was investigated.

Although an extension of the proposed method to an octree grid on a cubical
domain was presented, the full three-dimensional extension of the method to
irregular domains was not studied. The main difficulty in three dimensions is
the discretization of the divergence-free equation using the cut-cell
approach. As reported by Cheny et al. [9], there are many different cases, and
developing consistent discretization on irregular domains may require surface
integration on general polygons and curved elements. As a future work, we will
extend the proposed method on an octree grid to three-dimensional irregular
domains. Furthermore, the extension of the proposed method to a multiphase
flow problem and fluid–structure interaction will also be an interesting
future work.

## References

  * [1] Arakawa, A., Lamb, V.R.: Computational design of the basic dynamical processes of the ucla general circulation model.  General circulation models of the atmosphere 17(Supplement C), 173–265 (1977)
  * [2] Barrett, A., Fogelson, A.L., Griffith, B.E.: A hybrid semi-lagrangian cut cell method for advection-diffusion problems with robin boundary conditions in moving domains.  arXiv preprint arXiv:2102.07966 (2021)
  * [3] Bedrossian, J., Von Brecht, J.H., Zhu, S., Sifakis, E., Teran, J.M.: A second order virtual node method for elliptic problems with interfaces and irregular domains.  Journal of Computational Physics 229(18), 6405–6426 (2010)
  * [4] Blasco, J., Codina, R.: Space and time error estimates for a first order, pressure stabilized finite element method for the incompressible navier–stokes equations.  Applied Numerical Mathematics 38(4), 475–497 (2001)
  * [5] Bochev, P.B., Dohrmann, C.R., Gunzburger, M.D.: Stabilization of low-order mixed finite elements for the stokes equations.  SIAM Journal on Numerical Analysis 44(1), 82–101 (2006)
  * [6] Bochkov, D., Gibou, F.: Solving poisson-type equations with robin boundary conditions on piecewise smooth interfaces.  Journal of Computational Physics 376, 1156–1198 (2019)
  * [7] Brown, D.L., Cortez, R., Minion, M.L.: Accurate projection methods for the incompressible navier–stokes equations.  Journal of computational physics 168(2), 464–499 (2001)
  * [8] Chen, X., Li, Z., Álvarez, J.R.: A direct iim approach for two-phase stokes equations with discontinuous viscosity on staggered grids.  Computers & Fluids 172, 549–563 (2018)
  * [9] Cheny, Y., Botella, O.: The ls-stag method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties.  Journal of Computational Physics 229(4), 1043–1076 (2010)
  * [10] Cho, H., Kang, M.: Fully implicit and accurate treatment of jump conditions for two-phase incompressible navier-stokes equation.  arXiv preprint arXiv:2005.13724 (2020)
  * [11] Chorin, A.J.: A numerical method for solving incompressible viscous flow problems.  Journal of computational physics 135(2), 118–125 (1997)
  * [12] Coco, A.: A multigrid ghost-point level-set method for incompressible navier-stokes equations on moving domains with curved boundaries.  Journal of Computational Physics 418, 109623 (2020)
  * [13] Coco, A., Russo, G.: Finite-difference ghost-point multigrid methods on cartesian grids for elliptic problems in arbitrary domains.  Journal of Computational Physics 241, 464–501 (2013)
  * [14] Coco, A., Russo, G.: Second order finite-difference ghost-point multigrid methods for elliptic problems with discontinuous coefficients on an arbitrary interface.  Journal of Computational Physics 361, 299–330 (2018)
  * [15] Codina, R., Blasco, J.: Stabilized finite element method for the transient navier–stokes equations based on a pressure gradient projection.  Computer Methods in Applied Mechanics and Engineering 182(3-4), 277–300 (2000)
  * [16] Fedkiw, R.P., Aslam, T., Merriman, B., Osher, S.: A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method).  Journal of computational physics 152(2), 457–492 (1999)
  * [17] Ghia, U., Ghia, K.N., Shin, C.: High-re solutions for incompressible flow using the navier-stokes equations and a multigrid method.  Journal of computational physics 48(3), 387–411 (1982)
  * [18] Gibou, F., Fedkiw, R.P., Cheng, L.T., Kang, M.: A second-order-accurate symmetric discretization of the poisson equation on irregular domains.  Journal of Computational Physics 176(1), 205–227 (2002)
  * [19] Guittet, A., Theillard, M., Gibou, F.: A stable projection method for the incompressible navier–stokes equations on arbitrary geometries and adaptive quad/octrees.  Journal of Computational Physics 292, 215–238 (2015)
  * [20] Harlow, F.H., Welch, J.E.: Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface.  The physics of fluids 8(12), 2182–2189 (1965)
  * [21] Hellrung Jr, J.L., Wang, L., Sifakis, E., Teran, J.M.: A second order virtual node method for elliptic problems with interfaces and irregular domains in three dimensions.  Journal of Computational Physics 231(4), 2015–2048 (2012)
  * [22] Hughes, T.J., Liu, W.K., Brooks, A.: Finite element analysis of incompressible viscous flows by the penalty function formulation.  Journal of computational physics 30(1), 1–60 (1979)
  * [23] Kang, M., Fedkiw, R.P., Liu, X.D.: A boundary condition capturing method for multiphase incompressible flow.  Journal of Scientific Computing 15(3), 323–360 (2000)
  * [24] Lancaster, P., Salkauskas, K.: Surfaces generated by moving least squares methods.  Mathematics of computation 37(155), 141–158 (1981)
  * [25] LeVeque, R.J., Li, Z.: The immersed interface method for elliptic equations with discontinuous coefficients and singular sources.  SIAM Journal on Numerical Analysis 31(4), 1019–1044 (1994)
  * [26] Levin, D.: The approximation power of moving least-squares.  Mathematics of computation 67(224), 1517–1531 (1998)
  * [27] Liu, X.D., Fedkiw, R.P., Kang, M.: A boundary condition capturing method for poisson’s equation on irregular domains.  Journal of computational Physics 160(1), 151–178 (2000)
  * [28] Losasso, F., Gibou, F., Fedkiw, R.: Simulating water and smoke with an octree data structure.  In: ACM SIGGRAPH 2004 Papers, pp. 457–462 (2004)
  * [29] Min, C., Gibou, F.: A second order accurate projection method for the incompressible navier–stokes equations on non-graded adaptive grids.  Journal of computational physics 219(2), 912–929 (2006)
  * [30] Min, C., Gibou, F.: A second order accurate level set method on non-graded adaptive cartesian grids.  Journal of Computational Physics 225(1), 300–321 (2007)
  * [31] Min, C., Gibou, F., Ceniceros, H.D.: A supra-convergent finite difference scheme for the variable coefficient poisson equation on non-graded grids.  Journal of Computational Physics 218(1), 123–140 (2006)
  * [32] Ng, Y.T., Min, C., Gibou, F.: An efficient fluid–solid coupling algorithm for single-phase flows.  Journal of Computational Physics 228(23), 8807–8829 (2009)
  * [33] Olshanskii, M.A., Terekhov, K.M., Vassilevski, Y.V.: An octree-based solver for the incompressible navier–stokes equations with enhanced stability and low dissipation.  Computers & Fluids 84, 231–246 (2013)
  * [34] Osher, S., Sethian, J.A.: Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations.  Journal of computational physics 79(1), 12–49 (1988)
  * [35] Peskin, C.S.: Numerical analysis of blood flow in the heart.  Journal of computational physics 25(3), 220–252 (1977)
  * [36] Popinet, S.: Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries.  Journal of Computational Physics 190(2), 572–600 (2003)
  * [37] Purvis, J.W., Burkhalter, J.E.: Prediction of critical mach number for store configurations.  AIAA Journal 17(11), 1170–1177 (1979)
  * [38] Russo, G., Smereka, P.: A remark on computing distance functions.  Journal of computational physics 163(1), 51–67 (2000)
  * [39] Schroeder, C., Stomakhin, A., Howes, R., Teran, J.M.: A second order virtual node algorithm for navier–stokes flow problems with interfacial forces and discontinuous material properties.  Journal of Computational Physics 265, 221–245 (2014)
  * [40] Silvester, D.J., Kechkar, N.: Stabilised bilinear-constant velocity-pressure finite elements for the conjugate gradient solution of the stokes problem.  Computer Methods in Applied Mechanics and Engineering 79(1), 71–86 (1990)
  * [41] Theillard, M., Gibou, F., Saintillan, D.: Sharp numerical simulation of incompressible two-phase flows.  Journal of Computational Physics 391, 91–118 (2019)
  * [42] Xiu, D., Karniadakis, G.E.: A semi-lagrangian high-order method for navier–stokes equations.  Journal of computational physics 172(2), 658–684 (2001)
  * [43] Yoon, M., Yoon, G., Min, C.: On solving the singular system arisen from poisson equation with neumann boundary condition.  Journal of Scientific Computing 69(1), 391–405 (2016)