Update SUNDIALS to 7.8.0

src/OMSimulatorLib/AlgLoop.cpp

@@ -62,57 +62,49 @@ inline bool checkFlag(int flag, std::string functionName)
 }
 
 /**
- * @brief Error handler function given to KINSOL.
+ * @brief Error handler function of type SUNErrHandlerFn.
  *
- * @param errorCode   Error code from KINSOL
- * @param module      Name of the module reporting the error.
- * @param function    Name of the function in which the error occurred.
- * @param msg         Error Message.
- * @param user_data   Pointer to user data. Unused.
- */
-void oms::KinsolSolver::sundialsErrorHandlerFunction(int errorCode, const char *module,
-                                                     const char *function, char *msg,
-                                                     void *user_data)
-{
-  KINSOL_USER_DATA* kinsolUserData;
-  std::string systNum = "unknown";
-  std::string mod = module;
-  std::string func = function;
-
-  if (user_data != NULL)
-  {
-    kinsolUserData = (KINSOL_USER_DATA *)user_data;
-    systNum = std::to_string(kinsolUserData->algLoopNumber);
-  }
-
-  logError("SUNDIALS_ERROR: [system] " + systNum + " [module] " + mod + " | [function] " + func
-           + " | [error_code] " + std::to_string(errorCode) + "\n" + std::string(msg));
-}
-
-/**
- * @brief Info handler function given to KINSOL.
- *
- * Will only print information when debug loging is active.
+ * Pushed onto the error handler stack of the SUNDIALS context, so it is called
+ * for errors reported by KINSOL and by any other SUNDIALS module created with
+ * that context.
  *
- * @param module      Name of the module reporting the information.
- * @param function    Name of the function reporting the information.
- * @param msg         Message.
- * @param user_data   Pointer to user data. Unused.
+ * @param line            Line number in the SUNDIALS source file where the
+ *                        error occurred.
+ * @param func_name       Name of the SUNDIALS function in which the error
+ *                        occurred.
+ * @param file            Name of the SUNDIALS source file where the error
+ *                        occurred.
+ * @param msg             Error message. Can be NULL, in which case the message
+ *                        belonging to `err_code` is used.
+ * @param err_code        Error code. Either a SUNErrCode or a package level
+ *                        code, e.g. one of the KIN_* constants.
+ * @param err_user_data   Pointer to user data given to
+ *                        SUNContext_PushErrHandler(), a KINSOL_USER_DATA
+ *                        pointer. Can be NULL.
+ * @param sunctx          SUNDIALS context reporting the error. Unused.
  */
