# On Solving the Poisson Equation with Discontinuities on Irregular Interfaces: GFM and VIM

Authors: Ásdís Helgadóttir; Arthur Guittet; Frédéric Gibou
Journal: International Journal of Differential Equations
 (2018)

---

## Abstract

We analyze the accuracy of two numerical methods for the variable coefficient Poisson equation with discontinuities at an irregular interface. Solving the Poisson equation with discontinuities at an irregular interface is an essential part of solving many physical phenomena such as multiphase flows with and without phase change, in heat transfer, in electrokinetics, and in the modeling of biomolecules’ electrostatics. The first method, considered for the problem, is the widely known Ghost-Fluid Method (GFM) and the second method is the recently introduced Voronoi Interface Method (VIM). The VIM method uses Voronoi partitions near the interface to construct local configurations that enable the use of the Ghost-Fluid philosophy in one dimension. Both methods lead to symmetric positive definite linear systems. The Ghost-Fluid Method is generally first-order accurate, except in the case of both a constant discontinuity in the solution and a constant diffusion coefficient, while the Voronoi Interface Method is second-order accurate in theL∞-norm. Therefore, the Voronoi Interface Method generally outweighs the Ghost-Fluid Method except in special case of both a constant discontinuity in the solution and a constant diffusion coefficient, where the Ghost-Fluid Method performs better than the Voronoi Interface Method. The paper includes numerical examples displaying this fact clearly and its findings can be used to determine which approach to choose based on the properties of the real life problem in hand.


---

## Body

## 1. Introduction

The Poisson equations with discontinuities across irregular interfaces emerge in applications such as multiphase flows with and without phase change, in heat transfer, in electrokinetics, or in the modeling of biomolecules’ electrostatics. Several numerical methods have been proposed to solve this system, each with their own advantages and disadvantages. One approach is in the context of Discontinuous Galerkin methods, an extension of the finite element method (FEM); e.g., see [1–12] and the references therein. Finite element methods lead to symmetric positive definite linear systems that can be efficiently solved with fast iterative solvers [13]. In addition, FEM-type approaches can derive and usea priori error estimates to refine the mesh where higher resolution is needed. However, FEM-type methods rely on the quality of the underlying mesh, which is often difficult to obtain in cases where the irregular domain undergoes large deformations. In this case, it is challenging to generate a mesh with elements that pass a quality measure needed to ensure accurate solutions.Differential quadrature method is worth mentioning as presented in [14, 15]. There Lagrange interpolation and/or modified cubic B-splines are used, depending on boundary conditions, to approximate the solution to two-dimensional nonlinear hyperbolic partial differential equations. Irregular interfaces and discontinuities were, however, not studied there.Another approach is within the context of finite difference methods (FDM); e.g., see [16–25] and the references therein. For finite difference methods the grid is Cartesian (uniform or adaptive), which leads to a straightforward grid generation process. However, interfaces must be represented by other means and the treatment of boundary conditions requires additional considerations.In order to capture the interface and enforce the correct boundary conditions, several approaches have been explored. The immersed boundary method (for example, [26–30]) uses the δ-formulation that smears out the solution profile across the interface. This produces algorithms that are straightforward to implement since they are similar to solving the same equations on a regular domain. However, the smearing of the solution introduces O(1) errors near the interface. The immerse interface method (IIM) (see, for example, [31–37]) is a sharp interface method that leads to second-order accurate solution, albeit it is not a robust second-order method since it reaches its order by minimizing the truncation error. The method leads to neither symmetric nor positive definite linear systems and the application of the immerse interface method may be difficult in three spatial dimensions.In [20], Liuet al. introduced a Ghost-Fluid Methodology (a finite difference approach) that treats the jump condition in a dimension-by-dimension framework. This leads to a linear system that is symmetric positive definite and the jump conditions only affect the right-hand side of the linear system, leading to an easy-to-implement method. However, the GFM suffers a loss in accuracy due to smearing of the tangential derivative of the discontinuity at the interface caused by this dimension-by-dimension framework, as indicated in [20].The Voronoi Interface Method (VIM) [38] avoids the loss of accuracy of the GFM of [20] by constructing local Voronoi partitioning near the interface. The cells adjacent to the interface therefore have their faces orthogonal to the fluxes of the solution, hence providing a configuration that can leverage the Ghost-Fluid Methodology to its fullest. The advantage is therefore that the solution is second-order accurate. A drawback is that the solution is computed at the cells’ center of the Voronoi partition and an additional interpolation step is required, if the solution is needed on the original Cartesian mesh. As is the case of GFM of [20], the linear system is symmetric positive definite and only its right-hand side is modified by the jump conditions. Finally, we note that even though the construction of a global Voronoi partition may be difficult and costly, the construction of a local partition, i.e., only for cells that are adjacent to the interface, is straightforward. In [38], the authors used the Voro++ library of [39].In this paper we highlight the performance of both GFM and VIM. It is worth noting that both methods are unconditionally stable since they are implicit. Both methods have also shown that they can handle high ratios of discontinuity in their coefficients.


## 2. Governing Equations and Numerical Methods

We consider the Poisson equation with variable coefficient and discontinuities (jumps) across an irregular interface,Γ, which splits the computational domain Ω into two domains Ω+ and Ω-, both in Rn, n∈N. The governing equation is (1)∇·β∇u=fforx_∈Ω+∪Ω-,where f=f(x_) and β=β(x_) are given. Here, β is bounded from below by a positive constant and f is in L2(Ω). This equation is supplemented by jump conditions on the irregular interface, (2)u=a,β∇u·n_=h,where u=uΓ+-uΓ- denotes a jump in u across Γ. The functions a=a(x_) and h=h(x_) are given. Either Dirichlet, Neumann, or Robin boundary conditions can be imposed at the boundaries of the computational domain, ∂Ω. In order to represent the irregular interface, we use the signed distance level-set function, ϕ, which is positive inside Ω+, negative inside Ω-, and zero on Γ. The outward normal to the irregular interface can be computed from the level-set function by (3)n_=n1,n2,n3=∇ϕ∇ϕ.


