Multiple Grid Cells#
This tutorial will focus on running multiple grid cells. Because the Rate Constants (except user defined ones) and User-Defined Rate Constants showed both configuring the mechanism by hand and building with a configuration file, we will only show building up the mechanism by hand for this tutorial. We will use a simple 3-reaction 3-species mechanism.
We will use three grid cells. The second grid cells will have concentrations twice as large as the first grid cell. The third grid cell will have concentrations half as large as the first grid cell.
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>
// Use our namespace so that this example is easier to read
using namespace micm;
int main()
{
auto a = micm::Species("A");
auto b = micm::Species("B");
auto c = micm::Species("C");
micm::Phase gas_phase{ "gas", std::vector<micm::PhaseSpecies>{ a, b, c } };
micm::Process r1 = micm::ChemicalReactionBuilder()
.SetReactants({ a })
.SetProducts({ micm::StoichSpecies(b, 1) })
.SetRateConstant(micm::UserDefinedRateConstantParameters{ .label_ = "r1" })
.SetPhase(gas_phase)
.Build();
micm::Process r2 = micm::ChemicalReactionBuilder()
.SetReactants({ b, b })
.SetProducts({ micm::StoichSpecies(b, 1), micm::StoichSpecies(c, 1) })
.SetRateConstant(micm::UserDefinedRateConstantParameters{ .label_ = "r2" })
.SetPhase(gas_phase)
.Build();
micm::Process r3 = micm::ChemicalReactionBuilder()
.SetReactants({ b, c })
.SetProducts({ micm::StoichSpecies(a, 1), micm::StoichSpecies(c, 1) })
.SetRateConstant(micm::UserDefinedRateConstantParameters{ .label_ = "r3" })
.SetPhase(gas_phase)
.Build();
const Index number_of_grid_cells = 3;
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(
micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(micm::System(gas_phase ))
.SetReactions({ r1, r2, r3 })
.Build();
auto state = solver.GetState(number_of_grid_cells);
// mol m-3
state.SetConcentration(a, std::vector<Real>{ 1, 2, 0.5 });
state.SetConcentration(b, std::vector<Real>(3, 0));
state.SetConcentration(c, std::vector<Real>(3, 0));
Real k1 = 0.04;
Real k2 = 3e7;
Real k3 = 1e4;
state.SetCustomRateParameter("r1", std::vector<Real>(3, k1));
state.SetCustomRateParameter("r2", std::vector<Real>(3, k2));
state.SetCustomRateParameter("r3", std::vector<Real>(3, k3));
Real temperature = 272.5; // [K]
Real pressure = 101253.3; // [Pa]
Real air_density = 1e6; // [mol m-3]
for (Index cell = 0; cell < number_of_grid_cells; ++cell)
{
state.conditions_[cell].temperature_ = temperature;
state.conditions_[cell].pressure_ = pressure;
state.conditions_[cell].air_density_ = air_density;
}
// choose a timestep and print the initial state
Real time_step = 200; // 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;
while (elapsed_solve_time < time_step)
{
solver.UpdateStateParameters(state);
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));
}
}
Line-by-line explanation#
This mechanism only needs the user defined rate constant and the rosenbrock solver.
#include <micm/CPU.hpp>
#include <micm/util/types.hpp>
#include <iomanip>
After that, we’ll use the micm namespace.
// Use our namespace so that this example is easier to read
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 as
well as a micm::Phase.
{
auto a = micm::Species("A");
auto b = micm::Species("B");
auto c = micm::Species("C");
With the species and gas phase, we can define all of our reactions
micm::Process r1 = micm::ChemicalReactionBuilder()
.SetReactants({ a })
.SetProducts({ micm::StoichSpecies(b, 1) })
.SetRateConstant(micm::UserDefinedRateConstantParameters{ .label_ = "r1" })
.SetPhase(gas_phase)
.Build();
micm::Process r2 = micm::ChemicalReactionBuilder()
.SetReactants({ b, b })
.SetProducts({ micm::StoichSpecies(b, 1), micm::StoichSpecies(c, 1) })
.SetRateConstant(micm::UserDefinedRateConstantParameters{ .label_ = "r2" })
.SetPhase(gas_phase)
.Build();
micm::Process r3 = micm::ChemicalReactionBuilder()
.SetReactants({ b, c })
.SetProducts({ micm::StoichSpecies(a, 1), micm::StoichSpecies(c, 1) })
.SetRateConstant(micm::UserDefinedRateConstantParameters{ .label_ = "r3" })
Now we can define our RosenbrockSolver. This time we’ll form the reactions and chemical system in place.
Also, notice the false in our micm::RosenbrockSolverParameters. This tells the solver
not to reorder the state variables. The reordering is an optimization that can minizie fill-in for some
of the linear algebra operations. For this example, it’s turned off so that the order of the state matches
the order the species are added to the gas phase.
const Index number_of_grid_cells = 3;
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(
micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(micm::System(gas_phase ))
Now we need to get a state and set the concentations of each of the species. In the
“Initializing the state” section of the rate constants tutorial,
we used a std::unordered_map<std::string, std::vector<double>>
to set the concentrations. Here we will set the concentations by providing the micm::Species objects.
.Build();
auto state = solver.GetState(number_of_grid_cells);
// mol m-3
state.SetConcentration(a, std::vector<Real>{ 1, 2, 0.5 });
Then we set the reaction rates by creating a vector that is 3 elements long, one for each grid cell.
state.SetConcentration(c, std::vector<Real>(3, 0));
Real k1 = 0.04;
Real k2 = 3e7;
Real k3 = 1e4;
state.SetCustomRateParameter("r1", std::vector<Real>(3, k1));
And lastly set the temperature, pressure, and air density for each grid cell.
state.SetCustomRateParameter("r3", std::vector<Real>(3, k3));
Real temperature = 272.5; // [K]
Real pressure = 101253.3; // [Pa]
Real air_density = 1e6; // [mol m-3]
for (Index cell = 0; cell < number_of_grid_cells; ++cell)
{
state.conditions_[cell].temperature_ = temperature;
state.conditions_[cell].pressure_ = pressure;
Now we are ready to run the simulation.
}
// choose a timestep and print the initial state
Real time_step = 200; // 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;
while (elapsed_solve_time < time_step)
{
solver.UpdateStateParameters(state);
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_;
}
And these are the results.
time |
grid |
A |
B |
C |
|---|---|---|---|---|
0 |
0 |
1.00e+00 |
0.00e+00 |
0.00e+00 |
0 |
1 |
2.00e+00 |
0.00e+00 |
0.00e+00 |
0 |
2 |
5.00e-01 |
0.00e+00 |
0.00e+00 |
200 |
0 |
5.35e-01 |
4.49e-06 |
4.65e-01 |
200 |
1 |
1.21e+00 |
6.01e-06 |
7.89e-01 |
200 |
2 |
2.30e-01 |
3.31e-06 |
2.70e-01 |
400 |
0 |
4.50e-01 |
3.23e-06 |
5.50e-01 |
400 |
1 |
1.05e+00 |
4.42e-06 |
9.45e-01 |
400 |
2 |
1.85e-01 |
2.32e-06 |
3.15e-01 |
600 |
0 |
4.00e-01 |
2.63e-06 |
6.00e-01 |
600 |
1 |
9.59e-01 |
3.65e-06 |
1.04e+00 |
600 |
2 |
1.60e-01 |
1.85e-06 |
3.40e-01 |
800 |
0 |
3.64e-01 |
2.27e-06 |
6.36e-01 |
800 |
1 |
8.90e-01 |
3.18e-06 |
1.11e+00 |
800 |
2 |
1.42e-01 |
1.57e-06 |
3.58e-01 |
1000 |
0 |
3.37e-01 |
2.01e-06 |
6.63e-01 |
1000 |
1 |
8.35e-01 |
2.85e-06 |
1.16e+00 |
1000 |
2 |
1.29e-01 |
1.38e-06 |
3.71e-01 |
1200 |
0 |
3.15e-01 |
1.82e-06 |
6.85e-01 |
1200 |
1 |
7.91e-01 |
2.60e-06 |
1.21e+00 |
1200 |
2 |
1.19e-01 |
1.24e-06 |
3.81e-01 |
1400 |
0 |
2.96e-01 |
1.67e-06 |
7.04e-01 |
1400 |
1 |
7.54e-01 |
2.40e-06 |
1.25e+00 |
1400 |
2 |
1.11e-01 |
1.13e-06 |
3.89e-01 |
1600 |
0 |
2.81e-01 |
1.55e-06 |
7.19e-01 |
1600 |
1 |
7.21e-01 |
2.24e-06 |
1.28e+00 |
1600 |
2 |
1.04e-01 |
1.04e-06 |
3.96e-01 |
1800 |
0 |
2.67e-01 |
1.45e-06 |
7.33e-01 |
1800 |
1 |
6.93e-01 |
2.11e-06 |
1.31e+00 |
1800 |
2 |
9.77e-02 |
9.65e-07 |
4.02e-01 |
2000 |
0 |
2.55e-01 |
1.37e-06 |
7.45e-01 |
2000 |
1 |
6.68e-01 |
2.00e-06 |
1.33e+00 |
2000 |
2 |
9.25e-02 |
9.02e-07 |
4.07e-01 |