The Voronoi finite volume method

Construction of control volumes

  • Start with a boundary conforming Delaunay triangulation of a polygonal domain (intervals in 1D, triangles in 2D, tetrahedra in 3D). Such a triangulation can be generated by e.g. by the mesh generators triangle and TetGen. These are available in Julia via Triangulate.jl and TetGen.jl. For simple geometries – tensor products of lower dimensional grids – such triangulation can be created more easily. The package ExtendableGrids.jl manages the grid data structure which is used in this package. SimplexGridFactory.jl interfaces this grid structure with Triangulate.jl and TetGen.jl and provides an API for incrementally setting up geometry descriptions.

  • Join triangle circumcenters by lines $\rightarrow$ create Voronoi cells which can serve as control volumes, akin to representative elementary volumes (REV) used to derive conservation laws.

  • Black + green: triangle nodes
  • Gray: triangle edges
  • Blue: triangle circumcenters
  • Red: Boundaries of Voronoi cells

In order to make this construction possible, the triangulation must have the boundary conforming Delaunay property:

  • The interior of any triangle circumcircle does not contain any other node of the triangulation
  • All circumcircle centers lay within the domain

In 2D, an equivalent condition is:

  • The sum of triangle angles opposite to a given interior edge is less than $\pi$
  • Triangle angles opposite to boundary edges are less than $\frac\pi2$.

As a consequence, there is a 1:1 incidence between triangulation nodes and Voronoi cells. Moreover, the angle between the interface between two neighboring Voronoi cells and the edge between their corresponding nodes is $\frac\pi2$. Therefore the edge direction is aligned with the normal direction with respect to the boundary of the Voronoi cell. This makes it easy to use these Voronoi cells as REVs aka control volumes aka finite volume cells and to derive a space discretization for a conservation law based on very same balance rules used to derive this conservation law.

The discretization approach

Given a continuity equation $\nabla\cdot \vec j=f$ in a domain $\Omega$, integrate it over a control volume $\omega_k$ with associated node $\vec x_k$ and apply Gauss theorem:

\[\begin{aligned} 0&=\int_{\omega_k} (\nabla\cdot \vec j -f )\ d\omega =\int_{\partial\omega_k} \vec j\cdot \vec n ds - \int_{\omega_k} f d\omega\\ &=\sum_{l\in N_k} \int_{\omega_k\cap \omega_l} \vec j\cdot \vec n ds + \int_{\partial\omega_k\cap \partial\Omega} \vec j\cdot \vec n ds - \int_{\omega_k} f d\omega \\ &\approx \sum_{l\in N_k} \frac{\sigma_{kl}}{h_{kl}}g(u_k, u_l) - |\omega_k| f_k + \text{boundary terms} \end{aligned}\]

Here, $N_k$ is the set of neighbor control volumes, $\sigma_{kl}=|\omega_k\cap \omega_l|$, $h_{kl}=|\vec x_k -\vec x_l|$, where $|\cdot|$ denotes the measure (length resp. area) of a geometrical entity. In the approximation step, we replaced the normal flux integral over the interface between two control volumes by the measure of this interface multiplied by a function depending on the unknowns $u_k, u_l$ associated to the respective nodes divided by the distance between these nodes. The flux function $g$ can be derived from usual finite difference formulas discretizing a particular flux law.

Flux laws

For instance, for the diffusion flux $\vec j=-D\vec\nabla u$, we use $g(u_k, u_l)=D(u_k -u_l)$.

For a convective diffusion flux $\vec j = -D\vec \nabla u + u \vec v$, one can chose the upwind flux

\[\begin{aligned} g(u_k, u_l)=D(u_k -u_l) + v_{kl}\begin{cases} u_k,& v_{kl}>0\\ u_l,& v_{kl}\leq 0, \end{cases} \end{aligned}\]

where $v_{kl}=\frac{h_{kl}}{\sigma_{kl}}\int_{\omega_k\cap \omega_l} \vec v \cdot \vec n_{kl} \ ds$ Fluxes also can depend nonlinearily on $u$.

Boundary conditions

To implement a Robin boundary condition on $\Gamma=\partial\Omega$

\[- \vec j \cdot \vec n + a u = b,\]

we note that by the very construction, the discretization nodes associated to control volumes adjacent to the domain boundary are located at the domain boundary, thus we can assume that the boundary condition is valid in the corresponding collocation node $u_k$. We assume that $\partial\omega_k\cap \partial_\Omega= \cup_{m\in\mathcal M_k} \gamma_{km}$ is the union of a finite number of line (plane) segments. For interior nodes, we set $\mathcal M_k = \emptyset$ . Thus, for the boundary terms in the above equation, we have

\[\begin{aligned} \text{boundary terms}&=\sum_{m\in\mathcal M_k} \int_{\gamma_{km}} \vec j \cdot \vec n d s &\approx \sum_{m\in\mathcal M_k} |\gamma_{km}| \vec j \cdot \vec n\\ &\approx\sum_{m\in\mathcal M_k} |\gamma_{km}| (au_k -b), \end{aligned}\]

We observe that for $\varepsilon\to 0$, the Robin boundary condition

\[- \vec j \cdot \vec n + \frac{1}{\varepsilon}u = \frac{1}{\varepsilon}g\]