### 2.1. The Ghost-Fluid Method

A detailed description of the Ghost-Fluid Method is given in [20], we will explore here the main aspects in two spatial dimensions. The discontinuity in βun, where un refers to the derivative in the normal direction to the interface, is βunΓ=βuxΓn1+βuyΓn2. The discontinuity in βut, where ut refers to the derivative in the tangential direction to the interface, is βutΓ=βuxΓn2-βuyΓn1. These equations lead to βuxΓ=βunΓn1+βutΓn2 and βuyΓ=βunΓn2-βutΓn1. However, in order to devise a method in a dimension-by-dimension framework, the Ghost-Fluid Method smears out the discontinuity in the tangential derivative leading to the simplification βuxΓ=βunΓn1 and βuyΓ=βunΓn2. The discretization at each grid point i,j is then given by (4)βi+1/2,jui+1,j-ui,j/Δx-βi-1/2,jui,j-ui-1,j/ΔxΔx+βi,j+1/2ui,j+1-ui,j/Δy-βi,j-1/2ui,j-ui,j-1/ΔyΔy=fi,j+Fx+Fy,where βi+1/2,j=βi+βi+1/2 if xi,j and xi+1,j are on the same side of the interface or by βi+1/2,j=β+β-ϕ-+ϕ+/β+ϕ-+β-ϕ+ otherwise, where ϕ+ is the value of ϕ of the node in Ω+ and ϕ- is the value of ϕ of the node in Ω-. Here β± refers to the value of β adjacent to the interface in the Ω± domain and Δx and Δy are the cells’ sizes in the x- and y-directions, respectively. The left-hand side thus gives the same symmetric positive definite matrix as the one generated by the standard five-point stencil discretization of the Poisson equation on regular domains and only the right-hand side is altered when a discontinuity occurs. Furthermore, Fx=FL+FR and Fy=FB+FT are only activated if there is a discontinuity present in the local five point stencil. Here, FL and FR and the contribution from the left and right grid points to the current grid point and FB and FT are the contributions from the bottom and top grid points.We give the details forFL and FR in the x-direction, referring the reader to the original paper [20] for the description for FB and FT in the y-direction, since they follow the same procedure. If ϕi-1,j and ϕi,j have opposite signs, define θ=|ϕi-1,j|/ϕi,j+ϕi-1,j, aΓ=ai,jϕi-1,j+ai-1,jϕi,j/ϕi,j+ϕi-1,j, where ai,j is the jump at the grid node and βΓ=bi,jni,j1ϕi-1,j+bi-1,jni-1,j1ϕi,j/ϕi,j+ϕi-1,j. FL is then defined as (5)FL=βi-1/2,jaΓΔx2-βi-1/2,jbΓθβ+Δxifϕi,j≤0andϕi-1,j>0,-βi-1/2,jaΓΔx2+βi-1/2,jbΓθβ-Δxifϕi,j>0andϕi-1,j≤0,0otherwise.Similarly, if ϕi+1,j and ϕi,j have opposite signs, define θ=|ϕi+1,j|/ϕi,j+ϕi+1,j, aΓ=ai,jϕi+1,j+ai+1,jϕi,j/ϕi,j+ϕi+1,j, and βΓ=bi,jni,j1ϕi+1,j+bi+1,jni+1,j1ϕi,j/ϕi,j+ϕi+1,j. FR is then defined as (6)FR=βi+1/2,jaΓΔx2+βi+1/2,jbΓθβ+Δxifϕi,j≤0andϕi+1,j>0,-βi+1/2,jaΓΔx2-βi+1/2,jbΓθβ-Δxifϕi,j>0andϕi+1,j≤0,0otherwise.The Ghost-Fluid Method leads to a symmetric positive definite linear systems that captures the discontinuity in the normal derivative while smearing out the discontinuity in the tangential direction to the interface.


### 2.2. The Voronoi Interface Method

We present a summary of the method and refer the reader to [38] for a detailed description. For a given set of seeds, the Voronoi cell of a seed is defined as the points of space that are closer to that seed than any other. The union of the Voronoi cells is a tessellation of space, and the first step for the Voronoi Interface Method is to generate the Voronoi mesh associated with the background mesh chosen, which in our case is a uniform Cartesian grid. The seeds are of two types,(i)
if a cell of the background mesh is not crossed by the irregular interfaceΓ, its center is a seed of the Voronoi mesh,(ii)
if a cell of the background mesh is crossed by the irregular interfaceΓ, we locate the projection of its center onto the interface and generate two seeds located on either side of the projected point at a distance δ=min(Δx,Δy)/5.The Voronoi partition associated with those seeds is constructed using a simple geometric algorithm. We refer the reader to theVoro++ library [39] for an efficient tool to compute Voronoi partitions in both two and three spatial dimensions.Equation (1) is then discretized on the Voronoi mesh using a finite volume approach. The complete derivation is presented in [38] and leads to the discretization of the interaction between point i and point j(7)β~ijsijuj-uidij=β~ijsijdij-sign⁡ϕiu+dij2βjβ∇u·n_Γ+VolCifi,where any variable of the form γk would be the quantity γ at point k, dij is the distance between i and j, sij is the length of the face between i and j, β~ij=2βiβj/βi+βj, and Vol(Ci) is the volume of the Voronoi cell associated with point i. The discontinuities at the interface only affect the right-hand side of the linear system where the irregular interface is located and the system is symmetric positive definite.We remark that the solution is provided at the center of the Voronoi cells. If the solution is needed on the original background mesh, then an interpolation step is required. This can be done, for example, with least square interpolations and it does not impact the order of accuracy of the solution as long as the order of the polynomial interpolants is high enough. In this article, we work with second-order polynomial interpolants.


