- updated non-linear systems solving algorithm

  - changed extrapolation data
  - added nominal values for scaling
  - adjusted the strategy of solving
- marked some more ThermoSysPro examples as working


git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@12061 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

Compiler/susan_codegen/SimCode/CodegenC.tpl

@@ -285,7 +285,21 @@ template functionInitializeDataStruc2(ModelInfo modelInfo, list<SimEqSystem> all
 ::=
 match modelInfo
 case MODELINFO(varInfo=VARINFO(__)) then
+   /*
+   This fragment produces some empty lines, since emptyCount is not implemented in susan!
+   */
+   let eqnsDefines = (allEquations |> eq hasindex i0 => match eq
+    case SES_ALGORITHM(__)
+    case SES_WHEN(__)
+    case SES_RESIDUAL(__)
+    case SES_ARRAY_CALL_ASSIGN(__)
+    case SES_SIMPLE_ASSIGN(__)  then ''
+    case SES_MIXED(__)
+    case SES_LINEAR(__)
+    case SES_NONLINEAR(__) then '#define SIM_PROF_EQ_<%index%> <%i0%>'; empty; separator="\n")
   <<
+  /* Some empty lines, since emptyCount is not implemented in susan! */
+  <%eqnsDefines%>
   void setupDataStruc2(DATA *data)
   {  
     <%globalDataFunctionInfoArray("function_names", functions)%>
@@ -400,14 +414,14 @@ template globalDataParDefine(SimVar simVar, String arrayName)
   case SIMVAR(arrayCref=SOME(c),aliasvar=NOALIAS()) then
     <<
     #define <%cref(c)%> data->simulationInfo.<%arrayName%>[<%index%>]
-    #define $P$START<%cref(name)%> data->modelData.<%arrayName%>Data[<%index%>].attribute.start
+    #define $P$ATTRIBUTE<%cref(name)%> data->modelData.<%arrayName%>Data[<%index%>].attribute
     #define <%cref(name)%> data->simulationInfo.<%arrayName%>[<%index%>]
     #define <%cref(name)%>__varInfo data->modelData.<%arrayName%>Data[<%index%>].info
     >>
   case SIMVAR(aliasvar=NOALIAS()) then
     <<
     #define <%cref(name)%> data->simulationInfo.<%arrayName%>[<%index%>]
-    #define $P$START<%cref(name)%> data->modelData.<%arrayName%>Data[<%index%>].attribute.start
+    #define $P$ATTRIBUTE<%cref(name)%> data->modelData.<%arrayName%>Data[<%index%>].attribute
     #define <%cref(name)%>__varInfo data->modelData.<%arrayName%>Data[<%index%>].info
     >>  
 end globalDataParDefine;
@@ -426,7 +440,7 @@ template globalDataVarDefine(SimVar simVar, String arrayName, Integer offset)
     #define <%cref(name)%> _<%cref(name)%>(0)
     #define $P$PRE<%cref(c)%> data->simulationInfo.<%arrayName%>Pre[<%intAdd(offset,index)%>]
     #define $P$PRE<%cref(name)%> data->simulationInfo.<%arrayName%>Pre[<%intAdd(offset,index)%>]
-    #define $P$START<%cref(name)%> data->modelData.<%arrayName%>Data[<%intAdd(offset,index)%>].attribute.start
+    #define $P$ATTRIBUTE<%cref(name)%> data->modelData.<%arrayName%>Data[<%intAdd(offset,index)%>].attribute
     #define <%cref(name)%>__varInfo data->modelData.<%arrayName%>Data[<%intAdd(offset,index)%>].info
     >>
   case SIMVAR(aliasvar=NOALIAS()) then
@@ -435,7 +449,7 @@ template globalDataVarDefine(SimVar simVar, String arrayName, Integer offset)
     #define _<%cref(name)%>(i) data->localData[i]-><%arrayName%>[<%intAdd(offset,index)%>]
     #define <%cref(name)%> _<%cref(name)%>(0)
     #define $P$PRE<%cref(name)%> data->simulationInfo.<%arrayName%>Pre[<%intAdd(offset,index)%>]
-    #define $P$START<%cref(name)%> data->modelData.<%arrayName%>Data[<%intAdd(offset,index)%>].attribute.start
+    #define $P$ATTRIBUTE<%cref(name)%> data->modelData.<%arrayName%>Data[<%intAdd(offset,index)%>].attribute
     #define <%cref(name)%>__varInfo data->modelData.<%arrayName%>Data[<%intAdd(offset,index)%>].info
     >>  
 end globalDataVarDefine;
@@ -675,7 +689,7 @@ template functionUpdateBoundStartValues(list<SimEqSystem> startValueEquations)
     DEBUG_INFO(LOG_INIT, "updating start-values:");
     <%startValueEquations |> SES_SIMPLE_ASSIGN(__) =>
     'DEBUG_INFO_AL2(LOG_INIT, "   %s(start=%f)", <%cref(cref)%>__varInfo.name, (<%crefType(cref)%>) <%cref(cref)%>);
-$P$START<%cref(cref)%> = <%cref(cref)%>;'
+$P$ATTRIBUTE<%cref(cref)%>.start = <%cref(cref)%>;'
     ;separator="\n"%>
   
     return 0;
@@ -1935,7 +1949,8 @@ case SES_NONLINEAR(__) then
   <%crefs |> name hasindex i0 =>
     let namestr = cref(name)
     <<
-    nls_x[<%i0%>] = extraPolate(<%namestr%>, _<%namestr%>(1) /*old*/, _<%namestr%>(2) /*old2*/);
+    nls_diag[<%i0%>] = $P$ATTRIBUTE<%namestr%>.nominal;
+    nls_xEx[<%i0%>] = nls_x[<%i0%>] = extraPolate(<%namestr%>, _<%namestr%>(1) /*old*/, _<%namestr%>(2) /*old2*/);
     nls_xold[<%i0%>] = _<%namestr%>(1) /*old*/;
     >>
   ;separator="\n"%>
@@ -5833,11 +5848,11 @@ template daeExpCallStart(Exp exp, Context context, Text &preExp /*BUFP*/,
 ::=
   match exp
   case cr as CREF(__) then
-    '$P$START<%cref(cr.componentRef)%>'
+    '$P$ATTRIBUTE<%cref(cr.componentRef)%>.start'
   case ASUB(exp = cr as CREF(__), sub = {sub_exp}) then
     let offset = daeExp(sub_exp, context, &preExp /*BUFC*/, &varDecls /*BUFD*/)
     let cref = cref(cr.componentRef)
-    '*(&$P$START<%cref%> + <%offset%>)'
+    '*(&$P$ATTRIBUTE<%cref(cr.componentRef)%>.start + <%offset%>)'
   else
     error(sourceInfo(), 'Code generation does not support start(<%printExpStr(exp)%>)')
 end daeExpCallStart;
@@ -6890,12 +6905,6 @@ template equationInfo(list<SimEqSystem> eqs)
     <<
     <%preBuf%>
     <%res%>
-    <% eqs |> eq hasindex i0 => match eq
-      case SES_MIXED(__)
-      case SES_LINEAR(__)
-      case SES_NONLINEAR(__) then '#define SIM_PROF_EQ_<%index%> <%i0%>'
-    ; separator="\n"
-    %>
     const int n_omc_equationInfo_reverse_prof_index = 0<% eqs |> eq hasindex i0 => match eq
         case SES_MIXED(__)
         case SES_LINEAR(__)

SimulationRuntime/c/math-support/matrix.h

@@ -131,36 +131,77 @@ void * _omc_hybrj_(void(*) (int*, double*, double*, double *, int*, int*, void*
    int giveUp = 0; \
    int retries = 0; \
    int retries2 = 0; \
+   int retries3 = 0; \
+   int i; \
+   for (i = 0; i < n; i++) { nls_diagSave[i] = nls_diag[i]; }; \
+   if (DEBUG_FLAG(LOG_NONLIN_SYS)){ \
+     INFO2("Start solving Non-Linear System %s at time %f", no.name, data->localData[0]->timeValue); \
+   } \
    while(!giveUp) { \
      giveUp = 1; \
      _omc_hybrd_(residual,&n, nls_x,nls_fvec,&xtol,&maxfev,&ml,&mu,&epsfcn, \
           nls_diag,&mode,&factor,&nprint,&info,&nfev,nls_fjac,&ldfjac, \
         nls_r,&lr,nls_qtf,nls_wa1,nls_wa2,nls_wa3,nls_wa4, userdata); \
       if (info == 0) \
           printErrorEqSyst(IMPROPER_INPUT, no, data->localData[0]->timeValue); \
+      if (info == 1){ \
+        if (DEBUG_FLAG(LOG_NONLIN_SYS)) { \
+          INFO_AL("### System solved! ###"); \
+          INFO_AL2("\tSolution with %d retries and %d restarts.",retries, retries2); \
+          INFO_AL2("\tinfo = %d\tnfunc = %d",info, nfev); \
+          if (DEBUG_FLAG(LOG_DEBUG)) { \
+            for (i = 0; i < n; i++) { \
+              INFO_AL7("%d. scale-factor[%d] = %f\tresidual[%d] = %f\tx[%d] = %f",i,i,nls_diag[i],i,nls_fvec[i],i,nls_x[i]); \
+            } \
+          } \
+        } \
+      } \
       if ((info == 4 || info == 5) && retries < 3) { /* first try to decrease factor*/ \
         retries++;  giveUp = 0; \
         factor = factor / 10.0; \
         if (DEBUG_FLAG(LOG_NONLIN_SYS))  \
-          printErrorEqSyst(NO_PROGRESS_FACTOR, no, factor); \
+          INFO_AL1(" - iteration making no progress:\tdecrease factor to %f",factor); \
       } else if ((info == 4 || info == 5) && retries < 5) { /* Then, try with different starting point*/  \
-        int i = 0; \
         for (i = 0; i < n; i++) { nls_x[i] += 0.1; }; \
         retries++;  giveUp=0; \
         if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
-          printErrorEqSyst(NO_PROGRESS_START_POINT, no, 1e-6); \
+          INFO_AL(" - iteration making no progress:\tvary initial point by +1\%"); \
+      } else if ((info == 4 || info == 5) && retries < 7) {   \
+         for (i = 0; i < n; i++) { nls_x[i] = nls_xEx[i]*1.01; }; \
+         retries++; giveUp=0; \
+         if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
+           INFO_AL(" - iteration making no progress:\tvary initial point by adding 1%"); \
+      } else if ((info == 4 || info == 5) && retries < 9) { \
+        for (i = 0; i < n; i++) { nls_x[i] = nls_xEx[i] * 0.99; }; \
+        retries++; giveUp=0; \
+        if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
+          INFO_AL(" - iteration making no progress:\tvary initial point by %"); \
       } else if ((info == 4 || info == 5) && retries2 < 1) { /*Then try with old values (instead of extrapolating )*/ \
-        retries = 0; retries2++; giveUp = 0; \
-        int i = 0; \
+        factor = initial_factor; retries = 0; retries2++; giveUp = 0; \
         for (i = 0; i < n; i++) { nls_x[i] = nls_xold[i]; } \
+        if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
+          INFO_AL(" - iteration making no progress:\tuse old values instead extrapolated"); \
+      } else if ((info == 4 || info == 5) && retries3 < 1) { \
+         for (i = 0; i < n; i++) { nls_diag[i] = nls_diagSave[i];} \
+         factor = initial_factor; retries = 0; retries2 = 0; mode = 2; retries3++; giveUp = 0;\
+         if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
+           INFO_AL(" - iteration making no progress:\tchange scaling factors"); \
+     } else if ((info == 4 || info == 5) && retries3 < 2) { \
+        for (i = 0; i < n; i++) { nls_x[i] = nls_xEx[i]; nls_diag[i] = fmax(1e-2,fabs(nls_xEx[i]));} \
+        factor = initial_factor; retries = 0; retries2 = 0; mode = 1; retries3++; giveUp = 0;\
+        if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
+          INFO_AL(" - iteration making no progress:\tchange scaling factors"); \
+      } else if ((info == 4 || info == 5) && retries3 < 3) {  \
+        for (i = 0; i < n; i++) { nls_x[i] = nls_xEx[i]; nls_diag[i] = 1.0; } \
+        factor = initial_factor; retries = 0; retries2 = 0; retries3++; mode = 2; giveUp = 0;\
+        if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
+        INFO_AL(" - iteration making no progress:\tremove scaling factor at all!"); \
       } else if (info >= 2 && info <= 5) { \
-        int i = 0; \
         modelErrorCode=ERROR_NONLINSYS; \
         printErrorEqSyst(ERROR_AT_TIME, no, data->localData[0]->timeValue); \
-        if (DEBUG_FLAG(LOG_NONLIN_SYS)) { \
+        if (DEBUG_FLAG(LOG_DEBUG)) { \
           for (i = 0; i < n; i++) { \
-             DEBUG_INFO_AL2(LOG_NONLIN_SYS," residual[%d] = %f",i,nls_fvec[i]); \
-             DEBUG_INFO_AL2(LOG_NONLIN_SYS," x[%d] = %f",i,nls_x[i]); \
+             INFO_AL7("\t%d. scale-factor[%d] = %f\tresidual[%d] = %f\tx[%d] = %f",i,i,nls_diag[i],i,nls_fvec[i],i,nls_x[i]); \
           } \
         } \
       }\
@@ -239,24 +280,27 @@ _omc_dgesv_(&n,&nrhs,&A[0],&lda,ipiv,&b[0],&ldb,&info); \
 } while (0) /* (no trailing ; ) */
 
 #define start_nonlinear_system(size) { double nls_x[size]; \
+double nls_xEx[size] = {0}; \
 double nls_xold[size] = {0}; \
 double nls_fvec[size] = {0}; \
 double nls_diag[size] = {0}; \
+double nls_diagSave[size] = {0}; \
 double nls_r[(size*(size + 1) / 2)] = {0}; \
 double nls_qtf[size] = {0}; \
 double nls_wa1[size] = {0}; \
 double nls_wa2[size] = {0}; \
 double nls_wa3[size] = {0}; \
 double nls_wa4[size] = {0}; \
-double xtol = 1e-12; \
+double xtol =  1e-12; \
 double epsfcn = 1e-12; \
-int maxfev = 8000; \
+int maxfev = size*10000; \
 int n = size; \
 int ml = size - 1; \
 int mu = size - 1; \
 int mode = 1; \
 int info = 0, nfev = 0, njev = 0; \
 double factor = 100.0; \
+double initial_factor = 100.0; \
 int nprint = 0; \
 int lr = (size*(size + 1)) / 2; \
 int ldfjac = size; \

SimulationRuntime/c/openmodelica.h

@@ -46,6 +46,7 @@ extern "C" {
 #include <stdio.h>
 #include <limits.h>
 #include <assert.h>
+#include <float.h>
 
 /* adrpo: extreme windows crap! */
 #if defined(__MINGW32__) || defined(_MSC_VER)

SimulationRuntime/c/simulation/solver/dassl.c

@@ -111,6 +111,7 @@ dasrt_initial(DATA* simData, SOLVER_INFO* solverInfo, DASSL_DATA *dasslData){
   dasslData->dasslStatisticsTmp = (unsigned int*) calloc(numStatistics, sizeof(unsigned int));
   ASSERT(dasslData->dasslStatisticsTmp,"out of memory");
 
+  dasslData->info[2] = 1;
   /*********************************************************************
    *info[2] = 1;  //intermediate-output mode
    *********************************************************************
@@ -180,31 +181,50 @@ int dasrt_step(DATA* simData, SOLVER_INFO* solverInfo)
 
   /* Calculate time steps until TOUT is reached
    * (DASSL calculates beyond TOUT unless info[6] is set to 1!) */
-  do
+  tout = solverInfo->currentTime + solverInfo->currentStepSize;
+  /* Check that tout is not less than timeValue
+   * else will dassl get in trouble. If that is the case we skip the current step. */
+  if (solverInfo->currentTime - tout >= -1e-13)
   {
+    DEBUG_INFO(LOG_SOLVER, "**Desired step to small try next one");
+    DEBUG_INFO(LOG_SOLVER, "**Interpolate linear");
 
-    tout = solverInfo->currentTime + solverInfo->currentStepSize;
-    /* Check that tout is not less than timeValue
-     * else will dassl get in trouble. If that is the case we skip the current step. */
-    if (solverInfo->currentTime - tout >= -1e-13)
+    for (i = 0; i < simData->modelData.nStates; i++)
     {
-      DEBUG_INFO(LOG_SOLVER, "**Desired step to small try next one");
-      DEBUG_INFO(LOG_SOLVER, "**Interpolate linear");
+      sData->realVars[i] = sDataOld->realVars[i]
+          + stateDer[i] * solverInfo->currentStepSize;
+    }
+    sData->timeValue = tout;
+    functionODE(simData);
+    solverInfo->currentTime = tout;
+    /*TODO: interpolate states and evaluate the system again */
+    return 2;
+  }
+
+  DEBUG_INFO2(LOG_SOLVER, "**Calling DDASRT from %g to %g",
+      solverInfo->currentTime, tout);
+  do
+  {
+    DEBUG_INFO2(LOG_SOLVER, "**Start step %g to %g", solverInfo->currentTime, tout);
+    if (dasslData->idid == 1){
+      DEBUG_INFO(LOG_SOLVER, "Rotate Ringbuffer!");
+
+      /* rotate RingBuffer before step is calculated */
+      rotateRingBuffer(simData->simulationData, 1, (void**) simData->localData);
+      sData = (SIMULATION_DATA*) simData->localData[0];
+      sDataOld = (SIMULATION_DATA*) simData->localData[1];
+      stateDer = sDataOld->realVars + mData->nStates;
+      sData->timeValue = solverInfo->currentTime;
 
-      for (i = 0; i < simData->modelData.nStates; i++)
-      {
-        sData->realVars[i] = sDataOld->realVars[i]
-            + stateDer[i] * solverInfo->currentStepSize;
-      }
-      sData->timeValue = tout;
-      functionODE(simData);
-      solverInfo->currentTime = tout;
-      /*TODO: interpolate states and evaluate the system again */
-      return 2;
     }
+    /* read input vars */
+    input_function(simData);
 
-    DEBUG_INFO2(LOG_SOLVER, "**Calling DDASRT from %g to %g",
-        solverInfo->currentTime, tout);
+    if (DEBUG_FLAG(LOG_DEBUG)){
+      for (i=0; i<3;i++){
+        printAllVars(simData,i);
+      }
+    }
 
     if (jac_flag)
     {
@@ -236,25 +256,6 @@ int dasrt_step(DATA* simData, SOLVER_INFO* solverInfo)
           dasslData->rpar, dasslData->ipar, dummy_Jacobian, dummy_zeroCrossing,
           &tmpZero, NULL);
     }
-     if (DEBUG_FLAG(LOG_SOLVER))
-     {
-       INFO1("DASSL call | value of idid: %d", dasslData->idid);
-       INFO1("DASSL call | current time value: %0.4g", solverInfo->currentTime);
-       INFO1("DASSL call | current integration time value: %0.4g", dasslData->rwork[3]);
-       INFO1("DASSL call | step size H to be attempted on next step: %0.4g", dasslData->rwork[2]);
-       INFO1("DASSL call | step size used on last successful step: %0.4g", dasslData->rwork[6]);
-       INFO1("DASSL call | number of steps taken so far: %d", dasslData->iwork[10]);
-       INFO1("DASSL call | number of calls of functionODE() : %d", dasslData->iwork[11]);
-       INFO1("DASSL call | number of calculation of jacobian : %d", dasslData->iwork[12]);
-       INFO1("DASSL call | total number of convergence test failures: %d", dasslData->iwork[13]);
-       INFO1("DASSL call | total number of error test failures: %d", dasslData->iwork[14]);
-     }
-    /* save dassl stats */
-    for (i = 0; i < numStatistics; i++)
-    {
-      assert(10 + i < dasslData->liw);
-      dasslData->dasslStatisticsTmp[i] = dasslData->iwork[10 + i];
-    }
 
     if (dasslData->idid == -1)
     {
@@ -275,15 +276,35 @@ int dasrt_step(DATA* simData, SOLVER_INFO* solverInfo)
       INFO1("DASRT can't continue. time = %f", sData->timeValue);
       return 1;
     }
-    sData->timeValue = solverInfo->currentTime;
-    functionODE(simData);
-  } while (dasslData->idid == -1
-      && solverInfo->currentTime <= simData->simulationInfo.stopTime);
 
-  if (simData->simulationInfo.stopTime >= solverInfo->currentTime)
+  } while (dasslData->idid == 1 || (dasslData->idid == -1
+      && solverInfo->currentTime <= simData->simulationInfo.stopTime));
+
+  if (DEBUG_FLAG(LOG_SOLVER))
   {
-    DEBUG_INFO(LOG_SOLVER, "DDASRT finished");
+    INFO1("DASSL call | value of idid: %d", dasslData->idid);
+    INFO1("DASSL call | current time value: %0.4g", solverInfo->currentTime);
+    INFO1("DASSL call | current integration time value: %0.4g", dasslData->rwork[3]);
+    INFO1("DASSL call | step size H to be attempted on next step: %0.4g", dasslData->rwork[2]);
+    INFO1("DASSL call | step size used on last successful step: %0.4g", dasslData->rwork[6]);
+    INFO1("DASSL call | number of steps taken so far: %d", dasslData->iwork[10]);
+    INFO1("DASSL call | number of calls of functionODE() : %d", dasslData->iwork[11]);
+    INFO1("DASSL call | number of calculation of jacobian : %d", dasslData->iwork[12]);
+    INFO1("DASSL call | total number of convergence test failures: %d", dasslData->iwork[13]);
+    INFO1("DASSL call | total number of error test failures: %d", dasslData->iwork[14]);
   }
+  /* save dassl stats */
+  for (i = 0; i < numStatistics; i++)
+  {
+   assert(10 + i < dasslData->liw);
+   dasslData->dasslStatisticsTmp[i] = dasslData->iwork[10 + i];
+  }
+
+  sData->timeValue = solverInfo->currentTime;
+  functionODE(simData);
+
+  DEBUG_INFO(LOG_SOLVER, "*** Finished DDASRT step! ***");
+
   return 0;
 }
 
@@ -351,16 +372,13 @@ int functionODE_residual(double *t, double *x, double *xd, double *delta,
                     fortran_integer *ires, double *rpar, fortran_integer *ipar)
 {
   DATA* data = (DATA*)(void*)rpar;
-  /*DASSL_DATA* dasslData = (DASSL_DATA*)(void*)rpar[1];*/
+
   double timeBackup;
-  double* statesBackup;
   long i;
 
   timeBackup = data->localData[0]->timeValue;
-  statesBackup = data->localData[0]->realVars;
 
   data->localData[0]->timeValue = *t;
-  data->localData[0]->realVars = x;
   functionODE(data);
 
   /* get the difference between the temp_xd(=localData->statesDerivatives)
@@ -369,7 +387,6 @@ int functionODE_residual(double *t, double *x, double *xd, double *delta,
     delta[i] = data->localData[0]->realVars[data->modelData.nStates + i] - xd[i];
   }
 
-  data->localData[0]->realVars = statesBackup;
   data->localData[0]->timeValue = timeBackup;
 
   return 0;