# High order finite-difference ghost-point methods
for elliptic problems in domains with curved boundaries

Armando Coco, Giovanni Russo

###### Abstract

In this paper a fourth order finite difference ghost point method for the
Poisson equation on regular Cartesian mesh is presented. The method can be
considered the high order extension of the second ghost method introduced
earlier by the authors. Three different discretizations are considered, which
differ in the stencil that discretizes the Laplacian and the source term. It
is shown that only two of them provide a stable method. The accuracy of such
stable methods are numerically verified on several test problems.

##  1 Introduction

Elliptic equations in two space dimensions are often encountered in several
physical and mathematical contexts, ranging from electrostatic (the Poisson
equation describes the electric potential associated to a distribution of
charges) to electromagnetic (the Helmoltz equation models electromagnetic wave
propagation in the frequency domain), from incompressible gas dynamics (the
pressure has to satisfy an elliptic equation to impose that the flow is
incompressible) to linear elasticity. Because of its ubiquitous applications,
solution of elliptic problems has always been of great interest, from
modeling, analytical and numerical points of view.

The most classical approach to the solution of such equations is usually based
on finite difference or finite element methods. Typically, finite difference
schemes are extremely simple to implement producing efficient codes when
discretizing the equations on a regular Cartesian grid on a rectangular
domain. For non rectangular domains, the finite element approach has the
advantage of a greater flexibility, since it can more easily treat complex
geometries by adopting, for example, a triangulation which well approximates
the boundary of the domain.

More recently, advancements have been made in the development of finite
difference methods on regular Cartesian grids. These methods discretize
elliptic equations within arbitrary domains implicitly defined, for instance,
by a level set function [10, 11, 18, 19, 20, 27, 4].

Such schemes have the advantage offered by regular Cartesian mesh, and are
very efficient for the treatment of problems with time dependent boundaries,
since the grid does not need to be re-computed at every time step. Shortley-
Weller method [25] for the Poisson equation is based on the use of additional
grid points at the intersection of the boundary with grid axis, and on the use
of such additional points for the discretization of the Laplacian operator
near the boundary (see left panel in Fig. 2). It can be shown (see [25]) that
even if the discrete operator obtained in this way is only first order
accurate, the resulting scheme based on the five-point approximation of the
Laplacian applied to a Dirichlet problem guarantees second order accuracy for
the solution.

A different alternative is offered by ghost point methods, which are based on
the following idea. The computational domain $\Omega$ is immersed in a
rectangular region $R$, which is discretized by regular Cartesian mesh. The
grid points on the mesh are classified in three disjoint sets: internal points
(points inside $\Omega$), ghost points (points outside $\Omega$ with a
neighbor inside $\Omega$) and inactive points (see right panel of Figure 2).
The unknown function is defined only on internal and ghost points. The
discrete equations are defined on the internal points, while boundary
conditions are assigned by imposing suitable conditions on the ghost points.
In [10], for each ghost point $G$ the boundary condition is assigned on the
closest point $B_{G}$ on the boundary, while in [18] the boundary conditions
are assigned at the intersection between the boundary of $\Omega$ and the grid
axis, giving rise to two different values of the function on the ghost point.

More recently, a different version of ghost based method defined on a regular
Cartesian mesh for elliptic problems has been proposed [2], which is second
order accurate and guarantees better theoretical properties (symmetry,
coercivity) leading to a convergence proof under mild regularity assumptions
on the boundary of the domain $\Omega$.

In this paper we propose a fourth order extension of the second order method
introduced in [10] for the Poisson equation in two dimensions. Other high
order finite-difference methods for Elliptic problems are [15, 3, 21, 16, 14].

High order methods are more cost-effective when accurate solutions are needed.
The benefit appears when the solution itself is sufficiently regular, which
implies regularity of the data (i.e. regularity of the boundary
$\partial\Omega$, the data on the boundary and the right hand side of the
equation). For free boundary problems, in which the boundary of the domain
$\partial\Omega$ evolves according to an assigned level set equation, there
has been an increase interest in fourth order schemes for the evolution of the
level set function and for the computation of the corresponding signed
distance function [13], which allows a second order accurate computation of
the surface curvature.

###  1.1 Model problem

We consider the following mixed Dirichlet-Neumann problem for for the Poisson
equation

$\begin{split}-\Delta u&=f\text{ in }\Omega\\\ u&=g_{D}\text{ on
}\Gamma_{D}\\\ u&=g_{N}\text{ on }\Gamma_{N}\\\ \end{split}$ (1)

The compact domain $\Omega$ is a subset of a rectangle $R$ of
$\mathbb{R}^{2}$, i.e. $\Omega\subset R\subset\mathbb{R}^{2}$, and it is
implicitly defined by a level set function $\phi:R\to\mathbb{R}$,
$\Omega=\\{(x,y)\in R|\phi(x,y)\leq 0\\}$. We assume that the boundary of the
domain $\Gamma=\partial\Omega$, defined by $\Gamma=\\{(x,y)\in
R\,|\,\phi(x,y)=0\\}$, is sufficiently smooth, which is guaranteed if the
function $\phi$ is sufficiently regular in a neighborhood of $\Gamma$, and its
gradient is bounded away from zero. The boundary $\Gamma$ is partitioned in a
Dirichlet boundary $\Gamma_{D}$ and Neumann boundary $\Gamma_{N}$:
$\Gamma=\Gamma_{D}\cup\Gamma_{N}$, $\Gamma_{D}\cap\Gamma_{N}=\emptyset$. As
usual, we denote by $\Delta\equiv\partial^{2}/\partial
x^{2}+\partial^{2}/\partial y^{2}$ the Laplacian operator in $\mathbb{R}^{2}$.
The method proposed in this paper is high order accurate for sufficiently
regular source term and boundary data, therefore we assume that $f$, $g_{D}$
and $g_{N}$ are $\mathcal{C}^{1}$ functions in $\Omega$. The problem setup is
illustrated in Fig. 1.

Fig. 1: Setup of problem (1): the domain $\Omega$, defined by $\phi(x,y)\leq
0$, is embedded into a rectangular domain $R$. Mixed boundary conditions are
imposed on the frontier $\Gamma$ of the domain. In each point of $\Gamma$ the
unit normal pointing outside of the domain can be expressed as
${\mathbf{n}}=\nabla\phi/|\nabla\phi|$.

###  1.2 Review of the second order method

