Skip to content

Commit

Permalink
[SC] remove the replacement of the DIV operator by a call
Browse files Browse the repository at this point in the history 
 - added thow execption to the optimization, since ipopt
   seems to be abel to handle it

Belonging to [master]:
  - OpenModelica/OMCompiler#2379
  - OpenModelica/OpenModelica-testsuite#927
  • Loading branch information
@OpenModelica-Hudson
Willi Braun authored and OpenModelica-Hudson committed Apr 18, 2018
1 parent de3a487 commit feb746e
Show file tree
Hide file tree
Showing 4 changed files with 33 additions and 296 deletions.
259 changes: 0 additions & 259 deletions Compiler/SimCode/SimCodeUtil.mo
View file
Expand Up @@ -449,21 +449,6 @@ algorithm
zeroCrossings := updateZeroCrossEqnIndex(zeroCrossings, eqBackendSimCodeMapping, BackendDAEUtil.equationArraySizeBDAE(dlow));
if debug then execStat("simCode: update zero crossing index"); end if;

// replace div operator with div operator with check of Division by zero
allEquations := List.map(allEquations, addDivExpErrorMsgtoSimEqSystem);
odeEquations := List.mapList(odeEquations, addDivExpErrorMsgtoSimEqSystem);
algebraicEquations := List.mapList(algebraicEquations, addDivExpErrorMsgtoSimEqSystem);
startValueEquations := List.map(startValueEquations, addDivExpErrorMsgtoSimEqSystem);
nominalValueEquations := List.map(nominalValueEquations, addDivExpErrorMsgtoSimEqSystem);
minValueEquations := List.map(minValueEquations, addDivExpErrorMsgtoSimEqSystem);
maxValueEquations := List.map(maxValueEquations, addDivExpErrorMsgtoSimEqSystem);
parameterEquations := List.map(parameterEquations, addDivExpErrorMsgtoSimEqSystem);
removedEquations := List.map(removedEquations, addDivExpErrorMsgtoSimEqSystem);
initialEquations := List.map(initialEquations, addDivExpErrorMsgtoSimEqSystem);
initialEquations_lambda0 := List.map(initialEquations_lambda0, addDivExpErrorMsgtoSimEqSystem);
removedInitialEquations := List.map(removedInitialEquations, addDivExpErrorMsgtoSimEqSystem);
if debug then execStat("simCode: addDivExpErrorMsgtoSimEqSystem"); end if;

// collect all LinearSystem and NonlinearSystem algebraic system in modelInfo and update
// the corresponding index (indexNonLinear, indexLinear) in SES_NONLINEAR and SES_LINEAR
// Also collect all jacobians
Expand Down Expand Up @@ -771,9 +756,6 @@ algorithm

(ouniqueEqIndex, removedEquations) := BackendEquation.traverseEquationArray(syst.removedEqs, traversedlowEqToSimEqSystem, (ouniqueEqIndex, {}));

equations := List.map(equations, addDivExpErrorMsgtoSimEqSystem);
removedEquations := List.map(removedEquations, addDivExpErrorMsgtoSimEqSystem);

prevClockedVars := {};
isPrevVar := arrayCreate(BackendVariable.varsSize(syst.orderedVars), false);

Expand Down Expand Up @@ -8521,247 +8503,6 @@ algorithm
(_, outExpLst) := Expression.traverseExpList(inExps, inFn, {});
end getMatchingExpsList;

protected function addDivExpErrorMsgtoExp "author: Frenkel TUD 2010-02
Adds the error msg to Expression.Div."
input DAE.Exp inExp;
input DAE.ElementSource inSource;
output DAE.Exp outExp;
algorithm
false := Expression.traverseCrefsFromExp(inExp, traversingXLOCExpFinder, false);
(outExp, _) := Expression.traverseExpBottomUp(inExp, traversingDivExpFinder, inSource);
end addDivExpErrorMsgtoExp;

protected function traversingXLOCExpFinder "author: Frenkel TUD 2010-02"
input DAE.ComponentRef inCref;
input Boolean inB;
output Boolean outB;
algorithm
outB := match inCref
case DAE.CREF_IDENT(ident="xloc", identType=DAE.T_ARRAY(dims={DAE.DIM_UNKNOWN()}))
then true;

