added improved symbolical inline solver

- symbolic inline integration is working by replacing the der operator with the
  difference quotient(forward/backward)
- solver supports explicit and implicit order 1 integration methods
  enabled by compiler flag --symSolver=(impEuler|expEuler)
- added solver methods to the runtime -s symSolver|symSolverSsc without and with step size control

Author: kbalzereit

Compiler/BackEnd/BackendDAE.mo

@@ -111,7 +111,6 @@ uniontype Shared "Data shared for all equation-systems"
     Variables externalObjects               "External object variables";
     Variables aliasVars                     "Data originating from removed simple equations needed to build
                                              variables' lookup table (in C output).
-
                                              In that way, double buffering of variables in pre()-buffer, extrapolation
                                              buffer and results caching, etc., is avoided, but in C-code output all the
                                              data about variables' names, comments, units, etc. is preserved as well as
@@ -132,6 +131,13 @@ uniontype Shared "Data shared for all equation-systems"
   end SHARED;
 end Shared;
 
+uniontype InlineData
+  record INLINE_DATA
+    EqSystems inlineSystems;
+    Variables knownVariables;
+  end INLINE_DATA;
+end InlineData;
+
 uniontype BasePartition
   record BASE_PARTITION
     .DAE.ClockKind clock;
@@ -170,6 +176,7 @@ uniontype BackendDAEType "BackendDAEType to indicate different types of BackendD
   record ARRAYSYSTEM     "Type for multi dim equation arrays BackendDAE.DAE" end ARRAYSYSTEM;
   record PARAMETERSYSTEM "Type for parameter system BackendDAE.DAE"          end PARAMETERSYSTEM;
   record INITIALSYSTEM   "Type for initial system BackendDAE.DAE"            end INITIALSYSTEM;
+  record INLINESYSTEM    "Type for inline system BackendDAE.DAE"             end INLINESYSTEM;
 end BackendDAEType;
 
 //
@@ -262,9 +269,8 @@ uniontype VarKind "variable kind"
   record OPT_LOOP_INPUT
     .DAE.ComponentRef replaceExp;
   end OPT_LOOP_INPUT;
-  record ALG_STATE "algebraic state"
-    VarKind oldKind;
-  end ALG_STATE;
+  record ALG_STATE  end ALG_STATE; // algebraic state used by inline solver
+  record ALG_STATE_OLD  end ALG_STATE_OLD; // algebraic state old value used by inline solver
   record DAE_RESIDUAL_VAR end DAE_RESIDUAL_VAR; // variable kind used for DAEmode
 end VarKind;
 
@@ -715,7 +721,7 @@ public constant String functionDerivativeNamePrefix = "$funDER";
 
 public constant String optimizationMayerTermName = "$OMC$objectMayerTerm";
 public constant String optimizationLagrangeTermName = "$OMC$objectLagrangeTerm";
-public constant String symEulerDT = "__OMC_DT";
+public constant String symSolverDT = "__OMC_DT";
 
 type FullJacobian = Option<list<tuple<Integer, Integer, Equation>>>;
 

Compiler/BackEnd/BackendDAECreate.mo

@@ -158,7 +158,8 @@ algorithm
                                                     extObjCls,
                                                     BackendDAE.SIMULATION(),
                                                     symjacs,inExtraInfo,
-                                                    BackendDAEUtil.emptyPartitionsInfo()));
+                                                    BackendDAEUtil.emptyPartitionsInfo()
+                                                    ));
   BackendDAEUtil.checkBackendDAEWithErrorMsg(outBackendDAE);
   varSize := BackendVariable.varsSize(vars_1);
   eqnSize := BackendDAEUtil.equationSize(eqnarr);

Compiler/BackEnd/BackendDAEOptimize.mo

@@ -5171,144 +5171,6 @@ algorithm
 
 end simplifyLoops_SplitFactors;
 