Here we give a brief review of the second order methods presented in [10]. The
rectangular region is meshed by a regular Cartesian square grid, and
throughout the paper we assume that the grid spacing $h=\Delta x=\Delta y$ is
considerably smaller than the smallest radius of curvature of the boundary
$\Gamma$, although the method can be extended to the case of Lipschitzian
boundaries by adopting the approach presented, for example, in [10, 27]. In
the internal points, the Laplacian is discretized by the classical second
order 5-point discrete Laplacian. For each ghost point, a condition is
assigned by imposing that a suitable reconstruction of the solution satisfies
the boundary condition on the projection of the ghost point on the boundary
(see right panel of Fig. 2. The simplest scheme is obtained by using bilinear
interpolation111Here by bilinear, biquadratic, bicubic, etc. interpolation we
mean interpolation obtained on a regular square grid, using multivariate
polynomials in two space variables, of degree at most $n$, with $n$
respectively, equal to 1, 2, 3, etc. in each variable. Such multivariate
polynomials will take the form
$p_{n}(x,y)=\sum_{i=0}^{n}\sum_{j=0}^{n}a_{ij}x^{i}y^{j}$. The name comes from
the consideration that each such interpolating polynomial can be obtained by
repeated interpolations in the coordinate directions by one dimensional
polynomials of degree at most $n$. For example, with reference to the right
panel of Fig. 2, assuming we have Dirichlet boundary conditions near ghost
point $G$, the condition is obtained by imposing that the bilinear
reconstructions from the values in the four points indexed by
$(i,j),(i,j+1),(i-1,j+1),(i-1,j)$, evaluated in $B_{G}$, is equal to
$g_{D}(B_{G})$. Notice that in the particular case in the figure, the
condition in point $G$ is independent on the conditions on the other ghost
points, so the value imposed in $G$ could be expressed in terms of the
internal points. However, in general this is not the case (as, for example,
for point $(i-1,j+2)$), so at the end one obtains a coupled system for the
internal and ghost values. It has been shown that this simple scheme provides
second order accuracy for the function $u$ in the case of Dirichlet boundary
conditions, and first order accuracy for Neumann conditions. Once the function
$u$ is known on internal and ghost points, the gradient can be easily computed
on the internal nodes by standard centered discretization. However the
technique based on bilinear interpolation only provides first order accuracy
for the evaluation of the gradient. A better accuracy is obtained by using
biquadratic interpolation, which is based on a 9-point stencil pointing in the
direction opposite to the normal, direction that we denote as upwind. In the
case of point $G$ in the figure, for example, the biquadratic interpolation
uses points $(i-\ell,j+1-k)$, $\ell=0,1,2$ and $k=0,1,2$. Such stencil
guarantees second order accuracy both for the function $u$ and its gradient.
In order to improve the efficiency in the solution of the resulting linear
system for the internal and ghost values, a geometric multigrid method has
been proposed, which is based on a suitable smoother constructed ad hoc for
the ghost point discretization.

The method has been generalized to the case in which the boundary $\Gamma$ is
only piecewise regular [10] and to the case of a more general elliptic
equation with piecewise continuous coefficients [11]. The method has been
successfully applied to problems in linear elasticity, with application to
volcanology [7]. The case of moving boundaries has been considered in [6]. An
improved multigrid applied to the ghost method has been proposed and analyzed
in [9], by optimizing the relaxation parameters associated with the boundary
condition to augment the method’s smoothing efficacy. Finally, such boundary
treatment has been adopted for hyperebolic problems, both in the case of fixed
and moving boundaries [5].

The convergence of the method has been studied in [8], in the particular case
of rectangular domains whose boundary does not fit with the Cartesian lines. A
general proof of convergence of the ghost method developed in [10] for a
general 2D domain has never been obtained.

A recent cut-FEM variant of the ghost method has been presented in [2]. The
method results in a symmetric, positive definite and coercive linear system,
for which it is possible to prove convergence for a general domain.

The rest of the paper is structured as follows. After this introduction, in
Section 2 we present the various space discretizations of the equation. The
third Section is devoted to the treatment of ghost points and boundary
conditions. Section 4 illustrates numerical results on two test problems. In
the final section we draw some conclusions.

Fig. 2: Left panel: Shortley-Weller stencil around grid point $(i,j)$. Right
panel: Ghost point discretization for the classic five-point discrete
Laplacian. Boundary condition on node $G$ are assigned by imposing that
polynomial that interpolated the solution on the (internal of ghost) nodes
sarisfied the boundary conditions on point $B_{G}$, as in [10], or on
intersection of the grid lines with the boundary (denoted by the asteriscs),
as in [17].

##  2 Numerical Method

The computational domain $\mathcal{D}=[-1,1]^{2}$ is discretized by a uniform
Cartesian grid of $(n_{x}+1)(n_{y}+1)$ equally spaced grid points
$(x_{i},y_{j})$, with $x_{i}=-1+ih_{x}$, $y_{j}=-1+jh_{y}$ for
$i=0,\ldots,n_{x}$, $j=0,\ldots,n_{y}$, where $h_{x}=2/n_{x}$ and
$h_{y}=2/n_{y}$ are the spatial steps in the $x$ and $y$ directions,
respectively. For simplicity, we assume that $n_{x}=n_{y}=N$, then
$h_{x}=h_{y}=h$, but the entire methodology can be easiliy extended to the
general case $h_{x}\neq h_{y}$. The set of grid points is denoted by
$\mathcal{D}_{h}$. Grid points that are inside $\Omega$ are called internal
grid points, forming the discrete version of $\Omega$, denoted by
$\Omega_{h}=\mathcal{D}_{h}\cap\Omega$. The remaining points are called
external grid points. We use the notation $\omega_{i,j}$ to denote the
approximation of $\omega(x_{i},y_{j})$ for any function
$\omega\colon\Omega\rightarrow\mathbb{R}$.

###  2.1 Discretization of the elliptic equation

Several fourth order finite difference discretization of (1) are possible [15,
3, 21, 16, 14]. In this paper, we examine three distinct discretizations based
on a fully coupled ghost-point approach, denoted as Method 1 , Method 2 and
Method 3 . We demonstrate that the first discretization is too ill-conditioned
to ensure convergence. In contrast, the other two exhibit the expected fourth-
order convergence rate.

The PDE (first equation of (1)) is approximated in each internal grid point
using a fourth-order discretization. We investigate two types of
discretizations: star-stencil and box-stencil discretizations. We use the
following matrix notation to denote the linear combination of grid values
$\omega_{i,j}$.

$S\omega_{i,j}=\sum_{p,q=-P}^{P}s_{p,q}\omega_{i+p,j+q}$ (2)

where $S\in{\mathbb{R}}^{(2P+1)\times(2P+1)}$, $P\in\\{0,1,2\\}$. With this
notation, the discretization of Eq. (1) takes the form

$S\,u_{i,j}=Mf_{i,j}+\mathcal{O}(h^{4})$

We consider the following two fourth order discretizations.

#### Star-stencil discretization

The star-stencil discretization consists of a natural higher order extension
of the classical second-order accurate central difference five-point star
stencil discretization. Second derivatives are discretized along the Cartesian
coordinates by using central difference formulas obtained by Taylor expansions
up to the desired order of accuracy. The elliptic equation
$-\nabla^{2}u(x_{i},y_{j})=f(x_{i},y_{j})$ is dscretized by:

$\frac{1}{12h^{2}}\begin{bmatrix}&&1&&\\\ &&-16&&\\\ 1&-16&60&-16&1\\\
&&-16&&\\\ &&1&&\\\ \end{bmatrix}u_{i,j}=f_{i,j}+\mathcal{O}(h^{4})$ (3)

for each internal grid point $(x_{i},y_{j})\in\Omega_{h}$, so $P=2$ for $S$
and $M=1$. This scheme is obtained by approximating the second derivative in
each coordinate direction by the corresponding fourth order five-point
discretization.

#### Box-stencil (Mehrstellen) discretization

The Mehrstellen discretization [12] uses a compact $3\times 3$ stencil around
the internal grid point $(x_{i},y_{j})$. The elliptic equation
$-\nabla^{2}u(x_{i},y_{j})=f(x_{i},y_{j})$ is discretized as:

$\frac{1}{6h^{2}}\begin{bmatrix}-1&-4&-1\\\ -4&20&-4\\\ -1&-4&-1\\\
\end{bmatrix}u_{i,j}=\frac{1}{12}\begin{bmatrix}&1&\\\ 1&8&1\\\ &1&\\\
\end{bmatrix}f_{i,j}+\mathcal{O}(h^{4})$ (4)

for each internal grid point $(x_{i},y_{j})\in\Omega_{h}$, so $P=1$ both for
$S$ and $M$. This scheme can be understood as follows. It can be shown that
matrix $S$ provides a fourth order approximation of the operator
$-\Delta-h^{2}\Delta^{2}$, i.e.

$Su(x_{i},y_{j})=-\Delta
u(x_{i},y_{j})-h^{2}\Delta^{2}u(x_{i},y_{j})+O(h^{4})$

while the matrix $M$ is a fourth order approximation of the operator
$I+h^{2}\Delta$, i.e.

$Mf(x_{i},y_{j})=f(x_{i},y_{j})+h^{2}Lf(x_{i},y_{j})=f(x_{i},y_{j})+h^{2}\Delta
f(x_{i},y_{j})+O(h^{4})$

where

$L=\frac{1}{h^{2}}\begin{bmatrix}&1&\\\ 1&-4&1\\\ &1&\\\ \end{bmatrix}$

denotes the matrix corresponding to the classic five point stencil Laplacian.
The approximation is then obtained by observing that

$-\Delta u-h^{2}\Delta^{2}u=f+h^{2}\Delta(\Delta u)=f+h^{2}\Delta f$

###  2.2 Domain representation by level-set

The physical domain $\Omega$ is implicitly described by a level-set function
$\phi:\mathcal{D}\rightarrow\mathbb{R}$ such that:

$\Omega=\left\\{(x,y)|\phi(x,y)<0\right\\}$
$\Gamma=\partial\Omega=\left\\{(x,y)|\phi(x,y)=0\right\\}.$

Level set methods are largely used to track interfaces and they have been
widely discussed in literature (see, for example, [20, 24]). Different level-
set functions can describe the same domain $\Omega$. For example, the disk
centered at $(x_{C},y_{C})$ with radius $r>0$ is described by both the level-
set functions $\phi_{1}(x,y)=(x-x_{C})^{2}+(y-y_{C})^{2}-r^{2}$ and
$\phi_{2}(x,y)=\sqrt{(x-x_{C})^{2}+(y-y_{C})^{2}}-r$. The latter is the so-
called signed-distance function:

$\phi(x,y)=\left\\{\begin{array}[]{rl}-d\left((x,y),\partial\Omega\right)&\mbox{
if }(x,y)\in\Omega,\\\ d\left((x,y),\partial\Omega\right)&\mbox{ if
}(x,y)\notin\Omega,\end{array}\right.$

where

$d\left((x,y),\partial\Omega\right)=\inf_{(\bar{x},\bar{y})\in\partial\Omega}d_{e}\left((x,y),(\bar{x},\bar{y})\right)$

is the distance between a point and the set $\partial\Omega$ ($d_{e}$ is the
Euclidean distance between points). For an arbitrary domain $\Omega$, the
signed-distance function can be obtained from a generic level-set function
$\phi_{0}(x,y)$, for example, by the reinitialization approach [26, 22, 13],
namely as the steady-state solution of the following evolutionary equation

$\frac{\partial\phi}{\partial
t}=\operatorname{sign}(\phi_{0})\left(1-\left|\nabla\phi\right|\right),\quad\phi=\phi_{0}\;\mbox{
at }\;t=0,$ (5)

where $t$ is a fictitious time that represents an iterative parameter. This
equation updates the signed distance function away from the zero level set
$\partial\Omega$, at a speed of one, conforming to the standard CFL-type
stability conditions for the time step. Consequently, in a single time step,
the signed distance function is updated by approximately one grid spacing. The
signed distance function is usually numerically more stable than a generic
level-set function, since the latter may develop steeper (or shallower)
gradients, while $|\nabla\phi|=1$ for the signed-distance function. We note
that it is sufficient to solve Eq. (5) only for a few time steps, since we
need to know the distance function only in the vicinity of the boundary.

The methods presented in this paper do not rely on the assumption that the
level-set is a signed distance function, therefore they can be applied with
generic level-set functions, provided the function is sufficiently smooth in a
neighborhood of $\partial\Omega$.

From the level-set function we can recover geometric properties at each
boundary point, such as the outward unit normal vector ${\mathbf{n}}$
(pointing from inside $\Omega$ towards outside) and the curvature $\kappa$:

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

##  3 Ghost-point method for boundary conditions

When the internal grid point $(x_{i},y_{j})$ is sufficiently close to the
boundary, the internal discretizations (3) and (4) require values of $u$ on
external grid points, since the stencil is intersected with $\partial\Omega$.
For example, in Fig. 3 (left panel), the star-stencil discretization of the
Laplace operator on $(x_{i},y_{j})$ requires the value of $u$ on purple grid
points, including the external grid points
$(x_{i+1},y_{j}),(x_{i+2},y_{j}),(x_{i},y_{j+1}),(x_{i},y_{j+2})$, while the
box-stencil discretization requires the value of $u$ on the points surrounded
by a boxed dashed line, including the external points
$(x_{i-1},y_{j+1}),(x_{i},y_{j+1}),(x_{i+1},y_{j+1}),(x_{i+1},y_{j})$. The aim
is to define a fictitious value of $u$ on external grid points (ghost value)
by enforcing the boundary conditions in such a way that the overall fourth-
order accuracy is not degraded.

The external grid points that are required by the discretization of the
Laplace operator are called ghost points. Depending on the type of the
internal discretization, namely star- or box- stencil, we have two different
sets of ghost points $\Gamma_{h}^{S}$ and $\Gamma_{h}^{B}$, respectively,
defined as follows. We first note that if a grid point $P$ is included in the
discretization stencil of the Laplacian operator centered at another grid
point $Q$, then $Q$ also appears in the stencil centered at $P$. Therefore, it
is possible to define the set of ghost points by considering the stencil
centered at the ghost point itself, rather than at an internal point. In
detail, for each external grid point $(x_{l},y_{m})$, we define the star-
shaped set of neighbor points with grid distance $k\in\\{1,2\\}$:

$\mathcal{N}_{l,m}^{S,(k)}=\left\\{(x_{l\pm k},y_{m}),(x_{l},y_{m\pm
k}),\right\\}\cap\mathcal{D}_{h}$ (6)

and the box-shaped set of neighbour points:

$\mathcal{N}_{l,m}^{B}=\left\\{(x_{l+k_{x}},y_{m+k_{y}}),\;k_{x},k_{y}=-1,0,1\right\\}\cap\mathcal{D}_{h}.$
(7)

The grid point $(x_{l},y_{m})$ is ghost if at least one of its neighbors is
internal. Formally, the sets of ghost points are defined as:

$(x_{l},y_{m})\in\Gamma_{h}^{S,(1)}\Longleftrightarrow\mathcal{N}_{l,m}^{S,(1)}\cap\Omega_{h}\neq\emptyset,$
$(x_{l},y_{m})\in\Gamma_{h}^{S,(2)}\Longleftrightarrow(x_{l},y_{m})\notin\Gamma_{h}^{S,(1)}\text{
and }\mathcal{N}_{l,m}^{S,(2)}\cap\Omega_{h}\neq\emptyset,$
$(x_{l},y_{m})\in\Gamma_{h}^{B}\Longleftrightarrow\mathcal{N}_{l,m}^{B}\cap\Omega_{h}\neq\emptyset.$

Sets $\Gamma_{h}^{S,(1)}$ and $\Gamma_{h}^{S,(2)}$ represent, respectively,
the first and second layer of ghost points. Finally, we define
$\Gamma_{h}^{S}=\Gamma_{h}^{S,(1)}\cup\Gamma_{h}^{S,(2)}$. See Fig. 3 (right
panel) for examples of neighbour grid points $\mathcal{N}_{l,m}^{B}$ and
$\mathcal{N}_{l,m}^{S,(k)}$ for $k=1,2$ for the ghost point $(x_{l},y_{m})$.

Fig. 3:  Left: Stencil adopted in the discretization of the Laplace operator
on a generic internal grid point $(x_{i},y_{j})$. For the star-stencil
discretization (3), the stencil is made by the purple grid points, while the
dashed line surrounds the stencil of the box-stencil discretization (4).
Right: Star-shaped set of neighbour points $\mathcal{N}_{l,m}^{S,(1)}$ (red
points) and $\mathcal{N}_{l,m}^{S,(2)}$ (blue points) for a generic ghost
point $(x_{l},y_{m})$, as defined in Eq. (6). The box-shaped set of neighbour
points $\mathcal{N}_{l,m}^{B}$ (defined in Eq. (7)) is surrounded by a dashed
line. The set of neighbors of external grid points is adopted to identify the
ghost points.

In the following, we refer to $\Gamma_{h}$ as either $\Gamma^{S}_{h}$ or
$\Gamma^{B}_{h}$, depending on the type of internal discretizations (3) or
(4), respectively. Grid points that are neither internal nor ghost points are
called inactive grid points. They form the set
$\mathcal{D}_{h}\backslash(\Omega_{h}\cup\Gamma_{h})$.

###  3.1 Boundary points

For each ghost point $(x_{l},y_{m})$ we determine the closest boundary point
$B_{l,m}$ by using the level-set function as follows 222The algorithm to
determine $B_{l,m}$ assumes that the level-set is known only on the grid
points $\phi_{i,j}$, namely the analytical expression is not available.
Firstly, we divide the space into four quadrants by the Cartesian lines
$x=x_{l}$ and $y=y_{m}$ and identify the quadrant $(s_{x},s_{y})=(\pm 1,\pm
1)$ in which we will search the boundary point, where $s_{x}=-1$ [$s_{x}=1$]
refers to the left [right] quadrant and $s_{y}=-1$ [$s_{y}=1$] refers to the
bottom [top] quadrant. For example, $(s_{x},s_{y})=(1,-1)$ means that the
boundary point must be found in the bottom-right quadrant. In Fig. 4 we have
$(s_{x},s_{y})=(-1,-1)$. The quadrant $(s_{x},s_{y})$ can be identified by
looking at the opposite direction of $\nabla\phi$, the latter being orthogonal
to the curve $\phi=\phi_{l,m}$. A second-order central difference
discretization of $\nabla\phi$ leads to:

$s_{x}=\operatorname{sign}(\phi_{l-1,m}-\phi_{l+1,m}),\quad
s_{y}=\operatorname{sign}(\phi_{l,m-1}-\phi_{l,m+1})$

where $\operatorname{sign}(x)$ is the sign function defined as
$\operatorname{sign}(x)=1$ if $x>0$, $\operatorname{sign}(x)=-1$ if $x<0$ and
$\operatorname{sign}(0)=0$.

Then, we define the $5\times 5-$point stencil with steps
$(r_{x},r_{y})\in\mathbb{N}^{2}$ in the quadrant $(s_{x},s_{y})$:

$\mathcal{S}_{l,m}^{(r_{x},r_{y})}=\left\\{(x_{l+s_{x}k_{x}r_{x}},y_{m+s_{y}k_{y}r_{y}}),\;k_{x},k_{y}=0,1,\ldots,4\right\\}.$
(8)

In this paper we will use $1\leq r_{x},r_{y}\leq 2$ and the specific choice of
the pair $(r_{x},r_{y})$ will be defined later for each test.

We say that the stencil is in the upwind direction, since the ghost point is
at one corner of the stencil rather than at the centre and the rest of the
stencil is taken in the opposite direction of the outer normal vector.

Fig. 4:  Examples of stencils $\mathcal{S}_{l,m}^{(r_{x},r_{y})}$ adopted in
the discretization of boundary conditions to determine the ghost point values.
Left: stencils adopted in Method 1 and Method 3 . For Method 1 , the green
point is a ghost point of $\Gamma^{S,(1)}_{h}$ and the purple point is a ghost
point of $\Gamma^{S,(2)}_{h}$. The stencils are surrounded by the dashed lines
of the same colour. Both stencils have $(r_{x},r_{y})=(1,1)$ in Eq. (8). The
respective boundary points are also represented. Right: stencils adopted in
Method 2 . The stencil for a ghost point of $\Gamma^{S,(2)}_{h}$ (purple
point) is enlarged, meaning that $(r_{x},r_{y})=(2,2)$ is adopted in Eq. (8).
The stencil is represented by the purple point (ghost point) and the purple
circles. The stencil for the green ghost point, that is a ghost point of
$\Gamma^{S,(2)}_{h}$, is unchanged with respect to the left plot.

####  3.1.1 Biquartic interpolation formulas

By using the $5\times 5-$point stencil $\mathcal{S}_{l,m}^{(r_{x},r_{y})}$, we
can interpolate a grid function $\omega_{ij}$ on an arbitrary point $(x,y)$ by
evaluating the interpolant
$\tilde{\omega}(x,y)=\sum_{p=0}^{4}\sum_{q=0}^{4}a_{pq}x^{p}y^{q}$ where the
coefficients $a_{pq}$ are determined by imposing
$\tilde{\omega}(x_{i},y_{j})=\omega_{ij}$ for each
$(x_{i},y_{j})\in\mathcal{S}_{l,m}^{(r_{x},r_{y})}$.

This bi-quartic interpolation can be easily computed by a tensor product of 1D
quartic interpolations. In detail, the 1D quartic interpolation
$\left.\tilde{\omega}(x)\right|_{x=x^{*}}$ of a grid function $\omega_{i}$ on
a 5-point stencil $x_{i}=x_{0}+ih$ for $i=0,\ldots,4$ is:

$\left.\tilde{\omega}(x)\right|_{x=x^{*}}=\sum_{p=0}^{4}c_{p}(\vartheta)\omega_{i+p}$

where $\vartheta=(x^{*}-x_{0})/h$ (see Fig. 5) and the interpolation
coefficients are:

$\displaystyle c_{0}(\vartheta)$
$\displaystyle=\frac{1}{24}(\vartheta-4)(\vartheta-3)(\vartheta-2)(\vartheta-1),$
$\displaystyle c_{1}(\vartheta)$
$\displaystyle=-\frac{1}{6}(\vartheta-4)(\vartheta-3)(\vartheta-2)\vartheta,$
$\displaystyle c_{2}(\vartheta)$
$\displaystyle=\frac{1}{4}(\vartheta-4)(\vartheta-3)(\vartheta-1)\vartheta,$
$\displaystyle c_{3}(\vartheta)$
$\displaystyle=-\frac{1}{6}(\vartheta-4)(\vartheta-2)(\vartheta-1)\vartheta,$
$\displaystyle c_{4}(\vartheta)$
$\displaystyle=\frac{1}{24}(\vartheta-3)(\vartheta-2)(\vartheta-1)\vartheta.$
(9)

Similarly, the first derivative is obtained by

$\left.\tilde{\omega}^{\prime}(x)\right|_{x=x^{*}}=\frac{1}{h}\sum_{p=0}^{4}c_{p}^{\prime}(\vartheta)\omega_{i+p}$

where $\vartheta=(x^{*}-x_{0})/h$ and the interpolation coefficients are:

$c_{0}^{\prime}(\vartheta)=\frac{1}{12}(2\vartheta-5)(5-5\vartheta+\vartheta^{2}),\quad
c_{1}^{\prime}(\vartheta)=\frac{1}{6}(24-52\vartheta+27\vartheta^{2}-4\vartheta^{3}),\quad$
$c_{2}^{\prime}(\vartheta)=\frac{1}{2}(\vartheta-2)(3-8\vartheta+2\vartheta^{2}),\quad
c_{3}^{\prime}(\vartheta)=\frac{1}{6}(8-28\vartheta+21\vartheta^{2}-4\vartheta^{3}),\quad$
$c_{4}^{\prime}(\vartheta)=\frac{1}{12}(2\vartheta-3)(1-3\vartheta+\vartheta^{2})$

$x_{0}$$x_{1}$$x_{2}$$x_{3}$$x_{4}$$x^{*}$$\vartheta\,h$$h$ Fig. 5: 1D grid
nodes for quartic interpolation at $x^{*}$.

In 2D, the bi-quartic interpolation formula is

$\tilde{\omega}(x,y)=\sum_{p,q=0}^{4}c_{p}(\vartheta_{x})c_{q}(\vartheta_{y})\omega_{l+s_{x}r_{x}p,m+s_{y}r_{y}q}$
(10)

and the first partial derivatives are interpolated as

$\displaystyle\frac{\partial\tilde{\omega}}{\partial x}(x,y)$
$\displaystyle=\frac{s_{x}}{h}\sum_{p,q=0}^{4}c^{\prime}_{p}(\vartheta_{x})c_{q}(\vartheta_{y})\omega_{l+s_{x}r_{x}p,m+s_{y}r_{y}q}$
(11a) $\displaystyle\frac{\partial\tilde{\omega}}{\partial y}(x,y)$
$\displaystyle=\frac{s_{y}}{h}\sum_{p,q=0}^{4}c_{p}(\vartheta_{x})c^{\prime}_{q}(\vartheta_{y})\omega_{l+s_{x}r_{x}p,m+s_{y}r_{y}q}$
(11b)

where

$\vartheta_{x}=s_{x}\frac{x-x_{l}}{r_{x}\,h},\quad\vartheta_{y}=s_{y}\frac{y-y_{m}}{r_{y}\,h}.$
(12)

####  3.1.2 Iteration method to determine the boundary points

The boundary point $B_{l,m}$ is determined by the following iterative method

${\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}+\alpha_{k}{\mathbf{p}}_{k},\quad{\mathbf{x}}^{(0)}=(x_{l},y_{m}).$

At each step we approximate $\phi({\mathbf{x}}^{(k)})$ and
$\nabla_{h}{\phi}({\mathbf{x}}^{(k)})$ by the bi-quartic interpolation
$\tilde{\phi}({\mathbf{x}}^{(k)})$ and
$\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(k)})$ as described in Sect. 3.1.1. The
search direction is an approximation of the unit normal vector

