Edit on GitHub

The Voronoi finite volume method

Construction of control volumes

  • Start with a 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, and - for simple geometries - from tensor products of lower dimensional grids.

  • 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

This construction requires (in 2D) that sums of angles opposite to triangle edges are less than $\pi$ and that angles opposite to boudary edges are less than $\frac\pi2$. This property is called the "boundary conforming Delaunay property". It has different equivalent definitions and analogues in 3D.

As a consequence, there is a 1:1 incidence between triangulation nodes and Voronoi cells. Moreover, the angle between the interface between two Voronoi cells and the edge between their corresponding nodes is $\frac\pi2$.

The idea is now to use these Voronoi cells as REVs aka control volumes aka finite volume cells.

The discretization approach

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

\[\begin{aligned} 0&=\int_{\omega_k} \nabla\cdot \vec j \ d\omega =\int_{\partial\omega_k} \vec j\cdot \vec n ds\\ &=\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\\ &\approx \sum_{l\in N_k} \frac{\sigma_{kl}}{h_{kl}}g(u_k, u_l)+ \gamma_k b(u_k) \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|$, $\gamma_k=|\partial\omega_k\cap \partial\Omega|$, 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 associated to the respective nodes divided by the distance between these nodes. The second integral has benn replaced by a boundary condition on the flux $\vec j\cdot \vec n + b(u)=0$. 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.

The flux function $g$ can be derived from usual finite difference formulas.

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

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)+ \gamma_k b(u_k) + |\omega_k| (r(u_k)- f(\vec x_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,b$ has a unique solution.

This approach generalizes to systems of $m$ partial differential equations, which formally can be written in the same way, but assuming that $u$ is an $m$-vector function of $\vec x,t$, and $r,g,b,s$ are $m$-vector funtions of their arguments.

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 boudaries, 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.

Why not this method ?

There are a number of cases where this method needs replaces 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 breaktrough 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 of the triangles
  • physical data: the number $m$ and the functions $s,g,r,b,f$ describing the particular problem.

This structure allows to describe the problem to be solved by data derived from the discretization grid and by the functions describing the physics, giving rise to a software API.

The solution of the nonlinear systems of equations can be performed by Newton's method combined with various direct and iterative linear solvers.

The generic programming capabilities of Julia allow for an implementation of the method which results in an API which consists in the implementation of functions $s,g,r,b,f$ without the need to write code for their derivatives.