Solver results#

This tutorial will focus on solver result. Solver result is made of a collection of information useful to understand convergence of solution. In a 3D model, if certain grid cells fail the solver convergence criteraia, solver statistics can be used to diagnose strange behavior. We will use a simple 3-reaction 3-species mechanism. The setup here is the same in Multiple Grid Cells. To understand the full setup, read that tutorial. Otherwise, we assume that configuring a rosenbrock solver is understood and instead we will focus on timing the solver.

A⟶B,k1,user defined2B⟶B+C,k2,user definedB+C⟶A+C,k3,user defined

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;

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

  const std::size_t number_of_grid_cells = 3;

  auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
                    .SetSystem(System(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;
  }
  solver.CalculateRateConstants(state);

  // choose a timestep and print the initial state
  double time_step = 200;  // s

  auto result = solver.Solve(time_step, state);
  std::cout << "Solver state: " << SolverStateToString(result.state_) << std::endl;
  std::cout << "accepted: " << result.stats_.accepted_ << std::endl;
  std::cout << "function_calls: " << result.stats_.function_calls_ << std::endl;
  std::cout << "jacobian_updates: " << result.stats_.jacobian_updates_ << std::endl;
  std::cout << "number_of_steps: " << result.stats_.number_of_steps_ << std::endl;
  std::cout << "accepted: " << result.stats_.accepted_ << std::endl;
  std::cout << "rejected: " << result.stats_.rejected_ << std::endl;
  std::cout << "decompositions: " << result.stats_.decompositions_ << std::endl;
  std::cout << "solves: " << result.stats_.solves_ << std::endl;
  std::cout << "final simulation time: " << result.final_time_ << std::endl;
}

Line-by-line explanation#

Up until now we have neglected to talk about what the solver returns, which is a micm::RosenbrockSolver::SolverResult.

There are four values returned.

  1. micm::RosenbrockSolver::SolverResult::final_time_

    • This is the final simulation time achieved by the solver. The micm::RosenbrockSolver::Solve() function attempts to integrate the passed in state forward a set number of seconds. Often, the solver is able to complete the integration. However, extremely stiff systems may only solve for a fraction of the time. It is imperative that the final_time_ value is checked. If it is not equal to the amount of time you intended to solve for, call Solve again as we do in the tutorials with the difference between what was solved and how long you intended to solve.

      Note

      This does not represent the amount of time taken by the solve routine.

  2. micm::RosenbrockSolver::SolverResult::result_

    • This contains the integrated state value; the concentrations reached at the end of Solve function after the amount of time specified by final_time_.

  3. micm::RosenbrockSolver::SolverResult::state_

    • There are many possible reasons for the solver to return. This value is one of the possible enum values defined on the micm::SolverState. Hopefully, you receive a micm::SolverState::Converged state. But, it is good to always check this to ensure the solver really did converge. You can print this value using the micm::StateToString() function.

  4. micm::RosenbrockSolver::SolverResult::stats_

    • This is an instance of a micm::RosenbrockSolver::SolverStats struct which contains information about the number of individual function calls during solving process and the number of accepted or rejected solutions at every time step.

Upon completion of the simulation, we can inspect the solver state and statistics including the number of individual function calls, acceptance/rejection of solutions and the existence of singular matrix.

  auto result = solver.Solve(time_step, state);
  std::cout << "Solver state: " << SolverStateToString(result.state_) << std::endl;
  std::cout << "accepted: " << result.stats_.accepted_ << std::endl;
  std::cout << "function_calls: " << result.stats_.function_calls_ << std::endl;
  std::cout << "jacobian_updates: " << result.stats_.jacobian_updates_ << std::endl;
  std::cout << "number_of_steps: " << result.stats_.number_of_steps_ << std::endl;
  std::cout << "accepted: " << result.stats_.accepted_ << std::endl;
  std::cout << "rejected: " << result.stats_.rejected_ << std::endl;
  std::cout << "decompositions: " << result.stats_.decompositions_ << std::endl;
  std::cout << "solves: " << result.stats_.solves_ << std::endl;
  std::cout << "final simulation time: " << result.final_time_ << std::endl;

Solver state: Converged
accepted: 20
function_calls: 40
jacobian_updates: 20
number_of_steps: 20
accepted: 20
rejected: 0
decompositions: 20
solves: 60
singular: 0
On this page
Show Source