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.

\[\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, copy below and be on your way! Otherwise, stick around for a line by line explanation.

OpenAtmos Configuration reading

Download zip configuration

#include <micm/configure/solver_config.hpp>
#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("PHOTO.r1", k1);
  state.SetCustomRateParameter("PHOTO.r2", k2);
  state.SetCustomRateParameter("PHOTO.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;

  SolverConfig solverConfig;

  std::string config_path = "./configs/robertson";
  solverConfig.ReadAndParse(config_path);

  micm::SolverParameters solver_params = solverConfig.GetSolverParams();

  auto chemical_system = solver_params.system_;
  auto reactions = solver_params.processes_;

  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

At present, the Rosenbrock class is not thread safe. Thread safety is being added and will be available soon in new release.

Until then, you can run one instance of the Rosenbrock solver on its own thread. Configuration data, at least, can be shared across the threads. This tutorial reads a configuraiton file and then configures three solvers, 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/configure/solver_config.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>

#include <omp.h>

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("PHOTO.r1", k1);
  state.SetCustomRateParameter("PHOTO.r2", k2);
  state.SetCustomRateParameter("PHOTO.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 reads the configuration file, 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;

  SolverConfig solverConfig;

  std::string config_path = "./configs/robertson";
  solverConfig.ReadAndParse(config_path);

  micm::SolverParameters solver_params = solverConfig.GetSolverParams();

  auto chemical_system = solver_params.system_;
  auto reactions = solver_params.processes_;

  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]);
}

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