Backware Euler with vectorizable matrix types (#596)

* starting to test all solver parameter types

* saving progress

* saving progress

* testing all stages analytically

* updating all interfaces

* correcting cuda build I hope

* testing jit against hires, e5, oregonator

* adding cuda solver builder test

* removing hires, e5, oregonator from cuda tests; they need their own kernels

* testing e5 from a configuration

* testing e5 jit integration

* testing e5 properly

* removing reset of L and U matrices (#594)

* oregonator from a configuration

* renaming things

* using different tolerances?

* moving state onto and off of host

* saving gpu changes

* updating cuda tests

* adding some better tolerances for cuda tests

* adding different tolerances for e5

* adding citation to e5

* thing

* formed hires equations

* using passing tolerances for cpu tests

* jit tolerances

* backward euler tests

* configuration for hires

* add AddToDiagonal function on sparse matrix

* use ForEach in Backward Euler

* add convergence check function to backward euler

* fix merge problems

* add vector matrix to analytical solver tests

* update JIT analytical tests

* set up general use analytical test function

* add general function for stiff analytical tests

* fix jit analytical tests

* update remaining analytical tests

* address review comments

* update cuda analytical tests

* update tolerances for cuda analytical tests

---------

Co-authored-by: Kyle Shores <kyle.shores44@gmail.com>
Co-authored-by: Kyle Shores <kshores@ucar.edu>

include/micm/configure/solver_config.hpp

@@ -514,7 +514,7 @@ namespace micm
  "Ea is specified when C is also specified for an Arrhenius reaction. Pick one." };
  }
  // Calculate 'C' using 'Ea'
- parameters.C_ = -1 * object["Ea"].get<double>() / BOLTZMANN_CONSTANT;
+ parameters.C_ = -1 * object["Ea"].get<double>() / constants::BOLTZMANN_CONSTANT;
  }
  arrhenius_rate_arr_.push_back(ArrheniusRateConstant(parameters));
  std::unique_ptr<ArrheniusRateConstant> rate_ptr = std::make_unique<ArrheniusRateConstant>(parameters);

include/micm/jit/solver/jit_linear_solver.inl

