Home / Learn / Essential CFD Libraries / Config files: the interface between you, your solver, and the user

Config files: the interface between you, your solver, and the user

Everyone believes reading input files is a simple task. Are you one of them? Look closer, and you realise that there are issues you haven't even thought about; issues that may come back to haunt you long after you have settled on a specific file type in your solver to use.  

CFD solvers need to read input files: solver settings, meshing parameters, boundary conditions, the list goes on. Once a simulation is done, we need to store not just data, but possibly also some metadata, e.g. the parameters used in a simulation over a parametric range.

Most people will either use the file types they already know or the first that they found through googling. If we all picked our life partner the way we pick file types, we'd all be stuck with somebody named .tmp; technically functional, but nobody really wants them around long-term.

There is a lot of depth to a seemingly simple problem, and so, in this article, we will explore the file types that exist (and are in use) for reading and writing input/output files which are either produced or consumed by CFD solvers.

By the end of this article, you will know more about the best file types to use in CFD than, let's face it, 99.999999% of the population (approximately). This also means that once you know all the pros and cons of each file type, you know when to use which, and your solver will never run into unforeseen issues without you knowing why.

So, if you are ready to become the file type guru in your university or place of work, or you are just trying to settle a domestic dispute with your significant other (well done for finding someone willing to debate file types with you!) about the best file type, read on, the answers are all here!

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

[custom_category_posts_list category_slug="essential-libraries-for-cfd-solver-development"]

In this article

Introduction

OK, so I have sold you on the idea that file types aren't as straight-forward as it may seem. Let me give you an example. For that, we have to go back to 2011 (ah, pre-covid, we were so naive back then). I was nearing the end of my undergraduate studies and decided it would be a good idea to embark on a 6 Mm (that's Megametres, or 6000 km, or, for my American friends (which, by now, are mostly LLMs based on traffic reports), 236220500 inches) car journey to the north cape in Norway.

That trip lasted about 2 weeks and I have learned that a car can be quite small. And, since we were on the road for most of the time, and the car was full of stuff for 2 weeks to eat and drink, there wasn't much space either for just taking a break from each other. I learned it's better to travel with someone who gets CFD to never run out of conversation topics. On long trips, my daughters, unfortunately, only want to talk about Peppa Pig and Teenieping, never about von Neumann stability analysis (sigh).

After a week and many, many comfort breaks (you'd be surprised about the outdoor toilet culture in Norway, see this book (a book I do actually own, I have a great uncle with a sense of humour (he gifted it to me (after the two of us went on our very own outdoor trekking adventures where I encountered many of them (story for another time))))), we made it:

the north cape and a slightly adventurous pilot
Figure 1

We got there, and it seems pilots do fly however they want around this part of the world, and I am sure this was an EASA approved manoeuvre. While it is advertised as the most northern point on mainland Europe, it isn't, the furthest point north is a 15km walk away from it (and, as pedantic German engineering students, of course we went there as well).

The north cape is located in Norway, and Norway has a fundamental problem that is relevant to our article, which is also known as the Norway problem. I am pretty sure no one in Norway is aware of it, but I will let you in on a little secret here, which is so nerdy that only you and me will appreciate!

Each country has a country code, which is standardised by the International Organization for Standardization, or ISO. The Norway problem stems from the fact that Norway's ISO code is NO. How does that relate to file types? In computer science, we have boolean values, with can take exactly two states: either 0 or 1 (well, 0, or any other value, but usually, we use 1). This fits perfectly with our binary number system our hardware runs on.

Since no one can remember what 0 and 1 means, we have given it special names: true and false. We have false == 0 and true == 1. Now, I'd argue this is a pretty self-explanatory system, rooted in logic, and I think I can speak for everyone if I say that we all understand this. Two states, true or false, let's move on.

But no, some file type designers were looking at that and thought: "How can we improve a system that does not need improving"? These are the same people who were asking "How can we improve pizza?" and then came up with this answer:

a pizza with burgers for crust, what an abomination
Figure 2: Figure reproduced from QSR media

Yes, a burger crust, you are seeing that correctly. Based on my experiences with how Italians react to pineapple on pizza, this picture will give them a cardiac arrest (presumably anyone eating this pizza, too). I should probably ban Italian IPs from now on. I cannot be held accountable for any adverse health reactions.

If we use the definitions of attempted, first-degree, and second-degree murder, I'd argue this is a second-degree pizza. Malicious but not premeditated, i.e. I find it a bit conspiratorial to believe that the burger crust was invented only as a vehicle to carry out an Italian genocide, but that is for a judge to decide.

Not everything in life needs to be improved, it is OK to plateau and enjoy life as it is. But no, while the burger crusts designers are pursuing a complete Mediterranean basin redesign, file type designers committed an equally atrocious act by widening the definitions of what 0 and 1 can mean.

So, instead of accepting true or false, certain file types allow you to also state yes and no. This may seem like a well-intended move, but it also means that this single move has put Norway (with its ISO abbreviation NO) on a blacklist it didn't know it was on.

If no is equivalent to false, then country: "NO" could equally be interpreted as country: false. In some programming languages, like Python, for example, an empty string is interpreted as false, while a non-empty string is interpreted as true:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
country = ""
if country:
    print("will never print")

country = "NO"
if country:
    print("this will print")

country = False
if country:
    print("this will not print")

Do you start to see the problem? If we make a check based on whether the country was provided, and the file type interprets NO as false, then the poor people of Norway will never be able to pass your if clause without using a VPN.

On a completely separate note, I am all in favour of renaming the country of Yemen to Yesmen, this will really mess with legacy web code overnight. Someone should look into this ...

Granted, this example may not be very relevant for CFD solvers and their input file types, but don't you worry, here is a CFD-specific example. Imagine you have something like this in your input file:

1
viscosity: 1e-5

Still with me? Now let's look at how OpenFOAM likes to read their input files. They provide wrapper functions that look something like this:

1
viscosity = lookupOrDefault<float>("viscosity", 0.001);

So, here we try to read the viscosity parameter. If it is provided, great, we take its value, if it is not, we take a value of 0.001, or 1e-3. What could be a problem here? There are two. The first is the dependence on strings to do parameter comparisons. If my input file was:

1
viscosty: 1e-5

Well, viscosity and viscosty look similar but are not the same. So the lookupOrDefault() function looks for viscosity but cannot find it (because of the typo in viscosty), and uses a default value of 0.001, while you believed all of your simulations were using a value of 0.00001. Your Reynolds number is now off by a factor of 100, and you only realised this after you submitted and published your papers. Whoopsie daisy!

What you now have to rely on (and hope) is that, whoever has implemented the lookupOrDefault() function will actually check that the parameter you are looking for is actually available in the input files. If it isn't, then it should, at the very least, complain loudly and let you know that it could not find the parameter you have specified.

If the person implementing this function is particularly kind to its users, they may implement a fuzzy search algorithm (matching strings that are approximately the same), and so they could log something to the console like "viscosty" was not found. Did you mean "viscosity" instead?.

OpenFOAM sits somewhere in the middle, i.e. if you specify something it does not recognise, it'll pull the emergency break, stop the solver, and tell you that your input is malformed. In my view, that is the correct response to an incorrect input. It will also tell you what options are available, which has popularised the banana test among users.

String comparisons are always tricky, and something best avoided if you can. But, this isn't the worst that can happen. Let's use the correctly spelled name of the viscosity now, i.e. we actually have now:

1
viscosity: 1e-5

Since there is no decimal place, some file type standards will read this number (given in perfectly valid exponential/scientific notation) as a string. If we had, instead:

1
viscosity: 1.0e-5

Then we would be getting a floating-point number. If you don't know this, and you are using a file type standard (or parser that implements that specific standard), then you may never actually read the correct viscosity value. Even worse, your lookupOrDefault() function will now find a correctly spelled parameter called viscosity and so happily accepts its value. So the only place that was designed to catch any issues is now silently accepting a number that is in reality a string.

This is a classical example of an undefined behaviour, and, in the worst case, we have seen that this can bring down multi-million projects, as in the example of the Ariane 5 disaster. Worst case for us: the lookupOrDefault() function accepts our input but then immediately discards the value and uses the default value. We do a parameter sweep over viscosity (thinking that we changed the Reynolds number) and then confidently claim that there is no Reynolds number dependence on our results.

I am sure this will now send a few of you frantically looking through your results to check if that is the case for you. And, if it is you, this is a specific issue with YAML and its 1.1. standard. Whatever you use should already depend on YAML 1.2 or higher. If not, good luck to you.

See, I promised you in the beginning that file types can be rather tricky to work with, and I hope these two examples bring that point home. So now that we are sensibilised, let us start to look at the zoo of file types you may find in someone else's (or your) solver.

Input and output file types you should know about

In this section, I want to build up a bird's-eye view of file types. In order to understand why we ended up with all of these different types (many of which I am sure you have never even heard of), we need to talk about serialisation, and why it is of such importance in CFD. Many people have no clue about it, but once you understand it, you'll wonder how you ever survived without it.

What is serialisation and why you should care about it

CFD has one big problem: complexity. If all you ever want to do is to run 2D channel flows for incompressible fluids at a Reynolds number of Re=100Re=100, you may skip this section (and maybe you should question your life choices). But, say you wish to run simulations of anything north of 10 million cells, perhaps an unsteady simulation, and you have no idea what your time step is (as it is computed adaptively from the CFL number), then how do you know for how long the simulation will run?

We simply don't know, is often the answer. Even if it is a steady state simulation, we have no idea how many iterations we need. Last time I checked, whenever I submit a simulation to our HPC cluster, there is no option to say run for as long as the simulation needs mate. I typically have to specify how long I want to use the cluster, with a hard limit set to 5 days (at least on our own cluster in Cranfield).

Or worse, I know exactly how long it takes, and for the sake of argument, let's say it takes 7 days. I cannot run this simulation in one go, I will need to split it across 2 sessions. This means: start the solver, read the mesh, start the simulation, iterate, write the solution, stop the solver, resubmit the same simulation, start the solver again, read the mesh and the previous solution, iterate, write the solution, hope the simulation finished.

While this work-around may be Cranfield-specific (although I am sure other researchers elsewhere will have similar time-limitations for how long they can use their cluster), the restart-from-previous-solution part isn't. Restarting a CFD simulation from a previous solution is a fairly common practice.

If you have ever written your own solver, there will be a point where you start to think about implementing a restart file. If you have reached that point, but have never heard about serialisation, wonderful hacks start to emerge in your solver because you are trying to reinvent the wheel, finding a solution to a problem that has already been solved, over and over and over again!

But if you haven't, you may ask, what's so difficult about it? And, yes, you are right, on the surface, this is an easy problem. Store the solution variables (pressure, velocity, etc.) and then read them back when you need them, right? It couldn't be simpler. OK, let's say you implement the time stepping with a second-order backward scheme. It is fairly popular since it is easy to implement and has second-order accuracy in time. It reads:

∂ϕ∂t≈3ϕn+1−4ϕn+ϕn−12Δt(1) \frac{\partial \phi}{\partial t} \approx \frac{3\phi^{n+1} - 4\phi^n + \phi^{n-1}}{2\Delta t} \tag{1}

This is correct as long as Δt\Delta t is constant. If you, instead, implement this scheme in a solver where the time step is changing at each time level, then you would have:

∂ϕ∂t≈3ϕn+1−4ϕn+ϕn−1Δtn−1,n+Δtn,n+1(2) \frac{\partial \phi}{\partial t} \approx \frac{3\phi^{n+1} - 4\phi^n + \phi^{n-1}}{\Delta t_{n-1,n} + \Delta t_{n,n+1}} \tag{2}

Here, Δta,b\Delta t_{a,b} is the time step between time levels aa and bb. So, if you now store only the solution of your velocity and pressure, but not the time step information, i.e. Δt\Delta t, how will you continue your time stepping once you restart your simulation? If you only store the flow field, but not the time step information, you will see that later in your result, e.g. looking at a time series of a variable.

So, you need to store the time step. Not just one, remember, the backward scheme needs both Δtn−1,n\Delta t_{n-1,n} and Δtn,n+1\Delta t_{n,n+1}. So store both.

Ok, so we got the time step stored. What about residuals? We typically compute them as:

res(ϕ)=∑iN(ϕin+1−ϕin)2(3) res(\phi) = \sqrt{\sum_i^{N}\left(\phi_i^{n+1}-\phi_i^n\right)^2} \tag{3}

This is the L2-norm of the generic flow variable ϕ\phi (i.e. this could be pressure, velocity, temperature, etc.). The problem is that pressure, for example, may have rather large differences. Using SI units and assuming our simulation takes place around normal mean sea level operating conditions, we have pressure values of the order of 10510^5 Pa. At the same time, the velocity may be of the order of 1 m/s or so.

If we assume that both the pressure and velocity change by exactly 1% between iterations, we have 0.01⋅105=103=10000.01 \cdot 10^5 = 10^3 = 1000 for the pressure, and 0.01⋅1=0.010.01 \cdot 1 = 0.01 for the velocity. This is exactly what the residual is, i.e. it characterises by how much the flow variable changes between time steps (or iterations), in this case 1%.

If we now tried to plot this residual over time/iteration, we have 5 orders of magnitude separating our pressure and velocity residuals. To make the comparison of residuals easier for us, we typically normalise the residuals by the residual that was computed at the first iteration.

If we assume that the example given above was for the first iteration, then we could state that our normalisation factor for pressure would be pres,0=1000p_{res,0}=1000, and for velocity it would be ures,0=0.01u_{res,0}=0.01. If we now divide both of the residuals by our normalisation factor, we have 1000/pres,0=1000/1000=11000 / p_{res,0} = 1000 / 1000 = 1 and 0.01/ures,0=0.01/0.01=10.01 / u_{res,0} = 0.01 / 0.01 = 1. In other words, our residuals now all start at 1, and so, we can easily see how they compare against each other.

So, we want to restart our simulation. If we want to have a clean continuation of our residual computation, we need to store the normalisation factors pres,0p_{res,0} and ures,0u_{res,0} as well.

I could go on; the bottom line is, there is always additional data, let's call it the solver metadata, that we need to store. And you may say, ok, easy, I got this, I will just write a text file alongside my simulation, where I store all of this data, and move on. Sure, you could, but have you considered how you would do that? Let's stick with the residual normalisation factors, for a moment, and think the process through.

How would you store that data? Plain ASCII, or would you go binary? If you store it as ASCII, how many digits do you write out? A fixed amount? Do you let C++ decide? Do you use floating point or scientific/exponential notation? If you use binary, are you using big-endian or little-endian notation (and do you know how to change that if you need to)? Do you know about endianness or are you happy for the catastrophic bugs that will come your way if you don't?

If you run in parallel, are you storing one file per core or per simulation? If you restart your simulation, but now using a different domain decomposition, will that affect your restart information? What happens when you need to add additional information to your restart file, but some people are still using your old version of the solver, will they now be unable to restart their simulation because their restart files are outdated?

You may think this is a bit over the top. If you do, try it for yourself. Write out the normalisation factors as plain ASCII, however you please, and see if you get a convergence history without any bumps. Sure, you can achieve it, but only if you really thought about it beforehand. I have because I spent too much time trying to figure out why my residuals were looking strange after restarting.

So, restart files aren't necessarily the best way to store information, instead, you want to serialise your solver. Serialisation is the process of looking at what is currently stored in memory (RAM) and writing that to disk, so that next time you are running your solver, you get the exact same state in RAM (byte-for-byte identical).

If we want to have byte-for-byte identical copies, then we really are talking about binary files, not ASCII files. And, writing binary files in C++ isn't fun. Let's look at how we could natively support serialisation for our residual example we have looked at before:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// The residual class that can be used per flow variable
class Residual {
  double normFactor; // the norm factor we have looked at, other variables omitted here

public:
  // writing the norm factor byte-for-byte identical to a file on disk
  void serialize(const std::string& filename) {
    // open the file to which to write to, make sure it is a binary file
    std::ofstream out(filename, std::ios::binary | std::ios::out);

    out.write(reinterpret_cast<char*>(&normFactor), sizeof(normFactor));
    // if we wanted, we could write other stuff to disk here as well, omitted here

    out.close();
  }

  // we can read the same file again to restore the norm factor
  // this is known as deserialisation
  void deserialize(const std::string& filename) {
    // open file from which to read the serialised information
    std::ifstream in(filename, std::ios::binary | std::ios::in);

    in.read(reinterpret_cast<char*>(&normFactor), sizeof(normFactor));
    // we could read other variables here, too
      
    in.close();
  }
};

Here we have the Residual class that contains our normFactor, which we can write to disk through the serialize() function. We can also read the file again through the deserialize() function. In each, we open a file in binary mode (opening through std::ofstream and reading through std::ifstream), as indicated by std::ios::binary.

When we write to or read from a binary file, (i.e. out.write() or in.read()), we are using reinterpret_cast<>() here. If you are not a hardcore C++ person, it may be the first time you see it, and that is a good thing. We should never use reinterpret_cast<>() UNLESS we know what we are doing (it took me very long before I felt comfortable using it, I am not sure if I actually do feel comfortable with it).

C++ is a type-safe language, that means, each variable has a type and that cannot change (if we can ignore std::any for a moment, which is an odd addition for a type-safe language (you can change types here)). Each type has a specific size in memory. For example, a double typically stores variables with 8 bytes, while an integer stores variables with only 4 bytes.

But, not only do different types require different amounts of storage space (you can think of each byte as being represented by one bucket of memory), but their type also signals to the compiler how these bytes should be arranged. For example, floating-point numbers are written using the IEEE 754 standard, which is very different from how an integer is stored, or a class, or a pixel of an image, or ...

This is where reinterpret_cast<>() is coming in. It essentially tells the compiler "shut up, I know what I'm doing" and the compiler will allow you to do whatever you want (with your types). Here, we are converting now a variable (normFactor), which is a double, to a char, i.e. raw characters. It doesn't make sense to write a double as text, but, we have to do that because the write() and read() functions only know how to write chars.

Furthermore, by writing everything as a char, we essentially tell C++ to write all of our bytes that our double holds as it is, it should not interpret this variable as a double (and potentially change bytes when writing to disk). This removes the issue I mentioned earlier (i.e. potentially losing accuracy when we read/write our variables to disk; after reading we have a byte-for-byte identical copy).

One more thing, we are not actually reading or writing chars, but, instead, char *, i.e. a pointer to a char. I know, I shouldn't bring up pointers, everyone loves them, but, in a nutshell, if we write to a char, we can store exactly one character. If we write to a char *, we can write several characters (think: a word, or a sentence, or an entire article on cfd.university). This is just a C++ quirk we have to live with.

So, since we write to a char *, we also have to use a pointer double * as our source. But, since a pointer is just pointing to a specific memory address, we can use the address of operator in C++ to get that from our double, and that is why we provide &normFactor, and not just normFactor to reinterpret_cast<>().

The second argument to read() is the number of bytes to write. For that, we use the sizeof() function, which you can pass either a variable or a type directly. So, both sizeof(normFactor) and sizeof(double) would work equally well here.

Therefore, this example shows how we can write and read, byte-for-byte identical copies of variables, to and from disk. I should add that I have structured this class here in this way just for this example, but for a real solver, you want to extract the serialisation and deserialisation logic into its own class.

If we didn't do that, we would have to write the same serialisation/deserialisation logic for each class, violating the DRY (don't repeat yourself) principle. Each class should also only have one functionality (single responsibility principle); doing both residual calculations and serialisation is a violation of that principle.

There still is one big problem, and we have not addressed that yet. Endianness. When we write a variable, say, a double, which consists of 8 bytes, we haven't actually defined which way around we want to write it. If we imagine these 8 bytes being neatly placed in memory into 8 buckets of memory storage, how should we write and read them? Have a look at this 4 byte (32 bit) integer and how it could be written to memory:

little and big endianness example
Figure 3: Figure reproduced from Aeroid

We can either write/read a variable from the front (big endian) or from the back (little endian). Different hardware will use different conventions, although little endian seems to be what you will most often find these days. Just imagine you run your solver on a cluster, get all the results back, copy them to your PC to look at the results, and they are all just numerical noise. It could be that your cluster used a different endianness to your PC. It's something I have never come across, but that is something that can happen.

So, we also need to take care of the endianness if we want to avoid any future bugs (or release our code to users who potentially run your code on all sorts of exotic hardware, did you know that AMD released CPUs with only 3 cores?! See AMD Phenom 8750B Triple-core for reference).

I hope you start to see that file writing does have numerous smaller issues you may have never thought of before. The good news is, there are people in the world who are even more boring than me, and they took it upon themselves to provide us, the world, with elegant file types we can use so we have to never worry about any of the above again. What a wonderful (-ly boring) world we live in!

We already looked at the CGNS file type for storing CFD simulation data and the mesh, as well as HDF5 for storing solution data (both in binary), so it is time to add an array of additional file types to that list, some of which you will already know, and some of which you had no idea existed. One day, I shall create https://filetype.university. One day ...

A comparative view on different file types