## 2.1. The Ghost-Fluid Method

A detailed description of the Ghost-Fluid Method is given in [20], we will explore here the main aspects in two spatial dimensions. The discontinuity in βun, where un refers to the derivative in the normal direction to the interface, is βunΓ=βuxΓn1+βuyΓn2. The discontinuity in βut, where ut refers to the derivative in the tangential direction to the interface, is βutΓ=βuxΓn2-βuyΓn1. These equations lead to βuxΓ=βunΓn1+βutΓn2 and βuyΓ=βunΓn2-βutΓn1. However, in order to devise a method in a dimension-by-dimension framework, the Ghost-Fluid Method smears out the discontinuity in the tangential derivative leading to the simplification βuxΓ=βunΓn1 and βuyΓ=βunΓn2. The discretization at each grid point i,j is then given by (4)βi+1/2,jui+1,j-ui,j/Δx-βi-1/2,jui,j-ui-1,j/ΔxΔx+βi,j+1/2ui,j+1-ui,j/Δy-βi,j-1/2ui,j-ui,j-1/ΔyΔy=fi,j+Fx+Fy,where βi+1/2,j=βi+βi+1/2 if xi,j and xi+1,j are on the same side of the interface or by βi+1/2,j=β+β-ϕ-+ϕ+/β+ϕ-+β-ϕ+ otherwise, where ϕ+ is the value of ϕ of the node in Ω+ and ϕ- is the value of ϕ of the node in Ω-. Here β± refers to the value of β adjacent to the interface in the Ω± domain and Δx and Δy are the cells’ sizes in the x- and y-directions, respectively. The left-hand side thus gives the same symmetric positive definite matrix as the one generated by the standard five-point stencil discretization of the Poisson equation on regular domains and only the right-hand side is altered when a discontinuity occurs. Furthermore, Fx=FL+FR and Fy=FB+FT are only activated if there is a discontinuity present in the local five point stencil. Here, FL and FR and the contribution from the left and right grid points to the current grid point and FB and FT are the contributions from the bottom and top grid points.We give the details forFL and FR in the x-direction, referring the reader to the original paper [20] for the description for FB and FT in the y-direction, since they follow the same procedure. If ϕi-1,j and ϕi,j have opposite signs, define θ=|ϕi-1,j|/ϕi,j+ϕi-1,j, aΓ=ai,jϕi-1,j+ai-1,jϕi,j/ϕi,j+ϕi-1,j, where ai,j is the jump at the grid node and βΓ=bi,jni,j1ϕi-1,j+bi-1,jni-1,j1ϕi,j/ϕi,j+ϕi-1,j. FL is then defined as (5)FL=βi-1/2,jaΓΔx2-βi-1/2,jbΓθβ+Δxifϕi,j≤0andϕi-1,j>0,-βi-1/2,jaΓΔx2+βi-1/2,jbΓθβ-Δxifϕi,j>0andϕi-1,j≤0,0otherwise.Similarly, if ϕi+1,j and ϕi,j have opposite signs, define θ=|ϕi+1,j|/ϕi,j+ϕi+1,j, aΓ=ai,jϕi+1,j+ai+1,jϕi,j/ϕi,j+ϕi+1,j, and βΓ=bi,jni,j1ϕi+1,j+bi+1,jni+1,j1ϕi,j/ϕi,j+ϕi+1,j. FR is then defined as (6)FR=βi+1/2,jaΓΔx2+βi+1/2,jbΓθβ+Δxifϕi,j≤0andϕi+1,j>0,-βi+1/2,jaΓΔx2-βi+1/2,jbΓθβ-Δxifϕi,j>0andϕi+1,j≤0,0otherwise.The Ghost-Fluid Method leads to a symmetric positive definite linear systems that captures the discontinuity in the normal derivative while smearing out the discontinuity in the tangential direction to the interface.


## 2.2. The Voronoi Interface Method

We present a summary of the method and refer the reader to [38] for a detailed description. For a given set of seeds, the Voronoi cell of a seed is defined as the points of space that are closer to that seed than any other. The union of the Voronoi cells is a tessellation of space, and the first step for the Voronoi Interface Method is to generate the Voronoi mesh associated with the background mesh chosen, which in our case is a uniform Cartesian grid. The seeds are of two types,(i)
if a cell of the background mesh is not crossed by the irregular interfaceΓ, its center is a seed of the Voronoi mesh,(ii)
if a cell of the background mesh is crossed by the irregular interfaceΓ, we locate the projection of its center onto the interface and generate two seeds located on either side of the projected point at a distance δ=min(Δx,Δy)/5.The Voronoi partition associated with those seeds is constructed using a simple geometric algorithm. We refer the reader to theVoro++ library [39] for an efficient tool to compute Voronoi partitions in both two and three spatial dimensions.Equation (1) is then discretized on the Voronoi mesh using a finite volume approach. The complete derivation is presented in [38] and leads to the discretization of the interaction between point i and point j(7)β~ijsijuj-uidij=β~ijsijdij-sign⁡ϕiu+dij2βjβ∇u·n_Γ+VolCifi,where any variable of the form γk would be the quantity γ at point k, dij is the distance between i and j, sij is the length of the face between i and j, β~ij=2βiβj/βi+βj, and Vol(Ci) is the volume of the Voronoi cell associated with point i. The discontinuities at the interface only affect the right-hand side of the linear system where the irregular interface is located and the system is symmetric positive definite.We remark that the solution is provided at the center of the Voronoi cells. If the solution is needed on the original background mesh, then an interpolation step is required. This can be done, for example, with least square interpolations and it does not impact the order of accuracy of the solution as long as the order of the polynomial interpolants is high enough. In this article, we work with second-order polynomial interpolants.


