#### 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 the -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.

#### 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 use* a 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 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], Liu* et 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 , . The governing equation is where and are given. Here, is bounded from below by a positive constant and is in . This equation is supplemented by jump conditions on the irregular interface, where denotes a jump in across . The functions and 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

##### 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 , where refers to the derivative in the normal direction to the interface, is . The discontinuity in , where refers to the derivative in the tangential direction to the interface, is . These equations lead to and . 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 and . The discretization at each grid point is then given by where if and are on the same side of the interface or by 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 and are the cells’ sizes in the - and -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, and are only activated if there is a discontinuity present in the local five point stencil. Here, and and the contribution from the left and right grid points to the current grid point and and are the contributions from the bottom and top grid points.

We give the details for and in the -direction, referring the reader to the original paper [20] for the description for and in the -direction, since they follow the same procedure. If and have opposite signs, define , , where is the jump at the grid node and . is then defined as Similarly, if and have opposite signs, define , , and . is then defined as

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 .

The Voronoi partition associated with those seeds is constructed using a simple geometric algorithm. We refer the reader to the Voro++ 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 and point where any variable of the form would be the quantity at point , is the distance between and , is the length of the face between and , , and is the volume of the Voronoi cell associated with point . 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 in two spatial dimensions in with the level set function . We take when and an exact solution of . For the region , we take and the exact solution to be . 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 - and -norms and the solution’s gradient in the -norm. However, the Voronoi method gives first-order accuracy for the gradient of the solution in the -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.

##### 3.2. Nonconstant -Coefficient

Let us consider in two spatial dimensions in with the level set function . We set and the exact solution to in the region where . In the region where , we set and the exact solution to . 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 - and the -norm, as well as in the -norm for the gradient of the solution. However, it is not consistent in the -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 - and -norms, and first-order accurate (resp., second-order accurate) gradients in the - (resp., -) 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.

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

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Disclosure

The present address of Arthur Guittet is Google LLC, 1600 Amphitheatre Parkway Mountain View, CA 94043.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

#### Acknowledgments

The research of Á. Helgadóttir was supported by the University of Iceland Research Fund 2015 under HI14090070. The researches of A. Guittet and F. Gibou were supported in part by the NSF under DMS-1412695 and DMREF-1534264.