Here we outline the mathematical and implementation details involved in implementing the finite volume method (FVM). We assume that our partial differential equation (PDE) is given by
together with some boundary conditions or internal conditions that we discuss later. We also discuss steady-state problems and systems of PDEs of the form .
Let us start by discretising inside . The first step in this discretisation is to compute a triangulation of , decomposing into a collection of disjoint triangles so that
This triangulation is typically a constrained Delaunay triangulation, denoted , with appropriate care taken if is multiply-connected; these triangulations can be computed using DelaunayTriangulation.jl. An example of such a domain and its triangulation is shown below, where we use a multiply-connected domain to emphasise that these details are not necessarily restricted to simple domains.

Key to the FVM are the control volumes, which are used to define volumes around individual vertices that we integrate the PDE over. To define these volumes, take , which is a vertex of , and take the set of triangles that have as a vertex. For each of these triangles , connect its centroid to the midpoints of the triangle's edges. Once this procedure is complete, we obtain a closed polygon around that we denote by and whose interior defines the control volume . We show the result of this procedure, applied to the above domain, below, where we show the centroids in red and the control volume boundaries in blue.

Let us now establish some notation for referring to these control volumes, using the figure below as a reference.
Symbol | Description | Example |
---|
| A vertex of | The blue point below |
| The control volume around | The green region below |
| The boundary of | The blue edges below |
| The volume of | The volume of the green region below |
| The set of edges of | The blue edges below |
| An edge | The magenta edge below. Note that |
| The midpoint of | The blue point below on |
| The outward unit normal vector to | The black arrow below |
| The set of triangles that have as a vertex | The black triangles surrounding below |
| The length of | The length of the magenta edge below |

Now that we have our concept of control volumes, we can discretise the PDE . We do this by considering each PDE inside each and integrating. For a given control volume , we can integrate the PDE to give
Using the divergence theorem, the second integral in becomes
where the last equality in follows from integrating over each individual line segment that defines , which is simply . We then define the control volume averages,
so that our integral formulation becomes
Note that is still an exact expression.
To proceed, we need to approximate the integrals . To accomplish this, we use a midpoint rule, writing
Then, replacing the control volume averages with their value at , we obtain the following approximation:
where and .
The final step in this part of the approximation is the evaluation of . To deal with this function, consider some so that is inside . We use a linear shape function inside to approximate , writing
where we have suppressed the dependence of the coefficients on . The vertices of are given by with corresponding coordinates , respectively. We can then solve for the coefficients in by requiring that is equal to the values at the vertices of , giving the system of equations
where . Note that the values on the left-hand side of are all known from either the initial condition or the previous time-step. Using Cramer's rule, we define
where
and thus we obtain
where and are the elements of . With and , we can approximate and thus obtain the approximation
where is the index of the triangle that contains . This linear shape function also allows us to compute gradients like , since .
Let us now discuss how boundary conditions (BCs) are handled. We assume that BCs take on any of the following forms:
where the functions are user-provided functions. The conditions – are called Neumann, time-dependent Dirichlet, and Dirichlet, respectively. We discuss how we handle incompatible BCs below, and then how each of these three types are implemented.
When we have a Dirichlet BC of the form , the implementation is simple: Rather than using , we instead leave and update the value of with at the end of the iteration using a callback; note that the expression is not an implicit equation for , rather it is simply a reassignment of to , i.e. .
For a time-dependent Dirichlet BC of the form , the implementation is again simple: Rather than using , simply compute instead.
Neumann boundary conditions are the most complex out of the three. Let us return to our integral formulation . Let be the set of edges in that have a Neumann BC associated with them, and . Then, also using , in the case of we can write
where is the BC function associated with . This integral is then approximated using a midpoint rule as done previously, giving
where .
We also allow for specifying internal conditions, meaning conditions of the form – that are applied away from the boundary. We do not currently allow for internal Neumann conditions directly. These conditions are handled in the same way as BCs, except that the user is to provide them per-vertex rather than per-edge.
We have now specified how we discretise the PDE itself, and how we handle both boundary and internal conditions. The remaining task is to actually discuss how we compute . As written, indicates that we loop over each vertex and, within each vertex, loop over each edge of its control volume. On average, (since, on average, a point in a Delaunay triangulation has six neighbours), and so computing for each requires loop iterates and many repeated computations (since each control volume edge appears in another control volume), where is the number of points in the mesh. An alternative approach is to instead loop over each triangle in and to then loop over each edge, adding the contributions from each to the associated vertices. This instead requires loop iterates, where is the number of triangles in , and we instead only need to compute the relevant quantities for each control volume edge a single time; note that by Euler's formula. This is the approach we take in our implementation.
Let us think about how we can instead loop over each triangle. Consider an interior control volume, shown below.

We denote the triangle in blue by , and refer to the blue, red, and green vertices by , , and , respectively. The relevant edges that contribute to , , and are , , and , as annotated above. In particular, contributes to both and , contributes to both and , and contributes to both and .
Let us focus on and . The contribution from to and is given by:
where
and is the midpoint of ; is the unit normal vector to the edge , where is the centroid of which should only be computed once for the current , which should point away from but towards ; , and are computed from using the vertices of . Notice that uses a minus sign for and a plus sign for , because we have brought the sum in to the right-hand side of the equation.
When we apply the procedure above to each triangle, we will have computed the contribution from each edge to each vertex - almost. The only issue is with boundary triangles, where the edges that lie on the boundary will not be iterated over as they not of the form (i.e., they are not connected to a centroid). There are two ways to handle this:
- For each triangle looped over, also check if it is a boundary triangle and then consider its boundary edges.
- After looping over all triangles, loop over all boundary edges to pick up the missing contributions.
The second approach is preferable, as we don't need to worry about needless checks for boundary triangles, the number of boundary edges it has, etc.
To understand how to pick up contributions from a single edge, consider the figure below which shows some control volumes in the corner of a domain:

Consider the edge shown in red. There two control volumes that lie on , the one for and the other for . We denote the midpoint of by , so that the two control volume edges are and for and , respectively. The contributions from the flux over each edge gives
where
where , , is the outward unit normal vector to , , , and are computed from using the vertices of the triangle that contains . If there is a Neumann boundary condition on , uses the boundary condition functions for computing the terms.
Now that we have looped over all triangles and also over all boundary edges, the final values for each is given by
Of course, if there is a Dirichlet boundary condition at we set , and if there is a boundary condition on we use that boundary condition instead.
We also provide support for systems of PDEs that take the form
where any of the divergences and source terms may also depend on the other variables . We can write this as,
where
and the divergence is defined as
The method we use for solving these equations is basically the same as what we do for a single PDE. Letting and , the analogous approximation to is
where , , and are the vectors of coefficients from , e.g.
where and are the elements of . The only difference between and is that we now have a vector of PDEs rather than a single PDE. The procedure for computing the contribution from each edge is the same as before, except that we now have to loop over each PDE in the system. In the code, we use a FVMSystem
type to represent a system of PDEs, constructed by providing a vector of FVMProblem
s. The FVMSystem
replaces the initial conditions with a matrix initial condition, where the i
th column refers to the i
th component of the system.
We provide support for steady-state problems, in which case (and similarly for ) becomes
This is solved in exactly the same way, except rootfinding is used and there is no timestepping.