## 3. Numerical Experiments: Results

We present a pair of two-dimensional numerical examples where the two methods are compared. Both examples have a star shaped irregular interface. The first example has a coefficientβ that is constant in the whole domain. This example has a constant discontinuity in the solution across the irregular interface but neither a discontinuity in the normal nor tangential derivatives at the interface. The second example has a β coefficient that is not constant and has a discontinuity across the irregular interface. This example has a nonconstant discontinuity in the solution and nonconstant discontinuities in the normal and tangential components of the gradient of the solution at the irregular interface.


### 3.1. Constantβ-Coefficient

Let us consider∇·β∇u=fx,y in two spatial dimensions in Ω=-1,12 with the level set function ϕ=-x2+y2+0.5+0.15cos5arctany/x. We take β=1 when ϕ≤0 and an exact solution of u=cos(x)cos(y). For the region ϕ>0, we take β=1 and the exact solution to be u=cos(x)cos(y)+1. A representation of the solution and the Voronoi mesh are given in Figure 1. The comparison of the Ghost-Fluid Method and the Voronoi Interface Method is shown in Table 1 for the maximum error in the solution, in Table 2 for the average error in the solution, in Table 3 for the maximum error of the gradient, and in Table 4 for the average error of the gradient. For the Voronoi Interface Method, we present the errors on both the Voronoi mesh and the Cartesian mesh. For the Voronoi mesh, the gradients are computed on the faces of the Voronoi cells. For example, the solution only experiences a constant discontinuity in its solution and no discontinuities in its normal or tangential derivatives nor is there a discontinuity in the β coefficient, which is constant in the entire domain. Both methods give second-order accuracy in the solution in the L∞- and L1-norms and the solution’s gradient in the L1-norm. However, the Voronoi method gives first-order accuracy for the gradient of the solution in the L∞-norm while the GFM gives second-order accuracy for such simplified problems. A likely explanation for this difference is that the fluxes in the volume of fluid derivation for the Voronoi Interface Method are not necessarily computed at the center of the faces between two points. We also observe that the interpolation from the Voronoi mesh to the Cartesian mesh does not impact the solution itself but affects its gradient. However, the order of accuracy is conserved. We conclude that for such simple problems the GFM is preferable.Table 1
Comparison of theL∞ error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u-uh∞ Order u-uh∞ Order u-uh∞ Order  0.25 2.627×10-3 — 3.709 × 1 0 - 3 — 4.098×10-1 —  0.125 6.587 × 1 0 - 4 2.00 6.172 × 1 0 - 4 2.59 5.901×10-2 2.80  0.0625 1.648 × 1 0 - 4 2.00 1.565 × 1 0 - 4 1.98 1.788×10-4 8.37  0.03125 4.121 × 1 0 - 5 2.00 4.016 × 1 0 - 5 1.96 4.030×10-5 2.15  0.015625 1.030 × 1 0 - 5 2.00 1.020 × 1 0 - 5 1.98 1.023×10-5 1.98  0.0078125 2.576 × 1 0 - 6 2.00 2.598 × 1 0 - 6 1.97 2.656×10-6 1.95  0.00390625 6.439 × 1 0 - 7 2.00 6.591 × 1 0 - 7 1.98 6.655×10-7 2.00Table 2
Comparison of theL1 error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u - u h 1 Order u - u h 1 Order u - u h 1 Order  0.25 9.010 × 1 0 - 4 — 1.790 × 1 0 - 2 — 4.694 × 1 0 - 2 —  0.125 2.597 × 1 0 - 4 1.79 2.918 × 1 0 - 3 5.04 1.417 × 1 0 - 3 5.04  0.0625 6.940 × 1 0 - 5 1.90 7.365 × 1 0 - 5 4.31 7.130 × 1 0 - 5 4.31  0.03125 1.792 × 1 0 - 5 1.95 1.848 × 1 0 - 5 2.00 1.780 × 1 0 - 5 2.00  0.015625 4.552 × 1 0 - 6 1.98 4.625 × 1 0 - 6 1.98 4.518 × 1 0 - 6 1.98  0.0078125 1.147 × 1 0 - 6 1.99 1.156 × 1 0 - 6 1.98 1.142 × 1 0 - 6 1.98  0.00390625 2.879 × 1 0 - 7 1.99 2.890 × 1 0 - 7 1.99 2.871 × 1 0 - 7 1.99Table 3
Comparison of theL∞ error in the gradient of the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order  0.25 4.829 × 1 0 - 3 — 8.160 × 1 0 - 3 — 1.872 × 1 0 - 0 —  0.125 1.398 × 1 0 - 3 1.79 1.550 × 1 0 - 3 2.40 3.666 × 1 0 - 1 2.35  0.0625 3.776 × 1 0 - 4 1.89 5.192 × 1 0 - 4 1.58 4.180 × 1 0 - 2 3.13  0.03125 9.866 × 1 0 - 5 1.94 4.367 × 1 0 - 4 0.25 2.082 × 1 0 - 2 1.01  0.015625 2.530 × 1 0 - 5 1.96 2.551 × 1 0 - 4 0.78 1.037 × 1 0 - 2 1.00  0.0078125 6.419 × 1 0 - 6 1.98 1.376 × 1 0 - 4 0.89 5.195 × 1 0 - 3 0.99  0.00390625 1.618 × 1 0 - 6 1.99 8.580 × 1 0 - 5 0.68 2.597 × 1 0 - 3 1.00Table 4
Comparison of theL1 error in the gradient of the solution for the Ghost-Fluid method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order  0.25 3.636 × 1 0 - 3 — 2.575 × 1 0 - 3 — 3.797 × 1 0 - 1 —  0.125 8.704 × 1 0 - 4 2.06 3.490 × 1 0 - 4 2.88 2.908 × 1 0 - 2 3.71  0.0625 2.155 × 1 0 - 4 2.01 7.553 × 1 0 - 5 2.21 5.617 × 1 0 - 3 2.37  0.03125 5.343 × 1 0 - 5 2.01 1.711 × 1 0 - 5 2.14 1.474 × 1 0 - 3 1.93  0.015625 1.331 × 1 0 - 5 2.01 4.126 × 1 0 - 6 2.05 3.707 × 1 0 - 4 1.99  0.0078125 3.322 × 1 0 - 6 2.00 9.930 × 1 0 - 7 2.05 9.345 × 1 0 - 5 1.99  0.00390625 8.298 × 1 0 - 7 2.00 2.465 × 1 0 - 7 2.01 2.350 × 1 0 - 5 1.99Figure 1
Left: visualization of the exact solution, for example, shown in Section3.1.  Right: example of a Voronoi mesh generated.


