# Chemical Reactions

MFIX-Exa categorizes chemical reactions into two classes:

- Eulerian reactions are homogeneous fluid reactions, reactions that only reference fluid species. These reactions are calculated by looping over grid cells, and are defined in units of moles per second per cubic meter \([mol \cdot s^{-1} \cdot m^{-3}]\).

- Lagrangian reactions include homogeneous particle reactions and heterogeneous fluid-particle reactions. These reactions are calculated by looping over each particle, and are defined in units of moles per second \([mol \cdot s^{-1}]\).

Users implement reaction rates by modifying the header file, `src/usr/mfix_usr_reactions_rates_K.H`.

The `operator()` method for the `EulerianReactionRates` class gets two arguments, `reactions` and `fluid`, which are two classes from which the user can extract the data needed to implement the corresponding rates.

    class EulerianReactionRates
    {
      public:
        /**
         *  Reactions class gives access to the following reactions quantities:
         *   - reaction_idx("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
         *   - rate(q), with q an integer from 0 to number of Eulerian reactions - 1
         *   - rate("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
         *   - time(), the current simulation time
         *   - dt(), the current simulation time step
         *   - cell_volume(), this cell geometric volume
         *   - vfrac(), this cell geometric volume fraction
         *
         *  Fluid class gives access to the following fluid quantities:
         *   - shear_viscosity()
         *   - volume_fraction()
         *   - density()
         *   - velocity() and velocity(dir), with dir an integer from 0 to 2
         *   - thermodynamic_pressure()
         *   - temperature()
         *   - specific_enthalpy()
         *   - species_idx("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - mass_fraction(n), with n an integer from 0 to number of fluid species - 1
         *   - mass_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - molar_concentration(n), with n an integer from 0 to number of fluid species - 1
         *   - molar_concentration("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - molar_mass(n), with n an integer from 0 to number of fluid species - 1
         *   - molar_mass("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - average_molar_mass()
         *   - mole_fraction(n), with n an integer from 0 to number of fluid species - 1
         *   - mole_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
         */
        template <class Reactions, class Fluid>
        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        void operator() (Reactions& reactions,
                         const Fluid& fluid) const
        {
          // User defined Eulerian reactions rates go here
        }
    };

The `operator()` of the `LagrangianReactionRates` class contains a third argument, `solids`, which is a class containing data for the particle involved in the reaction.

    class LagrangianReactionRates
    {
      public:
        /**
         *  Reactions class gives access to the following reactions quantities:
         *   - reaction_idx("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
         *   - rate(q), with q an integer from 0 to max number of Lagrangian reactions - 1
         *   - rate("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
         *   - time(), the current simulation time
         *   - dt(), the current simulation time step
         *   - cell_volume(), cell geometric volume of cell containing this particle
         *   - vfrac(), geometric volume fraction of cell containing this particle
         *
         *  Fluid class gives access to the following fluid quantities:
         *   - shear_viscosity()
         *   - volume_fraction()
         *   - density()
         *   - velocity() and velocity(dir), with dir an integer from 0 to 2
         *   - thermodynamic_pressure()
         *   - temperature()
         *   - specific_enthalpy()
         *   - species_idx("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - mass_fraction(n), with n an integer from 0 to max number of fluid species
         *   - mass_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - molar_concentration(n), with n an integer from 0 to max number of fluid species - 1
         *   - molar_concentration("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - molar_mass(n), with n an integer from 0 to max number of fluid species - 1
         *   - molar_mass("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *   - average_molar_mass()
         *   - mole_fraction(n), with n an integer from 0 to max number of fluid species - 1
         *   - mole_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
         *
         *  Solids class gives access to the following solids quantities:
         *   - id()
         *   - cpu()
         *   - position() and position(dir), with dir an integer from 0 to 2
         *   - radius()
         *   - volume()
         *   - density()
         *   - mass()
         *   - oneOverI(), inverse of the particle momentum of inertia
         *   - statistical_weight()
         *   - velocity() and velocity(dir), with dir an integer from 0 to 2
         *   - temperature()
         *   - specific_heat_capacity()
         *   - species_idx("species_name"), with species_name a string matching one of the solids species defined in inputs
         *   - mass_fraction(n), with n an integer from 0 to max number of solids species - 1
         *   - mass_fraction("species_name"), with species_name a string matching one of the solids species defined in inputs
         *   - molar_concentration(n), with n an integer from 0 to max number of solids species - 1
         *   - molar_concentration("species_name"), with species_name a string matching one of the solids species defined in inputs
         *   - molar_mass(n), with n an integer from 0 to max number of solids species - 1
         *   - molar_mass("species_name"), with species_name a string matching one of the solids species defined in inputs
         *   - average_molar_mass()
         *   - mole_fraction(n), with n an integer from 0 to max number of solids species - 1
         *   - mole_fraction("species_name"), with species_name a string matching one of the solids species defined in inputs
         */
        template <class Reactions, class Fluid, class Solids>
        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        void operator() (Reactions& reactions,
                         const Fluid& fluid,
                         const Solids& solids) const
        {
          // User defined Lagrangian reactions rates go here
        }
    };