We have now established some of the issues when we want to write data to disk. And, if you are one of the people who actually made it to this section in the article (presumably 95% of people starting to read this article have already left, it's just you and me now, what exclusivity, may I invite you to a Marks & Spencer pistachio Cookie on the house?), I have a treat for you (arguably one only you and I will appreciate!).

I have two tables for you. The first table shows you what features are supported by the file type, and the second table shows you how much it would cost you to integrate that into your existing/new solver (e.g. learning curve). Each column in the table uses a 3-star system, which is defined as:

Table 1
Symbol Meaning
★★★ Excellent: this is a reason to choose the format
★★ Adequate: works, with caveats you need to know about
★ Poor: either unsupported, or supported so badly that you want to give me your pistachio cookie back

The first table looks at the following properties to define its functionality:

  • Comments: Can the user provide comments/annotation to provide additional information for variables?
  • Standard: Is there a standard, or steering committee behind this file type, i.e. how confident can we be that the file type does not silently change overnight, and we have no idea that this has broken our code, without us having changed anything (besides updating to a new version of a file type parser library)?
  • Inheritance: Can we define a cascade of files that are all linked? E.g. can we define a base input file which is included in other input files so we only change variables that are changing from simulation to simulation, while keeping everything else that is constant in another base file that is included?
  • Expression & Units: Can I specify expressions (e.g. a time-varying boundary condition) and units (e.g. m/s, Pa, ...)?
  • Round-trip: Can I write data to a file, read it back, and get a byte-for-byte exact copy of what I just wrote?

There are quite a few options available to us in terms of file types, so, let us look at the functionality of the ones you would most likely find in the wild:

Table 2
Format Comments Standard Inheritance Expressions & units Round-trip
JSON ★ ★★★ ★ ★ ★★
YAML ★★★ ★★ ★★★ ★ ★
TOML ★★★ ★★ ★ ★ ★
INI / key–value ★★★ ★ ★ ★ ★★
XML ★★★ ★★★ ★★★ ★ ★★★
Embedded script (Lua / Python) ★★★ ★★ ★★★ ★★★ ★
HDF5 / CGNS / NetCDF ★★ ★★★ ★ ★ ★★★
CSV ★ ★ ★ ★ ★★
CBOR / MessagePack ★ ★★ ★★ ★ ★★★

I won't comment on any file type here and why I came up with this rating, we'll leave that for the individual file type sections. Instead, let us focus on what it will cost us to adopt any of these file types. For that, I have defined the following categories:

  • Error messages: Do you get clean, easy to understand and comprehend error messages when something goes wrong?
  • Build burden: If you want to support this file type, what dependencies do you have to bring in? Can you write a quick parser yourself, or include a single header file, or do you need to compile additional libraries, all with their own dependencies?
  • Learning curve: Can you pick up the file type in one afternoon or do you need to visit cfd.university to learn how to use it?
  • HPC / MPI Scaling: Does this file type scale to the most demanding simulations we may want to run in production? Does it work well for parallel simulation, does it have dedicated support for it, or will it introduce a bottleneck?

Based on these categories, I have come up with the following rating:

Table 3
Format Error messages Build burden Learning curve HPC / MPI scaling
JSON ★★ ★★★ ★★★ ★★★
YAML ★★ ★★ ★★ ★★★
TOML ★★★ ★★★ ★★★ ★★★
INI / key–value ★★ ★★★ ★★★ ★★★
XML ★★ ★★★ ★★ ★★
Embedded script (Lua / Python) ★★ ★★ ★★ ★
HDF5 / CGNS / NetCDF ★ ★ ★ ★★★
CSV ★ ★★★ ★★★ ★★
CBOR / MessagePack ★ ★★★ ★★★ ★★★

Again, we leave the detailed discussion to each file type's own section, which is up next.

A detailed look at CFD-specific file types

I have introduced quite a few file types above. In this section, I want to look at all of them and give you an idea for how to use them. Some of them are so rich that they deserve their own article (and CGNS and HDF5 already got their own mention here on the website), but we should at least get an idea how they work and how to adopt them in our code. Let's dive into the nitty-gritty part of this article.

JSON

If you ever had to work with an exchange format to store metadata, it is likely that you will have come across .json files already. I start with them because, in my view, it is likely the one easiest to get started with (though, spoiler alert, probably not the best one to stick with).

JSON stands for JavaScript Object Notation, and so the question is justified, why would we be dealing with something JavaScript specific if we are likely working with C/C++, Python, or dare I even say, Rust? Well, history (to the best of my knowledge anyway).

When I was in high-school, we had to learn about this new technology, called the internet (it had been around for some time at that point, but we are talking still pre-YouTube era). At that point, websites were mainly HTML for formatting, CSS for styling, and PHP for database integration. The look and feel of those websites really was something special, and this website has some screenshots that are pure nostalgia (but I am glad we have moved on).

While we had to learn how to code our own websites using these three frameworks (using notepad as our code editor), a new kid on the block was gaining in popularity: JavaScript. The history around JavaScript itself is fascinating, being created by a single person in 10 days, and that quick-and-dirty prototype developed in the 1990s has become one of the most used programming languages nowadays, which includes many of its original flaws (as you would expect from a quickly put together prototype).

JavaScript wasn't standardised, and we were living in a pre-google era, so Internet Explorer had a market share of 95% and, of course, Microsoft being Microsoft, saw JavaScript, adopted it, but changed it fundamentally. Because it had the highest market share, whatever Microsoft did was essentially the JavaScript standard (even though it originated from the Netscape browser (now Mozilla Firefox)).

It also meant that due to its incompatibility with pure JavaScript (and the rise of Chrome, which adopted pure JavaScript), the days of Internet Explorer were numbered. It's a bit ironic that Microsoft Edge is now essentially a Chrome clone and having no popularity whatsoever.

In any case, with the rise of JavaScript and application programming interfaces (APIs), we needed a mechanism to exchange data between servers. For example, your website wants to show the local weather for each user based on their region, well, then you first need to get the region from the user, send a request to a weather server, receive a response back, and then show it on your website.

This required some form of standardised format to send data, and, this is precisely what JavaScript offered. And so, the new standard in which information could be exchanged was the JavaScript Object Notation (JSON), which would later also become the basic building block for classes in JavaScript (not an original feature).

JSON has become a pretty standard exchange format, not just on the web, but in general for configuration files or information that needs to be stored in a structured way. JSON supports text (strings), numbers, boolean values, arrays, and key-value pairs (also known as hash tables (computer science), maps (C++), or dictionaries (Python), depending on your background (if you are a computer science person, I am probably insulting your intelligence by talking about JSON)).

An example JSON file may look like the following:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
{
  "meshing": {
    "numPointsX": 128,
    "numPointsY": 128
  },
  "boundaryConditions": {
    "north": "Dirichlet",
    "south": "Neumann",
    "east": "Neumann",
    "west": "Dirichlet"
  },
  "solver": {
    "turbulence": true,
    "model": "komega-SST",
    "CFL": [1e-6, 1, 10]
  }
}

Each JSON file starts with enclosing {} brackets. Then, we define everything as key-value pairs. In this example, we have 3 keys (meshing, boundaryConditions, and solver), and each of these keys has different values (meshing has two values numPointsX and numPointsY). We can have infinitely many nested key-value pairs. For example, we could have also specified the boundary conditions more explicitly as:

1
2
3
4
5
6
7
8
{
  "boundaryConditions": {
    "north": {
      "type": "Dirichlet",
      "value": 0.25
    }
  }
}

One thing that can catch you out is the requirement to always place a comma at the end of a key-value pair, except the last item. So, above, we have a comma after "Dirichlet", but not after "value". Equally, in this example we only have "boundaryConditions", which is the only outer key (and thus also last), so we don't have a comma here (unlike in the example above, where "boundaryConditions" is not the last item, so it does have a comma).

Thankfully, we don't really need to worry about it too much. If a library does the writing and reading for you, it will take care of that for you. And, if you have to write them manually yourself, you can always use a JSON validator, which will check your JSON file against the JSON standard and complain loudly if your file violates it. I don't have a favourite, but JSONlint seems capable enough.

Let's look at the solver settings for a moment:

1
2
3
4
5
6
7
{
  "solver": {
    "turbulence": true,
    "model": "komega-SST",
    "CFL": [1e-6, 1, 10]
  }
}

I have given you an example of how to use an array with the CFL number. It has three entries. What do these numbers mean? If I used this for my solver, and you had to use it to run simulations, and this is all you got, then you would either have to dig through the source code (if that is available) to figure out what these three values mean, look at the documentation (if it bothers to document what goes in here), use an LLM and hope it will not hallucinate, or, the approach 99% of people will take: guess and pray it is correct.

In this example, 1e-6 and 1 denote the lower and upper range for the min and max time step I want to allow, while the third argument is the CFL number itself. Thus, I don't want to allow time steps to be smaller than 1e-6 or larger than 1, regardless of what the CFL number-based computed time step is.

It would be really handy to have comments to communicate that with the user, but, since JSON was only ever designed to exchange information between servers (without humans ever reading them), comments were never part of the JSON standard. There are extensions (like JSONC (JSON with Comments)), but comments, by default, are not supported.

This brings us back to the tables we looked at before, we had for JSON's functionality:

Table 4
Format Comments Standard Inheritance Expressions & units Round-trip
JSON ★ ★★★ ★ ★ ★★

As we saw, comments are not a feature of JSON, hence 1 star only. It does have a solid standard behind it, so if you write JSON valid once, you can expect it to work across many different programming languages and operating systems. That's fantastic if we need robustness. There is no mechanism for inheritance. Say, for example, we have base.json with:

1
2
3
4
5
6
7
{
  "solver": {
    "turbulence": true,
    "model": "komega-SST",
    "CFL": [1e-6, 1, 10]
  }
}

If we now want to set up many different cases, in which this information is the same, we cannot include this base.json file. So, we cannot have channelFlow.json which includes base.json. If we wanted to do that, we would have to write our own code to support that (and, once we do, our JSON files become non-compliant and can only be understood and parsed by our solver). There are file types that allow for this natively, so there isn't any need to do that ourselves if we really need this feature.

There is no support for attaching units and expressions either. For example, the following would be really nice to have but does not work:

1
2
3
4
5
6
7
8
{
  "boundaryConditions": {
    "east": {
      "type": "Dirichlet",
      "expression": "min(1, t) * U_inlet m/s"
    }
  }
}

What I am doing here is to define a ramp-type boundary condition (which gradually brings in the inlet velocity U_inlet based on the current time t). It also states that the velocity should be expressed in units of m/s. Unit support, in general, is pretty poor in any programming language (and, by extension, file types), so don't get your hopes up.

Reading a JSON file and then writing it again usually works well, but floating-point numbers may not be exactly reproduced. Hence, only 2 stars out of 3.

Looking at the burden of adopting JSON for our code development, we had:

Table 5
Format Error messages Build burden Learning curve HPC / MPI scaling
JSON ★★ ★★★ ★★★ ★★★

As I mentioned previously, if you use a validator like JSONlint, it will give you error messages that will help you to identify what the issue is. They are helpful, but not perfect. The build burden is essentially non-existent. In Python, for example, JSON is part of the standard Python library, and you can use it without installing any dependencies.

In C++, there are pretty excellent libraries available to help you parse JSON files. The most popular one being nlohmann::json. You can download the entire git repository, or head over to the releases and download the single header json.hpp, place it in your project, and start using it. It also integrates well with Conan and CMake's fetchContent, so you should always be able to bring it into your environment.

It also has a pretty decent documentation, and going through the GitHub repository's first few examples is usually enough to know how to use it. Anyway, that's what I usually do.

Let's see how we can drop JSON into a C++ code. We'll write the same JSON file we have looked at before with the meshing, boundaryConditions, and solver key-value pairs. I'll concentrate on the code here only, but I have detailed build instructions in the zip file you can download at the beginning of this article on how to integrate nlohmann::json with your project using Conan and CMake.

First up, the following code creates the config.json file:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#include <fstream>
#include <iostream>
#include <nlohmann/json.hpp>

int main() {
  nlohmann::json config;

  config["meshing"]["numPointsX"] = 128;
  config["meshing"]["numPointsY"] = 128;

  config["boundaryConditions"]["north"] = "Dirichlet";
  config["boundaryConditions"]["south"] = "Neumann";
  config["boundaryConditions"]["east"]  = "Neumann";
  config["boundaryConditions"]["west"]  = "Dirichlet";

  config["solver"]["turbulence"] = true;
  config["solver"]["model"] = "komega-SST";
  config["solver"]["CFL"] = { 1e-6, 1, 10 };

  std::ofstream outFile("config.json");

  outFile << config.dump(4);
  outFile.close();

  return 0;
}

We first include the nlohmann/json header. And, thanks to this header name, I will always remember that Niels Lohmann wrote the library nlohmann::json. And now you will, too. I mean, look at his GitHub bio:

Niels Lohmann's bio page on GitHub
Figure 4

Yes, Niels, I do know you from your JSON library, thank you for reminding me. Just what we need, yet another solipsistic German on the internet ...

We declare an object of type JSON as nlohmann::json config. That's it, now we can add all sorts of key-value pairs to this config object. Let's look at the meshing key for a moment. In our JSON file, we had:

1
2
3
4
"meshing": {
  "numPointsX": 128,
  "numPointsY": 128
},

This maps to:

1
2
config["meshing"]["numPointsX"] = 128;
config["meshing"]["numPointsY"] = 128;

C++ isn't always the easiest, or prettiest to read, but in this instance, I find this mapping to be rather straightforward and easy to understand. If we have nested key-value pairs, these are just as easy to implement. For example, before we looked at this example:

1
2
3
4
5
6
7
8
{
  "boundaryConditions": {
    "north": {
      "type": "Dirichlet",
      "value": 0.25
    }
  }
}

In C++, we would express that as:

1
2
config["boundaryConditions"]["north"]["type"] = "Dirichlet";
config["boundaryConditions"]["north"]["value"] = 0.25;

Let's look at how we can read a JSON file. It's just as easy as writing a JSON file, but there is a bit of added effort to satisfy C++'s strongly typed language requirements:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <fstream>
#include <iostream>
#include <nlohmann/json.hpp>

int main() {
  std::ifstream inFile("config.json");
  nlohmann::json config;
  inFile >> config;

  int numPointsX = config["meshing"]["numPointsX"].get<int>();
  int numPointsY = config["meshing"]["numPointsY"].get<int>();

  std::string north = config["boundaryConditions"]["north"].get<std::string>();
  std::string south = config["boundaryConditions"]["south"].get<std::string>();
  std::string east  = config["boundaryConditions"]["east"].get<std::string>();
  std::string west  = config["boundaryConditions"]["west"].get<std::string>();

  bool turbulence = config["solver"]["turbulence"].get<bool>();
  std::string model = config["solver"]["model"].get<std::string>();
  std::vector<double> cfl = config["solver"]["CFL"].get<std::vector<double>>();

  return 0;
}

The slightly more challenging part here is that we don't have any native support in the nlohmann::json library to read JSON files directly (bad Niels). So, instead, we have to first open a file ourselves (std::ifstream inFile("config.json");), then create a JSON object (nlohmann::json config;), and then we can dump the content of inFile into our config object using the stream operator inFile >> config;.

If you have done C++ for long enough, this will feel natural and normal. With my 10,000+ hours of C++ development under my belt, I still find that syntax slightly irritating. Maybe one day I will feel comfortable with C++. Maybe ...

In any case, now that we have read the JSON file, we can access it mostly in the same way as we did before when we were writing the files.

The main difference now is that we have to declare the variable type which we want to store (e.g. int numPointsX), which needs to be the same type we read from the config file, e.g. config["meshing"]["numPointsX"].get<int>(). Here, we say that the information we want is located at config["meshing"]["numPointsX"], but we also have to specify that this variable ought to be read/interpreted as an integer (get<int>()).

Now, the C++ steering committee looked at that and said "We are defining the type twice here, what a waste!". Thanks Bjarne and friends. So, if we wanted, since 2011, with the introduction of the C++11 standard, we could replace all type information on the left with auto. So, the following C++ code is perfectly legal:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <fstream>
#include <iostream>
#include <nlohmann/json.hpp>

int main() {
  std::ifstream inFile("config.json");
  nlohmann::json config;
  inFile >> config;

  auto numPointsX = config["meshing"]["numPointsX"].get<int>();
  auto numPointsY = config["meshing"]["numPointsY"].get<int>();

  auto north = config["boundaryConditions"]["north"].get<std::string>();
  auto south = config["boundaryConditions"]["south"].get<std::string>();
  auto east  = config["boundaryConditions"]["east"].get<std::string>();
  auto west  = config["boundaryConditions"]["west"].get<std::string>();

  auto turbulence = config["solver"]["turbulence"].get<bool>();
  auto model = config["solver"]["model"].get<std::string>();
  auto cfl = config["solver"]["CFL"].get<std::vector<double>>();

  return 0;
}

I tend to prefer this way, not just because it is shorter (and, if you have a template heavy code, auto becomes the only sensible mechanism to declare your types!), but also because it is more precise. For example, the following line will compile and is perfectly valid:

1
double numPointsX = config["meshing"]["numPointsX"].get<int>();

It will compile, it won't even issue a warning during compilation (at least on Microslop's MSVC compiler). But, what did you actually mean? From context, this probably ought to be an integer, but internally, it is now handled and stored as a double, which may break your code later in ways you can't foresee now.

auto solves that issue by having one, and only one type definition. I'd say it is more mathematically exact.

I once got into an argument with a software engineer at the Mercedes-AMG F1 team over the use of using auto. I told him my view that I just shared with you, but he was having none of it. auto was bad and seemingly discouraged within the team.

You would think that a German car manufacturer would use the word auto wherever they can (auto in German translates to car in English), but, the engineer I spoke to was British, so I suppose he did not get that point. It was 2022, the beginning of the downfall of their dominance in F1, and I am pretty sure they lost the following championships because they banned auto.

And, while we are on the topic, shout out to my ex-neighbour James (name changed), who, at least in 2022, worked for the Mercedes-AMG F1 team. I moved to a new home in 2020 (just at the start of the first COVID lockdown in the UK) and so, like everyone else, was stuck at home, and so was James. He and his wife had, shall we say, some issues.

The only reason I know his name is James is that his wife used to shout his name from wherever she was outside the house. James would then magically appear. By the time we moved out in 2022 (just after COVID), James was busy building a home office in his garden. From all the verbal exchanges he and his wife had in the open, I believe it also served as his new bedroom. Poor James. (I just checked Google Maps, his second home is done, and it is taking 25% space of the garden. Fantastic!)

There was another couple across the street who kept a pet pig in their house, and I saw them walk the pig a few times. A bit too eccentric for my taste, but it was COVID, we all did silly things, right?

Moving on, the only other thing worth pointing out here is how we can read arrays. With the CFL number, we saw that we can read them through a std::vector:

1
2
3
4
5
6
7
8
9
// if you work at Mercedes F1
std::vector<double> cfl = config["solver"]["CFL"].get<std::vector<double>>();

// if you work anywhere else
auto cfl = config["solver"]["CFL"].get<std::vector<double>>();

// if you use claude/chatgpt/gemini
auto prompt = "";  // 3 hours later: 400 lines, 12 dependencies, still doesn't compile
auto prompt = "idk just make it work like James' marriage (barely, and in public)";

And, honestly, with what we have covered, you are now a certified black-belt JSON master. We haven't looked at the null type (because it is unlikely you'll ever need it) or at mixed variables within arrays. For clarity, probably not a good idea to mix types anyway, but if you have to, std::tuple is what you will need here.

But, with that, let's move on.

YAML

According to Wikipedia, a markup language is a text-encoding system which specifies the structure and formatting of a document and potentially the relationships among its parts. JSON satisfies this description, and we think of JSON as being a markup language. Do we need another one? Well, Clark Evans, Oren Ben-Kiki, and Ingy döt Net certainly thought so and joined forces to create Yet Another Markup Language, or YAML in short.

It seems, though, Evans, Ben-Kiki and döt Net (or ebk.net as I shall refer to them) seem to have since had a change of heart and did a 180-degree turn and rebranded to YAML Ain't Markup Language (still abbreviates to YAML, how convenient).

In 2001, British Petroleum (BP) rebranded to Beyond Petrol (25 years later, 60%+ of their revenue is still from fossil fuel related products, perhaps the next rebrand will see them adopt the name Back to Petroleum). This was after merging with other companies, but also as a result of some, shall we say, questionable business decisions, like creating the largest oil spill (on purpose).

It makes you wonder what catastrophic business decisions ebk.net took to require such a drastic U-turn. One can only speculate. YAML was introduced in 2001. BP rebranded in 2001. Coincidence?

Regardless of what they may, or may not, be on the run for, YAML, in a nutshell, is JSON without the braces, commas, and quotation marks. Our JSON example of a sample solver configuration is shown below in YAML format:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
meshing:
  numPointsX: 128
  numPointsY: 128

boundaryConditions:
  north: Dirichlet
  south: Neumann
  east: Neumann
  west: Dirichlet

solver:
  turbulence: true
  model: komega-SST
  CFL:
    - 1e-6
    - 1
    - 10

So, in JSON, we had:

1
2
3
4
5
6
{
  "meshing": {
    "numPointsX": 128,
    "numPointsY": 128
  }
}

Remove all {}, ,, and "", and you have YAML. Sort of. We see that the CFL number is given as:

1
2
3
4
CFL:
  - 1e-6
  - 1
  - 10

But, I could have also written that as:

1
CFL: [1e-6, 1, 10]

Of course, as we have discussed in the introduction of this article, we need to make sure we use the YAML 1.2 standard to avoid any nasty surprises with the 1e-6 value (or, to be safe, use 1.0e-6).

Nested key-value pairs are just as straightforward. We can nest additional key-value pairs inside an existing key. The example we looked at before with additional information for each boundary condition may be written as:

1
2
3
4
boundaryConditions:
  north: 
    type: Dirichlet
    value: 0.25

One of the strengths of YAML (and, where it probably makes sense to decide over which file type to support, beyond how JSON and YAML files look and feel) is that it supports a kind of inheritance. It does so by providing anchors and aliases. This is best shown through an example.

Let's use the same solver definition from before, but let's say we want to perform a grid dependency study. For that, we want to generate three grid levels using 64264^2, 1282128^2, and 2562256^2 cells. We can specify a base setup that is shared among all runs, and then only change the meshing parameters as shown in the following example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
case: &caseBase
  meshing:
    numPointsX: 128
    numPointsY: 128

  boundaryConditions:
    north: Dirichlet
    south: Neumann
    east: Neumann
    west: Dirichlet

  solver:
    turbulence: 1
    model: komega-SST
    CFL:
      - 1e-06
      - 1
      - 10

cases:
  coarse:
    <<: *caseBase
    meshing:
      numPointsX: 64
      numPointsY: 64
    
  medium:
    <<: *caseBase
    meshing:
      numPointsX: 128
      numPointsY: 128
  
  fine:
    <<: *caseBase
    meshing:
      numPointsX: 256
      numPointsY: 256

First, we define a common case setup. We give this case set up an anchor using the &caseBase identifier. This means we can now refer to case (and all of its sub-entries) through the variable (anchor) &caseBase. We define the case setup as before, but then, we also introduce a new key-value pair called cases.

In cases, we have three separate cases; coarse, medium, and fine. Each of these cases is inheriting the base case through <<: *caseBase. Here, *caseBase is an alias and references whatever was defined inside the anchor &caseBase. The <<: syntax simply states that we should take everything from the alias (caseBase) and dump it here.

Once we have done that, we overwrite the meshing node by specifying new values for numPointsX and numPointsY. We don't have to respecify everything that is stored in meshing, we can overwrite only those parameters we want, in this case it just so happens that we overwrite everything within meshing.

Also, I am being overly verbose here, the medium case could have also just been defined as:

1
2
medium:
  <<: *caseBase

The caseBase already defines the mesh to contain 128 points in each direction. I am including it here as you may want to use different spacing values for the grid dependency study.

It is a small feature, but if you want to run parametric sweeps and define that sensibly in an input file, this is one way to do it. Unfortunately, we can't split cases into their own separate YAML file and then include them in a base YAML file. If we wanted to do that, we would have to provide that logic in our code, but there is no native YAML support to include files inside another file.

And, just as with JSON, there isn't really much more to it. But let us look at the rating we gave YAML in our overview table and see how that rating came about:

Table 6
Format Comments Standard Inheritance Expressions & units Round-trip
YAML ★★★ ★★ ★★★ ★ ★

Yes, YAML provides native support for comments. The following is perfectly legal:

1
2
3
4
5
# boundary condition node
boundaryConditions:
  north: 
    type: Dirichlet
    value: 0.25