### 3.2. Nonconstantβ-Coefficient

Let us consider∇·β∇u=fx,y in two spatial dimensions in Ω=-1,12 with the level set function ϕ=-x2+y2+0.5+0.15cos5arctany/x. We set β=y2lnx+2+4 and the exact solution to u=ex in the region where ϕ≤0. In the region where ϕ>0, we set β=e-y and the exact solution to u=cosxsiny. Figure 2 provides a representation of the solution and of the diffusion coefficient. The comparison of the Ghost-Fluid Method and the Voronoi Interface Method is shown in Table 5 for the maximum error in the solution, in Table 6 for the average error in the solution, in Table 7 for the maximum error of the gradient, and in Table 8 for the average error of the gradient. Again, for the Voronoi Interface Method we present the results on both the Voronoi and the Cartesian meshes. For example, there is a nonconstant discontinuity in the solution, a nonconstant discontinuity in the normal derivative of the solution to the interface, a nonconstant discontinuity in the tangential derivative of the solution to the interface, a discontinuity in the β coefficient across the interface, and a β coefficient that is not constant in each domain as in the previous example. The Ghost-Fluid Method is only first-order accurate in both the L∞- and the L1-norm, as well as in the L1-norm for the gradient of the solution. However, it is not consistent in the L∞-norm for the gradient of the solution. This example is a case where the lack of orthogonality between the cells’ faces and the solution’s fluxes lowers the accuracy of the GFM’s dimension-by-dimension approach. The Voronoi Interface Method provides a solution to that drawback and, therefore, produces a second-order accurate solution in the L∞- and L1-norms, and first-order accurate (resp., second-order accurate) gradients in the L∞- (resp., L1-) norm. Similarly to the previous example, we observe a decrease in accuracy for the gradient computed on the Cartesian mesh after the interpolation step, but the order of accuracy is conserved. The solution itself does not suffer from the interpolation step. For this type of complex problems, the Voronoi Interface Method is able to provide a second-order accurate solution and consistent gradient that the Ghost-Fluid Method cannot produce. VIM is therefore the recommended approach.Table 5
Comparions of theL∞ error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u - u h ∞ Order u - u h ∞ Order u - u h ∞ Order  0.25 6.254 × 1 0 - 1 — 9.802 × 1 0 - 3 — 6.044 × 1 0 - 1 —  0.125 4.063 × 1 0 - 1 0.62 2.560 × 1 0 - 3 1.94 9.178 × 1 0 - 2 2.72  0.0625 2.493 × 1 0 - 1 0.70 7.709 × 1 0 - 4 1.73 7.041 × 1 0 - 4 7.03  0.03125 1.474 × 1 0 - 1 0.76 2.370 × 1 0 - 4 1.70 2.160 × 1 0 - 4 1.70  0.015625 8.488 × 1 0 - 2 0.80 5.560 × 1 0 - 5 2.09 4.929 × 1 0 - 5 2.13  0.0078125 4.802 × 1 0 - 2 0.82 1.413 × 1 0 - 5 1.98 1.293 × 1 0 - 5 1.93  0.00390625 1.481 × 1 0 - 2 0.84 3.618 × 1 0 - 6 1.97 3.377 × 1 0 - 6 1.94Table 6
Comparison of theL1 error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u - u h 1 Order u - u h 1 Order u - u h 1 Order  0.25 6.446 × 1 0 - 2 — 2.291 × 1 0 - 3 — 4.557 × 1 0 - 2 —  0.125 3.393 × 1 0 - 2 0.93 4.022 × 1 0 - 4 2.51 2.041 × 1 0 - 3 4.48  0.0625 1.692 × 1 0 - 2 1.00 1.038 × 1 0 - 4 1.95 9.891 × 1 0 - 5 4.37  0.03125 7.908 × 1 0 - 3 1.10 2.985 × 1 0 - 5 1.80 2.718 × 1 0 - 5 1.86  0.015625 3.982 × 1 0 - 3 0.99 7.243 × 1 0 - 6 2.04 6.729 × 1 0 - 6 2.01  0.0078125 2.155 × 1 0 - 3 0.89 1.790 × 1 0 - 6 2.02 1.717 × 1 0 - 6 1.97  0.00390625 1.035 × 1 0 - 3 1.06 4.509 × 1 0 - 7 1.99 4.401 × 1 0 - 7 1.96Table 7
Comparison of theL∞ error in the gradient of the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order  0.25 1.259 — 1.309 × 1 0 - 2 — 2.658 × 1 0 - 0 —  0.125 1.227 0.04 9.931 × 1 0 - 3 0.55 6.334 × 1 0 - 1 2.07  0.0625 1.139 0.11 2.954 × 1 0 - 3 1.60 8.950 × 1 0 - 2 2.82  0.03125 1.156 -0.02 1.589 × 1 0 - 3 0.89 4.373 × 1 0 - 2 1.03  0.015625 1.188 -0.04 1.038 × 1 0 - 3 0.62 2.155 × 1 0 - 2 1.02  0.0078125 1.205 -0.02 5.177 × 1 0 - 4 1.00 1.070 × 1 0 - 2 1.01  0.00390625 1.214 -0.01 3.194 × 1 0 - 4 0.70 5.329 × 1 0 - 3 1.01Table 8
Comparison of theL1 error in the gradient of the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order  0.25 2.611 × 1 0 - 1 — 2.126 × 1 0 - 3 — 3.895 × 1 0 - 1 —  0.125 1.438 × 1 0 - 1 0.86 6.395 × 1 0 - 4 1.73 3.400 × 1 0 - 2 3.52  0.0625 7.136 × 1 0 - 2 1.01 1.957 × 1 0 - 4 1.71 5.631 × 1 0 - 3 2.59  0.03125 3.922 × 1 0 - 2 0.86 6.409 × 1 0 - 5 1.61 1.464 × 1 0 - 3 1.94  0.015625 2.019 × 1 0 - 2 0.96 1.593 × 1 0 - 5 2.01 3.713 × 1 0 - 4 1.98  0.0078125 1.033 × 1 0 - 2 0.97 4.132 × 1 0 - 6 1.95 9.359 × 1 0 - 5 1.99  0.00390625 5.222 × 1 0 - 3 0.98 1.099 × 1 0 - 6 1.91 2.358 × 1 0 - 5 1.99Figure 2
Visualization of the exact solutionu (left) and of the diffusion coefficient β (right), for example, shown in Section 3.2.


