In last week’s article, we developed the discretised equation for the 1D heat diffusion equation using the finite volume method. We looked at why we need an implicit time integration technique and why this requires a linear system of equations solver. In this week’s article, we’ll explore how to set up a basic project structure that we’ll use to code our library. We’ll also write the code for the 1D heat equation solver, even though it won’t compile yet without the library, but this provides us a test for the library and allows us to design the classes’ interfaces. Let’s get started!
Having a common library structure will help others to work with your library more easily. By the end of this article, you’ll understand the basic library structure required to develop our (and, in fact, any other) library and we will build upon this structure in the next articles.
Download Resources
All developed code and resources in this article are available for download. If you are encountering issues running any of the scripts, please refer to the instructions for running scripts downloaded from this website.
In this series
- Part 1: Understanding static, dynamic, and header-only C++ libraries
- Part 2: How to write a CFD library: Discretising the model equation
- Part 3: How to write a CFD library: Basic library structure
- Part 4: How to write a CFD library: The Conjugate Gradient class
- Part 5: How to write a CFD library: The vector class
- Part 6: How to write a CFD library: The sparse matrix class
- Part 7: How to write a CFD library: Compiling and testing
- Part 8: How to integrate external CFD libraries in your code
- Part 9: How to handle C++ libraries with a package manager
In this article
Basic library structure
Let’s start by exploring the basic library structure that I’ll be using. We assume that we have some root folder in which all of our files and directories reside. Within that folder, we have two folders and two files, one folder to collect all build and executable files, one folder for the library code itself, and then one file for the heat equation solver, as well as a script to compile and execute our solver. This structure is shown below:
root
├── build
├── linearAlgebraLib
│ ├── linearAlgebraLib.hpp
│ └── src
│ ├── vector.(hpp|cpp)
│ ├── sparseMatrixCSR.(hpp|cpp)
│ ├── linearAlgebraSolverBase.(hpp|cpp)
│ └── conjugateGradient.(hpp|cpp)
├── main.cpp
└── buildAndRun.(sh|bat)The linearAlgebraLib folder contains the actual library code that we will develop. We have a header include file named linearAlgebraLib.hpp, which is the file we will later include in our code to make use of the library. Within the src folder, we have all required library classes.
The vector.(hpp|cpp) files contain a class Vector that provides some overloaded operators, predominantly, so that we can write easy vector-vector and matrix-vector multiplications. The sparseMatrixCSR.(hpp|cpp) files, on the other hand, define the SparseMatrixCSR class which we’ll be using to store our coefficients. It is a sparse matrix and uses the compressed row storage (CRS) data structure, which allows us to only store non-zero elements in the matrix to reduce memory requirements.
The linearAlgebraSolverBase.(hpp|cpp) files define an interface for a matrix solver to implement, and we define here just one such matrix solver based on the conjugate gradient method (within the conjugateGradient.(hpp|cpp) files). This is the file where we will solve the actual implementation of x=A\b.
Finally, within the main.cpp file, we will define a test that we will use to make sure the conjugate gradient algorithm has worked, which is the implementation of the heat equation that we have discretised in the previous article. Furthermore, the buildAndRun.(sh|bat) files provide scripts for either UNIX (Linux/macOS) or Windows to compile and run the code. Ideally, these should be replaced by a build script but we have not looked at them yet, so this is something we’ll return to in the future.
The “linearAlgebraLib.hpp” header include file
Each library comes with at least one header-include file. The job of this file is to provide an entry point into your library. Other programmers can use your library by including this header file in their code. If your library gets particularly big, then you may have more than one header file, but in our case, we are looking at a relatively small library, so we’ll stick with a single header file. The content of the linearAlgebraLib.hpp file is given below:
#include "src/sparseMatrixCSR.hpp"
#include "src/vector.hpp"
#include "src/linearAlgebraSolverBase.hpp"
#include "src/conjugateGradient.hpp"All that this file does is include the respective header files for each class that is located in the src folder, i.e. where we keep all of our source files (*.hpp and *.cpp files). Later, when we want to use the library, all we do is include this header file in our code and the remaining headers will be loaded. We need to compile all *.cpp files along with our project (for now), but later will see how we can separate the compilation process for our library from the compilation process of our main.cpp file (the way an actual library would be compiled and used).
Compiling and executing the library code
In our project structure above, we saw that we have a buildAndRun.sh and buildAndRun.bat file, which can be used to compile our library code. We use the buildAndRun.sh file on UNIX (Linux, macOS) and the buildAndRun.bat file for Windows. Let’s look into the these files in turn.
Unix (Linux and macOS)
The content of the buildAndRun.sh file looks as follows:
#!/bin/bash
rm -rf build
mkdir -p build
# compile
g++ -c -g -O0 -Wall -Wextra -I. linearAlgebraLib/src/sparseMatrixCSR.cpp -o build/sparseMatrixCSR.o
g++ -c -g -O0 -Wall -Wextra -I. linearAlgebraLib/src/vector.cpp -o build/vector.o
g++ -c -g -O0 -Wall -Wextra -I. linearAlgebraLib/src/linearAlgebraSolverBase.cpp -o build/linearAlgebraSolver.o
g++ -c -g -O0 -Wall -Wextra -I. linearAlgebraLib/src/conjugateGradient.cpp -o build/conjugateGradient.o
g++ -c -g -O0 -Wall -Wextra -I. main.cpp -o build/main.o
# link
g++ -o build/main.out build/sparseMatrixCSR.o build/vector.o build/linearAlgebraSolver.o build/conjugateGradient.o build/main.o
# run
./build/main.outWe create a build directory on line 1 and then compile all library code along with the main.cpp file on lines 7-11. On line 14, we link all code together and finally execute the compiled executable on line 17. Note, if you are using this file and don’t have g++ installed, install gcc first through your package manager (use sudo apt install gcc on Ubuntu and brew install gcc on macOS. For Ubuntu, I would recommend installing the build-essential package, which does contain gcc/g++ and a range of other useful developer tools and libraries, in this case, use sudo apt install build-essential).
Here, the compiler flags mean the following:
- –
c: Compile the current file only into an object, but don’t link it just yet. -g: Include debugging information, useful for when we want to step through our code to find bugs.-O0: Lowest level of code optimisation, i.e. level0(= no optimisation).-Wall: (W)arn all = Turn all (useful) compiler warnings on. It does not contain all compiler warnings that exist, but it will include the most commonly used one- –
Wextra: Turn on some extra warnings and compiler checks. Together with-Wall, the compiler will do a good job of warning us of any potential syntactic issues in our code. -I.: Include the root folder (indicated by.) as a search path for header files. The root folder is defined as the folder from where we call the compiler, which in this case is the root folder as this one contains ourbuildAndRun.(sh|bat)scripts. We do this so that we can import any header files relative to the root directory. If we do not include the-I.flag, then we need to specify all paths as absolute paths, which will break your code if you want to execute your code on a different machine.-o: Directory where to put the compiled code, i.e. in this case we want to move everything into thebuildfolder that we generated on line 1.
Windows
On Windows, we make use of the buildAndRun.bat file and its content looks as follows:
@echo off
del /q build
mkdir build
REM compile
cl /nologo /c /EHsc /Zi /Od /I. /Fobuild\sparseMatrixCSR.obj linearAlgebraLib\src\sparseMatrixCSR.cpp
cl /nologo /c /EHsc /Zi /Od /I. /Fobuild\vector.obj linearAlgebraLib\src\vector.cpp
cl /nologo /c /EHsc /Zi /Od /I. /Fobuild\linearAlgebraSolver.obj linearAlgebraLib\src\linearAlgebraSolverBase.cpp
cl /nologo /c /EHsc /Zi /Od /I. /Fobuild\conjugateGradient.obj linearAlgebraLib\src\conjugateGradient.cpp
cl /nologo /c /EHsc /Zi /Od /I. /Fobuild\main.obj main.cpp
REM link
cl /Fe:build\main.exe build\sparseMatrixCSR.obj build\vector.obj build\linearAlgebraSolver.obj build\conjugateGradient.obj build\main.obj
REM run
build\main.exeTo use that file, make sure you have Visual Studio installed (either the community, professional or enterprise edition) and that you have a suitable C++ development environment installed. Search on your Windows search bar for Visual Studio Installer, select the version of Visual Studio that you want to use (you probably only see one version) and select the Modify button. Next, select a package that contains MSVC and a C++ compiler as a minimum. In my case, I use the Desktop Development with C++ package, which contains the most useful C++ tools and is loosely related to the build-essential package on Ubuntu.
When you start a new PowerShell session, you want to start the Visual Studio developer version of the PowerShell (in my case, it is called Developer PowerShell for VS 2022), which will set up your terminal so that you can use the compiler straight away. Alternatively, you can start a normal PowerShell session and then source the files manually. Type the following command and ensure the file exists at the specified path, you may have to change the version and/or the flavour of your Visual Studio installation (Community, Professional, Enterprise)
cmd.exe /k "C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvarsall.bat" amd64
The amd64 at the end just specifies my PC’s architecture and in this day and age you should have a 64-bit processor as well (unless you insist on developing on a Surface Pro, in which case you may still have a 32-bit processor, and the command should be replaced with x86 (yes I know, makes perfect sense …)).
The actual commands are essentially the same on Windows as on UNIX, the syntax just changes, i.e. we have /c for compiling but not linking, /Zi for debugging information, /Od for no optimisation, /I. to include the root folder, and /Fo to specify the output file for the compiled code. On Windows, you get a compilation error if there is a space between /Fo and the location where you want to place that file, so make sure there is no space (see above as well for reference).
The additional /nologo flag is probably best described as a spam filter, i.e. you are not getting unnecessary messages telling you that you are using Microsoft’s compiler and the /EHsc flag is weirdly required for consistent exception-handling behaviour. We have also removed all compiler warnings, i.e. the equivalent of -Wall on UNIX can be achieved with /Wall on Windows. However, if you value your sanity, I’d strongly suggest to not use this flag. If you do, you’ll get a flood of warnings from the standard template library (which is implemented by Microsoft themselves, job well done!)
Both the buildAndRun.sh and buildAndRun.bat file will compile and run the main executable for us so that we can check any output printed to the console. With these files implemented, we are ready to look at our heat equation implementation.
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
A first solver implementation of the heat equation
If you want to become a serious programmer, you’ll have to embrace software testing and make it a habit to follow test-driven development (TDD) practices. I’ll cover that in a future series, but in short, TDD requires us programmers to write a test for each piece of code that we intend to write. We write the test before we write code that will make that test work, and then we start to implement the code. We implement and refine our code until our test starts working, and then move on to the next test/implementation.
There are a few philosophies about how we should structure our tests, but all what these tests are doing is to ensure that different components at different level work. For example, we provide small test pieces (called unit tests) that only test functions, but then also provide integration or system tests to test single classes and their combined behaviour, respectively.
I’ll skip the proper test set up for now, but what I do want to borrow from this concept is a system test approach to write our library. I want to write the code for the 1D heat equation upfront so that we can compile our library code and check if it is working for our intended use case. If it does, chances are it will work on other problems as well. This approach also allows us to design the interface to our classes, i.e. the function names, their arguments, and in which order they need to be called. This makes it easier to implement the classes later.
Having said that, the main.cpp function, which includes the entire discretised 1D heat equation example, is given below. The code won’t compile as of now, as we have not provided any implementation for the library code. Once that is there and correctly implemented, you should see this code running without throwing an error message at the end.
#include "linearAlgebraLib/linearAlgebraLib.hpp"
// solve the heat equation implicitly of the form dT/dt = gamma * d^2 T/ dx^2
// over a domain L using the conjugate gradient methdod
// initial condition: 0 everywhere
// boundary condition: T(0) = 0, T(L) = 1
int main() {
// input variables
const double gamma = 1.0;
const unsigned numberOfCells = 100;
const double domainLength = 1.0;
const double boundaryValueLeft = 0.0;
const double boundaryValueRight = 1.0;
const double dx = domainLength / (numberOfCells);
// vectors and matrices
linearAlgebraLib::Vector coordinateX(numberOfCells);
linearAlgebraLib::Vector temperature(numberOfCells);
linearAlgebraLib::Vector boundaryConditions(numberOfCells);
linearAlgebraLib::SparseMatrixCSR coefficientMatrix(numberOfCells, numberOfCells);
// initialise arrays and set-up 1D mesh
for (unsigned i = 0; i < numberOfCells; ++i) {
coordinateX[i] = i * dx + dx / 2.0;
temperature[i] = 0.0;
boundaryConditions[i] = 0.0;
}
// calculate individual matrix coefficients
const double aE = gamma / dx;
const double aW = gamma / dx;
const double aP = -1.0 * (aE + aW);
// set individual matrix coefficients
for (unsigned i = 1; i < numberOfCells - 1; ++i) {
coefficientMatrix.set(i, i, aP);
coefficientMatrix.set(i, i + 1, aE);
coefficientMatrix.set(i, i - 1, aW);
}
coefficientMatrix.set(0, 0, -1.0 * (aE + 2.0 * aW));
coefficientMatrix.set(0, 1, aE);
coefficientMatrix.set(numberOfCells - 1, numberOfCells - 2, aW);
coefficientMatrix.set(numberOfCells - 1, numberOfCells - 1, -1.0 * (2.0 * aE + aW));
// set boundary conditions
boundaryConditions[0] = -2.0 * aW * boundaryValueLeft;
boundaryConditions[numberOfCells - 1] = -2.0 * aE * boundaryValueRight;
// solve the linear system using the conjugate gradient method
linearAlgebraLib::ConjugateGradient CGSolver(numberOfCells);
CGSolver.setCoefficientMatrix(coefficientMatrix);
CGSolver.setRightHandSide(boundaryConditions);
temperature = CGSolver.solve(100, 1e-10);
// the obtain temperature profile is a linear one of the form T(x) = x. Thus, we can compare it directly against
// the coordinate vector (which in this case acts as an analytic solution)
linearAlgebraLib::Vector difference(numberOfCells);
for (unsigned i = 0; i < numberOfCells; ++i) {
difference[i] += std::fabs(temperature[i] - coordinateX[i]);
}
// ensure that temperature has converged to at least single precision
assert(difference.getL2Norm() < 1e-8);
std::cout << "Simulation successful, final residual: " << difference.getL2Norm() << std::endl;
return 0;
}Lines 9-14 define some input variables and on lines 17-21 we define the vectors and coefficient matrix that we will need later to solve our linear system of equations. We construct a 1D mesh on lines 24-28, as well as initialise our vectors here. For the mesh, the integration points are located in between vertices. Thus, we first find the leftmost vertex with i * dx (at the first iteration, i is zero and so the vertex is located at x=0 and then add half the cell length to that vertex location with dx / 2.0.
Returning back to the discretised equation we developed in last week’s article, repeated below for convenience:
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)=0we see that we can define our coefficients on lines 31-33 based in the above equation, i.e. with the coefficient for T_{i} being labelled a_P, and the coefficients for T_{i+1} and T_{i-1} being labelled as a_W and a_E, respectively. You’ll also notice that the coefficients on lines 31-33 do not contain the area of the cells A, since this is constant and thus can be divided in all terms so it will cancel out. We can also state that a_P= -(a_E + a_W), which is shown on line 33.
On lines 36-40, we set the coefficients of the matrix but we ignore the first and last row (for now), as these need different treatment. This is evident by looking at the matrix definition again from above:
\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}The first and last rows only contain two matrix coefficients, whereas internal rows contain three entries, and thus we treat them separately. Lines 42-43 add the missing first row coefficients, while lines 44-45 add the coefficients on the last row. Lines 48-49 set the right-hand side vector, called b in the above matrix notation.
Once the matrix, temperature, and right-hand side vector are all set, we can create our linear system of equation solver on line 52, set the coefficient matrix, as well as the right-hand side vector, and then solve this system for the temperature. We specify two arguments to solve() function, the first being the maximum number of iterations, and the second the convergence threshold below which residuals must fall.
On lines 59-67, we essentially check that the results are correct. The behaviour of the heat diffusion equation is such that given two boundary points, a linear and smooth profile will be developed between these two boundary points for a steady state case. If we constrain the left boundary to have a value of 0, and the right boundary to have a value of 1, and then further let the grid start at 0 and go to 1, then we can exploit the fact that at convergence, we should have the same temperature=coordianteX. This is shown in the figure below for different level of convergence.
Here, t_0 is our initial solution, and then each subsequent entry t_1, t_2, and t_3, represent some advancement in the solution. At t_3, we have found a steady state solution and we can see indeed how the temperature profile is a smooth variation between the left and right boundary value and that we have T=x. Thus, we are able to judge convergence in this way.
Summary
This article showed the basic library structure we’ll use for our library. Pretty much all libraries you’ll find out there will use a similar structure and so if we understand this week’s article, you’ll find yourself right at home once you start working with other CFD libraries. We briefly touched upon common software development practices in the form of test-driven development and software testing, and we provided our very own test that we’ll use to ensure our library is working as intended. With this prerequisite out of the way, we are ready to develop some library code!
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.