diff --git a/3rdParty b/3rdParty index 34b05ece4..714704485 160000 --- a/3rdParty +++ b/3rdParty @@ -1 +1 @@ -Subproject commit 34b05ece412f8dbf0376f7969e5404ae0242cf75 +Subproject commit 714704485cca375a76d3e5d27909d9efe0e1cb63 diff --git a/src/OMSimulatorLib/AlgLoop.cpp b/src/OMSimulatorLib/AlgLoop.cpp index 4587ebe88..10a44f85c 100644 --- a/src/OMSimulatorLib/AlgLoop.cpp +++ b/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(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); diff --git a/src/OMSimulatorLib/AlgLoop.h b/src/OMSimulatorLib/AlgLoop.h index 5c17241e8..78b50369f 100644 --- a/src/OMSimulatorLib/AlgLoop.h +++ b/src/OMSimulatorLib/AlgLoop.h @@ -42,9 +42,11 @@ #include "OMSimulator/Types.h" #include "DirectedGraph.h" +#include /* SUNContext */ +#include /* SUNLogger */ #include #include -#include /* Default dense linear solver */ +#include /* 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 diff --git a/src/OMSimulatorLib/CMakeLists.txt b/src/OMSimulatorLib/CMakeLists.txt index c7fae62f9..88af05d04 100644 --- a/src/OMSimulatorLib/CMakeLists.txt +++ b/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 diff --git a/src/OMSimulatorLib/SystemSC.cpp b/src/OMSimulatorLib/SystemSC.cpp index 306f4bac5..0c72d8f86 100644 --- a/src/OMSimulatorLib/SystemSC.cpp +++ b/src/OMSimulatorLib/SystemSC.cpp @@ -44,7 +44,6 @@ #include "ssd/Tags.h" #include -#include #include #include #include @@ -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(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(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(n_states)); + solverData.cvode.abstol = N_VNew_Serial(static_cast(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) { diff --git a/src/OMSimulatorLib/SystemSC.h b/src/OMSimulatorLib/SystemSC.h index 7452005ed..b1246a299 100644 --- a/src/OMSimulatorLib/SystemSC.h +++ b/src/OMSimulatorLib/SystemSC.h @@ -42,18 +42,20 @@ #include -#include /* prototypes for CVODE fcts., consts. */ -#include /* serial N_Vector types, fcts., macros */ -#include /* Default dense linear solver */ +#include /* SUNContext */ +#include /* SUNLogger */ +#include /* prototypes for CVODE fcts., consts. */ +#include /* serial N_Vector types, fcts., macros */ +#include /* Default dense linear solver */ namespace oms { class Model; class ComponentFMUME; class Component; - int cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data); - int cvode_rhs_algebraic(realtype t, N_Vector y, N_Vector ydot, void* user_data); - int cvode_roots(realtype t, N_Vector y, realtype *gout, void* user_data); + int cvode_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + int cvode_rhs_algebraic(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + int cvode_roots(sunrealtype t, N_Vector y, sunrealtype *gout, void* user_data); class SystemSC : public System { @@ -113,6 +115,7 @@ namespace oms struct SolverDataCVODE_t { + SUNContext sunctx; /* SUNDIALS simulation context */ void *mem; N_Vector y; SUNLinearSolver linSol; /* linear solver object */ @@ -127,9 +130,9 @@ namespace oms SolverDataCVODE_t cvode; } solverData; - friend int oms::cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data); - friend int oms::cvode_rhs_algebraic(realtype t, N_Vector y, N_Vector ydot, void* user_data); - friend int oms::cvode_roots(realtype t, N_Vector y, realtype *gout, void* user_data); + friend int oms::cvode_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + friend int oms::cvode_rhs_algebraic(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + friend int oms::cvode_roots(sunrealtype t, N_Vector y, sunrealtype *gout, void* user_data); }; }