How to write a CFD library: Discretising the model equation

In this article, I want to quick-start our discussion on libraries and introduce a model problem that is simple enough for us to code, yet can be exploited to show how we can develop a library that is solving a real-world CFD problem; solving the linear system of equations $Ax=b$ . While we are limiting ourselves here to the 1D heat diffusion equation, it can be applied to any other equation as well (well, we limit ourselves to symmetric matrices, which typically result from pressure Poisson or diffusion equation in general, but they could be extended, let me know if you are interested in that).

In this article, we’ll look at the 1D heat equation and discretise it using the finite volume method. We go through it step by step so that you understand exactly how we arrive at $Ax=b$ . We then discuss the different solution algorithms, i.e. implicit and explicit integration and see that in the absence of a time derivative, we need to solve the system implicitly. We discuss boundary conditions as well and see how the equations change. By the end of this article, you’ll have a firm grasp on our discretised equation and we can use that knowledge in the next articles to write our library. Let’s get started!

In this series

Motivation and background

In the next few articles, I want to walk you through the process of writing a library which has purpose and value for real CFD applications. I want to look at writing a linear algebra solver that provides an approximation for $x$ in the linear system $Ax=b$ . This is a common use case, and every CFD solver under the sun will solve this equation at some point, with the exception of purely explicit solvers.

By the end of this series, you will have a simple but effective library to solve the linear system of equations, which you can apply to any discretised equation, as long as an implicit solution is required. In this article, I want to focus on the maths, i.e. discretising our model equation so that we know what we are solving for and why we need a linear system of equation solver. In subsequent articles, we will look at the library structure and the actual implementation of the various classes.

Discretising the steady-state heat equation

The steady-state heat equation is a good example to study when developing a linear algebra library to solve the discretised equation $Ax=b$ . In the following, I want to walk through the discretisation process so that we know how we are solving the steady-state heat equation. When we discretise an equation, we first need to pick a suitable approximation theory, and we can choose from either the finite difference, finite volume, or finite element method.

Historically speaking, CFD applications were discretised using the finite difference method, which seeks to find approximations of derivatives through a Taylor series expansion. Trying to approximate derivatives in cases of discontinuities, i.e. shock-waves, proved to be a difficult task and so the finite volume method has taken over for mainstream CFD solver development. Within the finite volume method, derivatives are reformulated into integral equations, and instead of approximating derivatives, we are only interested in finding fluxes across cell faces, so shock waves are easier to resolve.

The finite element method can be used for CFD applications as well but has been applied to structural problems primarily. While there is a strong research effort ongoing to use the finite element method in CFD applications, applications of the finite element method in commercial CFD solvers remain, to the best of my knowledge, untested.

The steady-state heat diffusion equation is given below, and we seek to solve this equation in 1D using the finite volume method (which is what pretty much all commercial CFD solver would do):

\Gamma\frac{\partial^2 T}{\partial \mathbf{x}^2}=0

Before we get our hands dirty, though, and derive all the required equations, I want to look at the difference between implicit and explicit time integration techniques.

Understanding the difference between explicit and implicit time integration

Before we continue, I want to look at explicit and implicit numerical schemes first as we need to decide how we want to integrate our equation. While we can only choose implicit in this case, I want to expand on this a bit to show you why. This will also help us understand, in general, how explicit and implicit equations are implemented and solved and we will apply this knowledge shortly to our temperature equation.

To understand explicit and implicit methods, let’s look at an example to clarify this point. Consider the numerical evaluation of the simple advection equation:

\frac{\partial \phi}{\partial t}+a\frac{\partial \phi}{\partial x}=0

Let’s ignore the time derivative for the moment and concentrate on the space derivative, i.e. $a\partial\phi/\partial x$ . We can either calculate the derivative with values from the previous iteration (or initial solution if it is the first time step), or from the next iteration, for which the values of $\phi$ are still unknown. We denote the time level of the previous iteration (or initial solution) with the letter $n$ and the solution at the next time step, i.e. the time level at which we want to obtain a solution in the current iteration with time level $n+1$ .

