SUNDIALS Upgrades

CHANGELOG.md

@@ -73,7 +73,7 @@
 - Add new `LinearSolver` interface for linear solvers.
 - Added new `Rosenbrock` integrator.
 - Clarified naming conventions for macros.
-- Added `dt_fixed`, `rel_tol`, `abs_tol`, and `max_steps` options to phasor dynamics solver JSON files and renamed the `dt` option to `dt_monitor`.
+- Added `dt_fixed`, `rel_tol`, `abs_tol`, `max_steps`, and `max_order` options to phasor dynamics solver JSON files and renamed the `dt` option to `dt_monitor`.
 - Added EMT model and operator documentation.
 - Added `REGCA` converter model implementation for PhasorDynamics.
 - Remove unnecessary data copying while evaluating `PowerElectronics` models, speeding up large simulations by up to 3x

GridKit/Solver/Dynamic/Ida.cpp

@@ -96,7 +96,7 @@ namespace AnalysisManager
       retval = IDASetId(solver_, tag_);
       checkOutput(retval, "IDASetId");
 
-      setIDAOptions(solver_, time_step_, rel_tol_, abs_tol_override_, max_steps_, suppress_alg_);
+      setIDAOptions(solver_, time_step_, rel_tol_, abs_tol_override_, max_steps_, max_order_, suppress_alg_);
 
       // Set up linear solver
       return this->configureLinearSolver();
@@ -551,6 +551,7 @@ namespace AnalysisManager
                     backward_rel_tol_,
                     backward_abs_tol_override_,
                     backward_max_steps_,
+                    backward_max_order_,
                     backward_suppress_alg_);
 
       retval = IDASetUserDataB(solver_, backwardID_, model_);
@@ -1250,6 +1251,32 @@ namespace AnalysisManager
       backward_max_steps_ = max_steps;
     }
 
+    /**
+     * @brief Set the maximum integration method order
+     *
+     * @param max_order The maximum integration method order
+     * @tparam ScalarT Scalar data type
+     * @tparam IdxT Index data type
+     */
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::setMaxOrder(int max_order)
+    {
+      max_order_ = max_order;
+    }
+
+    /**
+     * @brief Set the maximum integration method order for the backward simulation
+     *
+     * @param max_order The maximum integration method order
+     * @tparam ScalarT Scalar data type
+     * @tparam IdxT Index data type
+     */
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::setBackwardMaxOrder(int max_order)
+    {
+      backward_max_order_ = max_order;
+    }
+
     /**
      * @brief A helper function to set common IDA options
      *
@@ -1260,6 +1287,7 @@ namespace AnalysisManager
      *        absolute tolerance for the nonlinear solver rather than the
      *        model's default absolute tolerance
      * @param max_steps The maximum number of steps
+     * @param max_order The maximum integration method order
      * @param suppress_alg If true, algebraic variables are excluded from IDA's
      *        local error test
      * @tparam ScalarT Scalar data type
@@ -1271,6 +1299,7 @@ namespace AnalysisManager
                                            ScalarT rel_tol,
                                            ScalarT abs_tol_override,
                                            IdxT    max_steps,
+                                           int     max_order,
                                            bool    suppress_alg)
     {
       int retval = 0;
@@ -1280,6 +1309,8 @@ namespace AnalysisManager
       checkOutput(retval, "IDASetMaxStep");
       retval = IDASetMaxNumSteps(mem, static_cast<long int>(max_steps));
       checkOutput(retval, "IDASetMaxNumSteps");
+      retval = IDASetMaxOrd(mem, time_step == 0 ? max_order : std::min(max_order, 2));
+      checkOutput(retval, "IDASetMaxOrd");
       retval = IDASetSuppressAlg(mem, suppress_alg ? SUNTRUE : SUNFALSE);
       checkOutput(retval, "IDASetSuppressAlg");
 
@@ -1289,11 +1320,6 @@ namespace AnalysisManager
       }
       else
       {
-        /* Since the starting procedure is first order, the maximum global order
-         * of convergence is two */
-        retval = IDASetMaxOrd(mem, 2);
-        checkOutput(retval, "IDASetMaxOrd");
-
         /* Enable more nonlinear iterations because a failed nonlinear solve
          * causes a failed integration with fixed steps */
         static constexpr int FIXED_STEP_NONLIN_ITRS = 16;
