Try without homotopy also for adaptive homotopy

- Adaptive homotopy is a fallback initialization method
- Use -homotopyOnFirstTry to directly activate homotopy

Belonging to [master]:
  - #2165
  - OpenModelica/OpenModelica-testsuite#844

SimulationRuntime/c/openmodelica_func.h

@@ -150,7 +150,7 @@ int (*functionInitialEquations)(DATA *data, threadData_t*);
  * 2: new global homotopy approach (adaptive lambda steps)
  * 3: new local homotopy approach (adaptive lambda steps)
  */
-const int useHomotopy;
+int useHomotopy;
 
 /*! \fn functionInitialEquations_lambda0
  *

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

@@ -187,6 +187,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
   MODEL_DATA *mData = data->modelData;
   int homotopySupport = 0;
   int solveWithGlobalHomotopy;
+  int adaptiveGlobal;
 
 #if !defined(OMC_NUM_NONLINEAR_SYSTEMS) || OMC_NUM_NONLINEAR_SYSTEMS>0
   for(i=0; i<mData->nNonLinearSystems; i++) {
@@ -196,7 +197,9 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
     }
   }
 #endif
-  solveWithGlobalHomotopy = homotopySupport && ((data->callback->useHomotopy == 1 && init_lambda_steps > 1) || data->callback->useHomotopy == 2);
+  adaptiveGlobal = data->callback->useHomotopy == 2;
+  solveWithGlobalHomotopy = homotopySupport
+                            && ((data->callback->useHomotopy == 1 && init_lambda_steps > 1) || adaptiveGlobal);
 
 #if !defined(OMC_NDELAY_EXPRESSIONS) || OMC_NDELAY_EXPRESSIONS>0
   /* initial sample and delay before initial the system */
