Mathematical and Implementation Details

• Mathematical and Implementation Details
• Interior Discretisation
• • Control volumes
  • Discretising the PDE
• Boundary Conditions
• • Dirichlet boundary conditions
  • Time-dependent Dirichlet boundary conditions
  • Neumann boundary conditions
• Internal Conditions
• Putting Everything Together
• Systems of Equations
• Steady-State Problems

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

$$\begin{array}{}\text{(1)}& \frac{\partial u(\mathbf{x},t)}{\partial t}+\mathbf{\nabla }\mathbf{\cdot }\mathbf{q}(\mathbf{x},t,u)=S(\mathbf{x},t,u),\mathbf{x}\in \mathrm{\Omega },\end{array}$$

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 $\text{(1)}$.

Interior Discretisation

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

$$\mathrm{\Omega }=\bigcup _k{T}_k.$$

This triangulation is typically a constrained Delaunay triangulation, denoted $\mathcal{T}(\mathrm{\Omega })$, with appropriate care taken if $\mathrm{\Omega }$ is multiply-connected; these triangulations can be computed using DelaunayTriangulation.jl. An example of such a domain $\mathrm{\Omega }$ and its triangulation $\mathcal{T}(\mathrm{\Omega })$ is shown below, where we use a multiply-connected domain to emphasise that these details are not necessarily restricted to simple domains.

Control volumes

Key to the FVM are the control volumes, which are used to define volumes ${\mathrm{\Omega }}_i$ around individual vertices ${\mathbf{x}}_i$ that we integrate the PDE over. To define these volumes, take ${\mathbf{x}}_i\in \mathrm{\Omega }$, which is a vertex of $\mathcal{T}(\mathrm{\Omega })$, and take the set of triangles ${\mathcal{T}}_i=\{T_k\}$ that have ${\mathbf{x}}_i$ as a vertex. For each of these triangles $T_{k^{\prime }}\in {\mathcal{T}}_i$, connect its centroid to the midpoints of the triangle's edges. Once this procedure is complete, we obtain a closed polygon around ${\mathbf{x}}_i$ that we denote by $\partial {\mathrm{\Omega }}_i$ and whose interior defines the control volume ${\mathrm{\Omega }}_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.

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

Discretising the PDE

Now that we have our concept of control volumes, we can discretise the PDE $\text{(1)}$. We do this by considering each PDE inside each ${\mathrm{\Omega }}_i$ and integrating. For a given control volume ${\mathrm{\Omega }}_i$, we can integrate the PDE to give

$$\begin{array}{}\text{(2)}& \frac{\mathrm{d}}{\mathrm{d}t}{\iint }_{{\mathrm{\Omega }}_i}u\mathrm{d}V+{\iint }_{{\mathrm{\Omega }}_i}\mathbf{\nabla }\mathbf{\cdot }\mathbf{q}\mathrm{d}V={\iint }_{{\mathrm{\Omega }}_i}S\mathrm{d}V.\end{array}$$

Using the divergence theorem, the second integral in $\text{(2)}$ becomes

$$\begin{array}{}\text{(3)}& {\iint }_{{\mathrm{\Omega }}_i}\mathbf{\nabla }\mathbf{\cdot }\mathbf{q}={\oint }_{\partial {\mathrm{\Omega }}_i}\mathbf{q}\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }\mathrm{d}s=\sum _{\sigma \in {\mathcal{E}}_i}{\int }_{\sigma }\mathbf{q}\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }\mathrm{d}s,\end{array}$$

where the last equality in $\text{(3)}$ follows from integrating over each individual line segment that defines $\partial {\mathrm{\Omega }}_i$, which is simply ${\mathcal{E}}_i$. We then define the control volume averages,

$$\begin{array}{}\text{(4)}& {\overline{u}}_i=\frac{1}{V_i}{\iint }_{{\mathrm{\Omega }}_i}u\mathrm{d}V,{\overline{S}}_i=\frac{1}{V_i}{\iint }_{{\mathrm{\Omega }}_i}S\mathrm{d}V,\end{array}$$

so that our integral formulation $\text{(2)}$ becomes

$$\begin{array}{}\text{(5)}& \frac{\mathrm{d}{\overline{u}}_i}{\mathrm{d}t}+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i}{\int }_{\sigma }\mathbf{q}\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }\mathrm{d}s={\overline{S}}_i.\end{array}$$

Note that $\text{(5)}$ is still an exact expression.

