[DAEmode] improve event handling and discrete loops

Belonging to [master]:
  - #2291
  - OpenModelica/OpenModelica-testsuite#884

Author: Willi Braun

Compiler/BackEnd/DAEMode.mo

@@ -46,6 +46,7 @@ import BackendDAE;
 protected
 
 import Absyn;
+import Array;
 import BackendDAEOptimize;
 import BackendDAEUtil;
 import BackendDAEFunc;
@@ -67,6 +68,7 @@ import Flags;
 import Global;
 import Initialization;
 import List;
+import Matching;
 import StackOverflow;
 import Util;
 
@@ -196,6 +198,7 @@ uniontype TraverseEqnAryFold
     BackendDAE.Variables systemVars;
     DAE.FunctionTree functionTree;
     Boolean recursiveStrongComponentRun;
+    BackendDAE.Shared shared;
   end TRAVERSER_CREATE_DAE;
 end TraverseEqnAryFold;
 
@@ -221,7 +224,7 @@ algorithm
   systemSize := BackendDAEUtil.systemSize(syst);
   newDAEVars := BackendVariable.emptyVars();
   newDAEEquations := BackendEquation.emptyEqnsSized(systemSize);
-  travArgs := TRAVERSER_CREATE_DAE(globalDAEData, newDAEVars, newDAEEquations, syst.orderedVars, shared.functionTree, false);
+  travArgs := TRAVERSER_CREATE_DAE(globalDAEData, newDAEVars, newDAEEquations, syst.orderedVars, shared.functionTree, false, shared);
   if debug then BackendDump.printEqSystem(syst); end if;
   // for every equation create corresponding residiual variable(s)
   travArgs := BackendDAEUtil.traverseEqSystemStrongComponents(syst, traverserStrongComponents, travArgs);
@@ -269,8 +272,11 @@ algorithm
   (traverserArgs) :=
   matchcontinue(inEqns, traverserArgs.recursiveStrongComponentRun, isStateVarInvoled)
     local
+      BackendDAE.EqSystem syst;
+      BackendDAE.IncidenceMatrix adjMatrix;
+
       list<BackendDAE.Equation> newResEqns;
-      list<BackendDAE.Var> newResVars, newAuxVars;
+      list<BackendDAE.Var> newResVars, newAuxVars, discVars, contVars;
       BackendDAE.Var var;
       BackendDAE.Variables systemVars;
       BackendDAE.Var dummyVar;
@@ -282,6 +288,8 @@ algorithm
       DAE.FunctionTree funcsTree;
       list<DAE.ComponentRef> crlst;
       Boolean b1, b2;
+      list<BackendDAE.Equation> discEqns;
+      list<BackendDAE.Equation> contEqns;
 
       DAE.Algorithm alg;
       DAE.ElementSource source;
@@ -482,14 +490,26 @@ algorithm
 
     case(_, false, _)
       algorithm
-        vars := inVars;
-        for e in inEqns loop
+
+        (discVars, contVars) := List.splitOnTrue(inVars, BackendVariable.isVarDiscrete);
+        (discEqns, contEqns) := getDiscAndContEqns(inVars, inEqns, discVars, contVars, traverserArgs.shared.functionTree);
+
+        // create discrete
+        for e in discEqns loop
+          size := BackendEquation.equationSize(e);
+          newAuxVars := List.firstN(discVars, size);
+          traverserArgs := traverserStrongComponents({e}, newAuxVars, {}, {}, traverserArgs);
+          discVars := List.stripN(discVars, size);
+        end for;
+
+        // create continuous
+        for e in contEqns loop
           size := BackendEquation.equationSize(e);
-          newAuxVars := List.firstN(vars, size);
+          newAuxVars := List.firstN(contVars, size);
           traverserArgs.recursiveStrongComponentRun := true;
           traverserArgs := traverserStrongComponents({e}, newAuxVars, {}, {}, traverserArgs);
           traverserArgs.recursiveStrongComponentRun := false;
-          vars := List.stripN(vars, size);
+          contVars := List.stripN(contVars, size);
         end for;
       then
         (traverserArgs);
