Backware Euler with vectorizable matrix types

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

@@ -75,7 +75,7 @@ namespace micm
         const std::function<LuDecompositionPolicy(const SparseMatrixPolicy&)> create_lu_decomp);
 
     /// @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
@@ -84,7 +84,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]";
+  }
+}