While this is legal, anyone writing comments like this should be referred to a mental hospital. So, support is there, three out of three stars (I sound like a participant on Dinner Date now ... what an absolute car crash, almost as good as Naked Attraction (I don't have the guts to explain that concept to you, perhaps Suzi Ruffell can do that more eloquently)).

We have already discussed the issues with the YAML standard in the introduction, I touched upon it here as well, so, while there is a standard which is solid these days, there may still be some arcane libraries which have implemented an older standard and thus may be prone to the issues discussed. So two out of three stars.

Inheritance works great out of the box using anchors and aliases, so it scores highly here, but just like JSON, there isn't any support for expressions and units. For round-tripping, it also scores low, same issues as JSON, but also inconsistent standards can pollute your files (writing and reading with different YAML parsers can give you different results). Moreover, there are multiple ways to write arrays, so you may end up with an entirely different structure. All perfectly valid YAML, but different files nonetheless.

In terms of the error message and learning curve, we had:

Table 7
Format Error messages Build burden Learning curve HPC / MPI scaling
YAML ★★ ★★ ★★ ★★★

Error messages are usually good, have a look at YAML Lint. Write a non-compliant YAML tree, and it will give you error messages which is the offending line. The build burden depends on the library you are using. Now, this is somewhat depending on the library you use for parsing YAML files, but I have found it to be slightly more complicated, still manageable though, to make it work.

If you have a look at Conan Central for available YAML libraries, you have some choice, but some libraries are better than others. I have picked rapidyaml to test C++ integration (it is my first time to actually parse YAML files myself; otherwise I have only used them as config files to be parsed by another program). It seems to have the least brain-dead implementation of YAML, but an array of other issues (which I will get on to).

Long story short, it does compile well with Conan, but I had to troubleshoot a few issues before it could work. Two out of three. In terms of learning curve, if you want to master inheritance, it takes you slightly longer to learn YAML compared to JSON, and also, the available YAML parser libraries aren't great, as discussed, so only two stars here. HPC and MPI scaling is the same as JSON, so no difference here.

OK, so let's look at how we can use rapidyaml in C++ to parse YAML files. We will generate the same solver input file again, and the C++ code to do that is shown below:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include <ryml.hpp>
#include <ryml_std.hpp>

#include <fstream>
#include <string>

int main() {
  c4::yml::Tree tree;
  auto root = tree.rootref();
  root |= c4::yml::MAP;

  // meshing
  auto meshing = root["meshing"];
  meshing |= c4::yml::MAP;
  meshing["numPointsX"] << 128;
  meshing["numPointsY"] << 128;

  // boundaryConditions
  auto bc = root["boundaryConditions"];
  bc |= c4::yml::MAP;
  bc["north"] << "Dirichlet";
  bc["south"] << "Neumann";
  bc["east"]  << "Neumann";
  bc["west"]  << "Dirichlet";

  // solver
  auto solver = root["solver"];
  solver |= c4::yml::MAP;

  solver["turbulence"] << true;
  solver["model"] << "komega-SST";

  auto cfl = solver["CFL"];
  cfl |= c4::yml::SEQ;
  cfl.append_child() << 1e-6;
  cfl.append_child() << 1.0;
  cfl.append_child() << 10.0;

  // Emit YAML to std::string
  std::string yaml = c4::yml::emitrs_yaml<std::string>(tree);

  std::ofstream out("config.yaml");
  out << yaml;

  return 0;
}

I am not happy about that. Not because of the abusive use of streaming operators (which is fine, albeit somewhat unorthodox, why can't we just use simple map assignments like we did in JSON?), or the definitely unseen use of the or-equal |= operator (which is the first time I have ever seen it in use), but because the above code is generated with ChatGPT (with edits on top from my side).

Now, when it comes to learning a new framework or library, I am old school. I still go through documentation, which, presumably, was written by a human for other humans. This may change in the future. And, for those of you that have been for long enough on my email list, you will know that I have a certain fetish for documentation. Heck, I have an entire series on how to document your code correctly, which almost nuked this website (as in, I lost half of my readers overnight when I started publishing that series).

Accepting the possibility that you will leave, die of boredom, or hit your monitor with your keyboard in rage, let me walk through the basics of documentation. This is what it should do:

  1. Tell me how I can install this framework or library.
  2. Give me the TL;DR quick-start guide I need to use the most common and useful features of this framework/library. It can be a series of small tutorials.
  3. An in-depth description of what the framework/library can do, with all its features. It can be source-code level documentation, I don't care, but just something that shows me what else I can do that isn't covered by the quick-start guide.

Let's look at rapidyaml's documentation, which is hosted on Read the Docs. This is the table of contents:

  • Doxygen docs
  • YAML playground
  • Quick links
  • Is it rapid?
  • YAML standard conformance
  • Using ryml in your project
  • CMake build settings
  • Other languages
  • Alternative libraries

Doxygen docs points me to the documentation (wait, what, I thought we are already at the documentation). Yes, the documentation points you to the real (?!) documentation. And, once you land on the Doxygen page, it points you back to the read the docs' documentation. Yes, we have a loop, and no documentation. I am not sure if this is some kind of insider joke, I am not getting it (but I am close to rage-quit myself!).

The YAML playground supposedly allows me to test different YAML libraries head on, rapidyaml being one of them. So you try to run it, only to realise that you have to spin up a sandbox server yourself, and then go to a different, locally hosted website to test it. I hope you have Docker installed, and know how to use it; otherwise you cannot use it. What is the point of having a website that tells you to run your own server only to see a clone of the exact website locally?

Quick links points you, again, to the Doxygen documentation, the GitHub's README.md file and the changelog, all things I do not care about in a documentation (unless the README.md file actually contains documentation (which it doesn't!)).

Is it rapid? is self-praising on the highest level that would make even Andrew Tate look like an amateur (that is also the same content you'll find in the README.md file). Again, I don't care how much faster this library is compared to others. I am parsing files of a few lines, maybe even a few hundred lines, but if that takes 10ms or 100ms, who cares?

YAML standard conformance, my favourite section, starts with "ryml is fully conformant and feature-complete with regards to the YAML specification". The heading which follows literally 3 sentences later is called: "Deliberate deviations". This section starts with "ryml deliberately makes no effort to follow the standard". I couldn't make this up myself.

Using ryml in your project's first heading is "Quick-start build samples". I read that with great joy, finally, some samples, something to show me how to use it. Then you look at the samples and realise, well, it really should be sample (singular), there is just one; the quickstart.cpp file. But the pattern has been established by now, so if I told you the quick-start file, which is supposed to show you a high-level overview of the library, has 7413 lines of code, I think we can agree that this is taking the piss. But it is on-brand!

CMake build settings is potentially helpful in getting this library integrated into your own CMake project (that is, if you enjoy the pain). Just use a package manager like Conan, and you won't have to look at that until something breaks (which it rarely does).

Other languages and Alternative libraries talk about rapidyaml's integration into Python and JavaScript, and why you shouldn't use other YAML libraries.

At no point did the documentation do its job and actually tell me how to use the library; all I want is to write and parse some YAML file. It shouldn't be that difficult.

But, alas, you click through the Doxygen documentation, you go through the nested headers, you find a section that seemingly does what you want, and there, hidden behind Doxygen's beautiful 1990s HTML table structure, some source code emerges, showing you how to use the library. Fantastic, my prayers have been heard. So, I located the part which allows me to create a tree, which then will be used to store all the entries we will later see in the YAML file, which is hidden in the samplecreatetree_style.

But, I don't really mind as much, I am ok to lose some time by going through the documentation, even if I cannot find what I am looking for immediately, I read other things and learn adjacent topics, so it all helps to get started with a library.

Finally, some code, the first three lines:

1
2
3
ryml::Tree tree;
ryml::NodeRef root = tree;
root.set_map(ryml::BLOCK);

I see, so this is how you are supposed to write a tree. And then, presumably, I can add to it with the following syntax (taken from the example before in the documentation, i.e. samplecreatetree):

1
root["pi"].set_serialized(ryml::fmt::real(3.141592654, 5));

OK, so let's make that a valid C++ program and try to get the tree to work:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
#include <ryml.hpp>
#include <ryml_std.hpp>

int main() {
  ryml::Tree tree;
  ryml::NodeRef root = tree;
  root.set_map(ryml::BLOCK);

  return 0;
}

Let's compile:

1
C:\essentialLibraries\03_IO\src\writeYAML.cpp(8,10): error C2039: 'set_map': is not a member of 'c4::yml::NodeRef' [C:\essentialLibraries\03_IO\build\writeYAML.vcxproj]

What little code and examples are there, does not work! Here is one caveat, though: The example is taken from the 7413-line quick-start guide. The quick-start guide defines all sorts of helper functions, potentially something which I need to do as well to make this work, though if that is the case, it is even worse! I am using data types that are supposed to be baked into the library, that should compile regardless. But it doesn't!

And this is why I am not happy about the C++ example. When I look at it, I cannot find anything resembling the syntax ChatGPT suggests I should be using in the documentation. But, the documentation is crap. Honestly, an LLM generated documentation at this point would be more helpful than the self-praising, North Korean self-propaganda that is the current state of the documentation.

I don't think we can use the word documentation here. It's more like the library's Tinder profile, boasting all of its qualities while there is no substance behind it. Call it what you want, but it ain't documentation, and that upsets me!

My feelings aside, that just points to a shift in how we program and I think it is worth pointing that out explicitly. We as humans go the path of least resistance. If the documentation is crap, and an LLM is waiting for your input, and it is free (well, currently it still is, let's see how that will change), users will resort to prompting their favourite LLM and take code unquestioned. How are they supposed to know if it is right or wrong if there is nothing to check it against?

This is a new reality which we have to accept will become the norm, and in the future, we will unlikely read documentation ourselves any more. I mean, MCP servers have only been recently introduced (and, if you read this in the future, when MCP is mainstream, you may wonder why I even bother to point this out) and they already read the latest documentation of the library we want to use, so that their code generation is grounded in the documentation.

So, perhaps, we will never read documentation ourselves any more, who knows. But, all I can say is that with the speed at which LLMs provide us with code, and the pace at which we have become dependent on them, it is fascinating to pick up a book about programming now and read it. All of a sudden, everything slows down, things start to make a lot more sense (if you have a good book!), and it almost feels like meditation.

I may sound weird, but I am riding the LLM wave to see how it can help to benefit me, while I also deeply care about programming as a craft. I don't just want to generate the code for the sake of it (unless it is a simple helper script or a quick demo to test a few things before committing to a larger design), but I want to look at the code, learn from it, and become a better programmer myself.

That may just be me, I mean, I have talked now for so long about documentation that, realistically speaking, no reader should have reached this point. If you have actually made it this far into this article, this is your Everest, you have achieved something great today. You deserve pudding tonight!

So then, let's try to make sense then of the C++ code given by ChatGPT and learn from it. Maybe Joao, the developer of rapidyaml, can look through this as well and then finally add some examples people can use.

To create a new YAML file, we write:

1
2
3
c4::yml::Tree tree;
auto root = tree.rootref();
root |= c4::yml::MAP;

The first question you may have, what is c4 here? Well, let's look at the documentation again, shall we?

ryml ships with parts of c4core, a small C++ utilities multiplatform library, but you can ignore the in-source version of c4core and use a custom or system-installed c4core version.

That's all we are getting. So, naively you may then look c4core up, find the GitHub page, and, wow, the first thing you are seeing in the table of contents of the README.md file is Quicktour. Fantastic, this developer gets me, thank you. Within each heading within the quick tour, there are sections for, presumably, all the things you can do with c4core.

But then I noticed that there was no heading that would help me how to understand how to use the Tree class in c4core, so I scrolled through the Quicktour section and, oh no ... there is no documentation, but instead, a shopping list of items this library can do. But the real shock hit once I noticed that the Quicktour was just full of performance plots. Hang on, who develops this library?

It is Joao ...

OK, great, so let's throw whatever documentation c4core may or may not have overboard, at this point, I can't be bothered to spend any more time with Joao's self-glorifying pamphlet, I mean, even for me, as a German, his focus and dedication to speed and performance is a bit much. I'm looking at the source code now instead.

Good news: the first line, i.e. c4::yml::Tree tree; makes you believe that there should be a Tree class in c4core, right? After all, c4core re-implements the C++ standard library (which does not need reimplementation) and a Tree would be a logical container to add to all the other reimplementation (a tree just being a std::map (ordered/unordered) in C++, so also not needed). So you search the c4core library for a Tree and you find nothing.

So where is the Tree class located? Obviously within the rapidyaml library, because that is where a sane person would put it. OK, so after this near mental breakdown on my side and finding the tree, here are the lines that are of importance:

1
2
3
4
NodeRef Tree::rootref()
{
  return NodeRef(this, root_id());
}

To be fair, rootref is fairly self-explanatory, so seeing this Getter here isn't a huge surprise. We get back a reference (pointer) to NodeRef, which holds a node in the tree (a key in a key-value entry). Then we saw the use of the operator|=() being used. So, searching for that provided me with this:

1
void operator|= (type_bits t) { create(); m_tree->_add_flags(m_id, t); }

The create() function essentially creates a new node and then writes the value t that is given to the function to it. So, when we are writing root |= c4::yml::MAP;, where c4::yml::MAP is t in the operator|=() function, we essentially create a new node in the tree, and it has a new map associated with it. We can achieve the same in C++ using:

1
std::map<Type1, Type2> tree;

We just have to define the types here, which is the tricky bit (i.e. Type1 and Type2 determine the types of the keys and values). To make this a concrete example, let's say we want to store the boundary conditions, then we would write:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
#include <iostream>
#include <string>
#include <map>

int main() {
  std::map<std::string, std::string> bc;

  bc["north"] = "Dirichlet";
  bc["south"] = "Neumann";
  bc["east"]  = "Neumann";
  bc["west"]  = "Dirichlet";

  return 0;
}

This works as long as we don't need Type2 to change. YAML, however, requires more flexibility. For example, when we looked at the solver entry, we had:

1
2
3
4
5
6
7
solver:
  turbulence: 1
  model: komega-SST
  CFL:
    - 1e-06
    - 1
    - 10

So now we are storing a boolean value, a string, and some numbers. So you may say, OK, Joao has a point, and we need to define a special tree of some sort which allows us to store any type. Nope. We can implement the solver tree as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
#include <iostream>
#include <string>
#include <map>
#include <any>
#include <vector>

int main() {
  std::map<std::string, std::any> solver;

  solver["turbulence"] = true;
  solver["model"] = std::string("komega-SST");
  solver["CFL"] = std::vector<double>{1e-6, 1.0, 10.0};

  return 0;
}

Since the C++ standard of 2017, we have the std::any class that allows us to store any data type. Problem solved and we have a much cleaner syntax. While I haven't checked, I have to assume that c4core's Tree class, which is stored in rapidyaml, is more performant than a standard std::map<T1, T2> in C++, but, to me, there is not just a performance cost but also a maintenance cost (or technical debt). The more I get into rapidyaml, the more of a liability it seems to become!

OK, let's continue, thankfully the rest isn't that much more of a headache. To create a new node within the tree and then write elements to it, we saw, for example, with the mesh parameters, that we can achieve that with the following syntax:

1
2
3
4
auto meshing = root["meshing"];
meshing |= c4::yml::MAP;
meshing["numPointsX"] << 128;
meshing["numPointsY"] << 128;

First, we create a new node meshing in the root of our Tree. Then, we use the same, awkward operator|=() to indicate that we are storing a map here. This map has two entries, the numPointsX and numPointsY entries, both storing 128 as a number. We use the operator<<() now to write into this node, by now, I have given up and have accepted that this is the syntax.

There is precedent for it, the stringstream class, defined in the C++ standard library, also makes liberal use of the operator<<() to write into a stream, so while I am not a big fan of it, at least this is a defensible proposition (I sound like Claude Code now, one day I shall release TomGPT ...).

The rest of the YAML tree is built up in the same way, except for the CFL number, where we are storing an array. This is defined as:

1
2
3
4
5
auto cfl = solver["CFL"];
cfl |= c4::yml::SEQ;
cfl.append_child() << 1e-6;
cfl.append_child() << 1.0;
cfl.append_child() << 10.0;

We create the CFL node within the solver node, and assign it a type of c4::yml::SEQ, i.e. a sequence (or, as any other programmer, computer scientist, mathematician, physicist, engineer (yes, even civil engineers), anyone graduating high-school, would call it, an array). Sequences have an additional method append_child(), which mostly does what it says on the tin. I guess for consistency reasons (keeping a consistent bad style?!), we write:

1
cfl.append_child() << 1e-6;

and not:

1
cfl.append_child(1e-6);

That's what every other sane developer would likely do, just my 2 cents. If you try to do that, you get a compiler error, but, this code is perfectly fine:

1
cfl.append_child();

There is not even a compiler warning. So, if you were seeing this code, you may think this would be fine, the compiler would agree, and if you ran the following code:

1
2
3
4
5
6
7
8
auto cfl = solver["CFL"];
cfl |= c4::yml::SEQ;
cfl.append_child() << 1e-6;
cfl.append_child() << 1.0;
cfl.append_child() << 10.0;
cfl.append_child();
cfl.append_child();
cfl.append_child();

This is what you would get:

1
2
3
4
5
6
7
CFL:
  - 1e-06
  - 1
  - 10
  -
  -
  -

Beautiful, if it ain't a bug, it's a feature. So we can write empty arrays (excuse me, I mean ghost sequences, as they are probably referred to in c4core, I mean rapidyaml, I mean, I don't know any more what I mean ...). Not good.

At least writing a YAML tree to disk restores some sanity (mostly because Joao now trusts in the C++ standard library again, and he is not implementing his own file writing capability, let's not give him any wrong ideas ...). We had:

1
2
3
std::string yaml = c4::yml::emitrs_yaml<std::string>(tree);
std::ofstream out("config.yaml");
out << yaml;

First we convert the tree to a string, then we open a file, and we dump all of the tree (now as a string) into that file. Hallelujah, we managed to write a YAML file. Our next goal is to read the YAML file again.

At this point, I'd recommend taking a break, getting a tea, or doing something to reset your mind. Believe me, you'll need it. If you want to continue, that is fine, but do so at your own risk. It is not getting better. You have been warned.

I'll give you the entire code first, and then we discuss it below.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include <ryml.hpp>
#include <ryml_std.hpp>
#include <c4/yml/file.hpp>

#include <fstream>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>

int main() {
  // Read entire file
  // std::ifstream in("config.yaml");
  // std::string contents(
  //     (std::istreambuf_iterator<char>(in)),
  //     std::istreambuf_iterator<char>());

  // // rapidyaml parses in-place, so use a mutable buffer
  // std::vector<char> buffer(contents.begin(), contents.end());

  // c4::yml::Tree tree = c4::yml::parse_in_place(
  //     ryml::to_substr(buffer));
  std::string contents = c4::yml::file_get_contents<std::string>("config.yaml");
  ryml::Tree tree = ryml::parse_in_arena(ryml::to_csubstr(contents));
  auto root = tree.rootref();

  // meshing
  int numPointsX, numPointsY;
  root["meshing"]["numPointsX"] >> numPointsX;
  root["meshing"]["numPointsY"] >> numPointsY;

  // boundary conditions
  std::string north, south, east, west;
  root["boundaryConditions"]["north"] >> north;
  root["boundaryConditions"]["south"] >> south;
  root["boundaryConditions"]["east"]  >> east;
  root["boundaryConditions"]["west"]  >> west;

  // solver
  bool turbulence;
  root["solver"]["turbulence"] >> turbulence;
  
  std::string model;
  root["solver"]["model"] >> model;

  auto cfl = root["solver"]["CFL"];
  double cfl0, cfl1, cfl2;
  cfl[0] >> cfl0;
  cfl[1] >> cfl1;
  cfl[2] >> cfl2;

  std::cout << "numPointsX = " << numPointsX << '\n';
  std::cout << "numPointsY = " << numPointsY << '\n';

  std::cout << "north = " << north << '\n';
  std::cout << "south = " << south << '\n';
  std::cout << "east  = " << east  << '\n';
  std::cout << "west  = " << west  << '\n';

  std::cout << "turbulence = "
            << std::boolalpha << turbulence << '\n';

  std::cout << "model = " << model << '\n';

  std::cout << "CFL = ["
            << cfl0 << ", "
            << cfl1 << ", "
            << cfl2 << "]\n";

  return 0;
}

You may have noticed some code that is commented out in the beginning. Let's look at that in more detail:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
std::ifstream in("config.yaml");
std::string contents(
  (std::istreambuf_iterator<char>(in)),
  std::istreambuf_iterator<char>());

// rapidyaml parses in-place, so use a mutable buffer
std::vector<char> buffer(contents.begin(), contents.end());

c4::yml::Tree tree = c4::yml::parse_in_place(
  ryml::to_substr(buffer));

First line, OK, we are just creating an input stream (std::istream) and point it at the file which should be read into this stream. I am OK with that, not with the rest. I mean, lines 2-4 create a std::string and then use a stream iterator to go over the entire file to read it into the contents variable.

As a brief reminder, an iterator allows you to go over a range of items, think of an iterator as allowing you to go over an array of entries, for example, and it will know where the first item is, how to go to the next, and where the last item is. If you have ever seen a for loop (i.e. for (i = 0; i < size; ++i)), an iterator is essentially that, just as an object. And C++ is big on objects, and this whole object-oriented programming.

So, lines 2-4 read the content of the config.yaml file into contents by reading the file character by character, and the first streambuf_iterator is the starting point of the file, while the second denotes the end of a file. Does this make sense? To someone out there in the universe, I am sure, not to me, but it is one way of constructing a std::string (see std::basic_string's constructors, you want to look at the fourth type).

Once we have that string, we convert it to an array (std::vector) of characters on line 7 using again the start/end iterator concept we saw in lines 2-4. However, std::string (as any other container in the standard template library (STL) which has an iterator) exposes all of its starting and ending iterators through the begin() and end() methods.

You'll notice the comment on line 6: rapidyaml parses in-place, so use a mutable buffer. Translating from brain-dead AI slop (well, it is technically a very accurate computer-science terminology that is used here, but I would never call a computer scientist brain-dead) to English, mutable just means we can change things while immutable means things cannot change, and that usually means we either use const or constexpr.

We have this comment here because there are also immutable (const) iterators, and they are accessed through the cbegin() and cend() methods.

Now that we have read an entire file into a string using an istreambuf_iterator and then converted that into an array of characters, we can now parse all of these characters as content that represents a YAML file structure. We do that on lines 9-10. We create a Tree object, and then use c4core to parse in place.

What are we parsing in place? Well, we call the ryml::to_substr() method on our buffer (the array of characters). If we look at the implementation of that function, we have:

1
2
3
4
5
c4::substr to_substr(std::vector<char, Alloc> &vec)
{
  char *data = vec.empty() ? nullptr : vec.data(); // data() may or may not return a null pointer.
  return c4::substr(data, vec.size());
}

OK, we pass it a std::vector<char, Alloc> &vec, that checks out. We check if the vector actually contains data or if it is empty (line 3), then we simply return the data as a c4core substr. So, what is this substring object? Well, looking at the substring class implementation, it becomes apparent that we are now dealing with essentially a copy of C++'s std::string.

Yes, it is different, and I am sure it is fast, but we have duplicated efforts here. Why use std::string when you can reinvent the wheel? At this point, I wouldn't be surprised if I found c0core, c1core, c2core, and c3core as well as libraries, with c0core implementing its own operating system, c1core its own C++ compiler, c2core its own version of C, and then c3core adding object orientated capabilities a la C++.

I hope you can see that this is indeed the brain-dead version of implementing things. But this is what ChatGPT gave me. I refused. While I think Joao is capable of a lot of crimes against humanity, I don't believe this is one of them. I don't think he thought of reading files in this way even in his wildest fantasies. So I dug a bit deeper through the documentation to find a sensible way of opening files. This is such a basic operation, there must be a better way. And there is.

What I eventually stumbled upon was a combination of examples that allowed me to write this code:

1
2
3
std::string contents = c4::yml::file_get_contents<std::string>("config.yaml");
ryml::Tree tree = ryml::parse_in_arena(ryml::to_csubstr(contents));
auto root = tree.rootref();

3 lines of code. It does the same as the explicit version before, but it hides the complexity from the end-user. Joao still gets to work his magic around string-to-character-back-to-a-sort-of-but-not-C++-compatible-string conversion, but at least we don't have to look at the code. We can hide that complexity and maintain a cleaner code base. Abstraction at last.

The rest of reading the file content is then actually fairly straightforward (I mean, after this discussion, anything is ...). We first declare a variable and its type, and then point at the YAML file from where we want to read it. For example, we had:

1
2
3
4
5
bool turbulence;
root["solver"]["turbulence"] >> turbulence;

std::string model;
root["solver"]["model"] >> model;

Declare turbulence and model first, and then get it from root as shown. Job done. Arrays require one more step, we first have to get the node under which the array is stored, and then read individual items by indexing into the array that should be read as:

1
2
3
4
5
auto cfl = root["solver"]["CFL"];
double cfl0, cfl1, cfl2;
cfl[0] >> cfl0;
cfl[1] >> cfl1;
cfl[2] >> cfl2;

At this point, we have spent so much time with YAML that I am ready to move on, but I want to show you how you can work with anchors as well. Earlier, we looked at the following example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
case: &caseBase
  meshing:
    numPointsX: 128
    numPointsY: 128

  boundaryConditions:
    north: Dirichlet
    south: Neumann
    east: Neumann
    west: Dirichlet

  solver:
    turbulence: 1
    model: komega-SST
    CFL:
      - 1e-06
      - 1
      - 10

cases:
  coarse:
    <<: *caseBase
    meshing:
      numPointsX: 64
      numPointsY: 64
    
  medium:
    <<: *caseBase
    meshing:
      numPointsX: 128
      numPointsY: 128
  
  fine:
    <<: *caseBase
    meshing:
      numPointsX: 256
      numPointsY: 256

Let's say we want to read now the coarse, medium, and fine case for a grid independence study. We can do that with the following code:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include <ryml.hpp>
#include <ryml_std.hpp>
#include <c4/yml/file.hpp>

#include <fstream>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>

int main() {
  std::string contents = c4::yml::file_get_contents<std::string>("configWithAnchors.yaml");
  ryml::Tree tree = ryml::parse_in_arena(ryml::to_csubstr(contents));
  auto root = tree.rootref();

  // loop over all cases
  for (const char *currentCase : {"coarse", "medium", "fine"}) {
        
    // meshing
    int numPointsX, numPointsY;
    root["cases"][currentCase]["meshing"]["numPointsX"] >> numPointsX;
    root["cases"][currentCase]["meshing"]["numPointsY"] >> numPointsY;

    std::cout << "currentCase = " << currentCase << std::endl;
    std::cout << "numPointsX  = " << numPointsX << std::endl;
    std::cout << "numPointsY  = " << numPointsY << std::endl;
    std::cout << std::endl;
    
  }

  return 0;
}

Here, I am reading only the number of grid points, but we could have read any of the solver parameters here as well. First, we loop over an array of strings holding the entries "coarse", "medium", and "fine". We use that when we want to get a specific entry from our root variable. If you run this example, you will get:

``` currentCase = coarse numPointsX = 64 numPointsY = 64

currentCase = medium numPointsX = 128 numPointsY = 128

currentCase = fine numPointsX = 256 numPointsY = 256 ```

This is a pretty amazing functionality that YAML offers us. Creating nested parameter files has never been easier. I really love YAML for it, but not as much as the next file type we have to look at!

TOML

Call me egocentric, but Tom's Obvious Minimal Language, or TOML in short, somehow speaks to me as a fellow Tom. And no, I did not create it. JSON is often criticised for having too many brackets and commas. YAML takes those away, essentially, but still resembles JSON a lot. It gives us anchors/aliases on top, which is nice, but it still feels like a file generated by a computer for another computer. At least that is another common opinion I have heard.

That is the criticism TOML has taken as its inspiration and hence the Minimal Language identifier in the name; it is designed to be so simple, not just a computer, but even a human can read and understand it. As of the time of writing, the TOML homepage literally states:

A config file format for humans. TOML aims to be a minimal configuration file format that's easy to read due to obvious semantics.

I have been, for the longest time, a JSON and (to some degree) YAML supporter. I have used them in my previous codes. I always knew TOML existed, but I have never test-driven it.

Until one day, in summer 2025, a Friday night where I must have deviated from my very healthy diet of 1 milk-cream bread, 1 royal tea, and 1 bag of Korean crisps (the identifier Korean here is important, as half of their crisps are just cereals. Honestly, I have eaten most with milk. Call me disgusting. I am German. We are not known for our great taste!), and I must have run low on sugar and had the crazy idea of replacing my JSON config files by TOML.

Before you call me crazy, search for pong crush, these pongs are sold as cereals in the western world. In South Korea, you will find them in the crisp section. And there are more. For a country where it is almost impossible to find chocolate because "we don't eat sweets", it is surprising how much of what would be savoury/salty for most is sold as a sweet in Korea.

Ever tried garlic jam? Head to your nearest South Korean Lotte Mart and you will find it sold next to Nutella. Sure, they would say this is for making garlic bread, but if you buy that, it is also sweet, not savoury. Or how about a delicious onion tea? Order here. Why get just 1 bottle when 40 bottles are on offer?

And you thought I was the crazy one ...

Be that as it may, while it is a great tragedy I never got to consume gallons of onion tea, I did get to enjoy the benefits of using TOML over JSON. Once I tried that, I changed it in all of the CFD solvers I was working on at that time. First gridlock (private repo, I am not ready for you to see the C++ mess I cooked up there), and then swiftcfd, which is a very slow Python solver against the naming.

What's swift about it is its equation prototyping capabilities. I wrote it mainly for one of my ongoing research projects and gave it to one of my PhD students to help her progress beyond the simple toy example we came up with before. She is now using a code by Imperial College London. The irony ...

And, I can now proudly announce that TOML has my seal of approval. It is a great config file format and, spoiler alert, the one I would recommend anyone to use for their CFD solver for steering input files. Let's look at the rating first and then we jump into an example. The capabilities of the TOML were rated as:

Table 8
Format Comments Standard Inheritance Expressions & units Round-trip
TOML ★★★ ★★ ★ ★ ★

Doesn't look great, does it? Well, it does support comments, good, and while there is no standard, it does have a full specification that is published on its website. That is the minimum, someone wanting to implement it, can now go to a specific specification (all are versioned) and implement a parser against it. But, there is no steering committee, so no official mechanism to introduce new features, but, in general, most file types are so small and simple that they don't need one.

TOML is one of them, with its focus on being a minimal (markup) language, it probably doesn't need a committee behind it. One Tom making the decisions is enough. There is no inheritance, like in most other files, as well as no support for expressions and physical units. When you write the same file type back, you may end up with a slightly different file, so one star here as well for that.

While I have never found that to be an issue, as I mostly use TOML files as input files (i.e. reading only from them, never reproducing them), if you think about it in terms of a serialisation format, you probably don't want to use it. In any case, for the learning and integration, we had:

Table 9
Format Error messages Build burden Learning curve HPC / MPI scaling
TOML ★★★ ★★★ ★★★ ★★★

Error messages are crisp and precise. Have a look at the TOML formatter & validator, you can place your TOML file here, or even type it out live, and you will get detailed and precise error messages. While it does depend on the parser, the language itself is so simple and minimal that it is easy to implement these rule checking systems.

The build-burden is low as well. If you use TOML11 (which I would recommend, that is the library I end up using myself), you are dealing with a header-only library, so drop it into your existing project and continue coding, no need to add additional dependencies. Or, as usual, use Conan Center and you don't have any build burden anyway.

Regarding the learning curve, you will find that if TOML challenges you, you should pick up Marketing, not CFD. It is dead simple. No issues with scaling here again, just as JSON and YAML, it scales well across cores.

So, let us look at our solver input file again, now formatted in TOML syntax. We have:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
[meshing]
numPointsX = 128
numPointsY = 128

[boundaryConditions]
north = "Dirichlet"
south = "Neumann"
east = "Neumann"
west = "Dirichlet"

[solver]
turbulence = true
model = "komega-SST"
CFL = [1e-6, 1, 10]

We always start with a header enclosed in square brackets. The nice thing about the headers are they can be extended. For example, in JSON we introduced nested key-value pairs for the boundaries as:

1
2
3
4
5
6
7
8
{
  "boundaryConditions": {
    "north": {
      "type": "Dirichlet",
      "value": 0.25
    }
  }
}

In TOML, we achieve that by chaining headers together. The following is valid TOML and equivalent to the above-shown JSON file content:

1
2
3
[boundaryConditions.north]
type  = "Dirichlet"
value = 0.25

So, chaining these nested key-value pairs is as simple as joining them via dots in the header. What took 6 lines of input for JSON (not counting the first and last line), took 3 lines with TOML. No braces, quotation marks and commas needed. Raw information with, well, minimal markup. You see why I love the format so much?

OK, let's write a parser that can write and read our input file. I will use TOML11, as advertised previously. It is dead-simple to use, and you will find that writing and reading a TOML file is really straightforward. I hope our friend Joao will read this. There is a lot to be learned about clean code in this section! First up, writing a TOML file:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include <fstream>
#include <toml.hpp>

int main()
{
  toml::value config;

  // [meshing]
  config["meshing"] = toml::table{
      {"numPointsX", 128},
      {"numPointsY", 128}
  };

  // [boundaryConditions]
  config["boundaryConditions"] = toml::table{
    {"north", "Dirichlet"},
    {"south", "Neumann"},
    {"east",  "Neumann"},
    {"west",  "Dirichlet"}
  };

  // [solver]
  config["solver"] = toml::table{
    {"turbulence", true},
    {"model", "komega-SST"},
    {"CFL", toml::array{1e-6, 1.0, 10.0}}
  };

  std::ofstream ofs("config.toml");
  ofs << toml::format(config);

  return 0;
}

Single header include, then we define a TOML object on line 6, and then we simply write to it. TOML11 exposes the special type toml::table, which allows us to easily write a key-value pair under a specific header, as seen on lines 9, 15, and 23. You can also use STL containers directly (e.g. std::vector or std::map), and the documentation shows you how, but I find the toml::table syntax extremely expressive and so prefer that.

And yes, TOML11 comes with documentation, and one that is useful. To be honest, you can probably learn all you need from the README file on GitHub, it shows the most common use cases and then a bit more.

We define headers by simply adding them through the config["headername"] syntax. We then assign the content of the toml::table, which itself just takes key-value pairs as inputs. Lines 29-30 then write the entire TOML object to disk. Even if you have zero C++ knowledge, I bet you will still understand this code example, especially if we contrast that with what we had to go through with ryml.

Reading a file is just as simple:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <iostream>
#include <vector>
#include <toml.hpp>

int main()
{
  const auto config = toml::parse("config.toml");

  // meshing node
  const auto numPointsX = toml::find<int>(config, "meshing", "numPointsX");
  const auto numPointsY = toml::find<int>(config, "meshing", "numPointsY");

  // boundaryConditions node
  const auto north = toml::find<std::string>(config, "boundaryConditions", "north");
  const auto south = toml::find<std::string>(config, "boundaryConditions", "south");
  const auto east  = toml::find<std::string>(config, "boundaryConditions", "east");
  const auto west  = toml::find<std::string>(config, "boundaryConditions", "west");

  // solver node
  const auto turbulence = toml::find<bool>(config, "solver", "turbulence");
  const auto model = toml::find<std::string>(config, "solver", "model");
  const auto CFL = toml::find<std::vector<double>>(config, "solver", "CFL");

  // output (sanity check)
  std::cout << "numPointsX = " << numPointsX << std::endl;
  std::cout << "numPointsY = " << numPointsY << std::endl;

  std::cout << "north = " << north << std::endl;
  std::cout << "south = " << south << std::endl;
  std::cout << "east  = " << east  << std::endl;
  std::cout << "west  = " << west  << std::endl;

  std::cout << "turbulence = " << std::boolalpha << turbulence << std::endl;
  std::cout << "model = " << model << std::endl;

  std::cout << "CFL = ";
  for(const auto c : CFL)
      std::cout << c << ' ';
  std::cout << std::endl;

  return 0;
}

First, we have a one line solution for reading a file on line 7. This should be the norm, the fact that this feels like a victory tells us something about the unnecessary complexity of C++ (and how some people abuse it).

Finding elements in a TOML file is as easy as using toml::find<TYPE>(fileHandler, "header", "key"). You have to learn one function to read a TOML file. That's it. The function is templated on the type you are trying to read (i.e. the <TYPE> identifier), which will convert whatever variable you are reading to the type you want it to be.

The 3 arguments that we pass are the object that holds the TOML file information itself, the header where to locate it, and then the key we want to look up. It returns the value of that key-value pair.

It took me 3 paragraphs to explain how to read a TOML file in C++. I invite you to read JSON and YAML again and see if you find TOML easier. For me, it is not just easier, but a clean file, with a well-written parser library behind it just makes me happy, and it puts me in a good mood for the rest of the day. Today, I'll meet some of the new MSc students joining us at Cranfield, so they will probably be overwhelmed by my positivity, and they have no idea it is due to dealing with TOML files the entire morning!

What about nested key-value pairs? Well, if we wanted to get something like:

1
2
3
[boundaryConditions.north]
type = "Dirichlet"
value = 1.0

Then we can write it like:

1
2
3
4
config["boundaryConditions"]["north"] = toml::table{
  {"type", "Dirichlet"},
  {"value", 1.0}
};

If we want to read it we can do that with:

1
const auto north = toml::find<std::string>(config, "boundaryConditions", "north", "type");

We simply append more strings to query the TOML file tree until we find what we are looking for. This is probably the best example you can find if you ever wondered how variadic templates work in C++. They exist for exactly this type of scenario, i.e. we don't know how many function arguments we have, as our file tree we are parsing may have many nested key-value pairs.

Before we leave TOML files behind, I want to highlight one more functionality that is TOML11 specific. Besides the already very useful toml::find function, there is one more that is extremely useful if you are dealing with option parameters. If you look up a value, and it does not exist in your TOML file, then TOML11 will throw an exception.

Imagine your solver can do both fixed time stepping, and adaptive time stepping based on a fixed CFL number. You can now read your time step or CFL number from the input file, but, if you do fixed time stepping, you don't need CFL, if you do variable time stepping, you don't need a fixed time step to be read from your config file.

Now, you can implement logic in your file reading routine and then only read the information you need, or you can ask TOML11 to help you a bit by using toml::find_or<TYPE>(fileHandler, header, key, fallBack). The toml::find_or function looks for the header and key in your TOML file and will only read it if it can find it. If it can't, it will use the fallBack value you have provided.

Specifically, in our example, we can do something like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
const auto dt  = toml::find_or<double>(config, "time", "dt",  -1.0);
const auto cfl = toml::find_or<double>(config, "time", "cfl", -1.0);

if (dt > 0.0 && cfl < 0.0) {
  // fixed time-stepping
} else if (cfl > 0.0 && dt < 0.0) {
  // cfl-based time stepping
} else {
  // throw exception if logic not satisfied
}

There are better ways to deal with this, for sure (it will throw an exception if your TOML file contains both cfl and dt as an argument), but it just shows you how you may want to use fallback names. Another use case may be the writing of an output file. There is nothing worse than running your simulation to completion, and then having your solver crash at the very end when it wants to write the output file because no filename for the output file was specified. A simple fix would be to use:

1
const auto outputName  = toml::find_or<std::string>(config, "output", "filename", "out/results.vtk");

There really is a lot to like about TOML, I hope to have converted you to become a follower of Tom's Outstanding Minimal Language. I think that was the acronym ...

INI files

Ah, INI files, the files that dominated my childhood. Believe it or not, I didn't always want to become a lecturer in CFD when I was a child, no, first I wanted to become an architect (thankfully only for a very short time), and then a pilot. The pilot thing lasted for a good 10-15 years, I always thought I would end up flying a B747 for Lufthansa (even though I felt very patriotic and more like an Airbus guy, but there you go).

So, naturally, while most of my friends were playing football outside, in shorts, in winter (idiots), I was sitting in the comfort of my room playing the Microsoft Flight Simulator 95/98/2000/... I may have boasted about that before.

At one point, I circumnavigated the world in a B747 (see, I told you), starting in Frankfurt, with stops in Tokyo and San Francisco. Now, the party-poopers among you will point out that for a true circumnavigation, I need to cross the equator, to those, I say, why do you want to pick a fight with a minor (which I was back then)?

Anyway, the trip took me 37 hours, and I am ashamed to admit that I may have fallen asleep at the controls at some point, but nobody noticed, and all passengers got home safely. As I learned later, many of the people I studied with in my aeronautical engineering course (and, by many, I mean one other person), made their own circumnavigation. Though, they (he) did it with a Cessna 172. Respect!

Once you play the flight simulator for long enough, you start to mod the game. From installing third-party planes to making your own ones. Which I did. And, of course, each aircraft needs its own performance data to drive the simulation, and this is where the world of flight simulations meets aerodynamics, CFD, and INI files.

The aircraft performance data was stored in config files using a .cfg extension, but these files were, essentially, just plain INI files. The [latest] flight simulator appears to have the same file type, with some additions bolted on top of it. So, naturally, I spent a lot of time checking how these parameters were influencing the simulation.

It was also the time I got more interested in operating systems and discovered the console/DOS, so naturally that is where I could be found when I was not playing, I mean, simulating flying. And yes, I started changing system INI files as well, wanting to see how that would change the performance of my PC. Most of the time, it required a fresh installation of Windows and formatting of the hard-drive. We learn from our mistakes, right?

OK, so if you have never worked with INI files, but you have read the previous section on TOML files, you'd be surprised how similar these are. But, there are small differences, and these are important for us in the context of CFD.

An INI file simply stores key-value pairs like TOML, and it also uses headers in the same way as TOML, but it does not support nested headers/key-value pairs. TOML has native support for types, so we can differentiate between a string and a float, for example. In INI files, everything is a string (raw text). There is also no support for arrays in INI files.

In a nutshell, if TOML were a child, then INI would be his toddler brother. Both humans, but very different skill sets. Now that we have gone through JSON, YAML, and TOML, it might be surprising to learn that most CFD solvers out there use INI files for their setup. Some accept that and state so plainly (PyFR), others are oblivious to the fact and believe they have invented their own syntax (SU2).

My very first functional CFD solver used INI files as well, but I had no clue I was using INI files. I re-invented the file format without knowing it. Probably, for the same reason, you will find many other CFD solvers using this format. Like me, they most likely just don't know that they have better alternatives available. If you know such a unicorn, direct them towards this article, and they will be forever grateful to you.  

An INI file for the same solver configuration we used before can be written out as:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
[meshing]
numpointsx = 128
numpointsy = 128

[boundaryconditions]
north = Dirichlet
south = Neumann
east = Neumann
west = Dirichlet

[solver]
turbulence = true
model = komega-SST
cfl = 1e-6, 1, 10

We see where TOML stole the header ideas from. Not all INI files may have headers. SU2, for example, does not use headers, instead, they prepend headers to variables. So, the above config file in SU2 syntax would look like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
meshing_numpointsx = 128
meshing_numpointsy = 128

boundaryconditions_north = Dirichlet
boundaryconditions_south = Neumann
boundaryconditions_east = Neumann
boundaryconditions_west = Dirichlet

solver_turbulence = true
solver_model = komega-SST
solver_cfl = 1e-6, 1, 10

One way of doing it, but why would you? In any case, we notice the absence of any quotation marks, commas, and the like, and also, that there are no capital letters. No arrays, (OK, we have commas here, but that is just to separate values and to simulate an array), and, if you want to be on the safe side, I would also not recommend any whitespaces. It would be advisable to have the config file instead as:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
[meshing]
numpointsx=128
numpointsy=128
[boundaryconditions]
north=Dirichlet
south=Neumann
east=Neumann
west=Dirichlet
[solver]
turbulence=true
model=komega-SST
cfl=1e-6,1,10

Beautiful, isn't it? It has this Windows 95 retro look to it, at least to me. Before I start looking for the perfect synth music that would make the nostalgia complete (and, knowing myself, I would end up on YouTube looking for the perfect tune and lose about 2 hours doing so), let's look at our functionality table. We had:

Table 10
Format Comments Standard Inheritance Expressions & units Round-trip
INI / key–value ★★★ ★ ★ ★ ★★

While not shown in our example, we can use comments with INI files. Comments start with a semicolon (;) and any text after it will be ignored. Most of the time anyway, sometimes comments start with a hash sign (#). It depends on the day of the week, who you ask, and, more importantly, on your operating system and the parser you use.

As you can see, there is no standard to speak of. INI files are mostly driven by convention. Microsoft may have popularised them, but I guess no one ever thought of INI files as their own dedicated file format, and so no one ever thought of standardising the design. The standards are now embodied in the parser you are using to read and produce INI files. It is the Markdown situation all over again!

It's not great, and you can see that with SU2 (again, sorry for picking on you, I actually quite like SU2 as a solver!). If you throw their input files through an INI validator, you'll see it is non-compliant. So, you could now argue that SU2 either uses non-compliant INI files, a superset of INI with their own extensions, or simply their own bespoke file format that just happens to overlap 95% with INI files. Your call.

The bottom line is, if you still want to use INI files in your solver, pick a standard (or parser) and stick with it. But, be ready that people on the internet with a website on CFD will look at your solver and scrutinise your use of INI files.

It will, then, come as no surprise that we have no inheritance as in YAML, or any support for expressions and units. For round-tripping, I give it 2 stars because we are dealing with strings, we can't really break anything. But, if we start converting strings into numbers and then write those numbers out again, we do get different results.

In terms of learning curve, we had:

Table 11
Format Error messages Build burden Learning curve HPC / MPI scaling
INI / key–value ★★ ★★★ ★★★ ★★★

In terms of error messages, I gave it a tentative 2-star rating. That is mostly my ignorance. I have not checked and used all INI file parsers out there. I am sure error handling is a thing for some, but not for others. So you may, or may not, get lucky with yours. But then again, INI files are pretty forgiving. Whatever you do to them, they mostly work fine, even if you put rubbish into them. For example, take the following INI file:

1
2
3
4
5
6
; this is a comment
# this may be a comment, too
solver=SIMPLE
scheme========upwind
!rans!
what is the purpose of this file format?

Using the INI validator mentioned above, this is a valid INI file. It is nonsensical, but it does produce the following converted JSON file:

1
2
3
4
5
6
{
  "solver": "SIMPLE",
  "scheme": "=======upwind",
  "!rans!": true,
  "what is the purpose of this file format?": true
}

So, yeah, a parser may tell you that your INI file is stupid, but it can still generate a valid output. In terms of build-burden, there is none, really. There are plenty of header-only libraries out there, one of which we will look at in just a second, and I'd also argue that there is no learning curve.

In fact, I should have given it 4 stars out of 3 stars. Why? Because it is so simple, even software engineers at Microsoft can comprehend it, and as we know, this is the gold standard.

For HPC and MPI-scaling, I give it 3 stars. It's text, it's simple, what can go wrong?

So, let's look at libraries that will help us to read and write INI files. I have settled on mINI. Why? It is header-only, and the README file does give me all the information I need to start using it. I think ryml gave me PTSD. It's good to work with a sensible library like this one for a change.  

It is available on Conan Center, so installation is as difficult as running one Conan command. After that, we can include it in our project with the #include "mini/ini.h". To write an INI file using mINI, this is how we can achieve it:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#include <iostream>
#include <sstream>
#include <vector>
#include <string>

#include "mini/ini.h"

int main()
{
  mINI::INIStructure config;

  // [meshing]
  config["meshing"]["numPointsX"] = "128";
  config["meshing"]["numPointsY"] = "128";

  // [boundaryConditions]
  config["boundaryConditions"]["north"] = "Dirichlet";
  config["boundaryConditions"]["south"] = "Neumann";
  config["boundaryConditions"]["east"]  = "Neumann";
  config["boundaryConditions"]["west"]  = "Dirichlet";

  // [solver]
  config["solver"]["turbulence"] = "true";
  config["solver"]["model"]      = "komega-SST";
  config["solver"]["CFL"]        = "1e-6,1,10";

  mINI::INIFile outFile("config.ini");
  outFile.generate(config);

  return 0;
}

Fantastic! We declare a mINI object as mINI::INIStructure config;, and then start dumping entries to it. Notice that I write strings only here, which is what INI files expect. Once all is inside the config object, we save it to a file called config.ini. As simple as that!

In fact, if you look at the source code of mINI, it only is about 800 lines of code! INI files are simple, and simple files don't require large libraries to make them work.

Reading this config.ini file now is just as simple:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <iostream>
#include <sstream>
#include <vector>
#include <string>

#include "mini/ini.h"

int main()
{
  mINI::INIStructure config;
  mINI::INIFile inFile("config.ini");
  inFile.read(config);

  int numPointsX = std::stoi(config["meshing"]["numPointsX"]);
  int numPointsY = std::stoi(config["meshing"]["numPointsY"]);

  std::string north = config["boundaryConditions"]["north"];
  std::string south = config["boundaryConditions"]["south"];
  std::string east  = config["boundaryConditions"]["east"];
  std::string west  = config["boundaryConditions"]["west"];

  bool turbulence = config["solver"]["turbulence"] == "true";

  std::string model = config["solver"]["model"];

  // Parse CFL array
  std::vector<double> cfl;
  {
    std::stringstream ss(config["solver"]["CFL"]);
    std::string argument;
    while (std::getline(ss, argument, ','))
    {
        cfl.push_back(std::stod(argument));
    }
  }

  // Print values
  std::cout << "numPointsX = " << numPointsX << '\n';
  std::cout << "numPointsY = " << numPointsY << '\n';

  std::cout << "north = " << north << '\n';
  std::cout << "south = " << south << '\n';
  std::cout << "east  = " << east  << '\n';
  std::cout << "west  = " << west  << '\n';

  std::cout << "turbulence = " << std::boolalpha
            << turbulence << '\n';

  std::cout << "model = " << model << '\n';

  std::cout << "CFL = ";
  for (double v : cfl)
    std::cout << v << " ";
  std::cout << '\n';

  return 0;
}

You'll notice a few things here. First, we can read an INI file again using built-in functionality:

1
2
3
mINI::INIStructure config;
mINI::INIFile inFile("config.ini");
inFile.read(config);

That's what I like to see. Since we store all of our information as text, we need to convert strings to the right data type now. For example, for the number of meshing nodes, we had:

1
int numPointsX = std::stoi(config["meshing"]["numPointsX"]);

We use std::stoi (string to int) to convert, well, a string to an int. For the boolean evaluation we used a ternary operator in C++ as:

1
bool turbulence = config["solver"]["turbulence"] == "true";

This is the shorthand version, but we could have written it as well in its more verbose form as:

1
bool turbulence = config["solver"]["turbulence"] == "true" ? true : false;

So, if config["solver"]["turbulence"] is equal to "true", then this expression will evaluate to true and that will be the value for our turbulence variable. If not, then we set it to false.

To convert the CFL number array into a std::vector, we first read the line using a stringstream (which is, unfortunately, usually abbreviated as ss). I am going with this convention here so I can point out that you may want to use a different name. It honestly gives me the creeps! We had:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
std::vector<double> cfl;
{
  std::stringstream ss(config["solver"]["CFL"]);
  std::string argument;

  while (std::getline(ss, argument, ','))
  {
    cfl.push_back(std::stod(argument));
  }
}

Stringstreams provide some favourable properties over a classical string, when it comes to manipulating them (that includes conversion). So, we initialise our stringstream as std::stringstream ss(config["solver"]["CFL"]);. We can then use std::getline, which will take a stringstream as an argument, a destination to write the content to (here the argument), and the delimiter by which we want to split the stringstream, here ,.

We then get each argument in our CFL array, and we simply have to convert that string to a double using std::stod. Furthermore, we could have also converted to a float with std::stof, but I would only recommend doing that if you are on pre-millennial hardware. In our day and age, we should be comfortable with double as the default floating-point representation.

The file then just prints the content it read from file so we can make sure we got everything correctly. And yes, I have tested it, it works like a charm. So, if you really want to use INI files, go ahead. The reason I have placed INI files after JSON, YAML, and TOML is so that you see that there are other, better, alternatives out there. I would encourage you to use them instead.

It's only going to get worse. Let's see, what is next? XML! Oh, no ...

PS: I did end up on YouTube again to look for that perfect synth music and I stumbled upon this infomercial?! I cannot stop watching it. I leave you with it, if you have any affinity towards synth music, this is a strange one ...

XML files

When you get to a certain age, it is recommended, as a male, to get a prostate exam. While I don't want to get into the nitty-gritty details of it, I had mine last year in December (an early Christmas present if you like ...), but I got it while living in South Korea, making good use of their excellent health care provision. South Korean prostate exams are performed while you are under full anaesthesia, so you don't actually notice any of the procedure.

The weird part about the "sleeping prostate exam", as they sell it, is that at some point you wake up and the nurse will help you walk out. They keep you in a dark room to wake up, and, as you leave, you walk past rows of beds with others waiting to wake up. I think this was the closest I, hopefully, ever got to organ farming. The setting was surreal, straight out of Resident Evil, and I think I still have all my organs. Or so at least I hope ...

It does sound like a nice, painless procedure, but it means that they ram a needle the size of a screwdriver into the back of your hand to administer the anaesthetic. I'd rather have two of those fitted to both of my hands now than having to write about XML. Honestly, the person who invented XML deserves to be shot!

Before writing this article, I had no idea who invented XML. So, I looked it up and found this golden nugget of information. I read it and then thought to myself who has a website and writes about file types, how pathetic! ... It took a second to realise that I am one of those people, apparently. But it is OK, I mean, we talk about file types. Realistically speaking, no one will ever reach this far into the article, so I am mostly entertaining myself now knowing that no one will mind.

The website writes how Goldfarb, one of the persons involved in getting XML off the ground, stated:

[XML is] the holy grail of computing, solving the problem of universal data interchange between dissimilar systems.

To his credit, he, together with his colleagues Mosher and Lorie, invented a first version of what would later become XML. Fearing being forgotten to history, they did what any person with a half-baked solution would do: name it after themselves, so people would have to mention them, at least in a footnote, when writing about XML. Mission accomplished.

So, Goldfarb, Mosher, and Lorie came up with, drum roll: GML. Brilliant! Later, apparently Goldfarb realised that GML could be so much more and the mistake of naming it after themselves. So, he took to ChatGPT (in the 1970s, or whatever equivalent method was around back then (probably Cleverbot (it seems to be running an LLM now as well, it used to be hilarious))) and rebranded to General Markup Language (still GML), or Standard General Markup Language (SGML).

As we all know if you include the name "Standard" in your name, people will automatically take you more serious (especially me). I, personally, would have loved to see Goldfarb rebranding to Goldfarb's Markup Language, and, this is how I will remember it because it does suggest a large falling out between the once 3 friends who now have decided to go their own way. Perhaps that is still the case, who knows, I just like to spice up a long forgotten drama of the 1970s that may not even have taken place (it probably hasn't).

From that standard, XML would later emerge (the eXtensible Markup Language). While I do think GML really is just a footnote in history nowadays, it is one worth mentioning as it popularised the word Markup in the field of computer science/software engineering. In fact, Goldfarb is widely seen as the inventor of Markup languages, not because he was first, but because he shouted loudest and people listened (that is my interpretation anyway of the Wikipedia summary).

Right, I cannot put it off any longer, it is time to face the demon that is XML. And, for context, I should add, XML is probably a useful thing in certain applications. I cannot deny that. But, what I am saying is that there is a right place for the right file type, and XML+CFD is as good as Garlic+Jam. Sure, you can mix it, but it is an explosive mixture!

So, why bother? If I hate it so much, why exacerbate my PTSD and go through all of the XML horror if I don't believe it has a place in CFD? Well, because that is just my opinion. Clearly others see it differently and believe XML+CFD is the way to go (code_saturne, for example, but they are French, so it is on-brand!).

Because you will have the (dis)pleasure of encountering XML in the wild, I'd be doing you a disservice by letting my opinions steer my articles (I'd never do that of course ...). Returning to our original solver input file, in XML, we would write it down like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
<?xml version="1.0" encoding="UTF-8" ?>
<configuration>
  <meshing>
    <numPointsX>128</numPointsX>
    <numPointsY>128</numPointsY>
  </meshing>
  <boundaryConditions>
    <north>Dirichlet</north>
    <south>Neumann</south>
    <east>Neumann</east>
    <west>Dirichlet</west>
  </boundaryConditions>
  <solver>
    <turbulence>true</turbulence>
    <model>komega-SST</model>
    <CFL>
      <Item1>0.000001</Item1>
      <Item2>1</Item2>
      <Item3>10</Item3>
    </CFL>
  </solver>
</configuration>

We also had nested key-value pair entries, and those can just as easily be represented by XML:

1
2
3
4
5
6
<boundaryConditions>
  <north>
    <type>Dirichlet</type>
    <value>0.25</value>
  </north>
</boundaryConditions>

Place your header or key in <> brackets, and then close them with </> again. If you want to store a value, place it between <NAME> and </NAME>. Conceptually it is just as easy as learning an INI file, while it allows for an arbitrarily deep nesting of key-value pairs. But, everything is a string again. We can see how we can store arbitrarily complex information with this simple structure, but for my taste, we can do exactly the same with JSON, or TOML, where we have a nice separation of headers and key-value pairs.

No support for arrays, so we need to place array-like information into a nested key-value pair entry and then because we have to give each entry a name, use something like Item1, Item2, Item3, etc. It works, but it is too much typing to write out an array that can be stored as [1e-6,1,10] in most other formats.

Of course, we could have also separated each entry by a comma and then produced only a single entry, like we did with INI, but XML favours having each item wrapped around in its own entry. Let's obey the wishes of the XML gods (I am sure God farb, I mean Goldfarb, sees himself as one of them, as of the writing of today, he is still alive).

XML has too much typing and too much information. It is alright for the simple configuration entry we are looking at right now, but we should never generalise from a simple toy example. If you have ever looked at longer XML files, you'll know what I mean. In fact, in our HDF5 article, we looked at how we can use XML to read HDF5 files in Paraview. Let's just review it one more time, and you will see what I mean:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
<?xml version="1.0" ?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="3.0">
  <Domain>
    <Grid Name="TaylorGreenVortex" GridType="Uniform">
      <!-- 3D rectilinear mesh: 32x32x32 nodes -->
      <!-- Dimensions order matches HDF5 storage: slowest-varying index first (x here) -->
      <Topology TopologyType="3DRectMesh" Dimensions="32 32 32"/>
      <!-- Separate 1D coordinate vectors for each axis -->
      <Geometry GeometryType="VXVYVZ">
        <DataItem Name="x" Dimensions="32" NumberType="Float" Precision="8" Format="HDF">
          tgv-highfive.h5:/mesh/x
        </DataItem>
        <DataItem Name="y" Dimensions="32" NumberType="Float" Precision="8" Format="HDF">
          tgv-highfive.h5:/mesh/y
        </DataItem>
        <DataItem Name="z" Dimensions="32" NumberType="Float" Precision="8" Format="HDF">
          tgv-highfive.h5:/mesh/z
        </DataItem>
      </Geometry>
      <!-- Velocity component u -->
      <Attribute Name="u" AttributeType="Scalar" Center="Node">
        <DataItem Dimensions="32 32 32" NumberType="Float" Precision="8" Format="HDF">
          tgv-highfive.h5:/solution/u
        </DataItem>
      </Attribute>
      <!-- Velocity component v -->
      <Attribute Name="v" AttributeType="Scalar" Center="Node">
        <DataItem Dimensions="32 32 32" NumberType="Float" Precision="8" Format="HDF">
          tgv-highfive.h5:/solution/v
        </DataItem>
      </Attribute>
      <!-- Velocity component w -->
      <Attribute Name="w" AttributeType="Scalar" Center="Node">
        <DataItem Dimensions="32 32 32" NumberType="Float" Precision="8" Format="HDF">
          tgv-highfive.h5:/solution/w
        </DataItem>
      </Attribute>
      <!-- Pressure p -->
      <Attribute Name="p" AttributeType="Scalar" Center="Node">
        <DataItem Dimensions="32 32 32" NumberType="Float" Precision="8" Format="HDF">
          tgv-highfive.h5:/solution/p
        </DataItem>
      </Attribute>
    </Grid>
  </Domain>
</Xdmf>

If you look at how much setup is required for each entry and compare that against the data we are actually trying to store, you see that it is quickly getting out of hand. Let's take one specific example to make that obvious. Let's take the last attribute (pressure) we are storing. In XML, we have:

1
2
3
4
5
<Attribute Name="p" AttributeType="Scalar" Center="Node">
  <DataItem Dimensions="32 32 32" NumberType="Float" Precision="8" Format="HDF">
    tgv-highfive.h5:/solution/p
  </DataItem>
</Attribute>

If we wanted to store all of the same information, but, say, with TOML, then we'd have:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
[p]
AttributeType = "Scalar"
Center        = "Node"

[p.DataItem]
Dimension     = [32, 32, 32]
NumberType    = "Float"
Precision     = 8
Format        = "HDF5"
source        = "tgv-highfive.h5:/solution/p"

Sure, it has more lines, and I have also taken some liberty here to introduce some whitespaces to align the values, I just personally prefer data being laid out this way than the dense XML way where we repeat each header and key. Well, as I said, that is my opinion, you are entitled to your own.

Trying to be objective now for a bit, this is the rating I came up with for its features:

Table 12
Format Comments Standard Inheritance Expressions & units Round-trip
XML ★★★ ★★★ ★★★ ★ ★★★

We can use comments as we have seen in the HDF5 example, there is a strong and mature standard behind it, and it does support some form of inheritance. While I am not a big fan of it, as it is trying to mash up object orientated programming with markup language, it is mature and well tested. For that, we need XML Schema Definition, or XSD. It is essentially an extension and this is where we can really feast on the ugliness of XML.

I couldn't do a better service than Anders Møller to explain XSD, mainly because the website proudly states that it has not been updated since 2003! I do like a good old HTML programmed website by academics from the 2000s, the only thing missing is 8-bit background music! So support is there, it isn't pretty, but it deserves 3 stars for it.

There isn't any way to express units or expressions, unless you add your own units to a quantity. For example, you could have:

1
2
3
4
5
<boundaryConditions>
  <north>
    <value units="m/s">1.0</value>
  </north>
</boundaryConditions>

This can be done, XML allows for it, but no XML parser will look for your unit argument. In the end, you would have to write your own XML parser, and that is true for all file types we looked at previously. We can always extend it ourselves, but there is no native support for this.

XML does offer one advantage over other file types and that is it can be stored in binary. As a result, it is trivial to write out data and read it back, and get bit-identical results. Round-tripping deserves 3 stars here for that reason.

Overall, it is a very capable file format. London is, I am sure, a beautiful city to some. I still wouldn't want to live in it and neither do I want to use XML in my life anywhere.

In terms of development burden, we had:

Table 13
Format Error messages Build burden Learning curve HPC / MPI scaling
XML ★★ ★★★ ★★ ★★

Parsers do provide sensible error messages most of the time, but XML is an old, mature standard that allows you to do mostly anything you want so it is a bit difficult to foresee every possible failure scenario and warn against it. Because it is so old and mature, many libraries exist that will help you read and write XML files. No need to do that yourself. XML does scale well on clusters if we use binary formats, but it could be better, so 2 stars only here.

After having learned my lesson with rapid YAML, I didn't jump onto the first library I could get my hands on. I did read the documentation a bit more carefully this time of a few libraries and I have settled on pugixml. It has its own website with a quick-start guide (always welcome!), it can be installed through Conan Center, the only downside, if anything, is that it is compiled, not header-only.

I have come to appreciate header-only libraries since working with Eigen extensively. It does slow down your compilation, but it does optimise your performance potentially more than using compiled libraries injected into your executable. It also removes the hassle of troubleshooting issues when you compile the library yourself and it does not work.

Pugixml compiled fine on my end, so far so good, so let's look at how we can write our solver configuration file with it in C++:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#include <pugixml.hpp>

int main()
{
  pugi::xml_document doc;

  // XML declaration
  pugi::xml_node header = doc.prepend_child(pugi::node_declaration);
  header.append_attribute("version") = "1.0";
  header.append_attribute("encoding") = "UTF-8";

  // Root element
  pugi::xml_node root = doc.append_child("configuration");

  // meshing
  pugi::xml_node meshing = root.append_child("meshing");

  meshing.append_child("numPointsX").text() = 128;
  meshing.append_child("numPointsY").text() = 128;

  // boundaryConditions
  pugi::xml_node bc = root.append_child("boundaryConditions");

  bc.append_child("north").text() = "Dirichlet";
  bc.append_child("south").text() = "Neumann";
  bc.append_child("east").text()  = "Neumann";
  bc.append_child("west").text()  = "Dirichlet";

  // solver
  pugi::xml_node solver = root.append_child("solver");

  solver.append_child("turbulence").text() = true;
  solver.append_child("model").text() = "komega-SST";

  pugi::xml_node cfl = solver.append_child("CFL");
  cfl.append_child("Item1").text() = 0.000001;
  cfl.append_child("Item2").text() = 1.0;
  cfl.append_child("Item3").text() = 10.0;

  // Save pretty-printed XML
  doc.save_file("config.xml");

  return 0;
}

Even though I am not the biggest fan of XML, at least the library to write it provides a clean interface. And you know what, if I never had to actually look at the darn XML file, I would be OK with using pugixml to read and write files, I do appreciate a clean and well-designed API.

Let's look at a few things here. First, we declared the XML document object as pugi::xml_document doc;. We start dumping data to it, and then we use it to save whatever state it contains. Then, when we open a file and write the XML header, the part that any XML file will need and have, we had:

1
2
3
pugi::xml_node header = doc.prepend_child(pugi::node_declaration);
header.append_attribute("version") = "1.0";
header.append_attribute("encoding") = "UTF-8";

This produces:

1
<?xml version="1.0" encoding="UTF-8" ?>

We use prepend_child() here to ensure that it is the first element in the XML file later. Also, we pass it the pugi::node_declaration argument, which is a special identifier to indicate we want to write the mandatory header information. And then, we also see the append_attribute() calls, which we perform on the variable header. An attribute is anything that appears inside the header, we see that we set version=1.0 and encoding=UTF-8.

It is the only place in the document where we add attributes. But let's look at it a bit further and see how we could exploit that. To write the meshing node, we had:

1
2
3
4
5
6
7
8
// Root element
pugi::xml_node root = doc.append_child("configuration");

// meshing
pugi::xml_node meshing = root.append_child("meshing");

meshing.append_child("numPointsX").text() = 128;
meshing.append_child("numPointsY").text() = 128;

First, we create the mandatory root node configuration that contains everything. This is our root. Then, we append a child node to the root, which is our meshing node, by calling root.append_child("meshing"). A reference to this newly created root is stored in the variable meshing, and now we can append to it. We do that when we set the number of points in X and Y as seen on lines 7-8.

Notice that we write the information here as text. XML only understands text. So, we have the same conversion issue as INI files again, and will need to later convert to an appropriate type. But, I wanted to show you here that we could also append attributes here as well. For example, we could also write:

1
2
3
pugi::xml_node numPointsX = meshing.append_child("numPointsX");
numPointsX.text() = 128;
numPointsX.append_attribute("unit") = "m";

This would produce:

1
2
3
<meshing>
  <numPointsX unit="m">128</numPointsX>
</meshing>

So, we could add information such as units to a quantity, XML provides us with a mechanism to store any type of attribute. If we now coupled that with our solver and used a library like mp-units, then we could have decent unit support. But it would not be a property of the file format, instead, it would be a property of our solver to convert any type of input into the correct units. So, XML provides us with some support, but it does not have in-built units.

The rest of the file follows a very similar pattern. So, instead of looking at each entry again, let's do something different here because, on my machine, pugixml finally gave me an example of what issues you may expect when we are trying to write a file and then try to read it back again (round-tripping). This is the file that is produced by the code above:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
<?xml version="1.0" encoding="UTF-8"?>
<configuration>
  <meshing>
    <numPointsX>128</numPointsX>
    <numPointsY>128</numPointsY>
  </meshing>
  <boundaryConditions>
    <north>Dirichlet</north>
    <south>Neumann</south>
    <east>Neumann</east>
    <west>Dirichlet</west>
  </boundaryConditions>
  <solver>
    <turbulence>true</turbulence>
    <model>komega-SST</model>
    <CFL>
      <Item1>9.9999999999999995e-07</Item1>
      <Item2>1</Item2>
      <Item3>10</Item3>
    </CFL>
  </solver>
</configuration>

Did you spot it? When we write the CFL number, we get:

1
<Item1>9.9999999999999995e-07</Item1>

However, we were writing to our XML file as:

1
2
pugi::xml_node cfl = solver.append_child("CFL");
cfl.append_child("Item1").text() = 0.000001;

This is an example of where we are not getting the same output of what we were trying to write, and this is a particular issue we face with floating-point numbers. We get some curious behaviour when we try to read it back, though. Let's look at that as well:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <iostream>
#include <vector>
#include <pugixml.hpp>

int main()
{
  // read XML file
  pugi::xml_document doc;
  pugi::xml_parse_result result = doc.load_file("config.xml");

  pugi::xml_node root = doc.child("configuration");

  // meshing
  pugi::xml_node meshing = root.child("meshing");

  int nx = meshing.child("numPointsX").text().as_int();
  int ny = meshing.child("numPointsY").text().as_int();

  // boundary conditions
  pugi::xml_node bc = root.child("boundaryConditions");

  std::string north = bc.child("north").text().as_string();
  std::string south = bc.child("south").text().as_string();
  std::string east  = bc.child("east").text().as_string();
  std::string west  = bc.child("west").text().as_string();

  // solver
  pugi::xml_node solver = root.child("solver");
  bool turbulence = solver.child("turbulence").text().as_bool();
  std::string model = solver.child("model").text().as_string();

  std::vector<double> cfl;
  pugi::xml_node cflNode = solver.child("CFL");
  cfl.push_back(cflNode.child("Item1").text().as_double());
  cfl.push_back(cflNode.child("Item2").text().as_double());
  cfl.push_back(cflNode.child("Item3").text().as_double());

  std::cout << "cfl[0]: " << cfl[0] << std::endl;
  std::cout << "cfl[1]: " << cfl[1] << std::endl;
  std::cout << "cfl[2]: " << cfl[2] << std::endl;

  return 0;
}

Loading an XML file is sensibly done with a single line of code. Accessing data is also done very similarly to how we were writing data before. Instead of using append_child() to add child elements to the XML object, we use child() to query the XML object to find a specific child. We still call the text() function on that child since we want to read that property (as opposed to an attribute, for example). Pugixml also provides us with convenient converters that we can chain into the call.

So, let's look at XML writing and reading side by side, using the meshing node again as an example:

1
2
3
4
5
6
7
// write number of Points in X under the meshing node
pugi::xml_node meshing = root.append_child("meshing");
meshing.append_child("numPointsX").text() = 128;

// read number of Points in X under the meshing node
pugi::xml_node meshing = root.child("meshing");
int nx = meshing.child("numPointsX").text().as_int();

This is likely one of the cleanest APIs I have seen, and just because of this API I really want to like XML enough to use it but still, I feel allergic to it. Such a shame!

The rest of the document follows again a similar convention, and we see there are converters for int, double, bool, and strings, which allows us to read the text in any type we want. Beautiful!

In any case, you see that I am also printing the CFL number entries at the end. Want to take a guess what numbers I am getting back? Here it is:

cfl[0]: 1e-06 cfl[1]: 1 cfl[2]: 10

This may seem surprising! Initially, we wanted to write cfl[0]=1e-6, but saw that we actually got cfl[0]=9.9999999999999995e-07. Now after reading this value back into XML we get again cfl[0]=1e-6. There are two things at play here we need to understand to make sense of the situation:

  • The IEEE-754 standard on how to store floating-point numbers (the standard essentially everyone has adopted for floating-point numbers)
  • How numbers are displayed in C++

The IEEE-754 floating-point standard sets out how to store floating-point numbers. We can't store them natively on our hardware, as everything is based on binary numbers, and we can only store bit information, i.e. 0 or 1 (true or false, on or off). The standard says that we are storing floating-point numbers in this format:

−1sign⋅(1+fraction)⋅2exponent−bias(4) -1^{sign}\cdot (1+fraction)\cdot 2^{exponent - bias} \tag{4}

The actual things we store are the sign, the fraction, and the exponent. These are integers and can be represented in binary numbers easily. For example, if we had sign=1sign = 1, exponent=1025exponent=1025, bias=1023bias = 1023 (this is always the case, defined in the standard), and fraction=0012=0/21+0/22+1/23=0/2+0/4+1/8fraction=001_2=0/2^1+0/2^2+1/2^3=0/2+0/4+1/8, we get:

−11⋅(1+02+04+18)⋅21025−1023−11⋅(1.125)⋅22−1⋅1.125⋅4−4.5(5) -1^{1}\cdot (1+\frac{0}{2}+\frac{0}{4}+\frac{1}{8})\cdot 2^{1025 - 1023}\\[1em] -1^{1}\cdot (1.125)\cdot 2^{2}\\[1em] -1\cdot 1.125\cdot 4\\[1em] -4.5 \tag{5}

However, one limitation of this format is that we cannot store any arbitrary floating-point number. Just like we cannot express 1/31/3 with a finite number of decimals, we cannot represent all floating-point numbers using IEEE-754. In the example above, the next largest number we can represent is by increasing the binary number from 0012001_2 to 0102010_2, which gives us 0/2+1/4+0/8=0.250/2+1/4+0/8=0.25. Carrying out the same computation as above, we get:

−11⋅(1+02+14+08)⋅21025−1023−11⋅(1.25)⋅22−1⋅1.25⋅4−5.0(6) -1^{1}\cdot (1+\frac{0}{2}+\frac{1}{4}+\frac{0}{8})\cdot 2^{1025 - 1023}\\[1em] -1^{1}\cdot (1.25)\cdot 2^{2}\\[1em] -1\cdot 1.25\cdot 4\\[1em] -5.0 \tag{6}

So, if we wanted to express, say, −4.75-4.75, we couldn't in this example because we have only used 3-bits to store the fraction. Now, a floating-point number that has a round-off error of 0.250.25, as shown in this example, would be pretty bad! And that's why we use more than 3-bits to store the fraction.

A single precision float has space to store 32-bits. It uses 1-bit for the sign, 8-bits for the exponent, and 23-bits for the fraction (compared to my 3-bit in the example above). That will make the gap between floating-point numbers smaller, but it does also show that we always have a finite representation of floats; we will never be able to represent all possible floating-point values!

However, just like you can seemingly get out of any trouble if you are rich enough, you can also solve the problem here by more advanced hardware. And, since computer hardware followed Moore's law, we have such an abundance of hardware resources available these days that a 32-bit float is rarely used. Instead, we use 64-bit doubles, where we still use 1-bit for the sign, 11-bits for the exponent, and 52-bits for the fraction. This will be enough for all practical cases!

But, even though the gaps between subsequent floating-point numbers that can be represented are getting rather small, they are still there, and we can explicitly show them in C++, which is hiding it from us as best as it can.

Take the following simple code, for example:

1
2
3
4
5
6
#include <iostream>

int main() {
  std::cout << 1e-6 << std::endl;
  return 0;
}

If you compile and run it, you will get 1e-06 printed as an output. Sensible. Now, let's change the formatting a bit:

1
2
3
4
5
6
7
#include <iostream>
#include <iomanip>

int main() {
  std::cout << std::setprecision(17) << 1e-6 << std::endl;
  return 0;
}

This forces C++ to print the ugly truth, i.e. 17 significant digits. The result? Well, it gives me 9.9999999999999995e-07. Isn't that interesting? It is the same number that was stored in the XML file! Whenever we print this number in C++, it does round it for us to the significant digits a double can represent (which is 16 digits), so by default, C++ prints 1e-06. This makes sense. But once we force it to print some more digits, we see the real number that is stored.

Fortran is another compiled language that usually isn't as kind. You will very often see floating-point numbers being represented as above, i.e. 1e-06 may be the input you prescribe, but if you print it to the console, it will give you 9.9999999999999995e-07. C++ is, in that respect, a bit more user-friendly, and it is very rare that we can say that about C++!

Hopefully, this gives you some idea about the round-tripping issues mentioned in a few places already, finally we were able to see it in one place. And also, after working with pugixml to survive this section, I have to say that I changed my view on XML. I still don't like it, but if I can put a firewall between me and XML files so that I never have to deal with them myself, but instead, use an API like the one from pugixml, it does appear to be a capable format.

In other words, I am glad I don't have a prostate exam scheduled for today!

Embedded scripting

At this point, I want to make a bit of a break here. So far, we have looked at different file types that, for the most part, have been pretty similar. Sure, they all have their own syntax, but they also all share, more or less, the same drawbacks. This is where embedded scripting comes into play, which is one alternative of dealing with the shortcomings discussed thus far, specifically the lack of expression and unit support.

By embedded scripting I mean using a scripting language (like Python or Lua) to drive the simulation setup. In a scripting language, it is straightforward to express a time varying inlet profile. For example, in Python, if I wanted to store a time-varying inlet boundary condition, where I set the inlet velocity to a sine wave (for simplicity, think of it as a pulsating inlet condition), I could declare it first and then store it in a list of solver parameters as:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import math

def inlet_profile(t, u0, u_mean):
    return u_mean + u0 * math.sin(t)

config = {
    'boundaryConditions': {
        'inlet': {
            'u': inlet_profile,
            'v': 0.0,
        }
    }
}

print(config['boundaryConditions']['inlet']['u'](1.0))

So, we can easily evaluate expressions and store them in our configuration parameters of the solver. Python also has third-party support to deal with units easily. For example, using pint, a unit conversion library, we could also store units and convert easily between different units as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
from pint import UnitRegistry
ureg = UnitRegistry()

config = {
    'boundaryConditions': {
        'inlet': {
            'u': 50 * ureg("km/h")
        }
    }
}

print(config['boundaryConditions']['inlet']['u'].to("m/s"))
# 13.88888888888889 meter / second ... oh no, not this again ...

Of course, inheritance becomes as simple as using Python (data)classes and inheriting from a base case. For example, earlier we looked at creating a base case in YAML and then deriving a coarse, medium, and fine grid representation of it. In Python, we could implement something like that as well:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
from dataclasses import dataclass, field
from enum import Enum

# helper class for boundary conditions
class BoundaryCondition(Enum):
    DIRICHLET = "Dirichlet"
    NEUMANN = "Neumann"

# meshing
@dataclass
class Meshing:
    numPointsX = 128
    numPointsY = 128

# boundary conditions
@dataclass
class BoundaryConditions:
    north = BoundaryCondition.DIRICHLET
    south = BoundaryCondition.NEUMANN
    east = BoundaryCondition.NEUMANN
    west = BoundaryCondition.DIRICHLET

# solver
@dataclass
class Solver:
    turbulence = True
    model = "komega-SST"
    CFL = [1e-6, 1.0, 10.0]

# define a base case from which to derive and then later modify
@dataclass
class Case:
    meshing = field(default_factory = Meshing)
    boundaryConditions = field(default_factory=BoundaryConditions)
    solver = field(default_factory = Solver)

# coarse mesh, keep the rest of the settings as is
@dataclass
class CoarseCase(Case):
    meshing: Meshing = field(
        default_factory = lambda: Meshing(
            numPointsX=64,
            numPointsY=64
        )
    )

# medium mesh, keep the rest of the settings as is
@dataclass
class MediumCase(Case):
    meshing: Meshing = field(
        default_factory = lambda: Meshing(
            numPointsX=128,
            numPointsY=128
        )
    )

# fine mesh, keep the rest of the settings as is
@dataclass
class FineCase(Case):
    meshing: Meshing = field(
        default_factory = lambda: Meshing(
            numPointsX=256,
            numPointsY=256
        )
    )

There is a bit more involved here, but we could potentially split this over several files or arrange it in a way that makes sense to us, so we only would have to make the changes we want in one file, or provide a nice user interface that then sets all the corresponding entries on the backend for us. While there is more typing required, it does bring us additional flexibility that we can use to drive our solver.

It does have its drawbacks, as you can probably imagine. While it was easy to inspect an input file using JSON, YAML, TOML, INI, and yes, even XML, using an embedded scripting approach means we are now mixing implementations with parameters. If we are not taking good care of what should go into the script itself and what should remain on the backend in the solver, we end up with endless configuration files that become difficult to read.

Mistakes are more difficult to spot and may silently break your simulation. For example, we may want to steer some mesh generator with a Python script. To compute the first cell height based on my y+ requirements, I may want to dynamically compute that instead and set it as a parameter. For example, we may have:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
def first_cell_height(target_y_plus, L, U, nu):
    Re = U * L / nu
    Cf = 0.0592 * Re**(1/5)
    u_tau = U * (Cf / 2)**0.5
    y = target_y_plus * nu / u_tau
    return y

config = {
    'meshing': {
        'wall_patch': {
            'first_cell_height': first_cell_height(1.0, 1.0, 1.0, 1e-6)
        }
    }
}

Looks good? Except, on line 3, it should be Re**(-1/5), not Re**(1/5)! Without noticing this, you may silently introduce wrong behaviour and then start to debug your mesher (which is perfectly fine), and you will probably go just a tiny bit insane trying to fix it.

So, there are excellent opportunities for embedded scripts, but they do come with their own warnings, and you must decide if you are willing to open yourself up to these types of issues, not to speak of the complexity of getting Python, or Lua, integrated with your main C++ code. There are ways, but that is an entire article in its own right.

When it comes to the tables that we looked at before, it is a bit difficult to score embedded scripts here, as it really depends on the scripting language you are using. Here is an attempt anyway. For features, we had:

Table 14
Format Comments Standard Inheritance Expressions & units Round-trip
Embedded script (Lua / Python) ★★★ ★★ ★★★ ★★★ ★

All programming languages support comments, no issues here. Most will have a good standard or steering committee behind it as well, smaller languages may not, but mature languages will not misbehave in unexpected ways. Inheritance and expressions/units can usually be integrated without any issues, so three stars here. Round-tripping does not make any sense here; we are not writing out a Python or Lua file just to read it back again, that only makes sense for static config files.

For the integration burden, we had:

Table 15
Format Error messages Build burden Learning curve HPC / MPI scaling
Embedded script (Lua / Python) ★★ ★★ ★★ ★

Error messages depend on the language again so I give it two stars here, sort of an average, same for the build burden and learning curve. It is difficult to generalise this. Since we are now bringing in an entire scripting language to drive the simulations, the HPC / MPI scaling becomes a question. Sure, there are ways for how we can integrate that with excellent scaling behaviour, but in general, this isn't the case, and we need to put in quite a bit of effort to get there. So, one star only.

As I have mentioned, this is likely its own article to see how we can embed a scripting language into C++ so I won't go into details here and try to cover something in a few paragraphs that requires a bespoke discussion. But suffice to say it is a powerful way to steer your simulation, and more and more CFD solvers are switching to this type of configuration.

SU2, for example, does have a Python driver, Fluent is switching as well to Python and you can set up simulations using PyAnsys and PyFluent (though there is very little reason to switch from the arcane lisp-based scripting to Python). I have offered my frustration to the developer and got a reply, but I don't understand the reply. Maybe you do? What does "Anyone writing small verticals can use settings API and expect some rewrite later." mean? I stopped responding out of confusion.

Be that as it may, the opportunities that a Python-based environment provides do outweigh the disadvantages and so, I believe we will see more and more examples in the future coming up where the backend may be C/C++, Fortran, Julia, dare I even say Rust(?!), while we have a Python front-end doing all the simulation setup for us.

HDF5 / CGNS / NetCDF

Cast your mind back to (this sentence autocompletes in my mind very differently) the time we talked about CGNS and HDF5. We have already looked at these two file types, but not yet NetCDF, which stands for Network Common Data Form. If you have not read either of these articles, or you have no idea about what either format does, let me give you a brief overview of each:

CGNS: The CFD General Notation System, or CGNS, is a data format that aims to be a portable and all-inclusive data format for CFD. By portable, I mean you should be able to write the file in one CFD software, and then read it back in another CFD software. The best analogy I always use is to think of CGNS as the PDF file of CFD simulations, i.e. it is a specific CFD file format that can be understood by any piece of software (like PDFs are universally accepted to display a document).

By all-inclusive I mean a CGNS file wants to store everything a CFD simulation may need. It mostly comes down to the mesh and the corresponding solution, so the dream is you create a CGNS file in your mesher, then read it into your solver to run the simulations on, then write the solution to the same CGNS file, and then look at it in your favourite post-processor.

That's the dream, I don't think it exists outside NASA. Every CFD solver I have ever used sees CGNS as an add-on, a nice-to-have, but not as a first-class citizen. At best, it is there to show support for the cause, but every software then falls back onto using its own file format. In my view, the CGNS file format is (too?!) ambitious, meaning that there are different ways how to interpret it, and that has caused incompatibility issues.

Maybe we'll fix them one day, but when I wrote my first CFD solver and I believed in the same dream, I wanted to adopt CGNS as the main file format and most of the time I spent debugging issues with getting CGNS to work with at least one mesher (ICEM-CFD, oh the joy). I had Pointwise available as well (now bought by Cadence, who rebranded it to the-super-high-fidelity-cadences-in-house-mesh-generator-formerly-known-as-pointwise, or something to that effect (which I refuse to call it)), but Pointwise and ICEM-CFD would write CGNS files differently.

In any case, CGNS works by storing data in a hierarchical, node-based form. Each node has a standard name/type, so it can be read by other software, so the structure of a CGNS file is locked in, you just have to find the right node to write your data into. For example, there is a node for the mesh and one for your solution. Each variable is already defined, so you have to write your data using the predefined variable identifier (VelocityX, VelocityY, and VelocityZ instead of u, v, and w).

In theory, you could also write additional, solver meta-data into the CGNS file if you wanted to store additional solver configurations. The argument is that if you already use CGNS for your mesh and solution data, you may as well store the solver inputs alongside so you can make sense of this later.

However, CGNS files are binary files, and we cannot simply open them, unless we use specialist CGNS inspection tools. So, it is impractical to write solver configurations into them. If you wanted to set up a new case, or change the settings of an existing one, you would have to first launch a different application to access your input file. This is too complicated, only the French would do it this way.

HDF5: This file format is best described as "a compression format for scientific data, on steroids". Virtually all CFD software that want to reduce file sizes use HDF5 these days. Even CGNS is based off HDF5. NetCDF (we are getting to it) also uses HDF5 in the backend.

Just like CGNS, HDF5 files can be thought of as a node-based representation, where we write data (like the mesh, or the solution arrays) under a specific node. Unlike a CGNS file format, an HDF5 file wants to support all types of scientific applications, and so there is no pre-determined structure, and we can nest nodes however we want, or name them however we want.

HDF5 files are also somewhat easier to write than CGNS files, as you can essentially come up with a structure for your data yourself and then just write it, while you have to identify the correct node in CGNS to write your data into.

NetCDF: We have not talked about NetCDF around here on this website. And for good reason. I made the comparison between TOML and INI files before, where I said that INI files are like the toddler brother of TOML files in terms of capabilities. That is, TOML has support for headers, nested key-value pairs, and even types, while INI files have single-header only, no data type support, and everything is a string.

The same analogy holds here where NetCDF is the incapable toddler brother of HDF5, the strong, alpha-male older brother. NetCDF (which always translates to NetCFD in my head, I cannot be the only one ...) is, in essence, a file format that allows us to store arrays. No nested nodes, just plain arrays named arrays. So, in theory, it is very simple to pick up, but just like HDF5, we have to go through a library to write these files and the library follows 1980s rules of using as many cryptic names as possible. It's no joy.

Under the hood, NetCDF has switched to storing its data using HDF5. Let's unpack that. HDF5 offers a wide range and flexibility in storing data. NetCDF offers a fraction of this, and only achieves comparable speeds because it uses HDF5 as well. Brilliant. It's like me releasing my own open-source CFD solver, which is just a bash script calling OpenFOAM routines and solvers, but I have removed all support for unstructured grids, parallelisation with MPI, and even turbulence modelling capabilities.

You may ask yourself, what is the point of NetCDF then? I don't have an answer.

What it does mean, though, is that I have made this section incredibly simple for myself. I don't think there is any point in showing you how to write a NetCDF file, as you will be able to do much better with HDF5 and CGNS. Though, I wouldn't recommend either for storing solver configurations. These file types exist for storing your solution (i.e. array-like structures), and so using these files to store solver configurations is a bit like eating popcorn with chopsticks. Sure, you can do it, but why would you?

In case you must, here is my verdict on its features:

Table 16
Format Comments Standard Inheritance Expressions & units Round-trip
HDF5 / CGNS / NetCDF ★★ ★★★ ★ ★ ★★★

Since everything is binary, comments don't make sense here. But, we can typically provide annotations that do the same job, but we need a specialist program to read them, so these comments may also get lost. All file types have a steering committee behind it as best as I can tell, so we can be confident that the file types themselves will be maintained and do not change overnight breaking our solvers.

Since we store mostly array-like structure, there is no concept of inheritance and expressions/units, though both HDF5 and CGNS have compositions, i.e. we can link several files together (e.g. one mesh file linked with many solution files to avoid storing the mesh for each time step). But, since we store all of that in binary, we get byte-for-byte identical results, so we can be confident in storing data and reading it back exactly as is, good for serialisation tasks (and that is what we would use it for).

In terms of integration into our solver, we had:

Table 17
Format Error messages Build burden Learning curve HPC / MPI scaling
HDF5 / CGNS / NetCDF ★ ★ ★ ★★★

The files don't have any in-built mechanism to support error messages, that would happen at the time of going through the library. As a developer, you may see errors, but not necessarily as a user. Binary files without any error messages are difficult to debug in case things go wrong. Regardless of which file type you go for, there will be a library you have to go through, and these can bring in heavy dependencies (well, mostly HDF5, but CGNS also builds with TK/TCL if you want GUIs for some useful CGNS inspection tools).

HDF5 got its own article, CGNS its own series with 6 articles in total. If I can fill an entire series on one file type, that probably means there is a lot to learn. So it is not something you can pick up like any of the other file types we have looked at thus far, but usually, you also learn these file types because they are foundational to your solver infrastructure.

What it does give you is that extra bit of performance, especially when we are looking at parallel performance. HDF5 has dedicated support for writing distributed data, and since CGNS and NetCDF are built on it, they inherit this at no extra cost to them. It is the reverse of Pareto's law. We say that 80% of the results come from 20% of the work. Here it is the reverse; to get the last 20% of performance out of your CFD solver, you have to spend 80% of your time to implement and fine-tune HDF5, CGNS, and if you must, NetCDF.

My recommendation, though, if you have any say in this: avoid NetCDF, adopt HDF5 for any data storage, not solver configurations, and use CGNS if you have to, or you have a specific use case (i.e. reading a mesh from a mesh generator and CGNS is the easiest file format you can get out of it).

CSV

After we have looked at arguably the most difficult file types (and no one noticed that I managed to go through the last section without writing a single piece of code), CSV files probably do not require any introduction and are so common that I can reasonably assume most, if not all, have worked with them at some point.

CSV, or comma separated values, is just that, columns that are comma separated and used to store arrays, or data, in plain ASCII. So, if we look at our original JSON solver configuration:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
{
  "meshing": {
    "numPointsX": 128,
    "numPointsY": 128
  },
  "boundaryConditions": {
    "north": "Dirichlet",
    "south": "Neumann",
    "east": "Neumann",
    "west": "Dirichlet"
  },
  "solver": {
    "turbulence": true,
    "model": "komega-SST",
    "CFL": [1e-6, 1, 10]
  }
}

we could represent that as a CSV file in the following form:

meshing.numPointsX,meshing.numPointsY,boundaryConditions.north,boundaryConditions.south,boundaryConditions.east,boundaryConditions.west,solver.turbulence,solver.model,solver.CFL 128,128,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]"

If you look at this and say "hang on, this doesn't look very nice", you would be right. CSV files don't look nice, and since we mostly store arrays in columns, reading an input file from left to right, rather than from top to bottom, becomes annoying, especially for a lot of input parameters.

However, one thing we can probably all agree on is that simply separating input values by commas is likely the simplest way of representing data. As Milo Edwards would put it, everything in life is more complicated than CSV files; if you cannot comprehend how CSV works, you may as well dignitas it now. Or at least, something to that effect. He may have used the driving theory test as an example, but I hope you get my point.

There is one big advantage, though, that CSV files provide us. And this is where the left to right parameter reading actually will help us. And that is parametric simulations. Let's say we want to do the same grid dependency study we have looked at in other examples, where the solver configurations stayed the same, but we changed the mesh size only. We could express that in our CSV file as:

meshing.numPointsX,meshing.numPointsY,boundaryConditions.north,boundaryConditions.south,boundaryConditions.east,boundaryConditions.west,solver.turbulence,solver.model,solver.CFL 64,64,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]" 128,128,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]" 256,256,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]"

If we wanted to be kind to ourselves, we could also include a case ID:

caseID,meshing.numPointsX,meshing.numPointsY,boundaryConditions.north,boundaryConditions.south,boundaryConditions.east,boundaryConditions.west,solver.turbulence,solver.model,solver.CFL 1,64,64,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]" 2,128,128,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]" 3,256,256,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]"