@@ -504,6 +524,54 @@ algorithm
   end matchcontinue;
 end traverserStrongComponents;
 
+function getDiscAndContEqns
+  input list<BackendDAE.Var> inAllVars;
+  input list<BackendDAE.Equation> inAllEqns;
+  input list<BackendDAE.Var> inDiscVars;
+  input list<BackendDAE.Var> inContVars;
+  input DAE.FunctionTree functionTree;
+  output list<BackendDAE.Equation> discEqns;
+  output list<BackendDAE.Equation> contEqns;
+protected
+  BackendDAE.EqSystem syst;
+  BackendDAE.IncidenceMatrix adjMatrix;
+
+  list<Integer> varsIndex, eqnIndex;
+  array<Integer> assignVarEqn "eqn := assignVarEqn[var]";
+  array<Integer> assignEqnVar "var := assignEqnVar[eqn]";
+
+  array<Integer> mapEqnScalarArray;
+  constant Boolean debug = false;
+algorithm
+  try
+    // create syst for a matching
+    syst := BackendDAEUtil.createEqSystem(BackendVariable.listVar1(inAllVars), BackendEquation.listEquation(inAllEqns) );
+    if debug then BackendDump.printEqSystem(syst); end if;
+    (adjMatrix, _, _, mapEqnScalarArray) := BackendDAEUtil.incidenceMatrixScalar(syst, BackendDAE.NORMAL(), SOME(functionTree));
+    if debug then BackendDump.dumpIncidenceMatrix(adjMatrix); end if;
+    (assignVarEqn, assignEqnVar, true) := Matching.RegularMatching(adjMatrix, BackendDAEUtil.systemSize(syst), BackendDAEUtil.systemSize(syst));
+    if debug then BackendDump.dumpMatching(assignVarEqn); end if;
+
+    // get discrete vars indexes and then the equations
+    varsIndex := BackendVariable.getVarIndexFromVars(inDiscVars, syst.orderedVars);
+    if debug then print("discVarsIndex: "); BackendDump.dumpIncidenceRow(varsIndex);  end if;
+    eqnIndex := List.map1(varsIndex, Array.getIndexFirst, assignVarEqn);
+    if debug then print("discEqnIndex: "); BackendDump.dumpIncidenceRow(eqnIndex);  end if;
+    eqnIndex := List.unique(list(mapEqnScalarArray[i] for i in eqnIndex));
+    discEqns := BackendEquation.getList(eqnIndex, syst.orderedEqs);
+    if debug then BackendDump.equationListString(discEqns, "Discrete Equations");  end if;
+
+    // get continuous equations
+    varsIndex := BackendVariable.getVarIndexFromVars(inContVars, syst.orderedVars);
+    eqnIndex := List.map1(varsIndex, Array.getIndexFirst, assignVarEqn);
+    eqnIndex := List.unique(list(mapEqnScalarArray[i] for i in eqnIndex));
+    if debug then print("contEqnIndex: "); BackendDump.dumpIncidenceRow(eqnIndex);  end if;
+    contEqns := BackendEquation.getList(eqnIndex, syst.orderedEqs);
+    if debug then BackendDump.equationListString(contEqns, "Continuous Equations");  end if;
+  else
+    fail();
+  end try;
+end getDiscAndContEqns;
 
 function addVarsGlobalData
   input output BackendDAE.BackendDAEModeData globalDAEData;
@@ -513,8 +581,6 @@ protected
 algorithm
   // prepare algebraic states
   vars := List.filterOnTrue(inVars, BackendVariable.isNonStateVar);
-  // check for discrete vars
-  {} := List.filterOnTrue(vars, BackendVariable.isVarDiscrete);
   vars := list(BackendVariable.setVarKind(v, BackendDAE.ALG_STATE()) for v in vars);
   //print("Alg vars: " + BackendDump.varListString(vars, "") + "\n");
   globalDAEData.algStateVars := listAppend(vars, globalDAEData.algStateVars);

SimulationRuntime/c/simulation/simulation_runtime.cpp

