User-Defined Rate Constants#
This tutorial extends the Rate Constants (except user defined ones) tutorial. That one showed how to use every type of rate constant except the user-defined ones. The difference is that a user-defined rate constant must be updated by the user, whereas the other rate constants update in response to a changing state.
Internal to MICM, user-defined rate constants provide the ability to represent processes like emissions, loss, and photolysis that have rate constants that are calculated externally from MICM. These can be static and the same for each time step, or dynamically updated during the simulation.
To show how this works, we’ll add one photolysis reaction, one first-order loss reaction, and one emission
\[\begin{split}A &\longrightarrow B, &k_{1, \mathrm{arrhenius}} \\
B &\longrightarrow C (\mathrm{alkoxy\ products}) + D (\mathrm{nitrate\ products}), &k_{2, \mathrm{branched}} \\
C &\longrightarrow E, &k_{3, \mathrm{surface}} \\
D &\longrightarrow 2F, &k_{4, \mathrm{ternary\ chemical\ activation}} \\
2E &\longrightarrow G, &k_{5, \mathrm{troe}} \\
F &\longrightarrow G, &k_{6, \mathrm{tunneling}} \\
C &\longrightarrow G, &k_{7, \mathrm{photolysis}} \\
&\longrightarrow A, &k_{8, \mathrm{emission}} \\
B &\longrightarrow, &k_{9, \mathrm{loss}} \\\end{split}\]
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.
// Each rate constant is in its own header file
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/process/branched_rate_constant.hpp>
#include <micm/process/surface_rate_constant.hpp>
#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
#include <micm/process/troe_rate_constant.hpp>
#include <micm/process/tunneling_rate_constant.hpp>
#include <micm/process/user_defined_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>
#include <iomanip>
#include <iostream>
#include <map>
// Use our namespace so that this example is easier to read
using namespace micm;
int main(const int argc, const char* argv[])
{
auto a = Species("A");
auto b = Species("B");
auto c = Species(
"C",
std::map<std::string, double>{ { "molecular weight [kg mol-1]", 0.025 },
{ "diffusion coefficient [m2 s-1]", 2.3e2 } });
auto d = Species("D");
auto e = Species("E");
auto f = Species("F");
auto g = Species("G");
Phase gas_phase{ std::vector<Species>{ a, b, c, d, e, f, g } };
Process r1 = Process::Create()
.SetReactants({ a })
.SetProducts({ Yields(b, 1) })
.SetRateConstant(ArrheniusRateConstant({ .A_ = 2.15e-1, .B_ = 0, .C_ = 110 }))
.SetPhase(gas_phase);
// 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 = Process::Create()
.SetReactants({ b })
.SetProducts({ Yields(c, 1) })
.SetRateConstant(BranchedRateConstant(branched_params))
.SetPhase(gas_phase);
branched_params.branch_ = BranchedRateConstantParameters::Branch::Nitrate;
Process r3 = Process::Create()
.SetReactants({ b })
.SetProducts({ Yields(d, 1) })
.SetRateConstant(BranchedRateConstant(branched_params))
.SetPhase(gas_phase);
// A surface rate constant also needs to know the effective radius and particle number concentration
// we will set those later
Process r4 = Process::Create()
.SetReactants({ c })
.SetProducts({ Yields(e, 1) })
.SetRateConstant(SurfaceRateConstant({ .label_ = "C", .species_ = c, .reaction_probability_ = 0.90 }))
.SetPhase(gas_phase);
Process r5 = Process::Create()
.SetReactants({ d })
.SetProducts({ Yields(f, 2) })
.SetRateConstant(TernaryChemicalActivationRateConstant({ .k0_A_ = 1.2,
.k0_B_ = 2.3,
.k0_C_ = 302.3,
.kinf_A_ = 2.6,
.kinf_B_ = -3.1,
.kinf_C_ = 402.1,
.Fc_ = 0.9,
.N_ = 1.2 }))
.SetPhase(gas_phase);
// to have a stoichiemetric coefficient of more than one for reactants,
// list the reactant that many times
Process r6 = Process::Create()
.SetReactants({ e, e })
.SetProducts({ Yields(g, 1) })
.SetRateConstant(TroeRateConstant({ .k0_A_ = 1.2e4,
.k0_B_ = 167.0,
.k0_C_ = 3.0,
.kinf_A_ = 136.0,
.kinf_B_ = 5.0,
.kinf_C_ = 24.0,
.Fc_ = 0.9,
.N_ = 0.8 }))
.SetPhase(gas_phase);
Process r7 = Process::Create()
.SetReactants({ f })
.SetProducts({ Yields(g, 1) })
.SetRateConstant(TunnelingRateConstant({ .A_ = 1.2, .B_ = 2.3, .C_ = 302.3 }))
.SetPhase(gas_phase);
Process r8 = Process::Create()
.SetReactants({ c })
.SetProducts({ Yields(g, 1) })
.SetRateConstant(UserDefinedRateConstant({ .label_ = "my photolysis rate" }))
.SetPhase(gas_phase);
Process r9 = Process::Create()
.SetProducts({ Yields(a, 1) })
.SetRateConstant(UserDefinedRateConstant({ .label_ = "my emission rate" }))
.SetPhase(gas_phase);
Process r10 = Process::Create()
.SetReactants({ b })
.SetRateConstant(UserDefinedRateConstant({ .label_ = "my loss rate" }))
.SetPhase(gas_phase);
auto chemical_system = System(micm::SystemParameters{ .gas_phase_ = gas_phase });
auto reactions = std::vector<micm::Process>{ r1, r2, r3, r4, r5, r6, r7, r8, r9, r10 };
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.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 and timestep a print the initial state
double time_step = 500; // s
state.PrintHeader();
state.PrintState(0);
double photo_rate = 1e-10;
double emission_rate = 1e-20;
double loss = emission_rate * 1e-3;
// these rates are constant through the simulation
state.SetCustomRateParameter("my emission rate", emission_rate);
state.SetCustomRateParameter("my loss rate", loss);
// solve for ten iterations
for (int 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
double elapsed_solve_time = 0;
// this rate is updated at each time step and would typically vary with time
state.SetCustomRateParameter("my photolysis rate", photo_rate);
solver.CalculateRateConstants(state);
while (elapsed_solve_time < time_step)
{
auto result = solver.Solve(time_step - elapsed_solve_time, state);
elapsed_solve_time = result.final_time_;
}
state.PrintState(time_step * (i + 1));
photo_rate *= 1.5;
}
return 0;
}
// Each rate constant is in its own header file
#include <micm/configure/solver_config.hpp>
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/process/branched_rate_constant.hpp>
#include <micm/process/surface_rate_constant.hpp>
#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
#include <micm/process/troe_rate_constant.hpp>
#include <micm/process/tunneling_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>
#include <iomanip>
#include <iostream>
#include <map>
// Use our namespace so that this example is easier to read
using namespace micm;
int main(const int argc, const char* argv[])
{
SolverConfig solverConfig;
std::string config_path = "./configs/rate_constants_user_defined";
try
{
solverConfig.ReadAndParse(config_path);
}
catch (const std::system_error& e)
{
std::cerr << "Error reading and parsing config file: " << e.what() << std::endl;
return 1;
}
micm::SolverParameters solver_params = solverConfig.GetSolverParams();
auto chemical_system = solver_params.system_;
auto reactions = solver_params.processes_;
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
std::unordered_map<std::string, std::vector<double>> intial_concentration = {
{ "A", { 1.0 } }, // mol m-3
{ "B", { 0.0 } }, // mol m-3
{ "C", { 0.0 } }, // mol m-3
{ "D", { 0.0 } }, // mol m-3
{ "E", { 0.0 } }, // mol m-3
{ "F", { 0.0 } }, // mol m-3
{ "G", { 0.0 } }, // mol m-3
};
state.SetConcentrations(intial_concentration);
state.SetCustomRateParameter("SURF.C surface.effective radius [m]", 1e-7);
state.SetCustomRateParameter("SURF.C surface.particle number concentration [# m-3]", 2.5e6);
// choose a timestep a print the initial state
double time_step = 500; // s
state.PrintHeader();
state.PrintState(0);
double photo_rate = 1e-10;
double emission_rate = 1e-20;
double loss = emission_rate * 1e-3;
// these rates are constant through the simulation
state.SetCustomRateParameter("EMIS.my emission rate", emission_rate);
state.SetCustomRateParameter("LOSS.my loss rate", loss);
// solve for ten iterations
for (int 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
double elapsed_solve_time = 0;
// This rate is updated at each time step and would typically vary with time
state.SetCustomRateParameter("PHOTO.my photolysis rate", photo_rate);
solver.CalculateRateConstants(state);
while (elapsed_solve_time < time_step)
{
auto result = solver.Solve(time_step - elapsed_solve_time, state);
elapsed_solve_time = result.final_time_;
}
state.PrintState(time_step * (i + 1));
photo_rate *= 1.5;
}
return 0;
}
Line-by-line explanation#
Adding the custom rate constant is quite simple. Include the header file:
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/process/branched_rate_constant.hpp>
#include <micm/process/surface_rate_constant.hpp>
#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
#include <micm/process/troe_rate_constant.hpp>
#include <micm/process/tunneling_rate_constant.hpp>
+ #include <micm/process/user_defined_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
Then setup the reaction which will use this rate constant:
Process r7 = Process::Create()
.SetReactants({ f })
.SetProducts({ Yields(g, 1) })
.SetRateConstant(TunnelingRateConstant({ .A_ = 1.2, .B_ = 2.3, .C_ = 302.3 }))
.SetPhase(gas_phase);
+ Process r8 = Process::Create()
+ .SetReactants({ c })
+ .SetProducts({ Yields(g, 1) })
+ .SetRateConstant(UserDefinedRateConstant({.label_="my rate"}))
+ .SetPhase(gas_phase);
+ Process r9 = Process::Create()
+ .SetProducts({ Yields(a, 1) })
+ .SetRateConstant(UserDefinedRateConstant({.label_="my emission rate"}))
+ .SetPhase(gas_phase);
+ Process r10 = Process::Create()
+ .SetReactants({ b })
+ .SetRateConstant(UserDefinedRateConstant({.label_="my loss rate"}))
+ .SetPhase(gas_phase);
auto chemical_system = System(micm::SystemParameters{ .gas_phase_ = gas_phase });
- auto reactions = std::vector<micm::Process>{ r1, r2, r3, r4, r5, r6, r7 };
+ auto reactions = std::vector<micm::Process>{ r1, r2, r3, r4, r5, r6, r7, r8, r9, r10 };
In this case, you only need to add the configuration to the reactions.json file in the configuration directory.
+ {
+ "type": "PHOTOLYSIS",
+ "reactants": {
+ "C": {}
+ },
+ "products": {
+ "G": {}
+ },
+ "MUSICA name": "my photolysis rate"
+ },
+ {
+ "type": "FIRST_ORDER_LOSS",
+ "species": "B",
+ "MUSICA name": "my loss rate"
+ },
+ {
+ "type": "EMISSION",
+ "species": "A",
+ "MUSICA name": "my emission rate"
+ }
Finally, set and upate the rate constants as needed:
+ double photo_rate = 1e-10;
+ double emission_rate = 1e-20;
+ double loss = emission_rate * 1e-3;
+ // these rates are constant through the simulation
+ state.SetCustomRateParameter("my emission rate", emission_rate);
+ state.SetCustomRateParameter("my loss rate", loss);
// solve for ten iterations
for (int 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
double elapsed_solve_time = 0;
+ state.SetCustomRateParameter("my photolysis rate", photo_rate);
while (elapsed_solve_time < time_step)
{
auto result = solver.Solve(time_step - elapsed_solve_time, state);
elapsed_solve_time = result.final_time_;
state.variables_[0] = result.result_.AsVector();
}
print_state(time_step * (i + 1), state);
+ photo_rate *= 1.5;
}
In this case, you only need to add the configuration to the reactions.json file in the configuration directory.
When reading in from a configuration file, the loss, emissions, and photolysis rates are prefixed with
LOSS., EMIS., and PHOTO.. This differs slightly from defining the API by hand.
+ double photo_rate = 1e-10;
+ double emission_rate = 1e-20;
+ double loss = emission_rate * 1e-3;
+ // these rates are constant through the simulation
+ state.SetCustomRateParameter("EMIS.my emission rate", emission_rate);
+ state.SetCustomRateParameter("LOSS.my loss rate", loss);
// solve for ten iterations
for (int 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
double elapsed_solve_time = 0;
+ state.SetCustomRateParameter("PHOTO.my photolysis rate", photo_rate);
while (elapsed_solve_time < time_step)
{
auto result = solver.Solve(time_step - elapsed_solve_time, state);
elapsed_solve_time = result.final_time_;
state.variables_[0] = result.result_.AsVector();
}
print_state(time_step * (i + 1), state);
+ photo_rate *= 1.5;
}
And this is final output. Notice that the concentration of G ends up much higher than in the Rate Constants (except user defined ones) tutorial’s result.
time |
A |
B |
C |
D |
E |
F |
G |
|---|---|---|---|---|---|---|---|
0 |
1.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
0.00e+00 |
500 |
3.18e-09 |
3.66e-09 |
9.83e-01 |
3.88e-14 |
1.41e-03 |
2.02e-13 |
7.92e-03 |
1000 |
1.14e-14 |
1.31e-14 |
9.66e-01 |
1.39e-19 |
1.40e-03 |
7.24e-19 |
1.64e-02 |
1500 |
7.27e-20 |
6.40e-20 |
9.49e-01 |
6.53e-25 |
1.39e-03 |
3.19e-24 |
2.48e-02 |
2000 |
3.17e-20 |
1.70e-20 |
9.33e-01 |
1.55e-25 |
1.38e-03 |
5.92e-25 |
3.30e-02 |
2500 |
3.17e-20 |
1.70e-20 |
9.17e-01 |
1.55e-25 |
1.37e-03 |
5.92e-25 |
4.11e-02 |
3000 |
3.17e-20 |
1.70e-20 |
9.01e-01 |
1.55e-25 |
1.36e-03 |
5.92e-25 |
4.90e-02 |
3500 |
3.17e-20 |
1.70e-20 |
8.85e-01 |
1.55e-25 |
1.34e-03 |
5.92e-25 |
5.68e-02 |
4000 |
3.17e-20 |
1.70e-20 |
8.70e-01 |
1.55e-25 |
1.33e-03 |
5.92e-25 |
6.44e-02 |
4500 |
3.17e-20 |
1.70e-20 |
8.55e-01 |
1.55e-25 |
1.32e-03 |
5.92e-25 |
7.20e-02 |
5000 |
3.17e-20 |
1.70e-20 |
8.40e-01 |
1.55e-25 |
1.31e-03 |
5.92e-25 |
7.94e-02 |