Skip to content

Commit

Permalink
added improved symbolical inline solver
Browse files Browse the repository at this point in the history 
- 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
  • Loading branch information
@OpenModelica-Hudson
kbalzereit authored and OpenModelica-Hudson committed Feb 2, 2017
1 parent 69c53b3 commit 0dbd1e9
Show file tree
Hide file tree
Showing 41 changed files with 1,307 additions and 621 deletions.
16 changes: 11 additions & 5 deletions Compiler/BackEnd/BackendDAE.mo
View file
Expand Up @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -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;

//
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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>>>;

Expand Down
3 changes: 2 additions & 1 deletion Compiler/BackEnd/BackendDAECreate.mo
View file
Expand Up @@ -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);
Expand Down
138 changes: 0 additions & 138 deletions Compiler/BackEnd/BackendDAEOptimize.mo
View file
Expand Up @@ -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
//
Expand Down
34 changes: 26 additions & 8 deletions Compiler/BackEnd/BackendDAEUtil.mo
View file
Expand Up @@ -104,6 +104,7 @@ import ResolveLoops;
import SCode;
import Sorting;
import StateMachineFeatures;
import SymbolicImplicitSolver;
import SymbolicJacobian;
import SynchronousFeatures;
import System;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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";
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand All @@ -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);
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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"),
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -8346,7 +8363,8 @@ algorithm
backendDAEType,
{},
ei,
emptyPartitionsInfo());
emptyPartitionsInfo()
);
end createEmptyShared;

public function emptyPartitionsInfo
Expand Down
6 changes: 6 additions & 0 deletions Compiler/BackEnd/BackendDump.mo
View file
Expand Up @@ -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);
Expand All @@ -279,6 +282,7 @@ algorithm
if Flags.isSet(Flags.DUMP_FUNCTIONS) then
DAEDump.dumpFunctionTree(inShared.functionTree, "Functions");
end if;

end printShared;

public function printBasePartitions
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;

Expand Down
2 changes: 1 addition & 1 deletion Compiler/BackEnd/BackendEquation.mo
View file
Expand Up @@ -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;
Expand Down
11 changes: 11 additions & 0 deletions Compiler/BackEnd/BackendVariable.mo
View file
Expand Up @@ -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"
Expand Down
8 changes: 0 additions & 8 deletions Compiler/BackEnd/Differentiate.mo
View file
Expand Up @@ -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;

Expand Down

0 comments on commit 0dbd1e9

Please sign in to comment.