## 3.1. Constantβ-Coefficient

Let us consider∇·β∇u=fx,y in two spatial dimensions in Ω=-1,12 with the level set function ϕ=-x2+y2+0.5+0.15cos5arctany/x. We take β=1 when ϕ≤0 and an exact solution of u=cos(x)cos(y). For the region ϕ>0, we take β=1 and the exact solution to be u=cos(x)cos(y)+1. A representation of the solution and the Voronoi mesh are given in Figure 1. The comparison of the Ghost-Fluid Method and the Voronoi Interface Method is shown in Table 1 for the maximum error in the solution, in Table 2 for the average error in the solution, in Table 3 for the maximum error of the gradient, and in Table 4 for the average error of the gradient. For the Voronoi Interface Method, we present the errors on both the Voronoi mesh and the Cartesian mesh. For the Voronoi mesh, the gradients are computed on the faces of the Voronoi cells. For example, the solution only experiences a constant discontinuity in its solution and no discontinuities in its normal or tangential derivatives nor is there a discontinuity in the β coefficient, which is constant in the entire domain. Both methods give second-order accuracy in the solution in the L∞- and L1-norms and the solution’s gradient in the L1-norm. However, the Voronoi method gives first-order accuracy for the gradient of the solution in the L∞-norm while the GFM gives second-order accuracy for such simplified problems. A likely explanation for this difference is that the fluxes in the volume of fluid derivation for the Voronoi Interface Method are not necessarily computed at the center of the faces between two points. We also observe that the interpolation from the Voronoi mesh to the Cartesian mesh does not impact the solution itself but affects its gradient. However, the order of accuracy is conserved. We conclude that for such simple problems the GFM is preferable.Table 1
Comparison of theL∞ error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u-uh∞ Order u-uh∞ Order u-uh∞ Order  0.25 2.627×10-3 — 3.709 × 1 0 - 3 — 4.098×10-1 —  0.125 6.587 × 1 0 - 4 2.00 6.172 × 1 0 - 4 2.59 5.901×10-2 2.80  0.0625 1.648 × 1 0 - 4 2.00 1.565 × 1 0 - 4 1.98 1.788×10-4 8.37  0.03125 4.121 × 1 0 - 5 2.00 4.016 × 1 0 - 5 1.96 4.030×10-5 2.15  0.015625 1.030 × 1 0 - 5 2.00 1.020 × 1 0 - 5 1.98 1.023×10-5 1.98  0.0078125 2.576 × 1 0 - 6 2.00 2.598 × 1 0 - 6 1.97 2.656×10-6 1.95  0.00390625 6.439 × 1 0 - 7 2.00 6.591 × 1 0 - 7 1.98 6.655×10-7 2.00Table 2
Comparison of theL1 error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u - u h 1 Order u - u h 1 Order u - u h 1 Order  0.25 9.010 × 1 0 - 4 — 1.790 × 1 0 - 2 — 4.694 × 1 0 - 2 —  0.125 2.597 × 1 0 - 4 1.79 2.918 × 1 0 - 3 5.04 1.417 × 1 0 - 3 5.04  0.0625 6.940 × 1 0 - 5 1.90 7.365 × 1 0 - 5 4.31 7.130 × 1 0 - 5 4.31  0.03125 1.792 × 1 0 - 5 1.95 1.848 × 1 0 - 5 2.00 1.780 × 1 0 - 5 2.00  0.015625 4.552 × 1 0 - 6 1.98 4.625 × 1 0 - 6 1.98 4.518 × 1 0 - 6 1.98  0.0078125 1.147 × 1 0 - 6 1.99 1.156 × 1 0 - 6 1.98 1.142 × 1 0 - 6 1.98  0.00390625 2.879 × 1 0 - 7 1.99 2.890 × 1 0 - 7 1.99 2.871 × 1 0 - 7 1.99Table 3
Comparison of theL∞ error in the gradient of the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order  0.25 4.829 × 1 0 - 3 — 8.160 × 1 0 - 3 — 1.872 × 1 0 - 0 —  0.125 1.398 × 1 0 - 3 1.79 1.550 × 1 0 - 3 2.40 3.666 × 1 0 - 1 2.35  0.0625 3.776 × 1 0 - 4 1.89 5.192 × 1 0 - 4 1.58 4.180 × 1 0 - 2 3.13  0.03125 9.866 × 1 0 - 5 1.94 4.367 × 1 0 - 4 0.25 2.082 × 1 0 - 2 1.01  0.015625 2.530 × 1 0 - 5 1.96 2.551 × 1 0 - 4 0.78 1.037 × 1 0 - 2 1.00  0.0078125 6.419 × 1 0 - 6 1.98 1.376 × 1 0 - 4 0.89 5.195 × 1 0 - 3 0.99  0.00390625 1.618 × 1 0 - 6 1.99 8.580 × 1 0 - 5 0.68 2.597 × 1 0 - 3 1.00Table 4
Comparison of theL1 error in the gradient of the solution for the Ghost-Fluid method and the Voronoi Interface Method, for example, shown in Section 3.1.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order  0.25 3.636 × 1 0 - 3 — 2.575 × 1 0 - 3 — 3.797 × 1 0 - 1 —  0.125 8.704 × 1 0 - 4 2.06 3.490 × 1 0 - 4 2.88 2.908 × 1 0 - 2 3.71  0.0625 2.155 × 1 0 - 4 2.01 7.553 × 1 0 - 5 2.21 5.617 × 1 0 - 3 2.37  0.03125 5.343 × 1 0 - 5 2.01 1.711 × 1 0 - 5 2.14 1.474 × 1 0 - 3 1.93  0.015625 1.331 × 1 0 - 5 2.01 4.126 × 1 0 - 6 2.05 3.707 × 1 0 - 4 1.99  0.0078125 3.322 × 1 0 - 6 2.00 9.930 × 1 0 - 7 2.05 9.345 × 1 0 - 5 1.99  0.00390625 8.298 × 1 0 - 7 2.00 2.465 × 1 0 - 7 2.01 2.350 × 1 0 - 5 1.99Figure 1
Left: visualization of the exact solution, for example, shown in Section3.1.  Right: example of a Voronoi mesh generated.


