Rate Constants (except user defined ones)#
MICM supports a subset of the rate constants defined as part of the OpenAtmos Mechanism Configuration We will be adding more in the future. The parameter structs below hold all configuration for each rate constant type. Please check the OpenAtmos standard for specifics on the algorithm implemented for each rate constant type. At present, supported rate constants are:
micm::ArrheniusRateConstantParameters, OpenAtmos Arrhenius descriptionmicm::BranchedRateConstantParameters, OpenAtmos Branched descriptionmicm::SurfaceRateConstantParameters, OpenAtmos Surface descriptionmicm::TernaryChemicalActivationRateConstantParameters, OpenAtmos Ternay chemical activiation descriptionmicm::TroeRateConstantParameters, OpenAtmos Troe descriptionmicm::TunnelingRateConstantParameters, OpenAtmos Tunneling descriptionmicm::UserDefinedRateConstantParameters, this is a custom type, but this is how photolysis, first order loss, and emission rates are setup OpenAtmos Photolysis description
This tutorial covers all but the last one. See the User-Defined Rate Constants tutorial for examples and use cases on that.
We’ll setup and solve a fake chemical system with 7 species and 6 reactions,
MICM can be configured in two ways. We can either build the mechanism up by hand with the micm API,
or parse a valid mechanism Configuration
in the OpenAtmos format. In this tutorial, we will do both.
If you’re looking for a copy and paste, choose the appropriate tab below and be on your way! Otherwise, stick around for a line by line explanation.
#include <micm/CPU.hpp>
#include <micm/util/types.hpp>
#include <iomanip>
#include <iostream>
#include <map>
// Use our namespace so that this example is easier to read
using namespace micm;
// Conversion factor from moles m-3 to molecules cm-3 for consistency
// with the configuraion file
constexpr Real MOLES_M3_TO_MOLECULES_CM3 = 1.0e-6 * constants::AVOGADRO_CONSTANT;
int main(const int argc, const char* argv[])
{
auto a = Species("A");
auto b = Species("B");
auto c = Species(
"C",
std::map<std::string, Real>{ { "molecular weight [kg mol-1]", 0.025 }});
auto d = Species("D");
auto e = Species("E");
auto f = Species("F");
auto g = Species("G");
Phase gas_phase{ "gas", std::vector<PhaseSpecies>{ a, b, c, d, e, f, g } };
auto& phase_species_list = gas_phase.phase_species_;
auto it = std::find_if(phase_species_list.begin(), phase_species_list.end(),
[&c](const PhaseSpecies& ps) {
return ps.species_.name_ == c.name_; });
if (it == phase_species_list.end())
{
std::cout << "Species not found\n";
return 1; // Failure
}
Real c_diffusion_coefficient = 2.3e2;
Index surface_c_index = std::distance(phase_species_list.begin(), it);
phase_species_list[surface_c_index].SetDiffusionCoefficient(c_diffusion_coefficient);
Process r1 = ChemicalReactionBuilder()
.SetReactants({ a })
.SetProducts({ StoichSpecies(b, 1) })
.SetRateConstant(ArrheniusRateConstantParameters{ .A_ = 2.15e-4, .B_ = 0, .C_ = 110 })
.SetPhase(gas_phase)
.Build();
// a branched reaction has two output pathways
// this is represnted internal to micm as two different reactions
auto branched_params = BranchedRateConstantParameters{ .X_ = 1.2, .Y_ = 204.3, .a0_ = 1.0e-3, .n_ = 2 };
branched_params.branch_ = BranchedRateConstantParameters::Branch::Alkoxy;
Process r2 = ChemicalReactionBuilder()
.SetReactants({ b })
.SetProducts({ StoichSpecies(c, 1) })
.SetRateConstant(branched_params)
.SetPhase(gas_phase)
.Build();
branched_params.branch_ = BranchedRateConstantParameters::Branch::Nitrate;
Process r3 = ChemicalReactionBuilder()
.SetReactants({ b })
.SetProducts({ StoichSpecies(d, 1) })
.SetRateConstant(branched_params)
.SetPhase(gas_phase)
.Build();
// A surface rate constant also needs to know the effective radius and particle number concentration
// we will set those later
Process r4 = ChemicalReactionBuilder()
.SetReactants({ c })
.SetProducts({ StoichSpecies(e, 1) })
.SetRateConstant(SurfaceRateConstantParameters{ .label_ = "C", .phase_species_ = phase_species_list[surface_c_index], .reaction_probability_ = 0.90 })
.SetPhase(gas_phase)
.Build();
Process r5 = ChemicalReactionBuilder()
.SetReactants({ d })
.SetProducts({ StoichSpecies(f, 2) })
.SetRateConstant(TernaryChemicalActivationRateConstantParameters{ .k0_A_ = 1.2,
.k0_B_ = 2.3,
.k0_C_ = 302.3,
.kinf_A_ = 2.6 / MOLES_M3_TO_MOLECULES_CM3,
.kinf_B_ = -3.1,
.kinf_C_ = 402.1,
.Fc_ = 0.9,
.N_ = 1.2 })
.SetPhase(gas_phase)
.Build();
// to have a stoichiemetric coefficient of more than one for reactants,
// list the reactant that many times
//
// The falloff parameters are those measured for OH + NO2 + M -> HNO3 + M (JPL 19-5, Table 2):
// k0 = 1.8e-30 (T/300)^-3.0 cm6 molecule-2 s-1 and kinf = 2.8e-11 cm3 molecule-1 s-1. Using
// realistic magnitudes matters here -- k0_A_ is a third-order rate constant, so the conversion
// to mole-based units squares MOLES_M3_TO_MOLECULES_CM3, and an unphysically large k0 overflows
// single precision.
Process r6 =
ChemicalReactionBuilder()
.SetReactants({ e, e })
.SetProducts({ StoichSpecies(g, 1) })
.SetRateConstant(TroeRateConstantParameters{ .k0_A_ = 1.8e-30 * MOLES_M3_TO_MOLECULES_CM3 * MOLES_M3_TO_MOLECULES_CM3,
.k0_B_ = -3.0,
.k0_C_ = 0.0,
.kinf_A_ = 2.8e-11 * MOLES_M3_TO_MOLECULES_CM3,
.kinf_B_ = 0.0,
.kinf_C_ = 0.0,
.Fc_ = 0.6,
.N_ = 1.0 })
.SetPhase(gas_phase)
.Build();
Process r7 = ChemicalReactionBuilder()
.SetReactants({ f })
.SetProducts({ StoichSpecies(g, 1) })
.SetRateConstant(TunnelingRateConstantParameters{ .A_ = 1.2, .B_ = 2.3, .C_ = 302.3 })
.SetPhase(gas_phase)
.Build();
auto chemical_system = System(gas_phase );
auto reactions = std::vector<micm::Process>{ r1, r2, r3, r4, r5, r6, r7 };
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(
micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(chemical_system)
.SetReactions(reactions)
.Build();
State state = solver.GetState();
state.conditions_[0].temperature_ = 287.45; // K
state.conditions_[0].pressure_ = 101319.9; // Pa
state.conditions_[0].CalculateIdealAirDensity();
state.SetConcentration(a, 1.0); // mol m-3
state.SetConcentration(b, 0.0); // mol m-3
state.SetConcentration(c, 0.0); // mol m-3
state.SetConcentration(d, 0.0); // mol m-3
state.SetConcentration(e, 0.0); // mol m-3
state.SetConcentration(f, 0.0); // mol m-3
state.SetConcentration(g, 0.0); // mol m-3
state.SetCustomRateParameter("C.effective radius [m]", 1e-7);
state.SetCustomRateParameter("C.particle number concentration [# m-3]", 2.5e6);
// choose a timestep and print the initial state
Real time_step = 500; // s
state.PrintHeader();
state.PrintState(0);
// solve for ten iterations
for (Index i = 0; i < 10; ++i)
{
// Depending on how stiff the system is
// the solver integration step may not be able to solve for the full time step
// so we need to track how much time the solver was able to integrate for and continue
// solving until we finish
Real elapsed_solve_time = 0;
solver.UpdateStateParameters(state);
while (elapsed_solve_time < time_step)
{
auto result = solver.Solve(time_step - elapsed_solve_time, state);
// A solve that advances no time never satisfies the loop condition, so bail out rather than
// spin forever on a system the solver cannot integrate
if (result.stats_.final_time_ <= 0)
{
std::cerr << "Solver made no progress: " << SolverStateToString(result.state_) << std::endl;
return 1;
}
elapsed_solve_time += result.stats_.final_time_;
}
state.PrintState(time_step * (i + 1));
}
return 0;
}
Line-by-line explanation#
To get started, we’ll need to include each rate constant type and the rosenbrock solver at the top of the file.
#include <micm/CPU.hpp> #include <micm/util/types.hpp> #include <iomanip> #include <iostream> #include <map> // Use our namespace so that this example is easier to read using namespace micm; // Conversion factor from moles m-3 to molecules cm-3 for consistency // with the configuraion file constexpr Real MOLES_M3_TO_MOLECULES_CM3 = 1.0e-6 * constants::AVOGADRO_CONSTANT;
After that, we’ll use the micm namespace so that we don’t have to repeat it everywhere we need it.
{ auto a = Species("A");
To create a micm::RosenbrockSolver, we have to define a chemical system (micm::System)
and our reactions, which will be a vector of micm::Process We will use the species to define these.
To do this by hand, we have to define all of the chemical species in the system. This allows us to set any properties of the species that may be necessary for rate constanta calculations, like molecular weights and diffusion coefficients for the surface reaction. We will also put these species into the gas phase.
auto g = Species("G");
Phase gas_phase{ "gas", std::vector<PhaseSpecies>{ a, b, c, d, e, f, g } };
auto& phase_species_list = gas_phase.phase_species_;
auto it = std::find_if(phase_species_list.begin(), phase_species_list.end(),
[&c](const PhaseSpecies& ps) {
return ps.species_.name_ == c.name_; });
if (it == phase_species_list.end())
{
std::cout << "Species not found\n";
return 1; // Failure
Now that we have a gas phase and our species, we can start building the reactions. Two things to note are that
stoichiemtric coefficients for reactants are represented by repeating that product as many times as you need.
To specify the yield of a product, we’ve created a struct micm::Yield. Note that we add a conversion
for some rate constant parameters to be consistent with the configuration file that expects rate constants
to be in cm^3/molecule/s. (All units will be mks in the next version of the configuration file format.)
Real c_diffusion_coefficient = 2.3e2;
Index surface_c_index = std::distance(phase_species_list.begin(), it);
phase_species_list[surface_c_index].SetDiffusionCoefficient(c_diffusion_coefficient);
Process r1 = ChemicalReactionBuilder()
.SetReactants({ a })
.SetProducts({ StoichSpecies(b, 1) })
.SetRateConstant(ArrheniusRateConstantParameters{ .A_ = 2.15e-4, .B_ = 0, .C_ = 110 })
.SetPhase(gas_phase)
.Build();
// a branched reaction has two output pathways
// this is represnted internal to micm as two different reactions
auto branched_params = BranchedRateConstantParameters{ .X_ = 1.2, .Y_ = 204.3, .a0_ = 1.0e-3, .n_ = 2 };
branched_params.branch_ = BranchedRateConstantParameters::Branch::Alkoxy;
Process r2 = ChemicalReactionBuilder()
.SetReactants({ b })
.SetProducts({ StoichSpecies(c, 1) })
.SetRateConstant(branched_params)
.SetPhase(gas_phase)
.Build();
branched_params.branch_ = BranchedRateConstantParameters::Branch::Nitrate;
Process r3 = ChemicalReactionBuilder()
.SetReactants({ b })
.SetProducts({ StoichSpecies(d, 1) })
.SetRateConstant(branched_params)
.SetPhase(gas_phase)
.Build();
// A surface rate constant also needs to know the effective radius and particle number concentration
// we will set those later
Process r4 = ChemicalReactionBuilder()
.SetReactants({ c })
.SetProducts({ StoichSpecies(e, 1) })
.SetRateConstant(SurfaceRateConstantParameters{ .label_ = "C", .phase_species_ = phase_species_list[surface_c_index], .reaction_probability_ = 0.90 })
.SetPhase(gas_phase)
.Build();
Process r5 = ChemicalReactionBuilder()
.SetReactants({ d })
.SetProducts({ StoichSpecies(f, 2) })
.SetRateConstant(TernaryChemicalActivationRateConstantParameters{ .k0_A_ = 1.2,
.k0_B_ = 2.3,
.k0_C_ = 302.3,
.kinf_A_ = 2.6 / MOLES_M3_TO_MOLECULES_CM3,
.kinf_B_ = -3.1,
.kinf_C_ = 402.1,
.Fc_ = 0.9,
.N_ = 1.2 })
.SetPhase(gas_phase)
.Build();
// to have a stoichiemetric coefficient of more than one for reactants,
// list the reactant that many times
//
// The falloff parameters are those measured for OH + NO2 + M -> HNO3 + M (JPL 19-5, Table 2):
// k0 = 1.8e-30 (T/300)^-3.0 cm6 molecule-2 s-1 and kinf = 2.8e-11 cm3 molecule-1 s-1. Using
// realistic magnitudes matters here -- k0_A_ is a third-order rate constant, so the conversion
// to mole-based units squares MOLES_M3_TO_MOLECULES_CM3, and an unphysically large k0 overflows
// single precision.
Process r6 =
ChemicalReactionBuilder()
And finally we define our chemical system and reactions
.SetProducts({ StoichSpecies(g, 1) })
.SetRateConstant(TroeRateConstantParameters{ .k0_A_ = 1.8e-30 * MOLES_M3_TO_MOLECULES_CM3 * MOLES_M3_TO_MOLECULES_CM3,
Now that we have a chemical system and a list of reactions, we can create the RosenbrockSolver.
There are several ways to configure the solver. Here we are using a three stage solver. More options
can be found in the micm::RosenbrockSolverParameters and in the Solver configurations tutorial.
.k0_C_ = 0.0, .kinf_A_ = 2.8e-11 * MOLES_M3_TO_MOLECULES_CM3, .kinf_B_ = 0.0, .kinf_C_ = 0.0,
The rosenbrock solver will provide us a state, which we can use to set the concentrations, custom rate parameters, and temperature and pressure. Note that setting the custom rate paramters is different depending on if you define the configuration by hand or read it in. The parser has defaults for the names of the custom parameters and when defined by hand we choose these.
Initializing the state#
.Fc_ = 0.6,
.N_ = 1.0 })
.SetPhase(gas_phase)
.Build();
Process r7 = ChemicalReactionBuilder()
.SetReactants({ f })
.SetProducts({ StoichSpecies(g, 1) })
.SetRateConstant(TunnelingRateConstantParameters{ .A_ = 1.2, .B_ = 2.3, .C_ = 302.3 })
.SetPhase(gas_phase)
.Build();
auto chemical_system = System(gas_phase );
auto reactions = std::vector<micm::Process>{ r1, r2, r3, r4, r5, r6, r7 };
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(
Finally, we are ready to pick a timestep and solve the system.
.SetReactions(reactions) .Build(); State state = solver.GetState(); state.conditions_[0].temperature_ = 287.45; // K state.conditions_[0].pressure_ = 101319.9; // Pa state.conditions_[0].CalculateIdealAirDensity(); state.SetConcentration(a, 1.0); // mol m-3 state.SetConcentration(b, 0.0); // mol m-3 state.SetConcentration(c, 0.0); // mol m-3 state.SetConcentration(d, 0.0); // mol m-3 state.SetConcentration(e, 0.0); // mol m-3 state.SetConcentration(f, 0.0); // mol m-3 state.SetConcentration(g, 0.0); // mol m-3 state.SetCustomRateParameter("C.effective radius [m]", 1e-7); state.SetCustomRateParameter("C.particle number concentration [# m-3]", 2.5e6); // choose a timestep and print the initial state Real time_step = 500; // s state.PrintHeader();
This is the output:
0 |
1.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
500 |
8.54e-01 |
4.57e-04 |
1.44e-01 |
1.55e-04 |
0.00e+00 |
1.23e-22 |
6.44e-04 |
1000 |
7.30e-01 |
3.90e-04 |
2.65e-01 |
2.89e-04 |
8.18e-07 |
2.28e-22 |
2.44e-03 |
1500 |
6.23e-01 |
3.33e-04 |
3.66e-01 |
4.02e-04 |
9.66e-07 |
3.18e-22 |
5.20e-03 |
2000 |
5.32e-01 |
2.85e-04 |
4.49e-01 |
5.00e-04 |
1.07e-06 |
3.95e-22 |
8.77e-03 |
2500 |
4.55e-01 |
2.43e-04 |
5.18e-01 |
5.83e-04 |
1.15e-06 |
4.61e-22 |
1.30e-02 |
3000 |
3.88e-01 |
2.08e-04 |
5.75e-01 |
6.54e-04 |
1.21e-06 |
5.17e-22 |
1.78e-02 |
3500 |
3.32e-01 |
1.77e-04 |
6.21e-01 |
7.14e-04 |
1.26e-06 |
5.65e-22 |
2.30e-02 |
4000 |
2.83e-01 |
1.52e-04 |
6.59e-01 |
7.66e-04 |
1.30e-06 |
6.06e-22 |
2.86e-02 |
4500 |
2.42e-01 |
1.29e-04 |
6.88e-01 |
8.10e-04 |
1.33e-06 |
6.41e-22 |
3.45e-02 |
5000 |
2.07e-01 |
1.11e-04 |
7.11e-01 |
8.48e-04 |
1.35e-06 |
6.71e-22 |
4.06e-02 |