Solver configurations#
This tutorial will focus on configuring the solver with different stages. We will use a simple 3-reaction 3-species mechanism. The setup here is the same in Multiple Grid Cells, except only one grid cell is used.
\[\begin{split}A &\longrightarrow B, &k_{1, \mathrm{user\ defined}} \\
2B &\longrightarrow B + C, &k_{2, \mathrm{user\ defined}} \\
B + C &\longrightarrow A + C, \qquad &k_{3, \mathrm{user\ defined}} \\\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.
#include <micm/process/user_defined_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>
#include <chrono>
#include <iomanip>
#include <iostream>
// Use our namespace so that this example is easier to read
using namespace micm;
void test_solver_type(auto& solver)
{
auto state = solver.GetState();
// mol m-3
state.variables_[0] = { 1, 0, 0 };
double k1 = 0.04;
double k2 = 3e7;
double k3 = 1e4;
state.SetCustomRateParameter("r1", k1);
state.SetCustomRateParameter("r2", k2);
state.SetCustomRateParameter("r3", k3);
double temperature = 272.5; // [K]
double pressure = 101253.3; // [Pa]
double air_density = 1e6; // [mol m-3]
state.conditions_[0].temperature_ = temperature;
state.conditions_[0].pressure_ = pressure;
state.conditions_[0].air_density_ = air_density;
// choose a timestep and print the initial state
double time_step = 200; // s
state.PrintHeader();
state.PrintState(0);
SolverStats total_stats;
std::chrono::duration<double, std::nano> total_solve_time = std::chrono::nanoseconds::zero();
// 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)
{
auto start = std::chrono::high_resolution_clock::now();
solver.CalculateRateConstants(state);
auto result = solver.Solve(time_step - elapsed_solve_time, state);
auto end = std::chrono::high_resolution_clock::now();
total_solve_time += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start);
total_stats.function_calls_ += result.stats_.function_calls_;
total_stats.jacobian_updates_ += result.stats_.jacobian_updates_;
total_stats.number_of_steps_ += result.stats_.number_of_steps_;
total_stats.accepted_ += result.stats_.accepted_;
total_stats.rejected_ += result.stats_.rejected_;
total_stats.decompositions_ += result.stats_.decompositions_;
total_stats.solves_ += result.stats_.solves_;
total_stats.singular_ += result.stats_.singular_;
elapsed_solve_time = result.final_time_;
}
state.PrintState(time_step * (i + 1));
}
std::cout << "Total solve time: " << total_solve_time.count() << " nanoseconds" << std::endl;
std::cout << "accepted: " << total_stats.accepted_ << std::endl;
std::cout << "function_calls: " << total_stats.function_calls_ << std::endl;
std::cout << "jacobian_updates: " << total_stats.jacobian_updates_ << std::endl;
std::cout << "number_of_steps: " << total_stats.number_of_steps_ << std::endl;
std::cout << "accepted: " << total_stats.accepted_ << std::endl;
std::cout << "rejected: " << total_stats.rejected_ << std::endl;
std::cout << "decompositions: " << total_stats.decompositions_ << std::endl;
std::cout << "solves: " << total_stats.solves_ << std::endl;
std::cout << "singular: " << total_stats.singular_ << std::endl;
}
int main()
{
auto a = Species("A");
auto b = Species("B");
auto c = Species("C");
Phase gas_phase{ std::vector<Species>{ a, b, c } };
Process r1 = Process::Create()
.SetReactants({ a })
.SetProducts({ Yields(b, 1) })
.SetRateConstant(UserDefinedRateConstant({ .label_ = "r1" }))
.SetPhase(gas_phase);
Process r2 = Process::Create()
.SetReactants({ b, b })
.SetProducts({ Yields(b, 1), Yields(c, 1) })
.SetRateConstant(UserDefinedRateConstant({ .label_ = "r2" }))
.SetPhase(gas_phase);
Process r3 = Process::Create()
.SetReactants({ b, c })
.SetProducts({ Yields(a, 1), Yields(c, 1) })
.SetRateConstant(UserDefinedRateConstant({ .label_ = "r3" }))
.SetPhase(gas_phase);
auto system = System(SystemParameters{ .gas_phase_ = gas_phase });
auto reactions = std::vector<Process>{ r1, r2, r3 };
auto two_stage = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters())
.SetSystem(system)
.SetReactions(reactions)
.Build();
auto three_stage = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(system)
.SetReactions(reactions)
.Build();
auto four_stage = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters())
.SetSystem(system)
.SetReactions(reactions)
.Build();
auto four_stage_da = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters())
.SetSystem(system)
.SetReactions(reactions)
.Build();
auto six_stage_da = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters())
.SetSystem(system)
.SetReactions(reactions)
.Build();
std::cout << "Two stages: " << std::endl;
test_solver_type(two_stage);
std::cout << std::endl << "Three stages: " << std::endl;
test_solver_type(three_stage);
std::cout << std::endl << "Four stages: " << std::endl;
test_solver_type(four_stage);
std::cout << std::endl << "Four stages differential algebraic: " << std::endl;
test_solver_type(four_stage_da);
std::cout << std::endl << "Six stages differential algebraic: " << std::endl;
test_solver_type(six_stage_da);
}
Line-by-line explanation#
There are a total of five different configurations for the Rosenbrock solver. Each corresponds to a different number of stages. What each configuration means is beyond the scope of this tutorial. However, references are included for further reading.
- Two stage
An L-stable method, 2 stages, order 2. [HW96]
- Three stage
An L-stable method, 3 stages, order 3. [SVB+97]
- Four stage
An L-stable method, 4 stages, order 4. [HW96]
- Four stage, differential algebraic
A stiffly stable method, 4 stages, order 3, [HW96]
- Six stage, differential algebraic
Stiffly stable rosenbrock method of order 4 with 6 stages, [HW96]
Configuring the rosenbrock solver is as easy as providing the solver with a set of micm::RosenbrockSolverParameters
.SetSystem(system)
.SetReactions(reactions)
.Build();
auto four_stage = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters())
.SetSystem(system)
.SetReactions(reactions)
.Build();
auto four_stage_da = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters())
.SetSystem(system)
.SetReactions(reactions)
.Build();
After that, the usage is the same as a regular solver. A templated method was used here to run the same solving code for each of the different solver configurations.
using namespace micm;
void test_solver_type(auto& solver)
{
auto state = solver.GetState();
// mol m-3
state.variables_[0] = { 1, 0, 0 };
double k1 = 0.04;
double k2 = 3e7;
double k3 = 1e4;
state.SetCustomRateParameter("r1", k1);
state.SetCustomRateParameter("r2", k2);
state.SetCustomRateParameter("r3", k3);
double temperature = 272.5; // [K]
double pressure = 101253.3; // [Pa]
double air_density = 1e6; // [mol m-3]
state.conditions_[0].temperature_ = temperature;
state.conditions_[0].pressure_ = pressure;
state.conditions_[0].air_density_ = air_density;
// choose a timestep and print the initial state
double time_step = 200; // s
state.PrintHeader();
state.PrintState(0);
SolverStats total_stats;
std::chrono::duration<double, std::nano> total_solve_time = std::chrono::nanoseconds::zero();
// 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)
{
auto start = std::chrono::high_resolution_clock::now();
solver.CalculateRateConstants(state);
auto result = solver.Solve(time_step - elapsed_solve_time, state);
auto end = std::chrono::high_resolution_clock::now();
total_solve_time += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start);
total_stats.function_calls_ += result.stats_.function_calls_;
total_stats.jacobian_updates_ += result.stats_.jacobian_updates_;
total_stats.number_of_steps_ += result.stats_.number_of_steps_;
total_stats.accepted_ += result.stats_.accepted_;
total_stats.rejected_ += result.stats_.rejected_;
total_stats.decompositions_ += result.stats_.decompositions_;
total_stats.solves_ += result.stats_.solves_;
total_stats.singular_ += result.stats_.singular_;
elapsed_solve_time = result.final_time_;
}
state.PrintState(time_step * (i + 1));
}
std::cout << "Total solve time: " << total_solve_time.count() << " nanoseconds" << std::endl;
std::cout << "accepted: " << total_stats.accepted_ << std::endl;
std::cout << "function_calls: " << total_stats.function_calls_ << std::endl;
std::cout << "jacobian_updates: " << total_stats.jacobian_updates_ << std::endl;
std::cout << "number_of_steps: " << total_stats.number_of_steps_ << std::endl;
std::cout << "accepted: " << total_stats.accepted_ << std::endl;
std::cout << "rejected: " << total_stats.rejected_ << std::endl;
std::cout << "decompositions: " << total_stats.decompositions_ << std::endl;
std::cout << "solves: " << total_stats.solves_ << std::endl;
std::cout << "singular: " << total_stats.singular_ << std::endl;
}
int main()
{
auto a = Species("A");
auto b = Species("B");
auto c = Species("C");
Running this program should give an output similar to this:
Two stages:
time, A, B, C
0, 1.00e+00, 0.00e+00, 0.00e+00
200, 5.37e-01, 4.50e-06, 4.63e-01
400, 4.51e-01, 3.23e-06, 5.49e-01
600, 4.01e-01, 2.64e-06, 5.99e-01
800, 3.65e-01, 2.28e-06, 6.35e-01
1000, 3.38e-01, 2.02e-06, 6.62e-01
1200, 3.16e-01, 1.83e-06, 6.84e-01
1400, 2.97e-01, 1.68e-06, 7.03e-01
1600, 2.82e-01, 1.56e-06, 7.18e-01
1800, 2.68e-01, 1.46e-06, 7.32e-01
2000, 2.56e-01, 1.37e-06, 7.44e-01
Total solve time: 135748 nanoseconds
accepted: 178
function_calls: 356
jacobian_updates: 178
number_of_steps: 178
accepted: 178
rejected: 0
decompositions: 178
solves: 356
singular: 0
Three stages:
time, A, B, C
0, 1.00e+00, 0.00e+00, 0.00e+00
200, 5.35e-01, 4.50e-06, 4.65e-01
400, 4.50e-01, 3.23e-06, 5.50e-01
600, 4.00e-01, 2.63e-06, 6.00e-01
800, 3.64e-01, 2.27e-06, 6.36e-01
1000, 3.37e-01, 2.01e-06, 6.63e-01
1200, 3.15e-01, 1.82e-06, 6.85e-01
1400, 2.96e-01, 1.67e-06, 7.04e-01
1600, 2.81e-01, 1.55e-06, 7.19e-01
1800, 2.67e-01, 1.45e-06, 7.33e-01
2000, 2.55e-01, 1.37e-06, 7.45e-01
Total solve time: 110002 nanoseconds
accepted: 127
function_calls: 254
jacobian_updates: 127
number_of_steps: 127
accepted: 127
rejected: 0
decompositions: 127
solves: 381
singular: 0
Four stages:
time, A, B, C
0, 1.00e+00, 0.00e+00, 0.00e+00
200, 5.36e-01, 4.49e-06, 4.64e-01
400, 4.50e-01, 3.20e-06, 5.50e-01
600, 4.00e-01, 2.62e-06, 6.00e-01
800, 3.64e-01, 2.26e-06, 6.36e-01
1000, 3.37e-01, 2.01e-06, 6.63e-01
1200, 3.15e-01, 1.82e-06, 6.85e-01
1400, 2.97e-01, 1.67e-06, 7.03e-01
1600, 2.81e-01, 1.55e-06, 7.19e-01
1800, 2.67e-01, 1.45e-06, 7.33e-01
2000, 2.56e-01, 1.37e-06, 7.44e-01
Total solve time: 127291 nanoseconds
accepted: 127
function_calls: 381
jacobian_updates: 127
number_of_steps: 127
accepted: 127
rejected: 0
decompositions: 127
solves: 508
singular: 0
Four stages differential algebraic:
time, A, B, C
0, 1.00e+00, 0.00e+00, 0.00e+00
200, 5.36e-01, 4.49e-06, 4.64e-01
400, 4.51e-01, 3.23e-06, 5.49e-01
600, 4.00e-01, 2.63e-06, 6.00e-01
800, 3.64e-01, 2.27e-06, 6.36e-01
1000, 3.37e-01, 2.01e-06, 6.63e-01
1200, 3.15e-01, 1.82e-06, 6.85e-01
1400, 2.97e-01, 1.68e-06, 7.03e-01
1600, 2.81e-01, 1.55e-06, 7.19e-01
1800, 2.68e-01, 1.45e-06, 7.32e-01
2000, 2.56e-01, 1.37e-06, 7.44e-01
Total solve time: 126334 nanoseconds
accepted: 128
function_calls: 384
jacobian_updates: 128
number_of_steps: 128
accepted: 128
rejected: 0
decompositions: 128
solves: 512
singular: 0
Six stages differential algebraic:
time, A, B, C
0, 1.00e+00, 0.00e+00, 0.00e+00
200, 5.36e-01, 4.49e-06, 4.64e-01
400, 4.51e-01, 3.22e-06, 5.49e-01
600, 4.00e-01, 2.63e-06, 6.00e-01
800, 3.64e-01, 2.27e-06, 6.36e-01
1000, 3.37e-01, 2.01e-06, 6.63e-01
1200, 3.15e-01, 1.82e-06, 6.85e-01
1400, 2.97e-01, 1.67e-06, 7.03e-01
1600, 2.81e-01, 1.55e-06, 7.19e-01
1800, 2.67e-01, 1.45e-06, 7.33e-01
2000, 2.56e-01, 1.37e-06, 7.44e-01
Total solve time: 174959 nanoseconds
accepted: 127
function_calls: 762
jacobian_updates: 127
number_of_steps: 127
accepted: 127
rejected: 0
decompositions: 127
solves: 762
singular: 0