else inB;
end match;
end traversingXLOCExpFinder;

protected function traversingDivExpFinder "author: Frenkel TUD 2010-02"
input DAE.Exp inExp;
input DAE.ElementSource inSource;
output DAE.Exp outExp;
output DAE.ElementSource outSource;
algorithm
(outExp,outSource) := matchcontinue (inExp,inSource)
local
DAE.Exp e, e1, e2;
DAE.Type ty;
String se;
DAE.ElementSource source;
case (e as DAE.BINARY(operator = DAE.DIV(_), exp2 = e2), source)
equation
true = Expression.isConst(e2);
false = Expression.isZero(e2);
then (e, source);

case (DAE.BINARY(exp1 = e1, operator = DAE.DIV(ty), exp2 = e2), source)
then (DAE.CALL(Absyn.IDENT("DIVISION"), {e1, e2}, DAE.CALL_ATTR(ty, false, true, false, false, DAE.NO_INLINE(), DAE.NO_TAIL())), source);

case (e as DAE.BINARY(operator = DAE.DIV_ARRAY_SCALAR(_), exp2 = e2), source)
equation
true = Expression.isConst(e2);
false = Expression.isZero(e2);
then (e, source);
case (DAE.BINARY(exp1 = e1, operator = DAE.DIV_ARRAY_SCALAR(ty), exp2 = e2), source)
then (DAE.CALL(Absyn.IDENT("DIVISION_ARRAY_SCALAR"), {e1, e2}, DAE.CALL_ATTR(ty, false, true, false, false, DAE.NO_INLINE(), DAE.NO_TAIL())), source);

case (e as DAE.BINARY(operator = DAE.DIV_SCALAR_ARRAY(_), exp2 = e2), source)
equation
true = Expression.isConst(e2);
false = Expression.isZero(e2);
then (e, source);
case (DAE.BINARY(exp1 = e1, operator = DAE.DIV_SCALAR_ARRAY(ty), exp2 = e2), source)
then (DAE.CALL(Absyn.IDENT("DIVISION_SCALAR_ARRAY"), {e1, e2}, DAE.CALL_ATTR(ty, false, true, false, false, DAE.NO_INLINE(), DAE.NO_TAIL())), source);
else (inExp,inSource);
end matchcontinue;
end traversingDivExpFinder;

protected function addDivExpErrorMsgtosimJac
input tuple<Integer, Integer, SimCode.SimEqSystem> inJac;
output tuple<Integer, Integer, SimCode.SimEqSystem> outJac;
protected
Integer a, b;
SimCode.SimEqSystem ses, ses_;
algorithm
((a, b, ses)) := inJac;
ses_ := addDivExpErrorMsgtoSimEqSystem(ses);
outJac := if referenceEq(ses, ses_) then inJac else (a, b, ses_);
end addDivExpErrorMsgtosimJac;

protected function addDivExpErrorMsgtosymJac "helper for addDivExpErrorMsgtoSimEqSystem."
input Option<SimCode.JacobianMatrix> inJac;
output Option<SimCode.JacobianMatrix> outJac = inJac;
algorithm
outJac := match inJac
local
list<SimCode.SimEqSystem> eqns;
SimCode.JacobianColumn jacColumn;
SimCode.JacobianMatrix tmpJac;
case SOME(tmpJac as SimCode.JAC_MATRIX(columns={jacColumn as SimCode.JAC_COLUMN()}))
equation
jacColumn.columnEqns = List.map(jacColumn.columnEqns, addDivExpErrorMsgtoSimEqSystem);
tmpJac.columns = {jacColumn};
then
SOME(tmpJac);
else
outJac;
end match;
end addDivExpErrorMsgtosymJac;

