Update sundials to v7.3.0

install_KLU_Sundials.py

@@ -19,10 +19,10 @@ def build_solvers():
     SUITESPARSE_CHECKSUM = (
         "ce39b28d4038a09c14f21e02c664401be73c0cb96a9198418d6a98a7db73a259"
     )
-    SUNDIALS_VERSION = "6.5.0"
+    SUNDIALS_VERSION = "7.3.0"
     SUNDIALS_URL = f"https://github.com/LLNL/sundials/releases/download/v{SUNDIALS_VERSION}/sundials-{SUNDIALS_VERSION}.tar.gz"
     SUNDIALS_CHECKSUM = (
-        "4e0b998dff292a2617e179609b539b511eb80836f5faacf800e688a886288502"
+        "697b7b0dbc229f149e39b293d1ab03d321d61adb6733ffb78c0ddbffaf73d839"
     )
     DEFAULT_INSTALL_DIR = str(pathlib.Path(__file__).parent.resolve() / ".idaklu")
 

install_sundials.sh

@@ -28,7 +28,7 @@ function extract {
 }
 
 SUITESPARSE_VERSION=7.8.3
-SUNDIALS_VERSION=6.5.0
+SUNDIALS_VERSION=7.3.0
 
 SUITESPARSE_ROOT_ADDR=https://github.com/DrTimothyAldenDavis/SuiteSparse/archive
 SUNDIALS_ROOT_ADDR=https://github.com/LLNL/sundials/releases/download/v$SUNDIALS_VERSION

src/pybammsolvers/idaklu_source/Expressions/Base/Expression.hpp

@@ -21,8 +21,8 @@ class Expression {
    * @brief Evaluation operator (supplying data references)
    */
   virtual void operator()(
-    const std::vector<realtype*>& inputs,
-    const std::vector<realtype*>& results) = 0;
+    const std::vector<sunrealtype*>& inputs,
+    const std::vector<sunrealtype*>& results) = 0;
 
   /**
    * @brief The maximum number of elements returned by the k'th output

src/pybammsolvers/idaklu_source/Expressions/Base/ExpressionSet.hpp

@@ -69,16 +69,16 @@ class ExpressionSet
 
   std::vector<int64_t> jac_times_cjmass_rowvals;  // cppcheck-suppress unusedStructMember
   std::vector<int64_t> jac_times_cjmass_colptrs;  // cppcheck-suppress unusedStructMember
-  std::vector<realtype> inputs;  // cppcheck-suppress unusedStructMember
+  std::vector<sunrealtype> inputs;  // cppcheck-suppress unusedStructMember
 
   SetupOptions setup_opts;
 
-  virtual realtype *get_tmp_state_vector() = 0;
-  virtual realtype *get_tmp_sparse_jacobian_data() = 0;
+  virtual sunrealtype *get_tmp_state_vector() = 0;
+  virtual sunrealtype *get_tmp_sparse_jacobian_data() = 0;
 
 protected:
-  std::vector<realtype> tmp_state_vector;
-  std::vector<realtype> tmp_sparse_jacobian_data;
+  std::vector<sunrealtype> tmp_state_vector;
+  std::vector<sunrealtype> tmp_sparse_jacobian_data;
 };
 
 #endif // PYBAMM_IDAKLU_EXPRESSION_SET_HPP

src/pybammsolvers/idaklu_source/Expressions/Casadi/CasadiFunctions.cpp

@@ -61,8 +61,8 @@ const std::vector<expr_int>& CasadiFunction::get_col() {
   return m_cols;
 }
 
-void CasadiFunction::operator()(const std::vector<realtype*>& inputs,
-                                const std::vector<realtype*>& results)
+void CasadiFunction::operator()(const std::vector<sunrealtype*>& inputs,
+                                const std::vector<sunrealtype*>& results)
 {
   DEBUG("CasadiFunction operator() with inputs and results: " << m_func.name());
 

src/pybammsolvers/idaklu_source/Expressions/Casadi/CasadiFunctions.hpp

@@ -22,8 +22,8 @@ class CasadiFunction : public Expression
 
   // Method overrides
   void operator()() override;
-  void operator()(const std::vector<realtype*>& inputs,
-                  const std::vector<realtype*>& results) override;
+  void operator()(const std::vector<sunrealtype*>& inputs,
+                  const std::vector<sunrealtype*>& results) override;
   expr_int out_shape(int k) override;
   expr_int nnz() override;
   expr_int nnz_out() override;
@@ -144,10 +144,10 @@ class CasadiFunctions : public ExpressionSet<CasadiFunction>
   std::vector<CasadiFunction> dvar_dy_fcns_casadi;
   std::vector<CasadiFunction> dvar_dp_fcns_casadi;
 
-  realtype* get_tmp_state_vector() override {
+  sunrealtype* get_tmp_state_vector() override {
     return tmp_state_vector.data();
   }
-  realtype* get_tmp_sparse_jacobian_data() override {
+  sunrealtype* get_tmp_sparse_jacobian_data() override {
     return tmp_sparse_jacobian_data.data();
   }
 };

src/pybammsolvers/idaklu_source/Expressions/IREE/IREEFunction.hpp

@@ -22,8 +22,8 @@ class IREEFunction : public Expression
 
   // Method overrides
   void operator()() override;
-  void operator()(const std::vector<realtype*>& inputs,
-                  const std::vector<realtype*>& results) override;
+  void operator()(const std::vector<sunrealtype*>& inputs,
+                  const std::vector<sunrealtype*>& results) override;
   expr_int out_shape(int k) override;
   expr_int nnz() override;
   expr_int nnz_out() override;

src/pybammsolvers/idaklu_source/Expressions/IREE/IREEFunctions.cpp

@@ -177,7 +177,7 @@ void IREEFunction::evaluate(int n_outputs) {
   // Copy results to output array
   for(size_t k=0; k<n_outputs; k++) {
     for(size_t j=0; j<result[k].size(); j++) {
-      m_res[k][j] = static_cast<realtype>(result[k][j]);
+      m_res[k][j] = static_cast<sunrealtype>(result[k][j]);
     }
   }
 
@@ -213,8 +213,8 @@ const std::vector<expr_int>& IREEFunction::get_col() {
   return m_func.col;
 }
 
-void IREEFunction::operator()(const std::vector<realtype*>& inputs,
-                                const std::vector<realtype*>& results)
+void IREEFunction::operator()(const std::vector<sunrealtype*>& inputs,
+                                const std::vector<sunrealtype*>& results)
 {
   DEBUG("IreeFunction operator() with inputs and results");
   // Set-up input arguments, provide result vector, then execute function

src/pybammsolvers/idaklu_source/Expressions/IREE/IREEFunctions.hpp

@@ -130,10 +130,10 @@ class IREEFunctions : public ExpressionSet<IREEFunction>
   std::vector<IREEFunction> dvar_dy_fcns_iree;
   std::vector<IREEFunction> dvar_dp_fcns_iree;
 
-  realtype* get_tmp_state_vector() override {
+  sunrealtype* get_tmp_state_vector() override {
     return tmp_state_vector.data();
   }
-  realtype* get_tmp_sparse_jacobian_data() override {
+  sunrealtype* get_tmp_sparse_jacobian_data() override {
     return tmp_sparse_jacobian_data.data();
   }
 

src/pybammsolvers/idaklu_source/IDAKLUSolver.hpp

@@ -28,11 +28,11 @@ class IDAKLUSolver
    * @brief Abstract solver method that executes the solver
    */
   virtual SolutionData solve(
-    const std::vector<realtype> &t_eval,
-    const std::vector<realtype> &t_interp,
-    const realtype *y0,
-    const realtype *yp0,
-    const realtype *inputs,
+    const std::vector<sunrealtype> &t_eval,
+    const std::vector<sunrealtype> &t_interp,
+    const sunrealtype *y0,
+    const sunrealtype *yp0,
+    const sunrealtype *inputs,
     bool save_adaptive_steps,
     bool save_interp_steps
   ) = 0;

