ORNL · pelesh · Dec 10, 2024 · Dec 10, 2024 · Dec 10, 2024 · Dec 10, 2024
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -124,7 +124,7 @@ if(GRIDKIT_ENABLE_IPOPT)
   include(FindIpopt)
 endif()
 if(GRIDKIT_ENABLE_SUNDIALS)
-  find_package(SUNDIALS 6.0.0 REQUIRED CONFIG 
+  find_package(SUNDIALS 7.0.0 REQUIRED CONFIG 
                PATHS ${SUNDIALS_DIR}
                      ${SUNDIALS_DIR}/lib/cmake/sundials)
   message(STATUS "SUNDIALS configuration found: ${SUNDIALS_CONFIG}")

diff --git a/README.md b/README.md
@@ -14,7 +14,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) >= 6.0.0
+	- [SUNDIALS](https://github.com/LLNL/sundials) >= 7.0.0
 	- [Suitesparse](https://github.com/DrTimothyAldenDavis/SuiteSparse) >= 5.x (optional)
 	- [Ipopt](https://github.com/coin-or/Ipopt) >= 3.x (optional)
 - [CMake](https://cmake.org/) >= 3.12

diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp
@@ -61,12 +61,8 @@
 #include <iostream>
 #include <iomanip>
 
-#include <idas/idas_direct.h>              /* access to IDADls interface           */
 #include <idas/idas.h>
-
-//Sundials Sparse KLU
-#include <sunmatrix/sunmatrix_sparse.h>
-#include <sunlinsol/sunlinsol_klu.h>   
+#include <idas/idas_ls.h>
 
 #include "ModelEvaluator.hpp"
 #include "Ida.hpp"
@@ -84,7 +80,7 @@ namespace Sundials
         int retval = 0;
 
         // Create the SUNDIALS context that all SUNDIALS objects require
-        retval = SUNContext_Create(NULL, &context_);
+        retval = SUNContext_Create(SUN_COMM_NULL, &context_);
         checkOutput(retval, "SUNContext");
         solver_ = IDACreate(context_);
         tag_ = NULL;
@@ -142,7 +138,7 @@ namespace Sundials
         checkAllocation((void*) yp0_, "N_VClone");
 
         // Dummy initial time; will be overridden.
-        const realtype t0 = RCONST(0.0);
+        const sunrealtype t0 = SUN_RCONST(0.0);
 
         // Allocate and initialize IDA workspace
         retval = IDAInit(solver_, this->Residual, t0, yy_, yp_);
@@ -153,8 +149,8 @@ namespace Sundials
         checkOutput(retval, "IDASetUserData");
 
         // Set tolerances
-        realtype rtol;
-        realtype atol;
+        sunrealtype rtol;
+        sunrealtype atol;
 
         model_->setTolerances(rtol, atol); ///< \todo Function name should be "getTolerances"!
         retval = IDASStolerances(solver_, rtol, atol);
@@ -340,7 +336,7 @@ namespace Sundials
         int retval = 0;
 
         // Set all quadratures to zero
-        N_VConst(RCONST(0.0), q_);
+        N_VConst(SUN_RCONST(0.0), q_);
 
         // Initialize quadratures
         retval = IDAQuadReInit(solver_, q_);
@@ -424,8 +420,8 @@ namespace Sundials
     int Ida<ScalarT, IdxT>::initializeBackwardSimulation(real_type tf)
     {
         int retval = 0;
-        realtype rtol;
-        realtype atol;
+        sunrealtype rtol;
+        sunrealtype atol;
 
         model_->initializeAdjoint();
 
@@ -558,7 +554,7 @@ namespace Sundials
     }
 
     template <class ScalarT, typename IdxT>
-    int Ida<ScalarT, IdxT>::Residual(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
+    int Ida<ScalarT, IdxT>::Residual(sunrealtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
     {
         ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
@@ -574,7 +570,7 @@ namespace Sundials
     }
 
     template <class ScalarT, typename IdxT>
-    int Ida<ScalarT, IdxT>::Jac(realtype t, realtype cj, N_Vector yy, N_Vector yp, N_Vector resvec, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+    int Ida<ScalarT, IdxT>::Jac(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector resvec, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
     {
 
         ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
@@ -604,7 +600,7 @@ namespace Sundials
         }
 
         sunindextype *colvals = SUNSparseMatrix_IndexValues(J);
-        realtype *data = SUNSparseMatrix_Data(J);
+        sunrealtype  *data = SUNSparseMatrix_Data(J);
         //Copy data from model jac to sundials
         for (unsigned int i = 0; i < c.size(); i++ )
         {
@@ -616,7 +612,7 @@ namespace Sundials
     }
 
     template <class ScalarT, typename IdxT>
-    int Ida<ScalarT, IdxT>::Integrand(realtype tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void *user_data)
+    int Ida<ScalarT, IdxT>::Integrand(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void *user_data)
     {
         ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
@@ -632,7 +628,7 @@ namespace Sundials
     }
 
     template <class ScalarT, typename IdxT>
-    int Ida<ScalarT, IdxT>::adjointResidual(realtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rrB, void *user_data)
+    int Ida<ScalarT, IdxT>::adjointResidual(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rrB, void *user_data)
     {
         ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
@@ -651,7 +647,7 @@ namespace Sundials
 
 
     template <class ScalarT, typename IdxT>
-    int Ida<ScalarT, IdxT>::adjointIntegrand(realtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rhsQB, void *user_data)
+    int Ida<ScalarT, IdxT>::adjointIntegrand(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rhsQB, void *user_data)
     {
         ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
@@ -705,10 +701,10 @@ namespace Sundials
 
 
     template <class ScalarT, typename IdxT>
-    void Ida<ScalarT, IdxT>::printOutput(realtype t)
+    void Ida<ScalarT, IdxT>::printOutput(sunrealtype t)
     {
-        realtype *yval  = N_VGetArrayPointer_Serial(yy_);
-        realtype *ypval = N_VGetArrayPointer_Serial(yp_);
+        sunrealtype *yval  = N_VGetArrayPointer_Serial(yy_);
+        sunrealtype *ypval = N_VGetArrayPointer_Serial(yp_);
 
         std::cout << std::setprecision(5) << std::setw(7) << t << " ";
         for (IdxT i = 0; i < model_->size(); ++i)
@@ -723,9 +719,9 @@ namespace Sundials
     }
 
     template <class ScalarT, typename IdxT>
-    void Ida<ScalarT, IdxT>::printSpecial(realtype t, N_Vector y)
+    void Ida<ScalarT, IdxT>::printSpecial(sunrealtype t, N_Vector y)
     {
-        realtype *yval = N_VGetArrayPointer_Serial(y);
+        sunrealtype *yval = N_VGetArrayPointer_Serial(y);
         IdxT N = static_cast<IdxT>(N_VGetLength_Serial(y));
         std::cout << "{";
         std::cout << std::setprecision(5) << std::setw(7) << t;
@@ -791,10 +787,10 @@ namespace Sundials
         }
     }
 
-    // Compiler will prevent building modules with data type incompatible with realtype
-    template class Ida<realtype, long int>;
-    template class Ida<realtype, int>;
-    template class Ida<realtype, size_t>;
+    // Compiler will prevent building modules with data type incompatible with sunrealtype
+    template class Ida<sunrealtype, long int>;
+    template class Ida<sunrealtype, int>;
+    template class Ida<sunrealtype, size_t>;
 
 } // namespace Sundials
 } // namespace AnalysisManager
diff --git a/Solver/Dynamic/Ida.hpp b/Solver/Dynamic/Ida.hpp
@@ -157,30 +157,30 @@ namespace AnalysisManager
                 return NV_DATA_S(qB_);
             }
 
-            void printOutput(realtype t);
-            void printSpecial(realtype t, N_Vector x);
+            void printOutput(sunrealtype t);
+            void printSpecial(sunrealtype t, N_Vector x);
             void printFinalStats();
 
         private:
-            static int Residual(realtype t,
+            static int Residual(sunrealtype t,
                                 N_Vector yy, N_Vector yp,
                                 N_Vector rr, void *user_data);
 
-            static int Jac(realtype t, realtype cj,
+            static int Jac(sunrealtype t, sunrealtype cj,
                N_Vector yy, N_Vector yp, N_Vector resvec,
                SUNMatrix J, void *user_data,
                N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
 
-            static int Integrand(realtype t,
+            static int Integrand(sunrealtype t,
                                  N_Vector yy,   N_Vector yp,
                                  N_Vector rhsQ, void *user_data);
 
-            static int adjointResidual(realtype t,
+            static int adjointResidual(sunrealtype t,
                                        N_Vector yy,  N_Vector yp,
                                        N_Vector yyB, N_Vector ypB,
                                        N_Vector rrB, void *user_data);
 
-            static int adjointIntegrand(realtype t,
+            static int adjointIntegrand(sunrealtype t,
                                         N_Vector yy,  N_Vector yp,
                                         N_Vector yyB, N_Vector ypB,
                                         N_Vector rhsQB, void *user_data);

diff --git a/Solver/Optimization/DynamicConstraint.cpp b/Solver/Optimization/DynamicConstraint.cpp
@@ -297,8 +297,8 @@ void DynamicConstraint<ScalarT, IdxT>::finalize_solution(SolverReturn status,
 
 
 
-template class DynamicConstraint<realtype, long int>;
-template class DynamicConstraint<realtype, size_t>;
+template class DynamicConstraint<sunrealtype, long int>;
+template class DynamicConstraint<sunrealtype, size_t>;
 
 } // namespace IpoptInterface
 } // namespace AnalysisManager
diff --git a/Solver/SteadyState/Kinsol.cpp b/Solver/SteadyState/Kinsol.cpp
@@ -91,7 +91,7 @@ namespace Sundials
         int retval = 0;
 
         // Create the SUNDIALS context that all SUNDIALS objects require
-        retval = SUNContext_Create(NULL, &context_);
+        retval = SUNContext_Create(SUN_COMM_NULL, &context_);
         checkOutput(retval, "SUNContext");
 
         solver_ = KINCreate(context_);
@@ -139,13 +139,9 @@ namespace Sundials
         retval = KINSetUserData(solver_, model_);
         checkOutput(retval, "KINSetUserData");
 
-        // Set output verbosity level
-        retval = KINSetPrintLevel(solver_, 0);
-        checkOutput(retval, "KINSetPrintLevel");
-
         // Set tolerances
-        realtype fnormtol;  ///< Residual tolerance
-        realtype scsteptol; ///< Scaled step tolerance
+        sunrealtype fnormtol;  ///< Residual tolerance
+        sunrealtype scsteptol; ///< Scaled step tolerance
 
         model_->setTolerances(fnormtol, scsteptol); ///< \todo Function name should be "getTolerances"!
         retval = KINSetFuncNormTol(solver_, fnormtol);
@@ -270,7 +266,7 @@ namespace Sundials
     template <class ScalarT, typename IdxT>
     void Kinsol<ScalarT, IdxT>::printOutput()
     {
-        realtype *yval  = N_VGetArrayPointer_Serial(yy_);
+        sunrealtype *yval  = N_VGetArrayPointer_Serial(yy_);
 
         std::cout << std::setprecision(5) << std::setw(7);
         for (IdxT i = 0; i < model_->size(); ++i)
@@ -281,9 +277,9 @@ namespace Sundials
     }
 
     template <class ScalarT, typename IdxT>
-    void Kinsol<ScalarT, IdxT>::printSpecial(realtype t, N_Vector y)
+    void Kinsol<ScalarT, IdxT>::printSpecial(sunrealtype t, N_Vector y)
     {
-        realtype *yval = N_VGetArrayPointer_Serial(y);
+        sunrealtype *yval = N_VGetArrayPointer_Serial(y);
         IdxT N = N_VGetLength_Serial(y);
         std::cout << "{";
         std::cout << std::setprecision(5) << std::setw(7) << t;
@@ -343,10 +339,10 @@ namespace Sundials
         }
     }
 
-    // Compiler will prevent building modules with data type incompatible with realtype
-    template class Kinsol<realtype, long int>;
-    template class Kinsol<realtype, int>;
-    template class Kinsol<realtype, size_t>;
+    // Compiler will prevent building modules with data type incompatible with sunrealtype
+    template class Kinsol<sunrealtype, long int>;
+    template class Kinsol<sunrealtype, int>;
+    template class Kinsol<sunrealtype, size_t>;
 
 } // namespace Sundials
 } // namespace AnalysisManager
diff --git a/Solver/SteadyState/Kinsol.hpp b/Solver/SteadyState/Kinsol.hpp
@@ -163,22 +163,22 @@ namespace AnalysisManager
             // }
 
             void printOutput();
-            void printSpecial(realtype t, N_Vector x);
+            void printSpecial(sunrealtype t, N_Vector x);
             void printFinalStats();
 
         private:
             static int Residual(N_Vector yy, N_Vector rr, void *user_data);
 
-            // static int Integrand(realtype t,
+            // static int Integrand(sunrealtype t,
             //                      N_Vector yy,   N_Vector yp,
             //                      N_Vector rhsQ, void *user_data);
 
-            // static int adjointResidual(realtype t,
+            // static int adjointResidual(sunrealtype t,
             //                            N_Vector yy,  N_Vector yp,
             //                            N_Vector yyB, N_Vector ypB,
             //                            N_Vector rrB, void *user_data);
 
-            // static int adjointIntegrand(realtype t,
+            // static int adjointIntegrand(sunrealtype t,
             //                             N_Vector yy,  N_Vector yp,
             //                             N_Vector yyB, N_Vector ypB,
             //                             N_Vector rhsQB, void *user_data);