@@ -37,7 +37,8 @@ namespace micm
  if (matrix.NumberOfBlocks() != L || matrix.GroupVectorSize() != L)
  {
  std::string msg =
- "JIT functions require the number of grid cells solved together to match the vector dimension template parameter, "
+ "JIT functions require the number of grid cells solved together (" +
+ std::to_string(matrix.NumberOfBlocks()) + ") to match the vector dimension template parameter, "
  "currently: " +
  std::to_string(L);
  throw std::system_error(make_error_code(MicmJitErrc::InvalidMatrix), msg);

include/micm/jit/solver/jit_lu_decomposition.inl

@@ -31,7 +31,8 @@ namespace micm
  if (matrix.NumberOfBlocks() > L)
  {
  std::string msg =
- "JIT functions require the number of grid cells solved together to match the vector dimension template parameter, "
+ "JIT functions require the number of grid cells solved together (" +
+ std::to_string(matrix.NumberOfBlocks()) + ") to match the vector dimension template parameter, "
  "currently: " +
  std::to_string(L);
  throw std::system_error(make_error_code(MicmJitErrc::InvalidMatrix), msg);

include/micm/jit/solver/jit_rosenbrock.hpp

@@ -101,7 +101,8 @@ namespace micm
  if (jacobian.GroupVectorSize() != jacobian.NumberOfBlocks())
  {
  std::string msg =
- "JIT functions require the number of grid cells solved together to match the vector dimension template "
+ "JIT functions require the number of grid cells solved together (" +
+ std::to_string(jacobian.NumberOfBlocks()) + ") to match the vector dimension template "
  "parameter, currently: " +
  std::to_string(jacobian.GroupVectorSize());
  throw std::system_error(make_error_code(MicmJitErrc::InvalidMatrix), msg);

include/micm/process/branched_rate_constant.hpp

@@ -79,8 +79,8 @@ namespace micm
 
  inline BranchedRateConstant::BranchedRateConstant(const BranchedRateConstantParameters& parameters)
  : parameters_(parameters),
- k0_(2.0e-22 * AVOGADRO_CONSTANT * 1.0e-6 * std::exp(parameters_.n_)),
- z_(A(293.0, 2.45e19 / AVOGADRO_CONSTANT * 1.0e6) * (1.0 - parameters_.a0_) / parameters_.a0_)
+ k0_(2.0e-22 * constants::AVOGADRO_CONSTANT * 1.0e-6 * std::exp(parameters_.n_)),
+ z_(A(293.0, 2.45e19 / constants::AVOGADRO_CONSTANT * 1.0e6) * (1.0 - parameters_.a0_) / parameters_.a0_)
  {
  }
 

include/micm/process/surface_rate_constant.hpp

@@ -63,8 +63,8 @@ namespace micm
 
  inline SurfaceRateConstant::SurfaceRateConstant(const SurfaceRateConstantParameters& parameters)
  : parameters_(parameters),
- diffusion_coefficient_(parameters.species_.GetProperty<double>(GAS_DIFFUSION_COEFFICIENT)),
- mean_free_speed_factor_(8.0 * GAS_CONSTANT / (M_PI * parameters.species_.GetProperty<double>(MOLECULAR_WEIGHT)))
+ diffusion_coefficient_(parameters.species_.GetProperty<double>(property_keys::GAS_DIFFUSION_COEFFICIENT)),
+ mean_free_speed_factor_(8.0 * constants::GAS_CONSTANT / (M_PI * parameters.species_.GetProperty<double>(property_keys::MOLECULAR_WEIGHT)))
  {
  }
 

include/micm/solver/backward_euler.hpp

@@ -9,6 +9,7 @@
 #include <micm/solver/state.hpp>
 #include <micm/system/system.hpp>
 #include <micm/util/jacobian.hpp>
+#include <micm/util/matrix.hpp>
 
 #include <algorithm>
 #include <array>
@@ -66,7 +67,21 @@ namespace micm
  /// @param time_step Time [s] to advance the state by
  /// @param state The state to advance
  /// @return result of the solver (success or failure, and statistics)
- SolverResult Solve(double time_step, auto& state);
+ SolverResult Solve(double time_step, auto& state) const;
+
+ /// @brief Determines whether the residual is small enough to stop the
+ /// internal solver iteration
+ /// @param parameters Solver parameters
+ /// @param residual The residual to check
+ /// @param state The current state being solved for
+ /// @return true if the residual is small enough to stop the iteration
+ template<class DenseMatrixPolicy>
+ static bool IsConverged(const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, const DenseMatrixPolicy& state)
+ requires(!VectorizableDense<DenseMatrixPolicy>);
+ template<class DenseMatrixPolicy>
+ static bool IsConverged(const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, const DenseMatrixPolicy& state)
+ requires(VectorizableDense<DenseMatrixPolicy>);
+
  };
 
 } // namespace micm

include/micm/solver/backward_euler.inl

@@ -43,25 +43,21 @@ inline std::error_code make_error_code(MicmBackwardEulerErrc e)
 namespace micm
 {
  template<class RatesPolicy, class LinearSolverPolicy>
- inline SolverResult BackwardEuler<RatesPolicy, LinearSolverPolicy>::Solve(double time_step, auto& state)
+ inline SolverResult BackwardEuler<RatesPolicy, LinearSolverPolicy>::Solve(double time_step, auto& state) const
  {
  // A fully implicit euler implementation is given by the following equation:
  // y_{n+1} = y_n + H * f(t_{n+1}, y_{n+1})
  // This is a root finding problem because you need to know y_{n+1} to compute f(t_{n+1}, y_{n+1})
  // you need to solve the equation y_{n+1} - y_n - H f(t_{n+1}, y_{n+1}) = 0
- // We will also use the same logic used by cam-chem to determine the time step
- // That scheme is this:
- // Start with H = time_step
- // if that fails, try H = H/2 several times
- // if that fails, try H = H/10 once
- // if that fails, accept the current H but do not update the Yn vector
- // the number of time step reduction is controlled by the time_step_reductions parameter
+ // A series of time step reductions are used after failed solves to try to find a solution
+ // These reductions are controlled by the time_step_reductions parameter in the solver parameters
+ // if the last attempt to reduce the timestep fails,
+ // accept the current H but do not update the Yn vector
 
  SolverResult result;
 
- double small = parameters_.small;
  std::size_t max_iter = parameters_.max_number_of_steps_;
- const auto time_step_reductions = parameters_.time_step_reductions;
+ const auto time_step_reductions = parameters_.time_step_reductions_;
 
  double H = time_step;
  double t = 0.0;
@@ -89,46 +85,30 @@ namespace micm
  // after the first iteration Yn1 is updated to the new solution
  // so we can use Yn1 to calculate the forcing and jacobian
  // calculate forcing
- std::fill(forcing.AsVector().begin(), forcing.AsVector().end(), 0.0);
+ forcing.Fill(0.0);
  rates_.AddForcingTerms(state.rate_constants_, Yn1, forcing);
  result.stats_.function_calls_++;
 
- // calculate jacobian
- std::fill(state.jacobian_.AsVector().begin(), state.jacobian_.AsVector().end(), 0.0);
+ // calculate the negative jacobian
+ state.jacobian_.Fill(0.0);
  rates_.SubtractJacobianTerms(state.rate_constants_, Yn1, state.jacobian_);
  result.stats_.jacobian_updates_++;
 
- // subtract the inverse of the time step from the diagonal
- // TODO: handle vectorized jacobian matrix
- for (auto& jac : state.jacobian_.AsVector())
- {
- jac *= -1;
- }
- for (std::size_t i_block = 0; i_block < state.jacobian_.NumberOfBlocks(); ++i_block)
- {
- auto jacobian_vector = std::next(state.jacobian_.AsVector().begin(), i_block * state.jacobian_.FlatBlockSize());
- for (const auto& i_elem : jacobian_diagonal_elements_)
- jacobian_vector[i_elem] -= 1 / H;
- }
-
+ // add the inverse of the time step from the diagonal
+ state.jacobian_.AddToDiagonal(1 / H);
+
  // We want to solve this equation for a zero
  // (y_{n+1} - y_n) / H = f(t_{n+1}, y_{n+1})
 
  // try to find the root by factoring and solving the linear system
  linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular);
  result.stats_.decompositions_++;
 
- auto yn1_iter = Yn1.begin();
- auto yn_iter = Yn.begin();
- auto forcing_iter = forcing.begin();
  // forcing_blk in camchem
- // residual = (Yn1 - Yn) / H - forcing;
+ // residual = forcing - (Yn1 - Yn) / H
  // since forcing is only used once, we can reuse it to store the residual
- for (; yn1_iter != Yn1.end(); ++yn1_iter, ++yn_iter, ++forcing_iter)
- {
- *forcing_iter = (*yn1_iter - *yn_iter) / H - *forcing_iter;
- }
-
+ forcing.ForEach([&](double& f, const double& yn1, const double& yn) { f -= (yn1 - yn) / H; }, Yn1, Yn);
+
  // the result of the linear solver will be stored in forcing
  // this represents the change in the solution
  linear_solver_.Solve(forcing, state.lower_matrix_, state.upper_matrix_);
@@ -137,31 +117,14 @@ namespace micm
  // solution_blk in camchem
  // Yn1 = Yn1 + residual;
  // always make sure the solution is positive regardless of which iteration we are on
- forcing_iter = forcing.begin();
- yn1_iter = Yn1.begin();
- for (; forcing_iter != forcing.end(); ++forcing_iter, ++yn1_iter)
- {
- *yn1_iter = std::max(0.0, *yn1_iter + *forcing_iter);
- }
+ Yn1.ForEach([&](double& yn1, const double& f) { yn1 = std::max( 0.0, yn1 + f ); }, forcing);
 
  // if this is the first iteration, we don't need to check for convergence
  if (iterations++ == 0)
  continue;
 
  // check for convergence
- forcing_iter = forcing.begin();
- yn1_iter = Yn1.begin();
-
- // convergence happens when the absolute value of the change to the solution
- // is less than a tolerance times the absolute value of the solution
- auto abs_tol_iter = parameters_.absolute_tolerance_.begin();
- do
- {
- // changes that are much smaller than the tolerance are negligible and we assume can be accepted
- converged = (std::abs(*forcing_iter) <= small) || (std::abs(*forcing_iter) <= *abs_tol_iter) ||
- (std::abs(*forcing_iter) <= parameters_.relative_tolerance_ * std::abs(*yn1_iter));
- ++forcing_iter, ++yn1_iter, ++abs_tol_iter;
- } while (converged && forcing_iter != forcing.end());
+ converged = IsConverged(parameters_, forcing, Yn1);
  } while (!converged && iterations < max_iter);
 
  if (!converged)
@@ -203,4 +166,74 @@ namespace micm
  result.final_time_ = t;
  return result;
  }