In this case, the header describes the parameters, and each row thereafter represents one simulation. Representing parametric studies, in a single file, becomes extremely compact. Not only that, it is also very easy to see how many simulations we have to perform by looking at the number of lines of the file. So there are pros and cons, and for the right solver, CSV files may actually be not as stupid as it may have seemed at first!

In terms of CSV features, we had:

Table 18
Format Comments Standard Inheritance Expressions & units Round-trip
CSV ★ ★ ★ ★ ★★

No comments, only plain columns separated by commas, though, it is so simple that comments are not required. There is no standard, but, again, what would you want to standardise here? The whole file consists of placing a comma at the right position, that's it. A standard would be overkill, just like a steering committee, but still, it has to score as 1 star. Or so I thought, at least.

Some more digging and I came up with RFC, which stands for request for comments and it describes itself as:

RFCs (Requests for Comments) are documents that describe Internet Standards, technical specifications, protocols, procedures, research, and other information related to the Internet and Internet-connected systems. The RFC Series is the archival publication series for Internet technical specifications and related documents.

So, searching for RFC and CSV, sure enough, you will find RFC4180 titled "Common Format and MIME Type for Comma-Separated Values (CSV) Files". So, is there a standard now? Well, the first line in this document states:

This memo provides information for the Internet community. It does not specify an Internet standard of any kind.

