Mathematical and Implementation Details

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

(1)u(x,t)t+q(x,t,u)=S(x,t,u),xΩ,

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 (1).

Interior Discretisation

Let us start by discretising (1) inside Ω. The first step in this discretisation is to compute a triangulation of Ω, decomposing Ω into a collection of disjoint triangles {Tk} so that

Ω=kTk.

This triangulation is typically a constrained Delaunay triangulation, denoted T(Ω), 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 T(Ω) is shown below, where we use a multiply-connected domain to emphasise that these details are not necessarily restricted to simple domains.

Example block output

Control volumes

Key to the FVM are the control volumes, which are used to define volumes Ωi around individual vertices xi that we integrate the PDE over. To define these volumes, take xiΩ, which is a vertex of T(Ω), and take the set of triangles Ti={Tk} that have xi as a vertex. For each of these triangles TkTi, connect its centroid to the midpoints of the triangle's edges. Once this procedure is complete, we obtain a closed polygon around xi that we denote by Ωi and whose interior defines the control volume Ωi. 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.

Example block output

Let us now establish some notation for referring to these control volumes, using the figure below as a reference.

SymbolDescriptionExample
xiA vertex of T(Ω)The blue point below
ΩiThe control volume around xiThe green region below
ΩiThe boundary of ΩiThe blue edges below
ViThe volume of ΩiThe volume of the green region below
EiThe set of edges of ΩiThe blue edges below
σAn edge σEiThe magenta edge below. Note that σEiσ=Ωi
xσThe midpoint of σEiThe blue point below on σ
n^σThe outward unit normal vector to σEiThe black arrow below
TiThe set of triangles that have xi as a vertexThe black triangles surrounding xi below
LσThe length of σEiThe length of the magenta edge below
Example block output

Discretising the PDE

Now that we have our concept of control volumes, we can discretise the PDE (1). We do this by considering each PDE inside each Ωi and integrating. For a given control volume Ωi, we can integrate the PDE to give

(2)ddtΩiudV+ΩiqdV=ΩiSdV.

Using the divergence theorem, the second integral in (2) becomes

(3)Ωiq=Ωiqn^σds=σEiσqn^σds,

where the last equality in (3) follows from integrating over each individual line segment that defines Ωi, which is simply Ei. We then define the control volume averages,

(4)u¯i=1ViΩiudV,S¯i=1ViΩiSdV,

so that our integral formulation (2) becomes

(5)du¯idt+1ViσEiσqn^σds=S¯i.

Note that (5) is still an exact expression.

To proceed, we need to approximate the integrals σqn^σds. To accomplish this, we use a midpoint rule, writing

(6)σqn^σds[q(xσ,t,u(xσ,t))n^σ]Lσ.

Then, replacing the control volume averages with their value at xi, we obtain the following approximation:

(7)duidt+1ViσEi[q(xσ,t,u(xσ,t))n^σ]Lσ=Si,

where ui=u(xi,t) and Si=S(xi,t).

The final step in this part of the approximation is the evaluation of q(xσ,t,u(xσ,t)). To deal with this function, consider some TkTi so that xσ is inside Tk. We use a linear shape function inside Tk to approximate u, writing

(8)u(x,t)=αkx+βky+γk,xTk,

where we have suppressed the dependence of the coefficients (αk,βk,γk) on t. The vertices of Tk are given by vk1,vk2,vk3 with corresponding coordinates xk1,xk2,xk3, respectively. We can then solve for the coefficients in (8) by requiring that u is equal to the values at the vertices of Tk, giving the system of equations

(9)u(xk1,t)=αkxk1+βkyk1+γk,u(xk2,t)=αkxk2+βkyk2+γk,u(xk3,t)=αkxk3+βkyk3+γk,

where xki=(xki,yki)T. Note that the values on the left-hand side of (9) are all known from either the initial condition or the previous time-step. Using Cramer's rule, we define

(10)Sk=1Δk[yk2yk3yk3yk1yk1yk2xk3xk2xk1xk3xk2xk1xk2yk3xk3yk2xk3yk1xk1yk3xk1yk2xk2yk1],

where