+
+ template<class RatesPolicy, class LinearSolverPolicy>
+ template<class DenseMatrixPolicy>
+ inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, const DenseMatrixPolicy& state)
+ requires(!VectorizableDense<DenseMatrixPolicy>)
+ {
+ double small = parameters.small_;
+ double rel_tol = parameters.relative_tolerance_;
+ auto& abs_tol = parameters.absolute_tolerance_;
+ auto residual_iter = residual.AsVector().begin();
+ auto state_iter = state.AsVector().begin();
+ const std::size_t n_elem = residual.NumRows() * residual.NumColumns();
+ const std::size_t n_vars = abs_tol.size();
+ for (std::size_t i = 0; i < n_elem; ++i)
+ {
+ if (std::abs(*residual_iter) > small && std::abs(*residual_iter) > abs_tol[i % n_vars] &&
+ std::abs(*residual_iter) > rel_tol * std::abs(*state_iter))
+ {
+ return false;
+ }
+ ++residual_iter, ++state_iter;
+ }
+ return true;
+ }
+
+ template<class RatesPolicy, class LinearSolverPolicy>
+ template<class DenseMatrixPolicy>
+ inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, const DenseMatrixPolicy& state)
+ requires(VectorizableDense<DenseMatrixPolicy>)
+ {
+ double small = parameters.small_;
+ double rel_tol = parameters.relative_tolerance_;
+ auto& abs_tol = parameters.absolute_tolerance_;
+ auto residual_iter = residual.AsVector().begin();
+ auto state_iter = state.AsVector().begin();
+ const std::size_t n_elem = residual.NumRows() * residual.NumColumns();
+ const std::size_t L = residual.GroupVectorSize();
+ const std::size_t n_vars = abs_tol.size();
+ const std::size_t whole_blocks = std::floor(residual.NumRows() / L) * residual.GroupSize();
+
+ // evaluate the rows that fit exactly into the vectorizable dimension (L)
+ for (std::size_t i = 0; i < whole_blocks; ++i)
+ {
+ if (std::abs(*residual_iter) > small && std::abs(*residual_iter) > abs_tol[(i / L) % n_vars] &&
+ std::abs(*residual_iter) > rel_tol * std::abs(*state_iter))
+ {
+ return false;
+ }
+ ++residual_iter, ++state_iter;
+ }
+
+ // evaluate the remaining rows
+ const std::size_t remaining_rows = residual.NumRows() % L;
+ if (remaining_rows > 0)
+ {
+ for (std::size_t y = 0; y < residual.NumColumns(); ++y)
+ {
+ const std::size_t offset = y * L;
+ for (std::size_t i = offset; i < offset + remaining_rows; ++i)
+ {
+ if (std::abs(residual_iter[i]) > small && std::abs(residual_iter[i]) > abs_tol[y] &&
+ std::abs(residual_iter[i]) > rel_tol * std::abs(state_iter[i]))
+ {
+ return false;
+ }
+ }
+ }
+ }
+ return true;
+ }
 } // namespace micm

