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

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|}.$$

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.

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.

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}$.

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

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}}.$$

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.

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}$.

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.

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.

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.

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.

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}.$

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

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

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.