OpenMP#
This tutorial will focus on running micm with OpenMP support. 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.
If you’re looking for a copy and paste, copy below and be on your way! Otherwise, stick around for a line by line explanation.
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>
#include <omp.h>
using namespace micm;
void print_header()
{
std::cout << std::setw(10) << "A"
<< "," << std::setw(10) << "B"
<< "," << std::setw(10) << "C" << std::endl;
}
void print_results(std::vector<double> results)
{
std::ios oldState(nullptr);
oldState.copyfmt(std::cout);
std::cout << std::scientific << std::setprecision(2) << std::setw(10) << results[0] << "," << std::setw(10) << results[1]
<< "," << std::setw(10) << results[2] << std::endl;
std::cout.copyfmt(oldState);
}
std::vector<double> run_solver_on_thread_with_own_state(auto& solver, auto& state)
{
std::cout << "Running solver on thread " << omp_get_thread_num() << std::endl;
// 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;
double time_step = 200; // s
for (int i = 0; i < 10; ++i)
{
double elapsed_solve_time = 0;
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_;
}
}
return state.variables_.AsVector();
}
int main()
{
constexpr size_t n_threads = 3;
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);
auto reactions = std::vector<micm::Process>{ r1, r2, r3 };
auto chemical_system = micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase });
std::vector<std::vector<double>> results(n_threads);
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(chemical_system)
.SetReactions(reactions)
.Build();
#pragma omp parallel num_threads(n_threads)
{
// each thread should use its own state
auto state = solver.GetState();
std::vector<double> result = run_solver_on_thread_with_own_state(solver, state);
results[omp_get_thread_num()] = result;
#pragma omp barrier
}
std::cout << "Thread 1" << std::endl;
print_header();
print_results(results[0]);
std::cout << "Thread 2" << std::endl;
print_header();
print_results(results[1]);
std::cout << "Thread 3" << std::endl;
print_header();
print_results(results[2]);
}
Line-by-line explanation#
The Rosenbrock class is threadsafe. Each thread, however, must have its own unique micm::State object.
This is because the state object is modified during the solver’s run. The state object is passed to the solver’s
run method. The solver will modify the state object in place.
This tutorial createas a mechanism, one rosenbrock method, and then three states, each on their own thread, with that same configuration information.
First we need to bring in the necessary imports and use the micm namespace.
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>
#include <omp.h>
using namespace micm;
Then we’ll define a funtion that we can call with the configuration information from each thread. This function will build and run the rosenbrock solver.
std::vector<double> run_solver_on_thread_with_own_state(auto& solver, auto& state)
{
std::cout << "Running solver on thread " << omp_get_thread_num() << std::endl;
// 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;
double time_step = 200; // s
for (int i = 0; i < 10; ++i)
{
double elapsed_solve_time = 0;
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_;
}
}
return state.variables_.AsVector();
}
The main function simply creates the mechansim, sets the number of threads, and then sets up an OpenMP blcok with three threads. The function defined above is called on each thread.
int main()
{
constexpr size_t n_threads = 3;
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);
auto reactions = std::vector<micm::Process>{ r1, r2, r3 };
auto chemical_system = micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase });
std::vector<std::vector<double>> results(n_threads);
auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
.SetSystem(chemical_system)
.SetReactions(reactions)
.Build();
#pragma omp parallel num_threads(n_threads)
{
// each thread should use its own state
auto state = solver.GetState();
std::vector<double> result = run_solver_on_thread_with_own_state(solver, state);
results[omp_get_thread_num()] = result;
#pragma omp barrier
}
Running this program should give an output similar to this:
Thread 1
A, B, C
2.55e-01, 1.37e-06, 7.45e-01
Thread 2
A, B, C
2.55e-01, 1.37e-06, 7.45e-01
Thread 3
A, B, C
2.55e-01, 1.37e-06, 7.45e-01