So maybe it is like 50% of a standard. There is no one who ratified this RFC, but again, what do you want to ratify in the first place? Put a comma, and that's it; what is surprising is that the RFC4180 achieves a total of 6 pages. The same page will also tell you that there is an updated RFC, RFC7111. That one manages to write 12 pages. 12 pages on comma separated values? Must be a page turner, but I just can't find the willpower to read through this ...

The usual issues with inheritance and expressions/units exist here as well, and since CSV is essentially stored as text (no data types are supported), we have the same round-tripping issue as described with INI files. If you don't touch the values between reading and writing, you get back the exact same document, but if you convert data to floats, you may end up losing precision here.

In terms of integration, we had:

Table 19
Format Error messages Build burden Learning curve HPC / MPI scaling
CSV ★ ★★★ ★★★ ★★

CSV files don't have any way to emit error messages, but when something goes wrong, it is usually because we placed our commas incorrectly. The most difficult to debug issue is likely when you start using whitespaces and expect them to be stripped. For example, if we stored the first two entries in our config file as:

meshing.numPointsX,meshing.numPointsY 128,128

We can read them in a program by querying the headers "meshing.numPointsX" and "meshing.numPointsY". However, we could also store it as:

meshing.numPointsX, meshing.numPointsY 128, 128

Here, I introduced a whitespace between each column. Reasonable, you may think, but reading "meshing.numPointsY" now would fail because we introduced the whitespace. So, either we strip them in our code, or we have to read " meshing.numPointsY". Either way is not great, and so, typically, we store CSV files without whitespaces, making them more difficult to read for humans, but less arbitrary for computers.