-
-// =============================================================================
-// section for symEuler
-//
-// replace der(x) with difference quotient
-// -->  implicit euler
-//
-// after removeSimpliEquation
-// before tearing
-// ToDo: not add to initial equation
-// author: Vitalij Ruge
-// =============================================================================
-
-public function symEuler
- input BackendDAE.BackendDAE inDAE;
- output BackendDAE.BackendDAE outDAE;
-algorithm
- outDAE := if Flags.getConfigBool(Flags.SYM_EULER) then symEulerWork(inDAE, true) else inDAE;
-end symEuler;
-
-protected function symEulerWork
-  input BackendDAE.BackendDAE inDAE;
-  input Boolean b " true => add, false => remove euler equation";
-  output BackendDAE.BackendDAE outDAE;
-protected
-  list<BackendDAE.EqSystem> osystlst = {};
-  BackendDAE.EqSystem syst_;
-  BackendDAE.Shared shared;
-  BackendDAE.Var tmpv;
-  DAE.ComponentRef cref;
-algorithm
-  // make dt
-  cref := ComponentReference.makeCrefIdent(BackendDAE.symEulerDT, DAE.T_REAL_DEFAULT, {});
-  tmpv := BackendVariable.makeVar(cref);
-  //tmpv := BackendVariable.setVarKind(tmpv, BackendDAE.PARAM());
-  tmpv := BackendVariable.setBindExp(tmpv, SOME(DAE.RCONST(0.0)));
-  shared := BackendVariable.addGlobalKnownVarDAE(tmpv, inDAE.shared);
-
-  for syst in inDAE.eqs loop
-   (syst_, shared) := symEulerUpdateSyst(syst, b, shared);
-   osystlst := syst_ :: osystlst;
-  end for;
-
-  outDAE := BackendDAE.DAE(osystlst, shared);
-  //BackendDump.bltdump("BackendDAEOptimize.removeSymEulerEquation", outDAE);
-end symEulerWork;
-
-protected function symEulerUpdateSyst
-  input BackendDAE.EqSystem iSyst;
-  input Boolean b;
-  input BackendDAE.Shared shared;
-  output BackendDAE.EqSystem oSyst;
-  output BackendDAE.Shared oShared = shared;
-protected
-  array<Option<BackendDAE.Equation>> equOptArr;
-  Option<BackendDAE.Equation> oeqn;
-  BackendDAE.Equation eqn;
-  BackendDAE.Variables vars;
-  BackendDAE.EquationArray eqns;
-  Integer n;
-  list<DAE.ComponentRef> crlst;
-
-algorithm
-  oSyst := match iSyst
-    local
-      BackendDAE.EqSystem syst;
-    case syst as BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns)
-      algorithm
-        BackendDAE.EQUATION_ARRAY(equOptArr=equOptArr) := eqns;
-        n := arrayLength(equOptArr);
-        crlst := {};
-        for i in 1:n loop
-          oeqn := arrayGet(equOptArr, i);
-          if isSome(oeqn) then
-            SOME(eqn) := oeqn;
-            (eqn, (_,crlst)) := BackendEquation.traverseExpsOfEquation(eqn, symEulerUpdateEqn, (b,crlst));
-            arrayUpdate(equOptArr, i, SOME(eqn));
-          end if;
-        end for;
-        // states -> vars
-        vars := symEulerState(vars, crlst, b);
-        syst.orderedVars := vars;
-        syst.orderedEqs := eqns;
-      then BackendDAEUtil.clearEqSyst(syst);
-  end match;
-end symEulerUpdateSyst;
-
-protected function symEulerState
-  input BackendDAE.Variables vars;
-  input list<DAE.ComponentRef> crlst;
-  input Boolean b;
-  output BackendDAE.Variables ovars = vars;
-
-protected
-  Integer idx;
-  BackendDAE.VarKind kind, oldKind;
-algorithm
-  for cref in crlst loop
-    (_, idx) := BackendVariable.getVar2(cref, ovars);
-    oldKind := BackendVariable.getVarKindForVar(idx,ovars);
-    if b then
-      kind := BackendDAE.ALG_STATE(oldKind);
-    else
-     BackendDAE.ALG_STATE(kind) := oldKind;
-    end if;
-    ovars :=  BackendVariable.setVarKindForVar(idx, kind, ovars);
-  end for;
-end symEulerState;
-
-protected function symEulerUpdateEqn
-  input DAE.Exp inExp;
-  input tuple<Boolean, list<DAE.ComponentRef>> inTpl;
-  output DAE.Exp outExp;
-  output tuple<Boolean, list<DAE.ComponentRef>> outTpl;
-algorithm
-  (outExp, outTpl) := Expression.traverseExpBottomUp(inExp, symEulerUpdateDer, inTpl);
-end symEulerUpdateEqn;
-
-protected function symEulerUpdateDer
-  input DAE.Exp inExp;
-  input tuple<Boolean, list<DAE.ComponentRef>> inTpl;
-  output DAE.Exp outExp;
-  output tuple<Boolean, list<DAE.ComponentRef>> outTpl;
-algorithm
-  (outExp, outTpl) := match (inTpl, inExp)
-            local DAE.Exp exp; DAE.Type tp; list<DAE.ComponentRef> cr_lst; DAE.ComponentRef cr;
-
-            case ((true,cr_lst), DAE.CALL(path=Absyn.IDENT(name="der"), expLst={exp as DAE.CREF(ty=tp, componentRef = cr)}))
-            then (Expression.makePureBuiltinCall("$_DF$DER", {exp}, tp), (true,List.unionElt(cr,cr_lst)));
-
-            case ((false,cr_lst), DAE.CALL(path=Absyn.IDENT(name="$_DF$DER"), expLst={exp as DAE.CREF(ty=tp, componentRef = cr)}))
-            then (Expression.makePureBuiltinCall("der", {exp}, tp),(false, List.unionElt(cr,cr_lst)));
-
-            else (inExp, inTpl);
-            end match;
-end symEulerUpdateDer;
-
-
 // =============================================================================
 // section for introduceDerAlias
 //

