[SC] remove the replacement of the DIV operator by a call

 - added thow execption to the optimization, since ipopt
   seems to be abel to handle it

Belonging to [master]:
  - #2379
  - OpenModelica/OpenModelica-testsuite#927

Author: Willi Braun

Compiler/SimCode/SimCodeUtil.mo

@@ -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
@@ -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);
 
@@ -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."

Compiler/Template/CodegenCFunctions.tpl

@@ -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)
@@ -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"
@@ -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

SimulationRuntime/c/optimization/optimizer_main.c

@@ -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;

SimulationRuntime/c/util/division.h

@@ -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;
 }