To proceed, we need to approximate the integrals ${\int }_{\sigma }\mathbf{q}\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }\mathrm{d}s$. To accomplish this, we use a midpoint rule, writing

$$\begin{array}{}\text{(6)}& {\int }_{\sigma }\mathbf{q}\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }\mathrm{d}s\approx [\mathbf{q}({\mathbf{x}}_{\sigma },t,u({\mathbf{x}}_{\sigma },t))\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }]L_{\sigma }.\end{array}$$

Then, replacing the control volume averages with their value at ${\mathbf{x}}_i$, we obtain the following approximation:

$$\begin{array}{}\text{(7)}& \frac{\mathrm{d}{u}_i}{\mathrm{d}t}+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i}[\mathbf{q}({\mathbf{x}}_{\sigma },t,u({\mathbf{x}}_{\sigma },t))\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }]L_{\sigma }=S_i,\end{array}$$

where $u_i=u({\mathbf{x}}_i,t)$ and $S_i=S({\mathbf{x}}_i,t)$.

The final step in this part of the approximation is the evaluation of $\mathbf{q}({\mathbf{x}}_{\sigma },t,u({\mathbf{x}}_{\sigma },t))$. To deal with this function, consider some $T_k\in {\mathcal{T}}_i$ so that ${\mathbf{x}}_{\sigma }$ is inside $T_k$. We use a linear shape function inside $T_k$ to approximate $u$, writing

$$\begin{array}{}\text{(8)}& u(\mathbf{x},t)={\alpha }_{k}x+{\beta }_{k}y+{\gamma }_k,\mathbf{x}\in T_k,\end{array}$$

where we have suppressed the dependence of the coefficients $({\alpha }_k,{\beta }_k,{\gamma }_k)$ on $t$. The vertices of $T_k$ are given by $v_{k1},v_{k2},v_{k3}$ with corresponding coordinates ${\mathbf{x}}_{k1},{\mathbf{x}}_{k2},{\mathbf{x}}_{k3}$, respectively. We can then solve for the coefficients in $\text{(8)}$ by requiring that $u$ is equal to the values at the vertices of $T_k$, giving the system of equations

$$\begin{array}{}\text{(9)}& \begin{array}{rl}u({\mathbf{x}}_{k1},t)& ={\alpha }_k{x}_{k1}+{\beta }_k{y}_{k1}+{\gamma }_k,\\ u({\mathbf{x}}_{k2},t)& ={\alpha }_k{x}_{k2}+{\beta }_k{y}_{k2}+{\gamma }_k,\\ u({\mathbf{x}}_{k3},t)& ={\alpha }_k{x}_{k3}+{\beta }_k{y}_{k3}+{\gamma }_k,\end{array}\end{array}$$

where ${\mathbf{x}}_{ki}=(x_{ki},y_{ki}{)}^{\mathsf{T}}$. Note that the values on the left-hand side of $\text{(9)}$ are all known from either the initial condition or the previous time-step. Using Cramer's rule, we define

$$\begin{array}{}\text{(10)}& {\mathbf{S}}_k=\frac{1}{{\mathrm{\Delta }}_k}\left[\begin{array}{ccc}{y}_{k2}-y_{k3}& y_{k3}-y_{k1}& y_{k1}-y_{k2}\\ x_{k3}-x_{k2}& x_{k1}-x_{k3}& x_{k2}-x_{k1}\\ x_{k2}{y}_{k3}-x_{k3}{y}_{k2}& x_{k3}{y}_{k1}-x_{k1}{y}_{k3}& x_{k1}{y}_{k2}-x_{k2}{y}_{k1}\end{array}\right],\end{array}$$

where

$$\begin{array}{}\text{(11)}& {\mathrm{\Delta }}_k=x_{k1}{y}_{k2}-x_{k2}{y}_{k1}-x_{k1}{y}_{k3}+x_{k3}{y}_{k1}+x_{k2}{y}_{k3}-x_{k3}{y}_{k2},\end{array}$$

and thus we obtain

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\alpha }_k& =s_{k,11}{u}_{k1}+s_{k,12}{u}_{k2}+s_{k,13}{u}_{k3},\\ {\beta }_k& =s_{k,21}{u}_{k1}+s_{k,22}{u}_{k2}+s_{k,23}{u}_{k3},\\ {\gamma }_k& =s_{k,31}{u}_{k1}+s_{k,32}{u}_{k2}+s_{k,33}{u}_{k3},\end{array}\end{array}$$

