fixes issues from ticket:4395

  - added runtime option newtonMaxStepFactor, which adjustes the maximum
    step size of a newton step in kinsol.
  - used as: mxnewtstep = maxStepFactor * norm2(xScaling)
  - added a much better initial value for the maximum step size

Author: Willi Braun

SimulationRuntime/c/simulation/simulation_runtime.cpp

@@ -862,6 +862,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_NEWTON_MAX_STEP_FACTOR]) {
+    maxStepFactor = atof(omc_flagValue[FLAG_NEWTON_MAX_STEP_FACTOR]);
+    infoStreamPrint(LOG_STDOUT, 0, "Maximum step size factor for a Newton step 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);

SimulationRuntime/c/simulation/solver/kinsolSolver.c

@@ -102,6 +102,7 @@ typedef struct NLS_KINSOL_DATA
   /* ### tolerances ### */
   double fnormtol;        /* function tolerance */
   double scsteptol;       /* step tolerance */
+  double maxstepfactor;   /* maximum newton step factor mxnewtstep = maxstepfactor * norm2(xScaling) */
 
   /* ### work arrays ### */
   N_Vector initialGuess;
@@ -187,6 +188,8 @@ int nlsKinsolAllocate(int size, NONLINEAR_SYSTEM_DATA *nlsData, int linearSolver
   kinsolData->fnormtol  = sqrt(newtonFTol);     /* function tolerance */
   kinsolData->scsteptol = sqrt(newtonXTol);     /* step tolerance */
 
+  kinsolData->maxstepfactor = maxStepFactor;     /* step tolerance */
+
   kinsolData->initialGuess = N_VNew_Serial(size);
   kinsolData->xScale = N_VNew_Serial(size);
   kinsolData->fScale = N_VNew_Serial(size);
@@ -652,6 +655,15 @@ void nlsKinsolInfoPrint(const char *module, const char *function, char *msg, voi
   }
 }
 
+static
+void nlsKinsolSetMaxNewtonStep(NLS_KINSOL_DATA *kinsolData, double maxstepfactor)
+{
+  /* set maximum step size */
+  N_VConst(maxstepfactor, kinsolData->fTmp);
+  KINSetMaxNewtonStep(kinsolData->kinsolMemory, N_VWL2Norm(kinsolData->xScale, kinsolData->fTmp));
+
+}
+
 static
 void nlsKinsolResetInitial(DATA* data, NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM_DATA* nlsData, int mode)
 {
@@ -775,6 +787,7 @@ int nlsKinsolConfigPrint(NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM_DATA *nls
   infoStreamPrint(LOG_NLS_V, 0, "KINSOL max iterations %d", 20*kinsolData->size);
   infoStreamPrint(LOG_NLS_V, 0, "KINSOL strategy %d", kinsolData->kinsolStrategy);
   infoStreamPrint(LOG_NLS_V, 0, "KINSOL current retry %d", kinsolData->retries);
+  infoStreamPrint(LOG_NLS_V, 0, "KINSOL max step %g", (*(KINMem)kinsolData->kinsolMemory).kin_mxnstepin);
 
   messageClose(LOG_NLS_V);
 
@@ -805,13 +818,16 @@ int nlsKinsolErrorHandler(int errorCode, DATA *data, NONLINEAR_SYSTEM_DATA *nlsD
     break;
   /* just retry with new initial guess */
   case KIN_MXNEWT_5X_EXCEEDED:
-    warningStreamPrint(LOG_NLS, 0, "initial guess was too far away from the solution. Try again.\n");
+    warningStreamPrint(LOG_NLS, 0, "Newton step exceed the maximum step size several times. Try again after increasing maximum step size.\n");
+    kinsolData->maxstepfactor *= 1e5;
+    nlsKinsolSetMaxNewtonStep(kinsolData, kinsolData->maxstepfactor);
     return 1;
     break;
   /* just retry without line search */
   case KIN_LINESEARCH_NONCONV:
     warningStreamPrint(LOG_NLS, 0, "kinsols line search did not convergence. Try without.\n");
     kinsolData->kinsolStrategy = KIN_NONE;
+    kinsolData->retries--;
     return 1;
   /* maybe happened because of an out-dated factorization, so just retry  */
   case KIN_LSETUP_FAIL:
@@ -927,6 +943,9 @@ int nlsKinsolSolve(DATA *data, threadData_t *threadData, int sysNumber)
   /* set f scaling */
   nlsKinsolFScaling(data, kinsolData, nlsData, SCALING_JACOBIAN);
 
+  /* set maximum step size */
+  nlsKinsolSetMaxNewtonStep(kinsolData, kinsolData->maxstepfactor);
+
   /*  */
   do{
     /* dump configuration */

SimulationRuntime/c/simulation/solver/model_help.c

@@ -55,6 +55,7 @@ double linearSparseSolverMaxDensity = 0.2;
 int linearSparseSolverMinSize = 4001;
 double nonlinearSparseSolverMaxDensity = 0.2;
 int nonlinearSparseSolverMinSize = 10001;
+double maxStepFactor = 1e12;
 double newtonXTol = 1e-12;
 double newtonFTol = 1e-12;
 double steadyStateTol = 1e-3;

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 maxStepFactor;
 extern double steadyStateTol;
 extern const size_t SIZERINGBUFFER;
 extern int compiledInDAEMode;

SimulationRuntime/c/util/simulation_options.c

@@ -85,6 +85,7 @@ const char *FLAG_NAME[FLAG_MAX+1] = {
   /* FLAG_MAX_STEP_SIZE */         "maxStepSize",
   /* FLAG_MEASURETIMEPLOTFORMAT */ "measureTimePlotFormat",
   /* FLAG_NEWTON_FTOL */           "newtonFTol",
+  /* FLAG_NEWTON_MAX_STEP_FACTOR */"newtonMaxStepFactor",
   /* FLAG_NEWTON_XTOL */           "newtonXTol",
   /* FLAG_NEWTON_STRATEGY */       "newton",
   /* FLAG_NLS */                   "nls",
@@ -174,6 +175,7 @@ const char *FLAG_DESC[FLAG_MAX+1] = {
   /* FLAG_MAX_STEP_SIZE */         "value specifies maximum absolute step size for supported solver",
   /* FLAG_MEASURETIMEPLOTFORMAT */ "value specifies the output format of the measure time functionality",
   /* FLAG_NEWTON_FTOL */           "[double (default 1e-12)] tolerance respecting residuals for updating solution vector in Newton solver",
+  /* FLAG_NEWTON_MAX_STEP_FACTOR */"[double (default 1e12)] maximum newton step factor mxnewtstep = maxStepFactor * norm2(xScaling). Used currently only by kinsol.",
   /* FLAG_NEWTON_XTOL */           "[double (default 1e-12)] tolerance respecting newton correction (delta_x) for updating solution vector in Newton solver",
   /* FLAG_NEWTON_STRATEGY */       "value specifies the damping strategy for the newton solver",
   /* FLAG_NLS */                   "value specifies the nonlinear solver",
@@ -353,6 +355,9 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
   "  Tolerance respecting residuals for updating solution vector in Newton solver."
   "  Solution is accepted if the (scaled) 2-norm of the residuals is smaller than the tolerance newtonFTol and the (scaled) newton correction (delta_x) is smaller than the tolerance newtonXTol."
   "  The value is a Double with default value 1e-12.",
+  /* FLAG_NEWTON_MAX_STEP_FACTOR */
+  "  Maximum newton step factor mxnewtstep = maxStepFactor * norm2(xScaling). "
+  "  Used currently only by kinsol.",
   /* FLAG_NEWTON_XTOL */
   "  Tolerance respecting newton correction (delta_x) for updating solution vector in Newton solver."
   "  Solution is accepted if the (scaled) 2-norm of the residuals is smaller than the tolerance newtonFTol and the (scaled) newton correction (delta_x) is smaller than the tolerance newtonXTol."
@@ -496,6 +501,7 @@ const int FLAG_TYPE[FLAG_MAX] = {
   /* FLAG_MAX_STEP_SIZE */         FLAG_TYPE_OPTION,
   /* FLAG_MEASURETIMEPLOTFORMAT */ FLAG_TYPE_OPTION,
   /* FLAG_NEWTON_FTOL */           FLAG_TYPE_OPTION,
+  /* FLAG_NEWTON_MAX_STEP_FACTOR */FLAG_TYPE_OPTION,
   /* FLAG_NEWTON_XTOL */           FLAG_TYPE_OPTION,
   /* FLAG_NEWTON_STRATEGY */       FLAG_TYPE_OPTION,
   /* FLAG_NLS */                   FLAG_TYPE_OPTION,

SimulationRuntime/c/util/simulation_options.h

@@ -93,6 +93,7 @@ enum _FLAG
   FLAG_MAX_STEP_SIZE,
   FLAG_MEASURETIMEPLOTFORMAT,
   FLAG_NEWTON_FTOL,
+  FLAG_NEWTON_MAX_STEP_FACTOR,
   FLAG_NEWTON_XTOL,
   FLAG_NEWTON_STRATEGY,
   FLAG_NLS,