To simply defining chemical reactions, fluid and solids species-related quantities can be accessed by calling the variable function with the species name specified in the inputs file. For example, if the fluid contains a species called `CO`, then the species mass fraction for `CO` can be accessed as follows:

    // Mass fraction of CO species
    amrex::Real const X_CO = fluid.mass_fraction("CO");

This feature allows users to write readable reaction rates, especially when the number of chemical species is large. However, at runtime this approach requires the solver to search for the species index associated with the name. The computational cost of the search can be non-trivial because it is performed for each iteration of the chemistry solver every time step, for all grid cells and all particles.

A python-script is provided that substitutes the string references with their corresponding indices to remove the runtime search. The tool is run by passing the `inputs` file and the reaction header, `mfix_usr_reactions_rates_K.H`. The tool creates a backup copy of the original reaction header, appending a time stamp for when the tool was run, before applying the optimizations.

Example usage:

    python tools/Chemistry/optimize_usr_rates.py inputs mfix/src/usr/mfix_usr_reactions_rates_K.H

# Example rates UDFs

The following code shows a simple reaction rate for methane combustion. The molar concentrations for oxygen (`O2`) and methane (`CH4`) are accessed using the names provided in the inputs file. Similarly, the reaction rate is stored using the name given in the inputs file.

Listing 19 Snippet of inputs defining fluid species and adiabatic flame reaction

    # Fluid model settings
    # -----------------------------------------------------------------------
    fluid.solve = fluid

    fluid.species = CH4 O2 CO2 H2O N2

    fluid.thermodynamic_pressure = 101325.
    fluid.reference_temperature = 298.15

    fluid.species = Air Vapor

    fluid.viscosity.molecular = constant
    fluid.viscosity.molecular.constant = 1.86e-5

    fluid.specific_heat = mixture

    fluid.thermal_conductivity          = constant
    fluid.thermal_conductivity.constant = 0.02662

    # Chemical reactions:
    # -----------------------------------------------------------------------
    chemistry.solve = AdiabaticFlame

    chemistry.AdiabaticFlame.reaction = CH4(g) + 2O2(g) --> CO2(g) + 2H2O(g)

Listing 20 Snippet of reaction rate header defining reaction rate for adiabatic flame

    class EulerianReactionRates
    {
      public:
        AMREX_GPU_HOST_DEVICE
        EulerianReactionRates () = default;

        template <class Reactions, class Fluid>
        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        void operator() (Reactions& reactions,
                         const Fluid& fluid) const
        {

          // Adiabatic flame:
          // CH4 + 2O2 --> CO2 + 2H2O  (mol/m^3.s)
          //---------------------------------------------------------------------//
          // Note: The CH4 combustion rate is artificial and used for the
          // adiabatic flame test case.

          amrex::Real const c_O2  = fluid.molar_concentration("O2");   // [mol/m^3]
          amrex::Real const c_CH4 = fluid.molar_concentration("CH4");  // [mol/m^3]

          amrex::Real k0(2.0e-1); // reaction rate coefficient

          reactions.rate("AdiabaticFlame") = k0 * c_O2 * c_CH4;

        }
    };