${\mathbf{p}}^{(k)}=\frac{\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(k)})}{|\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(k)})|}$
(13)

and the step length is proportional to the level-set function, normalized by
its gradient

$\alpha^{(k)}=\varepsilon\frac{\tilde{\phi}({\mathbf{x}}^{(k)})}{|\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(k)})|}$

Numerically, we directly compute the product:

$\alpha^{(k)}{\mathbf{p}}^{(k)}=\varepsilon\,\tilde{\phi}({\mathbf{x}}^{(k)})\frac{\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(k)})}{(\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(k)}))^{2}}.$

The constant of proportionality must be small respect to the spatial step, say
$\varepsilon=h/10$.

###### Remark 1

Selecting
${\mathbf{p}}^{(k)}=\frac{\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(0)})}{|\nabla_{h}\tilde{\phi}({\mathbf{x}}^{(0)})|}$
is sufficient for a signed distance function defined from the boundary to the
ghost point ${\mathbf{x}}^{(0)}$. However, in practice, complications may
arise. If the level-set function is not a signed distance function, if the
signed distance function has not been evolved up to ${\mathbf{x}}^{(0)}$, or
if the ghost point is adjacent to a concave boundary (where the gradient of
the signed distance function might be discontinuous), then the chosen search
direction might fail to intersect the boundary. Therefore, it is advisable to
update the search direction at each iteration point as in 13 to accommodate
these potential issues.