tends to the Dirichlet bundary condition

\[ u=g.\]

Therefore, a Dirichlet boundary condition can be approximated by choosing a small value of $\varepsilon$ and implying the aforementioned Robin boundary conditions. This approach called penalty method is chosen for the implementation of Dirichlet boundary conditions in this package.

Time dependent problems, reaction terms

This approach easily generalizes to time dependent nonlinear transport-reaction problems with storage terms $s(u)$, reaction terms $r(u)$ and source terms $f$:

\[\partial_t s(u) + \nabla \cdot \vec j + r(u) -f =0\]

Semidiscretization in time (for implicit Euler) leads to

\[\frac{s(u)-s(u^\flat)}{\tau} + \nabla \cdot \vec j + r(u) -f =0\]

where $\tau$ is the time step size and $u^\flat$ is the solution from the old timestep. The approximation approach then for each control volume gives

\[|\omega_k|\frac{s(u_k)-s(u_k^\flat)}{\tau} + \sum_{l\in N_k} \frac{\sigma_{kl}}{h_{kl}}g(u_k, u_l)+ \sum_{m\in\mathcal M_k} |\gamma_{km}| (au_k -b) + |\omega_k| (r(u_k)- f_k)=0\]

If $n$ is the number of discretization nodes, we get a system of $n$ equations with $n$ unknowns which under proper conditions on $r,g,s$ has a unique solution.

The implicit Euler method is the default solver for time dependent problems. Alternatively, ODE and DAE solvers from DifferentialEquations.jl can be used with an upcoming glue package.

Generalizations to systems

This approach generalizes to systems of partial differential equations, which formally can be written in the same way, but assuming that $u$ is a vector function of $\vec x,t$, and $r,g,s$ are vector functions of their arguments. The package allows to handle different sets of species in different subdomains of $\Omega$.

Boundary reactions, boundary species

In addition to $r,g,s$, the package allows to specify additional boundary species, boundary reaction, boundary flux and boundary storage terms.

Why this method ?

Independent of space dimension, the method (with properly chosen flux functions) is able to preserve a number of physical quantities if they are present on the continuous level:

  • local and global mass conservation
  • positivity of solutions
  • maximum principle: in the absence of source and reaction terms, local extrema of the stationary solution are located at the domain boundaries, never in the interior. For transient problems, local extrema in the interior can only come from the initial value.
  • Consistency to thermodynamics: entropy production etc.

Many of these properties are hard to prove for finite element methods, in particular for the convection-diffusion case.

Where is this method not appropriate ?

There are a number of cases where this method needs to be replaced by something else or at least to be applied with great care:

  • Anisotropic diffusion only works with proper mesh alignment
  • Strongly varying capacity (in the function $s$) at domain interfaces lead to inexact breakthrough curves.
  • Sharp moving convection fronts are smeared out too strongly

History and literature

The following list is work in progress and incomplete, but it references some sources behind the ideas in this package.

  • Macneal, R. H. (1953). An asymmetrical finite difference network. Quarterly of Applied Mathematics, 11(3), 295-310. (pdf via JSTOR). Perhaps this is the earliest mentioning of the method. Note that it was used on an electrical analog computer.
  • Gärtner, K., & Kamenski, L. (2019). Why do we need Voronoi cells and Delaunay meshes? arXiv preprint arXiv:1905.01738. A recent overview on the merits of the method. One of the authors belongs to the pioneers of its application in 3D.
  • Fuhrmann, J., & Langmach, H. (2001). Stability and existence of solutions of time-implicit finite volume schemes for viscous nonlinear conservation laws. Applied Numerical Mathematics, 37(1-2), 201-230. A discussion of the method applied to rather general nonlinear scalar problems.
  • Si, H., Gärtner, K., & Fuhrmann, J. (2010). Boundary conforming Delaunay mesh generation. Computational Mathematics and Mathematical Physics, 50(1), 38-53. Definition of the boundary conforming Delaunay property.
  • Eymard, R., Fuhrmann, J., & Gärtner, K. (2006). A finite volume scheme for nonlinear parabolic equations derived from one-dimensional local Dirichlet problems. Numerische Mathematik, 102(3), 463-495. General concept of the derivation of upwind fluxes for nonlinear problems.
  • Farrell, P., Rotundo, N., Doan, D. H., Kantner, M., Fuhrmann, J., & Koprucki, T. (2017). Drift-diffusion models. In Handbook of Optoelectronic Device Modeling and Simulation (pp. 733-772). CRC Press. Overview and introduction to the method applied to semiconductor device simulation. This problem class profits most from the desirable properties of the method.

Software API and implementation

The entities describing the discrete system can be subdivided into two categories:

  • Geometrical data: $|\omega_k|, \gamma_k, \sigma_{kl}, h_{kl}$ together with the connectivity information simplex grid. These data are calculated from the discretization grid.
  • Physics data: the number of species and the functions $s,g,r,f$ etc. describing the particular problem.

The solution of the nonlinear systems of equations is performed by Newton's method combined with various direct and iterative linear solvers. The Jacobi matrices used in Newton's method are assembled from the constitutive functions with the help of forward mode automatic differentiation implemented in ForwardDiff.jl.