Skip to content

Commit

Permalink
added symbolic implicit euler. Improvements will come!!
Browse files Browse the repository at this point in the history 
git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@25721 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Vitalij Ruge committed Apr 23, 2015
1 parent 681f96c commit 5147a69
Show file tree
Hide file tree
Showing 13 changed files with 280 additions and 9 deletions.
160 changes: 160 additions & 0 deletions Compiler/BackEnd/BackendDAEOptimize.mo
View file
Expand Up @@ -4103,6 +4103,166 @@ algorithm
end addedScaledVarsWork;


// =============================================================================
// 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;
public function symEulerInit
"
fix the difference quotient for initial equations[0/0]
ToDo: remove me
"
input BackendDAE.BackendDAE inDAE;
output BackendDAE.BackendDAE outDAE;
algorithm
outDAE := if Flags.getConfigBool(Flags.SYM_EULER) then symEulerWork(inDAE, false) else inDAE;

end symEulerInit;

protected function symEulerWork
input BackendDAE.BackendDAE inDAE;
input Boolean b " true => add, false => remove euler equation";
output BackendDAE.BackendDAE outDAE;
protected
list<BackendDAE.EqSystem> systlst, osystlst = {};
BackendDAE.EqSystem syst_;
BackendDAE.Shared shared;
BackendDAE.Var tmpv;
DAE.ComponentRef cref;
algorithm
BackendDAE.DAE(systlst, shared) := inDAE;