###  3.2 Enforcing boundary conditions

For each ghost point $(x_{l},y_{m})$, we write a linear equation involving
grid values $u_{i,j}$ on the stencil (8). The equation is obtained by
enforcing the boundary condition on $B_{l,m}$ up to the desired accuracy. Let
$\tilde{u}$ be the bi-quartic interpolant on the stencil (8) obtained as in
Sect. 3.1.1. If $B_{l,m}\in\Gamma_{D}$, then we use the Dirichlet boundary
condition and the equation is

$\tilde{u}(B_{l,m})=g_{D}(B_{l,m}),$ (14)

otherwise (namely $B_{l,m}\in\Gamma_{N}$) we enforce the Neumann boundary
condition:

$\nabla\tilde{u}(B_{l,m})\cdot\frac{\nabla\phi(B_{l,m})}{|\nabla\phi(B_{l,m})|}=g_{N}(B_{l,m}).$
(15)

###### Remark 2

The linear equation for a ghost point, either (14) or (15), may contain other
ghost values. For example, in Fig. 4 (left plot), the stencil for the purple
ghost point contains the green point, that is also a ghost point. Therefore,
the equations for ghost points are all coupled each other. Internal and ghost
equations are then solved simultaneously (see 1 to assemble the linear
system).

###### Remark 3