@@ -1302,10 +1328,7 @@ namespace AnalysisManager
 
         // Set a large tolerance so the error test will never fail
         static constexpr RealT FIXED_STEP_TOL_FAC = 1e100;
-        setTolerance(mem,
-                     FIXED_STEP_TOL_FAC * rel_tol,
-                     FIXED_STEP_TOL_FAC * abs_tol_override,
-                     FIXED_STEP_TOL_FAC);
+        setTolerance(mem, rel_tol, abs_tol_override, FIXED_STEP_TOL_FAC);
 
         /* We want the nonlinear solver tolerance to be ~rel_tol, but the with
          * the large tolerances set above, we need to choose this tolerance to
@@ -1329,22 +1352,21 @@ namespace AnalysisManager
      * @param abs_tol_override If positive, this value will be used as the
      *        absolute tolerance rather than the model's default absolute
      *        tolerance
-     * @param abs_tol_fac A factor to apply to the absolute tolerance if not
-     *        overridden
+     * @param tol_fac A factor to apply to the relative and absolute tolerances
      * @tparam ScalarT Scalar data type
      * @tparam IdxT Index data type
      */
     template <class ScalarT, typename IdxT>
     void Ida<ScalarT, IdxT>::setTolerance(void*   mem,
                                           ScalarT rel_tol,
                                           ScalarT abs_tol_override,
-                                          ScalarT abs_tol_fac)
+                                          ScalarT tol_fac)
     {
       int retval = 0;
 
       if (abs_tol_override > 0)
       {
-        retval = IDASStolerances(mem, rel_tol, abs_tol_override);
+        retval = IDASStolerances(mem, tol_fac * rel_tol, tol_fac * abs_tol_override);
         checkOutput(retval, "IDASStolerances");
         return;
       }
@@ -1353,9 +1375,9 @@ namespace AnalysisManager
       checkAllocation((void*) abs_tol_vec, "N_VClone");
       model_->setAbsoluteTolerance(rel_tol);
       copyVec(model_->absoluteTolerance(), abs_tol_vec);
-      N_VScale(abs_tol_fac, abs_tol_vec, abs_tol_vec);
+      N_VScale(tol_fac, abs_tol_vec, abs_tol_vec);
 
-      retval = IDASVtolerances(mem, rel_tol, abs_tol_vec);
+      retval = IDASVtolerances(mem, tol_fac * rel_tol, abs_tol_vec);
       checkOutput(retval, "IDASVtolerances");
 
       N_VDestroy(abs_tol_vec);

GridKit/Solver/Dynamic/Ida.hpp

@@ -142,6 +142,8 @@ namespace AnalysisManager
       void setConsistentICType(IdaConsistentICType consistent_ic_type);
       void setMaxSteps(IdxT maxSteps) override;
       void setBackwardMaxSteps(IdxT maxSteps);
+      void setMaxOrder(int max_order);
+      void setBackwardMaxOrder(int max_order);
 
       IdaStats getStats() const;
 
@@ -220,13 +222,15 @@ namespace AnalysisManager
       RealT               rel_tol_{DEFAULT_REL_TOL};
       RealT               abs_tol_override_{};
       IdxT                max_steps_{};
+      int                 max_order_{5};
       bool                suppress_alg_{false};
       IdaConsistentICType consistent_ic_type_{IdaConsistentICType::YA_YDP};
 
       RealT backward_time_step_{};
       RealT backward_rel_tol_{DEFAULT_REL_TOL};
       RealT backward_abs_tol_override_{};
       IdxT  backward_max_steps_{};
+      int   backward_max_order_{5};
       bool  backward_suppress_alg_{false};
 
       RealT quadrature_rel_tol_{0.1 * DEFAULT_REL_TOL};
@@ -250,11 +254,12 @@ namespace AnalysisManager
                          ScalarT rel_tol,
                          ScalarT abs_tol_override,
                          IdxT    max_steps,
+                         int     max_order,
                          bool    suppress_alg);
       void setTolerance(void*   mem,
                         ScalarT rel_tol,
                         ScalarT abs_tol_override,
-                        ScalarT abs_tol_fac = 1);
+                        ScalarT tol_fac = 1);
       void setQuadratureTolerance(void*   mem,
                                   ScalarT rel_tol,
                                   ScalarT abs_tol_override);

