.. _InputsChem:

Defining chemical reactions
===========================

.. |tol_eq| replace:: :math:`\left|\sum \mathrm{products} - \sum \mathrm{reactants}\right| < \mathrm{tolerance}`

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

+--------------------------+---------------------------------------------------------+----------+----------------+
|                          | Description                                             |   Type   | Default        |
+==========================+=========================================================+==========+================+
| solve                    | Specify the name(s) of the chemical reactions or set    |  Strings |                |
|                          | to None to disable the chemistry solver. The name(s)    |          |                |
|                          | assigned to the chemistry solver are used to specify    |          |                |
|                          | the chemical reactions equations.                       |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| [reaction_name].reaction | Chemical formula for the given reaction. The string     |  String  |                |
|                          | given as input must not contain white space and         |          |                |
|                          | the reaction direction must be specified as `\-\->`     |          |                |
|                          | or `<\-\-`. Chemical species phases must be defined as  |          |                |
|                          | '(g)' for the fluid phase or '(s)' for the solid phase. |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| mass_balance_tolerance   | Tolerance used to test chemical formula conserves mass. | Real     | 1e-12          |
|                          |                                                         |          |                |
|                          | |tol_eq|                                                |          |                |
|                          |                                                         |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| solids.update_type       | How to update solids quantities due to chemical         | String   | ConstantVolume |
|                          | reactions.                                              |          |                |
|                          |                                                         |          |                |
|                          | Options:                                                |          |                |
|                          |                                                         |          |                |
|                          | * ``ConstantVolume`` - keep particle volume constant    |          |                |
|                          | * ``ConstantDensity`` - keep particle density constant  |          |                |
|                          | * ``Constant_Species_Densities`` - compute volume and   |          |                |
|                          |   density from                                          |          |                |
|                          |   :ref:`solid species densities<InputsSpeciesDefs>`.    |          |                |
|                          |   Requires :ref:`solids density model<InputsSolids>`    |          |                |
|                          |   to be set to ``mixture``.                             |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| solids.mass_threshold    | Threshold value for the particle mass.  When            | Real     | 0              |
|                          | updating the solids quantities in chemistry, if         |          |                |
|                          | the particle solids mass is below the threshold,        |          |                |
|                          | the particle is marked as invalid and removed.          |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| solids.radius_threshold  | Threshold value for the particle radius.  When          | Real     | 0              |
|                          | updating the solids quantities in chemistry, if         |          |                |
|                          | the particle radius is below the threshold,             |          |                |
|                          | the particle is marked as invalid and removed.          |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| solids.density_threshold | Threshold value for the particle density. When          | Real     | 0              |
|                          | updating the solids quantities in chemistry, if         |          |                |
|                          | the solids density is below the threshold,              |          |                |
|                          | the particle is marked as invalid and removed.          |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+

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         | String   | ForwardEuler  |
|                         | (case-insensitive):                                                  |          |               |
|                         |                                                                      |          |               |
|                         | * ``ForwardEuler``                                                   |          |               |
|                         | * ``StiffSolver::ForwardEuler``                                      |          |               |
|                         | * ``StiffSolver::BackwardEuler``                                     |          |               |
|                         | * ``StiffSolver::VODE``                                              |          |               |
|                         |                                                                      |          |               |
+-------------------------+----------------------------------------------------------------------+----------+---------------+

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.    | String | Numerical |
|                        | 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``)       |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| broyden_type           | When using the Broyden method in the Backward Euler ODE integrator, | String | Jacobian  |
|                        | select whether to approximate the Jacobian or its inverse.          |        |           |
|                        | This option is available only for ``StiffSolver::BackwardEuler``.   |        |           |
|                        |                                                                     |        |           |
|                        | Options (case-insensitive):                                         |        |           |
|                        |                                                                     |        |           |
|                        | * ``Jacobian``                                                      |        |           |
|                        | * ``Inverse``                                                       |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| broyden_initialization | When using the Broyden method in the Backward Euler ODE integrator, | String | Identity  |
|                        | select whether to initialize the Jacobian or its inverse through the|        |           |
|                        | identity matrix or the numerical (finite-difference) approximation. |        |           |
|                        |                                                                     |        |           |
|                        | Options (case-insensitive):                                         |        |           |
|                        |                                                                     |        |           |
|                        | * ``Identity``                                                      |        |           |
|                        | * ``Numerical``                                                     |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| atol                   | Absolute tolerance for the ODE integration error test between       | Real   | 1.e-6     |
|                        | the solution and the fine-step solution.                            |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| rtol                   | Relative tolerance for the ODE integration error test between       | Real   | 1.e-6     |
|                        | the solution and the fine-step solution.                            |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| retry_atol             | Overwrites the absolute tolerance value atol in case the ODE        | Real   | -1        |
|                        | integration fails.                                                  |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| retry_rtol             | Overwrites the relative tolerance value rtol in case the ODE        | Real   | -1        |
|                        | integration fails.                                                  |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| retry_swap_jacobian    | Swaps the type of Jacobian (from 1 to 0 or vice versa) in           | Int    | 1         |
|                        | case the ODE integration fails.                                     |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| 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           | Real   | 1.001     |
|                        | for the change of the timestep for the internal substeps.           |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| max_iter               | Maximum number of iterations for the internal Newton solver.        | Int    | 100       |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| linalg_do_pivoting     | Only for ``StiffSolver::BackwardEuler`` and ``StiffSolver::VODE``,  | Int    | 1         |
|                        | enables pivoting when solving the linear algebra                    |        |           |
|                        | problem associated to the internal Newton solver.                   |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| rp_rtol                | Relative tolerance for the convergence test of the internal         | Real   | 1.e-6     |
|                        | Newton solver. Valid only for ``StiffSolver::BackwardEuler``.       |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| use_jacobian_caching   | Only for ``StiffSolver::VODE``, enable caching the numerical        | Int    | 1         |
|                        | Jacobian.                                                           |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+
| X_reject_buffer        | Only for ``StiffSolver::VODE``, constrain species abundances        | Real   | 1.0       |
|                        | such that they don't change by more than a certain factor in        |        |           |
|                        | a given step.                                                       |        |           |
+------------------------+---------------------------------------------------------------------+--------+-----------+


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 ``IdealGas`` 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.

.. code-block:: bash

  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.model     = constant
  species.diffusivity.constant  =  2.31000e-5  # m^2 / sec

  species.specific_heat.model = 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.molecular.model    = constant
  fluid.viscosity.molecular.constant = 4.73e-5  #air @ 1200K

  fluid.specific_heat = mixture

  fluid.thermal_conductivity.model    = 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.model = 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.distribution = constant
  ic.full-domain.solid0.density.distribution  = 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.