@@ -210,21 +213,21 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
      or homotopy is disabled by runtime flag '-ils=<lambda_steps>',
      solve WITHOUT HOMOTOPY or LOCAL HOMOTOPY. */
   if (!solveWithGlobalHomotopy){
-    if (data->callback->useHomotopy == 3 && !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
-      infoStreamPrint(LOG_INIT, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the adaptive local approach yet.");
     data->simulationInfo->lambda = 1.0;
     data->callback->functionInitialEquations(data, threadData);
 
   /* If there is homotopy in the model and global homotopy is activated
      and homotopy on first try is deactivated,
      TRY TO SOLVE WITHOUT HOMOTOPY FIRST.
-     To-Do: Activate trying without homotopy first also for the adaptive global homotopy approach */
-  } else if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY] && data->callback->useHomotopy != 2) {
+     TO-DO: For the adaptive global approach, provide a separate DAE with
+     the original unmanipulated systems for trying without homotopy */
+  } else if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY]) {
     /* try */
 #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);
@@ -234,18 +237,22 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
 #ifndef OMC_EMCC
   MMC_CATCH_INTERNAL(simulationJumpBuffer)
 #endif
+    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 there is homotopy in the model and global homotopy is activated
+  /* If there is homotopy in the model and the equidistant global homotopy approach is activated
      and solving without homotopy failed or is not wanted,
      use EQUIDISTANT GLOBAL HOMOTOPY METHOD. */
   if (data->callback->useHomotopy == 1 && solveWithGlobalHomotopy)
   {
     long step;
     char buffer[4096];
 
+    infoStreamPrint(LOG_INIT, 0, "Global homotopy with equidistant step size started.");
+
 #if !defined(OMC_NO_FILESYSTEM)
     const char sep[] = ",";
     if(ACTIVE_STREAM(LOG_INIT))
@@ -260,7 +267,6 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
     }
 #endif
 
-    infoStreamPrint(LOG_INIT, 0, "Global homotopy with equidistant step size started.");
     infoStreamPrint(LOG_INIT, 1, "homotopy process\n---------------------------");
     for(step=0; step<init_lambda_steps; ++step)
     {
@@ -302,9 +308,6 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
      use ADAPTIVE GLOBAL HOMOTOPY APPROACH. */
   if (data->callback->useHomotopy == 2 && solveWithGlobalHomotopy)
   {
-    if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
-      infoStreamPrint(LOG_INIT, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the adaptive global approach yet.");
-
     infoStreamPrint(LOG_INIT, 0, "Global homotopy with adaptive step size started.");
     infoStreamPrint(LOG_INIT, 1, "homotopy process\n---------------------------");
 

SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c

@@ -2106,13 +2106,12 @@ int solveHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
   int pos;
   int rank;
   int tries = 0;
-  /* Modelica homotopy operator could be used!! */
   int runHomotopy = 0;
   int skipNewton = 0;
   int numberOfFunctionEvaluationsOld = solverData->numberOfFunctionEvaluations;
   solverData->casualTearingSet = systemData->strictTearingFunctionCall != NULL;
   int constraintViolated;
-  solverData->initHomotopy = (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3) && systemData->homotopySupport;
+  solverData->initHomotopy = systemData->initHomotopy;
 
   modelica_boolean* relationsPreBackup;
   relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));

SimulationRuntime/c/simulation/solver/nonlinearSystem.c

@@ -809,12 +809,10 @@ int solveNLS(DATA *data, threadData_t *threadData, int sysNumber)
     mixedSolverData = nonlinsys->solverData;
     nonlinsys->solverData = mixedSolverData->newtonData;
 
-
     /* try to solve the system if it is the strict set, only try to solve the casual set if the constraints are satisfied */
     if ((!casualTearingSet) || constraintsSatisfied)
       success = solveHomotopy(data, threadData, sysNumber);
 
-
     /* check if solution process was successful, if not use alternative tearing set if available (dynamic tearing)*/
     if (!success && casualTearingSet){
       debugString(LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
@@ -855,6 +853,7 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
   int step;
   int equidistantHomotopy = 0;
   int solveWithHomotopySolver = 0;
+  int homotopyDeactivated = 0;
   int j;
   struct dataSolver *solverData;
   struct dataMixedSolver *mixedSolverData;
@@ -905,10 +904,34 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
                             && nonlinsys->homotopySupport
                             && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3);
 
+  homotopyDeactivated = !data->simulationInfo->initial           // Not an initialization system
+                        || !nonlinsys->homotopySupport           // There is no homotopy in this component
+                        || (data->callback->useHomotopy == 1)    // Equidistant homotopy is performed globally in symbolic_initialization()
+                        || (data->callback->useHomotopy == 0     // Equidistant local homotopy is selected but homotopy is deactivated ...
+                            && init_lambda_steps <= 1);          // ... by the number of steps
+
   nonlinsys->solved = 0;
-  /* If the adaptive local/global homotopy approach is activated and the component contains the homotopy operator
-     use the HOMOTOPY SOLVER, otherwise use the selected solver */
-  if (solveWithHomotopySolver) {
+  nonlinsys->initHomotopy = 0;
+
+  /* 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);
+
+    data->simulationInfo->lambda = 1.0;
+    /* SOLVE! */
+    nonlinsys->solved = solveNLS(data, threadData, sysNumber);
+  }
+
+  /* The following cases are only valid for initial systems with homotopy */
+  /* **********************************************************************/
+
+  /* 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])
+      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) {
       // First solve the lambda0-system separately
@@ -931,6 +954,7 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
         solverData = nonlinsys->solverData;
         nonlinsys->solverData = solverData->initHomotopyData;
       }
+      nonlinsys->initHomotopy = 1;
       nonlinsys->solved = solveHomotopy(data, threadData, sysNumber);
       if (data->simulationInfo->nlsMethod == NLS_MIXED) {
         nonlinsys->solverData = mixedSolverData;
@@ -939,22 +963,12 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
       }
     }
   }
-  else {
-    if (!(equidistantHomotopy && omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])) {
-      if (equidistantHomotopy)
-        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);
-    }
-  }
 
   /* If equidistant local homotopy is activated and trying without homotopy failed or is not wanted,
      use EQUIDISTANT LOCAL HOMOTOPY */
-  if (!nonlinsys->solved && equidistantHomotopy) {
+  if (equidistantHomotopy && !nonlinsys->solved) {
     if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
-      warningStreamPrint(LOG_ASSERT, 0, "Failed to solve initial system %d without homotopy method. The local homotopy method with equidistant step size is used now.", sysNumber);
+      warningStreamPrint(LOG_ASSERT, 0, "Failed to solve the initial system %d without homotopy method. The local homotopy method with equidistant step size is used now.", sysNumber);
     else
       infoStreamPrint(LOG_INIT, 0, "Local homotopy with equidistant step size started for nonlinear system %d.", sysNumber);
 #if !defined(OMC_NO_FILESYSTEM)

SimulationRuntime/c/simulation_data.h

@@ -257,6 +257,7 @@ typedef struct NONLINEAR_SYSTEM_DATA
   modelica_integer equationIndex;      /* index for EQUATION_INFO */
 
   modelica_boolean homotopySupport;    /* 1 if homotopy is available, 0 otherwise */
+  modelica_boolean initHomotopy;       /* 1 if the homotopy solver should be used to solve the initial system, 0 otherwise */
   modelica_boolean mixedSystem;        /* 1 if the system contains discrete variables, 0 otherwise */
   /* attributes of iteration variables */
   modelica_real *min;