There is no build burden, it is likely the only file type you can write your own parser for in seconds, and also, there is no learning curve, really. It scales just as well on HPC as any other ASCII-based file type.

And just as I was writing the last paragraph and while browsing the handful of CSV parsers available on Conan Central, I thought to myself, why don't we have a bit of fun and write our own CSV parser library? Your definition of fun may vary. But it is my website, and, as already pointed out before, no one will read beyond the introductory paragraph, so I can do whatever I want here.

So, let's do exactly that. I am going to call my library Exceln't, which we will pronounce excellent because it really is an excellent CSV parser. In fact, my goal is to make this parser so stupidly simple, that it would irritate Joao beyond belief (and, if you have fallen asleep halfway through this article and are back now, welcome, you may want to check out the section on YAML files (or maybe don't, probably best if you don't ...)).

3 days have passed since writing the last paragraph and oh boy did I have some fun writing my own CSV parser. It is live now in version 1.0.0 on GitHub. First I thought about showing how the code is written, but we are now at about 800 lines of code and there is no guarantee I will leave it at that length, the code may change, I already have improvement suggestions so by the time this article is published, it may have already changed.

So, instead of making the CSV section the longest, which really should be the simplest section to write about, let's look at the interface itself and see how we would use Exceln't in production.

We have our solver configuration given in CSV as:

meshing.numPointsX,meshing.numPointsY,boundaryConditions.north,boundaryConditions.south,boundaryConditions.east,boundaryConditions.west,solver.turbulence,solver.model,solver.CFL 128,128,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,"[1e-6,1,10]"

