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/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 DynamicSolver.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 DynamicSolver.hpp
     LINK_LIBRARIES
       PUBLIC
       SUNDIALS::nvecserial
@@ -42,3 +40,26 @@ else()
       PUBLIC
       GridKit::sparse_matrix)
 endif()
+
+if(GRIDKIT_ENABLE_RESOLVE)
+  gridkit_add_library(
+    rosenbrock
+    SOURCES Rosenbrock.cpp RosenbrockTableaus.cpp
+    HEADERS Rosenbrock.hpp
+            Native/ErrorNorm.hpp
+            Native/StepControl.hpp
+            Native/StepController.hpp
+    LINK_LIBRARIES
+      PRIVATE
+      ReSolve::ReSolve
+      PUBLIC
+      GridKit::dense_vector
+      PUBLIC
+      GridKit::definitions
+      PUBLIC
+      GridKit::utilities_logger
+      PUBLIC
+      GridKit::sparse_matrix)
+endif()
+
+add_subdirectory(Native)

GridKit/Solver/Dynamic/Native/AdaptiveStep.cpp

@@ -0,0 +1,33 @@
+#include "AdaptiveStep.hpp"
+
+#include <algorithm>
+#include <cmath>
+
+#include <GridKit/Constants.hpp>
+
+namespace AnalysisManager
+{
+  namespace NativeDynamicSolver
+  {
+    /**
+     * @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 NativeDynamicSolver
+} // namespace AnalysisManager

GridKit/Solver/Dynamic/Native/AdaptiveStep.hpp

@@ -0,0 +1,72 @@
+#pragma once
+
+#include <GridKit/Solver/Dynamic/Native/StepController.hpp>
+
+namespace AnalysisManager
+{
+  namespace NativeDynamicSolver
+  {
+
+    /**
+     * @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 NativeDynamicSolver
+} // namespace AnalysisManager

GridKit/Solver/Dynamic/Native/CMakeLists.txt

@@ -0,0 +1,21 @@
+if(GRIDKIT_ENABLE_RESOLVE)
+  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/Native/ErrorNorm.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include <GridKit/LinearAlgebra/Vector/Vector.hpp>
+#include <GridKit/LinearAlgebra/Vector/VectorHandler.hpp>
+#include <GridKit/ScalarTraits.hpp>
+
+namespace AnalysisManager
+{
+  namespace NativeDynamicSolver
+  {
+    /**
+     * @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 NativeDynamicSolver
+} // namespace AnalysisManager

GridKit/Solver/Dynamic/Native/FixedStep.cpp

@@ -0,0 +1,22 @@
+#include "FixedStep.hpp"
+
+namespace AnalysisManager
+{
+  namespace NativeDynamicSolver
+  {
+    /**
+     * @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 NativeDynamicSolver
+} // namespace AnalysisManager

GridKit/Solver/Dynamic/Native/FixedStep.hpp

@@ -0,0 +1,34 @@
+#pragma once
+
+#include <GridKit/Solver/Dynamic/Native/StepController.hpp>
+
+namespace AnalysisManager
+{
+  namespace NativeDynamicSolver
+  {
+    /**
+     * @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 NativeDynamicSolver
+} // namespace AnalysisManager

GridKit/Solver/Dynamic/Native/InfNorm.cpp

@@ -0,0 +1,55 @@
+#include "InfNorm.hpp"
+
+#include <format>
+
+#include <GridKit/Constants.hpp>
+
+namespace AnalysisManager
+{
+  namespace NativeDynamicSolver
+  {
+    /**
+     * @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 NativeDynamicSolver
+} // namespace AnalysisManager