Solve lambda0-system separately ...

... before applying adaptive local homotopy method.

Belonging to [master]:
  - #2038

Author: ptaeuber

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

@@ -218,7 +218,7 @@ static int symbolic_initialization(DATA *data, threadData_t *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 new global homotopy approach */
+     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) {
     /* try */
 #ifndef OMC_EMCC
@@ -240,7 +240,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
 
   /* If there is homotopy in the model and global homotopy is activated
      and solving without homotopy failed or is not wanted,
-     use GLOBAL HOMOTOPY METHOD. */
+     use EQUIDISTANT GLOBAL HOMOTOPY METHOD. */
   if (data->callback->useHomotopy == 1 && solveWithGlobalHomotopy)
   {
     long step;
@@ -250,7 +250,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
     const char sep[] = ",";
     if(ACTIVE_STREAM(LOG_INIT))
     {
-      sprintf(buffer, "%s_global_homotopy.csv", mData->modelFilePrefix);
+      sprintf(buffer, "%s_equidistant_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=%s\"\n%s", sep, "\"lambda\"");
@@ -297,9 +297,9 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
 #endif
   }
 
-  /* If there is homotopy in the model and the new global homotopy approach is activated
+  /* If there is homotopy in the model and the adaptive global homotopy approach is activated
      and solving without homotopy failed or is not wanted,
-     use NEW GLOBAL HOMOTOPY APPROACH. */
+     use ADAPTIVE GLOBAL HOMOTOPY APPROACH. */
   if (data->callback->useHomotopy == 2 && solveWithGlobalHomotopy)
   {
     if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
@@ -318,7 +318,6 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
     infoStreamPrint(LOG_INIT, 0, "run along the homotopy path and solve the actual system");
     data->callback->functionInitialEquations(data, threadData);
 
-    data->simulationInfo->homotopyUsed = 1;
     messageClose(LOG_INIT);
   }
 

SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c

@@ -1613,7 +1613,7 @@ static int homotopyAlgorithm(DATA_HOMOTOPY* solverData, double *x)
     const char sep[] = ",";
     if(solverData->initHomotopy && ACTIVE_STREAM(LOG_INIT))
     {
-      sprintf(buffer, "%s_syst%d_new_global_homotopy_%s.csv", solverData->data->modelData->modelFilePrefix, solverData->sysNumber, solverData->startDirection > 0 ? "pos" : "neg");
+      sprintf(buffer, "%s_nonlinsys%d_adaptive_%s_homotopy_%s.csv", solverData->data->modelData->modelFilePrefix, solverData->sysNumber, solverData->data->callback->useHomotopy == 2 ? "global" : "local", solverData->startDirection > 0 ? "pos" : "neg");
       infoStreamPrint(LOG_INIT, 0, "The homotopy path will be exported to %s.", buffer);
       pFile = fopen(buffer, "wt");
       fprintf(pFile, "\"sep=%s\"\n%s", sep, "\"lambda\"");
@@ -1792,6 +1792,7 @@ static int homotopyAlgorithm(DATA_HOMOTOPY* solverData, double *x)
     debugString(LOG_NLS_HOMOTOPY, "Newton iteration for corrector step begins!");
     for(j=0;j<maxiter;j++)
     {
+      debugInt(LOG_NLS_HOMOTOPY, "Iteration: ", j+1);
       if (vec2Norm(solverData->n, solverData->hvec)<hEps || vec2Norm(solverData->n, solverData->hvecScaled)<hEps)
       {
         debugString(LOG_NLS_HOMOTOPY, "step accepted!");
@@ -2244,24 +2245,31 @@ int solveHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
     runHomotopy++;
     /* debug output */
     debugString(LOG_NLS_HOMOTOPY, "======================================================");
-    debugInt(LOG_NLS_HOMOTOPY, "RUN HOMOTOPY",runHomotopy);
 
     if (solverData->initHomotopy) {
       if (runHomotopy == 1) {
         solverData->h_function = wrapper_fvec;
         solverData->hJac_dh = wrapper_fvec_der;
         solverData->startDirection = 1.0;
-        debugDouble(LOG_INIT,"STARTING HOMOTOPY METHOD; startDirection = ", solverData->startDirection);
+        debugInt(LOG_INIT, "Homotopy run: ", runHomotopy);
+        debugDouble(LOG_INIT,"startDirection = ", solverData->startDirection);
       }
 
       if (runHomotopy == 2) {
         solverData->h_function = wrapper_fvec;
         solverData->hJac_dh = wrapper_fvec_der;
         solverData->startDirection = -1.0;
-        debugDouble(LOG_INIT,"STARTING HOMOTOPY METHOD; startDirection = ", solverData->startDirection);
+        debugInt(LOG_INIT, "Homotopy run: ", runHomotopy);
+        debugDouble(LOG_INIT,"Try again with startDirection = ", solverData->startDirection);
+      }
+
+      if (runHomotopy == 3) {
+        success = 0;
+        break;
       }
     }
     else {
+      debugInt(LOG_NLS_HOMOTOPY, "Homotopy run: ", runHomotopy);
       if (runHomotopy == 1)
       {
         /* store x0 and calculate f(x0) -> newton homotopy, fJac(x0) -> taylor, affin homotopy */

SimulationRuntime/c/simulation/solver/nonlinearSystem.c

@@ -426,7 +426,7 @@ int initializeNonlinearSystems(DATA *data, threadData_t *threadData)
 #endif
 
     /* allocate solver data */
-    if((data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3) && nonlinsys[i].homotopySupport)
+    if(data->callback->useHomotopy == 2 && nonlinsys[i].homotopySupport)
       allocateHomotopyData(size-1, &nonlinsys[i].solverData);
     else{
       switch(data->simulationInfo->nlsMethod)
@@ -447,6 +447,8 @@ int initializeNonlinearSystems(DATA *data, threadData_t *threadData)
         break;
 #if !defined(OMC_MINIMAL_RUNTIME)
       case NLS_MIXED:
+        if (data->callback->useHomotopy == 3 && nonlinsys[i].homotopySupport)
+          size--;
         mixedSolverData = (struct dataNewtonAndHybrid*) malloc(sizeof(struct dataNewtonAndHybrid));
         allocateHomotopyData(size, &(mixedSolverData->newtonData));
 
@@ -530,7 +532,7 @@ int freeNonlinearSystems(DATA *data, threadData_t *threadData)
     }
 #endif
     /* free solver data */
-    if((data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3) && nonlinsys[i].homotopySupport)
+    if(data->callback->useHomotopy == 2 && nonlinsys[i].homotopySupport)
       freeHomotopyData(&nonlinsys[i].solverData);
     else{
       switch(data->simulationInfo->nlsMethod)
@@ -796,6 +798,7 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
   int equidistantHomotopy = 0;
   int solveWithHomotopySolver = 0;
   int j;
+  struct dataNewtonAndHybrid *mixedSolverData;
   char buffer[4096];
   FILE *pFile = NULL;
 
@@ -809,7 +812,7 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
   rt_ext_tp_tick(&nonlinsys->totalTimeClock);
 
   /* grab the initial guess */
-  infoStreamPrint(LOG_NLS_EXTRAPOLATE, 1, "############ Start new iteration for system %ld at time %g ############", nonlinsys->equationIndex, data->localData[0]->timeValue);
+  infoStreamPrint(LOG_NLS_EXTRAPOLATE, 1, "############ Start new iteration for nonlinear system %ld at time %g ############", nonlinsys->equationIndex, data->localData[0]->timeValue);
   /* if last solving is too long ago use just old values  */
   if (fabs(data->localData[0]->timeValue - nonlinsys->lastTimeSolved) < 5*data->simulationInfo->stepSize || casualTearingSet)
   {
@@ -847,11 +850,31 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
   /* 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) {
-    if (data->callback->useHomotopy == 3)
-      infoStreamPrint(LOG_INIT, 0, "Local homotopy with adaptive step size started for system %d.", sysNumber);
     data->simulationInfo->lambda = 0.0;
+    if (data->callback->useHomotopy == 3) {
+      // First solve the lambda0-system separately
+      infoStreamPrint(LOG_INIT, 0, "Local homotopy with adaptive step size started for nonlinear system %d.", sysNumber);
+      infoStreamPrint(LOG_INIT, 1, "homotopy process\n---------------------------");
+      infoStreamPrint(LOG_INIT, 0, "solve lambda0-system");
+      nonlinsys->homotopySupport = 0;
+      nonlinsys->solved = solveNLS(data, threadData, sysNumber);
+      nonlinsys->homotopySupport = 1;
+      infoStreamPrint(LOG_INIT, 0, "solving lambda0-system done\n---------------------------");
+      infoStreamPrint(LOG_INIT, 0, "run along the homotopy path and solve the actual system");
+      messageClose(LOG_INIT);
+    }
     /* SOLVE! */
-    nonlinsys->solved = solveHomotopy(data, threadData, sysNumber);
+    if (data->callback->useHomotopy == 2 || nonlinsys->solved) {
+      if (data->callback->useHomotopy == 3) {
+        mixedSolverData = nonlinsys->solverData;
+        nonlinsys->solverData = mixedSolverData->newtonData;
+      }
+      nonlinsys->solved = solveHomotopy(data, threadData, sysNumber);
+      if (data->callback->useHomotopy == 3) {
+        nonlinsys->solverData = mixedSolverData;
+      }
+    }
+    data->simulationInfo->homotopyUsed = 1;
   }
   else {
     if (!(equidistantHomotopy && omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])) {
@@ -870,12 +893,12 @@ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
     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);
     else
-      infoStreamPrint(LOG_INIT, 0, "Local homotopy with equidistant step size started for system %d.", sysNumber);
+      infoStreamPrint(LOG_INIT, 0, "Local homotopy with equidistant step size started for nonlinear system %d.", sysNumber);
 #if !defined(OMC_NO_FILESYSTEM)
     const char sep[] = ",";
     if(ACTIVE_STREAM(LOG_INIT))
     {
-      sprintf(buffer, "%s_homotopy_nls_%d.csv", data->modelData->modelFilePrefix, sysNumber);
+      sprintf(buffer, "%s_nonlinsys%d_equidistant_local_homotopy.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=%s\"\n%s", sep, "\"lambda\"");