I am going to make one change here. The CFL column is now packaging an array (which doesn't exist in CSV), so I will split the values up into what they actually mean: min_dt (which is set to 1e-6), max_dt (which is set to 1), and the CFL number itself, which has a value of 10. We can write this with Exceln't in the following C++ code:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include "excelnt.hpp"

int main() {

  excelnt::Table table;

  // meshing
  table["meshing.numPointsX"] = 128;
  table["meshing.numPointsY"] = 128;

  // boundaryConditions
  table["boundaryConditions.north"] = "Dirichlet";
  table["boundaryConditions.south"] = "Neumann";
  table["boundaryConditions.east"]  = "Neumann";
  table["boundaryConditions.west"]  = "Dirichlet";

  // solver
  table["solver.turbulence"] = true;
  table["solver.model"] = "komega-SST";
  table["solver.CFL.min_dt"] = "1e-6";
  table["solver.CFL.max_dt"] = "1.0";
  table["solver.CFL.CFL"] = "10.0";

  excelnt::write("config.csv", table);

  return 0;
}

You have now seen about 90% of the syntax that is supported. There is a bit more to come in the reading code of a CSV file, but it doesn't get more complicated than that. We define a table, which has a type of excelnt::Table. That is the only type that is exposed. Once a table is defined, we can start adding entries to it using the square brackets ["column_name"]. We then assign a value for that column. We use excelnt::write("csv_file.csv", table) to write out a table object to file.

In this example, we assign single scalars only, as this is the use case we have for our configuration, but we could also pass here a std::vector and store an array of values. For example, in the mesh dependency example we have looked at now a few times, we may change the number of points in x and y as:

1
2
3
// meshing
table["meshing.numPointsX"] = std::vector<int>{64, 128, 256};
table["meshing.numPointsY"] = std::vector<int>{64, 128, 256};

You may ask if we have to change the values of the other variables as well. For example, we had:

1
table["solver.turbulence"] = true;

You may reasonably expect that we now would have to write:

1
table["solver.turbulence"] = std::vector<bool>{true, true, true};

No! Exceln't got you covered and allows for broadcasting, a term I stole from NumPy. Broadcasting means here that you may have a column with a number of entries, as we do for the number of mesh points, for example, and then other columns with just a single scalar. Exceln't will then internally promote single scalars to arrays automatically which will have the same length as the other arrays.

So, if we defined the number of mesh points to have 3 entries and then we write:

1
table["solver.turbulence"] = true;

Exceln't will internally promote this to:

1
table["solver.turbulence"] = std::vector<bool>{true, true, true};

You don't have to do anything. For the sake of completeness, let's look at a complete script:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include "excelnt.hpp"

int main() {

  excelnt::Table table;

  // meshing
  table["meshing.numPointsX"] = std::vector<int>{64, 128, 256};
  table["meshing.numPointsY"] = std::vector<int>{64, 128, 256};

  // boundaryConditions
  table["boundaryConditions.north"] = "Dirichlet";
  table["boundaryConditions.south"] = "Neumann";
  table["boundaryConditions.east"]  = "Neumann";
  table["boundaryConditions.west"]  = "Dirichlet";

  // solver
  table["solver.turbulence"] = true;
  table["solver.model"] = "komega-SST";
  table["solver.CFL.min_dt"] = "1e-6";
  table["solver.CFL.max_dt"] = "1.0";
  table["solver.CFL.CFL"] = "10.0";

  excelnt::write("config.csv", table);

  return 0;
}

If you run this code, you will get:

meshing.numPointsX,meshing.numPointsY,boundaryConditions.north,boundaryConditions.south,boundaryConditions.east,boundaryConditions.west,solver.turbulence,solver.model,solver.CFL.min_dt,solver.CFL.max_dt,solver.CFL.CFL 64,64,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,1e-6,1.0,10.0 128,128,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,1e-6,1.0,10.0 256,256,Dirichlet,Neumann,Neumann,Dirichlet,true,komega-SST,1e-6,1.0,10.0

All entries are the same, and have 3 entries, except for the meshing entries, which all have their individual settings.

Reading a CSV file is just as painless:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <iostream>

#include "excelnt.hpp"

int main() {

  const auto table = excelnt::read("config.csv");

  // meshing
  const auto numPointsX = table["meshing.numPointsX"].asInt();
  const auto numPointsY = table["meshing.numPointsY"].asInt();

  // boundary conditions
  const auto north = table["boundaryConditions.north"].asString();
  const auto south = table["boundaryConditions.south"].asString();
  const auto east  = table["boundaryConditions.east"].asString();
  const auto west  = table["boundaryConditions.west"].asString();

  // solver
  const auto turbulence = table["solver.turbulence"].asBool();
  const auto model = table["solver.model"].asString();
  const auto minDt = table["solver.CFL.min_dt"].asDouble();
  const auto maxDt = table["solver.CFL.max_dt"].asDouble();
  const auto CFL = table["solver.CFL.CFL"].asDouble();

  // sanity check output (everything is read as a std::vector, so we print the first (and only) element)
  std::cout << "numPointsX = " << numPointsX[0] << std::endl;
  std::cout << "numPointsY = " << numPointsY[0] << std::endl;

  std::cout << "north = " << north[0] << std::endl;
  std::cout << "south = " << south[0] << std::endl;
  std::cout << "east  = " << east[0]  << std::endl;
  std::cout << "west  = " << west[0]  << std::endl;

  std::cout << "turbulence = " << std::boolalpha << turbulence[0] << std::endl;
  std::cout << "model = " << model[0] << std::endl;
  std::cout << "min_dt = " << minDt[0] << std::endl;
  std::cout << "max_dt = " << maxDt[0] << std::endl;
  std::cout << "CFL = " << CFL[0] << std::endl;

  return 0;
}

We read a file by calling const auto table = excelnt::read("file_name.csv");, and then we can access any key in that table by calling the square bracket [] operator again on the table object. When we read something from a CSV file, we need to specify the type we want to receive, and we have asInt(), asDouble(), asBool() and asString() available. Since everything is stored as a string in CSV, the asString() conversion is a no-op but is is there for completeness.

I also show how you would access each variable by printing it. Since everything is returned as a std::vector, we have to access each variable with the additional [0] to specify we want to get the first (and, in this case, only entry) of that vector. If we added more columns to our CSV file, e.g. we added the additional mesh points, then we could still use this code, it would then only print the first entry, though there may be more. If we had more than one entry per column, we could loop over all using, for example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
// get column with more than one element per column
const auto numPointsX = table["meshing.numPointsX"].asInt();

for (const auto point : numPointsX) {
  std::cout << "Number of mesh points in X: " << point << std::endl;
}

// in our example with 3 columns, this would print:
//
// Number of mesh points in X: 64
// Number of mesh points in X: 128
// Number of mesh points in X: 256

Congratulations, you are now a power user of Exceln't. You may now put your name on this certificate and frame it next to your other degrees on the wall!

Your Exceln't certificate
Figure 5

CBOR and MessagePack

Last week (as of the time of this writing anyway), I had to derive the Taylor Series in one of my lectures. And the name Taylor may mean different things to different people. I am living in a cave, apparently, as I only found out in about 2025 that Taylor Swift is a thing, and I did so by reading a book from someone who boasts they made a website for her.

I have since sampled a few tunes. It is OK-good (a throwback to the good old days of DigitalRev TV with Kai and Lok, anyone watched that? Just me?). Slightly above elevator music, slightly below Spotify advertisement.

Anyway, back to the main plot. As a child, I would look through the catalogue of my local music shop (which, I just learned, went bankrupt during covid, sad emoji), and for me, Taylor was synonymous with that guitar brand that was just miles out of reach. I mean, their cheapest guitars were like €1,000 even back then, but they were also selling rare and limited edition guitars anywhere north of €50,000. Makes the €1,000 guitars look like a steal!

Of course, these days, Taylor for me is synonymous with Taylor's of Harrogate, producer of Yorkshire Tea (which I just noticed is absent from the cfd.university creative department, aka the tea cupboard).

But, there is a new sheriff in town, and from now on, when we speak of Taylor, we are talking about Taylor from Taylor Town. I mean, that web design, already in my top 10, and we have only just met! If you think my writing is a bit much sometimes, then wait until you meet Taylor. I think his website deserves an 18+ age restriction, not for the nudity (I am sure I will find some of that eventually as well on there, I wouldn't be surprised anymore), but for his use of language.

In any case, as it so happens, our new friend Taylor here had to pick a serialisation format at some point, and he did all the research already, so I don't have to. In one of his articles, he compares CBor and MessagePack. So, if this is the first time you are hearing about these two serialisation formats, let me ask Taylor to introduce them for us.

Let's start with MessagePack:

Obviously MessagePack is what cool kids would use.

So far so good, what about CBor:

Everything about CBOR is uncool. It was designed by a committee. It reeks of RFCs. Acronyms are lame. Saying "SEE-BORE" is like licking a nickel. One of the authors is "Carsten Bormann", which makes the name feel [erm, Taylor, I don't think I can say this around here ...].

I was intrigued. Now, Taylor and I have different world views. A committee and RFCs are lame, sure, I get it, but without them, we get that one eccentric programmer who changes the interface between every release. I don't know his name, but I know he works at ANSYS, and he (or she, how politically correct of me) keeps messing up all of my Fluent journal files, and I'd rather have him (or her) steer a lame committee that provides stability 30+ years after the solver was first written.

So what are these formats? Well, in a nutshell, MessagePack and CBor allow us to represent JSON-like data as binary data. The advantage of that is that we no longer have any issues with round-tripping. All values are stored in binary exactly and so read back exactly in the same way. It also makes file sizes much smaller.

The example JSON file we were looking at is relatively small, but MessagePack and CBor are serialisation formats, meaning they target JSON files generated by a program, and these can get arbitrarily large. Reducing hundreds of Megabytes into only a few using binary will make a difference in the long run.

Unfortunately, it also means that it is more difficult to reason about these file formats, as I am not very good at writing raw 0s and 1s to binary files. But, this is where our good friend Niels is coming back to save us. You remember Niels, right? The JSON guy. Now, let's assume we have our JSON file on disk, stored as config.json, which, as a reminder, looked something like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
{
  "meshing": {
    "numPointsX": 128,
    "numPointsY": 128
  },
  "boundaryConditions": {
    "north": "Dirichlet",
    "south": "Neumann",
    "east": "Neumann",
    "west": "Dirichlet"
  },
  "solver": {
    "turbulence": true,
    "model": "komega-SST",
    "CFL": [1e-6, 1, 10]
  }
}

We can use nlohmann::json to write and read MessagePack and CBor data directly. Isn't that fantastic? So, let's have a look. The following program will first read the config.json file from disk and then convert the data to MessagePack and CBor, write it to disk as a binary representation, then read it back again, convert back to JSON, and finally, print one variable to sense check that all of them have the same result:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include <cstdint>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

#include <nlohmann/json.hpp>

int main()
{   
  // read the config file generated in the JSON section
  std::ifstream inFile("config.json");
  nlohmann::json config;
  inFile >> config;

  // convert the json object to cbor and messagepack
  auto cborData = nlohmann::json::to_cbor(config);
  auto msgpackData = nlohmann::json::to_msgpack(config);

  // write cbor and messagepack created objects to disk as a binary file
  std::ofstream cborWrite("config.cbor", std::ios::binary);
  cborWrite.write(reinterpret_cast<const char*>(cborData.data()), cborData.size());
  cborWrite.close();
  
  std::ofstream msgpackWrite("config.msgpack", std::ios::binary);
  msgpackWrite.write(reinterpret_cast<const char*>(msgpackData.data()), msgpackData.size());
  msgpackWrite.close();

  // reading the same binary files back from disk now
  std::ifstream cborConfig("config.cbor", std::ios::binary);
  std::vector<uint8_t> cborRead((std::istreambuf_iterator<char>(cborConfig)), std::istreambuf_iterator<char>());

  std::ifstream msgpackConfig("config.msgpack", std::ios::binary);
  std::vector<uint8_t> msgpackRead((std::istreambuf_iterator<char>(msgpackConfig)), std::istreambuf_iterator<char>());

  // convert their content to json
  auto cbor2json = nlohmann::json::from_cbor(cborRead);
  auto msgpack2json = nlohmann::json::from_msgpack(msgpackRead);

  // sanity check, did we get the same number of points back after serialisation? 
  std::cout << "Property   | JSON | CBOR | MessagePack" << std::endl;
  std::cout << "numPointsX | " << config["meshing"]["numPointsX"] << "  | ";
  std::cout << cbor2json["meshing"]["numPointsX"] << "  | ";
  std::cout << msgpack2json["meshing"]["numPointsX"] << std::endl;

  return 0;
}

Let's go through this block by block. In the first block, we simply read the JSON file and make it available as a config object. Next, we write the data to MessagePack and CBor using the inbuilt converters:

1
2
3
// convert the json object to cbor and messagepack
auto cborData = nlohmann::json::to_cbor(config);
auto msgpackData = nlohmann::json::to_msgpack(config);

In the next block, we write the results as binary files to disk. We create first an std::ofstream object for writing files, which we provide with the std::ios::binary flag. This allows us to write in binary instead of ASCII text. We write both CBor and MessagePack data in the same way. Binary files need to be converted (reinterpret_cast) to a char* or const char* before writing.

The way I remember it, and it may not be the most scientific explanation, we convert to char because the size of a char is 1 byte, and in binary, we write anything as a sequence of 1 byte chunks. So, if I had been around during the design of the C++ interface, I would have gone with something like binary_cast(cborData.data()), but Bjarne (the C++ inventor, I think we are on a first name basis) obviously thought reinterpret_cast<const char*>(cborData.data()) rolls better off the tongue. Who am I to question his authority and expertise?

We then also state how many entries we want to write by calling the size() function on each object.

1
2
3
4
5
6
7
8
// write cbor and messagepack created objects to disk as a binary file
std::ofstream cborWrite("config.cbor", std::ios::binary);
cborWrite.write(reinterpret_cast<const char*>(cborData.data()), cborData.size());
cborWrite.close();
    
std::ofstream msgpackWrite("config.msgpack", std::ios::binary);
msgpackWrite.write(reinterpret_cast<const char*>(msgpackData.data()), msgpackData.size());
msgpackWrite.close();

OK, now that we have written the data, we need to read it again. This is somewhat more involved and, as it is not every day that I am writing binary files, especially unsigned integers specifically, I needed to look that up as well and it took some time. But, thankfully there was a very helpful GitHub Gist that allowed me to come up with the following:

1
2
3
4
5
6
// reading the same binary files back from disk now
std::ifstream cborConfig("config.cbor", std::ios::binary);
std::vector<uint8_t> cborRead((std::istreambuf_iterator<char>(cborConfig)), std::istreambuf_iterator<char>());

std::ifstream msgpackConfig("config.msgpack", std::ios::binary);
std::vector<uint8_t> msgpackRead((std::istreambuf_iterator<char>(msgpackConfig)), std::istreambuf_iterator<char>());

Here, the std::istreambuf_iterator reads binary inputs (again, using char here to indicate that we are reading chunks of 1 byte-sized data). We see that we are using an iterator here, and if we look at the constructors of std::vector, we can identify that we are using the fifth defined constructor.

This constructor requires us to provide a start and end iterator. If you are familiar with the standard template library in C++, i.e. the STL, then you may be familiar with the begin() and end() functions that provide exactly that for us. So, a std::vector<double> data has data.begin() and data.end() defined, and we can always iterate over all of its entries using these two iterators.

Using a binary representation complicates things a bit (conceptually), and the std::streambuf_iterator works somewhat differently. First, we point it at the first entry, which we do with the std::istreambuf_iterator<char>(cborConfig) entry. By default, if we just write std::istreambuf_iterator<char>(), this is the end-of-stream iterator, so this becomes our second argument.

There are no clever ways to remember this. Either you spend the next 3 years writing all your outputs in binary and develop a level of acceptance (but never a level of agreement) that we are writing binary files the way we do, or you just live with the fact that binary numbers work differently to the rest of C++. In many ways, binary numbers are the UK of Europe, a reference you may only get if you have lived in either (or both) part(s) of the world.

The UK is a weird place. I love it here!

With the binary blobs read back into C++, we can now pass all of that back to our JSON library and forget about binary numbers, touch grass, and hope it has not reduced our life expectancy too much.

1
2
3
// convert their content to json
auto cbor2json = nlohmann::json::from_cbor(cborRead);
auto msgpack2json = nlohmann::json::from_msgpack(msgpackRead);

I am printing back values in the code above to sense check that the same data has arrived again. And it does. We should remember at this point that we wouldn't necessarily store our input files in either CBor or MessagePack format, that would be the literal definition of insanity (or being French? okok, I stop now ...).

The reasons I am showing this here are two-fold:

  1. I am a little bit in love with Taylor and needed a reason to talk about Taylor (I let you decide which Taylor. Yesterday I read a book with my daughter and the author's name was also Taylor. Taylor is everywhere ...)
  2. We talked about serialisation at the beginning of this article but thus far have not really looked at any suitable option. In my view, CBor and MessagePack are the options you want to look at.

I have looked at libraries outside nlohmann::json to work with CBor and MessagePack, but either couldn't find anything I liked, or perhaps I am just getting library-hopping-fatigued, which is a real medical condition I just made up (I have a PhD, I can call myself a Doctor, right?). So, my recommendation is to stick with this library. Once we have mastered the art of reading and writing binary data that either format emits, we have all we need to serialise data.

Moving on, in terms of our summary table, for the features, we had:

Table 20
Format Comments Standard Inheritance Expressions & units Round-trip
CBOR / MessagePack ★ ★★ ★★ ★ ★★★

Since it is binary, no comments, which is by design. For the standard, CBor has one, MessagePack doesn't, but it does have a specification. So perhaps 2.5 stars would be appropriate here, but since I can only print integer representation of stars, I have settled for 2 (after all, this isn't Dinner Date!). Inheritance works no better than JSON, no expression and units as usual, but the round tripping now works, thanks to the binary data representation.

In terms of the learning curve, well, this is what we had:

Table 21
Format Error messages Build burden Learning curve HPC / MPI scaling
CBOR / MessagePack ★ ★★★ ★★★ ★★★

Error messages aren't really part of the standard, but, if we reuse nlohmann::json, as I have done, then there is no build-burden or learning curve involved, we simply dump our JSON data from within a JSON library we already know. Binary files scale better than ASCII data, so HPC / MPI scaling scores well here.

So, if we want to serialise data (e.g. the content of a class to be read back when we start a solver again and want to resume a previous simulation), then CBor and MessagePack are both excellent options we get for free if we use nlohmann::json.

For completeness, I should add that nlohmann::json supports other serialisation formats as well (see a list of them here). So, we could use any other format as well. I haven't checked their syntax, nor the file formats themselves, but would assume they work in a similar manner to how we read and wrote CBor and MessagePack data. If anything, it will not be much more complicated.

And with that, we have reached the end of all the file types I wanted to cover. There is just one section left, which lists some honourable mentions here.

Honourable mentions

I could go on, talk about more file formats, more libraries, and more ways to input and output data until even some of the most loyal readers of cfd.university will report me to the authorities for mental abuse, so I'd rather stop here. But there are a few formats I want to mention before we move on, as I think they deserve at least a mention to make this list complete (it can never be complete, but as close to complete as possible to satisfy our OCD).

Let's start with OpenFOAM, the elephant in the room. OpenFOAM could only have been conceived in the UK, it wants to be different, and every user new to OpenFOAM has to learn, essentially, an entire new way of setting up simulations that feels unlike anything else they have done before. OpenFOAM's input file (or, shall I say, input file system?) is a collection of folders, files, bash scripts, and references that resolve to other folders, files, and references.

In many ways, if we see OpenFOAM as its own operating system, I think we have an accurate description of what it actually is (too complicated for what it needs to be).

The inputs to OpenFOAM are complicated. If you have never used it, I compared it once against ANSYS Fluent, comparing how many clicks I needed to do in Fluent to set up a case in its simplest form, and then checking how many lines of input I need in all of my OpenFOAM input files. In Fluent, I needed 20 clicks to get the simulation started, which translated to 700 lines of OpenFOAM input, a lot of which was just boilerplate.

However, hidden in that complexity is, I have to admit, a beautiful expression syntax that allows you to do all the things we were unable to achieve before with any of the other file formats. For example, unit (checking) is done by default. You cannot run a simulation with wrong units. While that is a nice feature, it only checks units, it does not offer conversion, which would have been a really nice feature.

Equally, expressions are pretty painless and extraordinarily powerful. Instead of writing your input, for example, a freestream velocity of (50 0 0), you could define a variable as #define AOA 5, for example, which represents the angle of attack (AOA), and then calculate the velocity component based on value. To give you a full example, let's define the inlet velocity as:

1
2
3
4
5
inlet
{
    type            patch;
    value           uniform (50 0 0);
}

Then, we can transform that into:

1
2
3
4
5
6
7
#define AOA 5

inlet
{
    type            fixedValue;
    value           uniform (#calc "50*cos($AOA*M_PI/180)" #calc "50*sin($AOA*M_PI/180)" 0);
}

Simply throwing a #calc "" in there means we can now evaluate an expression in place, and I always admired OpenFOAM for that. That is an elegant extension to your input script machinery. In the background, each #calc "" statement compiles into a shared library and is then hooked into the simulation. What I most admire about this is that I have used this type of customisation in a lot of my OpenFOAM simulations, and the compilation and linking has worked flawlessly every time.

To be fair, though, that isn't really a challenge if you lock your users into using UNIX-only, if only they also attempted to build on Windows, then I would have high respect for them.

OpenFOAM has, somewhat accidentally, developed a domain-specific language (DSL), at least sort of. This is the learning curve that has plagued now generations of OpenFOAM users where every new user has to learn a new programming language, essentially, where learning the syntax isn't enough, you also need to be a CFD wizard in order to set up simulations and understand what is going on.

I hope, one day, I will be competent enough to use OpenFOAM. 언젠가 ...

While Sam Altman is busy building the ultimate Armageddon machine, we need to learn to submit to our AI overlords and accept them in our lives with open hands (all of a sudden, this has become very North-Korean). So why not embrace the powerhouse that is Google and start to consume all of their great developments?

As it turns out, Google has given us both ProtoBuf and FlatBuf (well, FlatBuffers, but I like these two names more). While buf, I am sure, refers to "Buffers", I like to think of them as two buff mates, "Proto" and "Flat", who are not just very muscular, but also enjoy drinking tea from cups that are just comically too small for them.

Proto and Flat are two serialisation formats. Google describes them as:

Protobuf is a language-neutral, platform-neutral extensible mechanism for serializing structured data.

And

FlatBuffers is an efficient cross platform serialization library for C++, C#, C, Go, Java, Kotlin, JavaScript, Lobster, Lua, TypeScript, PHP, Python, Rust and (Taylor, is it you?) Swift.

I wonder if developers at Google actually talk to each other or are allowed to peek into what other people are doing. Do we really need two? In any case, memory efficient or not, cross-platform or not, there is only so much space and attention I can devote to different kinds of serialisation formats. I looked through their documentation and the tutorials (so you don't have to) and decided, for you, that you don't need (sticking to the North Korean theme ...).

We have looked at CBor and MessagePack, they offer a different syntax, one that we can master easily, so we don't need yet more serialisation formats. That is, of course, only true for now. Once Gemini has overtaken ChatGPT's market share, we shall all switch to Google only products.

Side note, to lighten the mood just a tad, while looking for market shares, I found this article with some numbers. I love that Copilot is still listed on place 7. The table also mysteriously cuts off at position 7. It almost feels that whoever collected the data and wanted to share this wanted to show that Copilot is at the bottom of the table. They could have only shown the top 5 or top 10, but decided, purposefully, that the table should have 7 entries. Respect!

So, the ideal file format here is one that is both HDF5-based (storage optimised) and understands CFD. In my view, somewhat unfortunately, CGNS is the only file type that offers us this, while also being open-source and well documented.

Now, if you wonder why I say unfortunately, it is due to my love-hate relationship with it. I really want to like it, but it has become a file format that is too complex for its own good, and it tries to cater for every possible CFD solver in the world, meaning there is just a little bit too much flexibility in the standard. It is not a serialisation format like Proto and Flat, but more like HDF5 and NetCDF. In fact, it states on its GitHub page that it uses HDF5 underneath the hood. Any binary file format that stores scientific data and uses HDF5 under the hood is just more wrapping around a perfectly fine, albeit version-confusing, file format.

I didn't even bother to read the documentation this time, I think I already saw enough. But, I mention it here as yet another HDF5 imposter so you can carry forward my opinionated views at the next cocktail party where, inevitably, the conversation will switch to file formats at some point (mostly because you have a great urge to discuss it with people who do not share the same obsessions as you).

I mean, we are now, what, 2 hours into this article and you are still reading. I know there is something wrong with me for writing this long about file formats, but for you, there is still hope, you should get checked out by a professional before it is too late (and CFD starts to consume your life).

Now that I have shown you all there is to know about file formats, would it surprise you that I actually once wrote a standard to write my very own file format? No? Good, because that's what I did. Another one of my weird passions I developed during Covid. Now, I thought, for a bit of fun, I would actually write a parser library for this file format, but the standard sits on a private GitHub repo and was never completed (it is like 80% written). I am not sure if I have it in me to complete the remaining 20%, so let's not bother.

The core idea was something like JSON or TOML, but with native inheritance a la YAML, and an easy way to express parametric studies. So, for example, we could define an input parameter like "AOA = [0, 10, 11]", and it would now create 11 different input files internally (within the parser library) and then return one each with an AOA value of 0, 1, 2, ..., 10 to whoever would consume this library. This was the main selling point and reason to develop it, but yeah, never got around to finishing it.

One of the days when I was working on it, a friend of mine came around (has the statute of limitations expired for friends visiting you during Covid?) and I told him all about my new standard. He was confused. Maybe that is the reason I stopped working on it, don't know. We studied together, we both got our PhD in CFD, I thought he would appreciate. Apparently not.

But I never gave it a good name, so I thought of proposing some here, seems only logical. Given that it is essentially a superset of TOML, I first thought of calling it TOML++. Clever, an extended version of TOML and expressing that the same way C++ differentiates itself from C. But then I thought, why go for TOML++ if Tom's Integrated Templating Syntax is on the table as well?

The templating referring to the fact that we can write an array template which will internally be expanded into the parameters needed for a parametric sweep. What's the acronym? TITS. Oh. Maybe it is time to let go of the idea and the standard for good ...

Recommendation: What should you use for CFD applications?

Now that we have looked at the different file formats and their pros and cons, I thought it would be good to come up with a recommendation table of which file format to use and in what situation. If you like, we have been building up to this section over the past few hours of pure reading joy.

I have split my recommendation into what file format may be best used for different situations, split by typical file sizes and how it is used. Here it is:

Table 22
Case setup Mesh Restart / checkpoint
Size kilobytes megabytes–gigabytes gigabytes–terabytes
Written by a human a mesher a solver
Read by a human, often your solver your solver only
Needs comments yes no no
Needs to be diffable yes no no
Needs parallel I/O no yes absolutely
Recommendation TOML / YAML / JSON CGNS / HDF5 / other formats HDF5 / CBor / MessagePack

So, case setup files are typically of the order of kilobytes, and so well within the size we humans can still browse, modify, and write. Once we get to files generated by a meshing application or our CFD solver, they are getting too large, and so we really want to start to think about binary files, i.e. optimised formats for large data storage.

In this category, we probably would like to use something that uses HDF5 under the hood, or straight up HDF5 itself. The issue with HDF5 itself is that it does not provide any standard for how to name things, and that is on purpose; it wants to be application-independent. That's good, it means we in CFD can use it, but the other scientific computing niches can use it for their purposes, too.

So, the ideal file format here is one that is both HDF5-based (storage optimised) and understands CFD. In my view, somewhat unfortunately, CGNS is the only file type that offers us this, while also being open-source and well documented.

Now, if you wonder why I say unfortunately, it is due to my love-hate relationship with it. I really want to like it, but it has become a file format that is too complex for its own good. It tries to cater for every possible CFD solver in the world, resulting in too much flexibility in the standard. It is possible to implement the same thing in different ways, making it a nightmare trying to read in CGNS files generated from different mesh generators.

Granted, that applies only for very specific things and, for the most part, the standard is tight enough that you can actually use it, but it is these little inconsistencies that really bug me out. I have spent too much time going through their website, browsing the high-level API and low-level SIDS commands trying to figure out how to achieve what I want. And all of that frustration probably still sits with me.

What can I say, I'm getting old, I am German, I have opinions ...

There is no good alternative, though. OpenFOAM writes everything in ASCII by default (though it does support binary), Fluent uses their own HDF5-based format, and Star-CCM+ decided they don't need the cfd.university wisdom and crammed everything into their own proprietary .sim file format. The mesh, the configuration, the solution, and probably also the lost sanity of the person implementing that file format.

Converge CFD has taken the hilarious decision that it doesn't need any post-processing capabilities and so every Converge CFD licence now comes with Teapot 360 to handle all of their post-processing. You have to use Teapot with Converge, whether you like it or not. In many ways, that is a better definition for a shotgun wedding, and we shall petition for the meaning to be changed.

I know it is Tecplot 360, not Teapot 360, but my new grammar tool does not know Tecplot yet and so suggests Teapot, which I find, in many ways, much better.

We should also not forget that we have primarily focused our attention here on writing input files, i.e. those we humans would most likely author and use to communicate with our solver. That is, the interface between us and the machine. Sure, we also looked briefly at CBor and MessagePack for serialisation, mostly because I mentioned it in the introduction, but also because once we have data in JSON, serialising our data is done with ease.

I have, on purpose, not touched any other file formats like writing mesh or solution files. That would be an article nobody would be willing to read. Paraview can write 25,000 file formats, at least that's a rough estimation on my behalf, of which, if you are lucky, half of them work or are implemented to a point where it supports the full feature set of that file format. The other half are variants of .vtk files, and I am struggling to remember which one to write for which type of scenario/application.

So, we are left with this wonderful mess of a scattered landscape of different file formats and no consistency. We tried (CGNS) and we failed. I am sure the good folks over at NASA (CGNS inventors) disagree. I don't care anymore.

So, gun to my head, which file type should you use? Start with TOML for all input/output. It has been developed specifically for humans and does not have any missing features compared to the other file formats. Well, it does not offer inheritance, that is a trade-off I can accept.

For anything else, I would probably still go with CGNS, or, HDF5 but with an additional XDMF3 layer which provides the structure an HDF5 file lacks, meaning a post-processor like Paraview can understand your HDF5 file. But that locks you in to use Paraview, Tecplot has no support for it, and that isn't a great strategy if you want to adopt a file format anyone can process.

Paraview and Tecplot both support CGNS, but Paraview struggles with the more advanced features like linking. It's simply not implemented. That, at least, was the status from 2013 when I tested it, and I did test it a few years later, it was still not implemented. Maybe it will change one day and the CGNS standard will be supported in its entirety.

Funnily enough, the only data types that both Tecplot and Paraview support here are either OpenFOAM simulations (and OpenFOAM itself does not have a file format, their solution files are written in their DSL-like layer) or legacy ASCII documents like Paraview's .vtk and Tecplot's .dat files. Well, I find that funny anyway, I probably just have a very specific sense of humour.

File format and configuration antipatterns

This section is my experience working with different CFD solvers over the last decade. I have been abused by software engineers that had to meet a deadline and were cutting corners when implementing a flawed, rushed, and sometimes broken user interface. I suffer from Stockholm syndrome because I still hang out with my captors (Fluent, OpenFOAM, SU2, Star-CCM+, Converge CFD, Paraview, Tecplot, Pointwise, just ANSYS (in general), ...).

In many ways, this is the #MeToo of CFD campaign, and it is time we talk about all of the nasty practices that should stop. You may not have thought about them, you may have never experienced any of the below. But, I can promise you, once you get intimately familiar with any piece of CFD software, you too will suffer. So let's talk about them.

Sin #1: Parsing configuration files mid-solver run

Imagine you are on an airplane and all of a sudden, there is an unexpected change to the flight plan. If you wanted to casually walk up to the pilot and politely ask them to divert, your request would probably be politely declined. But, of course, if you bring explosives to the party and threaten to use them, you'd be surprised how much attention and cooperation you will receive.

In general, I would always advise you to travel with an explosive in your hand luggage. It is one of the best insurance policies you can have because the likelihood of having two individuals on board with a bomb is vanishingly small. It might be difficult to argue your case at security, but if you give them my explanation, I am sure they see the point.

Be that as it may, I think we can all agree, we'd all be much rather on a plane without any diversions and threats of explosives. If I board an easyJet flight to Luton airport, the pilot better takes me there. But then again, if I had a choice between landing in Luton or being diverted anywhere else, I am not sure what sounds better. I might go with the latter.

In our CFD solver, it is the same. When I tell it to run for a certain number of iterations, or time steps, and use a specific scheme or model, I would, ideally, like it to do that and not switch from an incompressible to a compressible solver, or from a k−ωk-\omega SST to a Spalart-Allmaras model mid-solver run. If I change to compressible flows, my solver will crash because it doesn't have any idea how to treat density, if I switch to Spalart-Allmaras, then my solver doesn't have any boundary conditions for ν~\tilde{\nu}.

The end result? Best case, this request is ignored, worst case, the solver which was running for a few days or weeks now crashes without having stored any intermediate solution files.

Now, you may say that this is a bit of a farfetched example, why would that change mid-solver run? Because in OpenFOAM, for example, this isn't a bug, it is an explicit feature. They probably even have it as one of their USPs (unique selling points). You can go into any config files, change it, and OpenFOAM will read it mid-run, parse it, and make the necessary modifications.

It is absolutely possible to make changes mid run, but I ask, what is the point? Is it so we can start a simulation and, once it is converging, change it to a higher-order numerical scheme to reduce the numerical dissipation? If so, why not just run the simulation for a set number of iterations, make the change after we finished this initial run and then restart with the new scheme selection?

I think it is important we are able to have a deterministic way of running our simulations so they can be reproduced. You may then argue this is more for dynamic testing, but that practice just teaches the wrong behavioural patterns. We shouldn't be making any adjustments mid-simulation because we introduce randomness at best, and crash our simulation at worst.

Thankfully, we can switch this behaviour off, and there should never be a reason to switch it on in my view.

Sin #2: Silent configuration defaults that aren't documented

Many solvers have a function that is called something like lookupOrDefault(). I have used that type of function in my own solvers, I am not a saint. This function looks for a value in an input file and, if it cannot find one, assumes a default value. Now, let me give you an example that happened to one of my students, I cannot remember the solver we used, but the issue was not solver-specific, it can happen to anyone with any solver.

The student was investigating a case for which we had different results available for comparison at different Reynolds numbers. Then, instead of changing the density or the inlet velocity, the student opted to change the viscosity instead to change the Reynolds number. When the student did that, they mistyped the variable name for the viscosity. It wasn't a typo (syntax error), but a flawed logic (semantic error).

The student changed the parameter mu (μ\mu), but the solver defined the Reynolds number using the kinematic viscosity nu (ν\nu). It may have been the other way around, I don't recall. The behaviour of the solver was, if nu was not provided, rather than crashing and reporting that as an error, it would silently adopt a default value.

The problem was that we did not just change the Reynolds number between runs, but also the Mach number. And so, there were differences between runs, but the differences were entirely due to the Mach number, and not capturing the Reynolds number effects. We had to do that because the reference data was provided in this way, and we wanted to match that data first before moving on.

Suffice it to say, it took the student some days to track down the issue. If the solver was just honest about a missing viscosity value, it could have informed the user about it. This is bad software design, never assume default physical quantities that fundamentally change a simulation if nothing is specified. Or, if you want to, make that clear in the config file by setting a flag like useDefault = true. Now the user can look up what these defaults are and have a chance to catch physics bugs.

Sin #3: Too many positional arguments, which change from release to release

This is an ANSYS Fluent-only issue, as far as I can tell, and I hate that this is a feature. I find myself moving away more and more from Fluent, though I can't really avoid it as it is the university's default CFD solver, so I have to use it in teaching.

We can write journal files in Fluent, that is, text files that instruct the solver what to do, one line at a time. The following journal file is a valid file that will run in Fluent, assuming we have a mesh ready and it is called mesh.cas.h5.

1
2
3
4
/file/read-case "mesh.cas.h5"
/solve/initialize initialize-flow
/solve/iterate 100
/file/write-data "solution.dat.h5"

It probably won't do anything, but notice the absence of any solution-specific inputs? No boundary conditions were specified, no viscosity, nothing. Fluent will happily run your simulation, it will assume defaults for everything unless you change it.

But that is not my point. Let's say you actually want to change the material and add water as a liquid. Then you add the following command:

1
2
/define/materials/copy fluid water-liquid
/define/boundary-conditions/fluid water yes water-liquid n n n n 0 n 0 n 0 n 0 n 0 n 0 n n n n n

Do you have any idea what the 22 positional arguments do? No? Neither do I. And you know what, there is also no point in learning what they do because Fluent is dropping two releases per year, and they constantly change the arguments when they add a feature. There is no backwards compatibility.

Now, imagine what happens in practice. You have crafted a journal file that works with your current Fluent release. Then you update to the next release. Your journal file now uses an interface that is no longer supported and it will break. If you got an error message, it would be annoying, but you would know at least that something is wrong. But this is not what Fluent does. It will look at your command, see that it is incorrect, print an error message, and then continue using the default values.

This error message then gets buried in a long log file, and you can later inspect it yourself if you notice something went wrong. Because this is not an acceptable default behaviour, Fluent offers its very own hack. You can specify the version your journal file was set up with so that a newer version of Fluent can map your old command to the latest interface.

Sounds like a good solution at first (I am sure the engineers at ANSYS are still tapping themselves on the shoulders for a job well-done), but they violated the open-closed principle with their interface design and thus deserve to burn in hell. The open-closed principle, in object-orientated programming, states that a class should be open for extension but closed for modification.

Open for extension means that we should be able to add functionality to a class, like adding a new method (function) to it. Closed for modifications means that whatever changes I do should not affect the existing class design. This is very abstract, let me demonstrate with a simple example, we don't even need classes for that, just a simple function.

Let's imagine we want to compute the time step from a CFL number (and, for simplicity, I'll use Python here), we may reasonably do something like:

1
2
3
# version 1
def compute_dt(CFL, velocity, spacing):
    return CFL * spacing / velocity

Later, we realise that the computed time steps may, on some occasions, become too large, or too small. So we may want to clamp them to a min and max time step size. So we modify our function to:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# version 2
def compute_dt(CFL, velocity, spacing):
    dt = CFL * spacing / velocity
    max_dt = 1
    min_dt = 1e-6

    if dt > max_dt:
        dt = max_dt
    elif dt < min_dt:
        dt = min_dt

    return dt 

So far, so good. You may get an uncomfortable code review with a senior developer for not packaging lines 7-10 into its own function (essentially we have re-invented the clamp() function), but it could pass. So far, anyone who has written their scripts using version 1 can still execute version 2 as well. Then, we have a round table discussion, you proudly present your function as the implementation of the month, but get harsh comments for defining min_dt and max_dt in code, that should be something the user sets. So you write:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# version 3
def compute_dt(CFL, velocity, spacing, min_dt = 1e-6, max_dt = 1):
    dt = CFL * spacing / velocity

    if dt > max_dt:
        dt = max_dt
    elif dt < min_dt:
        dt = min_dt

    return dt 

The function signature grows, but since we specify default values here for min_dt and max_dt, we can still call version 3 as we did in version 1. So, we have rewritten the function now twice already (we modified it twice), but both times the user was able to use the same function signature and gained our modifications for free. They did not have to change anything on their end. We are open for extending our function's functionality, but closed to changing the interface we use to call this function.

Let's break this. Later we realise that the above definition only works for high Reynolds number flows, and we really ought to check the viscosity as well, so we differentiate now between the inviscid and viscous CFL number, where the viscous CFL number may dominate for low Reynolds number simulations, especially in the near-wall region. We now need the viscosity as well, and so we write:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# version 4
def compute_dt(CFL, velocity, spacing, nu, min_dt = 1e-6, max_dt = 1):
    inviscid_dt = CFL * spacing / velocity
    viscous_dt  = CFL * spacing * spacing / nu

    dt = min(inviscid_dt, viscous_dt)

    if dt > max_dt:
        dt = max_dt
    elif dt < min_dt:
        dt = min_dt

    return dt 

Now we have introduced nu as a fourth positional argument, without a default value, meaning that anyone calling this function now with only 3 arguments will get an error message. A sensible software engineer would stop here and examine the interface to see how this can be fixed. The fine people at ANSYS decided to show us the finger and came up with this implementation instead:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
# version 5
def compute_dt_implementation(CFL, velocity, spacing, nu, min_dt = 1e-6, max_dt = 1):
    inviscid_dt = CFL * spacing / velocity
    viscous_dt  = CFL * spacing * spacing / nu

    dt = min(inviscid_dt, viscous_dt)

    if dt > max_dt:
        dt = max_dt
    elif dt < min_dt:
        dt = min_dt

    return dt

def compute_dt(*args):
    try:
        return compute_dt_implementation(*args)
    except:
        return 3.141592

First, we rename the function so it stays a hidden implementation secret. Then, we create a new function that users were already calling before, and we no longer care what arguments the users provide. We don't even sanity check them. It is the software engineering equivalent of calling the customer service department and then waiting for an hour while repeatedly hearing the words "Your call is very important to us, please hold".

We then try to call the original function and just pass all the arguments the user provided us, you know, the arguments that are very important to us, so important that we can't even bother checking if the arguments are correct, or make sense. But, we know there is a decent chance that the function call will fail now, and so, to make sure a simulation never stops but always proceeds, we return a magic number.

It can be anything, the simulation doesn't have to be right, or the user doesn't need to know their time step computation failed. As long as the simulation progresses, users are happy. While you could argue this is still adhering to the open-closed principle (we can still call the function we wanted, and it does something), we have now turned a working function into a broken promise. The function no longer works as expected. There are no points for being technically correct when the end result is unusable.

So, we have already opened the door for the devil, why stop there? If we really want to implement this damn time step computation properly, we also need to differentiate between the inviscid and viscous CFL number. Each time step ought to be computed at its own correct CFL number, so we may extend this now as:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
# version 7
def compute_dt_implementation(inviscid_CFL, viscous_CFL, velocity, spacing, nu, min_dt = 1e-6, max_dt = 1):
    inviscid_dt = inviscid_CFL * spacing / velocity
    viscous_dt  = viscous_CFL * spacing * spacing / nu

    dt = min(inviscid_dt, viscous_dt)

    if dt > max_dt:
        dt = max_dt
    elif dt < min_dt:
        dt = min_dt

    return dt

def compute_dt(*args):
    try:
        return compute_dt_implementation(*args)
    except:
        return 2.71

Great, so now we have two different CFL numbers, we also changed the default time step that is returned without telling anyone because we can (and so we do), and also we jumped from version 5 to version 7, but that is an HR issue, as we just hired someone from the HDF5 steering committee, they are somewhat challenged when it comes to version numbering. A throwback to the HDF5 article ...

In version 7.7, we may then also remove the defaults for min_dt and max_dt, rearrange the positional arguments, or introduce additional time step computations (for example, based on a period, where we subdivide a period into N segments (think: vortex shedding and how many time steps we want to evaluate for the dominant frequency)). Our original compute_dt function may no longer exist.

You may think this is a bit over the top, or for comical effect, but I am being serious this time. What I have described above is the story of every ANSYS Fluent release. You know the saying "Nothing can be certain, except death and taxes"? We shall, from now on, use the phrase:

Nothing can be certain, except Fluent breaking the interface in their next release

I am being harsh here, yes, but for the past 5 years or so I have been battling with all of these issues and have become somewhat of a Fluent debugging expert. Don't ask ChatGPT, I have seen it all! But then again, please don't ask me, I am really sick and tired of their version incompatibility issues. Fluent used to follow semantic versioning until 2019 (I believe) when they switched to calendar versioning.

Semantic versioning communicates whether an interface breaks (as we have seen with our compute_dt() function). So, when we introduced nu as an additional parameter, we broke the interface. Semantic versioning communicates that explicitly, i.e. if we had, say, version 12.1.0 and then released 13.0.0, it communicates "We introduced breaking changes, use with caution".

Fluent used semantic versioning, but with their practice of breaking everything, 2 a year, we would be now on version 50.0.0 I reckon, which is a tad silly. So, for those software engineers who are a bit mentally challenged, we have introduced calendar versioning, where we simply slap on the year, and possibly a release number after it, like 2026.1 and 2026.2. It states nothing about backwards compatibility. Calendar versioning itself states:

CalVer is a versioning convention based on your project's release calendar. Versioning gets better with time.

Sure, and Western Union is a service exclusively used to send money to friends worldwide. It is not used by scammers to receive payments from victims worldwide.

I think you get my point. Enough Fluent for the day, I only just got up an hour ago and my watch is already telling me to do some breathing exercises ...

Sin #4: Consuming proprietary file formats like petrol station sushi

Sure, I enjoy sushi as much as the next guy (that is, not at all). Honestly, raw fish? I mean, just imagine German sushi: a raw potato. Or Italian sushi: uncooked pasta. Or British sushi: a raw (Cadbury cream-) egg. Would you eat that? Probably not, so why are you still eating food that someone forgot to prepare? Actually, I do like British sushi, thinking about it ... Ah well, I know what I'll be getting from Tesco Express for breakfast ... A nice scrambled Cadbury Egg. Delicious!

But going back to the analogy, let's say you actually do enjoy sushi, would you, when going to the gas station at midnight on a Sunday evening buy the last pack of sushi that is on sale?

CFD solvers, as well as pre-processors (meshing) and post-processors all use their own proprietary formats. ANSYS has their .cas.h5, Star-CCM+ their .sim, Tecplot their .plt, and the list goes on. As long as there is a common standard that everyone is implementing to the dot, there is no issue. But that isn't a realistic reality to believe in.

Every software will implement what it believes it needs. Not what it believes you need. Locking yourself into using only proprietary formats means you don't really own your own data. For example, imagine you buy a Fluent licence for a single project, but only for a year. You store all the output in Fluent's .dat.h5 format. Then, you open that in Tecplot, which you also got a licence for, for the year. You do your post-processing and everything is good. Then you return that licence as well.

A year later you need to check something in your data, but realise that you can only really read .dat.h5 files with Fluent or Tecplot. Paraview is open source, and that is what you turn to first, but it only supports the old .cas reader, and even that is flaky. You can convert Fluent data into a format that Paraview understands, but you need to modify the exported file itself. In any case, you need Fluent for that again.

This may be a very niche scenario, but one that could have been entirely avoided by simply adopting a more open and supported file format. CGNS is the obvious choice here, but again, Paraview doesn't like it that much. But with CGNS you do at least have the option of reading and converting it yourself, if you absolutely had to, so it is still a better long-term storage solution than storing your results in a proprietary format you lose access to the day you walk away from that licence.

The other issue is that you have no idea what is in the file. If there is a bug, and it is setting a solver setting you did not expect to be set, that could, in the worst case, lead to an incorrect simulation. Tracing that behaviour will become a nightmare, as it is pretty much impossible to figure that out on your own without the help of someone on the inside who is also willing to figure out whether there is a bug or not in the file format parsing.

The software community has embraced the open-source model, and it is a fantastic model; anyone can contribute, anyone can inspect source code, anyone can find issues on code and flag these. If we have a choice, we should always go for something that is easy to work with in the long term, even if software licences or workflow changes.

I'll stop now, I am as surprised as anyone that I managed to fill this article with that much content, apparently I just had a lot on my chest I needed to share.

Summary

If you have made it this far, respect! Let me check if you believe my first sentence to be correct. I started this article with:

Everyone believes reading input files is a simple task. Are you one of them? Look closer, and you realise that there are issues you haven't even thought about; issues that may come back to haunt you long after you have settled on a specific file type in your solver to use.

Do you agree now? And let's not forget, we've been looking at input files for the most part. I introduced data serialisation as a related concept (as it is mostly stored as JSON data, and we already talked a great bit about JSON), as well as to express my new love for Taylor (your choice which Taylor I mean, I love them all the same). We have left data sharing formats mostly untouched (e.g. to read/write meshes and solution data), as I have covered them already in great detail previously.

If you really sat through the entire article and read it word by word, you will have, inevitably, lost some sanity, and that is normal. When I started with this website, I wanted to cover the parts that never get any mention in CFD textbooks. You know the equations, and the numerical schemes, you may even know one or two things about parallelisation, and then you sit down and start typing your code. And then the floodgate of questions opens.

One question is how you get data from disk into your solver and yes, it is a somewhat trivial exercise, and it can be achieved with a simple INI file, but I think something as fundamental as the interface between your CFD solver and the human sitting in front of it and using it should get its own discussion, with a clear list of pros and cons, as well as failed designs from which we can learn.

Granted, file types aren't as sexy as, say, the sort of calendar your health and safety department would quietly confiscate (mine all show airplanes), but a wrong decision here will snowball into frustration later. I think there is no point in blowing things out of proportion (he said and then writes a 3-hour article on input files ...), so let's summarise it in the following way:

If you want to be kind to yourself, your fellow humans, and do some good in the world: use TOML. Not because I adore anyone whose name is Tom and has an obsessive need to let anyone know about themselves (that's not why I created this website), but because it is designed by a human for humans. And it shows in the design.

JSON and YAML or OK, you treat sanity here a bit for usability, but either choice is still far better than plain INI files. Of course, if you want to use YAML, you may want to look for a different parser than the one shown in this article. But there is no good alternative. If there is, please keep it to yourself. I think I am permanently damaged and just thinking about the YAML section gives me PTSD. We shall not speak of YAML around here anymore!

I can't change the past, or the reluctant professor with their 30-year-old legacy code, written in Fortran 66, which parse their own input file format, that is probably even more primitive than INI files and somehow still manages to be more buggy. Not that I speak of experience.

No, the target audience of this article are those who learn CFD fresh, who are willing to question if we really should be using INI for everything, or if there is, perhaps, a better alternative. Sure, you can prompt your favourite LLM, and it will give you a list of formats, and it gives you a list of pros and cons if you ask it, but what an LLM can never give you is experience.

I have worked with different solvers, different ways of setting up simulations. As the saying goes: I have seen things (and I can't unsee them!). So, this article is my contribution to let you learn from that shared experience. My hope with this article is not that we start to rewrite the ANSYS user interface, or make OpenFOAM accept non SI-units in their input file system (it won't happen), but instead, that we make informed decisions on how we want to present a human with a way to interact with our solvers.

You have the knowledge now, use it. And, of course, if you ever need to parse a CSV file, there should be no question, Exceln't is the way to go, of course!

Tom-Robin Teschner

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.

signature
gwanghwamun
Popup Image

Get started today!

Do you want to learn how to put the Navier-Stokes equations into code and write your own solver?

Sign up for my newsletter, and you get my eBook Write your first CFD solver - From theory to implemented CFD solver in less than a weekend for free!