Add Rosenbrock integrators

CHANGELOG.md

@@ -69,6 +69,8 @@
 - Removed `COO_Matrix` class.
 - Added portable `Vector` class and policy-based memory utilities.
 - Add new `LinearSolver` interface for linear solvers.
+- Added new `Rosenbrock` integrator.
+- Clarified naming conventions for macros.
 
 ## v0.1
 

CONTRIBUTING.md

@@ -285,6 +285,10 @@ name. Use all caps (screaming snake case).
    constexpr double EXP = 2.7183      // Yes
 ```
 
+### Macros
+
+Macros should use all caps, same as constants.
+
 ### Enums (enumerated types)
 
 Always define `enum`s inside `GridKit` namespace. The `enum` name should

GridKit/LinearAlgebra/CMakeLists.txt

@@ -1,4 +1,4 @@
 add_subdirectory(SparseMatrix)
 add_subdirectory(DenseMatrix)
 add_subdirectory(Vector)
-add_subdirectory(Solver)
+add_subdirectory(Solver)

GridKit/LinearAlgebra/Solver/CMakeLists.txt

@@ -3,8 +3,5 @@ if(GRIDKIT_ENABLE_RESOLVE)
     resolve_system_solver
     SOURCES ResolveSystemSolver.cpp
     HEADERS ResolveSystemSolver.hpp
-    LINK_LIBRARIES
-      PRIVATE
-      ReSolve::ReSolve
-    )
-endif()
+    LINK_LIBRARIES PRIVATE ReSolve::ReSolve)
+endif()

GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp

@@ -25,7 +25,8 @@ namespace GridKit
     }
 
     template <class ScalarT, typename IdxT>
-    ResolveSystemSolver<ScalarT, IdxT>::ResolveSystemSolver(ReSolve::SystemSolver& lin_solver, GridKit::memory::MemorySpace memspace) : lin_solver_(lin_solver), memspace_(memorySpaceAsResolve(memspace))
+    ResolveSystemSolver<ScalarT, IdxT>::ResolveSystemSolver(ReSolve::SystemSolver& lin_solver, GridKit::memory::MemorySpace memspace)
+      : lin_solver_(lin_solver), memspace_(memorySpaceAsResolve(memspace))
     {
     }
 
@@ -101,4 +102,4 @@ namespace GridKit
 
     template class ResolveSystemSolver<double, int>;
   } // namespace LinearAlgebra
-} // namespace GridKit
+} // namespace GridKit

GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp

@@ -42,4 +42,4 @@ namespace GridKit
       ReSolve::memory::MemorySpace          memspace_;
     };
   } // namespace LinearAlgebra
-} // namespace GridKit
+} // namespace GridKit

GridKit/MemoryUtilities/ResolveMemoryUtils.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include <GridKit/MemoryUtilities/MemoryUtils.hpp>
+
+#include <resolve/MemoryUtils.hpp>
+
+namespace GridKit
+{
+  namespace memory
+  {
+    /**
+     * @brief Converts a GridKit \ref MemorySpace to its corresponding Re::Solve MemorySpace
+     *
+     */
+    inline ReSolve::memory::MemorySpace memorySpaceAsResolve(MemorySpace memspace)
+    {
+      switch (memspace)
+      {
+      case HOST:
+        return ReSolve::memory::HOST;
+      case DEVICE:
+        return ReSolve::memory::DEVICE;
+      default:
+        throw "Memory space not supported";
+      }
+    }
+  } // namespace memory
+} // namespace GridKit

GridKit/Solver/Dynamic/AdaptiveStep.cpp

@@ -0,0 +1,29 @@
+#include "AdaptiveStep.hpp"
+
+#include <cmath>
+
+#include <GridKit/Constants.hpp>
+
+namespace Integrator
+{
+  /**
+   * @brief Standard textbook adaptive controller. Accept if `err <= 1` and use
+   *
+   * \f[h_{new} = h * \min \left\{fac_{max}, \max\left\{fac_{min}, fac_{scale} \cdot e ^{-1/p}\right\}\right\}.\f]
+   *
+   */
+  template <typename RealT>
+  StepControl<RealT> AdaptiveStep<RealT>::nextStep(RealT err, StepControl<RealT> prev_step, uint8_t method_order)
+  {
+    StepControl<RealT> next_step = prev_step;
+
+    double h_mult = std::min(params_.fac_max_, std::max(params_.fac_scale_ * std::pow(err, GridKit::MINUS_ONE<RealT> / method_order), params_.fac_min_));
+
+    next_step.accept_     = err <= 1;
+    next_step.step_size_ *= h_mult;
+
+    return next_step;
+  }
+
+  template class AdaptiveStep<double>;
+} // namespace Integrator

GridKit/Solver/Dynamic/AdaptiveStep.hpp

@@ -0,0 +1,69 @@
+#pragma once
+
+#include <GridKit/Solver/Dynamic/StepController.hpp>
+
+namespace Integrator
+{
+
+  /**
+   * @brief A simple textbook adaptive `StepController` which seeks to meet a relative and absolute tolerance
+   *        based on an error estimate.
+   *
+   */
+  template <typename RealT>
+  class AdaptiveStep : public StepController<RealT>
+  {
+    /**
+     * @brief Parameters for the step controller.
+     *
+     */
+    struct Parameters
+    {
+      /**
+       * @brief The minimum multiple by which the step size can be multiplied to obtain the new step size.
+       *        Increasing this can allow the integrator to be slightly more conservative in selecting the step size
+       *        - decreasing the number of steps taken but increasing the risk of failing the next step.
+       *
+       * @note Should be between 0 and 1.
+       *
+       */
+      RealT fac_min_   = 0.2;
+      /**
+       * @brief The maximum multiple by which the step size can be multiplied to obtain the new step size.
+       *        Decreasing this will make the integrator more conservative in selecting the step size -
+       *        increasing the number of steps taken but decreasing the risk of failing the next step.
+       *
+       * @note Should be greater than 1.
+       *
+       */
+      RealT fac_max_   = 5.0;
+      /**
+       * @brief A "fudge factor" introduced to decrease risk of failing a step. The larger the fudge factor,
+       *        the more likely steps will fail, but fewer steps will be taken.
+       *
+       * @note Should be between 0 and 1.
+       *
+       */
+      RealT fac_scale_ = 0.9;
+    } params_;
+
+  public:
+    AdaptiveStep(const Parameters& params)
+      : params_(params)
+    {
+    }
+
+    StepControl<RealT> nextStep(RealT err, StepControl<RealT> prev_step, uint8_t method_order) final;
+
+    /**
+     * @brief This controller uses error estimates.
+     *
+     * @see `nextStep()`
+     *
+     */
+    constexpr bool usesError() const final
+    {
+      return true;
+    }
+  };
+} // namespace Integrator

GridKit/Solver/Dynamic/CMakeLists.txt

@@ -5,13 +5,11 @@
 #    - Cameron Rutherford <cameron.rutherford@pnnl.gov>
 #]]
 
-set(_install_headers DynamicSolver.hpp Ida.hpp)
-
 if(GRIDKIT_ENABLE_SUNDIALS_SPARSE)
   gridkit_add_library(
     solvers_dyn
     SOURCES Ida.cpp
-    HEADERS ${_install_headers}
+    HEADERS Ida.hpp
     LINK_LIBRARIES
       PUBLIC
       SUNDIALS::nvecserial
@@ -29,7 +27,7 @@ else()
   gridkit_add_library(
     solvers_dyn
     SOURCES Ida.cpp
-    HEADERS ${_install_headers}
+    HEADERS Ida.hpp
     LINK_LIBRARIES
       PUBLIC
       SUNDIALS::nvecserial
@@ -42,3 +40,41 @@ else()
       PUBLIC
       GridKit::sparse_matrix)
 endif()
+
+if(GRIDKIT_ENABLE_RESOLVE)
+  gridkit_add_library(
+    rosenbrock
+    SOURCES Rosenbrock.cpp RosenbrockTableaus.cpp
+    HEADERS Rosenbrock.hpp
+    LINK_LIBRARIES
+      PRIVATE
+      ReSolve::ReSolve
+      PUBLIC
+      GridKit::dense_vector
+      PUBLIC
+      GridKit::definitions
+      PUBLIC
+      GridKit::utilities_logger
+      PUBLIC
+      GridKit::sparse_matrix)
+
+  gridkit_add_library(
+    inf_norm
+    SOURCES InfNorm.cpp
+    HEADERS InfNorm.hpp
+    LINK_LIBRARIES
+      PRIVATE
+      ReSolve::ReSolve
+      PUBLIC
+      GridKit::dense_vector)
+endif()
+
+gridkit_add_library(
+  adaptive_step
+  SOURCES AdaptiveStep.cpp
+  HEADERS AdaptiveStep.hpp)
+
+gridkit_add_library(
+  fixed_step
+  SOURCES FixedStep.cpp
+  HEADERS FixedStep.hpp)

GridKit/Solver/Dynamic/ErrorNorm.hpp

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <GridKit/LinearAlgebra/Vector/Vector.hpp>
+#include <GridKit/LinearAlgebra/Vector/VectorHandler.hpp>
+#include <GridKit/ScalarTraits.hpp>
+
+namespace Integrator
+{
+  /**
+   * @brief Interface for error norms. Used to calculate the `err` parameter in `StepController::nextStep` based on a residual state error vector.
+   *
+   */
+  template <class ScalarT, typename IdxT>
+  class ErrorNorm
+  {
+    using State = GridKit::LinearAlgebra::Vector<ScalarT, IdxT>;
+    using RealT = typename GridKit::ScalarTraits<ScalarT>::RealT;
+
+  public:
+    /**
+     * @brief Calculate an error to be used by a step controller. Typically, an error > 1 indicates an error which does not meet tolerances, while
+     *        an error < 1 indicates an error which meets tolerances. For that reason, tolerances should be included in the calculation of the error.
+     *
+     * @param err The state error residual being measured.
+     * @param y The state that the error was calculated from. Can be used for proper relative error normalization.
+     * @param yprev The state from the previous step. Can be used for proper relative error normalization.
+     * @param handler A vector handler which can be used to facilitate vector operations.
+     * @param memspace The memory space which vector operations should be performed in/.
+     * @return The error.
+     *
+     * @todo Allow this method to fail, since it will likely involve linear algebra calls.
+     */
+    virtual RealT errorNorm(State& err, State& y, State& yprev, GridKit::LinearAlgebra::VectorHandler<ScalarT, IdxT>& handler, GridKit::memory::MemorySpace memspace) const = 0;
+  };
+} // namespace Integrator

GridKit/Solver/Dynamic/FixedStep.cpp

@@ -0,0 +1,19 @@
+#include "FixedStep.hpp"
+
+namespace Integrator
+{
+  /**
+   * @brief Fixed step - accept every step, no matter the error, and keep the step size the same.
+   *
+   */
+  template <typename RealT>
+  StepControl<RealT> FixedStep<RealT>::nextStep([[maybe_unused]] RealT err, StepControl<RealT> prev_step, [[maybe_unused]] uint8_t method_order)
+  {
+    return StepControl<RealT>{
+        .accept_    = true,
+        .step_size_ = prev_step.step_size_,
+    };
+  }
+
+  template class FixedStep<double>;
+} // namespace Integrator

GridKit/Solver/Dynamic/FixedStep.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <GridKit/Solver/Dynamic/StepController.hpp>
+
+namespace Integrator
+{
+  /**
+   * @brief A fixed step controller which doesn't change the step size and accepts every step.
+   *        Useful if you know what time scale your simulation operates on apriori and you're
+   *        using a method without an embedded error controller.
+   *
+   *        To set the fixed size, set the `Rosenbrock::Parameters::starting_step` parameter.
+   *
+   */
+  template <typename RealT>
+  class FixedStep : public StepController<RealT>
+  {
+    StepControl<RealT> nextStep(RealT err, StepControl<RealT> prev_step, uint8_t method_order) final;
+
+    /**
+     * @brief This controller does not use error estimates.
+     *
+     * @see `nextStep()`
+     *
+     */
+    constexpr bool usesError() const final
+    {
+      return false;
+    }
+  };
+} // namespace Integrator

GridKit/Solver/Dynamic/InfNorm.cpp

@@ -0,0 +1,52 @@
+#include "InfNorm.hpp"
+
+#include <format>
+
+#include <GridKit/Constants.hpp>
+
+namespace Integrator
+{
+  /**
+   * @brief Calculate the infinity error norm as
+   *
+   * \f[e = \max_i\frac{|\hat{e}_i|}{Atol_i + Rtol \cdot \max\{|y_{0i}|, |y_{1i}|\}},\f]
+   *
+   * where \f(y_0\f) is the initial state, \f(y_1\f) is the next state, and \f(\hat{e}\f) is the estimated error made in calculating
+   * the next state (typically \f(\hat{e} = y_1 - \hat{y}_1\f) for some different-order approximation \f(\hat{y}_1\f)).
+   *
+   * @param err \f(\hat{e}\f) in the above formula.
+   * @param y \f(y_1\f) in the above formula.
+   * @param yprev \f(y_0\f) in the above formula.
+   * @param handler The handler to be used for performing linear algebra operations.
+   * @param memspace The memory space to be used for performing linear lagebra operations.
+   * @see `Rosenbrock::errorEstimate()`
+   */
+  template <class ScalarT, typename IdxT>
+  InfNorm<ScalarT, IdxT>::RealT InfNorm<ScalarT, IdxT>::errorNorm(State& err, State& y, State& yprev, GridKit::LinearAlgebra::VectorHandler<ScalarT, IdxT>& handler, GridKit::memory::MemorySpace memspace) const
+  {
+    if (int err_code = workspace_.out_->copyFromExternal(&err, memspace, memspace))
+    {
+      throw std::format("GridKit::LinearAlgebra::Vector::copyFromExternal failed with error code {}", err_code);
+    }
+    if (int err_code = workspace_.scale_->copyFromExternal(&y, memspace, memspace))
+    {
+      throw std::format("GridKit::LinearAlgebra::Vector::copyFromExternal failed with error code {}", err_code);
+    }
+    if (int err_code = workspace_.yprev_abs_->copyFromExternal(&yprev, memspace, memspace))
+    {
+      throw std::format("GridKit::LinearAlgebra::Vector::copyFromExternal failed with error code {}", err_code);
+    }
+
+    handler.abs(workspace_.scale_.get(), workspace_.scale_.get(), memspace);
+    handler.abs(workspace_.yprev_abs_.get(), workspace_.scale_.get(), memspace);
+    handler.max(workspace_.yprev_abs_.get(), workspace_.scale_.get(), workspace_.scale_.get(), memspace);
+
+    // TODO: This scal shouldn't be necessary, but axpy doesn't support scaling the y parameter. In the future,
+    // the scaling should be able to be put on the next axpy.
+    handler.scal(params_.rel_tol_, workspace_.scale_.get(), memspace);
+    handler.axpy(GridKit::ONE<RealT>, params_.abs_tol_.get(), workspace_.scale_.get(), memspace);
+    handler.diagSolve(workspace_.scale_.get(), workspace_.out_.get(), memspace);
+
+    return handler.amax(workspace_.out_.get(), memspace);
+  }
+} // namespace Integrator