protected function addDivExpErrorMsgtoSimEqSystem "Traverses all subexpressions of an expression of an equation."
input SimCode.SimEqSystem inSES;
output SimCode.SimEqSystem outSES;
algorithm
outSES:=
matchcontinue (inSES)
local
DAE.Exp e,left, e_;
DAE.ComponentRef cr;
Boolean partOfMixed,partOfMixed1;
list<SimCodeVar.SimVar> vars,vars1;
list<DAE.Exp> elst, elst1;
list<tuple<Integer, Integer, SimCode.SimEqSystem>> simJac, simJac1;
Integer index, index1, indexSys, indexSys1, nUnknowns, nUnknowns1;
list<DAE.ComponentRef> crefs, crefs1;
SimCode.SimEqSystem cont, cont1, elseWhenEq, elseWhen;
list<SimCode.SimEqSystem> discEqs, discEqs1;
list<DAE.ComponentRef> conditions;
Boolean initialCall;
DAE.ElementSource source;
Option<SimCode.JacobianMatrix> symJac, symJac1;
list<DAE.ElementSource> sources,sources1;
Boolean homotopySupport, homotopySupport1;
Boolean mixedSystem, mixedSystem1, tornSystem, tornSystem1;
list<BackendDAE.WhenOperator> whenStmtLst, whenOps;
BackendDAE.WhenOperator whenOpNew;
Option<SimCode.SimEqSystem> oelsewe;
BackendDAE.Constraints cons;
BackendDAE.EquationAttributes eqAttr;

case SimCode.SES_RESIDUAL(index= index, exp = e, source = source, eqAttr=eqAttr)
equation
e = addDivExpErrorMsgtoExp(e, source);
then
SimCode.SES_RESIDUAL(index, e, source, eqAttr);

case SimCode.SES_SIMPLE_ASSIGN(index= index, cref = cr, exp = e, source = source, eqAttr=eqAttr)
equation
e_ = addDivExpErrorMsgtoExp(e, source);
then
if referenceEq(e,e_) then inSES else SimCode.SES_SIMPLE_ASSIGN(index, cr, e_, source, eqAttr);

case SimCode.SES_SIMPLE_ASSIGN_CONSTRAINTS(index= index, cref = cr, exp = e, source = source, cons = cons, eqAttr=eqAttr)
equation
e_ = addDivExpErrorMsgtoExp(e, source);
then
if referenceEq(e,e_) then inSES else SimCode.SES_SIMPLE_ASSIGN_CONSTRAINTS(index, cr, e_, source, cons, eqAttr);

case SimCode.SES_ARRAY_CALL_ASSIGN(index = index, lhs = left, exp = e, source = source, eqAttr=eqAttr)
equation
e = addDivExpErrorMsgtoExp(e, source);
then
SimCode.SES_ARRAY_CALL_ASSIGN(index, left, e, source, eqAttr);
/*
case (SimCode.SES_ALGORITHM(), inDlowMode)
equation
e = addDivExpErrorMsgtoExp(e, (source, inDlowMode));
then
SimCode.SES_ALGORITHM();
*/

// no dynamic tearing
case SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(index, partOfMixed, tornSystem, vars, elst, simJac, discEqs, symJac, sources, indexSys, nUnknowns), NONE(), eqAttr=eqAttr)
equation
simJac1 = List.map(simJac, addDivExpErrorMsgtosimJac);
symJac = addDivExpErrorMsgtosymJac(symJac);
elst1 = List.map1(elst, addDivExpErrorMsgtoExp, DAE.emptyElementSource);
discEqs1 = List.map(discEqs, addDivExpErrorMsgtoSimEqSystem);
then
SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(index, partOfMixed, tornSystem, vars, elst1, simJac1, discEqs1, symJac, sources, indexSys, nUnknowns), NONE(), eqAttr);

case SimCode.SES_NONLINEAR(SimCode.NONLINEARSYSTEM(index=index, eqs=discEqs, crefs=crefs, indexNonLinearSystem=indexSys, nUnknowns=nUnknowns, jacobianMatrix=symJac, homotopySupport=homotopySupport, mixedSystem=mixedSystem, tornSystem=tornSystem), NONE(), eqAttr=eqAttr)
equation
discEqs = List.map(discEqs, addDivExpErrorMsgtoSimEqSystem);
symJac = addDivExpErrorMsgtosymJac(symJac);
then
SimCode.SES_NONLINEAR(SimCode.NONLINEARSYSTEM(index, discEqs, crefs, indexSys, nUnknowns, symJac, homotopySupport, mixedSystem, tornSystem), NONE(), eqAttr);