include/micm/solver/backward_euler_solver_parameters.hpp

@@ -21,9 +21,11 @@ namespace micm
 
  std::vector<double> absolute_tolerance_;
  double relative_tolerance_{ 1.0e-8 };
- double small{ 1.0e-40 };
+ double small_{ 1.0e-40 };
  size_t max_number_of_steps_{ 11 };
- std::array<double, 5> time_step_reductions{ 0.5, 0.5, 0.5, 0.5, 0.1 };
+ // The time step reductions are used to determine the time step after a failed solve
+ // This default set of time step reductions is used by CAM-Chem
+ std::array<double, 5> time_step_reductions_{ 0.5, 0.5, 0.5, 0.5, 0.1 };
  };
 
 } // namespace micm

include/micm/solver/linear_solver.hpp

@@ -77,7 +77,7 @@ namespace micm
  virtual ~LinearSolver() = default;
 
  /// @brief Decompose the matrix into upper and lower triangular matrices
- void Factor(const SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix);
+ void Factor(const SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix) const;
 
  /// @brief Decompose the matrix into upper and lower triangular matrices
  /// @param matrix Matrix to decompose into lower and upper triangular matrices
@@ -86,7 +86,7 @@ namespace micm
  const SparseMatrixPolicy& matrix,
  SparseMatrixPolicy& lower_matrix,
  SparseMatrixPolicy& upper_matrix,