Explicit schemes make use of the first definition, i.e. all values are known and we evaluate the derivative at time level $n$ . For implicit schemes, we use values at the next time level $n+1$ , even though we don’t know the solution for $\phi$ (yet). Thus, we obtain the following discretised equations, assuming a first-order upwind discretisation in space and time (and assuming $a>0$ :

Explicit time integration:

\frac{\phi^{n+1}_i-\phi^{n}_i}{\Delta t}+a\frac{\phi^{n}_{i}-\phi^{n}_{i-1}}{\Delta x}=0

Implicit time integration:

\frac{\phi^{n+1}_i-\phi^{n}_i}{\Delta t}+a\frac{\phi^{n+1}_{i}-\phi^{n+1}_{i-1}}{\Delta x}=0

Let’s look at this graphically. The following image shows our (1D) mesh on the x-axis, and we have picked an arbitrary point as our current point in the space loop. We denote this point with the letter $i$ and it provides us with a convenient way of referencing neighbouring nodes as $i+1$ and $i-1$ . The y-axis represents time and we have drawn discrete time levels $n$ and $n+1$ . These are the discrete points in time for which we want to obtain a solution.

We can now look at our explicit and implicit equations and draw arrows from values that occur at time level $n$ to $n+1$ . We see the explicit discretisation uses values of $\phi$ at $\phi_{i}^{n}$ , $\phi_{i-1}^{n}$ , $\phi_{i}^{n+1}$ , which are connected in the image above with the orange arrows. The implicit discretisation, on the other hand, has values of $\phi$ at $\phi_{i}^{n}$ , $\phi_{i-1}^{n+1}$ , $\phi_{i}^{n+1}$ , shown by the green arrows.

Let’s also look at this from an equation point of view. It is customary to collect all known quantities (i.e. those at time level $n$ ) on the right-hand side of the equation, and then all unknown quantities on the left-hand side of the equations. If we do that, we end up with the following discretised equations:

Explicit time integration:

\phi_i^{n+1}=\phi_i^n-a\frac{\Delta t}{\Delta x}\left(\phi_i^{n}-\phi_{i-1}^{n}\right)=0

Implicit time integration:

\phi_i^{n+1}\left(\frac{1}{\Delta t}+\frac{a}{\Delta x}\right)+\phi_{i-1}^{n+1}\left(\frac{-a}{\Delta x}\right)=\phi_i^n\left(\frac{1}{\Delta t}\right)

In the case of the explicit treatment, we see that there is only one unknown quantity on the left-hand side so we can solve this equation directly (explicitly). In the figure above, we only have one quantity at time level $n+1$ , respectively. For the implicit treatment, on the other hand, we have now two unknown quantities on the left-hand side and so need to first solve a system of linear equations, and thus we solve the equation indirectly (implicitly).

We commonly express our implicit system in the following form

a_e\phi_{i+1}^{n+1}+a_p\phi_i^{n+1}+a_w\phi_{i-1}^{n+1}=b_i

Comparing terms, we can see that $a_e=0$ , $a_p=(1/\Delta t + a/\Delta x)$ , $a_w=(-a/\Delta x)$ , and $b_i=\phi_i^n (1/\Delta t)$ . But, we could also be writing this in the form $Ax=b$ as

\begin{bmatrix}a_0 & a_1 & & & \dots & 0 \\ a_0 & a_1 & a_2 & & & \vdots \\ & a_1 &a_2 & a_3 & & \\ & & & \ddots & &\\ \vdots & & & a_{k-2} & a_{k-1} & a_{k} \\ 0 & \dots & & & a_{k-1} & a_{k} \end{bmatrix} \begin{bmatrix}\phi_0 \\ \phi_1 \\ \phi_2 \\ \vdots \\ \phi_{k-1} \\ \phi_{k}\end{bmatrix} = \begin{bmatrix}b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_{k-1} \\ b_{k}\end{bmatrix}

In the linear system above, looking at the second row, we have $a_0=a_w$ , $a_1=a_p$ , and $a_2=a_e$ . Solving this linear system of equations in any way, i.e. $x=A^{-1}\b$ will provide a solution for the implicitly discretised equation

\phi_i^{n+1}\left(\frac{1}{\Delta t}+\frac{a}{\Delta x}\right)+\phi_{i-1}^{n+1}\left(\frac{-a}{\Delta x}\right)=\phi_i^n\left(\frac{1}{\Delta t}\right)

For completeness, we could also write our explicit equation in matrix form. Let’s look at the equation again:

\phi_i^{n+1}=\phi_i^n-a\frac{\Delta t}{\Delta x}\left(\phi_i^{n}-\phi_{i-1}^{n}\right)=0

If we want to bring it into the form

a_e\phi_{i+1}^{n+1}+a_p\phi_i^{n+1}+a_w\phi_{i-1}^{n+1}=b_i

then we see that we have $a_p=1$ , while $a_e=0$ and $a_w=0$ as we do not have any neighbouring cell being evaluated at time level $n+1$ . The right-hand side becomes $b_i=\phi_i^n-(a\Delta t/\Delta x)(\phi_i^{n}-\phi_{i-1}^{n})$ and so we would have a trivial vector matrix problem to solve in the form of:

\begin{bmatrix}1 &0 & & & \dots & 0 \\ 0 & 1 & 0 & & & \vdots \\ & 0 &1 & 0 & & \\ & & & \ddots & &\\ \vdots & & & 0 & 1 & 0 \\ 0 & \dots & & & 0 & 1 \end{bmatrix} \begin{bmatrix}\phi_0 \\ \phi_1 \\ \phi_2 \\ \vdots \\ \phi_{k-1} \\ \phi_{k}\end{bmatrix} = \begin{bmatrix}\phi_i^n-(a\Delta t/\Delta x)(\phi_i^{n}-\phi_{boundary}^{n}) \\ \phi_i^n-(a\Delta t/\Delta x)(\phi_i^{n}-\phi_{i-1}^{n}) \\ \phi_i^n-(a\Delta t/\Delta x)(\phi_i^{n}-\phi_{i-1}^{n}) \\ \vdots \\ \phi_i^n-(a\Delta t/\Delta x)(\phi_{k-1}^{n}-\phi_{k-2}^{n}) \\ \phi_i^n-(a\Delta t/\Delta x)(\phi_k^{n}-\phi_{k-1}^{n})\end{bmatrix}

Thus, explicit systems always result in the identity matrix which is trivial to invert. This just means it is easy to compute $A^{-1}[/katex from [katex]A$ and so we can compute $x=A^{-1}\b$ . In the case of the identity matrix, we further have $A=A^{-1}$ and it follows that we can write $x=A\b$ ). But as we have already seen, there is no need to write all of this out as a matrix and we can just directly solve the equation as it is, however, this notation makes sense if you want to support both explicit and implicit methods and want to switch between them seamlessly.

With this discussion out of the way, we are ready to discretise our heat diffusion equation.

This website exists to create a community of like-minded CFD enthusiasts and I'd love to start a discussion with you. If you would like to be part of it, sign up using the link below and you will receive my OpenFOAM quick reference guide, as well as my guide on Tools every CFD developer needs for free.
Join now

Step 1: Finite volume discretisation

The first step involves an integration of the equation over small (finite) volumes, which will later represent the cells of our mesh. The equation is written as

\int_V\Gamma\frac{\partial^2 T}{\partial \mathbf{x}^2}\mathrm{d}V=0

Using the Gauss (or divergence) theorem, we can transform any volume integral into a surface integral:

\int_V\Gamma\frac{\partial^2 T}{\partial \mathbf{x}^2}\mathrm{d}V=\int_V\Gamma\frac{\partial}{\partial \mathbf{x}}\frac{\partial T}{\partial \mathbf{x}}\mathrm{d}V=\int_A\mathbf{n}\cdot\Gamma\frac{\partial T}{\partial \mathbf{x}}\mathrm{d}A=0

The Gauss theorem simply states that instead of checking how much a quantity is changing over a given volume, you can check what goes through its boundaries. Think about a balloon, instead of measuring the volume of the balloon directly (integrate over the shape of the balloon), you could also simply measure the amount of air entering the balloon over the inlet area, and use the density to determine the volume. This is, in a nutshell, what the Gauss theorem allows us to do.

The next step involves turning our integrals into sums. For that, we are going to apply a second-order accurate Gaussian Quadrature rule to our integration. This step results in our area integral being replaced by a summation, where we are now seeking the solution at the cell's faces, specifically at their mid-point. This is illustrated in the following image:

Let's say we are in the centre cell with centroid P. Instead of integrating over the cell volume of P and measuring what the rate of change of temperature is, we apply the Gauss theorem and are only interested in how much temperature enters the cell from the east (e) and west (w) face. We can further state that the normal vector always points outwards of a given cell, so it will point to the right for the east (e) face and the left for the west (w) face. Assuming the coordinate direction is going from the left to the right, then we can write the discretised form of the integral equation above as:

\int_A\mathbf{n}\cdot\Gamma\frac{\partial T}{\partial \mathbf{x}}\mathrm{d}A\approx\sum^{faces} \mathbf{n}\cdot\Gamma\frac{\partial T}{\partial \mathbf{x}}A =\Gamma A\frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{e} - \Gamma A\frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{w}=\Gamma A\left[\frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{e} - \frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{w} \right]=0

We have transformed our second-order derivative into a first-order derivative, so our next task involves finding an approximation for the derivatives at the faces (e) and (w). To numerically approximate a derivative, you may recall, either from high school or undergraduate studies, that we can approximate a derivative with a limit of the form

\frac{\partial T}{\partial x}=\lim_{h\rightarrow 0}\frac{T(x+h)-T(x)}{h}

We can take this approximation and set $h=\mathrm{d}x$ , where $\mathrm{d}x$ is now the size/length of the cell as shown in the Figure above, and obtain

\frac{\partial T}{\partial x}\approx\frac{T(x+\mathrm{d}x)-T(x)}{\mathrm{d}x}

We can now find the derivatives for the east (e) and west (w) faces, respectively, as

\frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{e}\approx\frac{T_E - T_P}{\mathrm{d}x}

and

\frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{w}\approx\frac{T_P - T_W}{\mathrm{d}x}

All of these quantities are known, so we can directly evaluate the derivatives at the face. Inserting both definitions into our discretised equation (and using the more conventional notation of $\mathrm{d}x=\Delta x$ results in the following form:

\Gamma A\left[\frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{e} - \frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{w} \right]\approx\Gamma A\left[\frac{T_E - T_P}{\Delta x} - \frac{T_P - T_W}{\Delta x}\right]=0

Step 2: Rearranging equations into explicit and implicit contributions

To arrive at a linear system of the form $Ax=b$ , we need to rearrange the equations in terms of the temperature values. To begin with, we expand the brackets and fractions as

\Gamma A\frac{T_E}{\Delta x}-\Gamma A\frac{T_P}{\Delta x}-\Gamma A\frac{T_P}{\Delta x}+\Gamma A\frac{T_W}{\Delta x}=0

We can now reverse the fraction as

T_E\frac{\Gamma A}{\Delta x}-T_P\frac{\Gamma A}{\Delta x}-T_P\frac{\Gamma A}{\Delta x}+T_W\frac{\Gamma A}{\Delta x}=0

Grouping terms now together

T_E\left(\frac{\Gamma A}{\Delta x}\right)+T_P\left(\frac{-2\Gamma A}{\Delta x}\right)+T_W\left(\frac{\Gamma A}{\Delta x}\right)=0

I am using here the common finite volume notation of east cell (E), west cell (W) and the cell between those two points (P). But, we could have also written the equation in the following form to bring it in line with our advection equation discretisation seen above as

T_{i+1}\left(\frac{\Gamma A}{\Delta x}\right)+T_i\left(\frac{-2\Gamma A}{\Delta x}\right)+T_{i-1}\left(\frac{\Gamma A}{\Delta x}\right)=0

Now we decide if we want to solve this equation implicitly or explicitly. Well, there is only one choice here. We have to solve the equation implicitly. Why? If we tried to solve this equation explicitly, all temperature values would be known and at time level $n$ , and so we would have no means to advance this equation to the next time level $n+1$ . Thus, we need to solve this equation implicitly. If you do that, then all of a sudden we only have unknown temperatures at time level $n+1$ , and at first it may seem that we have not gained anything. How can we solve this equation now?

We have to look at boundary conditions. Let's set $i=1$ and use an implicit treatment, in this case, we could also write

T_{2}^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)+T_1^{n+1}\left(\frac{-2\Gamma A}{\Delta x}\right)+T_{0}^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)=0