// dynamic tearing
case SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(index, partOfMixed, tornSystem, vars, elst, simJac, discEqs, symJac, sources, indexSys, nUnknowns), SOME(SimCode.LINEARSYSTEM(index1, partOfMixed1, tornSystem1, vars1, elst1, simJac1, discEqs1, symJac1, sources1, indexSys1, nUnknowns1)), eqAttr=eqAttr)
equation
simJac = List.map(simJac, addDivExpErrorMsgtosimJac);
symJac = addDivExpErrorMsgtosymJac(symJac);
elst = List.map1(elst, addDivExpErrorMsgtoExp, DAE.emptyElementSource);
discEqs = List.map(discEqs, addDivExpErrorMsgtoSimEqSystem);
simJac1 = List.map(simJac1, addDivExpErrorMsgtosimJac);
symJac1 = addDivExpErrorMsgtosymJac(symJac1);
elst1 = List.map1(elst1, addDivExpErrorMsgtoExp, DAE.emptyElementSource);
discEqs1 = List.map(discEqs1, addDivExpErrorMsgtoSimEqSystem);
then
SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(index, partOfMixed, tornSystem, vars, elst, simJac, discEqs, symJac, sources, indexSys, nUnknowns), SOME(SimCode.LINEARSYSTEM(index1, partOfMixed1, tornSystem1, vars1, elst1, simJac1, discEqs1, symJac1, sources1, indexSys1, nUnknowns1)), eqAttr);

case SimCode.SES_NONLINEAR(SimCode.NONLINEARSYSTEM(index=index, eqs=discEqs, crefs=crefs, indexNonLinearSystem=indexSys, nUnknowns=nUnknowns, jacobianMatrix=symJac, homotopySupport=homotopySupport, mixedSystem=mixedSystem, tornSystem=tornSystem), SOME(SimCode.NONLINEARSYSTEM(index=index1, eqs=discEqs1, crefs=crefs1, indexNonLinearSystem=indexSys1, nUnknowns=nUnknowns1, jacobianMatrix=symJac1, homotopySupport=homotopySupport1, mixedSystem=mixedSystem1, tornSystem=tornSystem1)), eqAttr=eqAttr)
equation
discEqs = List.map(discEqs, addDivExpErrorMsgtoSimEqSystem);
symJac = addDivExpErrorMsgtosymJac(symJac);
discEqs1 = List.map(discEqs1, addDivExpErrorMsgtoSimEqSystem);
symJac1 = addDivExpErrorMsgtosymJac(symJac1);
then
SimCode.SES_NONLINEAR(SimCode.NONLINEARSYSTEM(index, discEqs, crefs, indexSys, nUnknowns, symJac, homotopySupport, mixedSystem, tornSystem), SOME(SimCode.NONLINEARSYSTEM(index1, discEqs1, crefs1, indexSys1, nUnknowns1, symJac1, homotopySupport1, mixedSystem1, tornSystem1)), eqAttr);

case SimCode.SES_MIXED(index, cont, vars, discEqs, indexSys, eqAttr=eqAttr)
equation
cont1 = addDivExpErrorMsgtoSimEqSystem(cont);
discEqs1 = List.map(discEqs, addDivExpErrorMsgtoSimEqSystem);
then
SimCode.SES_MIXED(index, cont1, vars, discEqs1, indexSys, eqAttr);