- bool& is_singular);
+ bool& is_singular) const;
 
  /// @brief Solve for x in Ax = b. x should be a copy of b and after Solve finishes x will contain the result
  template<class MatrixPolicy>

include/micm/solver/linear_solver.inl

@@ -111,7 +111,7 @@ namespace micm
  inline void LinearSolver<SparseMatrixPolicy, LuDecompositionPolicy>::Factor(
  const SparseMatrixPolicy& matrix,
  SparseMatrixPolicy& lower_matrix,
- SparseMatrixPolicy& upper_matrix)
+ SparseMatrixPolicy& upper_matrix) const
  {
  MICM_PROFILE_FUNCTION();
 
@@ -123,7 +123,7 @@ namespace micm
  const SparseMatrixPolicy& matrix,
  SparseMatrixPolicy& lower_matrix,
  SparseMatrixPolicy& upper_matrix,
- bool& is_singular)
+ bool& is_singular) const
  {
  MICM_PROFILE_FUNCTION();
 

include/micm/util/constants.hpp

@@ -2,6 +2,12 @@
 // SPDX-License-Identifier: Apache-2.0
 #pragma once
 
-static constexpr double BOLTZMANN_CONSTANT = 1.380649e-23; // J K^{-1}
-static constexpr double AVOGADRO_CONSTANT = 6.02214076e23; // # mol^{-1}
-static constexpr double GAS_CONSTANT = BOLTZMANN_CONSTANT * AVOGADRO_CONSTANT; // J K^{-1} mol^{-1}
+namespace micm
+{
+ namespace constants
+ {
+ static constexpr double BOLTZMANN_CONSTANT = 1.380649e-23; // J K^{-1}
+ static constexpr double AVOGADRO_CONSTANT = 6.02214076e23; // # mol^{-1}
+ static constexpr double GAS_CONSTANT = BOLTZMANN_CONSTANT * AVOGADRO_CONSTANT; // J K^{-1} mol^{-1}
+ } // namespace constants
+} // namespace micm

include/micm/util/property_keys.hpp

@@ -4,5 +4,11 @@
 
 #include <string>
 
-static const std::string GAS_DIFFUSION_COEFFICIENT = "diffusion coefficient [m2 s-1]";
-static const std::string MOLECULAR_WEIGHT = "molecular weight [kg mol-1]";
+namespace micm
+{
+ namespace property_keys
+ {
+ static const std::string GAS_DIFFUSION_COEFFICIENT = "diffusion coefficient [m2 s-1]";
+ static const std::string MOLECULAR_WEIGHT = "molecular weight [kg mol-1]";
+ }
+}