src/pybammsolvers/idaklu_source/IDAKLUSolverGroup.cpp

@@ -13,10 +13,10 @@ std::vector<Solution> IDAKLUSolverGroup::solve(
   // If t_interp is empty, save all adaptive steps
   bool save_adaptive_steps =  t_interp_np.size() == 0;
 
-  const realtype* t_eval_begin = t_eval_np.data();
-  const realtype* t_eval_end = t_eval_begin + t_eval_np.size();
-  const realtype* t_interp_begin = t_interp_np.data();
-  const realtype* t_interp_end = t_interp_begin + t_interp_np.size();
+  const sunrealtype* t_eval_begin = t_eval_np.data();
+  const sunrealtype* t_eval_end = t_eval_begin + t_eval_np.size();
+  const sunrealtype* t_interp_begin = t_interp_np.data();
+  const sunrealtype* t_interp_end = t_interp_begin + t_interp_np.size();
 
   // Process the time inputs
   // 1. Get the sorted and unique t_eval vector
@@ -97,9 +97,9 @@ std::vector<Solution> IDAKLUSolverGroup::solve(
   const std::size_t solves_per_thread = number_of_groups / m_solvers.size();
   const std::size_t remainder_solves = number_of_groups % m_solvers.size();
 
-  const realtype *y0 = y0_np.data();
-  const realtype *yp0 = yp0_np.data();
-  const realtype *inputs_data = inputs.data();
+  const sunrealtype *y0 = y0_np.data();
+  const sunrealtype *yp0 = yp0_np.data();
+  const sunrealtype *inputs_data = inputs.data();
 
   std::vector<SolutionData> results(number_of_groups);
 
@@ -111,9 +111,9 @@ std::vector<Solution> IDAKLUSolverGroup::solve(
     try {
       for (int j = 0; j < solves_per_thread; j++) {
         const std::size_t index = i * solves_per_thread + j;
-        const realtype *y = y0 + index * y0_np.shape(1);
-        const realtype *yp = yp0 + index * yp0_np.shape(1);
-        const realtype *input = inputs_data + index * inputs.shape(1);
+        const sunrealtype *y = y0 + index * y0_np.shape(1);
+        const sunrealtype *yp = yp0 + index * yp0_np.shape(1);
+        const sunrealtype *input = inputs_data + index * inputs.shape(1);
         results[index] = m_solvers[i]->solve(t_eval, t_interp, y, yp, input, save_adaptive_steps, save_interp_steps);
       }
     } catch (std::exception &e) {
@@ -132,9 +132,9 @@ std::vector<Solution> IDAKLUSolverGroup::solve(
 
   for (int i = 0; i < remainder_solves; i++) {
     const std::size_t index = number_of_groups - remainder_solves + i;
-    const realtype *y = y0 + index * y0_np.shape(1);
-    const realtype *yp = yp0 + index * yp0_np.shape(1);
-    const realtype *input = inputs_data + index * inputs.shape(1);
+    const sunrealtype *y = y0 + index * y0_np.shape(1);
+    const sunrealtype *yp = yp0 + index * yp0_np.shape(1);
+    const sunrealtype *input = inputs_data + index * inputs.shape(1);
     results[index] = m_solvers[i]->solve(t_eval, t_interp, y, yp, input, save_adaptive_steps, save_interp_steps);
   }