Make KINSOL compatible with adaptive homotopy

This workaround sets -homotopyOnFirstTry, so trying without homotopy
first is not possible (for kinsol with adaptive homotopy).

Belonging to [master]:
  - #2185
  - OpenModelica/OpenModelica-testsuite#851

Author: ptaeuber

SimulationRuntime/c/simulation/solver/initialization/initialization.c

@@ -189,6 +189,11 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
   int homotopySupport = 0;
   int solveWithGlobalHomotopy;
   int adaptiveGlobal;
+  int kinsol = 0;
+
+#if !defined(OMC_MINIMAL_RUNTIME)
+  kinsol = (data->simulationInfo->nlsMethod == NLS_KINSOL);
+#endif
 
 #if !defined(OMC_NUM_NONLINEAR_SYSTEMS) || OMC_NUM_NONLINEAR_SYSTEMS>0
   for(i=0; i<mData->nNonLinearSystems; i++) {
@@ -227,12 +232,16 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
 #ifndef OMC_EMCC
   MMC_TRY_INTERNAL(simulationJumpBuffer)
 #endif
-    if (adaptiveGlobal)
-      data->callback->useHomotopy = 1;
-    data->simulationInfo->lambda = 1.0;
-    infoStreamPrint(LOG_INIT, 0, "Try to solve the initialization problem without homotopy first.");
-    data->callback->functionInitialEquations(data, threadData);
-    solveWithGlobalHomotopy = 0;
+    if (adaptiveGlobal && kinsol) {
+      infoStreamPrint(LOG_INIT, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the adaptive global approach in combination with KINSOL.");
+    } else {
+      if (adaptiveGlobal)
+        data->callback->useHomotopy = 1;
+      data->simulationInfo->lambda = 1.0;
+      infoStreamPrint(LOG_INIT, 0, "Try to solve the initialization problem without homotopy first.");
+      data->callback->functionInitialEquations(data, threadData);
+      solveWithGlobalHomotopy = 0;
+    }
 
     /* catch */
 #ifndef OMC_EMCC
@@ -241,7 +250,8 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
     if (adaptiveGlobal)
       data->callback->useHomotopy = 2;
     if(solveWithGlobalHomotopy) {
-      warningStreamPrint(LOG_ASSERT, 0, "Failed to solve the initialization problem without homotopy method. If homotopy is available the homotopy method is used now.");
+      if (!kinsol)
+        warningStreamPrint(LOG_ASSERT, 0, "Failed to solve the initialization problem without homotopy method. If homotopy is available the homotopy method is used now.");
       omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY] = 1;
       setAllParamsToStart(data);
       setAllVarsToStart(data);

SimulationRuntime/c/simulation/solver/nonlinearSystem.c

@@ -448,10 +448,10 @@ int initializeNonlinearSystems(DATA *data, threadData_t *threadData)
     case NLS_KINSOL:
       solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
       if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
+        // Try without homotopy not supported for kinsol
         // nlsKinsolAllocate(size-1, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
         // solverData->ordinaryData = nonlinsys[i].solverData;
-        // allocateHomotopyData(size-1, &(solverData->initHomotopyData));
-        throwStreamPrint(threadData, "kinsol solver not compatible with adaptive homotopy approaches");
+        allocateHomotopyData(size-1, &(solverData->initHomotopyData));
       } else {
         nlsKinsolAllocate(size, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
         solverData->ordinaryData = nonlinsys[i].solverData;
@@ -571,9 +571,10 @@ int freeNonlinearSystems(DATA *data, threadData_t *threadData)
       free(nonlinsys[i].solverData);
       break;
     case NLS_KINSOL:
-      nlsKinsolFree(&((struct dataSolver*) nonlinsys[i].solverData)->ordinaryData);
       if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
         freeHomotopyData(&((struct dataSolver*) nonlinsys[i].solverData)->initHomotopyData);
+      } else {
+        nlsKinsolFree(&((struct dataSolver*) nonlinsys[i].solverData)->ordinaryData);
       }
       free(nonlinsys[i].solverData);
       break;
@@ -890,11 +891,17 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
   int solveWithHomotopySolver = 0;
   int homotopyDeactivated = 0;
   int j;
+  int nlsLs;
+  int kinsol = 0;
   struct dataSolver *solverData;
   struct dataMixedSolver *mixedSolverData;
   char buffer[4096];
   FILE *pFile = NULL;
 
+#if !defined(OMC_MINIMAL_RUNTIME)
+  kinsol = (data->simulationInfo->nlsMethod == NLS_KINSOL);
+#endif
+
   data->simulationInfo->currentNonlinearSystemIndex = sysNumber;
 
   /* enable to avoid division by zero */
@@ -951,12 +958,16 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
   /* If homotopy is deactivated in this place or flag homotopyOnFirstTry is not set,
      solve the system with the selected solver */
   if (homotopyDeactivated || !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY]) {
-    if (!homotopyDeactivated && !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
-      infoStreamPrint(LOG_INIT, 0, "Try to solve nonlinear initial system %d without homotopy first.", sysNumber);
+    if (solveWithHomotopySolver && kinsol) {
+      infoStreamPrint(LOG_INIT, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the local global approach in combination with KINSOL.");
+    } else {
+      if (!homotopyDeactivated && !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
+        infoStreamPrint(LOG_INIT, 0, "Try to solve nonlinear initial system %d without homotopy first.", sysNumber);
 
-    data->simulationInfo->lambda = 1.0;
-    /* SOLVE! */
-    nonlinsys->solved = solveNLS(data, threadData, sysNumber);
+      data->simulationInfo->lambda = 1.0;
+      /* SOLVE! */
+      nonlinsys->solved = solveNLS(data, threadData, sysNumber);
+    }
   }
 
   /* The following cases are only valid for initial systems with homotopy */
@@ -965,7 +976,7 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
   /* If the adaptive local/global homotopy approach is activated and trying without homotopy failed or is not wanted,
      use the HOMOTOPY SOLVER */
   if (solveWithHomotopySolver && !nonlinsys->solved) {
-    if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
+    if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY] && !kinsol)
       warningStreamPrint(LOG_ASSERT, 0, "Failed to solve the initial system %d without homotopy method.", sysNumber);
     data->simulationInfo->lambda = 0.0;
     if (data->callback->useHomotopy == 3) {
@@ -974,7 +985,14 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
       infoStreamPrint(LOG_INIT, 1, "homotopy process\n---------------------------");
       infoStreamPrint(LOG_INIT, 0, "solve lambda0-system");
       nonlinsys->homotopySupport = 0;
-      nonlinsys->solved = solveNLS(data, threadData, sysNumber);
+      if (!kinsol) {
+        nonlinsys->solved = solveNLS(data, threadData, sysNumber);
+      } else {
+        nlsLs = data->simulationInfo->nlsLinearSolver;
+        data->simulationInfo->nlsLinearSolver = NLS_LS_LAPACK;
+        nonlinsys->solved = solveWithInitHomotopy(data, threadData, sysNumber);
+        data->simulationInfo->nlsLinearSolver = nlsLs;
+      }
       nonlinsys->homotopySupport = 1;
       infoStreamPrint(LOG_INIT, 0, "solving lambda0-system done with%s success\n---------------------------", nonlinsys->solved ? "" : " no");
       messageClose(LOG_INIT);