The value of $T_i^{n+1}$ is the value of the boundary, and so in reality we have $T_0^n$ , i.e. since this is a boundary value, we know this value (either through a Dirichlet or Neumann-type boundary condition). We can look at this also from a graphical point of view, in the following image, we show which nodes are known and unknown. If we have 3 vertices or cells, i.e. 2 on the boundaries and one internal, then even though our matrix will still be a 3 by 3 matrix, the first and last row can be directly solved and so we only have a single equation for a single unknown value of $T$ .

If we add a second internal cell, then we have a 4 by 4 matrix. We can't solve rows 1 and 4 directly anymore (as either will now contain only one known and two unknown quantities), but we can write out all 4 equations for $T$ . The same observations hold for three or more cells.

It is worth pointing out at this point that whenever we have a steady state system, i.e. an equation without a time derivative, we have to solve it implicitly. If we want to solve the equation explicitly for whatever reason, we add a pseudo time derivative to the equation, which, as we reach a steady state solution, vanishes and becomes zero (definition of steady state, i.e. no derivative in time).

So, our final equation looks like the following for internal nodes

T_{i+1}^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)+T_i^{n+1}\left(\frac{-2\Gamma A}{\Delta x}\right)+T_{i-1}^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)=0

Step 3: Boundary condition treatment