where $u_{ki}=u({\mathbf{x}}_{ki},t)$ and $s_{k,ij}$ are the elements of ${\mathbf{S}}_k$. With $\text{(8)}$ and $\text{(12)}$, we can approximate $\mathbf{q}({\mathbf{x}}_{\sigma },t,u({\mathbf{x}}_{\sigma },t))$ and thus obtain the approximation

$$\begin{array}{}\text{(13)}& \frac{\mathrm{d}{u}_i}{\mathrm{d}t}+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i}[\mathbf{q}({\mathbf{x}}_{\sigma },t,{\alpha }_{k(\sigma )}{x}_{\sigma }+{\beta }_{k(\sigma )}{y}_{\sigma }+{\gamma }_{k(\sigma )})\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }]L_{\sigma }=S_i,\end{array}$$

where $k(\sigma )$ is the index of the triangle that contains ${\mathbf{x}}_{\sigma }$. This linear shape function also allows us to compute gradients like $\mathbf{\nabla }u({\mathbf{x}}_{\sigma },t)$, since $\mathbf{\nabla }u({\mathbf{x}}_{\sigma },t)=({\alpha }_{k(\sigma )},{\beta }_{k(\sigma )}{)}^{\mathsf{T}}$.

Boundary Conditions

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

$$\begin{array}{}\text{(14)}& \mathbf{q}(\mathbf{x},t,u)\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }& =a(\mathbf{x},t,u)& \mathbf{x}\in \mathcal{B}\subseteq \partial \mathrm{\Omega },\text{(15)}& \frac{\mathrm{d}u(\mathbf{x},t)}{\mathrm{d}t}& =a(\mathbf{x},t,u)& \mathbf{x}\in \mathcal{B}\subseteq \partial \mathrm{\Omega },\text{(16)}& u(\mathbf{x},t)& =a(\mathbf{x},t,u)& \mathbf{x}\in \mathcal{B}\subseteq \partial \mathrm{\Omega },\end{array}$$

where the functions $a$ are user-provided functions. The conditions $\text{(14)}$–$\text{(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 $\text{(16)}$, the implementation is simple: Rather than using $\text{(13)}$, we instead leave $\mathrm{d}{u}_i/\mathrm{d}t=0$ and update the value of $u_i$ with $a({\mathbf{x}}_i,t,u_i)$ at the end of the iteration using a callback; note that the expression $u_i=a({\mathbf{x}}_i,t,u_i)$ is not an implicit equation for $u_i$, rather it is simply a reassignment of $u_i$ to $a({\mathbf{x}}_i,t,u_i)$, i.e. $u_i\leftarrow a({\mathbf{x}}_i,t,u_i)$.

Time-dependent Dirichlet boundary conditions

For a time-dependent Dirichlet BC of the form $\text{(15)}$, the implementation is again simple: Rather than using $\text{(13)}$, simply compute $\mathrm{d}{u}_i/\mathrm{d}t=a({\mathbf{x}}_i,t,u_i)$ instead.

Neumann boundary conditions

Neumann boundary conditions $\text{(14)}$ are the most complex out of the three. Let us return to our integral formulation $\text{(5)}$. Let ${\mathcal{E}}_i^n$ be the set of edges in ${\mathcal{E}}_i$ that have a Neumann BC associated with them, and ${\mathcal{E}}_i^c={\mathcal{E}}_i\setminus {\mathcal{E}}_i^n$. Then, also using $\text{(13)}$, in the case of $\text{(14)}$ we can write

$$\begin{array}{}\text{(17)}& \frac{\mathrm{d}{u}_i}{\mathrm{d}t}+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i^c}[\mathbf{q}({\mathbf{x}}_{\sigma },t,{\alpha }_{k(\sigma )}{x}_{\sigma }+{\beta }_{k(\sigma )}{y}_{\sigma }+{\gamma }_{k(\sigma )})\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }]L_{\sigma }+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i^n}{\int }_{\sigma }{a}_{\sigma }(\mathbf{x},t,u)\mathrm{d}s=S_i,\end{array}$$

where $a_{\sigma }$ is the BC function associated with $\sigma$. This integral is then approximated using a midpoint rule as done previously, giving

$$\begin{array}{}\text{(18)}& \frac{\mathrm{d}{u}_i}{\mathrm{d}t}+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i^c}[\mathbf{q}({\mathbf{x}}_{\sigma },t,u({\mathbf{x}}_{\sigma },t))\mathbf{\cdot }{\hat{\mathbf{n}}}_{\sigma }]L_{\sigma }+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i^n}[a_{\sigma }({\mathbf{x}}_{\sigma },t,u({\mathbf{x}}_{\sigma },t))]L_{\sigma }=S_i,\end{array}$$