## 3.2. Nonconstantβ-Coefficient

Let us consider∇·β∇u=fx,y in two spatial dimensions in Ω=-1,12 with the level set function ϕ=-x2+y2+0.5+0.15cos5arctany/x. We set β=y2lnx+2+4 and the exact solution to u=ex in the region where ϕ≤0. In the region where ϕ>0, we set β=e-y and the exact solution to u=cosxsiny. Figure 2 provides a representation of the solution and of the diffusion coefficient. The comparison of the Ghost-Fluid Method and the Voronoi Interface Method is shown in Table 5 for the maximum error in the solution, in Table 6 for the average error in the solution, in Table 7 for the maximum error of the gradient, and in Table 8 for the average error of the gradient. Again, for the Voronoi Interface Method we present the results on both the Voronoi and the Cartesian meshes. For example, there is a nonconstant discontinuity in the solution, a nonconstant discontinuity in the normal derivative of the solution to the interface, a nonconstant discontinuity in the tangential derivative of the solution to the interface, a discontinuity in the β coefficient across the interface, and a β coefficient that is not constant in each domain as in the previous example. The Ghost-Fluid Method is only first-order accurate in both the L∞- and the L1-norm, as well as in the L1-norm for the gradient of the solution. However, it is not consistent in the L∞-norm for the gradient of the solution. This example is a case where the lack of orthogonality between the cells’ faces and the solution’s fluxes lowers the accuracy of the GFM’s dimension-by-dimension approach. The Voronoi Interface Method provides a solution to that drawback and, therefore, produces a second-order accurate solution in the L∞- and L1-norms, and first-order accurate (resp., second-order accurate) gradients in the L∞- (resp., L1-) norm. Similarly to the previous example, we observe a decrease in accuracy for the gradient computed on the Cartesian mesh after the interpolation step, but the order of accuracy is conserved. The solution itself does not suffer from the interpolation step. For this type of complex problems, the Voronoi Interface Method is able to provide a second-order accurate solution and consistent gradient that the Ghost-Fluid Method cannot produce. VIM is therefore the recommended approach.Table 5
Comparions of theL∞ error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u - u h ∞ Order u - u h ∞ Order u - u h ∞ Order  0.25 6.254 × 1 0 - 1 — 9.802 × 1 0 - 3 — 6.044 × 1 0 - 1 —  0.125 4.063 × 1 0 - 1 0.62 2.560 × 1 0 - 3 1.94 9.178 × 1 0 - 2 2.72  0.0625 2.493 × 1 0 - 1 0.70 7.709 × 1 0 - 4 1.73 7.041 × 1 0 - 4 7.03  0.03125 1.474 × 1 0 - 1 0.76 2.370 × 1 0 - 4 1.70 2.160 × 1 0 - 4 1.70  0.015625 8.488 × 1 0 - 2 0.80 5.560 × 1 0 - 5 2.09 4.929 × 1 0 - 5 2.13  0.0078125 4.802 × 1 0 - 2 0.82 1.413 × 1 0 - 5 1.98 1.293 × 1 0 - 5 1.93  0.00390625 1.481 × 1 0 - 2 0.84 3.618 × 1 0 - 6 1.97 3.377 × 1 0 - 6 1.94Table 6
Comparison of theL1 error in the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx u - u h 1 Order u - u h 1 Order u - u h 1 Order  0.25 6.446 × 1 0 - 2 — 2.291 × 1 0 - 3 — 4.557 × 1 0 - 2 —  0.125 3.393 × 1 0 - 2 0.93 4.022 × 1 0 - 4 2.51 2.041 × 1 0 - 3 4.48  0.0625 1.692 × 1 0 - 2 1.00 1.038 × 1 0 - 4 1.95 9.891 × 1 0 - 5 4.37  0.03125 7.908 × 1 0 - 3 1.10 2.985 × 1 0 - 5 1.80 2.718 × 1 0 - 5 1.86  0.015625 3.982 × 1 0 - 3 0.99 7.243 × 1 0 - 6 2.04 6.729 × 1 0 - 6 2.01  0.0078125 2.155 × 1 0 - 3 0.89 1.790 × 1 0 - 6 2.02 1.717 × 1 0 - 6 1.97  0.00390625 1.035 × 1 0 - 3 1.06 4.509 × 1 0 - 7 1.99 4.401 × 1 0 - 7 1.96Table 7
Comparison of theL∞ error in the gradient of the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order ∇ u - ∇ u h ∞ Order  0.25 1.259 — 1.309 × 1 0 - 2 — 2.658 × 1 0 - 0 —  0.125 1.227 0.04 9.931 × 1 0 - 3 0.55 6.334 × 1 0 - 1 2.07  0.0625 1.139 0.11 2.954 × 1 0 - 3 1.60 8.950 × 1 0 - 2 2.82  0.03125 1.156 -0.02 1.589 × 1 0 - 3 0.89 4.373 × 1 0 - 2 1.03  0.015625 1.188 -0.04 1.038 × 1 0 - 3 0.62 2.155 × 1 0 - 2 1.02  0.0078125 1.205 -0.02 5.177 × 1 0 - 4 1.00 1.070 × 1 0 - 2 1.01  0.00390625 1.214 -0.01 3.194 × 1 0 - 4 0.70 5.329 × 1 0 - 3 1.01Table 8
Comparison of theL1 error in the gradient of the solution for the Ghost-Fluid Method and the Voronoi Interface Method, for example, shown in Section 3.2.
Ghost-Fluid Method Voronoi Interface Method VIM Interpolated dx ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order ∇ u - ∇ u h 1 Order  0.25 2.611 × 1 0 - 1 — 2.126 × 1 0 - 3 — 3.895 × 1 0 - 1 —  0.125 1.438 × 1 0 - 1 0.86 6.395 × 1 0 - 4 1.73 3.400 × 1 0 - 2 3.52  0.0625 7.136 × 1 0 - 2 1.01 1.957 × 1 0 - 4 1.71 5.631 × 1 0 - 3 2.59  0.03125 3.922 × 1 0 - 2 0.86 6.409 × 1 0 - 5 1.61 1.464 × 1 0 - 3 1.94  0.015625 2.019 × 1 0 - 2 0.96 1.593 × 1 0 - 5 2.01 3.713 × 1 0 - 4 1.98  0.0078125 1.033 × 1 0 - 2 0.97 4.132 × 1 0 - 6 1.95 9.359 × 1 0 - 5 1.99  0.00390625 5.222 × 1 0 - 3 0.98 1.099 × 1 0 - 6 1.91 2.358 × 1 0 - 5 1.99Figure 2
Visualization of the exact solutionu (left) and of the diffusion coefficient β (right), for example, shown in Section 3.2.


