Regression Tests Not Portable Between Different Compilers and Systems

[jtramm]

The Problem

The regression testing strategy in OpenMC compares results for each test against previously generated result files. This strategy is simple and highly effective at ensuring PRs don't break working code. However, the downside is that it fundamentally assumes that different compilers and systems will produce the same outputs given the same version of the code. I've found that I often get test failures (e.g., different test results) when compiling on both my local macbook as well as linux cluster systems as compared to what OpenMC's github CI produces. Ideally - we would like our tests to run cleanly everywhere.

Causes

We can get different results from OpenMC on different compilers and systems for at least three reasons:

1. Floating point non-associativity. This is not a big deal for serial operations, but when running with multiple OpenMP threads (or sometimes in MPI depending on if the code requires strict rank ordering), certain floating point operations (like tallying) can happen in inconsistent order depending on which system conditions causing threads to execute in different orders. This problem isn't typically very serious though, as it only shows up in a few places, and usually doesn't cause macro-scale divergence.

2. Differences in system math transcendental operations (e.g., std::sin(), std::exp()) etc. In the C++ regions of OpenMC, these functions map to the system's C standard library implementation (libm). On linux systems, this is typically the GNU C Library (glibc), on macs this is the apple libSystem (e.g., libSystem.B.dylib), and on windows this is provided by the microsoft C runtime. All of these libraries have independent implementations. Furthermore - transcendental functions are not governed by IEEE 754 rules, such that each implementation may be rounded differently. As such, something like std::exp(-0.1) may give a slightly different result, even just when running on linux but using two different versions of glibc.

It's easy to see how differences in log() might affect transport and cause macro-scale divergence over time (e.g., in Particle::event_advance()):

    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;

Similarly for scattering operations, sin/cos etc will adjust outgoing angles slightly, causing a particle that previously nicked a surface to miss it, and causing macro-scale divergence.

3. Differences in compiler FMA contraction. Modern CPUs and GPUs typically have specialized floating point instructions for fused multiply-add (FMA) operations. Consider a line of code like a = b + c*d. Normally a compiler would treat this as a multiply operation and then an addition operation, with an implicit rounding operation happening after each operation (i.e., two FP operations + 2 FP roundings). This can be "contracted" into a single FMA instruction. The FMA instruction is notable as it is just one FP operation + 1 FP rounding. As such, it can be more accurate due to only having to round once at the end, which can be quite valuable in some cases. Additionally, depending on the architecture, it may be faster as 1 FMA might just take 1 clock cycle as compared to needing 2 for two separate multiple then add instructions.

Compilers have different default configurations for how they decide whether or not to "contract" a candidate operation into an FMA operation. By default they may have different levels of aggressiveness. Even under the same flags, it's potentially possible that one may limit its analysis to single line contractions, whereas another might do more analysis to determine if two lines of C++ code might be contracted into a single FMA. As such, when any level of FMA contraction is enabled, you can potentially get different instructions being generated between compilers and compiler versions, even on the same system. As FMA operations affect rounding, the results of a single line of code like a = b + c*d may give different results for a.

This problem is notable in that it occurs right out of the gate in OpenMC. For example, in material.cpp in Material::normalize_density(), as we are initializing material data right after reading it in, we are already going to be rounding things differently depending on FMA treatment:

  for (int i = 0; i < nuclide_.size(); ++i) {
    int i_nuc = nuclide_[i];
    double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_ : 1.0;
    int z = settings::run_CE ? data::nuclides[i_nuc]->Z_ : 0.0;
    density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO;  # Potential FMA
    charge_density_ += atom_density_(i) * z; # Potential FMA
  }

Thus, before we've even sampled the first particle, our underlying material/XS data is a little different. It won't take long before you get macro-scale particle divergence.

4. Misc. There are a few other potential differences between compilers that might cause problems, though these all seem less severe than what is listed above. For instance, some systems might treat subnormal floating point values different (e.g., flush to zero). Some compilers might auto-vectorize a reduction loop, in which case the floating point ordering may change significantly.

Implications

The implications of the above three issues are that you may see both minor differences in simulation outputs (e.g., tally results) caused by minor changes in FP rounding. By itself, such small changes in rounding are not such a big deal, as we can make tests tolerant of tiny rounding differences. However, the bigger issues at play is that tiny differences in XS data, sampling, and scatter operations can rapidly cause macro-scale divergence. A particle that previously just barely nicks a curved surface may instead miss it, potentially causing it to sample its next reaction in a different material altogether, and thus causing it to scatter to a completely different angle or sample a totally different sort of reaction. This can cause more significant changes in tally results for fixed source problems (especially true if we are only running a few hundred particles, as is typically the case for fast regression testing). Furthermore, it can cause complete divergence in eigenvalue mode, where a minor change in one history in the first batch can affect which particles get sampled in the next batch, causing the simulations to be be running completely different particles after just a few batches.

Severity and Potential Fixes

Each of the problems listed in the "Causes" section above have different severities and potential remedies.

1. Parallel Floating point non-associativity. For most cases I can think of I don't see this causing macro-scale divergence, rather just minor rounding differences in tallies. As such, I don't see that we are sensitive to these types of differences yet, as our testing framework seems pretty reproducible (at least for a given compiler + system) regardless of how many threads are used. Thus, no remedy may even be needed.

