Lambda expressions are quite powerful in C++, but it is not immediately clear why we should use them over similar constructs such as function pointers. In this article, we look at what lambda expressions are, how they are related to function pointers, and specifically how we can apply them in some real CFD code examples. We end the article by looking back on the previous article, in which we discussed the standard template library (STL) and look at an example of how we can combine the STL and lambda expression. We will realise with this example how powerful lambda expressions really are.
In this series
- Part 1: Choosing the right programming language for CFD development
- Part 2: Why you should use C++ for CFD development
- Part 3: The complete guide to memory management in C++ for CFD
- Part 4: Object-orientated programming in CFD
- Part 5: How to handle inheritance and class hierarchies in C++
- Part 6: Templates in C++: Boost your CFD solver performance
- Part 7: Enhance readability with operator overloading in C++
- Part 8: The power of the standard template library (STL) in C++
- Part 9: Understanding Lambda expressions and how to use them in C++
- Part 10: Reduce memory bugs with smart pointers in C++
In this article
Function pointers: Assigning functions instead of values to variables
At this point, you know what a variable is. We use them all the time in our code, and, without them, programming as we know it would not be feasible. You have also learned about templates, which allow us to change the type of a variable at compile time, and we have looked at a few examples in the previously linked article, so give it a read if you haven’t already.
There is one more concept we need to understand, which is a powerful technique, yet it is heavily under-utilised, at least in the CFD community. I am talking about function pointers. That’s right, we can not just make variables hold different values, but also point to different functions. Let’s look at an example (which will not be very helpful to begin with, but stick with me, it’ll become clear later why we picked this example):
#include <iostream>
#include <functional>
double add(double a, double b) { return a + b; }
double subtract(double a, double b) { return a - b; }
int main()
{
// older C syntax
double (*cFunction)(double, double) = &add;
std::cout << cFunction(1, 2) << std::endl;
cFunction = &subtract;
std::cout << cFunction(1, 2) << std::endl;
// C++ syntax
std::function<double(double, double)> cppFunction = &add;
std::cout << cppFunction(1, 2) << std::endl;
cppFunction = &subtract;
std::cout << cppFunction(1, 2) << std::endl;
return 0;
}The code above will print the following content:
3
-1
3
-1On lines 4–5, we define two functions, one to add and one to subtract numbers. Then, in the main body, I have provided two implementations, one using the older C syntax, and the other using C++ syntax. On line 10, we define a pointer called cFunction, which is expecting a single return value of type double, and two function arguments, also of type double. We then set this function pointer equal to the add function, and on line 11, we see how we can now call cFunction, instead of add. On line 13, we see that we can change the function from add() to subtract(), and then when we call it again on line 14, the function call hasn’t changed, but the output has (second line of the output). On line 11, we see a result of 3 being printed, whereas on line 14, we see -1 being printed. Everything on lines 17–21 is exactly the same, except that we are using the C++ syntax here, which is less commonly seen in the wild, the reason being that function pointers are mostly replaced by lambda expressions (we’ll get to those in a second), or, because people are still programming in C and pretending it to be C++ code because they use a C++ compiler. Let’s not get me started on this one! … (but do feel free to give this short 5-minute video a try, this summarises the issue rather well)
Looking back at our function pointer example, we can see that the syntax is very similar between the older C syntax and the C++ syntax, i.e., on line 17, we define a function pointer cppFunction which has, again, a single return value of type double and expects two double arguments to the function call, and we set this equal to the add() function. On line 20, we set it to the subtract() function, and on lines 18 and 21, we call the cppFunction, again printing 3 and -1, respectively, as shown by the output.
Coming up with a CFD-relevant example: Boundary condition treatment
So, what’s the point of doing all this? Well, let me give you one good example where you may want to use this. Let’s look at the code first and then explore together what it is doing:
#include <iostream>
#include <functional>
#include <string>
void dirichlet(double value) {
std::cout << "dirichlet: " << value << std::endl;
}
void neumann(double value) {
std::cout << "neumann: " << value << std::endl;
}
void applyBoundaryCondition(std::string bcType, double value) {
std::function<void(double)> boundaryCondition;
if (bcType == "inlet")
boundaryCondition = &dirichlet;
else if (bcType == "outlet")
boundaryCondition = &neumann;
boundaryCondition(value);
}
int main()
{
applyBoundaryCondition("inlet", 0.57);
applyBoundaryCondition("outlet", 0.0);
return 0;
}This is an incomplete example but hopefully shows you enough to get an idea. The idea of this code is that we want to assign boundary conditions based on the type we have specified. For example, if we are dealing with the velocity vector, we may want to assign a Dirichlet-type boundary condition at the inlet and a Neumann-type boundary condition at the outlet. So, we define first the boundary condition types, i.e. dirichlet and neumann on lines 5-7 and 9-11, respectively. We also define a function called applyBoundaryConditions() on lines 13-22, and this function takes a std::string and a double as function arguments.
Within this function, we define a function pointer first on line 14, which expects a function with no return type (i.e. void) and a single function argument of type double. We then check which boundary condition we are using, e.g. inlet or outlet on lines 16-19, and then assign either the dirichlet() or neumann() function to our function pointer named boundaryCondition. On line 21, we call the right boundary condition based on the assigned type.
Within the main function, we call the applyBoundaryCondition() function on lines 26 and 27 with either an inlet or outlet condition, and executing this code will print the following in the console:
dirichlet: 0.57
neumann: 0We can see that the right function, either dirichlet() or neumann() was called in the background. Why is this a good example? Later when we are calling boundary conditions, we want to avoid code like:
int main()
{
std::string bcType("inlet");
if (bcType == "inlet")
dirichlet(1);
else if (bcType == "outlet")
neumann(1);
return 0;
}This type of code violates the open-closed principle again, as we would have to change the if/else statement as we add more boundary conditions. However, there is a more important reason than this; speed. If/else statements are slow, and we want to avoid them as much as possible inside loops unless they have a predictable pattern. For example, an if/else statement like the following may be tolerated:
int main()
{
std::string interpolationScheme("upwind");
// start loop over space and time here
if (interpolationScheme == "upwind")
upwind(1);
else if (interpolationScheme == "central")
central(1);
return 0;
}Again, we are violating the open-closed principle, but performance is not too badly affected. Why? Your runtime is smart and has a feature which is called branch prediction. This feature allows your runtime environment to guess which branch of the if/else statement is likely to be executed. Well, I say that the runtime environment is smart, in reality, it just checks which branch was used in the previous call, so if we called the first branch (e.g., here upwind) in the previous iteration, then your runtime environment will check the if condition for this branch first again. If it doesn’t change from iteration to iteration, as we would expect with something like a numerical integration or interpolation scheme, then we are essentially bypassing the entire if/else statement.
Back to our example, the boundary condition type does change, and we have no way of knowing which boundary condition is going to be executed next. Say we have a case with an inlet, outlet, and walls. Then, for each iteration, we have to loop over all three boundary conditions and apply them to each of our scalar and vector fields. The velocity will have a Dirichlet type for the inlet and wall and a Neumann type for the outlet. The pressure will have a Dirichlet type for the outlet and a Neumann type for the inlet and wall. So, the runtime would not be able to predict which function will be executed next (either dirichlet() or neumann(), it changes frequently) and so branch prediction will not help us if we try to solve this problem with an if/else statement.
However, if we assign boundary conditions before we start looping over them in space and time, then there is no more if/else statement inside the loop, and we just directly call the correct function when we get to the point of updating boundary conditions. The following pseudo-code should clarify this:
#include <functional>
#include <vector>
void dirichlet(double value) { /* ... implementation ... */ }
void neumann(double value) { /* ... implementation ... */ }
int main()
{
int numberOfBoundaryConditions = 3;
int numberOfIterations = 1000;
int numberOfCells = 5000;
std::vector<std::function<void(double)>> velocityBC(numberOfBoundaryConditions);
std::vector<std::function<void(double)>> pressureBC(numberOfBoundaryConditions);
std::vector<double> velocityBCValue(numberOfBoundaryConditions);
std::vector<double> pressureBCValue(numberOfBoundaryConditions);
// boundary condition 1: inlet
velocityBC[0] = &dirichlet;
pressureBC[0] = &neumann;
velocityBCValue[0] = 1.0; // inflow
pressureBCValue[0] = 0.0; // zero-gradient neumann condition
// boundary condition 2: outlet
velocityBC[1] = &neumann;
pressureBC[1] = &dirichlet;
velocityBCValue[1] = 0.0;
pressureBCValue[1] = 0.0;
// boundary condition 3: wall
velocityBC[2] = &dirichlet;
pressureBC[2] = &neumann;
velocityBCValue[2] = 0.0;
pressureBCValue[2] = 0.0;
// loop over time
for (int iteration = 0; iteration < numberOfIterations; ++iteration) {
// loop over space
for (int cell = 0; cell < numberOfIterations; ++cell) {
// calculate solution here ...
// at the end, update boundary conditions
for (int boundary = 0; boundary < numberOfBoundaryConditions; ++boundary) {
velocityBC[boundary](velocityBCValue[boundary]);
pressureBC[boundary](pressureBCValue[boundary]);
}
}
}
return 0;
}Lines 13 and 14 define a boundary condition vector that holds a function pointer to either dirichlet() or neumann(), which is set over lines 20 to 38. Also, we assign the value that is used at the boundary, which is stored in another vector defined in lines 16-17. There are better ways to define this, but I did not want to provide a full-blown class implementation of a BoundaryConditionManager class here. The important part is that on lines 48-51, we see how we execute the boundary condition treatment, i.e. we are within the time/iteration (line 41) and space loop (line 43), and our boundary condition treatment no longer requires if/else statements.
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
Towards lambda expressions: A for loop-based example
So the boundary condition treatment is just one example, but we can come up with an even better one. While the example with the boundary condition shows how we can get back some performance by avoiding if/else statements, we can also avoid code duplication, and this example will lead us to the main use cases of function pointers, and later lambda expressions.
One thing you will often find is the need to loop over your data in space. For example, if you calculate your gradients, you loop over space. Then you may assemble your right-hand-side vector, which requires an additional loop in space. When you do an explicit time integration, you loop again over it, and, if you want to output your data to a common file format, you loop over space again. See how code duplication will come up? Let’s put this into an example. For simplicity, let’s assume we have a three-dimensional, structured code, and that we are able to loop over space in the following way:
#include <iostream>
int main()
{
int nx, ny, nz = 10;
for (int i = 0; i < nx; ++i)
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k) {
// calcualte gradient, righ-hand-side, integrate, etc.
}
return 0;
}Chances are, if this is the code you have, you are going to copy and paste lines 7-9 a lot. And remember, every time you do copy and paste, regardless of how safe you feel doing so, you will likely introduce bugs in the future. So we want to isolate the looping part, have it defined in one place (and one place only!) and then reuse it whenever we need to loop over our data. By now, you should have all the knowledge you need to make this happen, so before reading on, look at the code above and try it yourself. If you have found a solution (or are impatient, I won’t judge), read on.
The answer is, unsurprisingly, function pointers. But now we have to expand our knowledge a bit. Keep in mind that function pointers are just regular variables, i.e. instead of pointing to a value in memory, they point to the entry point in a function. What that means for us is that we can pass function pointers to other functions as well. Let us write out the example and then look at what we have achieved:
#include <iostream>
#include <functional>
void loop(std::function<void(int, int, int)> function) {
int nx = 10; int ny = 10; int nz = 10;
for (int i = 0; i < nx; ++i)
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
function(i, j, k);
}
void gradient(int i, int j, int k) {
// calculate gradient, access cells with i, j, k indices
}
void integrate(int i, j, k) {
// integrate values in time for each cell
}
int main()
{
loop(gradient);
loop(integrate);
return 0;
}On lines 4-11, we have defined a function called loop which loops over all vertices or cells in our mesh. We would define line 5 probably somewhere else, but for the sake of simplicity, let’s keep it here. We also have two functions, called gradient and integrate, both of which require looping over the mesh and calculating values at each grid point. Instead of copying the loop construct on lines 7-9, what we do is specify that these functions take integer arguments i, j, k, and these correspond to the same i, j, k indices we see on lines 7-9.
For the loop function itself, we allow for a single function argument, which is defined as a function pointer, and then on line 10, we call this function with the i, j, k indices. What we have achieved with this is shown in lines 23-24, where we are now able to call the gradient and integrate the function in an abstract manner. No need to pass any indices to the function; the loop function will do that for us on line 10. And we, as the implementers of the gradient and integrate functions, only have to worry about how to calculate the gradient or integrate for a single cell, the location of which is coming from the i, j, k indices we received from the loop() function itself.
If you have understood this, then you are ready for lambda expressions, which are essentially the same, but with a twist. There are two differences between lambda expressions and the previous example. Lambda expressions are anonymous, and they are classes, thus allowing for an internal state. What does this mean? Well, look at the example above again, we had to specify the gradient() and integrate() function before we could use them, which makes sense. But what if you only ever use a function a single time?
Let’s say you want to have 2 loops in your gradient function, in the first loop you resize a std::vector to fit all your gradients and in the second you do the actual calculation. Would you now want to define a dedicated function for gradientSetup() and gradientCalculation()? If you push this further, you will come into situations where you have to loop over space several times within a single function, and so you may wonder if it is a good idea to write out function pointers for each individual loop. Instead, you may just want to provide an implementation for a function where it is used (i.e. in place and not somewhere else) and then pass it to another function. This is exactly what lambda expressions are doing and what is referred to as being anonymous.
The second part is that they capture a state, i.e. they are essentially classes. This is very useful and will become clearer once we look at more examples, but what this does allow you is to pass additional data to your lambda function without having to specify them as function arguments. Since lambdas are classes, passing additional data into the class that you can then use in the function is handled through the constructor of the lambda.
OK, enough talking, let’s look at our first lambda expression, one which is pretty useless but one that is so simple that we can concentrate on the syntax and how to call it.
#include <iostream>
int main()
{
auto myLambda = []() {
std::cout << "calling lambda" << std::endl;
};
myLambda();
return 0;
}We define our lambda expression on line 5 and note that the type here is auto. We set it to auto because we can’t know the type; this is only known to the compiler. You don’t have to worry about this part, just use auto, always, and you will be fine. A lambda expression consists of three parts, indicated by different braces; we have [], then (), and finally {}. The first part allows us to capture state, i.e. anything within the [] braces will be passed to the constructor of the class that is constructed by the compiler and then it will be made available within the lambda expression itself (now you can also understand why we have to use auto to store the lambda expression, as the compiler gets to choose the name for the lambda, we have no saying in it and so only the compiler knows the type, but not we (and we don’t need to)).
When the lambda expression gets called, it may have certain arguments passed to it (similar to our looping example, where we always passed the triplets of i, j, k to our functions). This part is done within the () braces. This is similar to any function call arguments and probably feels most natural about lambda expressions. Finally, we have the {} braces and this is where our actual implementation goes. In the example above, we are simply printing some text. We can then call the lambda expression as we do on line 9, but even if we can do that, this is actually not how we would use them, necessarily.
To understand this better, let’s make this example somewhat more real and look at an implementation you may find in an actual CFD solver. We’ll reuse the example of looping over the mesh, which now becomes:
#include <iostream>
#include <vector>
template<class FunctionType>
void loop(FunctionType function) {
int nx = 10; int ny = 10; int nz = 10;
for (int i = 0; i < nx; ++i)
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
function(i, j, k);
}
int main()
{
std::vector<std::vector<std::vector<double>>> velocityX;
velocityX.resize(10);
for (int i = 0; i < 10; ++i) {
velocityX[i].resize(10);
for (int j = 0; j < 10; ++j)
velocityX[i][j].resize(10);
}
loop([&velocityX](int i, int j, int k){
velocityX[i][j][k] = i * j * k;
});
return 0;
}Since we are now using a lambda expression, we can no longer know the type, so we had to remove our function pointer definition on line 5 and replace that with a template. The rest of the loop function remains the same. On lines 15-21, we define a 3D velocity vector, which we resize here to 10 in each direction. But on line 23, we are now using a lambda expression as the argument to the function loop. If you look carefully, you see that we define the lambda in-place, i.e. we are not even bothering to store the lambda expression somewhere and then pass the pointer to that lambda expression to our loop function, instead we define the lambda expression directly as the function argument of the loop function.
This is probably the most common use case, where we want to pass an arbitrary function to a function or class which then does some common operation first and then calls the lambda. In our case, the common operation is looping over the mesh in space, after which the lambda gets called. You’ll have noticed that we have passed the velocityX vector, which we defined on lines 15-21, to the lambda within the [] brackets (i.e. its constructor). We pass this vector by reference, indicated by the ampersand (&), i.e. the same rules for function arguments hold here as well. Within the lambda, we are modifying the velocityX vector and, after we are done and exiting this function, the velocity vector will have changed its values, and we can continue to work with the updated vector. What are some additional examples you may ask?
Catching the bad guy: lambda expressions in try/catch blocks
How about a class that checks for exceptions and automatically catches them? If you need a refresher or are completely new to exception handling, you may want to read up on them over at learncpp. In a nutshell, the runtime will throw an error if we are trying to do something that results in an undefined behaviour. For example, if we have a std::vector with 10 elements and then try to access its 11th element, that is undefined behaviour. We have no idea what is passed the memory of the std::vector and so you run into a runtime exception. Here is an example of a class that checks whether an exception is thrown:
#include <iostream>
#include <stdexcept>
class ExceptionChecker {
public:
template<class FunctionType>
static void check(FunctionType function) {
try {
function();
} catch (const std::exception &e) {
std::cerr << e.what() << std::endl;
throw;
} catch (...) {
std::cerr << "Caught exception: Terminating" << std::endl;
throw;
}
}
};
int main()
{
ExceptionChecker::check([](){
// do something silly here
throw std::runtime_error("being silly");
});
return 0;
}We are declaring a class ExceptionChecker, which only contains a single static function (if it is static, we do not have to generate an object from the class, but instead, we can use it directly. Instead of using the dot operator to call a specific function within the class, we use the double colon operator, e.g., ::). The function itself just checks if the function we pass into it will throw an error. If it does, we can check here for all different sorts of errors, but we are limiting ourselves here to only checking for standard exceptions (and then having a catch all block, indicated by the ...). In the main function, we call the check function directly on line 22, this time without capturing any state and using no function arguments. We are simulating here a silly mistake or bug by just throwing a runtime error, which is then picked up on lines 10-12, with the error message being printed to the screen. For completeness, this is the error message that you would see:
being silly
terminate called after throwing an instance of 'std::runtime_error'
what(): being silly
Program terminated with signal: SIGSEGVMind you, though, exception checking is very expensive, so much so that some companies, writing time-critical code, forbid the use of exceptions in their code. There is a lively debate within the community (here is a good discussion of why people wouldn’t want to use exceptions), and everyone seems to be glad to offer their opinion. I am not sure if this topic or tabs vs. spaces is more polarizing, but I digress…
A match made in heaven: Lambda expressions and the STL
There is one final thing I want to get to, and for that, we have to go back to the Standard Template Library (STL). If you have missed it, we discussed the STL in the last article, so give it a read if you haven’t to be up to speed. We needed to cover lambda expressions first because there are a range of functions implemented in the STL which require you to pass in a function or lambda expression. If you are unaware of it, then looking up the reference may be somewhat confusing.
Say you have generated an adaptive mesh refinement library, and you are constantly resizing your mesh. Perhaps every now and then you go through your memory and want to free up some of that. Since you are using a data structure from the STL, specifically a vector, you can easily use any of the defined algorithms in the STL. One of them is called std::remove_if, and as the name somewhat suggests, we want to remove an element from our data structure if it fits a certain criterion. To check if it meets the criterion, we need to pass a lambda expression to the std::remove_if function. Let’s look at this in code:
#include <iostream>
#include <vector>
#include <algorithm>
int main()
{
std::vector<int> meshID = {0, 1, 2, 3, 4, 5, 6, 7};
int largestID = 5;
meshID.erase(std::remove_if(meshID.begin(), meshID.end(), [largestID](int ID){
if (ID > largestID) return true;
else return false;
}), meshID.end());
meshID.shrink_to_fit();
for (const auto &ID : meshID)
std::cout << ID << " ";
std::cout << std::endl;
return 0;
}The output of this code is as follows:
0 1 2 3 4 5 We define some std::vector on line 7, here called meshID, which holds all indices to our cells. After resizing, we may have decided that we have lost a few and that the new, largest ID now is 5, as defined on line 8. On line 9, we make a call to the erase function, which is part of the std::vector data structure, and we pass the std::remove_if function as an argument, which will return iterators to the values that should be removed. We go over the entire range of the meshID vector, and the third argument is the lambda. It captures the largestID parameter, and itself gets passed each entry, one by one, as a function argument, here indicated by int ID within the () brackets of the lambda. We check if the ID we pass to the function is larger than the largestID, and if so, return true, which indicates to the std::remove_if function to return the iterator to that value, which then, in turn, gets deleted through the erase function. Printing the vector on line 16 will reveal that the largest value in the vector is indeed now 5, as we can see from the output above.
Summary
We looked at some use cases, but most commonly you will use it to either wrap it around some common code (as in the examples we looked at) or to call functions within the C++’s standard template library (STL). Lambda expressions may seem innocuous, but they are very powerful constructs, and understanding them and being able to use them will make your code better. The art is discovering in which places to use them. But my rule always is, that if I find myself tempted to copy and paste code, I try to think about how I can avoid this code duplication. Lambda expressions are just one more tool in our arsenal that will help us to avoid exactly that. If you want to get lost in all of the nitty-gritty details of lambdas, then I highly recommend Jason Turner’s playlist on lambda expressions.
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.