case SimCode.SES_WHEN(index=index, conditions=conditions, initialCall=initialCall, whenStmtLst=whenStmtLst, elseWhen=oelsewe, source=source, eqAttr=eqAttr)
algorithm
whenOps := {};
for whenOp in whenStmtLst loop
whenOpNew := match whenOp
case BackendDAE.ASSIGN(left, e, source) then BackendDAE.ASSIGN(left, addDivExpErrorMsgtoExp(e, source), source);
else whenOp;
end match;
whenOps := whenOpNew::whenOps;
end for;
if isSome(oelsewe) then
SOME(elseWhen) := oelsewe;
elseWhenEq := addDivExpErrorMsgtoSimEqSystem(elseWhen);
oelsewe := SOME(elseWhenEq);
else
oelsewe := NONE();
end if;
then
SimCode.SES_WHEN(index, conditions, initialCall, whenOps, oelsewe, source, eqAttr);
else inSES;
end matchcontinue;
end addDivExpErrorMsgtoSimEqSystem;

protected function addDivExpErrorMsgtoSimEqSystemTuple
input tuple<SimCode.SimEqSystem,Integer> inSES;
output tuple<SimCode.SimEqSystem, Integer> outSES;

protected
Integer sccIdx;
SimCode.SimEqSystem eqSyst;

algorithm
(eqSyst,sccIdx) := inSES;
eqSyst := addDivExpErrorMsgtoSimEqSystem(eqSyst);
outSES := (eqSyst,sccIdx);

end addDivExpErrorMsgtoSimEqSystemTuple;

protected function extractVarUnit "author: asodja, 2010-03-11
Extract variable's unit and displayUnit as strings from
DAE.VariablesAttributes structures."
Expand Down
59 changes: 24 additions & 35 deletions Compiler/Template/CodegenCFunctions.tpl
View file
Expand Up @@ -4899,13 +4899,21 @@ case BINARY(__) then
else '<%e1%> - (<%e2%>)'
case MUL(__) then '(<%e1%>) * (<%e2%>)'
case DIV(ty = ty) then
let tvar = tempDecl(expTypeModelica(ty),&varDecls)
let &preExp += '<%tvar%> = <%e2%>;<%\n%>'
let &preExp +=
if acceptMetaModelicaGrammar()
then 'if (<%tvar%> == 0) {<%generateThrow()%>;}<%\n%>'
else 'if (<%tvar%> == 0) {throwStreamPrint(threadData, "Division by zero %s", "<%Util.escapeModelicaStringToCString(ExpressionDumpTpl.dumpExp(exp,"\""))%>");}<%\n%>'
'(<%e1%>) / <%tvar%>'
(match context
case FUNCTION_CONTEXT(__)
case PARALLEL_FUNCTION_CONTEXT(__) then
let tvar = tempDecl(expTypeModelica(ty),&varDecls)
let &preExp += '<%tvar%> = <%e2%>;<%\n%>'
let &preExp += if acceptMetaModelicaGrammar() then 'if (<%tvar%> == 0) {<%generateThrow()%>;}<%\n%>'
'(<%e1%>) / <%tvar%>'
case SIMULATION_CONTEXT() then
let e2str = Util.escapeModelicaStringToCString(ExpressionDumpTpl.dumpExp(exp2,"\""))
'DIVISION_SIM(<%e1%>,<%e2%>,"<%e2str%>",equationIndexes)'
else
let e2str = Util.escapeModelicaStringToCString(ExpressionDumpTpl.dumpExp(exp2,"\""))
'DIVISION(<%e1%>,<%e2%>,"<%e2str%>")'
)

case POW(__) then
if isHalf(exp2) then
let tmp = tempDecl(expTypeFromExpModelica(exp1),&varDecls)
Expand Down Expand Up @@ -5054,7 +5062,15 @@ case BINARY(__) then
let type = match ty case T_ARRAY(ty=T_INTEGER(__)) then "integer_array"
case T_ARRAY(ty=T_ENUMERATION(__)) then "integer_array"
else "real_array"
'div_alloc_<%type%>_scalar(<%e1%>, <%e2%>)'
let e2str = Util.escapeModelicaStringToCString(ExpressionDumpTpl.dumpExp(exp2,"\""))
(match context
case FUNCTION_CONTEXT(__)
case PARALLEL_FUNCTION_CONTEXT(__) then
'div_alloc_<%type%>_scalar(<%e1%>, <%e2%>)'
else
'division_alloc_<%type%>_scalar(threadData,<%e1%>,<%e2%>,"<%e2str%>")'
)

