.. _sec-developer: Developer manual ================ Developer tips -------------- Compiling Hermes-3 quickly ~~~~~~~~~~~~~~~~~~~~~~~~~~ After compiling Hermes-3 and changing something, you can avoid unnecessary recompilation of unchanged files. Simply enter the build directory and do: .. code-block:: bash make -j 4 Where ``4`` will make it compile using 4 cores. You can increase or decrease this as necessary. Sometimes, CMake will think that it needs to recompile parts of BOUT++ even if it wasn't changed. To stop this try adding the flag ``-DHERMES_UPDATE_GIT_SUBMODULE=OFF`` in your CMake command. This and other useful options can be found in ``hermes-3/CMakeLists.txt``. Compiling documentation ~~~~~~~~~~~~~~~~~~~~~~~ The Hermes-3 documentation is built using `Sphinx `__ and `Doxygen `__. It's written in `ReStructuredText (RST) `__, which is a markup language similar to Markdown. Doxygen generates automatic documentation based on the C++ code, while Sphinx handles everything else. Editing documentation is much easier if you can compile it locally using the following steps: 1. **Install Sphinx** and our theme in your Python environment: .. code-block:: bash pip install sphinx sphinx_book_theme 2. **Install Doxygen** (modify as necessary for your OS) and Breathe, the package that connects it to Sphinx: .. code-block:: bash sudo apt install doxygen 3. **Install Breathe** (modify as necessary for your OS): .. code-block:: bash pip install breathe 4. **Run Doxygen** - this will parse the C++ code: .. code-block:: bash cd hermes-3/docs/doxygen doxygen Doxyfile 5. **Run Sphinx** - this will parse the RST files and generate the documentation. ``sphinx`` and ``build`` are the source and build directories, respectively. .. code-block:: bash cd hermes-3/docs sphinx-build sphinx build 6. **Open the generated HTML files**, either by double clicking on the file in your browser, or some other way. If you use VS Code locally or on a remote machine through SSH, you can use the extension `Live Preview `__ which can stream it to your browser. Debugging: running for one iteration ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Any BOUT++ code can be run for just one right-hand side (RHS) iteration. This will run instantly for any simulation and not need any kind of solver convergence, making it ideal for debugging. This can be done by setting the following in the input file: .. code-block:: ini nout = 0 The ``timestep`` setting will be ignored. Debugging: printing values ~~~~~~~~~~~~~~~~~~~~~~~~~~ When debugging, it can be useful to print things out. This is simple in C++. For example, if you want to print the value of the variable ``particle_flow``, do: .. code-block:: ini output << "\n*******************************\n" output << "particle_flow: " << particle_flow << "\n" output << "*******************************\n" Which will result in: .. code-block:: ini ******************************* particle_flow: 0.123456 ******************************* Sometimes you need to deal with printing whole fields. In this case, it is useful to print out every cell value along with its coordinates. The below code snippet will do this for the variable ``Tn``: .. code-block:: ini or(int ix=0; ix < mesh->LocalNx ; ix++){ for(int iy=0; iy < mesh->LocalNy ; iy++){ for(int iz=0; iz < mesh->LocalNz; iz++){ std::string string_count = std::string("(") + std::to_string(ix) + std::string(",") + std::to_string(iy)+ std::string(",") + std::to_string(iz) + std::string(")"); output << string_count + std::string(": ") + std::to_string(Tn(ix,iy,iz)) + std::string("; "); } } output << "\n"; } The output will look something like: .. code-block:: bash (0,0,0): 0.123456; (0,0,1): 0.654321; (0,0,2): 0.987654; (0,1,0): 0.123456; (0,1,1): 0.654321; (0,1,2): 0.987654; (1,0,0): 0.123456; (1,0,1): 0.654321; (1,0,2): 0.987654; (1,1,0): 0.123456; (1,1,1): 0.654321; (1,1,2): 0.987654; ... Note that there are multiple ways to print out values in C++, and both ``std::string(text)`` and ``"text"`` are valid. Debugging: segmentation faults ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ `Segmentation faults `__ can be frustrating because they give very little verbosity. In practice, the most common cause is trying to access a variable that hasn't been initialised yet. The easiest way to debug this is to carefully review the new lines of code to make sure all variables exist and have been declared and initialised. If this is tricky, another simple way is to comment out large parts of the code until the segmentation fault disappears, helping to narrow down its location. While the above methods are very simple and can be effective, debugging tools such as `gdb `__ or `valgrind `__ can be used to find the segmentation fault as well. Debugging: compiling in debug mode ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ This can be a useful way to catch errors. Please see the `relevant page `__ in the BOUT++ documentation. Header vs. implementation files ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ allows you to split code into implementation (.cxx) and header (.hxx) files. The convention of what should be in each one is not consistent in Hermes-3 at the moment. The most common standard is for the header file to contain all of the variable and class declarations and the implementation file to contain the rest. An example of this is in the ``evolve_density`` component - see the `implementation file `__ and the `header file `__. Data types ~~~~~~~~~~~~~~ Hermes-3 casts its variables in a variety of BOUT++ classes. Floats are usually represented as ``BoutReal``, and fields as ``Field3D``. Note that Hermes-3 always runs "in 3D" - when configured in 1D, the x and z dimensions are of unit length. See relevant `BOUT++ docs `__ for more info. There is also a data type called ``Options`` which is equivalent to a Python dictionary with extra functionality, and is used to store input options, the entire simulation state and many other data. Finally, there is the ``GuardedOptions`` datatype, which wraps an ``Options`` object and controls access to its contents. Adding new settings ~~~~~~~~~~~~~~~~~~~ This is simple and uses the following syntax: .. code-block:: ini bndry_flux = options["bndry_flux"] .doc("Allow flows through radial boundaries") .withDefault(true); See `this real world example `__. The variable must also be declared in the corresponding header file. Adding new diagnostics ~~~~~~~~~~~~~~~~~~~~~~ Adding new diagnostics is also simple, provided there is already an ``outputVars`` function set up in your component. This is usually located at the end of the implementation file. See this example from ``evolve_momentum.cxx``: .. code-block:: ini void EvolveMomentum::outputVars(Options &state) { AUTO_TRACE(); // Normalisations auto Nnorm = get(state["Nnorm"]); auto Omega_ci = get(state["Omega_ci"]); auto Cs0 = get(state["Cs0"]); state[std::string("NV") + name].setAttributes( {{"time_dimension", "t"}, {"units", "kg / m^2 / s"}, {"conversion", SI::Mp * Nnorm * Cs0}, {"standard_name", "momentum"}, {"long_name", name + " parallel momentum"}, {"species", name}, {"source", "evolve_momentum"}}); Each diagnostic is saved with metadata, which is available in xHermes, e.g. .. code-block:: ini ds["Nd"].attrs() time_dimension This indicates that the variable is time evolving. If you remove this line, it will be saved only at the first RHS evaluation. units A string showing the units for post-processing. xHermes picks this up. conversion A float representing the normalisation factor. xHermes picks this up to do automatic conversion to SI units. standard_name, long_name We aren't consistent on what should be in each, but they are meant to describe the variables in post-processing. species, source The relevant species and component that the diagnostic is coming from What is "Options"? ~~~~~~~~~~~~~~~~~~ Options is a dictionary-like class originally developed for parsing BOUT++ options. In Hermes-3, it is used as a general purpose dictionary. Getting/setting values ~~~~~~~~~~~~~~~~~~~~~~ Hermes-3 has a system to prevent quantities from being modified after they are used. This is important as it uses a single dictionary-like ``state`` class to hold all of the variables in one place, which could allow some components to overwrite others. In ``component.hxx`` there is the function ``get``, which once called sets the "final" and "final-domain" attributes: .. code-block:: ini T get(const Options& option, const std::string& location = "") { #if CHECKLEVEL >= 1 // Mark option as final, both inside the domain and the boundary const_cast(option).attributes["final"] = location; const_cast(option).attributes["final-domain"] = location; #endif return getNonFinal(option); } When you call ``set``, these attributes are checked for, so that if they have already been "gotten", they can't be set again: .. code-block:: ini Options& set(Options& option, T value) { // Check that the value has not already been used #if CHECKLEVEL >= 1 if (option.hasAttribute("final")) { throw BoutException("Setting value of {} but it has already been used in {}.", option.name(), option.attributes["final"].as()); } if (option.hasAttribute("final-domain")) { throw BoutException("Setting value of {} but it has already been used in {}.", option.name(), option.attributes["final-domain"].as()); } if (hermesDataInvalid(value)) { throw BoutException("Setting invalid value for '{}'", option.str()); } #endif There is a special use case which allows you to use this "locking" scheme for only the domain cells, leaving the guard cells to be settable using ``getNoBoundary``: .. code-block:: ini T getNoBoundary(const Options& option, const std::string& location = "") { #if CHECKLEVEL >= 1 // Mark option as final inside the domain const_cast(option).attributes["final-domain"] = location; #endif return getNonFinal(option); } And there is a corresponding ``setBoundary`` that can be used for BC operations: .. code-block:: ini Options& setBoundary(Options& option, T value) { // Check that the value has not already been used #if CHECKLEVEL >= 1 if (option.hasAttribute("final")) { throw BoutException("Setting boundary of {} but it has already been used in {}.", option.name(), option.attributes["final"].as()); } #endif option.force(std::move(value)); return option; } All of these functions are overloaded to accept both `Options` and `GuardedOptions` objects. These functions take a second argument which tells you where they were set, which is easier for debugging. They are wrapped into additional functions, ``GET_VALUE`` and ``GET_NOBOUNDARY`` which automatically include this argument. Please review `component.hxx `__ for more details. Looping over cells ~~~~~~~~~~~~~~~~~~ BOUT++ provides a really easy way to loop over the domain using ``BOUT_FOR`` and similar loops, see `BOUT++ docs `__. There is a way to way to tell if you are in the core or not. The ``mesh`` object has a function to indicate if the coordinate is in a periodic region or not. Only the core is periodic. See below for an example from ``evolve_pressure.cxx`` which makes sure a pressure source is set to zero outside of the core: .. code-block:: ini if (p_options["source_only_in_core"] .doc("Zero the source outside the closed field-line region?") .withDefault(false)) { for (int x = mesh->xstart; x <= mesh->xend; x++) { if (!mesh->periodicY(x)) { // Not periodic, so not in core for (int y = mesh->ystart; y <= mesh->yend; y++) { for (int z = mesh->zstart; z <= mesh->zend; z++) { source(x, y, z) = 0.0; } } } } .. _sec-code_structure: Code structure -------------- A hermes-3 model, like all `BOUT++ models `__, is an implementation of a set of Ordinary Differential Equations (ODEs). The time integration solver drives the simulation, calling the `Hermes::rhs` function to calculate the time-derivatives of all the evolving variables. The calculation of the time derivatives is coordinated by passing a state object between components. The state is a nested tree, and can have values inserted and retrieved by the components. The components are created and then run by a scheduler, based on settings in the input (BOUT.inp) file. For example a transport simulation with deuterium and tritium ions and atoms has an input file specifying the components .. code-block:: ini [hermes] components = (d+, d, t+, t, e, braginskii_collisions, braginskii_friction, braginskii_heat_exchange, sheath_boundary, recycling, reactions, braginskii_conduction) The governing equations for each species are specified e.g. .. code-block:: ini [d+] type = evolve_density, evolve_momentum, evolve_pressure, anomalous_diffusion AA = 2 # Atomic mass charge = 1 and other components have their configuration options e.g. for reactions: .. code-block:: ini [reactions] type = ( d + e -> d+ + 2e, # Deuterium ionisation t + e -> t+ + 2e, # Tritium ionisation ) In terms of design patterns, the method used here is essentially a combination of the `Encapsulate Context `__ and `Command `__ patterns. Simulation state ~~~~~~~~~~~~~~~~ The simulation state is passed between components, and is a tree of objects (Options objects). At the start of each iteration (rhs call) a new state is created and contains: * `time` BoutReal, the current simulation time * `linear` bool. True if the time integrator expects a linear response. * `units` * `seconds` Multiply by this to get units of seconds * `eV` Temperature normalisation * `Tesla` Magnetic field normalisation * `meters` Length normalisation * `inv_meters_cubed` Density normalisation so the temperature normalisation can be extracted using:: BoutReal Tnorm = state["units"]["eV"]; As the components of a model are run, they set, modify and use values stored in this state. To ensure that components use consistent names for their input and output variables, a set of conventions are used for new variables which are added to the state: * `species` Plasma species * `e` Electron species * `species1` Example "h", "he+2" * `AA` Atomic mass, proton = 1 * `charge` Charge, in units of proton charge (i.e. electron=-1) * `density` * `momentum` Parallel momentum * `pressure` * `velocity` Parallel velocity * `temperature` * `collision_frequency` Normalised collision frequency * `density_source` Normalised particle source * `momentum_source` Normalised momentum source * `energy_source` Normalised energy source * `particle_flow_xlow` Normalised particle flow through lower X cell face * `particle_flow_ylow` Normalised particle flow through lower Y cell face * `momentum_flow_xlow` Normalised momentum flow through lower X cell face * `momentum_flow_ylow` Normalised momentum flow through lower Y cell face * `energy_flow_xlow` Normalised energy flow through lower X cell face * `energy_flow_ylow` Normalised energy flow through lower Y cell face * `fields` * `vorticity` * `phi` Electrostatic potential * `Apar` Electromagnetic potential b dot A in induction terms * `Apar_flutter` The electromagnetic potential (b dot A) in flutter terms * `DivJdia` Divergence of diamagnetic current * `DivJcol` Divergence of collisional current * `DivJextra` Divergence of current, including 2D parallel current closures. Not including diamagnetic, parallel current due to flows, or polarisation currents For example to get the electron density:: Field3D ne = state["species"]["e"]["density"]; This way of extracting values from the state will print the value to the log file, and is intended mainly for initialisation. In `Component::transform` and `Component::finally` functions which run frequently, faster access methods are used which don't print to the log. To get a value:: Field3D ne = get(state["species"]["e"]["density"]); If the value isn't set, or can't be converted to the given type, then a `BoutException` will be thrown. To set a value in the state, there is the `set` function:: set(state["species"]["h"]["density"], ne); A common need is to add or subtract values from fields, such as density sources:: add(state["species"]["h"]["density_source"], recombination_rate); subtract(state["species"]["h+"]["density_source"], recombination_rate); Notes: - When checking if a subsection exists, use ``option.isSection``, since ``option.isSet`` is false if it is a section and not a value. - The species name convention is that the charge state is last, after the ``+`` or ``-`` sign: ``n2+`` is a singly charged nitrogen molecule, while ``n+2`` is a +2 charged nitrogen atom. Components ~~~~~~~~~~~~~~ The basic building block of all Hermes-3 models is the `Component`. This defines an interface to a class which takes a state (a tree of dictionaries/maps) and transforms (modifies) it. This is done by calling the public `Component::transform` method. This will call the private `Component::transform_impl` method, which must be overriden for each Component implementation. After all components have modified the state in turn, all components may then implement a `finally` method to take the final state but not modify it. This allows two components to depend on each other, but makes debugging and testing easier by limiting the places where the state can be modified. .. doxygenstruct:: Component :members: :protected-members: :private-members: Components are usually defined in separate files; sometimes multiple components in one file if they are small and related to each other (e.g. atomic rates for the same species). To be able to create components, they need to be registered in the factory. This is done in the header file using a code like:: #include "component.hxx" struct MyComponent : public Component { MyComponent(const std::string &name, Options &options, Solver *solver); ... }; namespace { RegisterComponent registercomponentmine("mycomponent"); } where `MyComponent` is the component class, and "mycomponent" is the name that can be used in the BOUT.inp settings file to create a component of this type. Note that the name can be any string except it can't contain commas or brackets, and shouldn't start or end with whitespace. Inputs to the component constructors are: * `name` * `alloptions` * `solver` The `name` is a string labelling the instance. The `alloptions` tree contains at least: * ``alloptions[name]`` options for this instance * ``alloptions['units']`` Component Permissions ````````````````````` All component constructors must pass a `Permissions` object (see below) to the constructor on the `Component::Component` base class. This specifies which variables will be read/written by the `Component::transform` method and will be used to construct a `GuardedOptions` object to be passed into `Component::transform_impl`. The `Permissions` object can be further updated in the body of the constructor of your component using the `Component::setPermissions` and `Component::substitutePermissions` methods. You should give read and write permissions to the minimum number of variables necessary, to avoid circular dependencies arising among components. A number of substitutions will automatically be performed on your permissions (see `Permission Substitution`_), so that you can specify permissions for some variables for each species. For example, the following permissions would give read access to pressure for all species and density of ions:: MyComponent::MyComponent(const std::string &name, Options &options, Solver *solver) : Component({readOnly("species:{all_species}:pressure"), readOnly("species:{ions}:density")}) {} See the documentation for `Component::declareAllSpecies` for a list of all substitutions that will be performed. Component scheduler ~~~~~~~~~~~~~~~~~~~ The simulation model is created in `Hermes::init` by a call to the `ComponentScheduler`:: scheduler = ComponentScheduler::create(options, Options::root(), solver); and then in `Hermes::rhs` the components are run by a call:: scheduler->transform(state); The call to `ComponentScheduler::create` treats the "components" option as a comma-separated list of names. The order of the components is the order that they are run in. For each name in the list, the scheduler looks up the options under the section of that name. .. code-block:: ini [hermes] components = component1, component2 [component1] # options to control component1 [component2] # options to control component2 This would create two `Component` objects, of type `component1` and `component2`. Each time `Hermes::rhs` is run, the `transform` functions of `component1` and then `component2` will be called, followed by their `finally` functions. It is often useful to group components together, for example to define the governing equations for different species. A `type` setting in the option section overrides the name of the section, and can be another list of components .. code-block:: ini [hermes] components = group1, component3 [group1] type = component1, component2 # options to control component1 and component2 [component3] # options to control component3 This will create three components, which will be run in the order `component1`, `component2`, `component3`: First all the components in `group1`, and then `component3`. .. doxygenclass:: ComponentScheduler :members: .. _sec-permissions: Permissions ~~~~~~~~~~~~~~ The ``Permissions`` class can be used to store information about which variables within an ``Options`` object are allowed to be accessed and for what purpose. This is used to control the variables used by a ``Component``. There is a hierarchy of four types of increasing permission. These are expressed using the `PermissionTypes` `enum `__: #. **ReadIfSet:** Only allowed to read variable if it is already set. #. **Read:** Can read the contents of the variable. Assumes it has already been set. #. **Write:** Can write variable. Makes no assumption about whether it has already been written or will be written again in future. #. **Final:** This will be the last component to write to the variable. Only one component may have ``Final`` permission for a given variable. The order these per permissions are listed in is significant: each higher permission implies a component also has all lower permissions. E.g., write permission implies read permission as well. Declaring Permissions for Particular Variables `````````````````````````````````````````````` Permission information for a variable is stored in a `Permissions::VarRights` object. The overwhelming majority of the permissions you would want to create can be constructed using one of the provided convenience-functions. For example:: Permissions::VarRights read_e_pressure = readOnly("species:e:pressure"); Permissions::VarRights write_d_density = readWrite("species:d:density"); Permissions::VarRights read_e_velocity_in_interior_if_set = readIfSet("species:e:velocity", Regions::Interior); Permissions can be set to apply only to a particular region of the domain (e.g., the boundary or the interior) using a `Regions` enum (see `Specifying a Region`_). Creating Permissions Objects ```````````````````````````` Permission data like that created in the previous example can be used to construct a ``Permissions`` object. These objects describe the permissions for multiple variables.:: Permissions p({readOnly("time"), readOnly("species:e:pressure"), readWrite("species:e:momentum", Regions::Interior)}); A permission applied to a section of an ``Options`` object will apply to all variables contained within that section, unless a more specific permission is also set. Therefore, if we have a state with variables ``species:e:pressure``, ``species:e:density``, ``species:e:velocity``, and ``species:e:momentum``, then the following are equivalent:: Permissions p({readOnly("species:e"), readWrite("species:e:momentum")}); Permissions p({readOnly("species:e:pressure"), readOnly("species:e:density"), readOnly("species:e:velocity"), readWrite("species:e:momentum")}); Specifying a Region ``````````````````` The `PermissionTypes` are applied to particular regions of the domain. This allows, e.g., for there to be read permissions for the interior of the domain but write permissions for the boundaries. Regions are expressed using the `Permissions::Regions` enum, which functions as a `bitset `__. You can combine regions using bitwise logical operators. .. doxygengroup:: RegionsGroup :members: Permission Substitution ``````````````````````` Variable names can include labels, marked in curly-braces, that will later be substituted (using `Permissions::substitute` and `Component::substitutePermissions`). Substitutions are necessary because, when declaring permissions for a `Component`, you may need to express that it can access some variable for all species (or all ions, all neutrals, etc.), but you won't yet know the names of all the species. For example, if you need to read the density of all species and write the collision frequency of all ions then you would write:: Permissions p({readOnly("species:{all_spcies}:density"), readWrite("species:{ions}:collision_frequency"}); If there are species e, d, d+, h, and h+ then the above will be equivalent to:: Permissions p({readOnly("species:e:density"), readOnly("species:d:density"), readOnly("species:d+:density"), readOnly("species:h:density"), readOnly("species:h+:density"), readWrite("species:d+:collision_frequency"), readWrite("species:h+:collision_frequency")}); These substitutions will be performed in `Component::declareAllSpecies`. See the documentation for that method for a full list of the substitutions which it can perform. It can also be useful to define your own substitutions, to save repetitive declarations. For example, you could declare read permissions for electron density, pressure, temperature, velocity, and momentum as follows:: Permissions p({readOnly("species:e:{inputs}"); p.substitute({"density", "pressure", "temperature", "velocity", "momentum"}); This is equivalent to having written:: Permissions p({readOnly("species:e:density"), readOnly("species:e:pressure")}, readOnly("species:e:temperature")}, readOnly("species:e:velocity")}, readOnly("species:e:momentum")}); Permission Factory Functions ```````````````````````````` .. doxygengroup:: PermissionFactories :members: Permissions Class ````````````````` .. doxygenclass:: Permissions :members: Further Implementation Details `````````````````````````````` The above information should be sufficient for users that are developing or modifying components. The following explains in more detail how permission data is stored and should be read by anyone looking to modify the `Permissions` or `GuardedOptions` classes. Permission information for a variable gets stored in `Permissions::AccessRights` objects, which are arrays of `Regions`. Each element of the array corresponds to information about a permission level: ``{read_if_set, read, write, final}``. To access the element for a desired permission level, you can index the array with the corresponding member of the `PermissionTypes` enum:: Permissions::AccessRights rights; Regions read_regions = rights[PermissionTypes::Read]; Regions write_regions = rights[PermissionTypes::Write]; The contents of each element of an `Permissions::AccessRights` array is the set of regions for which the permissions apply. For example:: Permissions::AccessRights read_boundaries_if_set = {Regions::Boundaries, Regions::Nowhere, Regions::Nowhere, Regions::Nowhere}; Permissions::AccessRights read_interior_write_boundaries = {Regions::Nowhere, Regions::Interior, Regions::Boundaries, Regions::Nowhere}; Permissions::AccessRights final_write_all_regions = {Regions::Nowhere, Regions::Nowhere, Regions::Nowhere, Regions::All}; The `Permissions::VarRights` struct is used to pair a variable name with a `Permissions::AccessRights` array containing the permission information for that variable. GuardedOptions ~~~~~~~~~~~~~~ ``GuardedOptions`` objects combine a `Permissions` object and an `Options` object. They can be indexed just like normal ``Options`` objects but will return another ``GuardedOptions``, wrapping the result. In order to read or write the contents of a ``GuardedOptions`` object you must use the ``get()`` or ``getWritable()`` methods, respectively. These will return the underlying (const) ``Options`` object, if you have the necessary permissions to access it. Otherwise, they will raise an exception. If ``CHECKLEVEL`` is 1 or above, then the ``GuardedOptions`` will track which variables have actually been accessed. Lists of unread/unwritten variables can be returned with the ``unreadItems()`` and ``unwrittenItems()`` methods. If ``CHECKLEVEL`` is zero then calling these methods will raise an exception. .. doxygenclass:: GuardedOptions :members: .. note:: When indexing a ``GuardedOptions`` object, it will create a new ``GuardedOptions`` on-demand. This is unlike with a normal ``Options`` object which returns a reference to a preexisting child ``Options`` object. You generally should not store ``GuardedOptions`` by reference. You may be able to pass them by reference, but this requires you to think carefully about whether the argument is going to be an r-value or an l-value. .. _sec-reactions-dev: Reactions ~~~~~~~~~ The code for handling reactions was significantly refactored after version 1.4.0. Currently, reactions involving Hydrogen isotopes and Helium use the new framework, whilst reactions involving impurity isotopes use their own set of classes. The intent is that all reaction classes in the code will eventually be migrated to the new framework. Executive summary ````````````````` The hierarchy of classes in the framework is: - `Reaction` - `CXReaction` - `IznRecReaction` - `IznReaction` - `RecReaction` Most reaction source term contributions are handled in the the :code:`transform_impl` method of the base :code:`Reaction` class. The source of the reaction rate data can be configured via the input file; if no reaction data type and data ID is set, a sensible default will be set on a per-reaction basis (e.g. Amjuel, H.2_3.1.8 for Hydrogen isotope charge exchange, :code:`h + d+ -> d + h+`). Subclasses of :code:`Reaction` set up appropriate diagnostics and :ref:`Permissions` for particular reactions in their constructors and/or implement a :code:`transform_additional` method to include source term contributions that aren't captured by :code:`Reaction::transform_impl`. The following subsections describe the general approach used to compute reaction sources and provide guidance to developers on adding new subclasses to handle particular reactions. The stoichiometry matrix ```````````````````````` A core part of the approach is to convert reaction strings into a stoichiometry matrix. stoichiometry matrices specify the population changes associated with a *mechanism* or system of reactions, with each column corresponding to a reaction and each row to a species. For a configuration .. _code-reactions_config_eg: .. code-block:: ini [reactions] type = ( h + e -> h+ + 2e, # R1 h+ + e -> h, # R2 he + e -> he+ + 2e, # R3 he+ + e -> he, # R4 ) the stoichiometry matrix would be: .. _fig-stoichiometry_matrix_eg: .. figure:: figs/stoichiometry_matrix_eg.* :alt: Example of a stoichiometry matrix for the reactions configured in :numref:code-reactions_config_eg In the `Reaction` base class, the constructor instantiates a `ReactionParser`, which assembles one column of this matrix from the reaction string. .. note:: The reaction string is extracted from the comma-separated list of strings provided as the 'type' of the reactions component in the input file. All reaction classes (including those that do not currently inherit from `Reaction`) inherit from `ReactionBase`, which ensures that subclasses are paired up with reaction strings correctly. It is anticipated that changes to the input file format and/or the migration of all reaction classes to the new framework will make :code:`ReactionBase` unnecessary in the future. .. tip:: In addition to computing population changes, `ReactionParser` can be used to identify different subsets of species via the `species_filter` enum. For instance, to obtain the names of all 'heavy' (non-electron) reactants: .. code-block:: cpp auto heavy_reactants = this->parser.get_species(species_filter::reactants, species_filter::heavy); Source term equations ````````````````````` Source terms related to the population changes are calculated from the stoichiometric coefficients in `Reaction::transform_impl`. The reaction rate is calculated as .. math:: K &= M \langle\sigma_v\rangle \\ &= \prod_{r} n_r \langle\sigma_v\rangle where :math:`\langle\sigma_v\rangle` is the reaction cross section and :math:`M` is called the *mass action factor*. The density source for a species :math:`s` follows directly from the rate: .. _fig-reactions_dens_src_eqn: .. math:: S = C_s K where :math:`C_s` is the population change for species :math:`s` (from the stoichiometry matrix) and :math:`K` is the reaction rate calculated above. The momentum source (:math:`F`) and energy source (:math:`E`) associated with a species population change depend on whether the reaction in question is a net consumer (:math:`C_s < 0`) or a net producer (:math:`C_s > 0`) of the species: .. math:: F = \begin{cases} C_s K (NV)_s & \textrm{for}\ C_s < 0 \\ f_{NV,s} K \sum_{r'} -C_{r'}(NV)_{r'} & \textrm{for}\ C_s > 0 \\ 0 & \textrm{otherwise} \\ \end{cases} .. math:: E = \begin{cases} C_s K \mathcal{E}_s & \textrm{for}\ C_s < 0 \\ f_{\mathcal{E},s} K \sum_{r'} -C_{r'}\mathcal{E}_{r'} & \textrm{for}\ C_s > 0 \\ 0 & \textrm{otherwise} \\ \end{cases} If the species is consumed, it has negative sources proportional to its current momentum (:math:`{NV}_s`) and energy (:math:`\mathcal{E}_s`). If the species is produced, it has positive sources, which are proportional to the total momentum and energy of all consumed reactants (:math:`r'`; implying :math:`C_r < 0`). :math:`f_{NV,s}` and :math:`f_{\mathcal{E},s}` are weights / splitting factors used to distribute momentum and energy between products. The default is to weight by mass for momentum, and by population change (number) for energy, that is .. math:: \begin{aligned} f_{NV,s} &= \frac{\Delta M_s}{\Delta M} \\ &= \frac{C_s m_s}{\sum_{p'}C_{p'} m_{p'}} \end{aligned} .. math:: \begin{aligned} f_{\mathcal{E},s} &= \frac{\Delta N_s}{\Delta N} \\ &= \frac{C_s}{\sum_{p'}C_{p'}} \end{aligned} where :math:`m` are species masses, :math:`C` are population changes and subscripts :math:`s` and :math:`p'` refer to the target species and to a *produced* species (:math:`C > 0`) respectively. :math:`\Delta M` and :math:`\Delta N` are the change in mass and particle number associated with one instance of the reaction. These default factors can be overridden by calling `set_energy_channel_weight` and `set_momentum_channel_weight` in the constructor of a reaction subclass. For example, for charge exchange, momentum and energy need to be transferred 1. from a reactant that donates an electron ('r1') to the resulting product ('p1') and 2. from a reactant that accepts an electron ('r2') to the resulting product ('p2'). The following code sets channels to enforce those requirements: .. code-block:: cpp // Energy exchange between r1 and p1 set_energy_channel_weight(this->r1, this->p1, 1.0); set_energy_channel_weight(this->r1, this->p2, 0.0); // Energy exchange between r2 and p2 set_energy_channel_weight(this->r2, this->p1, 0.0); set_energy_channel_weight(this->r2, this->p2, 1.0); // Momentum exchange between r1 and p1 set_momentum_channel_weight(this->r1, this->p1, 1.0); set_momentum_channel_weight(this->r1, this->p2, 0.0); // Momentum exchange between r2 and p2 // (Can be turned off for neutral-ion CX using config option) if (this->has_neutral_reactant && this->no_neutral_cx_mom_gain) { set_momentum_channel_weight(this->r2, this->p2, 0.0); } else { set_momentum_channel_weight(this->r2, this->p2, 1.0); } set_momentum_channel_weight(this->r2, this->p1, 0.0); .. note:: So-called `participation factors` are included in all source term calculations in the code. For now, these are all set to 1, reproducing the equations above exactly, but in future they could be made configurable to allow users to experiment with tuning the contribution of particular species to reaction source terms. Source term calculation ``````````````````````` Density, momentum and energy source fields are constructed in `Reaction::transform_impl` using the `RateHelper` class to calculate reaction rates and collision frequencies. On instantiation, :code:`RateHelper` assembles maps containing (pointers to) all reactant densities and any other fields required as inputs. `RateHelper::calc_rates` is "run-time templated", allowing it to accept rate functions with 1D and 2D parameterisations. It also includes an optional flag which causes reaction input parameters to be averaged in the parallel direction. Averaging is turned on by default, but can be disabled in a subclass by setting :code:`do_parallel_averaging = false`. Having calculated the reaction rate, :code:`Reaction::transform_impl` constructs sources by looping over the species population changes provided by :code:`ReactionParser`. .. note:: The loop to compute momentum and energy sources includes a subtle workaround for reactions like symmetric charge exchange. Rather than using the regular stoichiometric coefficients, which are all zero in such cases, the loop uses the result of `ReactionParser::get_mom_energy_pop_changes`. For symmetric reactions, this returns a :code:`std::multimap` in which each species appears twice; once as a reactant, with a negative population change and once as a product, with a positive population change. This effectively turns reactions like :math:`\textrm{a} + \textrm{b} \rightarrow \textrm{b} + \textrm{a}` into :math:`\textrm{a} + \textrm{b} \rightarrow \textrm{b}' + \textrm{a}'`, allowing the momentum and energy transfer to be computed as usual, before mapping sources back to the correct species. Adding new reaction subclasses `````````````````````````````` New reaction subclasses should inherit directly or indirectly from `Reaction`. They are registered in the code with the reaction string as an argument, e.g. .. code-block:: cpp struct MyReaction : public Reaction { ... } RegisterComponent register_myreaction("a + b -> c + d"); Any sources associated with the reaction that are not captured by `transform_impl` should be added by overriding `transform_additional` and using `update_source` to set source terms and associated diagnostics. Adding new reaction diagnostics ``````````````````````````````` New diagnostics are configured by making calls to `add_diagnostic` in the subclass constructor; e.g. for a density source associated with the neutral reactant in charge exchange: .. code-block:: cpp :caption: Example of adding a reaction diagnostic. :code:`r1` is the name of a reactant that charge-exchanges with reactant :code:`r2` creating product :code:`p1`. add_diagnostic( this->r1, fmt::format("S{:s}{:s}_cx", this->r1, this->r2), fmt::format("Particle transfer to {:s} from {:s} due to CX with {:s}", this->r1, this->p1, this->r2), ReactionDiagnosticType::density_src, standard_name, identity, "particle transfer"); An optional function pointer can be used to set a transformation that will be applied when the diagnostic is updated. The value passed here (also the default) is :code:`identity`, meaning that the diagnostic and associated source term will be updated in the same way. The fourth argument is an :code:`enum` value that determines, among other things, which metadata will be set when the diagnostic is written out. Each call adds an instance of `ReactionDiagnostic` to a `diagnostics map`, indexed by species name and a `ReactionDiagnosticType` enum. The map entries are updated via `update_source()`, which wraps the `add`, `subtract`, and `set` operations, applying them to both a source term and to an associated diagnostic at the same time. For example, this code .. code-block:: cpp update_source>(state, "h+", ReactionDiagnosticType::momentum_src, momentum_source); would add the :code:`momentum_source` field to :code:`state["species"]["h+"]["momentum_source"]`, then apply the diagnostic transform (which defaults to an identity function) to :code:`momentum_source` and add the result to any diagnostics of type :code:`momentum_src` associated with species :code:`h+`. .. note:: `add`, `subtract`, and `set` can still be used directly to update the state, as usual, if there is no associated diagnostic to update. Finally, diagnostics are written out in `outputVars()`, which simply iterates over the `diagnostics` map, copying fields into the output state. Appropriate metadata is set for each diagnostic depending on the associated :code:`ReactionDiagnosticType`. Adding new reaction data ```````````````````````` If the data being added is of a new type (i.e. one not handled by an existing subclass of `ReactionData`), developers will need to: 1. Create a new class that inherits from `ReactionData` (or from `ReactionDataWithCoeffs` if the rate calculation computes a fit value from coefficients). 2. Implement all pure virtual functions, including those that calculate rates with different parameterisations. 3. Add code to the constructor of the new subclass to read data into the `coeffs` data member, for file-based data sources. 4. Add a new value to the `ReactionDataType` enum. 5. Register the data class in a new header file. .. code-block:: cpp RegisterReactionData register_myclass("mylabel"); where "mylabel" must match the value of the `ReactionDataType` enum added in step 4, converted to lowercase. 6. Add an include statement for your class's header in :file:`hermes-3.hxx`. For both new and existing data types, developers need to 1. Add any associated json data files to the :file:`./json_database` directory. 2. Register a `Reaction` subclass for the reaction string. 3. Add any default data IDs to the `generate_default_data_ids_map` function in :file:`reaction_settings.hxx`. Step 3. is optional, but without it, users will always need to specify the data type and ID/label for the new reaction in their input file. Example: To add He+ ionisation: :code:`he+ + e -> he+2 + 2e` one could 1. Create :file:`./json_database/AMJUEL_H.4_2.2C.json` containing the relevant coefficients, matching the format of existing Amjuel json files. 2. Register the reaction with the standard ionisation reaction class: :code:`RegisterComponent register_izn_he2p("he+ + e -> he+2 + 2e");` 3. Set the new data to be the default source for He+ ionisation in :file:`include/reaction_settings.hxx` .. code-block:: cpp static inline void generate_default_data_ids_map() { ... add_default_id(ReactionDataTypes::Amjuel, "he+ + e -> he+2 + 2e", ReactionCoeffTypes::sigma_v, "H.4_2.2C"); ... } .. note:: For ionisation/recombination one would also need to add data to compute the losses associated with the ionisation potential energy cost and photon emission. At time of writing, appropriate Amjuel data is not available for this reaction, which is why He+ ionisation hasn't already been added. .. _sec-tests: Tests -------------- The specification of the Toro tests used here is taken from `Walker (2012) `__, originally from Toro's book `Riemann Solvers and Numerical Methods for Fluid Dynamics `__. 1D-recycling-dthe ~~~~~~~~~~~~~~~~~ This is a comprehensive 1D test featuring three species (deuterium tritium and helium) as well as all of the parallel closure terms apart from electron viscosity. It includes ionisation and recombination reactions for all species as well as charge exchange for D-D, T-T, D-T and T-D species pairs. The test checks the values of charge exchange channels in the final domain cell against a reference. The test file can be used to generate the test data if `gen_data` is set to `True` in the beginning of the script. 2D-production ~~~~~~~~~~~~~~ This is a test representing the most common production fidelity simulations at the time of implementation. It features a double null 2D tokamak from an upcoming publication with deuterium only, a neutral puff, H-mode profiles of anomalous diffusion coefficients, a neutral pump, fast/thermal recycling, fast/thermal reflection as well as decay length boundary conditions on the SOL and PFR. Viscosity, impurities, currents and drifts are neglected. The test downloads the restart files and the grid from Zenodo and runs a few very short timesteps. It compares electron pressure and ion momentum at the four targets, outer midplane, SOL boundary and core/pfr boundaries. If the test fails, it will print the degree of error at each boundary. The test file can be used to generate the test data if `gen_data` is set to `True` in the beginning of the script. 2D-recycling ~~~~~~~~~~~~~~ This test uses the same restart file as 2D-production, but tests target recycling by reproducing the relevant recycling.cxx section in Python. It compares the source of recycled neutral density and energy as well as the pumped neutral density and energy sources between the calculation and the reference simulation. In the future, this test will be extended to radial recycling. This test is a reproduction of C++ code, and therefore there is no golden answer reference to update. This test plots results and can be used to help with developing the recycling component. There is a `plot` flag near the beginning of the file. 1D fluid (MMS) ~~~~~~~~~~~~~~ ``tests/integrated/1D-fluid`` This convergence test using the Method of Manufactured Solutions (MMS) solves fluid equations in the pressure form: .. math:: \begin{aligned} \frac{\partial n}{\partial t} &= -\nabla\cdot\left(n\mathbf{b}v_{||}\right) \\ \frac{\partial p}{\partial t} &= -\nabla\cdot\left(p\mathbf{b}v_{||}\right) - \frac{2}{3}p\nabla\cdot\left(\mathbf{b}v_{||}\right) \\ \frac{\partial}{\partial t}\left(mnv_{||}\right) &= -\nabla\cdot\left(nv_{||}\mathbf{b}v_{||}\right) - \partial_{||}p \end{aligned} .. figure:: figs/fluid_norm.png :name: fluid_norm :alt: :width: 60% Sod shock ~~~~~~~~~~~~~~ ``tests/integrated/sod-shock`` and ``tests/integrated/sod-shock-energy`` Euler equations in 1D. Starting from a state with a jump at the middle of the domain. Left state density, velocity and pressure are :math:`\left(\rho_L, u_L, p_L\right) = \left(1.0, 0, 1.0\right)` Right state :math:`\left(\rho_R, u_R, p_R\right) = \left(0.125, 0, 0.1\right)`. The result is shown in figure below at time :math:`t = 0.2` for different resolutions in a domain of length 1. The solid black line is the analytic solution. .. figure:: figs/sod_shock.png :name: sod_shock :alt: :width: 60% When evolving pressure the position of the shock front lags the analytic solution, with the pressure behind the front slightly too high. This is a known consequence of solving the Euler equations in non-conservative form. If instead we evolve energy (internal + kinetic) then the result is much closer to the analytic solution. .. figure:: figs/sod_shock_energy.png :name: sod_shock_energy :alt: :width: 60% Toro test 1 ~~~~~~~~~~~~~~ ``tests/integrated/toro-1`` Toro's test problem #1, from `Riemann Solvers and Numerical Methods for Fluid Dynamics `__ is a variation of Sod's shock tube problem. The left state is moving into the right, increasing the speed of the resulting shock. Left state :math:`\left(\rho_L, u_L, p_L\right) = \left(1.0, 0.75, 1.0\right)` Right state :math:`\left(\rho_R, u_R, p_R\right) = \left(0.125, 0, 0.1\right)`. The size of the domain is 5, and the reference result is given at time :math:`t = 0.8`. Toro test 2 ~~~~~~~~~~~~~~ ``tests/integrated/toro-2`` and ``tests/integrated/toro-2-energy`` Toro's test problem #2 tests robustness to diverging flows and near-zero densities. The initial state has constant density and temperature, but a jump in velocity. Left state :math:`\left(\rho_L, u_L, p_L\right) = \left(1.0, -2.0, 0.4\right)` Right state :math:`\left(\rho_R, u_R, p_R\right) = \left(1.0, 2.0, 0.4\right)`. The result in a domain of length 5 at time :math:`t=0.6` is shown below. .. figure:: figs/toro-2.png :name: toro-2 :alt: :width: 60% Toro test 3 ~~~~~~~~~~~~~~ ``tests/integrated/toro-3`` and ``tests/integrated/toro-3-energy`` Toro's test problem #3 contains a strong shock close to a contact discontinuity. Left initial state :math:`\left(\rho_L, u_L, p_L\right) = \left(1.0, 0, 1000.0\right)` Right state :math:`\left(\rho_R, u_R, p_R\right) = \left(1.0, 0, 0.01\right)`. Time :math:`t = 0.04`. When evolving pressure, the simulation is robust but the density peak does not converge to the analytic solution (solid black line): .. figure:: figs/toro-3.png :name: toro-3 :alt: :width: 60% However by evolving energy the result converges towards the analytic solution: .. figure:: figs/toro-3-energy.png :name: toro-3-energy :alt: :width: 60% Toro test 4 ~~~~~~~~~~~~~~ ``tests/integrated/toro-4`` and ``tests/integrated/toro-4-energy`` Toro's test problem #4 produces two right-going shocks with a contact between them. Left state :math:`\left(\rho_L, u_L, p_L\right) = \left(5.99924, 19.5975, 460.894\right)` Right state :math:`\left(\rho_R, u_R, p_R\right) = \left(5.99242, -6.19633, 46.0950\right)`. Result at time :math:`t = 0.15`. Toro test 5 ~~~~~~~~~~~~~~ ``tests/integrated/toro-5`` and ``tests/integrated/toro-5-energy`` The initial conditions for Toro's test problem #5 are the same as test #3, but the whole system is moving to the left at a uniform speed. The velocity is chosen so that the contact discontinuity remains almost stationary at the initial jump location. Left state :math:`\left(\rho_L, u_L, p_L\right) = \left(1, -19.59745, 1000.0\right)` Right state :math:`\left(\rho_R, u_R, p_R\right) = \left(1, -19.59745, 0.01\right)`. Result at time :math:`t = 0.03`. Drift wave ~~~~~~~~~~~~~~ ``tests/integrated/drift-wave`` This calculates the growth rate and frequency of a resistive drift wave with finite electron mass. The equations solved are: .. math:: \begin{aligned} \frac{\partial n_i}{\partial t} =& -\nabla\cdot\left(n_i\mathbf{v}_{E\times B}\right) \\ n_e =& n_i \\ \frac{\partial}{\partial t}\nabla\cdot\left(\frac{n_0 m_i}{B^2}\nabla_\perp\phi\right) =& \nabla_{||}J_{||} = -\nabla_{||}\left(en_ev_{||e}\right) \\ \frac{\partial}{\partial t}\left(m_en_ev_{||e}\right) =& -\nabla\cdot\left(m_en_ev_{||e} \mathbf{b}v_{||e}\right) + en_e\partial_{||}\phi - \partial_{||}p_e - 0.51\nu_{ei}n_im_ev_{||e} \end{aligned} Linearising around a stationary background with constant density :math:`n_0` and temperature :math:`T_0`, using :math:`\frac{\partial}{\partial t}\rightarrow -i\omega` gives: .. math:: \begin{aligned} \tilde{n} =& \frac{k_\perp}{\omega}\frac{n_0}{BL_n}\tilde{\phi} \\ \tilde{\phi} =& -\frac{k_{||}}{\omega k_\perp^2}\frac{eB^2}{m_i}\tilde{v_{||e}} \\ \omega m_e \tilde{v_{||e}} =& -ek_{||}\tilde{\phi} + ek_{||}\frac{T_o}{n_0}\tilde{n} - i0.51\nu_{ei}m_e\tilde{v_{||e}} \end{aligned} where the radial density length scale coming from the radial :math:`E\times B` advection of density is defined as .. math:: \frac{1}{L_n} \equiv \frac{1}{n_0}\frac{\partial n_0}{\partial r} Substituting and rearranging gives: .. math:: i\left(\frac{\omega}{\omega*}\right)^3 \frac{\omega_*}{0.51\nu_{ei}} = \left(\frac{\omega}{\omega_*} - 1\right)\frac{i\sigma_{||}}{\omega_*} + \left(\frac{\omega}{\omega*}\right)^2 or .. math:: \frac{\omega_*}{0.51\nu_{ei}}\left(\frac{\omega}{\omega_*}\right)^3 + i\left(\frac{\omega}{\omega_*}\right)^2 - \frac{\sigma_{||}}{\omega_*}\left(\frac{\omega}{\omega_*}\right) + \frac{\sigma_{||}}{\omega_*} = 0 where .. math:: \begin{aligned} \omega_* =& \frac{k_\perp T_0}{BL_n} \\ \sigma_{||} =& \frac{k_{||}^2}{k_\perp^2}\frac{\Omega_i\Omega_e}{0.51\nu_{ei}} \\ \Omega_s =& eB / m_s \end{aligned} This is a cubic dispersion relation, so we find the three roots (using NumPy), and choose the root with the most positive growth rate (imaginary component of :math:`\omega`). .. figure:: figs/drift-wave.png :name: drift-wave :alt: Comparison of drift-wave growth rate (top) and frequency (bottom) :width: 60% Alfven wave ~~~~~~~~~~~~~~ The equations solved are .. math:: \begin{aligned} \frac{\partial}{\partial t}\nabla\cdot\left(\frac{n_0 m_i}{B^2}\nabla_\perp\phi\right) =& \nabla_{||}J_{||} = -\nabla_{||}\left(en_ev_{||e}\right) \\ \frac{\partial}{\partial t}\left(m_en_ev_{||e} - en_eA_{||}\right) =& -\nabla\cdot\left(m_en_ev_{||e} \mathbf{b}v_{||e}\right) + en_e\partial_{||}\phi - 0.51\nu_{ei}n_im_ev_{||e} \\ J_{||} =& \frac{1}{\mu_0}\nabla_\perp^2 A_{||} \end{aligned} Linearising around a stationary background with constant density :math:`n_0` and temperature :math:`T_0`, using :math:`\frac{\partial}{\partial t}\rightarrow -i\omega` gives: .. math:: \begin{aligned} \tilde{\phi} =& -\frac{k_{||}}{\omega k_\perp^2}\frac{eB^2}{m_i}\tilde{v_{||e}} \\ \omega \left( m_e \tilde{v_{||e}} - e\tilde{A}_{||}\right) =& -ek_{||}\tilde{\phi} - i0.51\nu_{ei}m_e\tilde{v_{||e}} \\ en_0\tilde{v_{||e}} =& -\frac{k_\perp^2}{\mu_0}\tilde{A}_{||} \end{aligned} Rearranging results in a quadratic dispersion relation: .. math:: \omega^2\left(1 + \frac{k_\perp^2 c^2}{\omega_{pe}^2}\right) + i 0.51\nu_{ei}\frac{k_\perp^2 c^2}{\omega_{pe}^2}\omega - k_{||}^2V_A^2 = 0 where :math:`V_A = B / \sqrt{\mu_0 n_0 m_i}` is the Alfven speed, and :math:`c / \omega_{pe} = \sqrt{m_e / \left(\mu_0 n_0 e^2\right)}` is the electron skin depth. When collisions are neglected, we obtain the result .. math:: \omega^2 = V_A^2\frac{k_{||}^2}{1 + k_\perp^2 c^2 / \omega_{pe}^2} .. figure:: figs/alfven-wave.png :name: alfven-wave :alt: Alfven wave speed, as function of parallel and perpendicular wavenumbers :width: 60% Collision frequency selection ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ There are two simple integrated tests to make sure that the collision frequency selection is correct across `neutral_mixed`, `evolve_pressure`, `ion_viscosity` and `neutral_parallel_diffusion`. A minimal 3D geometry is run for one RHS evaluation, and the test checks the log file to make sure the correct collisionalities were selected. One of the tests is for the `multispecies` mode across all components, while the other is for `braginskii` for plasma and `afn` for neutrals.