Try without homotopy also for local homotopy method

Author: ptaeuber

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

@@ -185,7 +185,8 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
   FILE *pFile = NULL;
   long i;
   MODEL_DATA *mData = data->modelData;
-  int solvedWithHomotopy = 0;
+
+  data->simulationInfo->homotopyUsed = 0;
 
 #if !defined(OMC_NDELAY_EXPRESSIONS) || OMC_NDELAY_EXPRESSIONS>0
   /* initial sample and delay before initial the system */
@@ -221,7 +222,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
   MMC_CATCH_INTERNAL(simulationJumpBuffer)
 #endif
     if(init_lambda_steps > 0)
-      warningStreamPrint(LOG_ASSERT, 0, "Failed to solve the initial system without homotopy method. If homotopy is available the homotopy method is used now.");
+      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
@@ -236,6 +237,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
     if(ACTIVE_STREAM(LOG_INIT))
     {
       sprintf(buffer, "%s_global_homotopy.csv", mData->modelFilePrefix);
+      infoStreamPrint(LOG_INIT, 0, "The homotopy path will be exported to %s.", buffer);
       pFile = fopen(buffer, "wt");
       fprintf(pFile, "\"sep=,\"\n%s", "lambda");
       for(i=0; i<mData->nVariablesReal; ++i)
@@ -276,7 +278,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
           check_mixed_solutions(data, 0))
         break;
     }
-    solvedWithHomotopy = 1;
+    data->simulationInfo->homotopyUsed = 1;
     messageClose(LOG_INIT);
   }
 
@@ -291,7 +293,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
   retVal = data->callback->functionRemovedInitialEquations(data, threadData);
 
   if (!retVal)
-    infoStreamPrint(LOG_STDOUT, 0, "The initialization finished successfully %s homotopy method.", solvedWithHomotopy ? "with" : "without");
+    infoStreamPrint(LOG_STDOUT, 0, "The initialization finished successfully %s homotopy method.", data->simulationInfo->homotopyUsed ? "with" : "without");
 
   TRACE_POP
   return retVal;

SimulationRuntime/c/simulation/solver/nonlinearSystem.c

@@ -44,6 +44,7 @@
 #include "newtonIteration.h"
 #endif
 #include "nonlinearSolverHomotopy.h"
+#include "simulation/options.h"
 #include "simulation/simulation_info_json.h"
 #include "simulation/simulation_runtime.h"
 #include "simulation/solver/model_help.h"
@@ -817,54 +818,67 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
      (data->callback->useHomotopy == 0) && init_lambda_steps > 1)
     lambda_steps = init_lambda_steps;
 
-#if !defined(OMC_NO_FILESYSTEM)
-  if(lambda_steps > 1 && ACTIVE_STREAM(LOG_INIT))
-  {
-    sprintf(buffer, "%s_homotopy_nls_%d.csv", data->modelData->modelFilePrefix, sysNumber);
-    infoStreamPrint(LOG_INIT, 0, "The homotopy path of system %d will be exported to %s.", sysNumber, buffer);
-    pFile = fopen(buffer, "wt");
-    fprintf(pFile, "\"sep=,\"\n%s", "lambda");
-    for(j=0; j<nonlinsys->size; ++j)
-      fprintf(pFile, ",%s", modelInfoGetEquation(&data->modelData->modelDataXml, nonlinsys->equationIndex).vars[j]);
-    fprintf(pFile, "\n");
+  nonlinsys->solved = 0;
+  if (lambda_steps == 1 || !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY]) {
+
+    if(data->callback->useHomotopy == 0)
+      data->simulationInfo->lambda = 1.0;
+
+    if (lambda_steps > 1)
+      infoStreamPrint(LOG_INIT, 0, "Try to solve initial system %d without homotopy first.", sysNumber);
+
+    /* SOLVE! */
+    nonlinsys->solved = solveNLS(data, threadData, sysNumber);
   }
-#endif
 
