Fixed control flags in the test files of GBODE (#9215)

* Fixed control flags in the test files of GBODE

* Several fixes to the handling of the interpolation error
-   Added new hermite interpolation using just the left hand derivative
-   Error handling for too small step size
-   Use interpolation error for step size selection, if corresponding
    error control is activated

* Improved the activity diagram and step size control
- Added Informationen on rejection type
- Better initial step size of the inner integrator
- Improved handling of the interpolation error

* Several improvements:
-   Improved logging of fast states
-   Solved synchronization of equidistant time grid and event detection
-   Corrected statistics

Co-authored-by: AnHeuermann <AnHeuermann@users.noreply.github.com>

OMCompiler/SimulationRuntime/c/simulation/solver/epsilon.h

@@ -52,6 +52,7 @@ static const double DASSL_STEP_EPS = 1e-13;
  * defines the minimal step size
  */
 static const double MINIMAL_STEP_SIZE = 1e-12;
+static const double GB_MINIMAL_STEP_SIZE = 1e-20;
 
 /*
  * used in model_help.c for function setZCtol

OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.c

@@ -58,6 +58,7 @@ const char *GB_INTERPOL_METHOD_NAME[GB_INTERPOL_MAX] = {
   /* GB_INTERPOL_UNKNOWN */           "unknown",
   /* GB_INTERPOL_LIN */               "linear",
   /* GB_INTERPOL_HERMITE */           "hermite",
+  /* GB_INTERPOL_HERMITE_a */         "hermite_a",
   /* GB_INTERPOL_HERMITE_b */         "hermite_b",
   /* GB_INTERPOL_HERMITE_ERRCTRL */   "hermite_errctrl",
   /* GB_DENSE_OUTPUT */               "dense_output",
@@ -68,7 +69,8 @@ const char *GB_INTERPOL_METHOD_DESC[GB_INTERPOL_MAX] = {
   /* GB_INTERPOL_UNKNOWN */         "unknown",
   /* GB_INTERPOL_LIN */             "Linear interpolation (1st order)",
   /* GB_INTERPOL_HERMITE */         "Hermite interpolation (3rd order)",
-  /* GB_INTERPOL_HERMITE_b */       "Hermite interpolation (only for left hand side)",
+  /* GB_INTERPOL_HERMITE_a */       "Hermite interpolation (only for left hand side)",
+  /* GB_INTERPOL_HERMITE_b */       "Hermite interpolation (only for right hand side)",
   /* GB_INTERPOL_HERMITE_ERRCTRL */ "Hermite interpolation with error control",
   /* GB_DENSE_OUTPUT */             "use dense output formula for interpolation",
   /* GB_DENSE_OUTPUT_ERRCTRL */     "use dense output fomular with error control"
@@ -77,9 +79,9 @@ const char *GB_INTERPOL_METHOD_DESC[GB_INTERPOL_MAX] = {
 /**
  * @brief Get Runge-Kutta method from simulation flag FLAG_SR or FLAG_MR.
  *
- * Defaults to method RK_LOBA_IIIB_4 for single-rate.
+ * Defaults to method RK_ESDIRK4 for single-rate.
  *
- * Defaults to method RK_SDIRK2 for multi-rate method, if single-rate method is implicit.
+ * Defaults to method RK_ESDIRK4 for multi-rate method, if single-rate method is implicit.
  * Otherwise us same method as single-rate method.
  *
  * Returns GB_UNKNOWN if flag is not known.
@@ -131,15 +133,15 @@ enum GB_METHOD getGB_method(enum _FLAG flag) {
     case RK_LOBA_IIIC_4:
       // Default value for inner integration method
       // if the outer integration method is full implicit
-      return RK_ESDIRK3;
+      return RK_ESDIRK4;
     default:
       return singleRateMethod;
     }
   }
 
   // Default value for single-rate method
-  infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode method: esdirk3 [default]");
-  return RK_ESDIRK3;
+  infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode method: esdirk4 [default]");
+  return RK_ESDIRK4;
 }
 
 /**

OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.h

@@ -56,14 +56,14 @@ extern const char *GB_CTRL_METHOD_NAME[GB_CTRL_MAX];
 extern const char *GB_CTRL_METHOD_DESC[GB_CTRL_MAX];
 
 enum GB_INTERPOL_METHOD {
-  GB_INTERPOL_UNKNOWN = 0,          /* Unknown interpolation method */
-  GB_INTERPOL_LIN = 1,              /* Linear interpolation */
-  GB_INTERPOL_HERMITE = 2,          /* Hermite interpolation */
-  GB_INTERPOL_HERMITE_b = 3,        /* Hermite interpolation (only for left hand side)*/
-  GB_INTERPOL_HERMITE_ERRCTRL = 4,  /* Hermite interpolation with error control */
-  GB_DENSE_OUTPUT = 5,              /* Dense output, if available else hermite */
-  GB_DENSE_OUTPUT_ERRCTRL = 6,       /* Dense output, if available else hermite with error control */
-
+  GB_INTERPOL_UNKNOWN = 0,      /* Unknown interpolation method */
+  GB_INTERPOL_LIN,              /* Linear interpolation */
+  GB_INTERPOL_HERMITE,          /* Hermite interpolation */
+  GB_INTERPOL_HERMITE_a,        /* Hermite interpolation (only for left hand side)*/
+  GB_INTERPOL_HERMITE_b,        /* Hermite interpolation (only for right hand side)*/
+  GB_INTERPOL_HERMITE_ERRCTRL,  /* Hermite interpolation with error control */
+  GB_DENSE_OUTPUT,              /* Dense output, if available else hermite */
+  GB_DENSE_OUTPUT_ERRCTRL,      /* Dense output, if available else hermite with error control */
 
   GB_INTERPOL_MAX
 };

OMCompiler/SimulationRuntime/c/simulation/solver/gbode_ctrl.c

@@ -227,7 +227,7 @@ void getInitStepSize(DATA* data, threadData_t* threadData, DATA_GBODE* gbData)
     h1 = fmax(1e-6, h0*1e-3);
   }
 
-  gbData->stepSize = 0.5*fmin(100*h0,h1)*50;
+  gbData->stepSize = 0.5*fmin(100*h0,h1);
   gbData->lastStepSize = 0.0;
 
   sData->timeValue = gbData->time;

OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.c

@@ -46,7 +46,21 @@
 #include "nonlinearSystem.h"
 
 #include "util/jacobian_util.h"
+#include "util/rtclock.h"
 
+/**
+ * @brief Specific error handling of kinsol for gbode
+ *
+ * @param error_code  Reported error code
+ * @param module      Module of failure
+ * @param function    Nonlinear function
+ * @param msg         Message of failure
+ * @param data        Pointer to userData
+ */
+void GB_KINErrHandler(int error_code, const char *module, const char *function, char *msg, void *data) {
+// Preparation for specific error handling of the solution process of kinsol for gbode
+// This is needed to speed up simulation in case of failure
+}
 
 /**
  * @brief Initialize static data of non-linear system for DIRK.
@@ -264,8 +278,14 @@ NONLINEAR_SYSTEM_DATA* initRK_NLS_DATA(DATA* data, threadData_t* threadData, DAT
     solverData->initHomotopyData = NULL;
     nlsData->solverData = solverData;
 
-    int flag = KINSetNumMaxIters(((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory, nlsData->size * 10);
+    int flag;
+    NLS_KINSOL_DATA* kin_mem = ((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory;
+    flag = KINSetNumMaxIters(kin_mem, nlsData->size * 4);
     checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNumMaxIters");
+    flag = KINSetMaxSetupCalls(kin_mem, 10);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetMaxSetupCalls");
+    flag = KINSetErrHandlerFn(kin_mem, GB_KINErrHandler, NULL);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetErrHandlerFn");
     break;
   default:
     throwStreamPrint(NULL, "Memory allocation for NLS method %s not yet implemented.", GB_NLS_METHOD_NAME[gbData->nlsSolverMethod]);
@@ -389,6 +409,96 @@ void freeRK_NLS_DATA(NONLINEAR_SYSTEM_DATA* nlsData) {
   return;
 }
 
+/**
+ * @brief Set kinsol parameters
+ *
+ * @param kin_mem       Pointer to kinsol data object
+ * @param numIter       Number of nonlinear iterations
+ * @param jacUpdate     Update of jacobian necessary (SUNFALSE => yes)
+ * @param maxJacUpdate  Maximal number of constant jacobian
+ */
+void set_kinsol_parameters(NLS_KINSOL_DATA* kin_mem, int numIter, int jacUpdate, int maxJacUpdate) {
+    int flag;
+
+    flag = KINSetNumMaxIters(kin_mem, numIter);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNumMaxIters");
+    flag = KINSetNoInitSetup(kin_mem, jacUpdate);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNoInitSetup");
+    flag = KINSetMaxSetupCalls(kin_mem, maxJacUpdate);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetMaxSetupCalls");
+
+}
+
+/**
+ * @brief Get the kinsol statistics object
+ *
+ * @param kin_mem Pointer to kinsol data object
+ */
+void get_kinsol_statistics(NLS_KINSOL_DATA* kin_mem) {
+  int flag;
+  long int nIters, nFuncEvals, nJacEvals;
+  double fnorm;
+
+  // Get number of nonlinear iteration steps
+  flag = KINGetNumNonlinSolvIters(kin_mem, &nIters);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetNumNonlinSolvIters");
+
+  // Get the error of the residual function
+  flag = KINGetFuncNorm(kin_mem, &fnorm);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetFuncNorm");
+
+  // Get the number of jacobian evaluation
+  flag = KINGetNumJacEvals(kin_mem, &nJacEvals);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetNumJacEvals");
+
+  // Get the number of function evaluation
+  flag = KINGetNumFuncEvals(kin_mem, &nFuncEvals);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetNumFuncEvals");
+
+  // Report numbers
+  infoStreamPrint(LOG_GBODE_NLS, 0, "Kinsol statistics: nIters = %ld, nFuncEvals = %ld, nJacEvals = %ld,  fnorm:  %14.12g", nIters, nFuncEvals, nJacEvals, fnorm);
+}
+/**
+ * @brief Special treatment when solving non linear systems of equations
+ *
+ *        Will be described, when it is ready
+ *
+ * @param data                Pointer to runtime data struct.
+ * @param threadData          Thread data for error handling.
+ * @param nlsData             Non-linear solver data.
+ * @param gbData              Runge-Kutta method.
+ * @return NLS_SOLVER_STATUS  Return NLS_SOLVED on success and NLS_FAILED otherwise.
+ */
+NLS_SOLVER_STATUS solveNLS_gb(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nlsData, DATA_GBODE* gbData) {
+  struct dataSolver * solverData = (struct dataSolver *)nlsData->solverData;
+  NLS_SOLVER_STATUS solved;
+
+  // Debug nonlinear solution process
+  rtclock_t clock;
+  double cpu_time_used;
+
+  if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
+    rt_ext_tp_tick(&clock);
+  }
+
+  if (gbData->nlsSolverMethod == GB_NLS_KINSOL) {
+    // Get kinsol data object
+    NLS_KINSOL_DATA* kin_mem = ((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory;
+
+    set_kinsol_parameters(kin_mem, nlsData->size * 4, SUNFALSE, 10);
+    solved = solveNLS(data, threadData, nlsData);
+    if (ACTIVE_STREAM(LOG_GBODE_NLS)) get_kinsol_statistics(kin_mem);
+  } else {
+    solved = solveNLS(data, threadData, nlsData);
+  }
+
+  if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
+      cpu_time_used = rt_ext_tp_tock(&clock);
+      infoStreamPrint(LOG_GBODE_NLS, 0, "Time needed for solving the NLS:  %20.16g", cpu_time_used);
+  }
+
+  return solved;
+}
 
 /**
  * @brief Residual function for non-linear system of generic multi-step methods.

OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.h

@@ -49,6 +49,9 @@ NONLINEAR_SYSTEM_DATA* initRK_NLS_DATA(DATA* data, threadData_t* threadData, DAT
 NONLINEAR_SYSTEM_DATA* initRK_NLS_DATA_MR(DATA* data, threadData_t* threadData, DATA_GBODEF* gbfData);
 void freeRK_NLS_DATA( NONLINEAR_SYSTEM_DATA* nlsData);
 
+//Specific treatment of NLS within gbode
+NLS_SOLVER_STATUS solveNLS_gb(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nonlinsys, DATA_GBODE* gbData);
+
 // Residuum and Jacobian functions for diagonal implicit (DIRK) and implicit (IRK) Runge-Kutta methods.
 void residual_MS(RESIDUAL_USERDATA* userData, const double *xloc, double *res, const int *iflag);
 void residual_DIRK(RESIDUAL_USERDATA* userData, const double *xloc, double *res, const int *iflag);