-void oms::KinsolSolver::sundialsInfoHandlerFunction(const char *module, const char *function,
-                                                    char *msg, void *user_data)
+void oms::KinsolSolver::sundialsErrorHandlerFunction(int line, const char *func_name,
+                                                     const char *file, const char *msg,
+                                                     SUNErrCode err_code,
+                                                     void *err_user_data, SUNContext sunctx)
 {
   KINSOL_USER_DATA* kinsolUserData;
   std::string systNum = "unknown";
-  std::string mod = module;
-  std::string func = function;
+  std::string file_location = std::string(file) + ":" + std::to_string(line);
+  std::string func = func_name;
 
-  if (user_data != NULL) {
-    kinsolUserData = (KINSOL_USER_DATA *)user_data;
+  if (err_user_data != NULL)
+  {
+    kinsolUserData = (KINSOL_USER_DATA *)err_user_data;
     systNum = std::to_string(kinsolUserData->algLoopNumber);
   }
 
-  logDebug("SUNDIALS_INFO: [system] " + systNum + " [module] " + mod + " | [function] " + func + "\n" + std::string(msg));
+  /* Package level codes (KIN_* and friends) are not SUNErrCodes, so SUNGetErrMsg()
+     only makes sense when SUNDIALS did not supply a message - same rule as
+     SUNDIALS' own default handler. */
+  logError("SUNDIALS_ERROR: [system] " + systNum + " [at] " + file_location + " | [function] " + func
+           + " | [error_code] " + std::to_string(err_code)
+           + "\n" + std::string(msg ? msg : SUNGetErrMsg(err_code)));
 }
 
 /**
@@ -285,6 +277,8 @@ oms::KinsolSolver::~KinsolSolver()
   SUNMatDestroy(this->J);
   N_VDestroy_Serial(this->y);
 
+  SUNContext_Free(&(this->sunctx));
+
   delete((KINSOL_USER_DATA*)(this->user_data));
 }
 
@@ -299,7 +293,6 @@ oms::KinsolSolver::~KinsolSolver()
 oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double relativeTolerance, const bool useDirectionalDerivative)
 {
   int flag;
-  int printLevel;
   /* Held by a unique_ptr until it is complete, so that the error returns below
    * do not leak the object and the SUNDIALS memory it allocated so far. */
   std::unique_ptr<KinsolSolver> kinsolSolver(new KinsolSolver());
@@ -308,16 +301,41 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
 
   kinsolSolver->size = size;
 
+  /* Create the SUNDIALS context every other SUNDIALS object is created with */
+  if (SUNContext_Create(SUN_COMM_NULL, &kinsolSolver->sunctx) != SUN_SUCCESS)
+  {
+    logError("SUNDIALS_ERROR: SUNContext_Create() failed");
+    return NULL;
+  }
+
+  /* Mute SUNDIALS' own output, use OMSimulator's logger */
+  {
+    SUNLogger logger = NULL;
+    if (SUNContext_GetLogger(kinsolSolver->sunctx, &logger) == SUN_SUCCESS && logger != NULL)
+    {
+      SUNLogger_SetErrorFilename(logger, "");
+      SUNLogger_SetWarningFilename(logger, "");
+      SUNLogger_SetInfoFilename(logger, "");
+      SUNLogger_SetDebugFilename(logger, "");
+    }
+  }
+
   /* Allocate memory */
-  kinsolSolver->initialGuess = N_VNew_Serial(kinsolSolver->size);
-  kinsolSolver->uScale = N_VNew_Serial(kinsolSolver->size);
-  kinsolSolver->fScale = N_VNew_Serial(kinsolSolver->size);
-  kinsolSolver->fTmp = N_VNew_Serial(kinsolSolver->size);
-  kinsolSolver->y = N_VNew_Serial(kinsolSolver->size);
+  kinsolSolver->initialGuess = N_VNew_Serial(kinsolSolver->size, kinsolSolver->sunctx);
+  kinsolSolver->uScale = N_VNew_Serial(kinsolSolver->size, kinsolSolver->sunctx);
+  kinsolSolver->fScale = N_VNew_Serial(kinsolSolver->size, kinsolSolver->sunctx);
+  kinsolSolver->fTmp = N_VNew_Serial(kinsolSolver->size, kinsolSolver->sunctx);
+  kinsolSolver->y = N_VNew_Serial(kinsolSolver->size, kinsolSolver->sunctx);
+  if (!kinsolSolver->initialGuess || !kinsolSolver->uScale || !kinsolSolver->fScale
+      || !kinsolSolver->fTmp || !kinsolSolver->y)
+  {
+    logError("SUNDIALS_ERROR: N_VNew_Serial() failed");
+    return NULL;
+  }
   kinsolSolver->kinsolMemory = NULL;
 
   /* Create KINSOL memory block */
-  kinsolSolver->kinsolMemory = KINCreate();
+  kinsolSolver->kinsolMemory = KINCreate(kinsolSolver->sunctx);
   if (kinsolSolver->kinsolMemory == NULL)
   {
     logError("SUNDIALS_ERROR: KINCreate() failed");
@@ -329,34 +347,29 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
   flag = KINSetUserData(kinsolSolver->kinsolMemory, kinsolSolver->user_data);
   if (!checkFlag(flag, "KINSetUserData")) return NULL;
 
-  /* Set error handler and print level */
-  if (logDebugEnabled())
-  {
-    logDebug("SUNDIALS KINSOL: Set print level to maximum.");
-    printLevel = 3;
-  }
-  else
-  {
-    printLevel = 0;
-  }
-  flag = KINSetPrintLevel(kinsolSolver->kinsolMemory, printLevel);
-  if (!checkFlag(flag, "KINSetPrintLevel")) return NULL;
-
-  flag = KINSetErrHandlerFn(kinsolSolver->kinsolMemory, sundialsErrorHandlerFunction, kinsolSolver->user_data);
-  if (!checkFlag(flag, "KINSetErrHandlerFn")) return NULL;
-
-  flag = KINSetInfoHandlerFn(kinsolSolver->kinsolMemory, sundialsInfoHandlerFunction, kinsolSolver->user_data);
-  if (!checkFlag(flag, "KINSetInfoHandlerFn")) return NULL;
+  /* Set error handler. KINSOL's progress output goes through the SUNLogger */
+  flag = SUNContext_PushErrHandler(kinsolSolver->sunctx, sundialsErrorHandlerFunction, kinsolSolver->user_data);
+  if (!checkFlag(flag, "SUNContext_PushErrHandler")) return NULL;
 
   /* Initialize KINSOL object */
   flag = KINInit(kinsolSolver->kinsolMemory, nlsKinsolResiduals, kinsolSolver->initialGuess);
   if (!checkFlag(flag, "KINInit")) return NULL;
 
   /* Create matrix object */
-  kinsolSolver->J = SUNDenseMatrix(kinsolSolver->size, kinsolSolver->size);
+  kinsolSolver->J = SUNDenseMatrix(kinsolSolver->size, kinsolSolver->size, kinsolSolver->sunctx);
+  if (kinsolSolver->J == NULL)
+  {
+    logError("SUNDIALS_ERROR: SUNDenseMatrix() failed");
+    return NULL;
+  }
 
   /* Create linear solver object */
-  kinsolSolver->linSol = SUNLinSol_Dense(kinsolSolver->y, kinsolSolver->J);
+  kinsolSolver->linSol = SUNLinSol_Dense(kinsolSolver->y, kinsolSolver->J, kinsolSolver->sunctx);
+  if (kinsolSolver->linSol == NULL)
+  {
+    logError("SUNDIALS_ERROR: SUNLinSol_Dense() failed");
+    return NULL;
+  }
 
   /* Set linear solver */
   flag = KINSetLinearSolver(kinsolSolver->kinsolMemory, kinsolSolver->linSol, kinsolSolver->J);

src/OMSimulatorLib/AlgLoop.h

@@ -42,9 +42,11 @@
 #include "OMSimulator/Types.h"
 #include "DirectedGraph.h"
 
+#include <sundials/sundials_context.h>  /* SUNContext */
+#include <sundials/sundials_logger.h>   /* SUNLogger */
 #include <kinsol/kinsol.h>
 #include <nvector/nvector_serial.h>
-#include <sunlinsol/sunlinsol_dense.h>       /* Default dense linear solver */
+#include <sunlinsol/sunlinsol_dense.h>  /* Default dense linear solver */
 
 namespace oms
 {
@@ -79,6 +81,7 @@ namespace oms
     N_Vector fTmp = nullptr;          /* Vector used for tmp computations */
 
     /* kinsol internal data */
+    SUNContext sunctx = nullptr;      /* SUNDIALS simulation context */
     void* kinsolMemory = nullptr;
     void* user_data = nullptr;
     int size = 0;
@@ -91,8 +94,8 @@ namespace oms
     /* member function */
     static int nlsKinsolJac(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2);
     static int nlsKinsolResiduals(N_Vector u, N_Vector fval, void *user_data);
-    static void sundialsErrorHandlerFunction(int error_code, const char *module, const char *function, char *msg, void *user_data);
-    static void sundialsInfoHandlerFunction(const char *module, const char *function, char *msg, void *user_data);
+    static void sundialsErrorHandlerFunction(int line, const char *func, const char *file, const char *msg,
+                                             SUNErrCode err_code, void *err_user_data, SUNContext sunctx);
   };
 
   class AlgLoop

src/OMSimulatorLib/CMakeLists.txt

@@ -98,8 +98,8 @@ target_include_directories(OMSimulatorLib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 target_link_libraries(OMSimulatorLib
                         PUBLIC
                           oms::public_includes
-                          oms::3rd::kinsol
-                          oms::3rd::cvode
+                          oms::3rd::sundials::kinsol
+                          oms::3rd::sundials::cvode
                           oms::3rd::fmi4c
                           oms::3rd::minizip
                           oms::3rd::zlib
@@ -129,8 +129,8 @@ target_include_directories(OMSimulatorLib_static PUBLIC ${CMAKE_CURRENT_SOURCE_D
 target_link_libraries(OMSimulatorLib_static
                         PUBLIC
                           oms::public_includes
-                          oms::3rd::kinsol
-                          oms::3rd::cvode
+                          oms::3rd::sundials::kinsol
+                          oms::3rd::sundials::cvode
                           oms::3rd::fmi4c
                           oms::3rd::minizip
                           oms::3rd::zlib

src/OMSimulatorLib/SystemSC.cpp

@@ -44,7 +44,6 @@
 #include "ssd/Tags.h"
 
 #include <algorithm>
-#include <cmath>
 #include <cstring>
 #include <sstream>
 #include <iostream>
@@ -63,7 +62,7 @@ namespace
   }
 }
 
-int oms::cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data)
+int oms::cvode_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data)
 {
   SystemSC* system = (SystemSC*)user_data;
   oms_status_enu_t status;
@@ -111,7 +110,7 @@ int oms::cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data)
   return 0;
 }
 
-int oms::cvode_rhs_algebraic(realtype t, N_Vector y, N_Vector ydot, void* user_data)
+int oms::cvode_rhs_algebraic(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data)
 {
   SystemSC* system = (SystemSC*)user_data;
 
@@ -123,7 +122,7 @@ int oms::cvode_rhs_algebraic(realtype t, N_Vector y, N_Vector ydot, void* user_d
   return 0;
 }
 
-int oms::cvode_roots(realtype t, N_Vector y, realtype *gout, void *user_data)
+int oms::cvode_roots(sunrealtype t, N_Vector y, sunrealtype *gout, void *user_data)
 {
   logDebug("cvode_roots at time " + std::to_string(t));
   SystemSC* system = (SystemSC*)user_data;
@@ -377,7 +376,23 @@ oms_status_enu_t oms::SystemSC::initialize()
     if (algebraic)
       n_states = 1;
 
-    solverData.cvode.y = N_VNew_Serial(static_cast<long>(n_states));
+    /* Create the SUNDIALS context every other SUNDIALS object is created with */
+    if (SUNContext_Create(SUN_COMM_NULL, &solverData.cvode.sunctx) != SUN_SUCCESS)
+      logError("SUNDIALS_ERROR: SUNContext_Create() failed");
+
+    /* Mute SUNDIALS' own output, use OMSimulator's logger */
+    {
+      SUNLogger logger = NULL;
+      if (SUNContext_GetLogger(solverData.cvode.sunctx, &logger) == SUN_SUCCESS && logger != NULL)
+      {
+        SUNLogger_SetErrorFilename(logger, "");
+        SUNLogger_SetWarningFilename(logger, "");
+        SUNLogger_SetInfoFilename(logger, "");
+        SUNLogger_SetDebugFilename(logger, "");
+      }
+    }
+
+    solverData.cvode.y = N_VNew_Serial(static_cast<long>(n_states), solverData.cvode.sunctx);
     if (!solverData.cvode.y) logError("SUNDIALS_ERROR: N_VNew_Serial() failed - returned NULL pointer");
 
     if (algebraic)
@@ -388,7 +403,7 @@ oms_status_enu_t oms::SystemSC::initialize()
           NV_Ith_S(solverData.cvode.y, k) = states[j][i];
     //N_VPrint_Serial(solverData.cvode.y);
 
-    solverData.cvode.abstol = N_VNew_Serial(static_cast<long>(n_states));
+    solverData.cvode.abstol = N_VNew_Serial(static_cast<long>(n_states), solverData.cvode.sunctx);
     if (!solverData.cvode.abstol) logError("SUNDIALS_ERROR: N_VNew_Serial() failed - returned NULL pointer");
 
     if (algebraic)
@@ -401,7 +416,7 @@ oms_status_enu_t oms::SystemSC::initialize()
 
     // Call CVodeCreate to create the solver memory and specify the
     // Backward Differentiation Formula and the use of a Newton iteration
-    solverData.cvode.mem = CVodeCreate(CV_BDF);
+    solverData.cvode.mem = CVodeCreate(CV_BDF, solverData.cvode.sunctx);
     if (!solverData.cvode.mem) logError("SUNDIALS_ERROR: CVodeCreate() failed - returned NULL pointer");
 
     int flag = CVodeSetUserData(solverData.cvode.mem, (void*)this);
@@ -422,13 +437,13 @@ oms_status_enu_t oms::SystemSC::initialize()
     if (flag < 0) logError("SUNDIALS_ERROR: CVodeSVtolerances() failed with flag = " + std::to_string(flag));
 
     // Call N_VNew_Serial and SUNDenseMatrix to generate dense vector abd natrix for lin. solver module
-    solverData.cvode.liny = N_VNew_Serial(n_states);
+    solverData.cvode.liny = N_VNew_Serial(n_states, solverData.cvode.sunctx);
     if (solverData.cvode.liny == NULL) logError("SUNDIALS_ERROR: N_VNew_Serial() failed");
-    solverData.cvode.J = SUNDenseMatrix(n_states, n_states);
-    if (solverData.cvode.J == NULL) logError("SUNDIALS_ERROR: N_VNew_Serial() failed");
+    solverData.cvode.J = SUNDenseMatrix(n_states, n_states, solverData.cvode.sunctx);
+    if (solverData.cvode.J == NULL) logError("SUNDIALS_ERROR: SUNDenseMatrix() failed");
 
     // Call SUNLinSol_Dense to creat linear solver object
-    solverData.cvode.linSol = SUNLinSol_Dense(solverData.cvode.liny, solverData.cvode.J);
+    solverData.cvode.linSol = SUNLinSol_Dense(solverData.cvode.liny, solverData.cvode.J, solverData.cvode.sunctx);
     if (solverData.cvode.linSol == NULL) logError("SUNDIALS_ERROR: SUNLinSol_Dense() failed");
 
     // Call CVodeSetLinearSolver to set the dense linear solver */
@@ -503,6 +518,7 @@ oms_status_enu_t oms::SystemSC::terminate()
     N_VDestroy_Serial(solverData.cvode.y);
     N_VDestroy_Serial(solverData.cvode.abstol);
     CVodeFree(&(solverData.cvode.mem));
+    SUNContext_Free(&(solverData.cvode.sunctx));
     solverData.cvode.mem = NULL;
   }
 
@@ -570,6 +586,7 @@ oms_status_enu_t oms::SystemSC::reset()
     N_VDestroy_Serial(solverData.cvode.y);
     N_VDestroy_Serial(solverData.cvode.abstol);
     CVodeFree(&(solverData.cvode.mem));
+    SUNContext_Free(&(solverData.cvode.sunctx));
     solverData.cvode.mem = nullptr;
   }
 
@@ -862,12 +879,14 @@ oms_status_enu_t oms::SystemSC::doStepCVODE()
 
   while (time < end_time)
   {
-    logDebug("CVode: " + std::to_string(time) + " -> " + std::to_string(end_time));
+    const fmi3Float64 tout = std::min(tnext, end_time);
+
+    logDebug("CVode: " + std::to_string(time) + " -> " + std::to_string(tout));
     for (size_t j=0, k=0; j < fmus.size(); ++j)
       for (size_t i=0; i < nStates[j]; ++i, ++k)
         NV_Ith_S(solverData.cvode.y, k) = states[j][i];
 
-    flag = CVode(solverData.cvode.mem, std::min(tnext, end_time), solverData.cvode.y, &time, CV_NORMAL);
+    flag = CVode(solverData.cvode.mem, tout, solverData.cvode.y, &time, CV_NORMAL);
 
     for (size_t i = 0, j=0; i < fmus.size(); ++i)
     {