-  for(step=0; step<lambda_steps-1; ++step)
-  {
-    data->simulationInfo->lambda = ((double)step)/(lambda_steps-1);
-    infoStreamPrint(LOG_INIT, 0, "[system %d] homotopy parameter lambda = %g", sysNumber, data->simulationInfo->lambda);
-    solveNLS(data, threadData, sysNumber);
+  if (lambda_steps > 1 && !nonlinsys->solved) {
+    if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
+      warningStreamPrint(LOG_ASSERT, 0, "Failed to solve initial system %d without homotopy method. The homotopy method is used now.", sysNumber);
 
 #if !defined(OMC_NO_FILESYSTEM)
     if(ACTIVE_STREAM(LOG_INIT))
     {
-      infoStreamPrint(LOG_INIT, 0, "[system %d] homotopy parameter lambda = %g done\n---------------------------", sysNumber, data->simulationInfo->lambda);
-      fprintf(pFile, "%.16g", data->simulationInfo->lambda);
+      sprintf(buffer, "%s_homotopy_nls_%d.csv", data->modelData->modelFilePrefix, sysNumber);
+      infoStreamPrint(LOG_INIT, 0, "The homotopy path of system %d will be exported to %s.", sysNumber, buffer);
+      pFile = fopen(buffer, "wt");
+      fprintf(pFile, "\"sep=,\"\n%s", "lambda");
       for(j=0; j<nonlinsys->size; ++j)
-        fprintf(pFile, ",%.16g", nonlinsys->nlsx[j]);
+        fprintf(pFile, ",%s", modelInfoGetEquation(&data->modelData->modelDataXml, nonlinsys->equationIndex).vars[j]);
       fprintf(pFile, "\n");
     }
 #endif
-  }
 
-  if(data->callback->useHomotopy == 0)
-    data->simulationInfo->lambda = 1.0;
+    for(step=0; step<lambda_steps; ++step)
+    {
+      data->simulationInfo->lambda = ((double)step)/(lambda_steps-1);
+
+      if(data->simulationInfo->lambda > 1.0) {
+          data->simulationInfo->lambda = 1.0;
+        }
+
+      infoStreamPrint(LOG_INIT, 0, "[system %d] homotopy parameter lambda = %g", sysNumber, data->simulationInfo->lambda);
+      nonlinsys->solved = solveNLS(data, threadData, sysNumber);
 
-  /* SOLVE! */
-  nonlinsys->solved = solveNLS(data, threadData, sysNumber);
+#if !defined(OMC_NO_FILESYSTEM)
+      if(ACTIVE_STREAM(LOG_INIT))
+      {
+        infoStreamPrint(LOG_INIT, 0, "[system %d] homotopy parameter lambda = %g done\n---------------------------", sysNumber, data->simulationInfo->lambda);
+        fprintf(pFile, "%.16g", data->simulationInfo->lambda);
+        for(j=0; j<nonlinsys->size; ++j)
+          fprintf(pFile, ",%.16g", nonlinsys->nlsx[j]);
+        fprintf(pFile, "\n");
+      }
+#endif
+    }
 
 #if !defined(OMC_NO_FILESYSTEM)
-  if(lambda_steps > 1 && ACTIVE_STREAM(LOG_INIT))
-  {
-    infoStreamPrint(LOG_INIT, 0, "[system %d] homotopy parameter lambda = %g done\n---------------------------", sysNumber, data->simulationInfo->lambda);
-    fprintf(pFile, "%.16g", data->simulationInfo->lambda);
-    for(j=0; j<nonlinsys->size; ++j)
-      fprintf(pFile, ",%.16g", nonlinsys->nlsx[j]);
-    fprintf(pFile, "\n");
-    fclose(pFile);
-  }
+    if(lambda_steps > 1 && ACTIVE_STREAM(LOG_INIT))
+    {
+      fclose(pFile);
+    }
 #endif
+    data->simulationInfo->homotopyUsed = 1;
+  }
 
   /* handle asserts */
   threadData->currentErrorStage = saveJumpState;

SimulationRuntime/c/simulation_data.h

@@ -590,6 +590,7 @@ typedef struct SIMULATION_INFO
   int jacobianEvals;                   /* number of different columns to evaluate functionODE */
   int currentJacobianEval;             /* current column to evaluate functionODE for Jacobian*/
 
+  int homotopyUsed;                    /* =1 the initialization problem was solved with a homotopy method, =0 otherwise */
   double lambda;                       /* homotopy parameter E [0, 1.0] */
 
   /* indicators for simulations state */