INSTALL.md

@@ -27,7 +27,7 @@ MSVC toolchain) are possible but not regularly tested.
 |---|---|---|
 | CMake | >= 3.13 | |
 | C++ compiler | C++20 | Clang or GCC |
-| SUNDIALS | `develop` branch | Optional; disabled by default |
+| SUNDIALS | >= 7.8.0 | Optional; disabled by default |
 | SuiteSparse (KLU) | >= 7.x | Optional; needed for sparse solvers in SUNDIALS |
 | Ipopt | >= 3.14 | Optional; disabled by default |
 | HSL | >= 2015 | Optional; required by Ipopt for efficient linear solvers |
@@ -353,7 +353,7 @@ spack develop --path=$(pwd) gridkit@develop
 spack compiler find
 
 # Add GridKit with desired variants, then build
-spack add gridkit+sundials+ipopt+klu ^sundials@develop
+spack add gridkit+sundials+ipopt+klu ^sundials@7.8.0
 spack concretize -f
 spack install
 spack env deactivate

README.md

@@ -21,7 +21,7 @@ Before installing GridKit™ make sure you have all needed dependencies.
 ### Dependencies
 You should have all of the following installed before installing GridKit™
 - A version of
-    - [SUNDIALS](https://github.com/LLNL/sundials) `develop` branch (optional)
+    - [SUNDIALS](https://github.com/LLNL/sundials) >= 7.8.0 (optional)
         - To support sparse linear algebra, SUNDIALS must also be built with [KLU support](https://sundials.readthedocs.io/en/latest/sundials/Install_link.html#cmakeoption-ENABLE_KLU). You most likely want this.
     - [Ipopt](https://github.com/coin-or/Ipopt) >= 3.x (optional)
     - [Enzyme](https://github.com/EnzymeAD/Enzyme) >=0.0.206 (optional). Note

application/PhasorDynamics/AnalysisUtilities.hpp

@@ -64,6 +64,8 @@ namespace GridKit
       double                                         dt_fixed;
       /// maximum number of solver time steps, or 0 for the IDA default
       std::size_t                                    max_steps;
+      /// maximum IDA integration method order
+      int                                            max_order;
       /// IDA consistent initial condition calculation type
       AnalysisManager::Sundials::IdaConsistentICType consistent_ic_type;
       /// set of system events
@@ -103,6 +105,7 @@ namespace GridKit
       c.abs_tol            = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL);
       c.dt_fixed           = j.value("dt_fixed", 0.0);
       c.max_steps          = j.value("max_steps", std::size_t{0});
+      c.max_order          = j.value("max_order", 5);
       c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP;
       if (j.contains("consistent_ic_type"))
       {

application/PhasorDynamics/ContingencyAnalysis.cpp

@@ -1,4 +1,5 @@
 #include <chrono>
+#include <exception>
 #include <filesystem>
 #include <future>
 
@@ -38,6 +39,7 @@ TestStatus runStudy(StudyData study_data)
   ida.setTolerance(study_data.rel_tol, study_data.abs_tol);
   ida.setFixedStep(study_data.dt_fixed);
   ida.setMaxSteps(study_data.max_steps);
+  ida.setMaxOrder(study_data.max_order);
   ida.setConsistentICType(study_data.consistent_ic_type);
   ida.configureSimulation();
 
@@ -148,7 +150,7 @@ void runStudyOpenMP(const StudyData& study_data, std::vector<TestStatus>& stat_v
 }
 #endif
 
-int main(int argc, const char* argv[])
+int runApplication(int argc, const char* argv[])
 {
   // Study file
   checkCommandLine(argc, "ContingencyAnalysis");
@@ -187,3 +189,17 @@ int main(int argc, const char* argv[])
 
   return status.get();
 }
+
+int main(int argc, const char* argv[])
+{
+  try
+  {
+    return runApplication(argc, argv);
+  }
+  catch (const std::exception& error)
+  {
+    Log::error() << "ContingencyAnalysis failed: " << error.what() << '\n';
+  }
+
+  return 1;
+}

application/PhasorDynamics/DynamicSimulation.cpp

@@ -1,3 +1,4 @@
+#include <exception>
 #include <filesystem>
 #include <fstream>
 
@@ -17,7 +18,7 @@ using scalar_type = double;
 using real_type   = double;
 using index_type  = size_t;
 
-int main(int argc, const char* argv[])
+int runApplication(int argc, const char* argv[])
 {
   // Study file
   checkCommandLine(argc, "DynamicSimulation");
@@ -32,6 +33,7 @@ int main(int argc, const char* argv[])
   ida.setTolerance(study.rel_tol, study.abs_tol);
   ida.setFixedStep(study.dt_fixed);
   ida.setMaxSteps(study.max_steps);
+  ida.setMaxOrder(study.max_order);
   ida.setConsistentICType(study.consistent_ic_type);
   ida.configureSimulation();
 
@@ -80,3 +82,17 @@ int main(int argc, const char* argv[])
 
   return status.get();
 }
+
+int main(int argc, const char* argv[])
+{
+  try
+  {
+    return runApplication(argc, argv);
+  }
+  catch (const std::exception& error)
+  {
+    Log::error() << "DynamicSimulation failed: " << error.what() << '\n';
+  }
+
+  return 1;
+}

application/PhasorDynamics/README.md

@@ -11,6 +11,7 @@
   `abs_tol`            | Absolute solver tolerance override (default: 1.0e-9)
   `dt_fixed`           | Fixed solver time step size, or 0 for adaptive stepping (default: 0)
   `max_steps`          | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0)
+  `max_order`          | Maximum IDA integration method order from 1 to 5 (default: 5; fixed stepping is capped at 2)
   `consistent_ic_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: "ya_ydp")
   `events`             | An array of event groups (see [Events](#events) below)
   `output_file`        | Path to output (CSV) file (optional)

buildsystem/README.md

@@ -11,7 +11,7 @@ spack env activate -p gridkit
 spack repo add buildsystem/spack_repo/gridkit
 spack develop --path=$(pwd) gridkit@develop
 spack compiler find
-spack add gridkit+enzyme+ipopt+klu+sundials ^sundials@develop
+spack add gridkit+enzyme+ipopt+klu+sundials ^sundials@7.8.0
 spack concretize -f
 spack install
 spack env deactivate

buildsystem/spack_repo/gridkit/packages/gridkit/package.py

@@ -33,8 +33,8 @@ class Gridkit(CMakePackage):
     depends_on("ipopt", when="+ipopt")
     depends_on("resolve@gridkit-pinned+klu", when="+resolve+klu")
     depends_on("resolve@gridkit-pinned~klu", when="+resolve~klu")
-    depends_on("sundials@develop+klu~mpi", when="+sundials+klu")
-    depends_on("sundials@develop~klu~mpi", when="+sundials~klu")
+    depends_on("sundials@7.8.0+klu~mpi", when="+sundials+klu")
+    depends_on("sundials@7.8.0~klu~mpi", when="+sundials~klu")
 
     def cmake_args(self):
         args = []

buildsystem/spack_repo/gridkit/packages/sundials/package.py

@@ -0,0 +1,7 @@
+from spack_repo.builtin.packages.sundials.package import Sundials as BuiltinSundials
+
+from spack.package import *
+
+class Sundials(BuiltinSundials):
+    version("7.9.0", tag="v7.9.0", commit="312fc0f3684f27209ca9dc9249d194436eb41a7a")
+    version("7.8.0", tag="v7.8.0", commit="aedc088437064dd55b35c000145f7f5db6ee49e3")

tests/UnitTests/Solver/Dynamic/IdaTests.hpp

@@ -526,15 +526,20 @@ namespace GridKit
       TestOutcome fixedStep()
       {
         const unsigned n_steps = 32;
+        const double   tol     = 1.0e-6;
         TestStatus     success = true;
 
         Model::NullEvaluator<ScalarT, IdxT> model;
 
         Ida<double, size_t> ida(&model);
         ida.setFixedStep(1.0 / n_steps);
-        ida.setTolerance(1.0e-6);
+        ida.setTolerance(tol);
         ida.configureSimulation();
 
+        // Fixed-step error-test scaling must not affect the tolerance used by
+        // the model to construct its absolute-tolerance vector.
+        success *= (model.absoluteTolerance().getData()[0] == tol);
+
         ida.initializeSimulation(0.0, false);
         ida.runSimulation(1.0);
         auto stats = ida.getStats();