Merge pull request #270 from NCAR/264-add-tutorial-on-timing-the-solv…

…er-and-investigating-solver-statistics

264 add tutorial on timing the solver and investigating solver statistics

docs/CMakeLists.txt

@@ -23,7 +23,7 @@ configure_file(${DOXYFILE_IN} ${DOXYFILE_OUT} @ONLY)
 file(MAKE_DIRECTORY ${DOXYGEN_OUTPUT_DIR}) #Doxygen won't create this for us
 
 add_custom_command(
-  OUTPUT ${DOXYGEN_INDEX_FILE}
+  OUTPUT ${DOXYGEN_INDEX_FILE} __fake_file_to_ensure_this_always_run
   DEPENDS ${PUBLIC_HEADERS}
   COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYFILE_OUT}
   MAIN_DEPENDENCY ${DOXYFILE_OUT} ${DOXYFILE_IN}

docs/source/user_guide/but_how_fast_is_it.rst

@@ -0,0 +1,90 @@
+.. _But how fast is it:
+
+But how fast is it?
+===================
+
+This tutorial will focus on timing the solver to show how you can measure performance.
+We will use a simple 3-reaction 3-species mechanism. The setup here is the same in
+:ref:`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.
+
+.. math::
+
+  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}} \\
+
+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.
+
+.. tabs::
+
+    .. tab:: Build the Mechanism with the API
+
+      .. literalinclude:: ../../../test/tutorial/test_but_how_fast_is_it.cpp
+        :language: cpp
+
+Line-by-line explanation
+------------------------
+
+Up until now we have neglected to talk about what the solver returns, which is a :cpp:class:`micm::RosenbrockSolver::SolverResult`.
+
+There are four values returned.
+
+#. :cpp:member:`micm::RosenbrockSolver::SolverResult::final_time_`
+
+    * This is the final simulation time achieved by the solver. The :cpp:func:`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. You must measure that yourself(shown below). The ``final_time_`` is simulation time.
+
+#. :cpp:member:`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_``.
+
+#. :cpp:member:`micm::RosenbrockSolver::SolverResult::state_`
+
+    * There are many possible reasons for the solver to return. This value is one of the possible enum values define on the :cpp:enum:`micm::SolverState`. Hopefully, you receive a :cpp:enumerator:`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 :cpp:func:`micm::StateToString` function.
+
+#. :cpp:member:`micm::RosenbrockSolver::SolverResult::stats_`
+
+    * This is an instance of a :cpp:class:`micm::RosenbrockSolver::SolverStats` struct which contains information about the number of function calls and optionally the total cumulative amount of time spent calling each function. For the time to be collected, you must call the ``Solve`` function with a ``true`` templated parameter. Please see the example below.
+
+First, let's run the simulation but without collecting the solve time. We'll inspect the solver state and look at what's collected
+in the stats object. 
+
+.. literalinclude:: ../../../test/tutorial/test_but_how_fast_is_it.cpp
+  :language: cpp
+  :lines: 75-86
+
+.. code-block:: console
+
+  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
+
+To get the total accumulated time of each function call, you need to specify the templated boolean argument to turn the timing on.
+We can also record the total runtime of the ``Solve`` function.
+
+.. literalinclude:: ../../../test/tutorial/test_but_how_fast_is_it.cpp
+  :language: cpp
+  :lines: 88-96
+
+.. code-block:: console
+
+  Total solve time: 24416 nanoseconds
+  total_forcing_time: 3167 nanoseconds
+  total_jacobian_time: 1710 nanoseconds
+  total_linear_factor_time: 4584 nanoseconds
+  total_linear_solve_time: 3290 nanoseconds
+
+.. note::
+  Your systems clock may not report the same values depending on how accurate your system clock is. 

docs/source/user_guide/index.rst

@@ -43,4 +43,5 @@ If you would like to include the json examples, you must configure micm to build
    rate_constant_tutorial
    user_defined_rate_constant_tutorial
    multiple_grid_cells
+   but_how_fast_is_it
    openmp

include/micm/solver/rosenbrock.hpp

@@ -129,12 +129,19 @@ namespace micm
   /// @brief The final state the solver was in after the Solve function finishes
   enum class SolverState
   {
+    /// @brief This is the initial value at the start of the Solve function
     NotYetCalled,
+    /// @brief This is only used for control flow in the Solve function
     Running,
+    /// @brief A successful integration will have this value
     Converged,
+    /// @brief If the number of steps exceeds the maximum value on the solver parameter, this value will be returned
     ConvergenceExceededMaxSteps,
+    /// @brief Very stiff systems will likely result in a step size that is not useable for the solver
     StepSizeTooSmall,
+    /// @brief Matrices that are singular more than once will set this value. At present, this should never be returned
     RepeatedlySingularMatrix,
+    /// @brief Mostly this value is returned by systems that tend toward chemical explosions
     NaNDetected
   };
 
@@ -149,29 +156,40 @@ namespace micm
    public:
     struct SolverStats
     {
-      uint64_t function_calls{};    // Nfun
-      uint64_t jacobian_updates{};  // Njac
-      uint64_t number_of_steps{};   // Nstp
-      uint64_t accepted{};          // Nacc
-      uint64_t rejected{};          // Nrej
-      uint64_t decompositions{};    // Ndec
-      uint64_t solves{};            // Nsol
-      uint64_t singular{};          // Nsng
-      uint64_t total_steps{};       // Ntotstp
-
+      /// @brief The number of forcing function calls
+      uint64_t function_calls{};
+      /// @brief The number of jacobian function calls
+      uint64_t jacobian_updates{};
+      /// @brief The total number of internal time steps taken
+      uint64_t number_of_steps{};
+      /// @brief The number of accepted integrations
+      uint64_t accepted{};
+      /// @brief The number of rejected integrations
+      uint64_t rejected{};
+      /// @brief The number of LU decompositions
+      uint64_t decompositions{};
+      /// @brief The number of linear solvers
+      uint64_t solves{};
+      /// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix is never singular
+      uint64_t singular{};
+      /// @brief The cumulative amount of time spent calculating the forcing function
       std::chrono::duration<double, std::nano> total_forcing_time{};
+      /// @brief The cumulative amount of time spent calculating the jacobian
       std::chrono::duration<double, std::nano> total_jacobian_time{};
+      /// @brief The cumulative amount of time spent calculating the linear factorization
       std::chrono::duration<double, std::nano> total_linear_factor_time{};
+      /// @brief The cumulative amount of time spent calculating the linear solve
       std::chrono::duration<double, std::nano> total_linear_solve_time{};
 
+      /// @brief Set all member variables to zero
       void Reset();
     };
 
     struct [[nodiscard]] SolverResult
     {
       /// @brief The new state computed by the solver
       MatrixPolicy<double> result_{};
-      /// @brief The finals state the solver was in
+      /// @brief The final state the solver was in
       SolverState state_ = SolverState::NotYetCalled;
       /// @brief A collection of runtime state for this call of the solver
       SolverStats stats_{};

include/micm/solver/rosenbrock.inl

@@ -4,9 +4,9 @@
 #define TIMED_METHOD(assigned_increment, time_it, method, ...) \
     { \
       if constexpr (time_it) { \
-          auto start = std::chrono::steady_clock::now(); \
+          auto start = std::chrono::high_resolution_clock::now(); \
           method(__VA_ARGS__); \
-          auto end = std::chrono::steady_clock::now(); \
+          auto end = std::chrono::high_resolution_clock::now(); \
           assigned_increment += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start); \
       } else { \
           method(__VA_ARGS__); \
@@ -389,7 +389,6 @@ namespace micm
     decompositions = 0;
     solves = 0;
     singular = 0;
-    total_steps = 0;
     total_forcing_time = std::chrono::nanoseconds::zero();
     total_jacobian_time = std::chrono::nanoseconds::zero();
     total_linear_factor_time = std::chrono::nanoseconds::zero();
@@ -608,10 +607,9 @@ namespace micm
                 parameters_.safety_factor_ / std::pow(error, 1 / parameters_.estimator_of_local_order_)));
         double Hnew = H * fac;
 
-        // Check the error magnitude and adjust step size
         stats_.number_of_steps += 1;
-        stats_.total_steps += 1;
 
+        // Check the error magnitude and adjust step size
         if (std::isnan(error))
         {
           Y.AsVector().assign(Ynew.AsVector().begin(), Ynew.AsVector().end());

test/tutorial/CMakeLists.txt

@@ -8,6 +8,7 @@ include(test_util)
 
 create_standard_test(NAME README_example SOURCES test_README_example.cpp)
 
+create_standard_test(NAME but_how_fast_is_it SOURCES test_but_how_fast_is_it.cpp)
 create_standard_test(NAME multiple_grid_cells SOURCES test_multiple_grid_cells.cpp)
 create_standard_test(NAME rate_constants_no_user_defined_example_by_hand SOURCES test_rate_constants_no_user_defined_by_hand.cpp)
 create_standard_test(NAME rate_constants_user_defined_example_by_hand SOURCES test_rate_constants_user_defined_by_hand.cpp)

test/tutorial/test_but_how_fast_is_it.cpp

@@ -0,0 +1,97 @@
+#include <iomanip>
+#include <iostream>
+#include <micm/process/user_defined_rate_constant.hpp>
+#include <micm/solver/rosenbrock.hpp>
+#include <chrono>
+
+// Use our namespace so that this example is easier to read
+using namespace micm;
+
+// The Rosenbrock solver can use many matrix ordering types
+// Here, we use the default ordering, but we still need to provide a templated
+// Arguent to the solver so it can use the proper ordering with any data type
+template<class T>
+using SparseMatrixPolicy = SparseMatrix<T>;
+
+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()
+                   .reactants({ a })
+                   .products({ yields(b, 1) })
+                   .rate_constant(UserDefinedRateConstant({ .label_ = "r1" }))
+                   .phase(gas_phase);
+
+  Process r2 = Process::create()
+                   .reactants({ b, b })
+                   .products({ yields(b, 1), yields(c, 1) })
+                   .rate_constant(UserDefinedRateConstant({ .label_ = "r2" }))
+                   .phase(gas_phase);
+
+  Process r3 = Process::create()
+                   .reactants({ b, c })
+                   .products({ yields(a, 1), yields(c, 1) })
+                   .rate_constant(UserDefinedRateConstant({ .label_ = "r3" }))
+                   .phase(gas_phase);
+
+  RosenbrockSolver<Matrix, SparseMatrixPolicy> solver{ System(SystemParameters{ .gas_phase_ = gas_phase }),
+                                                       std::vector<Process>{ r1, r2, r3 },
+                                                       RosenbrockSolverParameters::three_stage_rosenbrock_parameters(
+                                                           3, false) };
+
+  State<Matrix> 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 < solver.parameters_.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
+
+  auto result = solver.Solve(time_step, state);
+  std::cout << "Solver state: " << StateToString(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 << "singular: " << result.stats_.singular << std::endl;
+  std::cout << "final simulation time: " << result.final_time_ << std::endl;
+
+  auto start = std::chrono::high_resolution_clock::now();
+  result = solver.Solve<true>(time_step, state);
+  auto end = std::chrono::high_resolution_clock::now();
+  auto solve_time = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start);
+  std::cout << "Total solve time: " << solve_time.count() << " nanoseconds" << std::endl;
+  std::cout << "total_forcing_time: " << result.stats_.total_forcing_time.count() << " nanoseconds" << std::endl;
+  std::cout << "total_jacobian_time: " << result.stats_.total_jacobian_time.count() << " nanoseconds" << std::endl;
+  std::cout << "total_linear_factor_time: " << result.stats_.total_linear_factor_time.count() << " nanoseconds" << std::endl;
+  std::cout << "total_linear_solve_time: " << result.stats_.total_linear_solve_time.count() << " nanoseconds" << std::endl;
+}