further fixes to ticket:4395

 - adjust default linear solver to klu
 - fix tolerances
 - fix scaling factors

Author: Willi Braun

SimulationRuntime/c/simulation/simulation_runtime.cpp

@@ -306,14 +306,23 @@ static int getNewtonStrategy()
   return NEWTON_NONE;
 }
 
-static int getNlsLSSolver()
+static int getNlsLSSolver(int nlsSolver)
 {
   int i;
   const char *cflags = omc_flagValue[FLAG_NLS_LS];
   const string *method = cflags ? new string(cflags) : NULL;
 
   if(!method)
-    return NLS_LS_LAPACK; /* default method */
+  {
+    if (nlsSolver == NLS_KINSOL)
+    {
+      return NLS_LS_KLU; /* default kinsol linear solver method */
+    }
+    else
+    {
+      return NLS_LS_LAPACK; /* default method */
+    }
+  }
 
   for(i=1; i<NLS_LS_MAX; ++i)
     if(*method == NLS_LS_METHOD[i])
@@ -834,7 +843,7 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data, threadData_t *thr
   data->simulationInfo->lssMethod = getlinearSparseSolverMethod();
   data->simulationInfo->newtonStrategy = getNewtonStrategy();
   data->simulationInfo->nlsCsvInfomation = omc_flag[FLAG_NLS_INFO];
-  data->simulationInfo->nlsLinearSolver = getNlsLSSolver();
+  data->simulationInfo->nlsLinearSolver = getNlsLSSolver(data->simulationInfo->nlsMethod);
 
   if(omc_flag[FLAG_LSS_MAX_DENSITY]) {
     linearSparseSolverMaxDensity = atof(omc_flagValue[FLAG_LSS_MAX_DENSITY]);

SimulationRuntime/c/simulation/solver/kinsolSolver.c

@@ -185,8 +185,8 @@ int nlsKinsolAllocate(int size, NONLINEAR_SYSTEM_DATA *nlsData, int linearSolver
   kinsolData->linearSolverMethod = linearSolverMethod;
   kinsolData->solved = 0;
 
-  kinsolData->fnormtol  = sqrt(newtonFTol);     /* function tolerance */
-  kinsolData->scsteptol = sqrt(newtonXTol);     /* step tolerance */
+  kinsolData->fnormtol  = newtonFTol;     /* function tolerance */
+  kinsolData->scsteptol = newtonXTol;     /* step tolerance */
 
   kinsolData->maxstepfactor = maxStepFactor;     /* step tolerance */
 
@@ -317,9 +317,6 @@ static int nlsKinsolResiduals(N_Vector x, N_Vector f, void *userData)
   data->simulationInfo->nonlinearSystemData[sysNumber].residualFunc(dataAndThreadData, xdata,  fdata, (const int*) &iflag);
   iflag = 0;
 
-  _omc_printVectorWithEquationInfo(_omc_createVector(kinsolData->size, fdata),
-        "fResiduals", LOG_NLS_V, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber));
-
 #ifndef OMC_EMCC
     MMC_CATCH_INTERNAL(simulationJumpBuffer)
 #endif
@@ -695,7 +692,7 @@ void nlsKinsolXScaling(DATA* data, NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM
   {
     case SCALING_NOMINALSTART:
       for (i=0;i<nlsData->size;i++){
-        xScaling[i] = fmax(nlsData->nominal[i],fabs(xStart[i]));
+        xScaling[i] = 1/fmax(nlsData->nominal[i],fabs(xStart[i]));
       }
       break;
     case SCALING_ONES:
@@ -720,7 +717,6 @@ void nlsKinsolFScaling(DATA* data, NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM
       N_Vector tmp2 = N_VNew_Serial(kinsolData->size);
 
       /* calculate the right jacobian */
-      nlsKinsolResiduals(x, kinsolData->fTmp, &kinsolData->userData);
       if(nlsData->isPatternAvailable && kinsolData->linearSolverMethod == 3)
       {
         spJac = NewSparseMat(kinsolData->size,kinsolData->size,kinsolData->nnz);
@@ -750,6 +746,7 @@ void nlsKinsolFScaling(DATA* data, NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM
           fScaling[spJac->rowvals[i]] = fabs(spJac->data[i]);
         }
       }
+      N_VInv(kinsolData->fScale, kinsolData->fScale);
 
       N_VDestroy_Serial(tmp1);
       N_VDestroy_Serial(tmp2);
@@ -788,6 +785,7 @@ int nlsKinsolConfigPrint(NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM_DATA *nls
   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);
+  infoStreamPrint(LOG_NLS_V, 0, "KINSOL linear solver %d", kinsolData->linearSolverMethod);
 
   messageClose(LOG_NLS_V);
 
@@ -797,7 +795,7 @@ int nlsKinsolConfigPrint(NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM_DATA *nls
 static
 int nlsKinsolErrorHandler(int errorCode, DATA *data, NONLINEAR_SYSTEM_DATA *nlsData, NLS_KINSOL_DATA *kinsolData)
 {
-  int retValue, i, retValue2=0;
+  int retValue, i, retValue2=0, flag;
   double fNorm;
   double *xStart = NV_DATA_S(kinsolData->initialGuess);
   double *xScaling = NV_DATA_S(kinsolData->xScale);
@@ -830,7 +828,6 @@ int nlsKinsolErrorHandler(int errorCode, DATA *data, NONLINEAR_SYSTEM_DATA *nlsD
     kinsolData->retries--;
     return 1;
   /* maybe happened because of an out-dated factorization, so just retry  */
-  case KIN_LSETUP_FAIL:
   case KIN_LSOLVE_FAIL:
     warningStreamPrint(LOG_NLS, 0, "kinsols matrix need new factorization. Try again.\n");
     if (nlsData->isPatternAvailable){
@@ -842,6 +839,20 @@ int nlsKinsolErrorHandler(int errorCode, DATA *data, NONLINEAR_SYSTEM_DATA *nlsD
     warningStreamPrint(LOG_NLS, 0, "kinsols runs into issues retry with different configuration.\n");
     retValue = 1;
     break;
+  case KIN_LSETUP_FAIL:
+    /* in case of something goes wrong with the symbolic jacobian try the numerical */
+    if ( kinsolData->linearSolverMethod == 3 && nlsData->isPatternAvailable && nlsData->analyticalJacobianColumn != NULL){
+      flag = KINSlsSetSparseJacFn(kinsolData->kinsolMemory, nlsSparseJac);
+    }
+    if (checkReturnFlag(flag)){
+      errorStreamPrint(LOG_STDOUT, 0, "##KINSOL## Something goes wrong while initialize KINSOL Sparse Solver!");
+      return flag;
+    }
+    else
+    {
+      retValue = 1;
+    }
+    break;
   default:
     errorStreamPrint(LOG_STDOUT, 0, "kinsol has a serious solving issue ERROR %d\n", errorCode);
     return errorCode;