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/process/user_defined_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.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{ std::vector<micm::Species>{ a, b, c } };
micm::Process r1 = micm::Process::Create()
.SetReactants({ a })
.SetProducts({ Yields(b, 1) })
.SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r1" }))
.SetPhase(gas_phase);
micm::Process r2 = micm::Process::Create()
.SetReactants({ b, b })
.SetProducts({ Yields(b, 1), Yields(c, 1) })
.SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r2" }))
.SetPhase(gas_phase);
micm::Process r3 = micm::Process::Create()
.SetReactants({ b, c })
.SetProducts({ Yields(a, 1), Yields(c, 1) })
.SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r3" }))
.SetPhase(gas_phase);
const std::size_t number_of_grid_cells = 3;
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }))
.SetReactions({ r1, r2, r3 })
.SetNumberOfGridCells(number_of_grid_cells)
.Build();
auto state = solver.GetState();
// mol m-3
state.SetConcentration(a, std::vector<double>{ 1, 2, 0.5 });
state.SetConcentration(b, std::vector<double>(3, 0));
state.SetConcentration(c, std::vector<double>(3, 0));
double k1 = 0.04;
double k2 = 3e7;
double k3 = 1e4;
state.SetCustomRateParameter("r1", std::vector<double>(3, k1));
state.SetCustomRateParameter("r2", std::vector<double>(3, k2));
state.SetCustomRateParameter("r3", std::vector<double>(3, k3));
double temperature = 272.5; // [K]
double pressure = 101253.3; // [Pa]
double air_density = 1e6; // [mol m-3]
for (size_t 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
double time_step = 200; // s
state.PrintHeader();
state.PrintState(0);
// 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;
while (elapsed_solve_time < time_step)
{
solver.CalculateRateConstants(state);
auto result = solver.Solve(time_step - elapsed_solve_time, state);
elapsed_solve_time = result.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/process/user_defined_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>
After that, we’ll use the micm namespace.
#include <iostream>
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.
int main()
{
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::Phase gas_phase{ std::vector<micm::Species>{ a, b, c } };
micm::Process r1 = micm::Process::Create()
.SetReactants({ a })
.SetProducts({ Yields(b, 1) })
.SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r1" }))
.SetPhase(gas_phase);
micm::Process r2 = micm::Process::Create()
.SetReactants({ b, b })
.SetProducts({ Yields(b, 1), Yields(c, 1) })
.SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r2" }))
.SetPhase(gas_phase);
micm::Process r3 = micm::Process::Create()
.SetReactants({ b, c })
.SetProducts({ Yields(a, 1), Yields(c, 1) })
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.
.SetPhase(gas_phase);
const std::size_t number_of_grid_cells = 3;
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.
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }))
.SetReactions({ r1, r2, r3 })
.SetNumberOfGridCells(number_of_grid_cells)
.Build();
Then we set the reaction rates by creating a vector that is 3 elements long, one for each grid cell.
// mol m-3
state.SetConcentration(a, std::vector<double>{ 1, 2, 0.5 });
state.SetConcentration(b, std::vector<double>(3, 0));
state.SetConcentration(c, std::vector<double>(3, 0));
And lastly set the temperature, pressure, and air density for each grid cell.
double k2 = 3e7;
double k3 = 1e4;
state.SetCustomRateParameter("r1", std::vector<double>(3, k1));
state.SetCustomRateParameter("r2", std::vector<double>(3, k2));
state.SetCustomRateParameter("r3", std::vector<double>(3, k3));
double temperature = 272.5; // [K]
double pressure = 101253.3; // [Pa]
double air_density = 1e6; // [mol m-3]
Now we are ready to run the simulation.
{
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
double time_step = 200; // s
state.PrintHeader();
state.PrintState(0);
// 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;
while (elapsed_solve_time < time_step)
{
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 |