# Defining chemical reactions

The following inputs are defined using the prefix `chemistry`:

| Description | Type | Default
--- | --- | --- | ---
solve | Specify the name(s) of the chemical reactions or set to None to disable the chemistry solver. The name(s) assigned to the chemistry solver are used to specify the chemical reactions equations. | String |
[reaction0].reaction | Chemical formula for the given reaction. The string given as input must not contain white space and the reaction direction has to be specified as ‘-->’ or ‘<--’. Chemical species phases must be defined as ‘(g)’ for the fluid phase or ‘(s)’ for the solid phase. | String |
mass_balance_tolerance | Tolerance used to test chemical formula conserves mass. \(\left|\sum\text{products - \sum\text{reactants}\right| < \text{tolerance}\) | Real | 1e-12
solids.update_type | How to update solids quantities due to chemical reactions. Options: `ConstantVolume` - keep particle volume constant `ConstantDensity` - keep particle density constant |  | ConstantVolume
solids.mass_threshold | Sets a threshold value for the particles’ mass. When updating the solids quantities in chemistry, check whether the solids mass is below the threshold, and in that case set the particle as invalid and remove it. | Real | 0
solids.radius_threshold | Sets a threshold value for the particles’ radius. When updating the solids quantities in chemistry, check whether the solids radius is below the threshold, and in that case set the particle as invalid and remove it. | Real | 0
solids.density_threshold | Sets a threshold value for the particles’ density. When updating the solids quantities in chemistry, check whether the solids density is below the threshold, and in that case set the particle as invalid and remove it. | Real | 0

Choosing the integrator type for the chemistry ODE integration operation that is performed to determine the fluid and solids phases transfer quantities due to chemical reactions. The `StiffSolver` class of integrators is inspired by the class of integrators in AMReX-Astro/Microphysics (https://github.com/AMReX-Astro/Microphysics)

| Description | Type | Default
--- | --- | --- | ---
chemistry.integrator | Specified name of the chemistry ODE integrator type. Options (case-insensitive): `ForwardEuler` `StiffSolver::ForwardEuler` `StiffSolver::BackwardEuler` `StiffSolver::VODE` | String | ForwardEuler

One can define the environment variable `VODE_JACOBIAN_CACHING` at compile time to enable caching of the numerical approximation of the Jacobian matrix (the derivative of the ODE right-hand side), applicable when the VODE integrator is selected. The following inputs can be specified using the prefix `chemistry.integrator`:

| Description | Type | Default
--- | --- | --- | ---
burner_verbose | Enable printing of some integration statistics. | Int | 0
ode_max_steps | The maximum number of substeps for the ODE integration. | Int | 150000
jacobian_type | Select how to compute the Jacobian for the ODE iterative solver. This option is available for `StiffSolver::BackwardEuler` and `StiffSolver::VODE`. Options (case-insensitive): `Numerical` - first-order numerical approximation `Broyden` - Broyden’s method (only for `BackwardEuler`) | String | Numerical
atol | Absolute tolerance for the ODE integration error test between the solution and the fine-step solution. | Real | 1.e-6
rtol | Relative tolerance for the ODE integration error test between the solution and the fine-step solution. | Real | 1.e-6
retry_atol | Overwrites the absolute tolerance value atol in case the ODE integration fails. | Real | -1
retry_rtol | Overwrites the relative tolerance value rtol in case the ODE integration fails. | Real | -1
retry_swap_jacobian | Swaps the type of Jacobian (from 1 to 0 or vice versa) in case the ODE integration fails. | Int | 1
ode_max_dt | Maximum timestep size for the ODE integration substeps. | Real | 1.e30
max_dt_change_factor | Only for StiffSolver::ForwardEuler, sets a maximum factor for the change of the timestep for the internal substeps. | Real | 1.001
max_iter | Maximum number of iterations for the internal Newton solver. | Int | 100
linalg_do_pivoting | Only for `StiffSolver::BackwardEuler` and `StiffSolver::VODE`, enables pivoting when solving the linear algebra problem associated to the internal Newton solver. | Int | 1
rp_rtol | Relative tolerance for the convergence test of the internal Newton solver. Valid only for `StiffSolver::BackwardEuler`. | Real | 1.e-6
use_jacobian_caching | Only for `StiffSolver::VODE`, enable caching the numerical Jacobian. | Int | 1
X_reject_buffer | Only for `StiffSolver::VODE`, constrain species abundances such that they don’t change by more than a certain factor in a given step. | Real | 1.0

Below is an example for specifying chemical reactions for MFIX-Exa.

# Example inputs

The following setup shows how to define chemical reactions in the inputs file. Only `Ideal Gas` constraint types can be used when solving chemistry, as the `IncompressibleFluid` constraint would be inconsistent with the nature of a problem where the fluid density is not constant because of mass transfer due to chemical reactions.

    mfix.constraint_type  = IdealGasOpenSystem

    mfix.advect_density   = 1
    mfix.advect_enthalpy  = 1
    mfix.solve_species    = 1

    # Species model settings
    #----------------------------------------------------------------------
    species.solve  =  CH4  H2  CO  H2O  CO2  O2  N2  Fe2O3  FeO  inert

    species.diffusivity = constant
    species.diffusivity.constant  =  2.31000e-5  # m^2 / sec

    species.specific_heat = NASA7-poly

    # Methane
    species.CH4.molecular_weight  = 16.04276e-3
    species.CH4.specific_heat.NASA7.a0 =  5.14911468E+00    1.65326226E+00
    species.CH4.specific_heat.NASA7.a1 = -1.36622009E-02    1.00263099E-02
    species.CH4.specific_heat.NASA7.a2 =  4.91453921E-05   -3.31661238E-06
    species.CH4.specific_heat.NASA7.a3 = -4.84246767E-08    5.36483138E-10
    species.CH4.specific_heat.NASA7.a4 =  1.66603441E-11   -3.14696758E-14
    species.CH4.specific_heat.NASA7.a5 = -1.02465983E+04   -1.00095936E+04

    # Hydrogen
    species.H2.molecular_weight = 2.01588e-3
    species.H2.specific_heat.NASA7.a0 =  2.34433112E+00    2.93286575E+00
    species.H2.specific_heat.NASA7.a1 =  7.98052075E-03    8.26608026E-04
    species.H2.specific_heat.NASA7.a2 = -1.94781510E-05   -1.46402364E-07
    species.H2.specific_heat.NASA7.a3 =  2.01572094E-08    1.54100414E-11
    species.H2.specific_heat.NASA7.a4 = -7.37611761E-12   -6.88804800E-16
    species.H2.specific_heat.NASA7.a5 = -9.17935173E+02   -8.13065581E+02

    # Carbon monoxide
    species.CO.molecular_weight = 28.01040e-3
    species.CO.specific_heat.NASA7.a0 =  3.57953350E+00    3.04848590E+00
    species.CO.specific_heat.NASA7.a1 = -6.10353690E-04    1.35172810E-03
    species.CO.specific_heat.NASA7.a2 =  1.01681430E-06   -4.85794050E-07
    species.CO.specific_heat.NASA7.a3 =  9.07005860E-10    7.88536440E-11
    species.CO.specific_heat.NASA7.a4 = -9.04424490E-13   -4.69807460E-15
    species.CO.specific_heat.NASA7.a5 = -1.43440860E+04   -1.42661170E+04

    # Water vapor
    species.H2O.molecular_weight = 18.01528e-3
    species.H2O.specific_heat.NASA7.a0 =  0.41986352E+01    0.26770389E+01
    species.H2O.specific_heat.NASA7.a1 = -0.20364017E-02    0.29731816E-02
    species.H2O.specific_heat.NASA7.a2 =  0.65203416E-05   -0.77376889E-06
    species.H2O.specific_heat.NASA7.a3 = -0.54879269E-08    0.94433514E-10
    species.H2O.specific_heat.NASA7.a4 =  0.17719680E-11   -0.42689991E-14
    species.H2O.specific_heat.NASA7.a5 = -0.30293726E+05   -0.29885894E+05

    # Carbon dioxide
    species.CO2.molecular_weight = 44.00980e-3
    species.CO2.specific_heat.NASA7.a0 =  2.35681300E+00    4.63651110E+00
    species.CO2.specific_heat.NASA7.a1 =  8.98412990E-03    2.74145690E-03
    species.CO2.specific_heat.NASA7.a2 = -7.12206320E-06   -9.95897590E-07
    species.CO2.specific_heat.NASA7.a3 =  2.45730080E-09    1.60386660E-10
    species.CO2.specific_heat.NASA7.a4 = -1.42885480E-13   -9.16198570E-15
    species.CO2.specific_heat.NASA7.a5 = -4.83719710E+04   -4.90249040E+04

    # Oxygen
    species.O2.molecular_weight = 31.99880e-3
    species.O2.specific_heat.NASA7.a0 =  3.78245636E+00    3.66096065E+00
    species.O2.specific_heat.NASA7.a1 = -2.99673416E-03    6.56365811E-04
    species.O2.specific_heat.NASA7.a2 =  9.84730201E-06   -1.41149627E-07
    species.O2.specific_heat.NASA7.a3 = -9.68129509E-09    2.05797935E-11
    species.O2.specific_heat.NASA7.a4 =  3.24372837E-12   -1.29913436E-15
    species.O2.specific_heat.NASA7.a5 = -1.06394356E+03   -1.21597718E+03

    # Nitrogen
    species.N2.molecular_weight = 28.01340e-3
    species.N2.specific_heat.NASA7.a0 =  3.53100528E+00    2.95257637E+00
    species.N2.specific_heat.NASA7.a1 = -1.23660988E-04    1.39690040E-03
    species.N2.specific_heat.NASA7.a2 = -5.02999433E-07   -4.92631603E-07
    species.N2.specific_heat.NASA7.a3 =  2.43530612E-09    7.86010195E-11
    species.N2.specific_heat.NASA7.a4 = -1.40881235E-12   -4.60755204E-15
    species.N2.specific_heat.NASA7.a5 = -1.04697628E+03   -9.23948688E+02

    # Hematite
    species.Fe2O3.molecular_weight = 159.68820e-3
    species.Fe2O3.specific_heat.NASA7.a0 =  1.52218166E-01    2.09445369E+01
    species.Fe2O3.specific_heat.NASA7.a1 =  6.70757040E-02    0.00000000E+00
    species.Fe2O3.specific_heat.NASA7.a2 = -1.12860954E-04    0.00000000E+00
    species.Fe2O3.specific_heat.NASA7.a3 =  9.93356662E-08    0.00000000E+00
    species.Fe2O3.specific_heat.NASA7.a4 = -3.27580975E-11    0.00000000E+00
    species.Fe2O3.specific_heat.NASA7.a5 = -1.01344092E+05   -1.07936580E+05

    # Wustite
    species.FeO.molecular_weight = 71.84440e-3
    species.FeO.specific_heat.NASA7.a0 =  3.68765953E+00    1.81588527E+00
    species.FeO.specific_heat.NASA7.a1 =  1.09133433E-02    1.70742829E-02
    species.FeO.specific_heat.NASA7.a2 = -1.61179493E-05   -2.39919190E-05
    species.FeO.specific_heat.NASA7.a3 =  1.06449256E-08    1.53690046E-08
    species.FeO.specific_heat.NASA7.a4 = -2.39514915E-12   -3.53442390E-12
    species.FeO.specific_heat.NASA7.a5 = -3.34867527E+04   -3.30239565E+04

    # inert
    species.inert.molecular_weight = 92.82399e-3
    species.inert.specific_heat.NASA7.a0 =  1.52218166E-01    2.09445369E+01
    species.inert.specific_heat.NASA7.a1 =  6.70757040E-02    0.00000000E+00
    species.inert.specific_heat.NASA7.a2 = -1.12860954E-04    0.00000000E+00
    species.inert.specific_heat.NASA7.a3 =  9.93356662E-08    0.00000000E+00
    species.inert.specific_heat.NASA7.a4 = -3.27580975E-11    0.00000000E+00
    species.inert.specific_heat.NASA7.a5 = -1.01344092E+05   -1.07936580E+05

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

    fluid.thermodynamic_pressure = 101325.
    fluid.reference_temperature = 298.15

    fluid.species =  CH4  H2  CO  H2O  CO2  O2  N2

    fluid.viscosity          = constant
    fluid.viscosity.constant = 4.73e-5  #air @ 1200K

    fluid.specific_heat = mixture

    fluid.thermal_conductivity          = constant
    fluid.thermal_conductivity.constant = 0.024

    # Solids model settings
    #----------------------------------------------------------------------
    solids.types = solid0
    solids.reference_temperature = 298.15

    solids.species  =  Fe2O3  FeO  inert

    solids.specific_heat = mixture

    # Chemistry model settings
    #----------------------------------------------------------------------
    chemistry.solve  =  Hem_CH4  Hem_H2  Hem_CO  Wus_O2

    chemistry.Hem_CH4.reaction = Fe2O3(s)+0.25CH4(g)-->2.FeO(s)+0.25CO2(g)+0.5H2O(g)
    chemistry.Hem_H2.reaction  = Fe2O3(s)+H2(g)-->2.FeO(s)+H2O(g)
    chemistry.Hem_CO.reaction  = Fe2O3(s)+CO(g)-->2.FeO(s)+CO2(g)
    chemistry.Wus_O2.reaction  = FeO(s)+0.25O2(g)-->0.5Fe2O3(s)

    # Initial Conditions
    #----------------------------------------------------------------------
    ic.regions = full-domain

    ic.full-domain.fluid.volfrac       = 0.725
    ic.full-domain.fluid.velocity      = 0.015  0.00  0.00
    ic.full-domain.fluid.temperature   = 1273.
    ic.full-domain.fluid.species.CH4   = 0.1
    ic.full-domain.fluid.species.H2    = 0.2
    ic.full-domain.fluid.species.CO    = 0.1
    ic.full-domain.fluid.species.H2O   = 0.2
    ic.full-domain.fluid.species.CO2   = 0.1
    ic.full-domain.fluid.species.O2    = 0.1
    ic.full-domain.fluid.species.N2    = 0.2

    ic.full-domain.solids  = solid0
    ic.full-domain.packing = pseudo_random

    ic.full-domain.solid0.volfrac  =  0.275
    ic.full-domain.solid0.velocity =  0.00  0.00  0.00

    ic.full-domain.solid0.diameter = constant
    ic.full-domain.solid0.density  = constant

    ic.full-domain.solid0.diameter.constant =  100.0e-6
    ic.full-domain.solid0.density.constant  = 2000.0

    ic.full-domain.solid0.temperature = 1273.
    ic.full-domain.solid0.species.Fe2O3 = 0.3
    ic.full-domain.solid0.species.FeO   = 0.3
    ic.full-domain.solid0.species.inert = 0.4

    # Boundary Conditions
    #----------------------------------------------------------------------
    bc.regions = inflow  outflow

    bc.inflow = mi
    bc.inflow.fluid.volfrac       = 1.0
    bc.inflow.fluid.velocity      = 0.015  0.0  0.0
    bc.inflow.fluid.temperature   = 1273.
    bc.inflow.fluid.species.CH4   = 0.1
    bc.inflow.fluid.species.H2    = 0.2
    bc.inflow.fluid.species.CO    = 0.1
    bc.inflow.fluid.species.H2O   = 0.2
    bc.inflow.fluid.species.CO2   = 0.1
    bc.inflow.fluid.species.O2    = 0.1
    bc.inflow.fluid.species.N2    = 0.2

    bc.outflow = po
    bc.outflow.fluid.pressure =  101325.