// make dt
cref := ComponentReference.makeCrefIdent("$TMP$OMC$DT", 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.addKnVarDAE(tmpv, shared);

for syst in systlst 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, numberOfElement, size, arrSize;

Option<BackendDAE.IncidenceMatrix> m, mT;
BackendDAE.Matching matching;
BackendDAE.StateSets stateSets;
BackendDAE.BaseClockPartitionKind partitionKind;
list<DAE.ComponentRef> crlst;

algorithm
BackendDAE.EQSYSTEM(vars, eqns, m, mT, matching, stateSets, partitionKind) := iSyst;
BackendDAE.EQUATION_ARRAY(size, numberOfElement, arrSize, 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);
if b then // der(x)
oShared := symEulerDerVars(crlst, oShared);
end if;
oSyst := BackendDAE.EQSYSTEM(vars, eqns, NONE(), NONE(), BackendDAE.NO_MATCHING(), stateSets, partitionKind);
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;
algorithm
kind := if b then BackendDAE.VARIABLE() else BackendDAE.STATE(0,NONE());
for cref in crlst loop
(_, idx) := BackendVariable.getVar2(cref, ovars);
ovars := BackendVariable.setVarKindForVar(idx, kind, ovars);
end for;
end symEulerState;

protected function symEulerDerVars
input list<DAE.ComponentRef> crlst;
input BackendDAE.Shared shared;
output BackendDAE.Shared outShared = shared;
protected
DAE.ComponentRef cref;
BackendDAE.Var v;
algorithm
for cr in crlst loop
cref := ComponentReference.crefPrefixDer(cr);
v := BackendVariable.makeVar(cref);
outShared := BackendVariable.addKnVarDAE(v, outShared);
end for;
end symEulerDerVars;

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
2 changes: 2 additions & 0 deletions Compiler/BackEnd/BackendDAEUtil.mo
View file
Expand Up @@ -7717,6 +7717,8 @@ algorithm
(InlineArrayEquations.inlineArrayEqn, "inlineArrayEqn", false),
(BackendDAEOptimize.removeUnusedParameter, "removeUnusedParameter", false),
(BackendDAEOptimize.removeUnusedVariables, "removeUnusedVariables", false),
(BackendDAEOptimize.symEuler, "symEuler", false),
(BackendDAEOptimize.symEulerInit, "symEulerInit", false),
(SymbolicJacobian.constantLinearSystem, "constantLinearSystem", false),
(OnRelaxation.relaxSystem, "relaxSystem", false),
(BackendDAEOptimize.countOperations, "countOperations", false),
Expand Down
9 changes: 9 additions & 0 deletions Compiler/BackEnd/Differentiate.mo
View file
Expand Up @@ -1583,6 +1583,15 @@ algorithm
(exp_1, _) = Expression.makeZeroExpression(Expression.arrayDimension(tp));
then
(exp_1, inFuncs);

/*der(x)/dt*/
case ("$_DF$DER",_,_,_,_,_)
equation
(exp_1, funcs) = differentiateExp(exp, inDiffwrtCref, inInputData,inDiffType,inFuncs);
exp_2 = Expression.crefExp(ComponentReference.makeCrefIdent("$TMP$OMC$DT", DAE.T_REAL_DEFAULT, {}));
then
(Expression.expDiv(exp_1,exp_2), inFuncs);

end match;
end differentiateCallExp1Arg;

Expand Down
24 changes: 24 additions & 0 deletions Compiler/BackEnd/ExpressionSolve.mo
View file
Expand Up @@ -586,6 +586,7 @@ preprocessing for solve1,
con := true;
iter := 0;

x := unifyFunCalls(x, inExp3);
while con and iter < 1000 loop

(x, y, con) := preprocessingSolve2(x,y, inExp3);
Expand Down Expand Up @@ -1198,6 +1199,17 @@ protected function unifyFunCallsWork
e = DAE.IFEXP(DAE.RELATION(e1,DAE.GREATEREQ(tp), Expression.makeConstZero(tp),-1,NONE()),Expression.expMul(e1,e2), Expression.expMul(e1,e3));
then (e,true, iT);

// df_der(x) = (x-old(x))/dt
case(DAE.CALL(path = Absyn.IDENT(name = "$_DF$DER"),expLst = {e1}),X)
guard expHasCref(e1, X)
equation
tp = Expression.typeof(e1);
e2 = Expression.crefExp(ComponentReference.makeCrefIdent("$TMP$OMC$DT", DAE.T_REAL_DEFAULT, {}));
e3 = Expression.makePureBuiltinCall("$_old", {e1}, tp);
e3 = Expression.expSub(e1,e3);
e = Expression.expDiv(e3,e2);
then (e,true, iT);

else (inExp, true, iT);
end matchcontinue;

Expand Down Expand Up @@ -1598,6 +1610,18 @@ algorithm

then(e1, lhs, true, eqnForNewVars_, newVarsCrefs_, idepth + 1);

// $_DF$DER(x) =y -> (x-old(x))/dt = y -> x = y*dt + old(x)
case(DAE.CALL(path = Absyn.IDENT(name = "$_DF$DER"),expLst = {e1}), _)
equation
true = expHasCref(e1, inExp3);
false = expHasCref(inExp2, inExp3);
tp = Expression.typeof(e1);
e2 = Expression.crefExp(ComponentReference.makeCrefIdent("$TMP$OMC$DT", DAE.T_REAL_DEFAULT, {}));
lhs = Expression.makePureBuiltinCall("$_old", {e1}, tp);
lhs = Expression.expAdd(Expression.expMul(inExp2,e2), lhs);
then(e1, lhs, true, ieqnForNewVars, inewVarsCrefs, idepth + 1);


//QE
// a*x^n + b*x^m = c
// a*x^n - b*x^m = c
Expand Down
2 changes: 1 addition & 1 deletion Compiler/BackEnd/Initialization.mo
View file
Expand Up @@ -201,7 +201,7 @@ algorithm
initdae := BackendDAEUtil.mapEqSystem(initdae, solveInitialSystemEqSystem);

// transform and optimize DAE
pastOptModules := BackendDAEUtil.getPostOptModules(SOME({"constantLinearSystem", /* here we need a special case and remove only alias and constant (no variables of the system) variables "removeSimpleEquations", */ "tearingSystem","calculateStrongComponentJacobians", "solveSimpleEquations"}));
pastOptModules := BackendDAEUtil.getPostOptModules(SOME({"constantLinearSystem", /* here we need a special case and remove only alias and constant (no variables of the system) variables "removeSimpleEquations", */ "symEulerInit","tearingSystem","calculateStrongComponentJacobians", "solveSimpleEquations"}));
matchingAlgorithm := BackendDAEUtil.getMatchingAlgorithm(NONE());
daeHandler := BackendDAEUtil.getIndexReductionMethod(NONE());

Expand Down
3 changes: 2 additions & 1 deletion Compiler/FrontEnd/Expression.mo
View file
Expand Up @@ -6007,7 +6007,8 @@ algorithm
then (inExp,false,inTpl);
case (DAE.CALL(path = Absyn.IDENT(name = "$_round")), _)
then (inExp,false,inTpl);

case (DAE.CALL(path = Absyn.IDENT(name = "$_old")), _)
then (inExp,false,inTpl);

case (DAE.CREF(componentRef = cr1), (cr,false))
equation
Expand Down
11 changes: 7 additions & 4 deletions Compiler/FrontEnd/ExpressionSimplify.mo
View file
Expand Up @@ -1131,14 +1131,17 @@ algorithm

// smooth of constant expression
case DAE.CALL(path=Absyn.IDENT("smooth"),expLst={_,e1})
equation
true = Expression.isConst(e1);
guard Expression.isConst(e1)
then e1;

// df_der(const) --> 0
case DAE.CALL(path=Absyn.IDENT("$_DF$DER"),expLst={e1})
guard Expression.isConst(e1)
then Expression.makeConstZeroE(e1);

// delay of constant expression
case DAE.CALL(path=Absyn.IDENT("delay"),expLst={e1,_,_})
equation
true = Expression.isConst(e1);
guard Expression.isConst(e1)
then e1;

// delay of constant subexpression
Expand Down
35 changes: 34 additions & 1 deletion Compiler/Template/CodegenC.tpl
View file
Expand Up @@ -603,6 +603,8 @@ template simulationFile(SimCode simCode, String guid)
<%functionDAE(allEquations, whenClauses, modelNamePrefixStr)%>
<%functionSymEuler(modelInfo, modelNamePrefixStr)%>
<%functionODE(odeEquations,(match simulationSettingsOpt case SOME(settings as SIMULATION_SETTINGS(__)) then settings.method else ""), hpcomData.odeSchedule, modelNamePrefixStr)%>
/* forward the main in the simulation runtime */
Expand Down Expand Up @@ -694,7 +696,8 @@ template simulationFile(SimCode simCode, String guid)
<%symbolName(modelNamePrefixStr,"lagrange")%>,
<%symbolName(modelNamePrefixStr,"pickUpBoundsForInputsInOptimization")%>,
<%symbolName(modelNamePrefixStr,"setInputData")%>,
<%symbolName(modelNamePrefixStr,"getTimeGrid")%>
<%symbolName(modelNamePrefixStr,"getTimeGrid")%>,
<%symbolName(modelNamePrefixStr,"symEulerUpdate")%>
<%\n%>
};

