New feature: detect steady state

- New runtime flag -steadyState aborts simulation as soon as
  steady state is detected.
- New runtime flag -steadyStateTol=<value> overrides default
  tolerance to detect steady state.

Author: lochel

SimulationRuntime/c/simulation/simulation_runtime.cpp

@@ -863,6 +863,11 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data, threadData_t *thr
     infoStreamPrint(LOG_STDOUT, 0, "Tolerance for accepting accuracy in Newton solver changed to %g", newtonFTol);
   }
 
+  if(omc_flag[FLAG_STEADY_STATE_TOL]) {
+    steadyStateTol = atof(omc_flagValue[FLAG_STEADY_STATE_TOL]);
+    infoStreamPrint(LOG_STDOUT, 0, "Tolerance for steady state detection changed to %g", steadyStateTol);
+  }
+
   rt_tick(SIM_TIMER_INIT_XML);
   read_input_xml(data->modelData, data->simulationInfo);
   rt_accumulate(SIM_TIMER_INIT_XML);

SimulationRuntime/c/simulation/solver/model_help.c

@@ -57,6 +57,7 @@ double nonlinearSparseSolverMaxDensity = 0.2;
 int nonlinearSparseSolverMinSize = 10001;
 double newtonXTol = 1e-12;
 double newtonFTol = 1e-12;
+double steadyStateTol = 1e-3;
 const size_t SIZERINGBUFFER = 3;
 int compiledInDAEMode = 0;
 int compiledWithSymSolver = 0;

SimulationRuntime/c/simulation/solver/model_help.h

@@ -81,6 +81,7 @@ extern double nonlinearSparseSolverMaxDensity;
 extern int nonlinearSparseSolverMinSize;
 extern double newtonXTol;
 extern double newtonFTol;
+extern double steadyStateTol;
 extern const size_t SIZERINGBUFFER;
 extern int compiledInDAEMode;
 extern int compiledWithSymSolver;

SimulationRuntime/c/simulation/solver/perform_simulation.c

@@ -361,6 +361,27 @@ int prefixedName_performSimulation(DATA* data, threadData_t *threadData, SOLVER_
     }
 #endif
 
+    /* check for steady state */
+    if (omc_flag[FLAG_STEADY_STATE])
+    {
+      if (0 < data->modelData->nStates)
+      {
+        int i;
+        double maxDer = 0.0;
+        double currDer;
+        for(i=data->modelData->nStates; i<2*data->modelData->nStates; ++i)
+        {
+          currDer = fabs(data->localData[0]->realVars[i] / data->modelData->realVarsData[i].attribute.nominal);
+          if(maxDer < currDer)
+            maxDer = currDer;
+        }
+        if(maxDer < steadyStateTol)
+          throwStreamPrint(threadData, "steady state reached at time = %g\n  * max(|d(x_i)/dt|/nominal(x_i)) = %g\n  * relative tolerance = %g", solverInfo->currentTime, maxDer, steadyStateTol);
+      }
+      else
+        throwStreamPrint(threadData, "No states in model. Flag -steadyState can only be used if states are present.");
+    }
+
     omc_alloc_interface.collect_a_little();
 
     /* try */
@@ -472,6 +493,9 @@ int prefixedName_performSimulation(DATA* data, threadData_t *threadData, SOLVER_
 
   fmtClose(&fmt);
 
+  if (omc_flag[FLAG_STEADY_STATE])
+    warningStreamPrint(LOG_STDOUT, 0, "Steady state has not been reached.\nThis may be due to too restrictive relative tolerance (%g) or short stopTime (%g).", steadyStateTol, simInfo->stopTime);
+
   TRACE_POP
   return retValue;
 }

SimulationRuntime/c/util/simulation_options.c

@@ -109,6 +109,8 @@ const char *FLAG_NAME[FLAG_MAX+1] = {
   /* FLAG_RT */                    "rt",
   /* FLAG_S */                     "s",
   /* FLAG_SOLVER_STEPS */          "steps",
+  /* FLAG_STEADY_STATE */          "steadyState",
+  /* FLAG_STEADY_STATE_TOL */      "steadyStateTol",
   /* FLAG_UP_HESSIAN */            "keepHessian",
   /* FLAG_W */                     "w",
 
@@ -194,6 +196,8 @@ const char *FLAG_DESC[FLAG_MAX+1] = {
   /* FLAG_RT */                    "value specifies the scaling factor for real-time synchronization (0 disables)",
   /* FLAG_S */                     "value specifies the integration method",
   /* FLAG_SOLVER_STEPS */          "dumps the number of integration steps into the result file",
+  /* FLAG_STEADY_STATE */          "aborts if steady state is reached",
+  /* FLAG_STEADY_STATE_TOL */      "[double (default 1e-3)] This relative tolerance is used to detect steady state.",
   /* FLAG_UP_HESSIAN */            "value specifies the number of steps, which keep hessian matrix constant",
   /* FLAG_W */                     "shows all warnings even if a related log-stream is inactive",
 
@@ -416,6 +420,10 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
   "  Value specifies the integration method.",
   /* FLAG_SOLVER_STEPS */
   "  dumps the number of integration steps into the result file",
+  /* FLAG_STEADY_STATE */
+  "  Aborts the simulation if steady state is reached.",
+  /* FLAG_STEADY_STATE_TOL */
+  "This relative tolerance is used to detect steady state: max(|d(x_i)/dt|/nominal(x_i)) < steadyStateTol",
   /* FLAG_UP_HESSIAN */
   "  Value specifies the number of steps, which keep Hessian matrix constant.",
   /* FLAG_W */
@@ -503,6 +511,8 @@ const int FLAG_TYPE[FLAG_MAX] = {
   /* FLAG_RT */                    FLAG_TYPE_OPTION,
   /* FLAG_S */                     FLAG_TYPE_OPTION,
   /* FLAG_SOLVER_STEPS */          FLAG_TYPE_FLAG,
+  /* FLAG_STEADY_STATE */          FLAG_TYPE_FLAG,
+  /* FLAG_STEADY_STATE_TOL */      FLAG_TYPE_OPTION,
   /* FLAG_UP_HESSIAN */            FLAG_TYPE_OPTION,
   /* FLAG_W */                     FLAG_TYPE_FLAG
 };

SimulationRuntime/c/util/simulation_options.h

@@ -117,6 +117,8 @@ enum _FLAG
   FLAG_RT,
   FLAG_S,
   FLAG_SOLVER_STEPS,
+  FLAG_STEADY_STATE,
+  FLAG_STEADY_STATE_TOL,
   FLAG_UP_HESSIAN,
   FLAG_W,