#2178

- bug fixes for non-linear system. - make jacobian calculation continuous - added reset for non-contnuous solving case git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@16038 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
OpenModelica · May 14, 2013 · 9571a8a · 9571a8a
1 parent b1cec89
commit 9571a8a
Showing 1 changed file with 32 additions and 0 deletions.
diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c b/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c
@@ -79,6 +79,7 @@ typedef struct DATA_HYBRD
   double* wa2;
   double* wa3;
   double* wa4;
+
 } DATA_HYBRD;
 
 static int wrapper_fvec_hybrj(integer* n, double* x, double* f, double* fjac, integer* ldjac, integer* iflag, void* data);
@@ -325,6 +326,7 @@ static int wrapper_fvec_hybrj(integer* n, double* x, double* f, double* fjac, in
   int i, currentSys = ((DATA*)data)->simulationInfo.currentNonlinearSystemIndex;
   NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo.nonlinearSystemData[currentSys]);
   DATA_HYBRD* solverData = (DATA_HYBRD*)(systemData->solverData);
+  int continuous = ((DATA*)data)->simulationInfo.solveContinuous;
 
   /* debug output */
   INFO1(LOG_NLS_V, "# Call wrapper function with mode %d", (int)*iflag);
@@ -360,11 +362,22 @@ static int wrapper_fvec_hybrj(integer* n, double* x, double* f, double* fjac, in
       break;
 
     case 2:
+
+      /* set residual function continuous for jacobian calculation */
+      if(continuous)
+        ((DATA*)data)->simulationInfo.solveContinuous = 0;
+
       /* call apropreated jacobain function */
       if(systemData->jacobianIndex != -1)
         getAnalyticalJacobian(data, fjac);
       else
         getNumericalJacobian(data, fjac, x, f);
+
+      /* reset residual function again */
+      if(continuous)
+        ((DATA*)data)->simulationInfo.solveContinuous = 1;
+
+
       break;
 
     default:
@@ -400,12 +413,16 @@ int solveHybrd(DATA *data, int sysNumber)
   double initial_factor = solverData->factor;
   int nfunc_evals = 0;
   int continuous = 1;
+  int nonContinuousCase = 0;
 
   int giveUp = 0;
   int retries = 0;
   int retries2 = 0;
   int retries3 = 0;
 
+  modelica_boolean* relationsPreBackup;
+  relationsPreBackup = (modelica_boolean*) malloc(data->modelData.nRelations*sizeof(modelica_boolean));
+
   /* debug output */
   if(ACTIVE_STREAM(LOG_NLS))
   {
@@ -513,6 +530,7 @@ int solveHybrd(DATA *data, int sysNumber)
       storeRelations(data);
     }
 
+
     /* scaling residual vector */
     {
       int l=0;
@@ -566,6 +584,15 @@ int solveHybrd(DATA *data, int sysNumber)
     if(solverData->info == 1 && (xerror > local_tol && xerror_scaled > local_tol))
       solverData->info = 4;
 
+    /* reset non-contunuousCase */
+    if (nonContinuousCase && (solverData->info == 4 || solverData->info == 5)){
+      /* set x vector */
+      memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
+
+      memcpy(data->simulationInfo.relationsPre, relationsPreBackup, sizeof(modelica_boolean)*data->modelData.nRelations);
+      nonContinuousCase = 0;
+    }
+
     /* solution found */
     if(solverData->info == 1 || xerror <= local_tol || xerror_scaled <= local_tol)
     {
@@ -686,6 +713,9 @@ int solveHybrd(DATA *data, int sysNumber)
       /* try to solve a discontinuous system */
       continuous = 0;
 
+      nonContinuousCase = 1;
+      memcpy(relationsPreBackup, data->simulationInfo.relationsPre, sizeof(modelica_boolean)*data->modelData.nRelations);
+
       giveUp = 0;
       nfunc_evals += solverData->nfev;
       if(ACTIVE_STREAM(LOG_NLS)) {
@@ -853,5 +883,7 @@ int solveHybrd(DATA *data, int sysNumber)
   solverData->factor = initial_factor;
   solverData->mode = 1;
 
+  free(relationsPreBackup);
+
   return success;
 }