We observe that Eqs. (14) and (15) are linear equations of the grid values
$u_{i,j}$ in $\mathcal{S}_{l,m}^{(r_{x},r_{y})}\subseteq\mathcal{D}_{h}$. If
$\partial\Omega$ is sufficiently regular (say
$\partial\Omega\in\mathcal{C}^{1}$) and the spatial step is sufficiently small
with respect to the radius of curvature of $\partial\Omega$, then
$\mathcal{S}_{l,m}^{(r_{x},r_{y})}\subseteq\Omega_{h}\cup\Gamma_{h}$, where
$\Gamma_{h}$ is either $\Gamma^{S}_{h}$ or $\Gamma^{B}_{h}$ according to the
chosen internal discretization (3) or (4), respectively. We assume in this
paper that we are always in this case. However, if $\partial\Omega$ is only a
Lipschitzian boundary, in the vicinity of kink points it may happen that
$\mathcal{S}_{l,m}^{(r_{x},r_{y})}$ is not contained in
$\Omega_{h}\cup\Gamma_{h}$ regardless of the spatial step and one has to
modify the stencil accordingly, following a similar approach as in [10] for
the second order method.

We can label the inactive points of
$\mathcal{S}_{l,m}^{(r_{x},r_{y})}\backslash\left(\Omega_{h}\cup\Gamma_{h}\right)$
as ghost points and add them to $\Gamma_{h}$. For each additional ghost point
$(x_{l*},y_{m*})$, we compute again the stencil
$\mathcal{S}_{l*,m*}^{(r_{x},r_{y})}$ and check weather it is contained in
$\Omega_{h}\cup\Gamma_{h}$ or not. If not, we repeat the process again. After
a few iterations, we end up with the case where the stencil
$\mathcal{S}_{l,m}^{(r_{x},r_{y})}$ is contained in $\Omega_{h}\cup\Gamma_{h}$
for any ghost point $(x_{l},y_{m})$.

###### Remark 4

When the Mehrstellen discretization (4) is adopted, the grid values $f_{ij}$
are needed on a five-point stencil centred in $(x_{i},y_{j})$, not only on
$(x_{i},y_{j})$. This means that $f_{i,j}$ must be defined on some external
grid points. For example, in Fig. 3 the discretization of the Laplace equation
on the grid point $(x_{i},y_{j})$ requires the values of $f$ on the external
grid points $(x_{i+1},y_{j})$ and $(x_{i},y_{j+1})$. We do not have any
information on the boundary values of $f(x,y)$, meaning that we have to
extrapolate a grid function from internal to external grid points, without
using any boundary condition. The extrapolation must be done to the desired
accuracy, that is fourth-order accuracy in this paper. In the context of
arbitrary domains represented by level-set functions, the high-order
extrapolation can be obtained as in [1].

###  3.3 Linear Solver

The whole discretization process (at internal and ghost points) leads to a
linear system in the $u_{i,j}$ for $(x_{i},y_{j})\in\Omega_{h}\cup\Gamma_{h}$.
For simplicity, we assemble a linear system for the unknowns $u_{i,j}$ in all
grid points $(x_{i},y_{j})\in\mathcal{D}_{h}$ and keep in mind that the value
of $u_{i,j}$ on inactive points
$\mathcal{D}_{h}\backslash(\Omega_{h}\cup\Gamma_{h})$ is meaningless. We
investigate three methods in this paper:

  * •

Method 1 \- star-shaped internal discretization (3) and
$\Gamma_{h}=\Gamma^{S}_{h}$ with $(r_{x},r_{y})=(1,1)$ in (8) for all ghost
points

  * •

Method 2 \- star-shaped internal discretization (3) and
$\Gamma_{h}=\Gamma^{S}_{h}$ with $(r_{x},r_{y})=(k,k)$ in (8) for ghost points
in $\Gamma^{S,(k)}_{h}$ for $k=1,2$

  * •

Method 3 \- box-shaped internal discretization (Mehrstellen) (4) and
$\Gamma_{h}=\Gamma^{B}_{h}$ with $(r_{x},r_{y})=(1,1)$ for all ghost points

In summary, the linear system $A{\mathbf{u}}={\mathbf{b}}$ is assembled by the
Algorithm 1.

Initialize the matrix $A$ as the identity matrix of order $(n_{x}+1)(n_{y}+1)$
and ${\mathbf{b}}$ as the zero column vector

for $(x_{i},y_{j})\in\Omega_{h}$ do

$k=i+jn_{x}$

if Method 1 or Method 2 then

$A(k,[k,k-1,k+1,k-n_{x},k+n_{x},k-2,k+2,k-2n_{x},k+2n_{x}])$ ” ”
$=[60,-16,-16,-16,-16,1,1,1,1]/(12h^{2})$

$b(k)=f_{i,j}$

else if Method 3 then

$A(k,[k,k-1,k+1,k-n_{x},k+n_{x},k-1-n_{x},k-1+n_{x},k+1-n_{x},k+1+n_{x}])$ ” ”
$=[20,-4,-4,-4,-4,-1,-1,-1,-1]/(6h^{2})$

