class design for non-LTE physics modules with compile-time polymorphism

[BenWibking]

We want to be able to support a wide variety of non-LTE physics implementations that are problem-dependent.

An example of the "simplest" possible non-LTE model is a radiation model with independent gas, radiation, and dust temperatures. At the other extreme, one could imagine solving for the time-dependent non-equilibrium atomic level populations. Even for the latter, there is the choice of the number of atomic levels included in the model, which atoms to include, and what numerical method is used. There are easily thousands of combinatorial possibilies of non-LTE models users might want to use. A mature code might include several dozen distinct options intended to use across a variety of problem types.

In the long term, hard-coding these all of these models in a very long if or if constexpr statement is not a maintainable option. On the other hand, having the implementation of each model exist in the problem generator is not ideal -- although some users will want to, most users will not implement a new non-LTE for each new simulation type.

This requires a significant redesign of the class structure of RadiationSystem. The constraints on the design are:

1. It must be flexible enough for users to implement new models without changing the core radiation code.
2. It must allow compile-time optimizations so that only the radiation model that is used is compiled into the GPU kernel used for the radiation moment source terms, and
3. It must be maintainable (no spaghetti code, endless if statements, etc.).
4. It must not change the build process (no changes to CMake scripts)

One compelling option is to use the so-called "Factory" pattern in C++. A good introduction to this is this blog post. (This factory pattern is already used in Quokka for the in-situ diagnostics, such as the histogram and slice outputs. See src/Factory.H for how it works.)

The downside is that there is some boilerplate code that is not particular understandable to new code users, but it should never need to be changed once it is written. A more difficult issue is that we need to carefully design the Factory so that it allows for compile-time optimization of GPU kernels (i.e., only the user-selected non-LTE module should be compiled and run for a given simulation).

Are there other viable options for the class design that can accomplish the three goals listed above? For instance, is this possible to do with template specialization, like we do for problem_t? (For example, we currently do this when choosing the Riemann solver based on the options selected in Physics_Traits, but we have an explicit enumeration of all of the Riemann solvers with if constexpr for each option.) The issue I do not have an immediate solution for when using a pure template-based solution is how to avoid having a large if constexpr list of each radiation model somewhere in RadiationSystem.

cc @chongchonghe @markkrumholz

[markkrumholz]

I guess am I a bit less worried than you about the constexpr solution, partly because I do not see a realistic possibility that we are going to wind up with more than a handful of these implementations, mainly because I don't think all of the possible variations on implementations wind up in a single monolithic if block that will be unmaintainable. For example, if one wanted to implement non-LTE level populations, I imagine that the pattern would be a single constexpr choice for this, and that choices of numbers of atoms and levels would then be handled internal to a general non-LTE level populations model. That complexity won't all end up in RadiationSystem. So the way I see it we are probably looking at no more than 5 or so distinct treatments, which is fine for an if block.

[BenWibking]

Five if/else statements would be perfectly fine, agreed. If it gets to be much more than this, we can revisit this.