Compiler/BackEnd/BackendDAEUtil.mo

@@ -104,6 +104,7 @@ import ResolveLoops;
 import SCode;
 import Sorting;
 import StateMachineFeatures;
+import SymbolicImplicitSolver;
 import SymbolicJacobian;
 import SynchronousFeatures;
 import System;
@@ -486,6 +487,15 @@ algorithm
   outSystem := BackendDAE.EQSYSTEM(vars, eqns, m, mt, matching, inSystem.stateSets, inSystem.partitionKind, removedEqs);
 end copyEqSystem;
 
+public function copyEqSystems
+  input BackendDAE.EqSystems inSystems;
+  output BackendDAE.EqSystems outSystems = {};
+algorithm
+  for e in inSystems loop
+    outSystems := copyEqSystem(e) :: outSystems;
+  end for;
+end copyEqSystems;
+
 public function mergeEqSystems
   input BackendDAE.EqSystem System1;
   input output BackendDAE.EqSystem System2;
@@ -1305,7 +1315,7 @@ protected
   BackendDAE.Variables v;
 algorithm
   BackendDAE.EQSYSTEM(orderedVars = v,m=SOME(m)) := syst;
-  if Flags.getConfigBool(Flags.SYM_EULER) then
+  if (Flags.getConfigEnum(Flags.SYM_SOLVER) > 0) then
     (_,statevarindx_lst) := BackendVariable.getAllAlgStateVarIndexFromVariables(v);
   else
     (_,statevarindx_lst) := BackendVariable.getAllStateVarIndexFromVariables(v);
@@ -1675,7 +1685,11 @@ algorithm
     case syst as BackendDAE.EQSYSTEM( orderedEqs=ordererdEqs, orderedVars=v,
                                       matching=BackendDAE.MATCHING(ass1=ass1, ass2=ass2) )
       algorithm
