Improve event handling by using dense output instead of hermite (#10404)

OpenModelica · Mar 14, 2023 · 9d76c7b · 9d76c7b
1 parent 80e3f4b
commit 9d76c7b
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 7 deletions.
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_events.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_events.c
@@ -77,17 +77,17 @@ void bisection_gb(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo,
     c = 0.5 * (*a + *b);
 
     if (gbData->eventSearch == 0) {
-      /*calculates states at time c using hermite interpolation */
+      /*calculates states at time c using interpolation */
       if (isInnerIntegration) {
         gbfData = gbData->gbfData;
-        gb_interpolation(GB_INTERPOL_HERMITE,
+        gb_interpolation(gbfData->interpolation,
                     gbfData->timeLeft,  gbfData->yLeft,  gbfData->kLeft,
                     gbfData->timeRight, gbfData->yRight, gbfData->kRight,
                     c, gbfData->y1,
                     gbData->nStates, NULL,  gbData->nStates, gbfData->tableau, gbfData->x, gbfData->k);
         y = gbfData->y1;
       } else {
-        gb_interpolation(GB_INTERPOL_HERMITE,
+        gb_interpolation(gbData->interpolation,
                     gbData->timeLeft,  gbData->yLeft,  gbData->kLeft,
                     gbData->timeRight, gbData->yRight, gbData->kRight,
                     c, gbData->y1,
@@ -97,10 +97,16 @@ void bisection_gb(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo,
     } else {
       /*calculates states at time c using integration */
       if (isInnerIntegration) {
+        gbData->gbfData->time = gbData->gbfData->timeLeft;
+        memcpy(gbData->gbfData->yOld, gbData->gbfData->yLeft, gbData->nStates * sizeof(double));
+
         gbData->gbfData->stepSize = c - gbData->gbfData->time;
         gb_step_info = gbData->gbfData->step_fun(data, threadData, solverInfo);
         y = gbData->gbfData->y;
       } else {
+        gbData->time = gbData->timeLeft;
+        memcpy(gbData->yOld, gbData->yLeft, gbData->nStates * sizeof(double));
+
         gbData->stepSize = c - gbData->time;
         gb_step_info = gbData->step_fun(data, threadData, solverInfo);
         y = gbData->y;

diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.c
@@ -417,7 +417,7 @@ void freeRK_NLS_DATA(NONLINEAR_SYSTEM_DATA* nlsData) {
  * @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) {
+void set_kinsol_parameters(NLS_KINSOL_DATA* kin_mem, int numIter, int jacUpdate, int maxJacUpdate, double tolerance) {
     int flag;
 
     flag = KINSetNumMaxIters(kin_mem, numIter);
@@ -426,7 +426,7 @@ void set_kinsol_parameters(NLS_KINSOL_DATA* kin_mem, int numIter, int jacUpdate,
     checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNoInitSetup");
     flag = KINSetMaxSetupCalls(kin_mem, maxJacUpdate);
     checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetMaxSetupCalls");
-    flag = KINSetFuncNormTol(kin_mem, 100*DBL_EPSILON);
+    flag = KINSetFuncNormTol(kin_mem, tolerance);
     checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetFuncNormTol");
 }
 
@@ -486,12 +486,18 @@ NLS_SOLVER_STATUS solveNLS_gb(DATA *data, threadData_t *threadData, NONLINEAR_SY
     // Get kinsol data object
     NLS_KINSOL_DATA* kin_mem = ((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory;
 
-    set_kinsol_parameters(kin_mem, nlsData->size * 4, SUNTRUE, 10);
+    set_kinsol_parameters(kin_mem, nlsData->size * 4, SUNTRUE, 10, 100*DBL_EPSILON);
     solved = solveNLS(data, threadData, nlsData);
     /* Retry solution process with updated Jacobian */
     if (!solved) {
-      set_kinsol_parameters(kin_mem, nlsData->size * 4, SUNFALSE, 10);
+      infoStreamPrint(LOG_STDOUT, 0, "GBODE: Solution of NLS failed, Try with updated Jacobian.");
+      set_kinsol_parameters(kin_mem, nlsData->size * 4, SUNFALSE, 10, 100*DBL_EPSILON);
       solved = solveNLS(data, threadData, nlsData);
+      if (!solved) {
+        infoStreamPrint(LOG_STDOUT, 0, "GBODE: Solution of NLS failed, Try with less accuracy.");
+        set_kinsol_parameters(kin_mem, nlsData->size * 4, SUNFALSE, 10, 1000*DBL_EPSILON);
+        solved = solveNLS(data, threadData, nlsData);
+      }
     }
     if (ACTIVE_STREAM(LOG_GBODE_NLS)) get_kinsol_statistics(kin_mem);
   } else {