(11)Δk=xk1yk2xk2yk1xk1yk3+xk3yk1+xk2yk3xk3yk2,

and thus we obtain

(12)αk=sk,11uk1+sk,12uk2+sk,13uk3,βk=sk,21uk1+sk,22uk2+sk,23uk3,γk=sk,31uk1+sk,32uk2+sk,33uk3,

where uki=u(xki,t) and sk,ij are the elements of Sk. With (8) and (12), we can approximate q(xσ,t,u(xσ,t)) and thus obtain the approximation

(13)duidt+1ViσEi[q(xσ,t,αk(σ)xσ+βk(σ)yσ+γk(σ))n^σ]Lσ=Si,

where k(σ) is the index of the triangle that contains xσ. This linear shape function also allows us to compute gradients like u(xσ,t), since u(xσ,t)=(αk(σ),βk(σ))T.

Boundary Conditions

Let us now discuss how boundary conditions (BCs) are handled. We assume that BCs take on any of the following forms:

(14)q(x,t,u)n^σ=a(x,t,u)xBΩ,(15)du(x,t)dt=a(x,t,u)xBΩ,(16)u(x,t)=a(x,t,u)xBΩ,

where the functions a are user-provided functions. The conditions (14)(16) 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.

Dirichlet boundary conditions

When we have a Dirichlet BC of the form (16), the implementation is simple: Rather than using (13), we instead leave dui/dt=0 and update the value of ui with a(xi,t,ui) at the end of the iteration using a callback; note that the expression ui=a(xi,t,ui) is not an implicit equation for ui, rather it is simply a reassignment of ui to a(xi,t,ui), i.e. uia(xi,t,ui).

Time-dependent Dirichlet boundary conditions

For a time-dependent Dirichlet BC of the form (15), the implementation is again simple: Rather than using (13), simply compute dui/dt=a(xi,t,ui) instead.

Neumann boundary conditions

Neumann boundary conditions (14) are the most complex out of the three. Let us return to our integral formulation (5). Let Ein be the set of edges in Ei that have a Neumann BC associated with them, and Eic=EiEin. Then, also using (13), in the case of (14) we can write

(17)duidt+1ViσEic[q(xσ,t,αk(σ)xσ+βk(σ)yσ+γk(σ))n^σ]Lσ+1ViσEinσaσ(x,t,u)ds=Si,

where aσ is the BC function associated with σ. This integral is then approximated using a midpoint rule as done previously, giving

(18)duidt+1ViσEic[q(xσ,t,u(xσ,t))n^σ]Lσ+1ViσEin[aσ(xσ,t,u(xσ,t))]Lσ=Si,

where u(xσ,t)=αk(σ)xσ+βk(σ)yσ+γk(σ).

Internal Conditions

We also allow for specifying internal conditions, meaning conditions of the form (14)(16) that are applied away from the boundary. We do not currently allow for internal Neumann conditions directly.[1] These conditions are handled in the same way as BCs, except that the user is to provide them per-vertex rather than per-edge.

Putting Everything Together

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 dui/dt. As written, (13) indicates that we loop over each vertex and, within each vertex, loop over each edge of its control volume. On average, |Eic|=12 (since, on average, a point in a Delaunay triangulation has six neighbours), and so computing dui/dt for each i requires O(12n) loop iterates and many repeated computations (since each control volume edge appears in another control volume), where n is the number of points in the mesh. An alternative approach is to instead loop over each triangle in T(Ω) and to then loop over each edge, adding the contributions from each to the associated vertices. This instead requires O(3|T|) loop iterates, where |T| is the number of triangles in T(Ω), and we instead only need to compute the relevant quantities for each control volume edge a single time; note that |T|=O(n) 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.

Example block output

We denote the triangle in blue by T, and refer to the blue, red, and green vertices by vb, vr, and vg, respectively. The relevant edges that contribute to dub/dt, dur/dt, and dug/dt are σbr, σrg, and σgb, as annotated above. In particular, σbr contributes to both dub/dt and dur/dt, σrg contributes to both dur/dt and dug/dt, and σgb contributes to both dug/dt and dub/dt.

Let us focus on ub and ur. The contribution from ebr to dub/dt and dur/dt is given by:

(19)dubdtdubdtQ,durdtdurdt+Q,

where

Q=[q(xbr,t,αxbr+βybr+γ)n^br]Lbr,

and xbr is the midpoint of ebr; n^br is the unit normal vector to the edge σbr=xbrxT, where xT is the centroid of T which should only be computed once for the current T, which should point away from xb but towards xr; Lbr=xTxbr, and α,β,γ are computed from (12) using the vertices of T. Notice that (19) uses a minus sign for dub/dt and a plus sign for dur/dt, because we have brought the sum in (13) 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 xbrxT (i.e., they are not connected to a centroid). There are two ways to handle this:

  1. For each triangle looped over, also check if it is a boundary triangle and then consider its boundary edges.
  2. 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:

Example block output

Consider the edge eij shown in red. There two control volumes that lie on eij, the one for vi and the other for vj. We denote the midpoint of eij by xij=(xi+xj)/2, so that the two control volume edges are xixij and xijxj for vi and vj, respectively. The contributions from the flux over each edge gives

(20)duidtduidtQi,dujdtdujdtQj,

where

(21)Qi=[q(mi,t,αijmix+βijmiy+γij)n^ij]Li,Qj=[q(mj,t,αijmjx+βijmjy+γij)n^ij]Lj,

where mi=(xi+xij)/2=(mix,miy)T, mj=(xij+xj)/2=(mjx,mjy)T, nij is the outward unit normal vector to eij, Li=xijxi, Lj=xjxij, and αij,βij,γij are computed from (12) using the vertices of the triangle that contains eij. If there is a Neumann boundary condition on eij, (21) uses the boundary condition functions for computing the qn^ terms.

Now that we have looped over all triangles and also over all boundary edges, the final values for each dui/dt is given by

duidt1Viduidt+S(xi,t,ui).

Of course, if there is a Dirichlet boundary condition at ui we set dui/dt=0, and if there is a boundary condition on dui/dt we use that boundary condition instead.

Systems of Equations

We also provide support for systems of PDEs that take the form

(22)u1(x,t)t+q1(x,t,u1,,un)=S1(x,t,u1,,un),u2(x,t)t+q2(x,t,u1,,un)=S2(x,t,u1,,un),un(x,t)t+qn(x,t,u1,,un)=Sn(x,t,u1,,un),

where any of the divergences and source terms may also depend on the other variables u1,,un. We can write this as,

(23)u(x,t)t+Q(x,t,u)=S(x,t,u),

where

(24)u(x,t)=[u1(x,t)u2(x,t)un(x,t)]Rn×1,Q(x,t,u)=[q1(x,t,u)q2(x,t,u)qn(x,t,u)]R2×n,S(x,t,u)=[S1(x,t,u)S2(x,t,u)Sn(x,t,u)]Rn×1,

and the divergence Q is defined as

(25)Q(x,t,u)=[q1(x,t,u)q2(x,t,u)qn(x,t,u)]Rn×1.

The method we use for solving these equations is basically the same as what we do for a single PDE. Letting ui=u(xi,t) and Si=S(xi,t), the analogous approximation to (13) is

(26)duidt+1ViσEi[Q(xσ,t,t,αk(σ)xσ+βk(σ)yσ+γk)Tn^σ]Lσ=Si,

where αk(σ), βk(σ), and γk are the vectors of coefficients from (12), e.g.

αk(σ)=sk(σ),11uk1+sk(σ),12uk2+sk(σ),13uk3.

where uki=u(xki,t) and sk(σ),ij are the elements of Sk(σ). The only difference between (26) and (13) 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 FVMProblems. The FVMSystem replaces the initial conditions with a matrix initial condition, where the ith column refers to the ith component of the system.

Steady-State Problems

We provide support for steady-state problems, in which case (1) (and similarly for (22)) becomes

(27)q(x,t,u)=S(x,t,u).

This is solved in exactly the same way, except rootfinding is used and there is no timestepping.

  • 1This is a technical limitation due to how the control volumes are defined. For vertices away from the boundary, the control volume edges do not lie along any of the triangle's edges, which is where we would like to impose Neumann conditions.