## 4. Conclusion

This paper has considered the numerical solution of the Poisson equation with jump conditions across an irregular interface. In particular, we have compared the results obtained with the Ghost-Fluid Method and the Voronoi Interface Method. The Ghost-Fluid Method imposes the jump conditions in a dimension-by-dimension framework, leading to a linear system that is symmetric positive definite in which the jump conditions only affect the right-hand side. However, the dimension-by-dimension approach forces a smearing of the tangential quantities in the jump, leading to a loss of accuracy unless both the discontinuity in the solution and the variable coefficientβ are constant. The Voronoi Interface method solves that problem by constructing a Voronoi partition for cells adjacent to the interface. A finite volume discretization over those cells produces discretized fluxes that are orthogonal to the cells’ faces themselves and aligned with the normal direction to the interface. The Ghost-Fluid philosophy can therefore be readily applied, resulting in a linear system that is also symmetric positive definite with only its right-hand side affected by the jump conditions. The resulting solution is second-order accurate (versus first-order accurate in the general case of the Ghost-Fluid Method) and the solution’s gradient is first-order accurate (versus zeroth-order accurate in the case of the Ghost-Fluid Method). The Voronoi Interface Method can therefore be considered superior in general. In the particular case where both the discontinuity across the interface is constant and the variable coefficient β is constant over each subdomains and across the interface, the Ghost-Fluid Method gives a second-order accurate solution and also second-order accurate gradient, giving this approach an advantage over the Voronoi Interface Method. The likely reason is that the discrete fluxes for the Voronoi Interface Method are not located at the center of the faces between two points. Finally, the Voronoi Interface Method provides the solution at the center of the Voronoi cells. An interpolation step is required if the solution on the original Cartesian mesh is needed which can, for example, be carried out with least square interpolations. The order of accuracy is then preserved though the quality of the gradient of the solution is impacted.



---

*Source: 9216703-2018-10-17.xml*