For boundary nodes, we have to re-derive our equations. Consider the following arrangements of the boundary cells at the west boundary:

We can write the equations in the following form again

\Gamma A\left[\frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{e} - \frac{\partial T}{\partial \mathbf{x}}\Bigm\lvert_{w} \right]\approx\Gamma A\left[\frac{T_E - T_P}{\Delta x} - \frac{T_P - T_b}{\frac{1}{2}\Delta x}\right]=0

Expanding brackets we get (and keep in mind that $1/(1/2)=1/0.5=2/1=2$

\Gamma A\frac{T_E}{\Delta x}-\Gamma A\frac{T_P}{\Delta x}-2\Gamma A\frac{T_P}{\Delta x}+2\Gamma A\frac{T_b}{\Delta x}=0

Expressing the equations again in terms of temperatures, we get

T_E\frac{\Gamma A}{\Delta x}-T_P\frac{\Gamma A}{\Delta x}-T_P\frac{2\Gamma A}{\Delta x}+T_b\frac{2\Gamma A}{\Delta x}=0

Grouping terms results in

T_E\left(\frac{\Gamma A}{\Delta x}\right)+T_P\left(\frac{-3\Gamma A}{\Delta x}\right)+T_b\left(\frac{2\Gamma A}{\Delta x}\right)=0

Since $T_b$ is known, we can write the equation in terms of known and unknown contributions as

T_E^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)+T_P^{n+1}\left(\frac{-3\Gamma A}{\Delta x}\right)=T_b\left(\frac{-2\Gamma A}{\Delta x}\right)

