From 75621bbbdeb592866b3daa9ac30baaffe416450f Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 13 Jul 2026 15:53:54 -0400 Subject: [PATCH 1/6] Add LinearSolver interface --- CMakeLists.txt | 18 ++++ GridKit/LinearAlgebra/CMakeLists.txt | 1 + GridKit/LinearAlgebra/Solver/CMakeLists.txt | 10 +++ GridKit/LinearAlgebra/Solver/LinearSolver.hpp | 24 +++++ .../Solver/ResolveSystemSolver.cpp | 87 +++++++++++++++++++ .../Solver/ResolveSystemSolver.hpp | 33 +++++++ README.md | 1 + 7 files changed, 174 insertions(+) create mode 100644 GridKit/LinearAlgebra/Solver/CMakeLists.txt create mode 100644 GridKit/LinearAlgebra/Solver/LinearSolver.hpp create mode 100644 GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp create mode 100644 GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 729630e44..9240a41a4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,6 +33,9 @@ option(GridKit_ENABLE_SUNDIALS "Enable SUNDIALS support" OFF) # Enzyme support is disabled by default option(GridKit_ENABLE_ENZYME "Enable automatic differentiation with Enzyme" OFF) +# ReSolve support is disabled by default +option(GridKit_ENABLE_RESOLVE "Enable ReSolve support" OFF) + # Enable the address sanitizer option(GridKit_ENABLE_ASAN "Enable the address sanitizer" OFF) @@ -52,6 +55,7 @@ list( IPOPT SUNDIALS ENZYME + RESOLVE ASAN UBSAN OPENMP) @@ -147,6 +151,20 @@ if(GRIDKIT_ENABLE_THREADS) find_package(Threads REQUIRED) endif() +if(GRIDKIT_ENABLE_RESOLVE) + find_package( + ReSolve + CONFIG + PATHS + ${ReSolve_DIR} + ${ReSolve_DIR}/share/resolve/cmake + ENV + LD_LIBRARY_PATH + ENV + DYLD_LIBRARY_PATH + REQUIRED) +endif() + # Macro that adds libraries include(GridkitAddLibrary) diff --git a/GridKit/LinearAlgebra/CMakeLists.txt b/GridKit/LinearAlgebra/CMakeLists.txt index e645e23ed..0fe740023 100644 --- a/GridKit/LinearAlgebra/CMakeLists.txt +++ b/GridKit/LinearAlgebra/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory(SparseMatrix) add_subdirectory(DenseMatrix) add_subdirectory(Vector) +add_subdirectory(Solver) \ No newline at end of file diff --git a/GridKit/LinearAlgebra/Solver/CMakeLists.txt b/GridKit/LinearAlgebra/Solver/CMakeLists.txt new file mode 100644 index 000000000..356cc12e4 --- /dev/null +++ b/GridKit/LinearAlgebra/Solver/CMakeLists.txt @@ -0,0 +1,10 @@ +if(GRIDKIT_ENABLE_RESOLVE) + gridkit_add_library( + resolve_system_solver + SOURCES ResolveSystemSolver.cpp + HEADERS ResolveSystemSolver.hpp + LINK_LIBRARIES + PRIVATE + ReSolve::ReSolve + ) +endif() \ No newline at end of file diff --git a/GridKit/LinearAlgebra/Solver/LinearSolver.hpp b/GridKit/LinearAlgebra/Solver/LinearSolver.hpp new file mode 100644 index 000000000..1222f2543 --- /dev/null +++ b/GridKit/LinearAlgebra/Solver/LinearSolver.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include +#include + +#include "GridKit/LinearAlgebra/Vector/Vector.hpp" + +namespace GridKit +{ + namespace LinearAlgebra + { + + template + class LinearSolver + { + public: + using RealT = GridKit::ScalarTraits::RealT; + + virtual int configureSolver(GridKit::LinearAlgebra::CsrMatrix& matrix) = 0; + virtual int setupSolver(bool reuse_factors = false) = 0; + virtual int solve(GridKit::LinearAlgebra::Vector& rhs, GridKit::LinearAlgebra::Vector& lhs) = 0; + }; + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp new file mode 100644 index 000000000..eba405a48 --- /dev/null +++ b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp @@ -0,0 +1,87 @@ +#include "ResolveSystemSolver.hpp" + +#include "GridKit/MemoryUtilities/MemoryUtils.hpp" +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + ReSolve::memory::MemorySpace memorySpaceAsResolve(memory::MemorySpace memspace) + { + switch (memspace) + { + case memory::HOST: + return ReSolve::memory::HOST; + case memory::DEVICE: + return ReSolve::memory::DEVICE; + default: + throw "Memory space not supported"; + } + } + + template + ResolveSystemSolver::ResolveSystemSolver(ReSolve::SystemSolver& lin_solver, GridKit::memory::MemorySpace memspace) : lin_solver_(lin_solver), memspace_(memorySpaceAsResolve(memspace)) + { + } + + template + int ResolveSystemSolver::configureSolver(GridKit::LinearAlgebra::CsrMatrix& matrix) + { + matrix_ = std::make_unique(matrix.getNumRows(), matrix.getNumColumns(), matrix.getNnz()); + + if (int err = matrix_->setDataPointers(matrix.getRowData(), matrix.getColData(), matrix.getValues(), memspace_)) + return err; + + if (int err = lin_solver_.setMatrix(matrix_.get())) + return err; + + if (int err = lin_solver_.analyze()) + return err; + + // TODO: This function needs to be called to properly use a preconditioner in ReSolve (if there is any), but currently will error + // if there is no preconditioner configured. Once we can detect if a preconditioner is configured, we can restore this functionality. + // Also, we will always want to use *right* preconditioning. + // BUBBLE_FAIL(lin_solver_.preconditionerSetup("right")); + + return 0; + } + + template + int ResolveSystemSolver::setupSolver(bool reuse_factors) + { + if (reuse_factors) + { + if (int err = lin_solver_.refactorize()) + return err; + } + else + { + if (int err = lin_solver_.factorize()) + return err; + } + + return 0; + } + + template + int ResolveSystemSolver::solve(GridKit::LinearAlgebra::Vector& rhs, GridKit::LinearAlgebra::Vector& lhs) + { + ReSolve::vector::Vector resolve_rhs(rhs.getSize()); + ReSolve::vector::Vector resolve_lhs(lhs.getSize()); + + if (int err = resolve_rhs.setData(rhs.getData(), memspace_)) + return err; + + if (int err = resolve_lhs.setData(lhs.getData(), memspace_)) + return err; + + if (int err = lin_solver_.solve(&resolve_rhs, &resolve_lhs)) + return err; + + return 0; + } + + template class ResolveSystemSolver; + } // namespace LinearAlgebra +} // namespace GridKit \ No newline at end of file diff --git a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp new file mode 100644 index 000000000..3b529a9fa --- /dev/null +++ b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include + +#include + +#include +#include +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + template + class ResolveSystemSolver : public LinearSolver + { + public: + using RealT = LinearSolver::RealT; + + ResolveSystemSolver(ReSolve::SystemSolver& lin_solver, GridKit::memory::MemorySpace memspace = GridKit::memory::HOST); + + int configureSolver(GridKit::LinearAlgebra::CsrMatrix& matrix) final; + int setupSolver(bool reuse_factors = false) final; + int solve(GridKit::LinearAlgebra::Vector& rhs, GridKit::LinearAlgebra::Vector& lhs) final; + + private: + std::unique_ptr matrix_; + ReSolve::SystemSolver& lin_solver_; + ReSolve::memory::MemorySpace memspace_; + }; + } // namespace LinearAlgebra +} // namespace GridKit \ No newline at end of file diff --git a/README.md b/README.md index 16d5d946f..e80b68405 100644 --- a/README.md +++ b/README.md @@ -28,6 +28,7 @@ You should have all of the following installed before installing GridKit™ that the use of Enzyme is experimental, and some versions of it have been found to break GridKit code. - [LLVM](https://github.com/llvm/llvm-project) >= 15.x. GridKit is currently tested with LLVM 16. + - [ReSolve](https://github.com/ORNL/ReSolve) == commit `a93fb6571e3542b2fa779983fef0186a034c53c9` (optional) - [CMake](https://cmake.org/) >= 3.13 - C++ 20 compliant compiler From 259a3adfbc69be62a4d8b9fb5ba67f5a958d1674 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 13 Jul 2026 16:26:13 -0400 Subject: [PATCH 2/6] Documentation --- GridKit/LinearAlgebra/Solver/LinearSolver.hpp | 31 +++++++++++++++++++ .../Solver/ResolveSystemSolver.cpp | 17 ++++++++++ .../Solver/ResolveSystemSolver.hpp | 12 +++++++ 3 files changed, 60 insertions(+) diff --git a/GridKit/LinearAlgebra/Solver/LinearSolver.hpp b/GridKit/LinearAlgebra/Solver/LinearSolver.hpp index 1222f2543..5d06ce758 100644 --- a/GridKit/LinearAlgebra/Solver/LinearSolver.hpp +++ b/GridKit/LinearAlgebra/Solver/LinearSolver.hpp @@ -10,14 +10,45 @@ namespace GridKit namespace LinearAlgebra { + /** + * @brief An interface for linear solvers to be used in GridKit, such as in \ref Integrator::Rosenbrock. + * + */ template class LinearSolver { public: using RealT = GridKit::ScalarTraits::RealT; + /** + * @brief Configure the solver by analyzing the given matrix. The sparsity pattern of the given matrix must be set, and future calls + * to other methods will use this sparsity pattern, as well as the data pointer of the given matrix. + * + * @param matrix The matrix to configure the solver with. + * @return int An error code, or 0 if none. + */ virtual int configureSolver(GridKit::LinearAlgebra::CsrMatrix& matrix) = 0; + /** + * @brief Setup the solver with matrix data. + * + * @param reuse_factors If true, factors from a previous call to `setupSolver()` are reused. + * @pre \ref configureSolver() must be called first and every time the sparsity pattern of the matrix changes. + * @pre If `reuse_factors` is true, then `setupSolver` must have been called once before with the same sparsity pattern. + * @return int An error code, or 0 if none. + */ virtual int setupSolver(bool reuse_factors = false) = 0; + /** + * @brief Perform a linear solve using the configured matrix. + * + * @param rhs The right-hand-side vector (input). No data is changed. + * @param lhs The left-hand-side vector (output). + * @pre \ref setupSolver must be called first and every time the data of the matrix changes. + * @return int An error code, or 0 if none. + * + * @todo The data of the `rhs` parameter doesn't change, but we can't mark as cosntant because ReSolve currently requires + * a non-`const` pointer for vectors, even if it will only read from that vector. In the future this might change, so this + * can be updated to be `const`. + */ virtual int solve(GridKit::LinearAlgebra::Vector& rhs, GridKit::LinearAlgebra::Vector& lhs) = 0; }; } // namespace LinearAlgebra diff --git a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp index eba405a48..127a6119b 100644 --- a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp +++ b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp @@ -7,6 +7,10 @@ namespace GridKit { namespace LinearAlgebra { + /** + * @brief Helper function to convert GridKit `MemorySpace` into ReSolve `MemorySpace`. + * + */ ReSolve::memory::MemorySpace memorySpaceAsResolve(memory::MemorySpace memspace) { switch (memspace) @@ -25,6 +29,13 @@ namespace GridKit { } + /** + * @brief Creates a new ReSolve matrix with the dimensions of the given matrix, then sets its data pointers to those of the given matrix. + * Sets up the linear solver to use this new ReSolve matrix. + * + * @todo Right now preconditioning doesn't work. There should be a ReSolve PR soon for this. + * + */ template int ResolveSystemSolver::configureSolver(GridKit::LinearAlgebra::CsrMatrix& matrix) { @@ -64,6 +75,12 @@ namespace GridKit return 0; } + /** + * @brief Perform the solver by creating two new ReSolve vectors and setting their data to point to their GridKit counterparts. + * Then ask ReSolve to perform the solve. Since the matrix should be pointing at the correct GridKit buffers and the vectors + * as well, this solve should correctly fill `lhs`'s data buffer. + * + */ template int ResolveSystemSolver::solve(GridKit::LinearAlgebra::Vector& rhs, GridKit::LinearAlgebra::Vector& lhs) { diff --git a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp index 3b529a9fa..f666f9b4d 100644 --- a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp +++ b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp @@ -25,8 +25,20 @@ namespace GridKit int solve(GridKit::LinearAlgebra::Vector& rhs, GridKit::LinearAlgebra::Vector& lhs) final; private: + /** + * @brief The ReSolve matrix used for the linear solve. Its data pointers should just point to a \ref GridKit::LinearAlgebra::CsrMatrix. + * + */ std::unique_ptr matrix_; + /** + * @brief The ReSolve linear solver to use. + * + */ ReSolve::SystemSolver& lin_solver_; + /** + * @brief The memspace the data is stored in. + * + */ ReSolve::memory::MemorySpace memspace_; }; } // namespace LinearAlgebra From 36e2448bac1b03ec1529437aa11f82db6d7e13cc Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 13 Jul 2026 16:26:49 -0400 Subject: [PATCH 3/6] update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ec0df53b4..70e3416ef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -68,6 +68,7 @@ - Added IDA option to suppress algebraic variables in local error tests. - Removed `COO_Matrix` class. - Added portable `Vector` class and policy-based memory utilities. +- Add new `LinearSolver` interface for linear solvers. ## v0.1 From 1b80a8b83e8273d10fc1102a7916790043b8645e Mon Sep 17 00:00:00 2001 From: alexander-novo Date: Mon, 13 Jul 2026 20:32:30 +0000 Subject: [PATCH 4/6] Apply pre-commit fixes --- GridKit/LinearAlgebra/CMakeLists.txt | 2 +- GridKit/LinearAlgebra/Solver/CMakeLists.txt | 7 ++----- GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp | 5 +++-- GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp | 2 +- 4 files changed, 7 insertions(+), 9 deletions(-) diff --git a/GridKit/LinearAlgebra/CMakeLists.txt b/GridKit/LinearAlgebra/CMakeLists.txt index 0fe740023..a28565ed1 100644 --- a/GridKit/LinearAlgebra/CMakeLists.txt +++ b/GridKit/LinearAlgebra/CMakeLists.txt @@ -1,4 +1,4 @@ add_subdirectory(SparseMatrix) add_subdirectory(DenseMatrix) add_subdirectory(Vector) -add_subdirectory(Solver) \ No newline at end of file +add_subdirectory(Solver) diff --git a/GridKit/LinearAlgebra/Solver/CMakeLists.txt b/GridKit/LinearAlgebra/Solver/CMakeLists.txt index 356cc12e4..51da2bb1b 100644 --- a/GridKit/LinearAlgebra/Solver/CMakeLists.txt +++ b/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() \ No newline at end of file + LINK_LIBRARIES PRIVATE ReSolve::ReSolve) +endif() diff --git a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp index 127a6119b..4d6797147 100644 --- a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp +++ b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp @@ -25,7 +25,8 @@ namespace GridKit } template - ResolveSystemSolver::ResolveSystemSolver(ReSolve::SystemSolver& lin_solver, GridKit::memory::MemorySpace memspace) : lin_solver_(lin_solver), memspace_(memorySpaceAsResolve(memspace)) + ResolveSystemSolver::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; } // namespace LinearAlgebra -} // namespace GridKit \ No newline at end of file +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp index f666f9b4d..9ca55ef19 100644 --- a/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp +++ b/GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp @@ -42,4 +42,4 @@ namespace GridKit ReSolve::memory::MemorySpace memspace_; }; } // namespace LinearAlgebra -} // namespace GridKit \ No newline at end of file +} // namespace GridKit From 629c42a1ece517ac17242e12a52edefaeff1a816 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Wed, 15 Jul 2026 21:27:37 -0400 Subject: [PATCH 5/6] Add Rosenbrock integrators (#479) * Add rosenbrock * Add rosenbrock tests * Add RODAS5P tableau * Add rodas5p order test * Update Rosenbrock to use new LinearSolver interface * Moved Rosenbrock support classes to Native folder * Change Integrator namespace * Remove ResolveMemoryUtils.hpp --------- Co-authored-by: alexander-novo Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> Co-authored-by: pelesh --- CHANGELOG.md | 2 + CONTRIBUTING.md | 4 + GridKit/Solver/Dynamic/CMakeLists.txt | 29 +- .../Solver/Dynamic/Native/AdaptiveStep.cpp | 33 + .../Solver/Dynamic/Native/AdaptiveStep.hpp | 72 ++ GridKit/Solver/Dynamic/Native/CMakeLists.txt | 21 + GridKit/Solver/Dynamic/Native/ErrorNorm.hpp | 38 + GridKit/Solver/Dynamic/Native/FixedStep.cpp | 22 + GridKit/Solver/Dynamic/Native/FixedStep.hpp | 34 + GridKit/Solver/Dynamic/Native/InfNorm.cpp | 55 ++ GridKit/Solver/Dynamic/Native/InfNorm.hpp | 78 ++ GridKit/Solver/Dynamic/Native/StepControl.hpp | 27 + .../Solver/Dynamic/Native/StepController.hpp | 38 + GridKit/Solver/Dynamic/Rosenbrock.cpp | 912 ++++++++++++++++++ GridKit/Solver/Dynamic/Rosenbrock.hpp | 517 ++++++++++ GridKit/Solver/Dynamic/RosenbrockTableaus.cpp | 205 ++++ docs/CMakeLists.txt | 2 + docs/citations.bib | 10 + tests/UnitTests/Solver/CMakeLists.txt | 4 +- tests/UnitTests/Solver/Dynamic/CMakeLists.txt | 25 +- .../Solver/Dynamic/RosenbrockTests.hpp | 482 +++++++++ .../Solver/Dynamic/runRosenbrockTests.cpp | 18 + 22 files changed, 2616 insertions(+), 12 deletions(-) create mode 100644 GridKit/Solver/Dynamic/Native/AdaptiveStep.cpp create mode 100644 GridKit/Solver/Dynamic/Native/AdaptiveStep.hpp create mode 100644 GridKit/Solver/Dynamic/Native/CMakeLists.txt create mode 100644 GridKit/Solver/Dynamic/Native/ErrorNorm.hpp create mode 100644 GridKit/Solver/Dynamic/Native/FixedStep.cpp create mode 100644 GridKit/Solver/Dynamic/Native/FixedStep.hpp create mode 100644 GridKit/Solver/Dynamic/Native/InfNorm.cpp create mode 100644 GridKit/Solver/Dynamic/Native/InfNorm.hpp create mode 100644 GridKit/Solver/Dynamic/Native/StepControl.hpp create mode 100644 GridKit/Solver/Dynamic/Native/StepController.hpp create mode 100644 GridKit/Solver/Dynamic/Rosenbrock.cpp create mode 100644 GridKit/Solver/Dynamic/Rosenbrock.hpp create mode 100644 GridKit/Solver/Dynamic/RosenbrockTableaus.cpp create mode 100644 docs/citations.bib create mode 100644 tests/UnitTests/Solver/Dynamic/RosenbrockTests.hpp create mode 100644 tests/UnitTests/Solver/Dynamic/runRosenbrockTests.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 70e3416ef..c5d73ba7b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 1a446018e..00bece554 100644 --- a/CONTRIBUTING.md +++ b/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 diff --git a/GridKit/Solver/Dynamic/CMakeLists.txt b/GridKit/Solver/Dynamic/CMakeLists.txt index 2d7fa0152..e46292e41 100644 --- a/GridKit/Solver/Dynamic/CMakeLists.txt +++ b/GridKit/Solver/Dynamic/CMakeLists.txt @@ -5,13 +5,11 @@ # - Cameron Rutherford #]] -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) diff --git a/GridKit/Solver/Dynamic/Native/AdaptiveStep.cpp b/GridKit/Solver/Dynamic/Native/AdaptiveStep.cpp new file mode 100644 index 000000000..23cc0f306 --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/AdaptiveStep.cpp @@ -0,0 +1,33 @@ +#include "AdaptiveStep.hpp" + +#include +#include + +#include + +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 + StepControl AdaptiveStep::nextStep(RealT err, StepControl prev_step, uint8_t method_order) + { + StepControl next_step = prev_step; + + double h_mult = std::min(params_.fac_max_, std::max(params_.fac_scale_ * std::pow(err, GridKit::MINUS_ONE / method_order), params_.fac_min_)); + + next_step.accept_ = err <= 1; + next_step.step_size_ *= h_mult; + + return next_step; + } + + template class AdaptiveStep; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Native/AdaptiveStep.hpp b/GridKit/Solver/Dynamic/Native/AdaptiveStep.hpp new file mode 100644 index 000000000..11740696b --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/AdaptiveStep.hpp @@ -0,0 +1,72 @@ +#pragma once + +#include + +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 + class AdaptiveStep : public StepController + { + /** + * @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 nextStep(RealT err, StepControl 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 diff --git a/GridKit/Solver/Dynamic/Native/CMakeLists.txt b/GridKit/Solver/Dynamic/Native/CMakeLists.txt new file mode 100644 index 000000000..d42b58eb5 --- /dev/null +++ b/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) diff --git a/GridKit/Solver/Dynamic/Native/ErrorNorm.hpp b/GridKit/Solver/Dynamic/Native/ErrorNorm.hpp new file mode 100644 index 000000000..279e9255f --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/ErrorNorm.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include +#include +#include + +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 ErrorNorm + { + using State = GridKit::LinearAlgebra::Vector; + using RealT = typename GridKit::ScalarTraits::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& handler, GridKit::memory::MemorySpace memspace) const = 0; + }; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Native/FixedStep.cpp b/GridKit/Solver/Dynamic/Native/FixedStep.cpp new file mode 100644 index 000000000..965627ac0 --- /dev/null +++ b/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 + StepControl FixedStep::nextStep([[maybe_unused]] RealT err, StepControl prev_step, [[maybe_unused]] uint8_t method_order) + { + return StepControl{ + .accept_ = true, + .step_size_ = prev_step.step_size_, + }; + } + + template class FixedStep; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Native/FixedStep.hpp b/GridKit/Solver/Dynamic/Native/FixedStep.hpp new file mode 100644 index 000000000..82edda6cf --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/FixedStep.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include + +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 + class FixedStep : public StepController + { + StepControl nextStep(RealT err, StepControl 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 diff --git a/GridKit/Solver/Dynamic/Native/InfNorm.cpp b/GridKit/Solver/Dynamic/Native/InfNorm.cpp new file mode 100644 index 000000000..f1edc419e --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/InfNorm.cpp @@ -0,0 +1,55 @@ +#include "InfNorm.hpp" + +#include + +#include + +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 + InfNorm::RealT InfNorm::errorNorm(State& err, State& y, State& yprev, GridKit::LinearAlgebra::VectorHandler& 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, 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 diff --git a/GridKit/Solver/Dynamic/Native/InfNorm.hpp b/GridKit/Solver/Dynamic/Native/InfNorm.hpp new file mode 100644 index 000000000..6cc7a39a3 --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/InfNorm.hpp @@ -0,0 +1,78 @@ +#pragma once + +#include + +#include +#include + +namespace AnalysisManager +{ + namespace NativeDynamicSolver + { + + /** + * @brief The infinity error norm, which requires the error in every component in the system + * to meet tolerance. + * + */ + template + class InfNorm : public ErrorNorm + { + using State = GridKit::LinearAlgebra::Vector; + using RealT = ErrorNorm::RealT; + + /** + * @brief A workspace for the linear algebra operations required to calculate the norm. + * + */ + mutable struct + { + /** + * @brief The final vector which will have its norm taken. + * + */ + std::unique_ptr out_; + /** + * @brief The vector which will be used to scale the error in accordance with the tolerances. + * + */ + std::unique_ptr scale_; + /** + * @brief The absolute value of yprev. Used to calculate \ref scale_ + * + */ + std::unique_ptr yprev_abs_; + } workspace_; + + public: + /** + * @brief The configurable parameters of the error norm. + * + */ + struct Parameters + { + /** + * @brief A vector of absolute tolerances for each component. The norm will attempt to reject errors + * in each component above this tolerance. + * + */ + std::unique_ptr abs_tol_; + + /** + * @brief The relative tolerance. The norm will attempt to reject any error larger in percentage of + * the solution's maximum element than this. + * + */ + RealT rel_tol_; + } params_; + + InfNorm(Parameters&& params) + : params_(std::move(params)) + { + } + + RealT errorNorm(State& err, State& y, State& yprev, GridKit::LinearAlgebra::VectorHandler& handler, GridKit::memory::MemorySpace memspace) const final; + }; + + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Native/StepControl.hpp b/GridKit/Solver/Dynamic/Native/StepControl.hpp new file mode 100644 index 000000000..5dc678a4c --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/StepControl.hpp @@ -0,0 +1,27 @@ +#pragma once + +namespace AnalysisManager +{ + namespace NativeDynamicSolver + { + /** + * @brief Define control flow for `StepController`s to be able to control the step size of a `Rosenbrock` integrator. + * + */ + template + struct StepControl + { + /** + * @brief Whether or not the step is accepted. A rejected step will cause the time step controller to discard + * the next state and re-step with the new `step_size`. + * + */ + bool accept_; + /** + * @brief The step size the next step should take. + * + */ + RealT step_size_; + }; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Native/StepController.hpp b/GridKit/Solver/Dynamic/Native/StepController.hpp new file mode 100644 index 000000000..8c8abe5f8 --- /dev/null +++ b/GridKit/Solver/Dynamic/Native/StepController.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include + +#include + +namespace AnalysisManager +{ + namespace NativeDynamicSolver + { + /** + * @brief Interface for step size controllers. Used by `Rosenbrock` integrators to decide when to accept/reject steps and + * what size each step should be. + * + * @todo It may be best to have \ref usesError() return a reference to the \ref ErrorNorm that should be used. + */ + template + class StepController + { + public: + /** + * @brief Decide the control flow for the next step, based on information gathered by the integrator about the current step. + * + * @param err The estimated error made by the current step. Only calculated if `usesError()` returns true, otherwise it is assumed this + * value is not used. + * @param prev_step The control flow from before the current step was taken. Can be used to accurately update the step size. + * @param method_order The order of the method being used. + * @return StepControl + */ + virtual StepControl nextStep(RealT err, StepControl prev_step, uint8_t method_order) = 0; + /** + * @brief Return whether or not the `nextStep` method implementation uses the `err` parameter. If `false`, this parameter is not calculated. + * + */ + virtual bool usesError() const = 0; + }; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Rosenbrock.cpp b/GridKit/Solver/Dynamic/Rosenbrock.cpp new file mode 100644 index 000000000..8978bcd91 --- /dev/null +++ b/GridKit/Solver/Dynamic/Rosenbrock.cpp @@ -0,0 +1,912 @@ +#include "Rosenbrock.hpp" + +#include +#include +#include + +#include +#include + +/** + * @brief A small helper macro to "bubble" errors. The Rosenbrock implementations call many + * fallible model and linear algebra functions, and often the only thing that can be + * done when one of these fails is to fail the current function as well. This macro + * wraps a fallible function call and causes the current function to fail as well, + * returning the same error code as the called function. + * + */ +#define BUBBLE_FAIL(arg) \ + do \ + { \ + if (int err = (arg)) \ + return err; \ + } while (false) + +namespace AnalysisManager +{ + namespace NativeDynamicSolver + { + /** + * @brief Constructs a single-line report string out of the `StepInfo` object which is suitable to be included in a CSV file. + * + * Useful if you would like to plot step info. + * + */ + template + std::string Rosenbrock::StepInfo::csvReport() const + { + std::stringstream out; + out << std::scientific << std::setprecision(20) + << sim_time_ << ',' + << step_size_ << ',' + << next_step_size_ << ',' + << err_est_ << ',' + << step_no_ << ',' + << skip_lu_ << ',' + << skip_f_ << ',' + << accepted_; + + return out.str(); + } + + /** + * @brief Constructs a human-readable report string out of the `StepInfo` object. + * + * Useful for debug dumps. + * + */ + template + std::string Rosenbrock::StepInfo::report() const + { + std::stringstream out; + out << std::scientific << std::setprecision(20) + << "Simulation Time: " << sim_time_ << '\n' + << "Step Size: " << step_size_ << '\n' + << "Next Step Size: " << next_step_size_ << '\n' + << "Error Estimate: " << err_est_ << '\n' + << "Step Number: " << step_no_ << '\n' + << "Skip LU: " << skip_lu_ << '\n' + << "Skip F: " << skip_f_ << '\n' + << "Accepted: " << accepted_; + + return out.str(); + } + + /** + * @brief Constructs a human-readable report string out of the `Stats` object. + * + * Useful for reporting at the end of a simulation. + * + */ + template + std::string Rosenbrock::Stats::report() const + { + std::stringstream out; + out << "Rejections: " << rejections_.size() + << "\nSteps: " << num_steps_ + << "\nSkip LU Steps: " << skip_lu_steps_.size() + << "\nMin Step: " << min_step_ + << "\nMax Step: " << max_step_ + << "\nRHS Function Evals: " << f_evals_ + << "\nRHS Function Skipped: " << f_skipped_ + << "\nJacobian Evals: " << jac_evals_ + << "\nLinear Solves: " << decomp_solves_; + + return out.str(); + } + + /** + * @brief Accumulate statistics into one object. Useful to collect statistics over runs of several simulations, such as when + * there is a discrete event. + * + * @param other The other statistics to add to this one. + * + * @todo Right now, the step numbers for \ref rejections_ and \ref skip_lu_steps_ are impossible to tell apart from the different simulations. + */ + template + typename Rosenbrock::Stats& Rosenbrock::Stats::operator+=(const Stats& other) + { + rejections_.insert(rejections_.end(), other.rejections_.begin(), other.rejections_.end()); + skip_lu_steps_.insert(skip_lu_steps_.end(), other.skip_lu_steps_.begin(), other.skip_lu_steps_.end()); + + num_steps_ += other.num_steps_; + f_evals_ += other.f_evals_; + f_skipped_ += other.f_skipped_; + jac_evals_ += other.jac_evals_; + decomp_solves_ += other.decomp_solves_; + + min_step_ = std::min(min_step_, other.min_step_); + max_step_ = std::max(max_step_, other.max_step_); + + return *this; + } + + /** + * @brief Checks to see if the \ref asum_ variable can be re-used from the previous stage. + * + * Often, a row in the Rosenbrock A matrix is the exact same as the previous row, but with an additional + * non-zero element at the end. In this case, the \ref asum_ variable is the exact same as the previous stage, + * but with one extra additional term. The integrator can take advantage of this and reduce a matmul for computing + * \ref asum_ down to a single `axpy`. + * + * Typically, \ref asum_ is only initialized if it needs to be re-calculated from scratch. For this reason, this + * function will always return `false` for `stage == 0`, forcing \ref asum_ to be initialized for the first stage. + * + * @param stage The stage being checked. + */ + template + constexpr bool Rosenbrock::Tableau::canReuseAsum(size_t stage) const + { + assert(stage < num_stages_); + + if (stage == 0) + return false; + else + { + for (size_t j = 0; j < stage - 1; j++) + { + if (getA(stage, j) != getA(stage - 1, j)) + { + return false; + } + } + return true; + } + } + + /** + * @brief Checks to see if the \ref asum_ variable can be re-used from the last stage to compute the output state + * + * Often, the Rosenbrock m vector is the exact same as the last row in the A matrix, but with an additional + * non-zero element at the end. In this case, the output state for a step is equal to the \ref asum_ variables from + * the last stage but with a single weighted vector added to it. The integrator can take advantage of this and reduce + * a matmul for computing the final state down to a single `axpy`. + * + * If a method has only a single stage, then \ref asum_ will not have been initialized (as it will just be equal to y0). + * Therefore, this method will return `false`. + * + */ + template + constexpr bool Rosenbrock::Tableau::canReuseAsumForOut() const + { + if (num_stages_ == 1) + return false; + + for (size_t j = 0; j < num_stages_ - 1; j++) + { + if (getA(num_stages_ - 1, j) != m_[j]) + { + return false; + } + } + + return true; + } + + /** + * @brief Returns the index of a stage which can be used as an embedded error estimator, if that stage exists. + * + * Sometimes, Rosenbrock methods are designed in such a way where a stage can be used as an embedded error estimator. + * Typically, the embedded error estimator is a linear combination of the stages, so in this particular case the + * weights will be 1 on the estimator stage and 0 on all other stages. + * + * @pre This function is only valid to call if there is an embedded error estimator and its coefficients + * are included in this tableau. + */ + template + constexpr std::optional Rosenbrock::Tableau::errorEstimatorStage() const + { + assert(e_); + + std::optional re; + for (size_t j = 0; j < num_stages_; j++) + { + if (e_[j] == GridKit::ONE && !re) + { + re = j; + } + else if (e_[j] != GridKit::ZERO) + { + return {}; + } + } + return re; + } + + /** + * @brief Construct a new Rosenbrock integrator. + * + * @param tab The tableau to be used for this integrator. Since tableaus contain `std::unique_ptr`, it must be moved into the integrator. + * @param model The model to be simulated. Despite taking a pointer, this must be a valid pointer to an `Evaluator`. + * Must have \ref GridKit::Model::Evaluator::tag() set, and must be in Hessenberg form (\f(F(\dot y, y) = \dot y - f(y)\f)). + * @param lin_solver The linear solver to be used when constructing stages during simulation. The reference must remain valid for as long + * as the Rosenbrock integrator lives. + * @param vector_handler The vector handler to be used when simulating. The reference must remain valid for as long as the Rosenbrock + * integrator lives. + * @param err_norm The error norm to use in calculating error for the `StepController`. Will not be accessed if the `StepController` + * does not need error (such as `FixedStep`), so `nullptr` can be passed in that circumstance. + * @param memspace The memory space that linear algebra operations should be performed in. + */ + template + Rosenbrock::Rosenbrock(Tableau&& tab, + GridKit::Model::Evaluator* model, + GridKit::LinearAlgebra::LinearSolver& lin_solver, + GridKit::LinearAlgebra::VectorHandler& vector_handler, + const ErrorNorm* err_norm, + GridKit::memory::MemorySpace memspace) + : tab_(std::move(tab)), + model_(model), + lin_solver_(lin_solver), + vector_handler_(vector_handler), + err_norm_(err_norm), + memspace_(memspace) + { + } + + /** + * @brief Allocates memory based on the the size of \ref model_. Must be called before other methods. + * + * @note This method can fail. + * + * @pre Allocates the Jacobian matrix, which requires knowledge of the number of nonzero elements. This must be known at this point, + * so `allocate()` must be called beforehand on \ref model_ to count the nonzero elements and allocate the CSR Jacobian. + * + * @return An error code, with 0 as success. + */ + template + int Rosenbrock::allocate() + { + size_t size = static_cast(model_->size()); + + y_prev_ = std::make_unique(size); + y_cur_ = std::make_unique(size); + y_new_ = std::make_unique(size); + y_interp_ = std::make_unique(size); + workspace_.asum_ = std::make_unique(size); + workspace_.csum_ = std::make_unique(size); + workspace_.RHS_ = std::make_unique(size); + workspace_.RHS_first_stage_ = std::make_unique(size); + workspace_.dFdt_ = std::make_unique(size); + workspace_.mass_ = std::make_unique(size); + + BUBBLE_FAIL(y_prev_->allocate(memspace_)); + BUBBLE_FAIL(y_cur_->allocate(memspace_)); + BUBBLE_FAIL(y_new_->allocate(memspace_)); + BUBBLE_FAIL(y_interp_->allocate(memspace_)); + BUBBLE_FAIL(workspace_.asum_->allocate(memspace_)); + BUBBLE_FAIL(workspace_.csum_->allocate(memspace_)); + BUBBLE_FAIL(workspace_.RHS_->allocate(memspace_)); + BUBBLE_FAIL(workspace_.RHS_first_stage_->allocate(memspace_)); + BUBBLE_FAIL(workspace_.dFdt_->allocate(memspace_)); + BUBBLE_FAIL(workspace_.mass_->allocate(memspace_)); + + if (tab_.e_) + { + std::optional err_est_stage = tab_.errorEstimatorStage(); + if (!err_est_stage) + { + workspace_.err_est_ = std::make_unique(size); + BUBBLE_FAIL(workspace_.err_est_->allocate(memspace_)); + } + } + + workspace_.stages_ = std::make_unique[]>(tab_.num_stages_); + for (size_t i = 0; i < tab_.num_stages_; i++) + { + workspace_.stages_[i] = std::make_unique(size); + BUBBLE_FAIL(workspace_.stages_[i]->allocate(memspace_)); + workspace_.stages_[i]->setToZero(memspace_); + } + + if (tab_.order_ > 2) + { + dense_coeff_ = std::make_unique[]>(tab_.order_ - 2); + for (size_t i = 0; i < static_cast(tab_.order_ - 2); i++) + { + dense_coeff_[i] = std::make_unique(size); + BUBBLE_FAIL(dense_coeff_[i]->allocate(memspace_)); + } + } + + return 0; + } + + /** + * @brief Initializes the simulation. Must be called before \ref integrate() or \ref timeStep(). Must also be called after + * discrete events. + * + * - Sets the simulation time to `t0` and copies the initial condition from \ref model_. + * - Analyzes \ref model_ Jacobian sparsity and runs the preconditioner + * - Generates the mass matrix from \ref GridKit::Model::Evaluator::tag(). If the tag is not properly set, then initialization will fail. + * - Resets \ref stats_. + * + * @note This method can fail. + * + * @pre Must have called \ref allocate(). The `model.tag_` variable must be properly constructed. + * + * @todo Document mass matrix construction + * + * @param t0 The starting simulation time. + * @return An error code, with 0 as success. + */ + template + int Rosenbrock::initializeSimulation(RealT t0) + { + current_time_ = t0; + BUBBLE_FAIL(y_cur_->copyFromExternal(model_->y().data(), memspace_, memspace_)); + workspace_.jacobian_analyzed_ = false; + + BUBBLE_FAIL(lin_solver_.configureSolver(*model_->getCsrJacobian())); + + if (model_->tag().size() != static_cast(model_->size())) + { + std::cerr << "Model tag is either unset or does not match the size of the model\n"; + return 1; + } + + std::unique_ptr mass = std::make_unique(model_->tag().size()); + for (size_t i = 0; i < static_cast(model_->size()); i++) + { + mass[i] = model_->tag()[i] ? GridKit::ONE : GridKit::ZERO; + } + BUBBLE_FAIL(workspace_.mass_->copyFromExternal(mass.get(), memspace_, memspace_)); + + stats_ = Stats(); + + return 0; + } + + /** + * @brief Simulate the given model, producing output at the given output times. + * + * Implements a simple time stepping algorithm and facilitates output. + * 1. Calls \ref timeStep() to move the simulation forward by \ref step_size_ (starting at `params.starting_step`). + * 2. Consults the `step_controller` to see if the step is accepted, and what the next \ref step_size_ should be. + * 3. If accepted, advance simulation by adding \ref step_size_ to \ref current_time_ using Kahan summation and shifting + * \ref y_new_ -> \ref y_cur_ -> \ref y_prev_. + * 4. Any time we step over an output time, calculate dense coefficients if available. + * 5. For each output time stepped over in the last step, interpolate the output at that time and call `out_cb`. + * + * Time is advanced to at least the final output time. The simulation may step over the final time, so if you + * wish to restart simulation from the final output time, make sure to re-initialize the model's state and call + * \ref initializeSimulation(). + * + * @note This method can fail. + * + * @pre Must have called \ref initializeSimulation(). + * + * @todo A higher-order Hermite interpolation can be used as a fallback for interpolation when dense coefficients don't exist. + * Each step needs to sample the residual function at the beginning (and therefore end) of the step, which contains derivative + * information. This derivative information can be re-used for Hermite interpolation. + * + * @todo It doesn't really make sense to have the error estimator separate from the step controller, since your choice of one will + * affect your choice of the other. The error estimator should probably be inside the step controller, then accessed if needed. + * + * @todo Return an error when max steps is hit. + * + * @todo Check if current time is close enough to output time, and skip interpolation + * + * @todo Configure upper bound for skip lu step size increase + * + * @param out_times The times at which output is wanted. The simulation will stop once the final output time has been reached. + * @param step_controller The step size controller to use during the simulation. + * @param params The parameters to use during the simulation. + * @param out_cb An optional function which, if provided, will be called once for each time in `out_times`. The simulation time of the + * output is passed as the only argument. Before being called, the model will be updated with the output state and can be queried + * separately by the callback if needed. + * @param step_cb An optional function which, if provided, will be called once for each step the integrator takes. Information about the + * step which was just taken is provided as the argument. Before being called, the model will be updated with the current value of the state + * an can be queried separately by the callback if needed. Useful for debugging simulations. + * @return An error code, with 0 as success. + */ + template + int Rosenbrock::integrate(const std::vector& out_times, + StepController& step_controller, + Parameters params, + std::optional> out_cb, + std::optional> step_cb) + { + constexpr RealT ONE = GridKit::ONE; + constexpr RealT ZERO = GridKit::ZERO; + + skip_lu_ = false; + skip_f_ = false; + + bool prev_accept = true; + step_size_ = params.starting_step_; + + double next_step_size; + + // Kahan summation time buffer. The "leftover" time that was lost when trying to add h at some point that needs to be added + // later + double time_buffer = 0; + + // Generate output for each output time + for (double out_time : out_times) + { + while (current_time_ < out_time && stats_.num_steps_ < params.max_steps_) + { + BUBBLE_FAIL(timeStep(current_time_, step_size_)); + + double err = 0; + + if (step_controller.usesError()) + { + if (err_norm_ == nullptr) + { + std::cerr << "The provided step controller requires the use of an error norm, but none was provided!\n"; + + return -1; + } + + State& err_vec = errorEstimate(); + err = err_norm_->errorNorm(err_vec, *y_new_, *y_cur_, vector_handler_, memspace_); + } + + StepControl next_step = step_controller.nextStep(err, + StepControl{ + .accept_ = prev_accept, + .step_size_ = step_size_, + }, + tab_.order_); + prev_accept = next_step.accept_; + next_step_size = next_step.step_size_; + + if (prev_accept) + { + // Try to add the leftover time that we've stored up + double step_size_adj = step_size_ + time_buffer; + double next_time = current_time_ + step_size_adj; + + // Kahan summation - keep track of how much of step_size_adj we weren't able to add to current_time + // due to lack of precision + time_buffer = step_size_adj - (next_time - current_time_); + current_time_ = next_time; + + // Since time is advancing, we need to re-evaluate the residual function and dense coefficients + skip_f_ = false; + dense_coefficients_valid_ = false; + + stats_.num_steps_++; + if (skip_lu_) + { + stats_.skip_lu_steps_.push_back(StepInfo{ + .sim_time_ = current_time_, + .step_size_ = step_size_, + .next_step_size_ = next_step_size, + .err_est_ = err, + .step_no_ = stats_.num_steps_, + .skip_lu_ = skip_lu_, + .skip_f_ = skip_f_, + .accepted_ = prev_accept, + }); + } + stats_.min_step_ = std::min(stats_.min_step_, step_size_); + stats_.max_step_ = std::max(stats_.max_step_, step_size_); + + // Shift y_new_ -> y_cur_ -> y_prev_ + // Then y_new_ is free to be replaced (contains old y_prev_) + std::swap(y_prev_, y_cur_); + std::swap(y_cur_, y_new_); + } + else + { + skip_f_ = true; + + stats_.rejections_.push_back(StepInfo{ + .sim_time_ = current_time_, + .step_size_ = step_size_, + .next_step_size_ = next_step_size, + .err_est_ = err, + .step_no_ = stats_.num_steps_, + .skip_lu_ = skip_lu_, + .skip_f_ = skip_f_, + .accepted_ = prev_accept, + }); + } + + // Check if we can use time delay Jacobians. If we would have increased the step size (but not too much), + // instead keep the step size the same and use time-delay Jacobian. + // TODO: configure upper bound here + double step_gain = next_step_size / step_size_; + if (params.skip_lu_ && step_gain >= ONE && step_gain <= 1.2) + { + skip_lu_ = true; + } + else + { + skip_lu_ = false; + prev_step_size_ = step_size_; + step_size_ = next_step_size; + } + + // If there is a step_cb, update the model state and call it + if (step_cb) + { + BUBBLE_FAIL(y_cur_->copyToExternal(model_->y().data(), memspace_, memspace_)); + model_->updateTime(current_time_, ZERO); + + (*step_cb)(StepInfo{ + .sim_time_ = current_time_, + .step_size_ = step_size_, + .next_step_size_ = next_step_size, + .err_est_ = err, + .step_no_ = stats_.num_steps_, + .skip_lu_ = skip_lu_, + .skip_f_ = skip_f_, + .accepted_ = prev_accept, + }); + } + } + + // Check to make sure integration was paused because we reached an output time. + // Other reasons, like hitting max step count, shouldn't generate output. + if (current_time_ >= out_time) + { + // Theta = (t - t0) / h = (t - t1) / h + 1 + // current_time_ is t1 here + RealT theta = (out_time - current_time_) / prev_step_size_ + ONE; + + // Generate output at the appropriate time. + if (tab_.hasDenseOutput()) + { + if (!dense_coefficients_valid_) + { + BUBBLE_FAIL(calcDenseCoeff()); + dense_coefficients_valid_ = true; + } + + BUBBLE_FAIL(interpDense(theta)); + } + else + { + // TODO: Put code for alternative interpolation (Abdou) here + BUBBLE_FAIL(y_interp_->copyFromExternal(*y_prev_, memspace_, memspace_)); + vector_handler_.scal(1 - theta, y_interp_.get(), memspace_); + vector_handler_.axpy(theta, y_cur_.get(), y_interp_.get(), memspace_); + } + + // Update model with output and call out_cb if it exists + if (out_cb) + { + BUBBLE_FAIL(y_interp_->copyToExternal(model_->y().data(), memspace_, memspace_)); + model_->updateTime(out_time, ZERO); + + (*out_cb)(out_time); + } + } + else + { + BUBBLE_FAIL(y_interp_->copyFromExternal(*y_cur_, memspace_, memspace_)); + break; + } + } + + return 0; + } + + /** + * @brief Advance the simulation forward by one step, storing the new state in \ref y_new_. + * + * Apply the Rosenbrock scheme using the stored tableau. Each stage \f(u_i\f) is calculated as + * + * \f[\left(J - \frac{1}{h\gamma} M\right)u_i = -f \left(t_0 + \alpha_ih, y_0 + \sum_{j = 1}^{i - 1} a_{ij}u_j\right) - M \sum_{j=1}^{i-1} \left(\frac{c_{ij}}{h}\right)u_j,\f] + * + * and the next state \f(y_1\f) is calculated as + * + * \f[y_1 = y_0 + \sum_{j = 1}^sm_ju_j\f] + * + * where the coefficients \f(\gamma, \alpha_i, a_{ij}, c_{ij}, m_j\f) come from the tableau, and the mass matrix \f(M\f) and Jacobian \f(J\f) come from the model. + * + * This method uses some state which is maintained between calls for future calls to `timeStep()` and to communicate with \ref integrate(): + * - \ref y_cur_ is used as \f(y_0\f) and \ref y_new_ is used as \f(y_1\f) + * - \ref asum_ is used as \f(y_0 + \sum_{j = 1}^{i - 1} a_{ij}u_j\f) for every stage except the first. Often, \f(a_{ij} = a_{i-1,j}\f), so this variable + * can be re-used between stages to save computation. Similarly, often \f(a_{ij} = m_j\f), so this variable can be re-used for computing \ref y_new_. + * - \ref csum_ is used as \f(M \sum_{j=1}^{i-1} \left(\frac{c_{ij}}{h}\right)u_j\f) for every stage except the first + * - \ref RHS_ is used as \f(-f \left(t_0 + \alpha_ih, y_0 + \sum_{j = 1}^{i - 1} a_{ij}u_j\right) - M \sum_{j=1}^{i-1} \left(\frac{c_{ij}}{h}\right)u_j\f), i.e. + * the residual evaluated at \ref asum_ plus \ref csum_, for every stage but the first. This is used as the right-hand side vector for the linear solve to solve + * for the stage value \f(u_i\f). + * - \ref RHS_first_stage_ is used as a special \ref RHS_ for the first stage, since it may need to be saved for the next step for the \ref skip_f_ flag. + * Since \f(\alpha_1 = 0\f) always, this will have a value of \f(f \left(t_0, y_0\right)\f). + * - \ref stages_ stores all of the stages \f(u_i\f). These stages are necessary to be used in future calls to \ref errorEstimate() and \ref calcDenseCoeff(). + * - \ref skip_lu_ is used as a flag set by \ref integrate() to indicate that it is appropriate to use a time-delay Jacobian by re-using the factorization + * of the last step. + * - \ref skip_f_ is used as a flag set by \ref integrate() to indicate that \f(t_0, y_0\f) have not changed since the last time `timeStep()` was called, + * so \ref RHS_first_stage_ can be re-used from the previous step. This will only happen when a step was rejected, so \ref skip_lu_ should always be false, + * and the entire first stage can't be re-used. + * - \ref jacobian_analyzed_ keeps track of whether the Jacobian has been factored in a previous call to `timeStep()`. If so, then it can be re-factored + * in a faster way. The first factor must be done on actual data, so it cannot be performed pre-simulation. + * + * @pre Must call \ref initializeSimulation() beforehand. + * + * @note This method can fail. + * + * @todo This currently does not work with non-autonomous models. Some thought needs to be put in for how we want non-autonomous models to work. + * + * @todo It doesn't really make sense to pass `t0` here, which may be inconsistent with \ref current_time_. + * This function should just use \ref current_time_ in place of `t0`. + * + * @todo It doesn't really make sense for this method to be `public` since it relies on proper setup via \ref integrate(). + * + * @todo May be able to move the copying of model jacobian pointers to \ref allocate(). + * + * @todo Since \ref csum_ is multiplied by the mass matrix (which is currently diagonal with 0s for algebraic components), we can save some computation + * by only computing the differential parts of \ref csum_ and storing it in the same variable as \ref RHS_. + * + * @param t0 \f(t_0\f) in the above formula + * @param dt \f(h\f) in the above formula. The next state \f(y_1\f) will be an estimate of the state at \f(t_1 = t_0 + h\f). + * @return An error code, with 0 as success. + */ + template + int Rosenbrock::timeStep(RealT t0, RealT dt) + { + constexpr RealT ZERO = GridKit::ZERO; + constexpr RealT MINUS_ONE = GridKit::MINUS_ONE; + + // A flag to keep track of if y0 (stored in y_cur_) has been copied in to the model already, to avoid double-copying + // for evaluating the Jacobian and residual on stage 1 (both evaluated at y0). + bool y0_copied = false; + + // Form the left-hand side of the system. This is constant between stages. + // Can sometimes be skipped if the method allows for time-delay Jacobians (such as w-methods). + [[likely]] + if (!tab_.is_w_ || !skip_lu_) + { + BUBBLE_FAIL(y_cur_->copyToExternal(model_->y().data(), memspace_, memspace_)); + y0_copied = true; + + // GridKit, like IDA, expects to evaluate the Jacobian J = df/dy + alpha * df/dy', + // so we need a negative here since df/dy' = M. + model_->updateTime(t0, MINUS_ONE / (dt * tab_.gamma_)); + BUBBLE_FAIL(model_->evaluateJacobian()); + BUBBLE_FAIL(lin_solver_.setupSolver(workspace_.jacobian_analyzed_)); + workspace_.jacobian_analyzed_ = true; + + stats_.jac_evals_++; + } + + // First stage + [[unlikely]] + if (skip_f_) + { + stats_.f_skipped_++; + } + else + { + // TODO: non-autonomous model + if (!y0_copied) + { + BUBBLE_FAIL(y_cur_->copyToExternal(model_->y().data(), memspace_, memspace_)); + y0_copied = true; + } + model_->updateTime(t0, ZERO); + BUBBLE_FAIL(model_->evaluateResidual()); + BUBBLE_FAIL(workspace_.RHS_first_stage_->copyFromExternal(model_->getResidual().data(), memspace_, memspace_)); + vector_handler_.scal(-1, workspace_.RHS_first_stage_.get(), memspace_); + + stats_.f_evals_++; + } + BUBBLE_FAIL(lin_solver_.solve(*workspace_.RHS_first_stage_, *workspace_.stages_[0])); + stats_.decomp_solves_++; + + // Rest of stages + for (size_t i = 1; i < tab_.num_stages_; i++) + { + // Calculate asum + // We can sometimes reuse asum from the previous stage + if (i > 1 && tab_.canReuseAsum(i)) + { + if (tab_.A_[tab_.num_stages_ * i + i - 1] != ZERO) + vector_handler_.axpy(tab_.A_[tab_.num_stages_ * i + i - 1], workspace_.stages_[i - 1].get(), workspace_.asum_.get(), memspace_); + } + else + { + BUBBLE_FAIL(workspace_.asum_->copyFromExternal(*y_cur_, memspace_, memspace_)); + + for (size_t j = 0; j < i; j++) + { + if (tab_.A_[tab_.num_stages_ * i + j] != ZERO) + vector_handler_.axpy(tab_.A_[tab_.num_stages_ * i + j], workspace_.stages_[j].get(), workspace_.asum_.get(), memspace_); + } + } + + // Calculate csum + // TODO: Since csum is multiplied by the mass matrix, we can reduce calculations by just not calculating some indices + BUBBLE_FAIL(workspace_.csum_->setToZero(memspace_)); + for (size_t j = 0; j < i; j++) + { + if (tab_.C_[i * tab_.num_stages_ + j] != ZERO) + { + vector_handler_.axpy(tab_.C_[i * tab_.num_stages_ + j] / dt, workspace_.stages_[j].get(), workspace_.csum_.get(), memspace_); + } + } + + // TODO: non-autonomous model + BUBBLE_FAIL(workspace_.asum_->copyToExternal(model_->y().data(), memspace_, memspace_)); + model_->updateTime(t0 + tab_.alpha_sum_[i] * dt, ZERO); + BUBBLE_FAIL(model_->evaluateResidual()); + workspace_.RHS_->copyFromExternal(model_->getResidual().data(), memspace_, memspace_); + + vector_handler_.scal(MINUS_ONE, workspace_.RHS_.get(), memspace_); + vector_handler_.scal(workspace_.mass_.get(), workspace_.csum_.get(), memspace_); + vector_handler_.axpy(MINUS_ONE, workspace_.csum_.get(), workspace_.RHS_.get(), memspace_); + + BUBBLE_FAIL(lin_solver_.solve(*workspace_.RHS_, *workspace_.stages_[i])); + stats_.f_evals_++; + stats_.decomp_solves_++; + } + + // Compute the solution at time t + dt + // It happens often where the solution is just asum of the last stage + // plus some multiple of the last stage. In that case we can avoid a matmul + if (tab_.canReuseAsumForOut()) + { + std::swap(workspace_.asum_, y_new_); + vector_handler_.axpy(tab_.m_[tab_.num_stages_ - 1], workspace_.stages_[tab_.num_stages_ - 1].get(), y_new_.get(), memspace_); + } + else + { + BUBBLE_FAIL(y_new_->copyFromExternal(*y_cur_, memspace_, memspace_)); + + for (size_t j = 0; j < tab_.num_stages_; j++) + { + if (tab_.m_[j] != ZERO) + { + vector_handler_.axpy(tab_.m_[j], workspace_.stages_[j].get(), y_new_.get(), memspace_); + } + } + } + + return 0; + } + + /** + * @brief Calculates an estimation of the error produced by the last call to \ref timeStep() which can be used as the + * `err` argument for \ref ErrorNorm::errorNorm(). + * + * Calculate the embedded error as + * + * \f[\hat{e} = \sum_{j = 1}^s e_j u_j,\f] + * + * where \f(e_j\f) are tableau coefficients and \f(u_j\f) are stages computed by \ref timeStep(). It happens often that + * \f(e_j = 0\f) for all but one stage, where \f(e_j = 1\f). In that case, the stage itself is used as the error estimate, + * and extra calculation can be avoided. For this reason, a reference to the estimate is returned to avoid an unnecessary copy. + * + * @pre Must call \ref timeStep(), and tableau must have coefficients for an embedded error estimator. + * @note This method can fail. + * + * @todo This function is fallible, but the return type makes it difficult to return an error code. Right now it will throw an + * error, but it should be refactored to allow returning an error code, such as by returning `std::variant`. + * + * @return A reference to the estimated error. + */ + template + Rosenbrock::State& Rosenbrock::errorEstimate() const + { + // Test to see if the tableau allows us to use a stage as the error estimate, + // avoiding extra computation. + std::optional err_stage = tab_.errorEstimatorStage(); + + if (err_stage) + { + return *workspace_.stages_[*err_stage]; + } + else + { + // TODO: could make this function return recoverable errors by using std::variant + int err_code = workspace_.err_est_->copyFromExternal(*workspace_.stages_[0], memspace_, memspace_); + + if (err_code) + { + throw std::format("ReSolve::vector::Vector::copyFromExternal failed with error code {}", err_code); + } + + vector_handler_.scal(tab_.e_[0], workspace_.err_est_.get(), memspace_); + for (size_t j = 1; j < tab_.num_stages_; j++) + { + if (tab_.e_[j] != GridKit::ZERO) + { + vector_handler_.axpy(tab_.e_[j], workspace_.stages_[j].get(), workspace_.err_est_.get(), memspace_); + } + } + + return *workspace_.err_est_; + } + } + + /** + * @brief Calculate the interpolation nodes used by \ref interpDense() based on \ref stages_ computed by + * the last call to \ref timeStep(). Nodes are stored in \ref dense_coeff_. Only needs to be called once, + * but can be invalidated by future calls to \ref timeStep(). \ref integrate() keeps track of \ref dense_coefficients_valid_, + * which tells you if this function is needed to be called. + * + * @pre Must call \ref timeStep() first. + * + * @note This method can fail. + * + * @return An error code, with 0 as success. + */ + template + int Rosenbrock::calcDenseCoeff() + { + if (tab_.order_ > 2) + { + for (size_t j = 0; j < static_cast(tab_.order_ - 2); j++) + { + BUBBLE_FAIL(dense_coeff_[j]->setToZero(memspace_)); + } + + for (size_t i = 0; i < tab_.num_stages_; i++) + { + for (size_t j = 0; j < static_cast(tab_.order_ - 2); j++) + { + vector_handler_.axpy(tab_.H_[j * tab_.num_stages_ + i], workspace_.stages_[i].get(), dense_coeff_[j].get(), memspace_); + } + } + } + + return 0; + } + + /** + * @brief Calculate an interpolated state at \f(\theta = \frac{t - t_0}{h}\f) between the initial state + * \f(y_0\f) and final state \f(y_1\f) of the last step taken using dense interpolation nodes calculated + * by \ref calcDenseCoeff(). + * + * For a valid interpolation of appropriate order, \f(\theta\f) must be in \f([0, 1]\f), although values + * beyond 1 can be used for extrapolation if desired. + * + * Uses a number of interpolation nodes equal to the order. Since \f(y_0\f) and \f(y_1\f) are interpolation + * nodes, this method will only access \ref dense_coeff_ if the method's order is greater than 2. + * + * The inteporlation is calculated as + * + * \f[y(\theta) = (1 - \theta) y_0 + \theta \left(y_1 + (1 - \theta) \sum_{i = 1}^{p-2} \theta^{i-1} \hat{y}_i\right),\f] + * + * where \f(\hat{y}_i\f) are the dense interpolation nodes calculated in \ref calcDenseCoeff() and \f(p\f) is the order of the method. + * This calculation is carried out using synthetic division. + * + * @pre Must call \ref calcDenseCoeff() first. + * + * @note This method can fail. + * + * @note If `theta` is very close to 0 or very close to 1, it might be better to simply use \ref y_prev_ or \ref y_cur_ instead of + * calculating this interpolation. + * + * @param theta The fraction of time during the last step taken to calculate the interpolation at. \f(\theta = \frac{t - t_0}{h}\f) + * @return An error code, with 0 as success. + */ + template + int Rosenbrock::interpDense(RealT theta) + { + constexpr RealT ONE = GridKit::ONE; + + if (tab_.order_ > 2) + { + BUBBLE_FAIL(y_interp_->copyFromExternal(*dense_coeff_[tab_.order_ - 3], memspace_, memspace_)); + + for (size_t i = 1; i < static_cast(tab_.order_ - 2); i++) + { + vector_handler_.scal(theta, y_interp_.get(), memspace_); + vector_handler_.axpy(ONE, dense_coeff_[tab_.order_ - 3 - i].get(), y_interp_.get(), memspace_); + } + + // TODO: This scal can be removed and absorbed into the next axpy, except that it currently isn't possible to put a scalar + // multiple on the y term. + vector_handler_.scal(ONE - theta, y_interp_.get(), memspace_); + } + else + { + BUBBLE_FAIL(y_interp_->setToZero(memspace_)); + } + + vector_handler_.axpy(ONE, y_cur_.get(), y_interp_.get(), memspace_); + vector_handler_.scal(theta, y_interp_.get(), memspace_); + vector_handler_.axpy(ONE - theta, y_prev_.get(), y_interp_.get(), memspace_); + + return 0; + } + + template class Rosenbrock; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/Rosenbrock.hpp b/GridKit/Solver/Dynamic/Rosenbrock.hpp new file mode 100644 index 000000000..234d94fb3 --- /dev/null +++ b/GridKit/Solver/Dynamic/Rosenbrock.hpp @@ -0,0 +1,517 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace AnalysisManager +{ + namespace NativeDynamicSolver + { + /** + * @brief A Rosenbrock integrator designed to simulate a `GridKit::Model::Evaluator` model. + * Rosenbrock integrators are best for quick medium-accuracy simulation and + * boast a large amount of customization to simulating a given model. + * + * For the list of available Rosenbrock methods, see `Rosenbrock::Tableau`. + */ + template + class Rosenbrock + { + using RealT = typename GridKit::ScalarTraits::RealT; + using State = GridKit::LinearAlgebra::Vector; + + public: + /** + * @brief Keeps track of a variety of notable properties of a single step. + * + */ + struct StepInfo + { + /** + * @brief The simulation time at the beginning of the step + * + */ + RealT sim_time_; + /** + * @brief The size of the step. + * + */ + RealT step_size_; + /** + * @brief The size of the next step, as governed by the current `StepController` in use. + * + */ + RealT next_step_size_; + /** + * @brief The estimated error made by the step, as calculated by the current `ErrorNorm` in use. + * + */ + RealT err_est_; + /** + * @brief The step number, starting at 1. + * + */ + size_t step_no_; + /** + * @brief Whether or not the integrator decided to skip computing the decomposition of the Jacobian on this step. + * + */ + bool skip_lu_; + /** + * @brief Whether or not the integrator decided to skip evaluating the residual on the first stage on this step. + * + */ + bool skip_f_; + /** + * @brief Whether or not this step was accepted by the `StepController` in use. + * + */ + bool accepted_; + + std::string csvReport() const; + std::string report() const; + }; + + /** + * @brief Keeps track of running statistics of the integrator since the last `initializeSimulation()` call. + * + */ + struct Stats + { + /** + * @brief Information of each step which has been rejected. + * + */ + std::vector rejections_; + /** + * @brief Information of each step which the integrator decided to skip re-factoring the Jacobian. + * + */ + std::vector skip_lu_steps_; + /** + * @brief How many steps the integrator has taken. + * + */ + size_t num_steps_ = 0; + /** + * @brief Number of model residual function evaluations. + * + */ + size_t f_evals_ = 0; + /** + * @brief Number of model residual function evaluations which have been skipped by the integrator. + * + */ + size_t f_skipped_ = 0; + /** + * @brief Number of model Jacobian evaluations. + * + */ + size_t jac_evals_ = 0; + /** + * @brief Number of linear solves against the model Jacobian. + * + */ + size_t decomp_solves_ = 0; + /** + * @brief Minimum step size. + * + */ + RealT min_step_ = INFINITY; + /** + * @brief Maximum step size. + * + */ + RealT max_step_ = 0; + + std::string report() const; + Stats& operator+=(const Stats& other); + }; + + /** + * @brief Parameters of the integrator. + * + */ + struct Parameters + { + /** + * @brief What step size the first step should take. + * + * @todo Consider adding a starting step size selector to select this automatically. + */ + RealT starting_step_ = 1e-5; + /** + * @brief The maximum number of steps the integrator should take. If the integrator has not reached the final time before + * taking this many steps, then integration is stopped. For more details, see `integrate()`. + * + */ + size_t max_steps_ = 2000; + /** + * @brief Whether or not the integrator should attempt to skip Jacobian decompositions. + * + * @note This feature is only available if the underlying method is a Rosenbrock-W method and can speed up + * the time taken to compute each step. However, the overall number of steps taken will increase. + * + */ + bool skip_lu_ = false; + }; + + /** + * @brief A list of the coefficients needed to complete Rosenbrock integration in an accurate way. + * Each tableau can be considered a different method or scheme, and will have different properties, + * advantages, and disadvantages. + * + */ + struct Tableau + { + /** + * @brief The number of stages used by the method. Each stage requires one model residual evaluation. + * + */ + size_t num_stages_; + /** + * @brief The coefficient along the diagonal of the Gamma matrix. + * + */ + RealT gamma_; + /** + * @brief A vector of sums of rows of the alpha matrix. These are the classic + * Runge-Kutta 'c' coefficients, or abscissae. The size of this vector + * should be equal to `num_stages`. + * + */ + std::unique_ptr alpha_sum_; + /** + * @brief A vector of sums of rows of the Gamma matrix. The size of this vector + * should be equal to `num_stages`. + * + */ + std::unique_ptr gamma_sum_; + /** + * @brief A vector of weights for constructing the final solution from the stages. + * The size of this vector should be equal to `num_stages`. + * + */ + std::unique_ptr m_; + /** + * @brief OPTIONAL vector of coefficients for the embedded error method. If it exists, + * the size of this vector should be equal to `num_stages`. + * + */ + std::unique_ptr e_; + /** + * @brief The transformed A coefficient matrix. Strictly lower triangular and stored in dense row-major form. + * Upper triangular terms are not accessed. Should be `num_stages` by `num_stages` large. + * + */ + std::unique_ptr A_; + /** + * @brief The transformed C coefficient matrix. Strictly lower triangular and stored in dense row-major form. + * Upper triangular terms are not accessed. Should be `num_stages` by `num_stages` large. + * + */ + std::unique_ptr C_; + /** + * @brief OPTIONAL matrix of dense coefficients. Defines how the stages should be transformed into interpolant + * nodes for computing dense output. The interpolating polynomial has an order one less than the order of + * the method, and two interpolant nodes are already pre-computed, so if this matrix exists it should be + * `order` - 2 by `num_stages` large. + * + */ + std::unique_ptr H_; + /** + * @brief What ODE order these coefficients satisfy. If `is_dae` is true, then the coefficients must additionally satisfy + * DAE conditions up to this order. If `is_w` is true, then the coefficients must additionally satisfy ROW conditions + * up to this order. If `is_krylov` is true, then the coefficients must additionally satisfy ROK condition up to this order. + * + */ + + uint8_t order_; + /** + * @brief Whether or not these coefficients are appropriate to use in a Rosenbrock-Krylov (ROK) solver. + * + */ + bool is_krylov_; + /** + * @brief Whether or not these coefficients satisfy Rosenbrock-W (ROW) order conditions up to `order`. + * The integrator may take advantage of this fact by e.g. using time-delay Jacobians to speed up computation. + * + */ + bool is_w_; + /** + * @brief Whether or not these coefficients satisfy DAE order conditions up to `order`. If this is not true, + * these coefficients should not be used to solve models with algebraic conditions (indicated by a + * `Model::Evaluator::tag_` value of 0). + * + */ + bool is_dae_; + + /** + * @brief Whether or not this tableau contains an embedded error estimator method. + * + */ + constexpr bool hasEmbedded() const + { + return static_cast(e_); + } + + /** + * @brief Whether or not this tableau contains coefficients which can be used to construct dense output. + * + */ + constexpr bool hasDenseOutput() const + { + return static_cast(H_); + } + + /** + * @brief Helper function for accessing elements of `A` + * + */ + constexpr RealT getA(size_t row, size_t col) const + { + return A_[row * num_stages_ + col]; + } + + constexpr bool canReuseAsum(size_t stage) const; + constexpr bool canReuseAsumForOut() const; + constexpr std::optional errorEstimatorStage() const; + + static Tableau linImplicitEuler(); + static Tableau rodas5p(); + }; + + private: + /** + * @brief The current step size of the integrator. If integration is ever stopped and resumed, this value will be used for + * the initial step after resuming. + * + */ + RealT step_size_ = 0; + /** + * @brief The step size of the previous step. Used for operations which need to be done on the current step, but step size + * control for the next step has already been performed. + * + */ + RealT prev_step_size_ = 0; + /** + * @brief Whether or not the integrator should attempt to skip Jacobian decomposition on the next step. Controlled by the + * time stepping algorithm in \link integrate() \endlink . Generally, this should only be set if we suspect the Jacobian for the + * last step is a good enough approximation of the Jacobian on the next step. Non-ROW methods need exact Jacobians, + * so this should only be set for ROW methods, and when the step size for the next step is the same as the step size + * as the previous step. + * + */ + bool skip_lu_ = false; + /** + * @brief Whether or not the integrator should attempt to skip the residual function evaluation of the first stage on the + * next step. This should only be used when a step is rejected and the residual function is evaluated at the exact + * same arguments as the previous step. Then \ref RHS_first_stage_ can be re-used rather than re-calculated. + * + */ + bool skip_f_ = false; + /** + * @brief Keeps track of whether or not the integrator currently has valid dense coefficients. + * i.e. they have been computed and haven't been invalidated by taking another step. This can be used to avoid + * re-computing dense coefficients when interpolating states multiple times in one step. + * + */ + bool dense_coefficients_valid_ = false; + + /** + * @brief The tableau of Rosenbrock coefficients currently being used by the integrator. + * + */ + Tableau tab_; + /** + * @brief The model being simulated. + * + */ + GridKit::Model::Evaluator* model_; + /** + * @brief The linear solver to be used during integration in \link timeStep() \endlink. + * + */ + GridKit::LinearAlgebra::LinearSolver& lin_solver_; + /** + * @brief The vector handler to be used for vector operations by the integrator. + * + */ + GridKit::LinearAlgebra::VectorHandler& vector_handler_; + /** + * @brief The `ErrorNorm` to be used by the `StepController` in \link integrate() \endlink. + * + * @note Can be `nullptr` if no `ErrorNorm` is configured, in case the `StepController` does not need an error to be calculated. + * + * @todo Should be removed from the `Rosenbrock` class. Whether or not this is needed is dependent on the `StepController`, so it should be stored there. + */ + const ErrorNorm* err_norm_; + /** + * @brief The memory space where linear algebra operations hsould be done in. + * + */ + GridKit::memory::MemorySpace memspace_; + + /** + * @brief The current simulation time. + * + */ + RealT current_time_; + + /** + * @brief The state from last step. + * + */ + std::unique_ptr y_prev_; + /** + * @brief The state on the current step. + * + */ + std::unique_ptr y_cur_; + /** + * @brief The incoming state for the next step. + * + */ + std::unique_ptr y_new_; + /** + * @brief Used as output for dense output interpolation. + * + */ + std::unique_ptr y_interp_; + + /** + * @brief Configured parameters for the integrator. + * + */ + Parameters params_; + + /** + * @brief Running statistics of the integrator. Reset by \link initializeSimulation() \endlink. + * + */ + Stats stats_; + + public: + Rosenbrock(Tableau&& tab, + GridKit::Model::Evaluator* model, + GridKit::LinearAlgebra::LinearSolver& lin_solver, + GridKit::LinearAlgebra::VectorHandler& vector_handler, + const ErrorNorm* err_norm, + GridKit::memory::MemorySpace memspace = GridKit::memory::HOST); + + [[nodiscard("May fail. Check error code.")]] + int allocate(); + + [[nodiscard("May fail. Check error code.")]] + int initializeSimulation(RealT t0); + + [[nodiscard("May fail. Check error code.")]] + int integrate(const std::vector& out_times, + StepController& step_controller, + Parameters params = {}, + std::optional> out_cb = {}, + std::optional> step_cb = {}); + + private: + /** + * @brief A workspace used by `time_step` for the necessary linear algebra computations. + * + */ + mutable struct + { + /** + * @brief Sum of A-weighted states used for function evaluation in a stage. + * + */ + std::unique_ptr asum_; + /** + * @brief Sum of C-weighted states used in linearization in a stage. + * + */ + std::unique_ptr csum_; + /** + * @brief Right-hand side of linear solve used in a stage. + * + */ + std::unique_ptr RHS_; + /** + * @brief Right-hand side of linear solve used in the first stage. + * Stored separately from \ref RHS_ since it can be re-used by future steps if the current step is rejected. + * + * @see `skip_f_` + * + */ + std::unique_ptr RHS_first_stage_; + /** + * @brief Time Jacobian for non-autonomous models. Currently unused, since models are assumed to be autonomous. + * + */ + std::unique_ptr dFdt_; + /** + * @brief Vector representation of a diagonal mass matrix. + * + */ + std::unique_ptr mass_; + /** + * @brief Estimated error produced by a step in a method with an empbedded error estimator. + * + * @see `Tableau::e` + * + */ + std::unique_ptr err_est_; + + /** + * @brief Whether or not the Jacobian has been factorized. Initial factorization is required, + * but pivots can be re-used for a refactorization later. Re-factorization can introduce + * errors when the Jacobian differs significantly, so this flag should be reset when + * a significant event can happen, such as re-initializing the simulation. + * + * @see \link initializeSimulation() \endlink + * + */ + bool jacobian_analyzed_ = false; + + /** + * @brief Integration stages - each a low-order approximation of the state at the abscissae given by + * the tableau. + * + */ + std::unique_ptr[]> stages_; + } workspace_; + + /** + * @brief Dense output interpolation nodes. Used to generate output states in-between steps. + * + * @see \link calcDenseCoeff() \endlink, \link interpDense() \endlink + * + */ + std::unique_ptr[]> dense_coeff_; + + public: + [[nodiscard("May fail. Check error code.")]] + int timeStep(RealT t0, RealT dt); + + State& errorEstimate() const; + + [[nodiscard("May fail. Check error code.")]] + int calcDenseCoeff(); + + [[nodiscard("May fail. Check error code.")]] + int interpDense(RealT theta); + }; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/GridKit/Solver/Dynamic/RosenbrockTableaus.cpp b/GridKit/Solver/Dynamic/RosenbrockTableaus.cpp new file mode 100644 index 000000000..4fd952bc2 --- /dev/null +++ b/GridKit/Solver/Dynamic/RosenbrockTableaus.cpp @@ -0,0 +1,205 @@ +#include "Rosenbrock.hpp" + +namespace AnalysisManager +{ + namespace NativeDynamicSolver + { + template + Rosenbrock::Tableau Rosenbrock::Tableau::linImplicitEuler() + { + constexpr size_t num_stages = 1; + + Tableau re = { + .num_stages_ = num_stages, + .gamma_ = 1.0, + .alpha_sum_ = std::make_unique_for_overwrite(num_stages), + .gamma_sum_ = std::make_unique_for_overwrite(num_stages), + .m_ = std::make_unique_for_overwrite(num_stages), + .e_ = {}, + .A_ = std::make_unique_for_overwrite(num_stages * num_stages), + .C_ = std::make_unique_for_overwrite(num_stages * num_stages), + .H_ = {}, + .order_ = 1, + .is_krylov_ = true, + .is_w_ = false, + .is_dae_ = true, + }; + + re.alpha_sum_[0] = 0.0; + + re.gamma_sum_[0] = 1.0; + + re.m_[0] = 1.0; + + re.A_[0] = 0.0; + + re.C_[0] = 2.0; + + return re; + } + + /** + * @brief A 5th order Rosenbrock method suitable for solving parabolic PDEs. Created in \cite steinebach2023construction. + * + * + */ + template + Rosenbrock::Tableau Rosenbrock::Tableau::rodas5p() + { + constexpr size_t num_stages = 8; + + Tableau re = { + .num_stages_ = num_stages, + .gamma_ = 0.21193756319429014, + .alpha_sum_ = std::make_unique_for_overwrite(num_stages), + .gamma_sum_ = std::make_unique_for_overwrite(num_stages), + .m_ = std::make_unique_for_overwrite(num_stages), + .e_ = std::make_unique_for_overwrite(num_stages), + .A_ = std::make_unique_for_overwrite(num_stages * num_stages), + .C_ = std::make_unique_for_overwrite(num_stages * num_stages), + .H_ = std::make_unique_for_overwrite(num_stages * 3), + .order_ = 5, + .is_krylov_ = false, + .is_w_ = false, + .is_dae_ = true, + }; + + re.alpha_sum_[0] = 0.0; + re.alpha_sum_[1] = 0.6358126895828704; + re.alpha_sum_[2] = 0.4095798393397535; + re.alpha_sum_[3] = 0.9769306725060716; + re.alpha_sum_[4] = 0.4288403609558664; + re.alpha_sum_[5] = 1.0; + re.alpha_sum_[6] = 1.0; + re.alpha_sum_[7] = 1.0; + + re.gamma_sum_[0] = 0.21193756319429014; + re.gamma_sum_[1] = -0.42387512638858027; + re.gamma_sum_[2] = -0.3384627126235924; + re.gamma_sum_[3] = 1.8046452872882734; + re.gamma_sum_[4] = 2.325825639765069; + re.gamma_sum_[5] = 0.0; + re.gamma_sum_[6] = 0.0; + re.gamma_sum_[7] = 0.0; + + re.m_[0] = -7.502846399306121; + re.m_[1] = 2.561846144803919; + re.m_[2] = -11.627539656261098; + re.m_[3] = -0.18268767659942256; + re.m_[4] = 0.030198172008377946; + re.m_[5] = 1.0; + re.m_[6] = 1.0; + re.m_[7] = 1.0; + + re.e_[0] = 0.0; + re.e_[1] = 0.0; + re.e_[2] = 0.0; + re.e_[3] = 0.0; + re.e_[4] = 0.0; + re.e_[5] = 0.0; + re.e_[6] = 0.0; + re.e_[7] = 1.0; + + re.A_[1 * 8 + 0] = 3.0; + + re.A_[2 * 8 + 0] = 2.849394379747939; + re.A_[2 * 8 + 1] = 0.45842242204463923; + + re.A_[3 * 8 + 0] = -6.954028509809101; + re.A_[3 * 8 + 1] = 2.489845061869568; + re.A_[3 * 8 + 2] = -10.358996098473584; + + re.A_[4 * 8 + 0] = 2.8029986275628964; + re.A_[4 * 8 + 1] = 0.5072464736228206; + re.A_[4 * 8 + 2] = -0.3988312541770524; + re.A_[4 * 8 + 3] = -0.04721187230404641; + + re.A_[5 * 8 + 0] = -7.502846399306121; + re.A_[5 * 8 + 1] = 2.561846144803919; + re.A_[5 * 8 + 2] = -11.627539656261098; + re.A_[5 * 8 + 3] = -0.18268767659942256; + re.A_[5 * 8 + 4] = 0.030198172008377946; + + re.A_[6 * 8 + 0] = -7.502846399306121; + re.A_[6 * 8 + 1] = 2.561846144803919; + re.A_[6 * 8 + 2] = -11.627539656261098; + re.A_[6 * 8 + 3] = -0.18268767659942256; + re.A_[6 * 8 + 4] = 0.030198172008377946; + re.A_[6 * 8 + 5] = 1.0; + + re.A_[7 * 8 + 0] = -7.502846399306121; + re.A_[7 * 8 + 1] = 2.561846144803919; + re.A_[7 * 8 + 2] = -11.627539656261098; + re.A_[7 * 8 + 3] = -0.18268767659942256; + re.A_[7 * 8 + 4] = 0.030198172008377946; + re.A_[7 * 8 + 5] = 1.0; + re.A_[7 * 8 + 6] = 1.0; + + re.C_[1 * 8 + 0] = -14.155112264123755; + + re.C_[2 * 8 + 0] = -17.97296035885952; + re.C_[2 * 8 + 1] = -2.859693295451294; + + re.C_[3 * 8 + 0] = 147.12150275711716; + re.C_[3 * 8 + 1] = -1.41221402718213; + re.C_[3 * 8 + 2] = 71.68940251302358; + + re.C_[4 * 8 + 0] = 165.43517024871676; + re.C_[4 * 8 + 1] = -0.4592823456491126; + re.C_[4 * 8 + 2] = 42.90938336958603; + re.C_[4 * 8 + 3] = -5.961986721573306; + + re.C_[5 * 8 + 0] = 24.854864614690072; + re.C_[5 * 8 + 1] = -3.0009227002832186; + re.C_[5 * 8 + 2] = 47.4931110020768; + re.C_[5 * 8 + 3] = 5.5814197821558125; + re.C_[5 * 8 + 4] = -0.6610691825249471; + + re.C_[6 * 8 + 0] = 30.91273214028599; + re.C_[6 * 8 + 1] = -3.1208243349937974; + re.C_[6 * 8 + 2] = 77.79954646070892; + re.C_[6 * 8 + 3] = 34.28646028294783; + re.C_[6 * 8 + 4] = -19.097331116725623; + re.C_[6 * 8 + 5] = -28.087943162872662; + + re.C_[7 * 8 + 0] = 37.80277123390563; + re.C_[7 * 8 + 1] = -3.2571969029072276; + re.C_[7 * 8 + 2] = 112.26918849496327; + re.C_[7 * 8 + 3] = 66.9347231244047; + re.C_[7 * 8 + 4] = -40.06618937091002; + re.C_[7 * 8 + 5] = -54.66780262877968; + re.C_[7 * 8 + 6] = -9.48861652309627; + + re.H_[0 * 8 + 0] = 25.948786856663858; + re.H_[0 * 8 + 1] = -2.5579724845846235; + re.H_[0 * 8 + 2] = 10.433815404888879; + re.H_[0 * 8 + 3] = -2.3679251022685204; + re.H_[0 * 8 + 4] = 0.524948541321073; + re.H_[0 * 8 + 5] = 1.1241088310450404; + re.H_[0 * 8 + 6] = 0.4272876194431874; + re.H_[0 * 8 + 7] = -0.17202221070155493; + + re.H_[1 * 8 + 0] = -9.91568850695171; + re.H_[1 * 8 + 1] = -0.9689944594115154; + re.H_[1 * 8 + 2] = 3.0438037242978453; + re.H_[1 * 8 + 3] = -24.495224566215796; + re.H_[1 * 8 + 4] = 20.176138334709044; + re.H_[1 * 8 + 5] = 15.98066361424651; + re.H_[1 * 8 + 6] = -6.789040303419874; + re.H_[1 * 8 + 7] = -6.710236069923372; + + re.H_[2 * 8 + 0] = 11.419903575922262; + re.H_[2 * 8 + 1] = 2.8879645146136994; + re.H_[2 * 8 + 2] = 72.92137995996029; + re.H_[2 * 8 + 3] = 80.12511834622643; + re.H_[2 * 8 + 4] = -52.072871366152654; + re.H_[2 * 8 + 5] = -59.78993625266729; + re.H_[2 * 8 + 6] = -0.15582684282751913; + re.H_[2 * 8 + 7] = 4.883087185713722; + + return re; + } + + template class Rosenbrock; + } // namespace NativeDynamicSolver +} // namespace AnalysisManager diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index 0b0368213..38a7a9eba 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -9,6 +9,8 @@ if(${DOXYGEN_FOUND}) set(DOXYGEN_SOURCE_BROWSER YES) set(DOXYGEN_INTERACTIVE_SVG YES) set(DOXYGEN_DISTRIBUTE_GROUP_DOC YES) + set(DOXYGEN_USE_MATHJAX YES) + set(DOXYGEN_CITE_BIB_FILES ${CMAKE_SOURCE_DIR}/docs/citations.bib) doxygen_add_docs(GridKitDocs ${CMAKE_SOURCE_DIR}/GridKit ${CMAKE_SOURCE_DIR}/README.md) endif() diff --git a/docs/citations.bib b/docs/citations.bib new file mode 100644 index 000000000..739cd0a19 --- /dev/null +++ b/docs/citations.bib @@ -0,0 +1,10 @@ +@article{steinebach2023construction, + title={Construction of Rosenbrock--Wanner method Rodas5P and numerical benchmarks within the Julia Differential Equations package}, + author={Steinebach, Gerd}, + journal={BIT Numerical Mathematics}, + volume={63}, + number={2}, + pages={27}, + year={2023}, + publisher={Springer} +} diff --git a/tests/UnitTests/Solver/CMakeLists.txt b/tests/UnitTests/Solver/CMakeLists.txt index b6a9ff80a..5acde4132 100644 --- a/tests/UnitTests/Solver/CMakeLists.txt +++ b/tests/UnitTests/Solver/CMakeLists.txt @@ -3,6 +3,4 @@ # - Slaven Peles #]] -if(TARGET SUNDIALS::idas) - add_subdirectory(Dynamic) -endif() +add_subdirectory(Dynamic) diff --git a/tests/UnitTests/Solver/Dynamic/CMakeLists.txt b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt index 93b92a471..35d1f70e2 100644 --- a/tests/UnitTests/Solver/Dynamic/CMakeLists.txt +++ b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt @@ -1,9 +1,24 @@ # -add_executable(test_ida runIdaTests.cpp) -target_link_libraries( - test_ida GridKit::solvers_dyn GridKit::testing) +if(TARGET SUNDIALS::idas) + add_executable(test_ida runIdaTests.cpp) + target_link_libraries(test_ida GridKit::solvers_dyn GridKit::testing) -add_test(NAME IDATest COMMAND $) + add_test(NAME IDATest COMMAND $) -install(TARGETS test_ida RUNTIME DESTINATION bin) + install(TARGETS test_ida RUNTIME DESTINATION bin) +endif() + +if(GRIDKIT_ENABLE_RESOLVE) + add_executable(test_rosenbrock runRosenbrockTests.cpp) + target_link_libraries( + test_rosenbrock + GridKit::rosenbrock + GridKit::testing + GridKit::fixed_step + GridKit::resolve_system_solver + ReSolve::ReSolve) + add_test(NAME RosenbrockTest COMMAND $) + + install(TARGETS test_rosenbrock RUNTIME DESTINATION bin) +endif() diff --git a/tests/UnitTests/Solver/Dynamic/RosenbrockTests.hpp b/tests/UnitTests/Solver/Dynamic/RosenbrockTests.hpp new file mode 100644 index 000000000..49ac22649 --- /dev/null +++ b/tests/UnitTests/Solver/Dynamic/RosenbrockTests.hpp @@ -0,0 +1,482 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace Model + { + /** + * @brief A test DAE (referred to as the "trigonometric DAE") that is useful for evaluating order + * of DAE integrators. The system is 2-dimensional, with one differential variable and one algebraic variable. + * The DAE is index 1 and the derivatives of the model are non-vanishing. + * + * The valid simulation time interval is \f([0.5,2]\f) and the model is initialized at \f(t = 0.5\f). + * + */ + template + class TrigonometricDaeEvaluator : public Model::Evaluator + { + public: + using RealT = typename Model::Evaluator::RealT; + + constexpr static size_t SIZE = 2; + + TrigonometricDaeEvaluator() + { + } + + int allocate() override + { + constexpr size_t NNZ = SIZE * SIZE; + + IdxT* row_ptrs = new IdxT[SIZE + 1]; + IdxT* cols = new IdxT[NNZ]; + RealT* vals = new RealT[NNZ]; + + for (size_t i = 0; i < SIZE + 1; i++) + { + row_ptrs[i] = static_cast(i * SIZE); + } + + for (size_t i = 0; i < SIZE; i++) + { + for (size_t j = 0; j < SIZE; j++) + { + cols[i * SIZE + j] = static_cast(j); + } + } + + csr_jac_ = std::make_unique>(SIZE, SIZE, NNZ, &row_ptrs, &cols, &vals); + + return 0; + } + + int initialize() override + { + y_ = {sinh(0.5), tanh(0.5)}; + yp_ = {0.0, 0.0}; + + tag_ = {true, false}; + + f_ = {0.0, 0.0}; + g_ = {0.0, 0.0}; + return 0; + } + + IdxT size() override + { + return SIZE; + } + + IdxT nnz() override + { + return SIZE * SIZE; + } + + bool hasJacobian() override + { + return true; + } + + IdxT sizeQuadrature() override + { + return 0; + } + + IdxT sizeParams() override + { + return 0; + } + + int setAbsoluteTolerance([[maybe_unused]] RealT rel_tol) override + { + return 0; + } + + std::vector& absoluteTolerance() override + { + return abs_tol_; + } + + const std::vector& absoluteTolerance() const override + { + return abs_tol_; + } + + int tagDifferentiable() override + { + return 0; + } + + int evaluateResidual() override + { + double y02 = y_[0] * y_[0]; + double y12 = y_[1] * y_[1]; + + f_ = {y02 / (y_[1] * sqrt(std::pow(y_[0] / y_[1], 2) - 1)), // + y12 + 1 / (1 + y02) - (y02 / y12 - y02)}; + + return 0; + } + + int evaluateJacobian() override + { + RealT* vals = csr_jac_->getValues(); + + double y1 = y_[0]; + double y2 = y_[1]; + + double y12 = y1 * y1; + double y13 = y12 * y1; + double y22 = y2 * y2; + double y23 = y22 * y2; + double y24 = y22 * y22; + + double tmp = std::pow(y12 / y22 - 1, 1.5); + double tmp2 = std::pow(y12 + 1, 2); + + vals[0] = alpha_ + -(-y13 + 2 * y1 * y22) / (y23 * tmp); + vals[1] = y12 / (y22 * tmp); + vals[2] = -2 * y1 * (1 / y22 - 1) - (2 * y1) / tmp2; + vals[3] = (2 * (y12 + y24)) / y23; + + return 0; + } + + int evaluateIntegrand() override + { + return 0; + } + + int initializeAdjoint() override + { + return 0; + } + + int evaluateAdjointResidual() override + { + return 0; + } + + int evaluateAdjointIntegrand() override + { + return 0; + } + + void updateTime([[maybe_unused]] RealT t, RealT a) override + { + alpha_ = a; + } + + std::vector& y() override + { + return y_; + } + + const std::vector& y() const override + { + return y_; + } + + std::vector& yp() override + { + return yp_; + } + + const std::vector& yp() const override + { + return yp_; + } + + std::vector& tag() override + { + return tag_; + } + + const std::vector& tag() const override + { + return tag_; + } + + std::vector& yB() override + { + return yB_; + } + + const std::vector& yB() const override + { + return yB_; + } + + std::vector& ypB() override + { + return ypB_; + } + + const std::vector& ypB() const override + { + return ypB_; + } + + std::vector& param() override + { + return param_; + } + + const std::vector& param() const override + { + return param_; + } + + std::vector& param_up() override + { + return param_up_; + } + + const std::vector& param_up() const override + { + return param_up_; + } + + std::vector& param_lo() override + { + return param_lo_; + } + + const std::vector& param_lo() const override + { + return param_lo_; + } + + std::vector& getResidual() override + { + return f_; + } + + const std::vector& getResidual() const override + { + return f_; + } + + std::vector& getIntegrand() override + { + return g_; + } + + const std::vector& getIntegrand() const override + { + return g_; + } + + std::vector& getAdjointResidual() override + { + return fB_; + } + + const std::vector& getAdjointResidual() const override + { + return fB_; + } + + std::vector& getAdjointIntegrand() override + { + return gB_; + } + + const std::vector& getAdjointIntegrand() const override + { + return gB_; + } + + IdxT getIDcomponent() + { + return 0; + } + + GridKit::LinearAlgebra::CsrMatrix* getCsrJacobian() const override + { + return csr_jac_.get(); + } + + protected: + std::vector y_; + std::vector yp_; + std::vector tag_; + std::vector f_; + std::vector g_; + + std::vector yB_; + std::vector ypB_; + std::vector fB_; + std::vector gB_; + + std::vector abs_tol_; + + std::unique_ptr> csr_jac_; + + RealT alpha_; + + std::vector param_; + std::vector param_up_; + std::vector param_lo_; + }; + } // namespace Model + + namespace Testing + { + template + class RosenbrockTests + { + using Rosenbrock = AnalysisManager::NativeDynamicSolver::Rosenbrock; + using RealT = typename GridKit::ScalarTraits::RealT; + + public: + /** + * @brief Test a Rosenbrock tableau by verifying its order empirically. Empirical order is calculated + * by running a convergence test on the problem modeled in \ref Model::TrigonometricDaeEvaluator. 21 + * simulations are run on the problem with a fixed-step step size controller with a number of steps + * logarithmically distributed in \f(\left[10^a, 10^b\right]\f) where \f(a\f) is `step_exponent_lower` + * and \f(b\f) is `step_exponent_upper`. The average slope between pairs of simulations is taken to be + * the empirical order. The test succeeds if the empirical order is at least 85% of the expected theoretical order. + * + * @param tab The tableau to test + * @param step_exponent_lower The exponent describing the smallest number of steps to take during a simulation. + * @param step_exponent_upper The exponent describing the largest number of steps to take during a simulation. + */ + TestOutcome test_order(Rosenbrock::Tableau&& tab, double step_exponent_lower, double step_exponent_upper) + { + TestStatus success = true; + + // Tableaus keep track of their theoretical order. We will attempt to match the empirical order to this. + uint8_t expected_order = tab.order_; + + // Setup the model, linear solver, and integrator + Model::TrigonometricDaeEvaluator model; + model.allocate(); + model.initialize(); + + ReSolve::LinAlgWorkspaceCpu linear_workspace; + ReSolve::SystemSolver resolve_solver(&linear_workspace, "klu", "klu", "klu"); + GridKit::LinearAlgebra::VectorHandler vec_handler; + GridKit::LinearAlgebra::ResolveSystemSolver lin_solver(resolve_solver); + + resolve_solver.initialize(); + + Rosenbrock integrator(std::move(tab), &model, lin_solver, vec_handler, nullptr); + if (integrator.allocate()) + { + success = false; + return success.report(__func__); + } + + AnalysisManager::NativeDynamicSolver::FixedStep step_controller; + + // The number of simulations to run to calculate empirical order. Must be at least 2, + // since empirical order is calculated pairwise between simulations. + size_t num_samples = 21; + + // Output vectors to keep track of data used to calculate empirical order. + // Each simulation will record its step size and final error. + std::vector step_sizes; + std::vector errors; + + // An output callback to populate step_sizes and errors. Each simulation will call this callback + // at the end of the simulation (final_time). It will then calculate the final error by comparing the + // solution of the simulation to the true answer. + double final_time = 2.0; + std::vector out_times = {final_time}; + auto out_cb = [&]([[maybe_unused]] double t) + { + double error = 0.0; + double sol_norm = 0.0; + + // The final solution of the simulation + const std::vector& state = model.y(); + + // The difference from the simulated solution to the true solution + error += std::pow(state[0] - sinh(final_time), 2); + error += std::pow(state[1] - tanh(final_time), 2); + + sol_norm += std::pow(sinh(final_time), 2); + sol_norm += std::pow(tanh(final_time), 2); + + // Error relative to the true solution + errors.push_back(std::sqrt(error) / std::sqrt(sol_norm)); + }; + + // Perform all of the simulations and populate step_sizes and errors using out_cb + for (size_t i = 0; i < num_samples; i++) + { + // Logarithmically distribute step_size based on step_exponent_lower and step_exponent_upper. + // Round it to ensure num_steps is integral (do not invoke the dense output, which can add additional errors). + double step_size = std::pow(10, step_exponent_lower + static_cast(i) * (step_exponent_upper - step_exponent_lower) / static_cast(num_samples - 1)); + double num_steps = round((final_time - 0.5) / step_size); + step_size = (final_time - 0.5) / num_steps; + step_sizes.push_back(step_size); + + model.initialize(); + if (integrator.initializeSimulation(0.5)) + { + success = false; + return success.report(__func__); + } + + typename Rosenbrock::Parameters params; + params.starting_step_ = step_size; + params.max_steps_ = static_cast(ceil((final_time - 0.5) / step_size)) + 10; + if (integrator.integrate(out_times, step_controller, params, out_cb)) + { + success = false; + return success.report(__func__); + } + } + + // Print output data + std::cerr << "Step sizes\n"; + for (double step_size : step_sizes) + { + std::cerr << std::scientific << std::setprecision(20) << step_size << ' '; + } + + std::cerr << "\n\nErrors\n"; + for (double error : errors) + { + std::cerr << std::scientific << std::setprecision(20) << error << ' '; + } + std::cerr << "\n\n"; + + // Calculate empirical order. Each pairwise order is calculated, then averaged. + std::vector pairwise_orders; + double empirical_order = 0.0; + for (size_t i = 1; i < num_samples; i++) + { + pairwise_orders.push_back((log(errors[i]) - log(errors[i - 1])) / (log(step_sizes[i]) - log(step_sizes[i - 1]))); + empirical_order += pairwise_orders.back(); + } + empirical_order /= static_cast(num_samples - 1); + + // Print test result - observed empirical order and the expected order. + std::cerr << "Empirical order: " << std::fixed << std::setprecision(5) << empirical_order << " v.s. expected " << static_cast(expected_order) << '\n'; + success *= empirical_order > expected_order * 0.85; + + return success.report(__func__); + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/Solver/Dynamic/runRosenbrockTests.cpp b/tests/UnitTests/Solver/Dynamic/runRosenbrockTests.cpp new file mode 100644 index 000000000..c311c2e8e --- /dev/null +++ b/tests/UnitTests/Solver/Dynamic/runRosenbrockTests.cpp @@ -0,0 +1,18 @@ +#include + +#include "RosenbrockTests.hpp" + +int main() +{ + using namespace GridKit; + using namespace GridKit::Testing; + using namespace AnalysisManager::NativeDynamicSolver; + + GridKit::Testing::TestingResults result; + GridKit::Testing::RosenbrockTests test; + + result += test.test_order(Rosenbrock::Tableau::linImplicitEuler(), -1.0, -5.0); + result += test.test_order(Rosenbrock::Tableau::rodas5p(), -0.5, -2.5); + + return result.summary(); +} From bdec5dc2413bde84caffc7711386f14f75a625ba Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Fri, 17 Jul 2026 13:55:05 -0400 Subject: [PATCH 6/6] Update resolve dependency requirement Co-authored-by: pelesh --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index e80b68405..76d37a7bc 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,7 @@ You should have all of the following installed before installing GridKit™ that the use of Enzyme is experimental, and some versions of it have been found to break GridKit code. - [LLVM](https://github.com/llvm/llvm-project) >= 15.x. GridKit is currently tested with LLVM 16. - - [ReSolve](https://github.com/ORNL/ReSolve) == commit `a93fb6571e3542b2fa779983fef0186a034c53c9` (optional) + - [ReSolve](https://github.com/ORNL/ReSolve) >= 0.99.3 (optional) - [CMake](https://cmake.org/) >= 3.13 - C++ 20 compliant compiler