-        (_, statevarindx_lst) := BackendVariable.getAllStateVarIndexFromVariables(v);
+        if (Flags.getConfigEnum(Flags.SYM_SOLVER) > 0) then
+          (_,statevarindx_lst) := BackendVariable.getAllAlgStateVarIndexFromVariables(v);
+        else
+          (_,statevarindx_lst) := BackendVariable.getAllStateVarIndexFromVariables(v);
+        end if;
         indx_lst_v := BackendVariable.getVarIndexFromVariables(iVars, v);
 
         indx_lst_v := listAppend(indx_lst_v, statevarindx_lst) "overestimate";
@@ -2981,6 +2995,9 @@ algorithm
     case (BackendDAE.VAR(varKind = BackendDAE.VARIABLE())::rest,i::irest,_,_)
       guard not AvlSetInt.hasKey(vars, i)
       then incidenceRowExp1(rest,irest,AvlSetInt.add(vars, i),diffindex);
+    case (BackendDAE.VAR(varKind = BackendDAE.ALG_STATE())::rest,i::irest,_,_)
+      guard not AvlSetInt.hasKey(vars, i)
+      then incidenceRowExp1(rest,irest,AvlSetInt.add(vars, i),diffindex);
     case (BackendDAE.VAR(varKind = BackendDAE.DISCRETE())::rest,i::irest,_,_)
       guard not AvlSetInt.hasKey(vars, i)
       then incidenceRowExp1(rest,irest,AvlSetInt.add(vars, i),diffindex);
@@ -6770,6 +6787,7 @@ public function getSolvedSystem "Run the equation system pipeline."
   input Option<list<String>> strPostOptModules = NONE();
   output BackendDAE.BackendDAE outSimDAE;
   output BackendDAE.BackendDAE outInitDAE;
+  output Option<BackendDAE.InlineData > outInlineData;
   output Boolean outUseHomotopy "true if homotopy(...) is used during initialization";
   output Option<BackendDAE.BackendDAE> outInitDAE_lambda0;
   output list<BackendDAE.Equation> outRemovedInitialEquationLst;
@@ -6781,6 +6799,7 @@ protected
   list<tuple<BackendDAEFunc.optimizationModule, String>> postOptModules;
   tuple<BackendDAEFunc.StructurallySingularSystemHandlerFunc, String, BackendDAEFunc.stateDeselectionFunc, String> daeHandler;
   tuple<BackendDAEFunc.matchingAlgorithmFunc, String> matchingAlgorithm;
+  BackendDAE.InlineData inlineData;
 algorithm
   preOptModules := getPreOptModules(strPreOptModules);
   postOptModules := getPostOptModules(strPostOptModules);
@@ -6833,6 +6852,9 @@ algorithm
   simDAE := BackendDAEOptimize.addInitialStmtsToAlgorithms(simDAE);
   simDAE := Initialization.removeInitializationStuff(simDAE);
 
+  // generate inline solver
+  outInlineData := SymbolicImplicitSolver.symSolver(simDAE);
+
   // post-optimization phase
   outSimDAE := postOptimizeDAE(simDAE, postOptModules, matchingAlgorithm, daeHandler);
 
@@ -7398,7 +7420,6 @@ protected function allPostOptimizationModules
     (BackendDAEOptimize.inlineFunctionInLoops, "forceInlineFunctionInLoops"), // before simplifyComplexFunction
     (BackendDAEOptimize.simplifyComplexFunction, "simplifyComplexFunction"),
     (ExpressionSolve.solveSimpleEquations, "solveSimpleEquations"),
-    (BackendDAEOptimize.symEuler, "symEuler"),
     (ResolveLoops.reshuffling_post, "reshufflePost"),
     (DynamicOptimization.reduceDynamicOptimization, "reduceDynamicOptimization"), // before tearing
     (Tearing.tearingSystem, "tearingSystem"),
@@ -7573,10 +7594,6 @@ algorithm
     enabledModules := deprecatedDebugFlag(Flags.ADD_SCALED_VARS, enabledModules, "addScaledVars_states", "postOptModules+");
     enabledModules := deprecatedDebugFlag(Flags.ADD_SCALED_VARS_INPUT, enabledModules, "addScaledVars_inputs", "postOptModules+");
 