where $u({\mathbf{x}}_{\sigma },t)={\alpha }_{k(\sigma )}{x}_{\sigma }+{\beta }_{k(\sigma )}{y}_{\sigma }+{\gamma }_{k(\sigma )}$.

Internal Conditions

We also allow for specifying internal conditions, meaning conditions of the form $\text{(14)}$–$\text{(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 $\mathrm{d}{u}_i/\mathrm{d}t$. As written, $\text{(13)}$ indicates that we loop over each vertex and, within each vertex, loop over each edge of its control volume. On average, $|{\mathcal{E}}_i^c|=12$ (since, on average, a point in a Delaunay triangulation has six neighbours), and so computing $\mathrm{d}{u}_i/\mathrm{d}t$ for each $i$ requires $\mathcal{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 $\mathcal{T}(\mathrm{\Omega })$ and to then loop over each edge, adding the contributions from each to the associated vertices. This instead requires $\mathcal{O}(3|\mathcal{T}|)$ loop iterates, where $|\mathcal{T}|$ is the number of triangles in $\mathcal{T}(\mathrm{\Omega })$, and we instead only need to compute the relevant quantities for each control volume edge a single time; note that $|\mathcal{T}|=\mathcal{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.

We denote the triangle in blue by $T$, and refer to the blue, red, and green vertices by $v_b$, $v_r$, and $v_g$, respectively. The relevant edges that contribute to $\mathrm{d}{u}_b/\mathrm{d}t$, $\mathrm{d}{u}_r/\mathrm{d}t$, and $\mathrm{d}{u}_g/\mathrm{d}t$ are ${\sigma }_{br}$, ${\sigma }_{rg}$, and ${\sigma }_{gb}$, as annotated above. In particular, ${\sigma }_{br}$ contributes to both $\mathrm{d}{u}_b/\mathrm{d}t$ and $\mathrm{d}{u}_r/\mathrm{d}t$, ${\sigma }_{rg}$ contributes to both $\mathrm{d}{u}_r/\mathrm{d}t$ and $\mathrm{d}{u}_g/\mathrm{d}t$, and ${\sigma }_{gb}$ contributes to both $\mathrm{d}{u}_g/\mathrm{d}t$ and $\mathrm{d}{u}_b/\mathrm{d}t$.

Let us focus on $u_b$ and $u_r$. The contribution from $e_{br}$ to $\mathrm{d}{u}_b/\mathrm{d}t$ and $\mathrm{d}{u}_r/\mathrm{d}t$ is given by:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}\frac{\mathrm{d}{u}_b}{\mathrm{d}t}& \leftarrow \frac{\mathrm{d}{u}_b}{\mathrm{d}t}-Q,\\ \frac{\mathrm{d}{u}_r}{\mathrm{d}t}& \leftarrow \frac{\mathrm{d}{u}_r}{\mathrm{d}t}+Q,\end{array}\end{array}$$

where

$$Q=[\mathbf{q}({\mathbf{x}}_{br},t,\alpha x_{br}+\beta y_{br}+\gamma )\mathbf{\cdot }{\hat{\mathbf{n}}}_{br}]L_{br},$$

and ${\mathbf{x}}_{br}$ is the midpoint of $e_{br}$; ${\hat{\mathbf{n}}}_{br}$ is the unit normal vector to the edge ${\sigma }_{br}=\overrightarrow{{\mathbf{x}}_{br}{\mathbf{x}}_T}$, where ${\mathbf{x}}_T$ is the centroid of $T$ which should only be computed once for the current $T$, which should point away from ${\mathbf{x}}_b$ but towards ${\mathbf{x}}_r$; $L_{br}=\Vert {\mathbf{x}}_T-{\mathbf{x}}_{br}\Vert$, and $\alpha ,\beta ,\gamma$ are computed from $\text{(12)}$ using the vertices of $T$. Notice that $\text{(19)}$ uses a minus sign for $\mathrm{d}{u}_b/\mathrm{d}t$ and a plus sign for $\mathrm{d}{u}_r/\mathrm{d}t$, because we have brought the sum in $\text{(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 $\overrightarrow{{\mathbf{x}}_{br}{\mathbf{x}}_T}$ (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:

Consider the edge $e_{ij}$ shown in red. There two control volumes that lie on $e_{ij}$, the one for $v_i$ and the other for $v_j$. We denote the midpoint of $e_{ij}$ by ${\mathbf{x}}_{ij}=({\mathbf{x}}_i+{\mathbf{x}}_j)/2$, so that the two control volume edges are $\overrightarrow{{\mathbf{x}}_i{\mathbf{x}}_{ij}}$ and $\overrightarrow{{\mathbf{x}}_{ij}{\mathbf{x}}_j}$ for $v_i$ and $v_j$, respectively. The contributions from the flux over each edge gives

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\frac{\mathrm{d}{u}_i}{\mathrm{d}t}& \leftarrow \frac{\mathrm{d}{u}_i}{\mathrm{d}t}-Q_i,\\ \frac{\mathrm{d}{u}_j}{\mathrm{d}t}& \leftarrow \frac{\mathrm{d}{u}_j}{\mathrm{d}t}-Q_j,\end{array}\end{array}$$

where

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{Q}_i& =[\mathbf{q}({\mathbf{m}}_i,t,{\alpha }_{ij}{m}_{ix}+{\beta }_{ij}{m}_{iy}+{\gamma }_{ij})\mathbf{\cdot }{\hat{\mathbf{n}}}_{ij}]L_i,\\ Q_j& =[\mathbf{q}({\mathbf{m}}_j,t,{\alpha }_{ij}{m}_{jx}+{\beta }_{ij}{m}_{jy}+{\gamma }_{ij})\mathbf{\cdot }{\hat{\mathbf{n}}}_{ij}]L_j,\end{array}\end{array}$$

where ${\mathbf{m}}_i=({\mathbf{x}}_i+{\mathbf{x}}_{ij})/2=(m_{ix},m_{iy}{)}^{\mathsf{T}}$, ${\mathbf{m}}_j=({\mathbf{x}}_{ij}+{\mathbf{x}}_j)/2=(m_{jx},m_{jy}{)}^{\mathsf{T}}$, ${\mathbf{n}}_{ij}$ is the outward unit normal vector to $e_{ij}$, $L_i=\Vert {\mathbf{x}}_{ij}-{\mathbf{x}}_i\Vert$, $L_j=\Vert {\mathbf{x}}_j-{\mathbf{x}}_{ij}\Vert$, and ${\alpha }_{ij},{\beta }_{ij},{\gamma }_{ij}$ are computed from $\text{(12)}$ using the vertices of the triangle that contains $e_{ij}$. If there is a Neumann boundary condition on $e_{ij}$, $\text{(21)}$ uses the boundary condition functions for computing the $\mathbf{q}\mathbf{\cdot }\hat{\mathbf{n}}$ terms.

Now that we have looped over all triangles and also over all boundary edges, the final values for each $\mathrm{d}{u}_i/\mathrm{d}t$ is given by

$$\frac{\mathrm{d}{u}_i}{\mathrm{d}t}\leftarrow \frac{1}{V_i}\frac{\mathrm{d}{u}_i}{\mathrm{d}t}+S({\mathbf{x}}_i,t,u_i).$$

Of course, if there is a Dirichlet boundary condition at $u_i$ we set $\mathrm{d}{u}_i/\mathrm{d}t=0$, and if there is a boundary condition on $\mathrm{d}{u}_i/\mathrm{d}t$ we use that boundary condition instead.

Systems of Equations

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

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\frac{\partial u_1(\mathbf{x},t)}{\partial t}+\mathbf{\nabla }\mathbf{\cdot }{\mathbf{q}}_1(\mathbf{x},t,u_1,\dots ,u_n)& =S_1(\mathbf{x},t,u_1,\dots ,u_n),\\ \frac{\partial u_2(\mathbf{x},t)}{\partial t}+\mathbf{\nabla }\mathbf{\cdot }{\mathbf{q}}_2(\mathbf{x},t,u_1,\dots ,u_n)& =S_2(\mathbf{x},t,u_1,\dots ,u_n),\\ & ⋮\\ \frac{\partial u_n(\mathbf{x},t)}{\partial t}+\mathbf{\nabla }\mathbf{\cdot }{\mathbf{q}}_n(\mathbf{x},t,u_1,\dots ,u_n)& =S_n(\mathbf{x},t,u_1,\dots ,u_n),\end{array}\end{array}$$

where any of the divergences and source terms may also depend on the other variables $u_1,\dots ,u_n$. We can write this as,

$$\begin{array}{}\text{(23)}& \frac{\partial \mathbf{u}(\mathbf{x},t)}{\partial t}+\mathbf{\nabla }\mathbf{\cdot }\mathbf{Q}(\mathbf{x},t,\mathbf{u})=\mathbf{S}(\mathbf{x},t,\mathbf{u}),\end{array}$$

where

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\mathbf{u}(\mathbf{x},t)& =\left[\begin{array}{c}{u}_1(\mathbf{x},t)\\ u_2(\mathbf{x},t)\\ ⋮\\ u_n(\mathbf{x},t)\end{array}\right]\in {\mathbb{R}}^{n\times 1},\\ \mathbf{Q}(\mathbf{x},t,\mathbf{u})& =\left[\begin{array}{cccc}{\mathbf{q}}_1(\mathbf{x},t,\mathbf{u})& {\mathbf{q}}_2(\mathbf{x},t,\mathbf{u})& \cdots & {\mathbf{q}}_n(\mathbf{x},t,\mathbf{u})\end{array}\right]\in {\mathbb{R}}^{2\times n},\\ \mathbf{S}(\mathbf{x},t,\mathbf{u})& =\left[\begin{array}{c}{S}_1(\mathbf{x},t,\mathbf{u})\\ S_2(\mathbf{x},t,\mathbf{u})\\ ⋮\\ S_n(\mathbf{x},t,\mathbf{u})\end{array}\right]\in {\mathbb{R}}^{n\times 1},\end{array}\end{array}$$

and the divergence $\mathbf{\nabla }\mathbf{\cdot }\mathbf{Q}$ is defined as

$$\begin{array}{}\text{(25)}& \mathbf{\nabla }\mathbf{\cdot }\mathbf{Q}(\mathbf{x},t,\mathbf{u})=\left[\begin{array}{c}\mathbf{\nabla }\mathbf{\cdot }{\mathbf{q}}_1(\mathbf{x},t,\mathbf{u})\\ \mathbf{\nabla }\mathbf{\cdot }{\mathbf{q}}_2(\mathbf{x},t,\mathbf{u})\\ ⋮\\ \mathbf{\nabla }\mathbf{\cdot }{\mathbf{q}}_n(\mathbf{x},t,\mathbf{u})\end{array}\right]\in {\mathbb{R}}^{n\times 1}.\end{array}$$

The method we use for solving these equations is basically the same as what we do for a single PDE. Letting ${\mathbf{u}}_i=\mathbf{u}({\mathbf{x}}_i,t)$ and ${\mathbf{S}}_i=\mathbf{S}({\mathbf{x}}_i,t)$, the analogous approximation to $\text{(13)}$ is

$$\begin{array}{}\text{(26)}& \frac{\mathrm{d}{u}_i}{\mathrm{d}t}+\frac{1}{V_i}\sum _{\sigma \in {\mathcal{E}}_i}\left[\mathbf{Q}{({\mathbf{x}}_{\sigma },t,t,{\mathit{\alpha }}_{k(\sigma )}{x}_{\sigma }+{\mathit{\beta }}_{k(\sigma )}{y}_{\sigma }+{\mathit{\gamma }}_k)}^{\mathsf{T}}{\hat{\mathbf{n}}}_{\sigma }\right]L_{\sigma }={\mathbf{S}}_i,\end{array}$$

where ${\mathit{\alpha }}_{k(\sigma )}$, ${\mathit{\beta }}_{k(\sigma )}$, and ${\mathit{\gamma }}_k$ are the vectors of coefficients from $\text{(12)}$, e.g.

$${\mathit{\alpha }}_{k(\sigma )}=s_{k(\sigma ),11}{\mathbf{u}}_{k1}+s_{k(\sigma ),12}{\mathbf{u}}_{k2}+s_{k(\sigma ),13}{\mathbf{u}}_{k3}.$$

where ${\mathbf{u}}_{ki}=\mathbf{u}({\mathbf{x}}_{ki},t)$ and $s_{k(\sigma ),ij}$ are the elements of ${\mathbf{S}}_{k(\sigma )}$. The only difference between $\text{(26)}$ and $\text{(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 $\text{(1)}$ (and similarly for $\text{(22)}$) becomes

$$\begin{array}{}\text{(27)}& \mathbf{\nabla }\mathbf{\cdot }\mathbf{q}(\mathbf{x},t,u)=S(\mathbf{x},t,\mathbf{u}).\end{array}$$

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