2. Math Transcendentals. This is a much higher severity problem, and I believe does cause many test failures. Different linux ssytems may have different glibc versions. MacOS systems will use the apple standard library instead of the GNU C standard library as on linux, etc. Thus, the only way I can see to get around this is to use a fixed implementation of the transcendental functions instead of using the system library. There are actually good options here, like the core-math https://gitlab.inria.fr/core-math/core-math/ library. In theory this could be ingested as a vendor dependency. Or - as it is distributed under the MIT license, we could ingest it permanently into the OpenMC repo, which would make things like adding GPU flags much easier.

3. FMA Contraction. This is also a high severity problem, which is most noticeable when compiling with LLVM clang vs. GNU gcc. Thankfully, there are compiler flags that can globally disable/enable/set these policies such that consistent behavior can be achieved by simply altering the OpenMC CMakeLists.txt to be more explicit about what we want, and not leave it up to default CMake/compiler-specific choice.

4. Misc. I'm guessing these are all low severity and we aren't affected by this, but we may need to look more closely at these areas if the above items are all remedied and differences are still showing up in tests.

[jtramm]

So as seen in #3823, the fix for issue (3) about FMA differences is quite easy. The remaining major issue (2) above regarding differences in math transcendental functions still causes about 14 test failures on my mac, due to use of the apple math library instead of the GNU C math library glibc. It's entirely possible we may see a similar scale of divergence if we ran the tests on GPU (once those capabilities are in).

There may be other ways of making a fixe here, but to me, the most obvious way of fixing this is to use a consistent set of math operations in OpenMC instead of calling the system library ones (e.g., std::cos(), std::log(), etc.). There's basically two options for this that I can see:

1. Add the core-math library (https://gitlab.inria.fr/core-math/core-math/ ) (or something similar?) as a dependency. This implements correctly rounded operations for all of the math functions we need, and seems well validated (even moreso than system library!) so would be a reliable choice. The downside is that we're adding another submodule, and dealing with future dependency rot, possible API changes, inconsistent versions between submodule and something floating around on people's systems, etc. The impression I get is also that core-math ensures that rounding is correct, though at the cost of performance and larger code size etc. Comparatively, the apple and GNU libraries tend to balance accuracy with speed, with more of an emphasis on speed. Another potential problem is that, in looking through the code, the code is not very human readable, with a lot of platform specific IFDEFs, platform-specific intrinsics, etc. It's about 60k lines of this. If we go the submodule route, we would also potentially not be able to use the library on GPU, as we may not be able to easily apply #pragma omp device around function definitions, and there are also some 128 bit types and intrinsics that are common on CPU but not supported on GPU.

So: I think this route would work to close the gap between linux vs. mac, but may not be extensible to GPU. Ideally we'd pick a solution that would help universally.

2. We "eat" the core-math library (distributed under the MIT license, so this wouldn't be a problem). We can then strip out the functions we don't need (going from about 60k lines down to 20k lines) and simplify the rest (further reducing things from 20k lines down to about 6k lines) to be fully portable. I've actually tried this (branch is at: https://github.com/jtramm/openmc/tree/core-math-transcendentals) and it does appear to fix almost all the tests. The exception is that three of the tests that compile code (like the .cpp driver test), or use dlopen, don't seem to work still on macos, but that's a separate issue.

I'm curious what people think about option (2) above. On the one hand, it seems nuts to make fine-grained edits on 16 different math functions and then maintain them ourselves. However -- they are all derived from a really reputable library (core-math) (so we aren't cooking up all the methods for them from scratch or anything). Additionally, they are extremely easy to test empirically by feeding randomized inputs and comparing outputs against the std:: functions. In that light, there doesn't appear to be many opportunities to go wrong as long as we're able to empirically demonstrate accuracy against reference implementations. The AI agent was able to simplify the functions enough to be portable, and strip out a lot of the really expensive code that is required to adjudicate the "tie breaker" scenarios governing rounding of the last bit. In the first draft, measuring our simplified core-math library vs. glibc (what you would typically get on linux), in terms of ulp we get:

Function	Max ulp error
exp	1
expm1	1
log	1
log1p	1
sin	1
cos	2
asin	1
acos	1
atan	2
atan2	1
sinh	2
erf	1
erfc	4
lgamma	7
tgamma	7
pow	1

which seems pretty solid, given we are looking more for reproducibility than perfection. We could alternatively work towards keeping the accuracy of the original core-math library but just trying to make it more portable, if more accuracy seems useful. In any case, I would want to do a lot more analysis here before proposing a PR with this, but won't bother to dig any deeper if people aren't receptive to the idea of having our own math functions.

[HunterBelanger]

I would agree that this is definitely a problem, and it is one of the reasons that I never finished my PR on Windows support in #2919. I am fairly confident that OpenMC was working fine on Windows, but there were several tests that would fail, and I think it was due to floating point differences on Windows, leading to different results than what we had in the reference files.

I don't get to contribute to OpenMC enough now to have an opinion that matters much, but I agree with your statement that it is insane to maintain our own math function implementation. However, if it's really only 6k LOC and will work on CPU and GPU at that point, I think it could be highly valuable and potentially worth the headache / maintenance.