@@ -455,7 +455,7 @@ int startNonInteractiveSimulation(int argc, char**argv, DATA* data, threadData_t
   }
   /* if daeMode is turned on than use also the ida solver */
   if (omc_flag[FLAG_DAE_MODE] && std::string("ida") != data->simulationInfo->solverMethod) {
-    data->simulationInfo->solverMethod = std::string("ida").c_str();
+    data->simulationInfo->solverMethod = GC_strdup(std::string("ida").c_str());
     infoStreamPrint(LOG_SIMULATION, 0, "overwrite solver method: %s [DAEmode works only with IDA solver]", data->simulationInfo->solverMethod);
   }
 
@@ -644,7 +644,7 @@ static int callSolver(DATA* simData, threadData_t *threadData, string init_initM
   /* if no states are present, then we can
    * use euler method, since it does nothing.
    */
-  if (simData->modelData->nStates < 1 && solverID != S_OPTIMIZATION && solverID != S_SYM_SOLVER) {
+  if (simData->modelData->nStates < 1 && solverID != S_OPTIMIZATION && solverID != S_SYM_SOLVER && !compiledInDAEMode) {
     solverID = S_EULER;
   }
 

SimulationRuntime/c/simulation/solver/ida_solver.c

@@ -555,22 +555,21 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo
   /* configure algebraic variables as such */
   if (idaData->daeMode)
   {
-    if (!omc_flag[FLAG_NO_SUPPRESS_ALG])
+    if (omc_flag[FLAG_NO_SUPPRESS_ALG])
     {
       flag = IDASetSuppressAlg(idaData->ida_mem, TRUE);
       if (checkIDAflag(flag)){
         throwStreamPrint(threadData, "##IDA## Suppress algebraic variables in the local error test failed");
       }
+    }
+    for(i=0; i<idaData->N; ++i)
+    {
+      tmp[i] = (i<data->modelData->nStates)? 1.0: 0.0;
+    }
 
-      for(i=0; i<idaData->N; ++i)
-      {
-        tmp[i] = (i<data->modelData->nStates)? 1.0: 0.0;
-      }
-
-      flag = IDASetId(idaData->ida_mem, N_VMake_Serial(idaData->N,tmp));
-      if (checkIDAflag(flag)){
-        throwStreamPrint(threadData, "##IDA## Mark algebraic variables as such failed!");
-      }
+    flag = IDASetId(idaData->ida_mem, N_VMake_Serial(idaData->N,tmp));
+    if (checkIDAflag(flag)){
+      throwStreamPrint(threadData, "##IDA## Mark algebraic variables as such failed!");
     }
   }
 
@@ -703,6 +702,7 @@ ida_event_update(DATA* data, threadData_t *threadData)
   if (compiledInDAEMode == 3){
     if (initializedSolver){
 
+      data->simulationInfo->needToIterate = 0;
       /* get new values from data -> TODO: update all discrete */
       data->simulationInfo->discreteCall = 1;
 
@@ -731,9 +731,9 @@ ida_event_update(DATA* data, threadData_t *threadData)
       }
 
       /* increase limits of the non-linear solver */
-      IDASetMaxNumStepsIC(idaData->ida_mem, 2*data->modelData->nStates*10);
-      IDASetMaxNumJacsIC(idaData->ida_mem, 2*data->modelData->nStates*10);
-      IDASetMaxNumItersIC(idaData->ida_mem, 2*data->modelData->nStates*10);
+      IDASetMaxNumStepsIC(idaData->ida_mem, 2*idaData->N*10);
+      IDASetMaxNumJacsIC(idaData->ida_mem, 2*idaData->N*10);
+      IDASetMaxNumItersIC(idaData->ida_mem, 2*idaData->N*10);
       /* Calc Consistent y_algebraic and y_prime with current y */
       flag = IDACalcIC(idaData->ida_mem, IDA_YA_YDP_INIT, data->localData[0]->timeValue+init_h);
 
@@ -1029,11 +1029,12 @@ ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
   /* initialize states and der(states) */
   if (idaData->daeMode)
   {
-    idaData->residualFunction(data->localData[0]->timeValue, idaData->y, idaData->yp, idaData->newdelta, idaData);
     memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
     // and  also algebraic vars
     data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, idaData->states + data->modelData->nStates);
     memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
+    sData->timeValue = solverInfo->currentTime;
+    idaData->residualFunction(sData->timeValue, idaData->y , idaData->yp, idaData->newdelta, idaData);
   }
 
   /* sensitivity mode */
@@ -1149,7 +1150,6 @@ int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, voi
   {
     setContext(data, &time, CONTEXT_ODE);
   }
-  timeBackup = data->localData[0]->timeValue;
   data->localData[0]->timeValue = time;
 
   saveJumpState = threadData->currentErrorStage;
@@ -1233,8 +1233,6 @@ int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, voi
 
   threadData->currentErrorStage = saveJumpState;
 
-  data->localData[0]->timeValue = timeBackup;
-
   if (data->simulationInfo->currentContext == CONTEXT_ODE){
     unsetContext(data);
   }
@@ -1250,8 +1248,9 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void*
   IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
   DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
   threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->threadData);
+  double *states = N_VGetArrayPointer(yy);
+  double *statesDer = N_VGetArrayPointer(yp);
 
-  double timeBackup;
   int saveJumpState;
 
   infoStreamPrint(LOG_SOLVER_V, 1, "### eval rootsFunctionIDA ###");
@@ -1270,21 +1269,24 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void*
   saveJumpState = threadData->currentErrorStage;
   threadData->currentErrorStage = ERROR_EVENTSEARCH;
 
-  timeBackup = data->localData[0]->timeValue;
-  data->localData[0]->timeValue = time;
-
   if (idaData->daeMode)
   {
-    memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
-    data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, idaData->states + data->modelData->nStates);
-    memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
+    memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
+    data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, states + data->modelData->nStates);
+    memcpy(data->localData[0]->realVars + data->modelData->nStates, statesDer, sizeof(double)*data->modelData->nStates);
   }
 