Equation Summary

Let's summarise all equations to have them in one convenient place. For internal nodes, that do not contain any boundary nodes, then we have the following discretised equation:

T_{E}^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)+T_P^{n+1}\left(\frac{-2\Gamma A}{\Delta x}\right)+T_{W}^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)=0

On the boundary on the left side of the domain (west boundary), we have

T_E^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)+T_P^{n+1}\left(\frac{-3\Gamma A}{\Delta x}\right)=T_{W,boundary}\left(\frac{-2\Gamma A}{\Delta x}\right)

On the boundary on the right side of the domain (east boundary), we have

T_P^{n+1}\left(\frac{-3\Gamma A}{\Delta x}\right)+T_W^{n+1}\left(\frac{\Gamma A}{\Delta x}\right)=T_{E,boundary}\left(\frac{-2\Gamma A}{\Delta x}\right)

Summary

We looked at the model equation (heat diffusion) that we want to solve and how to discretise that using the finite volume method. We explored how we can use either explicit or implicit time integration techniques to advance a discretised equation and how the implicit discretisation results in a linear system of equations in the form of $Ax=b$ . We derived the discretised heat diffusion equation step by step and learned how to treat internal and boundary cells. With that knowledge in hand, we are ready to look at the library that will help us solve the linear system of equations $Ax=b$

Tom-Robin Teschner is a senior lecturer in computational fluid dynamics and course director for the MSc in computational fluid dynamics and the MSc in aerospace computational engineering at Cranfield University.