Expand Down Expand Up @@ -1294,6 +1297,29 @@ template functionInput(ModelInfo modelInfo, String modelNamePrefix)
end match
end functionInput;

template functionSymEuler(ModelInfo modelInfo, String modelNamePrefix)
"Generates function in simulation file."
::=
match modelInfo
case MODELINFO(vars=SIMVARS(__)) then
<<
int <%symbolName(modelNamePrefix,"symEulerUpdate")%>(DATA *data, modelica_real dt)
{
TRACE_PUSH
#ifdef $P$TMP$OMC$DT
$P$TMP$OMC$DT = dt;
#else
return -1;
#endif
TRACE_POP
return 0;
}

>>
end match
end functionSymEuler;

template functionOutput(ModelInfo modelInfo, String modelNamePrefix)
"Generates function in simulation file."
::=
Expand Down Expand Up @@ -8590,10 +8616,17 @@ template daeExpCall(Exp call, Context context, Text &preExp, Text &varDecls, Tex
case CALL(path=IDENT(name="$_initialGuess"), expLst={arg as CREF(__)}) then
let namestr = cref(arg.componentRef)
'( <%namestr%>)' //
case CALL(path=IDENT(name="$_old"), expLst={arg as CREF(__)}) then
let namestr = cref(arg.componentRef)
'( _<%namestr%>(1))' //
// if arg >= 0 then 1 else -1
case CALL(path=IDENT(name="$_signNoNull"), expLst={e1}) then
let var1 = daeExp(e1, context, &preExp, &varDecls, &auxFunction)
'(<%var1%> >= 0.0 ? 1.0:-1.0)'
// numerical der()
case CALL(path=IDENT(name="$_DF$DER"), expLst={arg as CREF(__)}) then
let namestr = cref(arg.componentRef)
'(initial() ? 0 : (($P$TMP$OMC$DT == 0.0) ? ($P$DER<%namestr%>) : ($P$DER<%namestr%> = ((<%namestr%> - _<%namestr%>(1))/$P$TMP$OMC$DT))))'
// round
case CALL(path=IDENT(name="$_round"), expLst={e1}) then
let var1 = daeExp(e1, context, &preExp, &varDecls, &auxFunction)
Expand Down
8 changes: 7 additions & 1 deletion Compiler/Util/Flags.mo
View file
Expand Up @@ -743,6 +743,7 @@ constant ConfigFlag POST_OPT_MODULES = CONFIG_FLAG(16, "postOptModules",
"solveLinearSystem",
"addScaledVars",
"removeSimpleEquations",
"symEuler",
"encapsulateWhenConditions", // must called after remove simple equations
"reshufflePost",
"reduceDynamicOptimization", // before tearing
Expand Down Expand Up @@ -1045,6 +1046,10 @@ constant ConfigFlag DYNAMIC_TEARING = CONFIG_FLAG(68, "dynamicTearing",
NONE(), EXTERNAL(), BOOL_FLAG(false), NONE(),
Util.gettext("Activates dynamic tearing (TearingSet can be changed automatically during runtime)"));

constant ConfigFlag SYM_EULER = CONFIG_FLAG(69, "symEuler",
NONE(), EXTERNAL(), BOOL_FLAG(false), NONE(),
Util.gettext("Rewritte the ode system for inplicit euler."));

protected
// This is a list of all configuration flags. A flag can not be used unless it's
// in this list, and the list is checked at initialization so that all flags are
Expand Down Expand Up @@ -1117,7 +1122,8 @@ constant list<ConfigFlag> allConfigFlags = {
MAX_SIZE_FOR_SOLVE_LINIEAR_SYSTEM,
CPP_FLAGS,
REMOVE_SIMPLE_EQUATIONS,
DYNAMIC_TEARING
DYNAMIC_TEARING,
SYM_EULER
};

public function new
Expand Down
4 changes: 4 additions & 0 deletions SimulationRuntime/c/openmodelica_func.h
View file
Expand Up @@ -287,7 +287,11 @@ int (*setInputData)(DATA* data, const modelica_boolean file);
int (*getTimeGrid)(DATA *data, modelica_integer * nsi, modelica_real**t);


/*
* update parameter for symEuler
*/

int (*symEulerUpdate)(DATA * data, modelica_real dt);

};

Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/c/simulation/simulation_runtime.cpp
View file
Expand Up @@ -659,7 +659,7 @@ int callSolver(DATA* simData, string init_initMethod, string init_file,
/* if no states are present, than we can
* use euler method, since it does nothing.
*/
if (simData->modelData.nStates < 1 && solverID != S_OPTIMIZATION) {
if (simData->modelData.nStates < 1 && solverID != S_OPTIMIZATION && solverID != S_SYM_EULER) {
solverID = S_EULER;
}

Expand Down

0 comments on commit 5147a69

Please sign in to comment.