+  data->localData[0]->timeValue = time;
+
   /* read input vars */
   externalInputUpdate(data);
   data->callback->input_function(data, threadData);
   /* eval needed equations*/
   if (idaData->daeMode && compiledInDAEMode == 1){}
+  else if (idaData->daeMode && compiledInDAEMode == 3)
+  {
+   idaData->residualFunction(time, yy, yp, idaData->newdelta, idaData);
+  }
   else
   {
     data->callback->function_ZeroCrossingsEquations(data, threadData);
@@ -1293,7 +1295,6 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void*
   data->callback->function_ZeroCrossings(data, threadData, gout);
 
   threadData->currentErrorStage = saveJumpState;
-  data->localData[0]->timeValue = timeBackup;
 
   /* scale data again */
   if (omc_flag[FLAG_IDA_SCALING])

SimulationRuntime/c/simulation/solver/solver_main.c

@@ -687,15 +687,19 @@ int solver_main(DATA* data, threadData_t *threadData, const char* init_initMetho
   externalInputallocate(data);
   /* set tolerance for ZeroCrossings */
   setZCtol(fmin(data->simulationInfo->stepSize, data->simulationInfo->tolerance));
-
-  omc_alloc_interface.collect_a_little();
-  /* initialize all parts of the model */
-  retVal = initializeModel(data, threadData, init_initMethod, init_file, init_time);
   omc_alloc_interface.collect_a_little();
 
-  if(0 == retVal) {
-    retVal = initializeSolverData(data, threadData, &solverInfo);
-    initSolverInfo = 1;
+  /* initialize solver data */
+  /* For the DAEmode we need to initialize solverData before the initialization,
+   * since the solver is used to obtain consistent values also via updateDiscreteSystem
+   */
+  retVal = initializeSolverData(data, threadData, &solverInfo);
+  initSolverInfo = 1;
+
+  /* initialize all parts of the model */
+  if (0 == retVal){
+    retVal = initializeModel(data, threadData, init_initMethod, init_file, init_time);
+    omc_alloc_interface.collect_a_little();
   }
 
 #if !defined(OMC_MINIMAL_RUNTIME)