$b(k)=(8f_{i,j}+f_{i+1,j}+f_{i-1,j}+f_{i,j+1}+f_{i,j-1})/12$

end if

end for

for $(x_{l},y_{m})\in\Gamma_{h}$ do

Find the boundary point $B_{l,m}\equiv(x_{B},y_{B})$ and $s_{x},s_{y}$ by the
method described in Sect. 3.1

$\vartheta_{x}=s_{x}(x_{B}-x_{l})/(r_{x}\,h)$

$\vartheta_{y}=s_{y}(y_{B}-y_{m})/(r_{y}\,h)$

$k=l+m\,n_{x}$

if Method 2 and $(x_{l},y_{m})\in\Gamma^{S,(2)}_{h}$ then

$(r_{x},r_{y})=(2,2)$

else

$(r_{x},r_{y})=(1,1)$

end if

for $l=1,\ldots,5$ do

for $m=1,\ldots,5$ do

$\texttt{stencil}(l+5(m-1))=k+s_{x}r_{x}(l-1)+s_{y}r_{y}(m-1)$

end for

end for

$\texttt{weights}_{x}=[c_{1}(\vartheta_{x}),\ldots,c_{5}(\vartheta_{x})],\quad\texttt{weights}_{y}=[c_{1}(\vartheta_{y}),\ldots,c_{5}(\vartheta_{y})]$

if $B_{l,m}\in\Gamma_{D}$ then

$b(k)=g_{D}(B_{l,m})$

${\tt coeffs}={\tt weights}_{x}\cdot({\tt weights}_{y})^{\top}$

else if $B_{l,m}\in\Gamma_{N}$ then

$b(k)=g_{N}(B_{l,m})$

$\texttt{weights}_{dx}=[c^{\prime}_{1}(\vartheta_{x}),\ldots,c^{\prime}_{5}(\vartheta_{x})]/r_{x},\quad\texttt{weights}_{dy}=[c^{\prime}_{1}(\vartheta_{y}),\ldots,c^{\prime}_{5}(\vartheta_{y})]/r_{y}$

${\tt coeffs}_{dx}=s_{x}{\tt weights}_{y}\cdot({\tt
weights}_{dx})^{\top},\quad{\tt coeffs}_{dy}=s_{y}{\tt weights}_{x}\cdot({\tt
weights}_{dy})^{\top}$

$n_{x}=\phi({\tt stencil})\cdot({\tt coeffs}_{dx})^{\top},\quad
n_{y}=\phi({\tt stencil})\cdot({\tt coeffs}_{dy})^{\top}$

$n_{x}=n_{x}/\sqrt{n_{x}^{2}+n_{y}^{2}},\quad
n_{y}=n_{y}/\sqrt{n_{x}^{2}+n_{y}^{2}}$

${\tt coeffs}=n_{x}{\tt coeffs}_{dx}+n_{y}{\tt coeffs}_{dy}$

end if

$A(k,{\tt stencil})={\tt coeffs}$

end for

Algorithm 1 Pseudocode to assemble the matrix $A$ and the right-hand side $b$
of the linear system (Sect. 3.3)

Finally, we adopt a multigrid method to solve $A{\mathbf{u}}={\mathbf{b}}$,
similarly to the one described in [10] for second order accuracy. In this
paper we are interersted in investigating the accuracy of the discretization
and the conditioning of the system. The efficiency of the solver is not within
the scope of this study, although we experienced a similar performance to what
is observed for the second order accurate method in [10, 11].

##  4 Numerical results

We perform numerical tests on two different domains (see Fig. 6)): a circle
and a flower-shaped domain. The level-set representations $\phi_{C}$ and
$\phi_{F}$ of the two domains are:

$\text{\sc{(Circle)}}\quad\phi_{C}(x,y)=\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}}-R\text{
with }\;x_{0}=\sqrt{2}/10,\;y_{0}=-\sqrt{3}/20,\;R=\sqrt{5}/3$ (16)
$\text{\sc{(Flower)}}\quad\phi_{F}(x,y)=r-r_{1}-\frac{r_{2}}{r^{5}}\left((y-y_{0})^{5}+5(x-x_{0})^{4}(y-y_{0})-10(x-x_{0})^{2}(y-y_{0})^{3}\right),\\\
\text{ with }r=\sqrt{x^{2}+y^{2}},\quad r_{1}=0.5,\quad r_{2}=0.2.$ (17)

The flower-shaped domain can also be described by the parametric
representation:

$\left\\{\begin{aligned} x(\vartheta)&=r(\gamma)\cos(\gamma)+x_{0}\\\
y(\vartheta)&=r(\gamma)\sin(\gamma)+y_{0}\end{aligned}\right.$

with $\gamma\in[0,2\pi]$ and $r(\gamma)=r_{1}+r_{2}\sin(5\gamma)$.

Fig. 6:  Domains used in numerical tests (Sect. 4).

For each test we choose an exact solution $u_{\text{EXA}}(x,y)$ and then
compute analitically $f$, $g_{D}$ and $g_{N}$ from (1). We observe that the
analytical expression of $f(x,y)$ can be used to define grid values outside
the domain. However, we assume that $f_{ij}$ is only known on internal grid
points and then we extrapolate it outside $\Omega$ according to Remark 4.

The partition of the boundary into $\Gamma_{D}$ and $\Gamma_{N}$ is performed
in the following way in all numerical tests: a boundary point of
$\partial\Omega$ is in $\Gamma_{D}$ if $x\leq 0$, otherwise it is in
$\Gamma_{N}$.

The relative numerical error is computed as

$e_{h}=\frac{\left\|u_{h}-u_{\text{EXA}}\right\|_{p}}{\left\|u_{\text{EXA}}\right\|_{p}}=\left(\frac{\sum_{(x_{i},y_{j})\in\Omega_{h}\cup\Gamma_{h}}\left|u_{i,j}-u_{\text{EXA}}(x_{i},y_{j})\right|^{p}}{\sum_{(x_{i},y_{j})\in\Omega_{h}\cup\Gamma_{h}}\left|u_{\text{EXA}}(x_{i},y_{j})\right|^{p}}\right)^{1/p},\quad\text{
with }p=1,\infty,$

where $u_{h}$ is the numerical solution defined on the grid nodes of
$\Omega_{h}\cup\Gamma_{h}$. We also approximate the gradient of the solution
by a fourth-order accurate finite-difference formula:

$\nabla u_{h}=\left(\frac{\partial u_{h}}{\partial x},\frac{\partial
u_{h}}{\partial
y}\right)=\left(\frac{-u_{i+2,j}+8u_{i+1,j}-8u_{i-1,j}+u_{i-2,j}}{12h},\frac{-u_{i,j+2}+8u_{i,j+1}-8u_{i,j-1}+u_{i,j-2}}{12h}\right)+O(h^{4}).$
(18)

Then, we compute the error on the gradient of the solution by:

$e^{\nabla u}_{h}=\frac{\left\|\nabla u_{h}-\nabla
u_{\text{EXA}}\right\|_{p}}{\left\|\nabla
u_{\text{EXA}}\right\|_{p}}=\left(\frac{\sum_{(x_{i},y_{j})\in\Omega_{h}}\left|\nabla
u_{i,j}-\nabla
u_{\text{EXA}}(x_{i},y_{j})\right|^{p}}{\sum_{(x_{i},y_{j})\in\Omega_{h}}\left|\nabla
u_{\text{EXA}}(x_{i},y_{j})\right|^{p}}\right)^{1/p},\quad\text{ with
}p=1,\infty,$

where $\left|(v_{1},v_{2})\right|=\sqrt{v_{1}^{2}+v_{2}^{2}}$ is the modulus
of a two-dimensional vector.

###  4.1  Method 1

We choose $u_{\text{EXA}}=\sin(2x)\sin(5y)$ for the circular domain and
$u_{\text{EXA}}=\log(1+3xy)$ for the flower-shaped domain. The errors are
plotted in Fig. 7, showing that the method is not consistent, namely the
numerical error does not even converge to the exact solution when $N$ goes to
infinity. Indeed, it is an increasing function of $N$. This behavior is
explained by the condition number of the matrix $A$, that is very large (close
to $10^{16}$ and then almost no digits are accurate in the numerical
solution), resulting in an ill-conditioned linear system, as shown in Fig. 8.
This is caused by the presence of possible ill-conditioned interpolations. For
example, in (10) for Dirichlet boundary conditions (14) related to the second
layer of ghost-points $\Gamma^{S,(2)}_{h}$, if we solve (14) for the ghost
value $u_{l,m}$, we obtain

$u_{l,m}=\frac{1}{c_{0}(\vartheta_{x})c_{0}(\vartheta_{y})}\left(g_{D}(B_{l,m})-\sum_{\begin{subarray}{c}(p,q)\in\left\\{0,\ldots,4\right\\}^{2}\\\
(p,q)\neq(0,0)\end{subarray}}c_{p}(\vartheta_{x})c_{q}(\vartheta_{y})u_{l+s_{x}r_{x}p,m+s_{y}r_{y}q}\right).$

This is ill-conditioned when $c_{0}(\vartheta_{x})c_{0}(\vartheta_{y})\approx
0$, namely when $\vartheta_{x}$ or $\vartheta_{y}$ are close to 1 (see (9)).
This is observed, for example, in Fig. 4 for the purple ghost point (that
belongs to $\Gamma^{S,(2)}_{h}$), since its boundary projection is close to
the grid point $(x_{l-1},y_{m-1})$. To overcome this issue, Method 2 enlarges
the stencil for ghost points of $\Gamma^{S,(2)}_{h}$ by using
$(r_{x},r_{y})=(2,2)$ in (8). In this way, $\vartheta_{x}$ and $\vartheta_{y}$
are smaller than $1/2$ (from (12)) and the interpolation is always well-
conditioned.

The ill-conditioned behavior can also be explained in terms of small diagonal
terms of the matrix $A$ (compared to off-diagonal terms of the same row). In
fact, in the case of Dirichlet boundary condition, the diagonal term is
$c_{0}(\vartheta_{x})c_{0}(\vartheta_{y})$ and the off-diagonal terms are
$c_{p}(\vartheta_{x})c_{q}(\vartheta_{y})$ for
$(p,q)\in\left\\{0,\ldots,4\right\\}^{2}$, $(p,q)\neq(0,0)$. If
$\vartheta_{x}$ or $\vartheta_{y}$ are close to 1, then the diagonal term
$c_{0}(\vartheta_{x})c_{0}(\vartheta_{y})\approx 0$.

We note that the issue is observed in several methods where the domain is
embedded in a fixed grid that does not conform the boundary, falling under the
category of unfitted boundary methods, such as other ghost point methods [18]
or cut-cell methods (where the issue is caused by the presence of small cut
cells [23]).

Fig. 7:  Method 1 : numerical error of the solution and the gradient in
$L^{1}$ and $L^{\infty}$ norms vs different values of $N$, for both the
circular (left) and the flower-shaped (right) domains. The slope of the
bestfit lines provide an approximation of the accuracy orders, shown in the
legend.

Fig. 8:  Method 1 : condition number of the matrix $A$ for different values of
$N$ for the circular (left) and the flower-shaped (right) domains. Reference
slopes of $O(N^{2})$ and $O(N^{8})$ are represented by dashed and dotted
lines, respectively.

###  4.2  Method 2

We perform the same tests as in 4.1, except that we use en enlarged stencil
for the ghost points of the second layer $\Gamma^{S,(2)}_{h}$, namely
$(r_{x},r_{y})=(2,2)$ in (8), so that $\vartheta_{x}$ and $\vartheta_{y}$ are
always smaller than $1/2$ and then $c_{0}(\vartheta_{x})c_{0}(\vartheta_{y})$
is not close to zero. We observe that the method is fourth-order accurate in
the solution and its gradient (Fig. 9) and the condition number is improved
with respect to 4.1, scaling as $\kappa(A)=O(N^{2})$ (Fig. 11), that is the
expected scaling of the discrete Laplacian matrix, meaning that the ghost
point approach does not affect the overall stability of the method.

Fig. 9:  Method 2 : numerical error of the solution and the gradient in
$L^{1}$ and $L^{\infty}$ norms vs different values of $N$, for both the
circular (left) and the flower-shaped (right) domains. The slope of the
bestfit lines provide an approximation of the accuracy orders, shown in the
legend.

###  4.3  Method 3

We implement the box-shaped internal discretization (Mehrstellen) and plot the
errors in Fig. 10. We observe that the central difference formula based on a
stencil centered in $(x_{i},y_{j})$ to approximate the gradient (18) cannot be
computed for those internal grid points that are close to the boundary, since
the second layer of external points is made by inactive points in this method.
Therefore, for those points we use a different fourth-order formula by
selecting a different stencil, maintaining the same accuracy. For example, the
discretization of ${\partial u_{h}}/{\partial x}$ on the grid point
$(x_{i},y_{j})$ of Fig. 3 is:

$\frac{\partial u_{h}}{\partial
x}=\frac{-u_{i-3,j}+6u_{i-2,j}-18u_{i-1,j}+10u_{i,j}+3u_{i+1,j}}{12h}+O(h^{4}).$

instead of the one used in (18). Like Method 2 of Sect. 4.2, Method 3 is
fourth-order accurate in the solution and its gradient. The linear system is
well conditioned and $\kappa(A)=O(N^{2})$ as shown in Fig. 11.

Fig. 10:  Method 3 : numerical error of the solution and the gradient in
$L^{1}$ and $L^{\infty}$ norms vs different values of $N$, for both the
circular (left) and the flower-shaped (right) domains. The slope of the
bestfit lines provide an approximation of the accuracy orders, shown in the
legend.

Fig. 11:  Method 2 and Method 3 : condition number of the matrix $A$ for
different values of $N$ for the circular (left) and the flower-shaped (right)
domains. A reference slope of $O(N^{2})$ is represented by a dashed line.

###  4.4 Comparison of the methods

Method 1 is not consistent due to the ill-conditioned matrix, and therefore is
not further discussed in this section. Method 2 and Method 3 are both fourth-
order accurate in the solution and its gradient, say that the errors are
asymptotically scaling as $C_{2}h^{4}$ and $C_{3}h^{4}$ for Method 2 and
Method 3 respectively. However, Method 3 is considerably more accurate than
Method 2 meaning that $C_{3}<C_{2}$. Condition number of the matrix $A$ is
very similar for both methods and scales as expected, $\kappa(A)=O(N^{2})$,
although Method 3 is slightly more stable.

The improved accuracy and condition of Method 3 are justified by the adoption
of a better stencil for the discretization of the Laplace operator, since the
box-stencil is more compact than the star-stencil. This is also preferred by
the ghost-point approach, since ghost values are only required in the vicinity
of the boundary, while the star-shaped stencil needs the definition of grid
values on a farther second layer of ghost points.

The main drawback of Method 3 lies in the requirement that the source term $f$
must be evaluated on a $5$-point star-stencil (Eq. (4)), as highlighted in
Remark 4. For internal grid points that are close to the boundary, this means
that $f$ must be evaluated also on external grid points, therefore a suitable
technique must be adopted to extrapolate $f$ from internal to external points.
The ghost-point approach adopted in this paper to define the ghost values of
the solution outside the domain cannot be employed to extrapolate $f$, since
no boundary conditions are available for the source term. In this paper, the
artificial solution $u_{\text{EXA}}$ automatically defines an analytical
expression for the source term $f$ that can be used to evaluate the source
terms in points outside of the domain. In real applications this might not
possible, since usually only internal grid values are available. Therefore,
one has to adopt some other techniques to extrapolate $f$ from internal to
external grid points to the desired accuracy order. To this purpose, a very
efficient approach is the one proposed by Aslam in [1], where grid values can
be extrapolated outside of a domain with arbitrary geometry (defined by a
level-set function) up to the desired accuracy.

##  5 Conclusion

In this paper fourth order extensions of the finite difference ghost-point
method introduced in [10] have been presented. The methods are based on two
main ingredients:

  * •

finite difference discretization on the internal nodes;

  * •

interpolation based on an upwind stencil to apply the boundary condition at
the orthogonal projection of each ghost point onto the boundary.

To attain fourth-order accuracy not only for the function $u$ but also for its
gradient, it is necessary to employ biquartic interpolation. This is because
bicubic interpolation only ensures fourth-order accuracy for the function and
third-order accuracy for the gradient.

Three methods, which are theoretically fourth-order accurate for both the
function and its gradient, are considered. However, one of these methods
exhibits such poor conditioning that it fails to converge. The remaining two
methods demonstrate fourth-order accuracy for both the function and the
gradient. Among these, the method utilizing Mehrstellen discretization
achieves superior accuracy and is more practical, as it requires only a single
layer of ghost points.

This paper primarily examines the accuracy of the method while setting aside
several critical aspects, such as efficiency and a detailed convergence
analysis of the schemes. Concerning efficiency, a key advantage highlighted in
[10] is the compatibility of the method with a geometric multigrid approach
for solving the linear system, which significantly enhances its efficiency. In
[9] the efficiency of the multigrid method has been improved by optimizing the
relaxation parameters used during ghost value iteration. This optimization
specifically aims to smooth the residual in the tangential direction along the
boundary.

The convergence analysis for rectangular domains whose boundaries do not align
with Cartesian grid lines has been investigated in [8].

We expect that similar strategies regarding efficiency and convergence
analysis can be effectively applied to the fourth-order ghost method discussed
in this paper. Efforts to pursue these enhancements are actively underway.

## Acknowledgement

This work has been supported by the Spoke 1 Future HPC & Big Data of the
Italian Research Center on High-Performance Computing, Big Data and Quantum
Computing (ICSC) funded by MUR Missione 4 Componente 2 Investimento 1.4:
Potenziamento strutture di ricerca e creazione di “campioni nazionali di R &S
(M4C2-19)” - Next Generation EU (NGEU).

The work of A.C. has been supported from Italian Ministerial grant PRIN 2022
“Efficient numerical schemes and optimal control methods for time-dependent
partial differential equations”, No. 2022N9BM3N - Finanziato dall’Unione
europea - Next Generation EU – CUP: E53D23005830006.

The work of A.C. and G.R. has been supported from Italian Ministerial grant
PRIN 2022 PNRR “FIN4GEO: Forward and Inverse Numerical Modeling of
hydrothermal systems in volcanic regions with application to geothermal energy
exploitation”, No. P2022BNB97 - Finanziato dall’Unione europea - Next
Generation EU – CUP: E53D23017960001.

A.C. and G.R. are members of the Gruppo Nazionale Calcolo Scientifico-Istituto
Nazionale di Alta Matematica (GNCS-INdAM).

We acknowledge the CINECA award under the ISCRA initiative (ISCRA C project
BCMG, code HP10C7YOPZ), for the availability of high-performance computing
resources and support.

## References

  * [1] T. D. Aslam.  A partial differential equation approach to multidimensional extrapolation.  Journal of Computational Physics, 193:349–355, 2003.
  * [2] C. Astuto, D. Boffi, G. Russo, and U. Zerbinati.  A nodal ghost method based on variational formulation and regular square grid for elliptic problems on arbitrary domains in two space dimensions.  arXiv preprint arXiv:2402.04048, 2024.
  * [3] A. Baeza, P. Mulet, and D. Zorío.  High order boundary extrapolation technique for finite difference methods on complex domains with cartesian meshes.  Journal of Scientific Computing, 66:761–791, 2016.
  * [4] F. Chantalat, C.-H. Bruneau, C. Galusinski, and A. Iollo.  Level-set, penalization and cartesian meshes: A paradigm for inverse problems and optimal design.  Journal of Computational Physics, 228:6291–6315, 2009.
  * [5] A. Chertock, A. Coco, A. Kurganov, and G. Russo.  A second-order finite-difference method for compressible fluids in domains with moving boundaries.  Commun Comput Phys, 23(1), 2018.
  * [6] A. Coco.  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.
  * [7] A. Coco, G. Currenti, C. Del Negro, and G. Russo.  A second order finite-difference ghost-point method for elasticity problems on unbounded domains with applications to volcanology.  Communications in Computational Physics, 16(4):983–1009, 2014.
  * [8] A. Coco, S.-E. Ekström, G. Russo, S. Serra-Capizzano, and S. C. Stissi.  Spectral and norm estimates for matrix-sequences arising from a finite difference approximation of elliptic operators.  Linear Algebra and its Applications, 667:10–43, 2023.
  * [9] A. Coco, M. Mazza, and M. Semplice.  A ghost-point smoothing strategy for geometric multigrid on curved boundaries.  Journal of Computational Physics, 478:111982, 2023.
  * [10] A. Coco and G. Russo.  Finite-Difference Ghost-Point Multigrid Methods on Cartesian Grids for Elliptic Problems in Arbitrary Domains.  Journal of Computational Physics, 241:464–501, 2013.
  * [11] A. Coco and G. Russo.  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.
  * [12] L. Collatz.  The numerical treatment of differential equations, volume 60.  Springer Science & Business Media, 2012.
  * [13] A. Du Chéné, C. Min, and F. Gibou.  Second-order accurate computation of curvatures in a level set framework using novel high-order reinitialization schemes.  Journal of Scientific Computing, 35:114–131, 2008.
  * [14] Q. Feng, B. Han, and P. Minev.  A high order compact finite difference scheme for elliptic interface problems with discontinuous and high-contrast coefficients.  Applied Mathematics and Computation, 431:127314, 2022.
  * [15] J. Fernández-Fidalgo, S. Clain, L. Ramírez, I. Colominas, and X. Nogueira.  Very high-order method on immersed curved domains for finite difference schemes with regular cartesian grids.  Computer Methods in Applied Mechanics and Engineering, 360:112782, 2020.
  * [16] J. Gabbard and W. M. van Rees.  A high-order finite difference method for moving immersed domain boundaries and material interfaces.  Journal of Computational Physics, page 112979, 2024.
  * [17] F. Gibou and R. Fedkiw.  A second-order-accurate symmetric discratization of the poisson equation on irregular domains.  Journal of Computational Physics, 176:205–227, 2002.
  * [18] F. Gibou, R. P. Fedkiw, L.-T. Cheng, and M. Kang.  A second-order-accurate symmetric discretization of the poisson equation on irregular domains.  Journal of Computational Physics, 176(1):205–227, 2002.
  * [19] Y. T. Ng, H. Chen, C. Min, and F. Gibou.  Guidelines for Poisson solvers on irregular domains with Dirichlet boundary conditions using the ghost fluid method.  J. Sci. Comput., 41:300–320, 2009.
  * [20] S. Osher and R. Fedkiw.  Level Set Methods and Dynamic Implicit Surfaces.  Springer-Verlag New York, Applied Mathematical Sciences, 2002.
  * [21] Y. Ren, H. Feng, and S. Zhao.  A fft accelerated high order finite difference method for elliptic boundary value problems over irregular domains.  Journal of Computational Physics, 448:110762, 2022.
  * [22] G. Russo and P. Smereka.  A remark on computing distance functions.  Journal of Computational Physics, 163:51–67, 2000.
  * [23] L. Schneiders, D. Hartmann, M. Meinke, and W. Schröder.  An accurate moving boundary formulation in cut-cell methods.  Journal of Computational Physics, 235:786–809, 2013.
  * [24] J. Sethian.  Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science.  Cambridge University Press, 1999.
  * [25] G. H. Shortley and R. Weller.  The numerical solution of laplace’s equation.  J. Appl. Phys., 9:334–348, 1938.
  * [26] M. Sussman, P. Smereka, and S. Osher.  A level set approach for computing solutions to incompressible 2-phase flow.  Journal of Computational Physics, 114:146–159, 1994.
  * [27] S. Yu, Y. Zhou, and G. Wei.  Matched Interface and Boundary (MIB) method for elliptic problems with sharp-edged interfaces.  Journal of Computational Physics, 224:729–756, 2007.