case DIV_SCALAR_ARRAY(__) then
let type = match ty case T_ARRAY(ty = T_INTEGER(__)) then "integer_array"
case T_ARRAY(ty = T_ENUMERATION(__)) then "integer_array"
Expand Down Expand Up @@ -5501,33 +5517,6 @@ template daeExpCall(Exp call, Context context, Text &preExp, Text &varDecls, Tex
let var2 = daeExp(e2, context, &preExp, &varDecls, &auxFunction)
var2

case CALL(path=IDENT(name="DIVISION"),
expLst={e1, e2}) then
let var1 = daeExp(e1, context, &preExp, &varDecls, &auxFunction)
let var2 = daeExp(e2, context, &preExp, &varDecls, &auxFunction)
let var3 = Util.escapeModelicaStringToCString(ExpressionDumpTpl.dumpExp(e2,"\""))
(match context
case FUNCTION_CONTEXT(__)
case PARALLEL_FUNCTION_CONTEXT(__) then
'DIVISION(<%var1%>,<%var2%>,"<%var3%>")'
else
'DIVISION_SIM(<%var1%>,<%var2%>,"<%var3%>",equationIndexes)'
)

case CALL(attr=CALL_ATTR(ty=ty),
path=IDENT(name="DIVISION_ARRAY_SCALAR"),
expLst={e1, e2}) then
let type = match ty case T_ARRAY(ty=T_INTEGER(__)) then "integer_array"
case T_ARRAY(ty=T_ENUMERATION(__)) then "integer_array"
else "real_array"
let var1 = daeExp(e1, context, &preExp, &varDecls, &auxFunction)
let var2 = daeExp(e2, context, &preExp, &varDecls, &auxFunction)
let var3 = Util.escapeModelicaStringToCString(ExpressionDumpTpl.dumpExp(e2,"\""))
'division_alloc_<%type%>_scalar(threadData,<%var1%>,<%var2%>,"<%var3%>")'

case exp as CALL(attr=CALL_ATTR(ty=ty), path=IDENT(name="DIVISION_ARRAY_SCALAR")) then
error(sourceInfo(),'Code generation does not support <%ExpressionDumpTpl.dumpExp(exp,"\"")%>')

case CALL(path=IDENT(name="der"), expLst={arg as CREF(__)}) then
cref(crefPrefixDer(arg.componentRef))
case CALL(path=IDENT(name="der"), expLst={exp}) then
Expand Down
1 change: 1 addition & 0 deletions SimulationRuntime/c/optimization/optimizer_main.c
View file
Expand Up @@ -44,6 +44,7 @@ int runOptimizer(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo){
OptData *optData, optData_;

solverInfo->solverData = &optData_;
data->simulationInfo->noThrowDivZero = 1;

pickUpModelData(data, threadData, solverInfo);
optData = (OptData*) solverInfo->solverData;
Expand Down
10 changes: 8 additions & 2 deletions SimulationRuntime/c/util/division.h
View file
Expand Up @@ -62,8 +62,14 @@ static inline modelica_real __OMC_DIV_SIM(threadData_t *threadData, const modeli
else
res = a / division_error_equation_time(threadData, a, b, msg, equationIndexes, time_, noThrowDivZero);

if(!valid_number(res))
throwStreamPrintWithEquationIndexes(threadData, equationIndexes, "division leads to inf or nan at time %g, (a=%g) / (b=%g), where divisor b is: %s", time_, a, b, msg);
if(!valid_number(res)){
if(noThrowDivZero) {
warningStreamPrintWithEquationIndexes(LOG_UTIL, 0, equationIndexes, "division leads to inf or nan at time %g, (a=%g) / (b=%g), where divisor b is: %s", time_, a, b, msg);
}
else {
throwStreamPrintWithEquationIndexes(threadData, equationIndexes, "division leads to inf or nan at time %g, (a=%g) / (b=%g), where divisor b is: %s", time_, a, b, msg);
}
}
return res;
}

Expand Down

0 comments on commit feb746e

Please sign in to comment.