Merge branch 'cvode_7.1.1' into SoybeanParameterization

yufenghe22 · yufenghe22 · commit c1a154ba7579 · 2024-12-17T15:51:16.000-06:00
diff --git a/bin/ePhotosynthesis.cpp b/bin/ePhotosynthesis.cpp
@@ -175,16 +175,7 @@ int main(int argc, const char* argv[]) {
           readFile(evn, inputs);
         if (fileProvided(atpcost, atpcost))
           readFile(atpcost, inputs);
-#ifdef SUNDIALS_CONTEXT_REQUIRED
-        SUNContext context;
-        if (SUNContext_Create(0, &context) < 0) {
-          std::cout << "SUNContext_Create failed" << std::endl;
-          throw std::runtime_error("SUNContext_Create");
-        }
-        Variables *theVars = new Variables(&context);
-#else // SUNDIALS_CONTEXT_REQUIRED
-	Variables *theVars = new Variables();
-#endif // SUNDIALS_CONTEXT_REQUIRED
+        Variables *theVars = new Variables();
         if (fileProvided(enzyme, enzymeFile)) {
             std::cerr << "ENZYME DATA PROVIDED" << std::endl;
             readFile(enzymeFile, theVars->EnzymeAct, true);
@@ -316,9 +307,6 @@ int main(int argc, const char* argv[]) {
             delete theVars;
         }
         delete maindriver;
-#ifdef SUNDIALS_CONTEXT_REQUIRED
-        SUNContext_Free(&context);
-#endif // SUNDIALS_CONTEXT_REQUIRED
         return (EXIT_SUCCESS);
     } catch (std::exception& e) {
         std::cout << "An error occurred: " << e.what() << std:: endl;
diff --git a/cmake/FindSUNDIALS.cmake b/cmake/FindSUNDIALS.cmake
@@ -98,8 +98,9 @@ find_path(SUNDIALS_INCLUDE_DIR sundials_types.h
         include
         include/sundials
 )
+cmake_path(GET SUNDIALS_INCLUDE_DIR PARENT_PATH SUNDIALS_INCLUDE_DIR_PARENT)
 set(SUNDIALS_INCLUDE_DIRS
-    "${SUNDIALS_INCLUDE_DIR}/.."
+    "${SUNDIALS_INCLUDE_DIR_PARENT}"
     "${SUNDIALS_INCLUDE_DIR}"
 )
 
diff --git a/include/Variables.hpp b/include/Variables.hpp
@@ -26,6 +26,7 @@
  *
  **********************************************************************************************************************************************/
 #include <map>
+#include <memory>
 #include "definitions.hpp"
 #include "vel/FIVel.hpp"
 #include "vel/BFVel.hpp"
@@ -56,17 +57,55 @@
 
 namespace ePhotosynthesis {
 
+#ifdef SUNDIALS_CONTEXT_REQUIRED
+extern std::shared_ptr<SUNContext> global_context; //!< Global context
+
+/**
+   Initialize the global SUNContext.
+*/
+void init_global_sundials_context();
+
+/**
+   Destroy the global SUNContext.
+ */
+void cleanup_global_sundials_context();
+  
+#endif // SUNDIALS_CONTEXT_REQUIRED
+    
 /**
   Structure to hold global variables
   */
 class Variables {
 public:
 #ifdef SUNDIALS_CONTEXT_REQUIRED
+
+private:
+    std::shared_ptr<SUNContext> _context; /**< Sundials context */
+    int _context_flags = 0; /**< Flags describing the context */
+    /** Flags describing context traits */
+    enum CONTEXT_FLAGS : int {
+        CONTEXT_FLAG_GLOBAL = 0x00000001, //< The context is global
+    };
+public:
     /**
        Construct a new variables instance for a given sundials context
        \param[in, out] ctx Sundials context.
     */
-    Variables(SUNContext* ctx) : context(ctx) {}
+    Variables(SUNContext* ctx = nullptr);
+    /**
+       Get the shared pointer to the Sundials context.
+       \return Sundials context shared pointer.
+     */
+    std::shared_ptr<SUNContext> context_ptr() {
+      return _context;
+    }
+    /**
+       Get the Sundials context.
+       \return Sundials context.
+     */
+    SUNContext& context() {
+      return *(_context.get());
+    }
 #else // SUNDIALS_CONTEXT_REQUIRED
     /**
        Construct a new variables instance for a given sundials context
@@ -86,9 +125,6 @@ class Variables {
     Variables& operator=(const Variables& other);
     friend std::ostream& operator<<(std::ostream &out, const Variables *in);
 
-#ifdef SUNDIALS_CONTEXT_REQUIRED
-    SUNContext* context = NULL; /**< Sundials context */
-#endif // SUNDIALS_CONTEXT_REQUIRED
     bool record = false;
     bool BF_FI_com = false;
     bool EPS_SUCS_com = false;
diff --git a/include/definitions.hpp b/include/definitions.hpp
@@ -47,7 +47,14 @@
 #ifndef sunrealtype
 #define sunrealtype double
 #endif
-#endif // 
+#endif //
+
+#if (SUNDIALS_VERSION_MAJOR >= 7)
+#include <sundials/sundials_errors.h>
+#else
+typedef int SUNErrCode;
+inline const char* SUNGetErrMsg(const SUNErrCode&) { return ""; }
+#endif
 
 #ifndef uint
 #define uint unsigned int
@@ -58,6 +65,9 @@
 #ifndef realtype
 #define realtype sunrealtype
 #endif
+#ifndef SUN_COMM_NULL
+#define SUN_COMM_NULL 0
+#endif
 
 namespace ePhotosynthesis {
 
diff --git a/include/drivers/driver.hpp b/include/drivers/driver.hpp
@@ -27,6 +27,7 @@
  **********************************************************************************************************************************************/
 
 #include <vector>
+#include <memory>
 
 #include "definitions.hpp"
 #include <nvector/nvector_serial.h>
@@ -146,7 +147,13 @@ class Driver {
     double smoothPenalty(const std::vector<double>& x, double threshold);
 
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-    SUNContext* context;
+private:
+    std::shared_ptr<SUNContext> _context;
+public:
+    /** Get the context */
+    SUNContext& context() {
+        return *(_context.get());
+    }
 #endif // SUNDIALS_CONTEXT_REQUIRED
     realtype abstol;            // absolute tolerance
     realtype reltol;            // relative tolerance
diff --git a/src/Variables.cpp b/src/Variables.cpp
@@ -1,12 +1,49 @@
+#include <cstdlib>
 #include "Variables.hpp"
 
+#ifdef SUNDIALS_CONTEXT_REQUIRED
+std::shared_ptr<SUNContext> ePhotosynthesis::global_context(NULL);
+
+void ePhotosynthesis::init_global_sundials_context() {
+    if (!global_context) {
+        SUNContext* context = new SUNContext();
+        SUNErrCode error_code = SUNContext_Create(SUN_COMM_NULL, context);
+        if (error_code) {
+            throw std::runtime_error("SUNContext_Create failed: " +
+                                     std::string(SUNGetErrMsg(error_code)));
+        }
+        global_context.reset(context);
+        std::atexit(cleanup_global_sundials_context);
+    }
+}
+
+void ePhotosynthesis::cleanup_global_sundials_context() {
+    if (global_context) {
+        global_context.reset();
+    }
+}
+
+#endif // SUNDIALS_CONTEXT_REQUIRED
+
 using namespace ePhotosynthesis;
 
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-Variables::Variables(const Variables& other) : Variables(other.context) {
+
+Variables::Variables(SUNContext* ctx) : _context(), _context_flags(0) {
+    if (_context == nullptr) {
+        init_global_sundials_context();
+        _context = global_context;
+        _context_flags |= CONTEXT_FLAG_GLOBAL;
+    } else {
+        _context.reset(ctx);
+    }
+}
+Variables::Variables(const Variables& other) :
+  _context(other._context), _context_flags(other._context_flags) {
     *this = other;
 }
-Variables::Variables(const Variables* other) : Variables(other->context) {
+Variables::Variables(const Variables* other) :
+  _context(other->_context), _context_flags(other->_context_flags) {
     *this = *other;
 }
 #else // SUNDIALS_CONTEXT_REQUIRED
@@ -58,7 +95,8 @@ Variables* Variables::deepcopy() const {
 
 Variables& Variables::operator=(const Variables &other) {
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-    context = other.context;
+    _context = other._context;
+    _context_flags = other._context_flags;
 #endif // SUNDIALS_CONTEXT_REQUIRED
     record = other.record;
     BF_FI_com = other.BF_FI_com;
diff --git a/src/drivers/CVodeMem.cpp b/src/drivers/CVodeMem.cpp
@@ -12,7 +12,7 @@ void CVodeMem::cvode_mem_init(Driver* driver, realtype t0, N_Vector y) {
     if (!initialized) {
         data = alloc_calc_data();
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-        _cvode_mem = CVodeCreate(CV_BDF, driver->context[0]);
+        _cvode_mem = CVodeCreate(CV_BDF, driver->context());
 #else // SUNDIALS_CONTEXT_REQUIRED
 	_cvode_mem = CVodeCreate(CV_BDF);
 #endif // SUNDIALS_CONTEXT_REQUIRED
diff --git a/src/drivers/driver.cpp b/src/drivers/driver.cpp
@@ -43,7 +43,7 @@ Driver::Driver(Variables *theVars, const double startTime, const double stepSize
 	       const int maxSubsteps, const double atol, const double rtol,
 	       const bool showWarn) {
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-  context = theVars->context;
+  _context = theVars->context_ptr();
 #endif // SUNDIALS_CONTEXT_REQUIRED
   inputVars = theVars;
   start = startTime;
@@ -72,7 +72,7 @@ arr Driver::run() {
         sunindextype N =  static_cast<long>(constraints.size());
         N_Vector y;
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-        y = N_VNew_Serial(N, *context);
+        y = N_VNew_Serial(N, context());
 #else // SUNDIALS_CONTEXT_REQUIRED
         y = N_VNew_Serial(N);
 #endif // SUNDIALS_CONTEXT_REQUIRED
@@ -95,9 +95,9 @@ arr Driver::run() {
         data->drv = this;
 
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-        SUNMatrix A = SUNDenseMatrix(N, N, *context);
-        SUNNonlinearSolver NLS = SUNNonlinSol_Newton(y, *context);
-        SUNLinearSolver LS = SUNLinSol_Dense(y, A, *context);
+        SUNMatrix A = SUNDenseMatrix(N, N, context());
+        SUNNonlinearSolver NLS = SUNNonlinSol_Newton(y, context());
+        SUNLinearSolver LS = SUNLinSol_Dense(y, A, context());
 #else // SUNDIALS_CONTEXT_REQUIRED
         SUNMatrix A = SUNDenseMatrix(N, N);
         SUNNonlinearSolver NLS = SUNNonlinSol_Newton(y);
diff --git a/src/rate/RedoxReg_Rate.cpp b/src/rate/RedoxReg_Rate.cpp
@@ -63,9 +63,9 @@ void RedoxReg::_Rate(const double t, const RedoxRegCondition* const RedoxReg_Con
 
         N_Vector y, constraints, scaling;
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-        y = N_VNew_Serial(1, theVars->context[0]);
-        constraints = N_VNew_Serial(1, theVars->context[0]);
-        scaling = N_VNew_Serial(1, theVars->context[0]);
+        y = N_VNew_Serial(1, theVars->context());
+        constraints = N_VNew_Serial(1, theVars->context());
+        scaling = N_VNew_Serial(1, theVars->context());
 #else // SUNDIALS_CONTEXT_REQUIRED
         y = N_VNew_Serial(1);
         constraints = N_VNew_Serial(1);
@@ -79,10 +79,10 @@ void RedoxReg::_Rate(const double t, const RedoxRegCondition* const RedoxReg_Con
         constraints_ptr[0] = ZERO;
         void *kmem = nullptr;
 #ifdef SUNDIALS_CONTEXT_REQUIRED
-        kmem = KINCreate(theVars->context[0]);
+        kmem = KINCreate(theVars->context());
 
-        SUNMatrix A = SUNDenseMatrix(1, 1, theVars->context[0]);
-        SUNLinearSolver LS = SUNLinSol_Dense(y, A, theVars->context[0]);
+        SUNMatrix A = SUNDenseMatrix(1, 1, theVars->context());
+        SUNLinearSolver LS = SUNLinSol_Dense(y, A, theVars->context());
 #else // SUNDIALS_CONTEXT_REQUIRED
         kmem = KINCreate();
 

Original file line number	Diff line number	Diff line change
`@@ -98,8 +98,9 @@ find_path(SUNDIALS_INCLUDE_DIR sundials_types.h`
`98`	`98`	`include`
`99`	`99`	`include/sundials`
`100`	`100`	`)`
	`101`	`+cmake_path(GET SUNDIALS_INCLUDE_DIR PARENT_PATH SUNDIALS_INCLUDE_DIR_PARENT)`
`101`	`102`	`set(SUNDIALS_INCLUDE_DIRS`
`102`		`- "${SUNDIALS_INCLUDE_DIR}/.."`
	`103`	`+ "${SUNDIALS_INCLUDE_DIR_PARENT}"`
`103`	`104`	`"${SUNDIALS_INCLUDE_DIR}"`
`104`	`105`	`)`
`105`	`106`