-    if Flags.getConfigBool(Flags.SYM_EULER) then
-      enabledModules := "symEuler"::enabledModules;
-    end if;
-
     if Flags.getConfigInt(Flags.SIMPLIFY_LOOPS) > 0 then
       enabledModules := "simplifyLoops"::enabledModules;
     end if;
@@ -8346,7 +8363,8 @@ algorithm
                               backendDAEType,
                               {},
                               ei,
-                              emptyPartitionsInfo());
+                              emptyPartitionsInfo()
+                              );
 end createEmptyShared;
 
 public function emptyPartitionsInfo

Compiler/BackEnd/BackendDump.mo

@@ -253,6 +253,9 @@ end printClassAttributes;
 
 public function printShared "This function dumps the BackendDAE.Shared representation to stdout."
   input BackendDAE.Shared inShared;
+protected
+  BackendDAE.EqSystems eqS;
+  BackendDAE.InlineData inlineData;
 algorithm
   print("\nBackendDAEType: ");
   printBackendDAEType(inShared.backendDAEType);
@@ -279,6 +282,7 @@ algorithm
   if Flags.isSet(Flags.DUMP_FUNCTIONS) then
     DAEDump.dumpFunctionTree(inShared.functionTree, "Functions");
   end if;
+
 end printShared;
 
 public function printBasePartitions
@@ -351,6 +355,7 @@ algorithm
     case (BackendDAE.ARRAYSYSTEM()) then "multidim equation arrays";
     case (BackendDAE.PARAMETERSYSTEM()) then "parameter system";
     case (BackendDAE.INITIALSYSTEM()) then "initialization";
+    case (BackendDAE.INLINESYSTEM()) then "inline system";
   end match;
 end printBackendDAEType2String;
 
@@ -2406,6 +2411,7 @@ algorithm
     case BackendDAE.OPT_TGRID()  then "OPT_TGRID";
     case BackendDAE.OPT_LOOP_INPUT()  then "OPT_LOOP_INPUT";
     case BackendDAE.ALG_STATE()  then "ALG_STATE";
+    case BackendDAE.ALG_STATE_OLD()  then "ALG_STATE_OLD";
   end match;
 end kindString;
 

Compiler/BackEnd/BackendEquation.mo

@@ -649,7 +649,7 @@ algorithm
 end traverseExpsOfEquationListList_WithStop;
 
 public function traverseExpsOfEquation<T> "author: Frenkel TUD 2010-11
-  Traverses all expressions of a equation.
+  Traverses all expressions of an equation.
   It is possible to change the equation."
   input BackendDAE.Equation inEquation;
   input FuncExpType inFunc;

Compiler/BackEnd/BackendVariable.mo

@@ -1184,6 +1184,17 @@ algorithm
   end match;
 end isRealOptimizeDerInput;
 
+public function isAlgebraicOldState
+"Return true if variable is old algebraic variable for inline integration."
+  input BackendDAE.Var inVar;
+  output Boolean outBoolean;
+algorithm
+  outBoolean := match (inVar)
+    case (BackendDAE.VAR(varKind =  BackendDAE.ALG_STATE_OLD())) then true;
+    else false;
+  end match;
+end isAlgebraicOldState;
+
 public function hasMayerTermAnno
 "author: Vitalij Ruge
  Return true if variable has isMayer=true annotation"

Compiler/BackEnd/Differentiate.mo

@@ -1651,14 +1651,6 @@ algorithm
       then
        (exp_1, inFuncs);
 
-    /*der(x)/dt*/
-    case ("$_DF$DER",_)
-      equation
-        (exp_1,_) = differentiateExp(exp, inDiffwrtCref, inInputData,inDiffType,inFuncs, maxIter, expStack);
-        exp_2 =  Expression.crefExp(ComponentReference.makeCrefIdent(BackendDAE.symEulerDT, DAE.T_REAL_DEFAULT, {}));
-      then
-       (Expression.expDiv(exp_1,exp_2), inFuncs);
-
   end match;
 end differentiateCallExp1Arg;