From ad95e1a338b6da41309d67738e5ffbfc46c17865 Mon Sep 17 00:00:00 2001 From: kabdelhak Date: Fri, 6 Sep 2019 10:41:23 +0200 Subject: [PATCH] [BE] implement sanity check for artificial states - ticket 5459 - revert some stateSelect.prefer based on partial derivatives - revert stateSelect.prefer if contained in smooth(0, cr) - force in stateSelect.never more rigorously - make function tree more accessible when stuff gets differentiated - add some convenience functions --- .../Compiler/BackEnd/AdjacencyMatrix.mo | 12 + .../Compiler/BackEnd/BackendDAECreate.mo | 10 +- .../Compiler/BackEnd/BackendDAEOptimize.mo | 6 +- OMCompiler/Compiler/BackEnd/BackendDAEUtil.mo | 117 +++---- OMCompiler/Compiler/BackEnd/BackendDump.mo | 34 ++ .../Compiler/BackEnd/BackendVariable.mo | 11 + .../Compiler/BackEnd/ExpressionSolve.mo | 3 +- OMCompiler/Compiler/BackEnd/HpcOmEqSystems.mo | 10 +- OMCompiler/Compiler/BackEnd/IndexReduction.mo | 303 +++++++++++++++--- OMCompiler/Compiler/BackEnd/Matching.mo | 204 +++++++++--- .../Compiler/BackEnd/SymbolicJacobian.mo | 2 + OMCompiler/Compiler/FrontEnd/Expression.mo | 33 ++ OMCompiler/Compiler/Util/Array.mo | 35 ++ OMCompiler/Compiler/Util/Error.mo | 32 +- OMCompiler/Compiler/Util/List.mo | 25 +- 15 files changed, 642 insertions(+), 195 deletions(-) diff --git a/OMCompiler/Compiler/BackEnd/AdjacencyMatrix.mo b/OMCompiler/Compiler/BackEnd/AdjacencyMatrix.mo index b078873433d..9ec27d6aeba 100644 --- a/OMCompiler/Compiler/BackEnd/AdjacencyMatrix.mo +++ b/OMCompiler/Compiler/BackEnd/AdjacencyMatrix.mo @@ -284,5 +284,17 @@ algorithm end for; end absAdjacencyMatrix; +public function isEmpty + input BackendDAE.AdjacencyMatrix m; + output Boolean b = true; +algorithm + for element in m loop + if not listEmpty(element) then + b := false; + return; + end if; + end for; +end isEmpty; + annotation(__OpenModelica_Interface="backend"); end AdjacencyMatrix; diff --git a/OMCompiler/Compiler/BackEnd/BackendDAECreate.mo b/OMCompiler/Compiler/BackEnd/BackendDAECreate.mo index 978d593ec29..6983946d257 100644 --- a/OMCompiler/Compiler/BackEnd/BackendDAECreate.mo +++ b/OMCompiler/Compiler/BackEnd/BackendDAECreate.mo @@ -1096,12 +1096,14 @@ protected function lowerVarkind output BackendDAE.VarKind outVarKind; algorithm outVarKind := match(inVarKind, daeAttr) - // variable -> state if have stateSelect = StateSelect.always + // variable -> artificial state if have stateSelect = StateSelect.always case (DAE.VARIABLE(), SOME(DAE.VAR_ATTR_REAL(stateSelectOption = SOME(DAE.ALWAYS())))) + guard(not Types.isDiscreteType(inType)) then BackendDAE.STATE(1, NONE(), false); - // variable -> state if have stateSelect = StateSelect.prefer + // variable -> artificial state if have stateSelect = StateSelect.prefer case (DAE.VARIABLE(), SOME(DAE.VAR_ATTR_REAL(stateSelectOption = SOME(DAE.PREFER())))) + guard(not Types.isDiscreteType(inType)) then BackendDAE.STATE(1, NONE(), false); else @@ -2022,7 +2024,7 @@ algorithm case DAE.DEFINE(componentRef = cr, exp = e, source = source)::xs equation - (e, _) = ExpressionSolve.solve(Expression.crefExp(cr), e, Expression.crefExp(cr)); + (e, _) = ExpressionSolve.solve(Expression.crefExp(cr), e, Expression.crefExp(cr), NONE()); (DAE.PARTIAL_EQUATION(e), source) = ExpressionSimplify.simplifyAddSymbolicOperation(DAE.PARTIAL_EQUATION(e),source); whenOp = BackendDAE.ASSIGN(Expression.crefExp(cr), e, source); whenEq = BackendDAE.WHEN_STMTS(inCond, {whenOp}, NONE()); @@ -2052,7 +2054,7 @@ algorithm case (el as DAE.EQUATION(exp = (cre as DAE.CREF()), scalar = e, source = source))::xs algorithm try - e := ExpressionSolve.solve(cre, e, cre); + e := ExpressionSolve.solve(cre, e, cre, NONE()); else Error.addCompilerError("Failed to solve " + DAEDump.dumpElementsStr({el})); fail(); diff --git a/OMCompiler/Compiler/BackEnd/BackendDAEOptimize.mo b/OMCompiler/Compiler/BackEnd/BackendDAEOptimize.mo index 0a1db50c390..bc29f94e2a8 100644 --- a/OMCompiler/Compiler/BackEnd/BackendDAEOptimize.mo +++ b/OMCompiler/Compiler/BackEnd/BackendDAEOptimize.mo @@ -600,7 +600,7 @@ algorithm (_,(false,_,_,_,_)) = Expression.traverseExpTopDown(e2, traversingTimeEqnsFinder, (false,vars,globalKnownVars,true,false)); cr = BackendVariable.varCref(var); cre = Expression.crefExp(cr); - (_,{}) = ExpressionSolve.solve(e1,e2,cre); + (_,{}) = ExpressionSolve.solve(e1,e2,cre, NONE()); then (); // a = const case ({i},_,_,_,_) @@ -619,7 +619,7 @@ algorithm (_,(false,_,_,_,_)) = Expression.traverseExpTopDown(e2, traversingTimeEqnsFinder, (false,vars,globalKnownVars,false,false)); cr = BackendVariable.varCref(var); cre = Expression.crefExp(cr); - (_,{}) = ExpressionSolve.solve(e1,e2,cre); + (_,{}) = ExpressionSolve.solve(e1,e2,cre, NONE()); then (); // a = der(b) case ({_,_},_,_,_,_) @@ -3654,7 +3654,7 @@ algorithm try BackendDAE.EQUATION(exp = lhs, scalar = rhs) := potentialGlobalKnownEquation; crefExp := BackendVariable.varExp(potentialLocalKnownVar); - (binding,_) := ExpressionSolve.solve(lhs,rhs,crefExp); + (binding,_) := ExpressionSolve.solve(lhs,rhs,crefExp, NONE()); potentialLocalKnownVar := BackendVariable.setBindExp(potentialLocalKnownVar, SOME(binding)); localKnownVars := vindex::localKnownVars; localKnownEqns := eindex::localKnownEqns; diff --git a/OMCompiler/Compiler/BackEnd/BackendDAEUtil.mo b/OMCompiler/Compiler/BackEnd/BackendDAEUtil.mo index d8e0d604f12..4d92e36884f 100644 --- a/OMCompiler/Compiler/BackEnd/BackendDAEUtil.mo +++ b/OMCompiler/Compiler/BackEnd/BackendDAEUtil.mo @@ -4007,7 +4007,6 @@ public function getAdjacencyMatrixEnhancedScalar input BackendDAE.EqSystem syst; input BackendDAE.Shared shared; input Boolean trytosolve "determine the solvability by solving for the variable instead of deriving, needed for 'Casual Tearing Set'"; - input Boolean staticStateSelection = false "Consider states with stateSelect.never as always solvable. Only for index reduction."; output BackendDAE.AdjacencyMatrixEnhanced outIncidenceMatrix; output BackendDAE.AdjacencyMatrixTEnhanced outIncidenceMatrixT; output array> outMapEqnIncRow; @@ -4036,7 +4035,7 @@ algorithm arrT = arrayCreate(numberofVars, {}); // create the array to mark if a variable is allready found in the equation rowmark = arrayCreate(numberofVars, 0); - (arr,arrT,mapEqnIncRow,mapIncRowEqn,varsSolvedInWhenEqnsTupleList) = adjacencyMatrixDispatchEnhancedScalar(vars, eqns,arrT, numberOfEqs, rowmark,globalKnownVars ,trytosolve, staticStateSelection); + (arr,arrT,mapEqnIncRow,mapIncRowEqn,varsSolvedInWhenEqnsTupleList) = adjacencyMatrixDispatchEnhancedScalar(vars, eqns,arrT, numberOfEqs, rowmark,globalKnownVars ,trytosolve, shared); then (arr,arrT,mapEqnIncRow,mapIncRowEqn); @@ -4117,7 +4116,7 @@ protected function adjacencyMatrixDispatchEnhancedScalar input array rowmark; input BackendDAE.Variables globalKnownVars; input Boolean trytosolve; - input Boolean staticStateSelection; + input BackendDAE.Shared shared; output array> omapEqnIncRow; output array omapIncRowEqn; output list>> varsSolvedInWhenEqnsTupleListOut = {}; @@ -4135,7 +4134,7 @@ algorithm // get the equation e := BackendEquation.get(eqArr, i1); // compute the row - (row,size,varsSolvedInWhenEqnsTuple) := adjacencyRowEnhanced(vars, e, i1, rowmark, globalKnownVars, trytosolve, staticStateSelection); + (row,size,varsSolvedInWhenEqnsTuple) := adjacencyRowEnhanced(vars, e, i1, rowmark, globalKnownVars, trytosolve, shared); rowindxs := List.intRange2(rowSize+1, rowSize+size); rowSize := rowSize + size; imapIncRowEqn := List.consN(size,i1,imapIncRowEqn); @@ -4180,7 +4179,7 @@ algorithm arrT = arrayCreate(numberofVars, {}); // create the array to mark if a variable is allready found in the equation rowmark = arrayCreate(numberofVars, 0); - (arr,arrT) = adjacencyMatrixDispatchEnhanced(vars, eqns, arr,arrT, 0, numberOfEqs, intLt(0, numberOfEqs),rowmark,globalKnownVars,trytosolve); + (arr,arrT) = adjacencyMatrixDispatchEnhanced(vars, eqns, arr,arrT, 0, numberOfEqs, intLt(0, numberOfEqs),rowmark,globalKnownVars,trytosolve,shared); then (arr,arrT); @@ -4206,7 +4205,7 @@ protected function adjacencyMatrixDispatchEnhanced input array rowmark; input BackendDAE.Variables globalKnownVars; input Boolean trytosolve; - input Boolean staticStateSelection = false; + input BackendDAE.Shared shared; output BackendDAE.AdjacencyMatrixEnhanced outIncidenceArray; output BackendDAE.AdjacencyMatrixTEnhanced outIncidenceArrayT; algorithm @@ -4231,11 +4230,11 @@ algorithm // get the equation e = BackendEquation.get(eqArr, i1); // compute the row - (row,_,_) = adjacencyRowEnhanced(vars, e, i1, rowmark, globalKnownVars, trytosolve,staticStateSelection); + (row,_,_) = adjacencyRowEnhanced(vars, e, i1, rowmark, globalKnownVars, trytosolve,shared); // put it in the arrays iArr = arrayUpdate(inIncidenceArray, i1, row); iArrT = fillincAdjacencyMatrixTEnhanced(row,{i1},inIncidenceArrayT); - (iArr,iArrT) = adjacencyMatrixDispatchEnhanced(vars, eqArr, iArr, iArrT, i1, numberOfEqs, intLt(i1, numberOfEqs), rowmark, globalKnownVars, trytosolve); + (iArr,iArrT) = adjacencyMatrixDispatchEnhanced(vars, eqArr, iArr, iArrT, i1, numberOfEqs, intLt(i1, numberOfEqs), rowmark, globalKnownVars, trytosolve, shared); then (iArr,iArrT); end match; @@ -4303,7 +4302,7 @@ protected function adjacencyRowEnhanced input array rowmark; input BackendDAE.Variables globalKnownVars; input Boolean trytosolve; - input Boolean staticStateSelection; + input BackendDAE.Shared shared; output BackendDAE.AdjacencyMatrixElementEnhanced outRow; output Integer size; output list>> varsSolvedInWhenEqnsTuple={}; @@ -4330,7 +4329,7 @@ algorithm equation lst = adjacencyRowExpEnhanced(e1, vars, mark, rowmark, {}); lst = adjacencyRowExpEnhanced(e2, vars, mark, rowmark, lst); - row = adjacencyRowEnhanced1(lst,e1,e2,vars,globalKnownVars,mark,rowmark,{},trytosolve,staticStateSelection); + row = adjacencyRowEnhanced1(lst,e1,e2,vars,globalKnownVars,mark,rowmark,{},trytosolve,shared); then (row,1); // COMPLEX_EQUATION @@ -4338,7 +4337,7 @@ algorithm equation lst = adjacencyRowExpEnhanced(e1, vars, mark, rowmark, {}); lst = adjacencyRowExpEnhanced(e2, vars, mark, rowmark, lst); - row = adjacencyRowEnhanced1(lst,e1,e2,vars,globalKnownVars,mark,rowmark,{},trytosolve,staticStateSelection); + row = adjacencyRowEnhanced1(lst,e1,e2,vars,globalKnownVars,mark,rowmark,{},trytosolve,shared); then (row,size); // ARRAY_EQUATION @@ -4346,7 +4345,7 @@ algorithm equation lst = adjacencyRowExpEnhanced(e1, vars, mark, rowmark, {}); lst = adjacencyRowExpEnhanced(e2, vars, mark, rowmark, lst); - row = adjacencyRowEnhanced1(lst,e1,e2,vars,globalKnownVars,mark,rowmark,{},trytosolve,staticStateSelection); + row = adjacencyRowEnhanced1(lst,e1,e2,vars,globalKnownVars,mark,rowmark,{},trytosolve,shared); size = List.fold(ds,intMul,1); then (row,size); @@ -4357,21 +4356,21 @@ algorithm expCref = Expression.crefExp(cr); lst = adjacencyRowExpEnhanced(expCref, vars, mark, rowmark, {}); lst = adjacencyRowExpEnhanced(e, vars, mark, rowmark, lst); - row = adjacencyRowEnhanced1(lst,expCref,e,vars,globalKnownVars,mark,rowmark,{},trytosolve,staticStateSelection); + row = adjacencyRowEnhanced1(lst,expCref,e,vars,globalKnownVars,mark,rowmark,{},trytosolve,shared); then (row,1); // RESIDUAL_EQUATION case (vars,BackendDAE.RESIDUAL_EQUATION(exp = e),_,_,_) equation lst = adjacencyRowExpEnhanced(e, vars, mark, rowmark, {}); - row = adjacencyRowEnhanced1(lst,e,DAE.RCONST(0.0),vars,globalKnownVars,mark,rowmark,{},trytosolve,staticStateSelection); + row = adjacencyRowEnhanced1(lst,e,DAE.RCONST(0.0),vars,globalKnownVars,mark,rowmark,{},trytosolve,shared); then (row,1); // WHEN_EQUATION case (vars,BackendDAE.WHEN_EQUATION(size=size,whenEquation = elsewe),_,_,_) equation - (row,varsSolvedInWhenEqns) = adjacencyRowWhenEnhanced(elsewe, mark, rowmark, vars, globalKnownVars, {}, {},staticStateSelection); + (row,varsSolvedInWhenEqns) = adjacencyRowWhenEnhanced(elsewe, mark, rowmark, vars, globalKnownVars, {}, {},shared); varsSolvedInWhenEqnsTuple = {(mark,varsSolvedInWhenEqns)}; then (row,size); @@ -4394,7 +4393,7 @@ algorithm // special case for it initial() then ... else ... end if; only else branch needs to be checked case(_,BackendDAE.IF_EQUATION(conditions={DAE.CALL(path=Absyn.IDENT("initial"))},eqnstrue={_},eqnsfalse=eqnselse),_,_,_) equation - (row,size) = adjacencyRowEnhancedEqnLst(eqnselse,inVariables,mark,rowmark,globalKnownVars,trytosolve,staticStateSelection); + (row,size) = adjacencyRowEnhancedEqnLst(eqnselse,inVariables,mark,rowmark,globalKnownVars,trytosolve,shared); then (row,size); @@ -4410,12 +4409,12 @@ algorithm // mark all negative because the when condition cannot used to solve a variable lst = List.fold3(expl, adjacencyRowExpEnhanced, vars, mark, rowmark, {}); _ = List.fold1(lst,markNegativ,rowmark,mark); - row1 = adjacencyRowEnhanced1(lst,DAE.RCONST(0.0),DAE.RCONST(0.0),vars,globalKnownVars,mark,rowmark,{},trytosolve,staticStateSelection); + row1 = adjacencyRowEnhanced1(lst,DAE.RCONST(0.0),DAE.RCONST(0.0),vars,globalKnownVars,mark,rowmark,{},trytosolve,shared); - (row, size) = adjacencyRowEnhancedEqnLst(eqnselse, vars, mark, rowmark, globalKnownVars, trytosolve,staticStateSelection); + (row, size) = adjacencyRowEnhancedEqnLst(eqnselse, vars, mark, rowmark, globalKnownVars, trytosolve,shared); lst = List.map(row,Util.tuple31); - (lst, row, size) = List.fold6(eqnslst, adjacencyRowEnhancedEqnLstIfBranches, vars, mark, rowmark, globalKnownVars, trytosolve, staticStateSelection, (lst, row, size)); + (lst, row, size) = List.fold6(eqnslst, adjacencyRowEnhancedEqnLstIfBranches, vars, mark, rowmark, globalKnownVars, trytosolve, shared, (lst, row, size)); lstall = List.map(row, Util.tuple31); (_, lst, _) = List.intersection1OnTrue(lstall, lst, intEq); @@ -4441,7 +4440,7 @@ protected function adjacencyRowEnhancedEqnLstIfBranches input array rowmark; input BackendDAE.Variables globalKnownVars; input Boolean trytosolve; - input Boolean staticStateSelection; + input BackendDAE.Shared shared; input tuple,BackendDAE.AdjacencyMatrixElementEnhanced,Integer> intpl; output tuple,BackendDAE.AdjacencyMatrixElementEnhanced,Integer> outtpl; protected @@ -4451,7 +4450,7 @@ protected algorithm (inLstAllBranch,iRow,iSize) := intpl; for eqn in iEqns loop - (row,size,_) := adjacencyRowEnhanced(inVariables, eqn, mark, rowmark, globalKnownVars, trytosolve,staticStateSelection); + (row,size,_) := adjacencyRowEnhanced(inVariables, eqn, mark, rowmark, globalKnownVars, trytosolve,shared); lst := List.map(row,Util.tuple31); inLstAllBranch := List.intersectionOnTrue(lst, inLstAllBranch,intEq); iSize := iSize + size; @@ -4467,7 +4466,7 @@ protected function adjacencyRowEnhancedEqnLst input array rowmark; input BackendDAE.Variables globalKnownVars; input Boolean trytosolve; - input Boolean staticStateSelection; + input BackendDAE.Shared shared; output BackendDAE.AdjacencyMatrixElementEnhanced outRow = {}; output Integer oSize = 0; protected @@ -4475,7 +4474,7 @@ protected Integer size; algorithm for eqn in iEqns loop - (row,size,_) := adjacencyRowEnhanced(inVariables,eqn,mark,rowmark,globalKnownVars,trytosolve,staticStateSelection); + (row,size,_) := adjacencyRowEnhanced(inVariables,eqn,mark,rowmark,globalKnownVars,trytosolve,shared); outRow := listAppend(row,outRow); oSize := oSize + size; end for; @@ -4598,7 +4597,7 @@ protected function adjacencyRowWhenEnhanced input BackendDAE.Variables globalKnownVars; input list iLst; input BackendDAE.AdjacencyMatrixElementEnhanced iRow; - input Boolean staticStateSelection; + input BackendDAE.Shared shared; output BackendDAE.AdjacencyMatrixElementEnhanced outRow = iRow; output list varsSolvedInWhenEqns={}; protected @@ -4628,7 +4627,7 @@ algorithm _ = List.fold1(lst,markNegativ,rowmark,mark); //leftexp = Expression.crefExp(left); lst = adjacencyRowExpEnhanced(leftexp, vars, mark, rowmark, lst); - outRow = adjacencyRowEnhanced1(lst,leftexp,right,vars,globalKnownVars,mark,rowmark,outRow,false,staticStateSelection); + outRow = adjacencyRowEnhanced1(lst,leftexp,right,vars,globalKnownVars,mark,rowmark,outRow,false,shared); then (); case BackendDAE.ASSIGN(leftexp, right) equation @@ -4639,7 +4638,7 @@ algorithm // mark all negative because the when condition cannot used to solve a variable _ = List.fold1(lst,markNegativ,rowmark,mark); lst = adjacencyRowExpEnhanced(leftexp, vars, mark, rowmark, lst); - outRow = adjacencyRowEnhanced1(lst,leftexp,right,vars,globalKnownVars,mark,rowmark,outRow,false,staticStateSelection); + outRow = adjacencyRowEnhanced1(lst,leftexp,right,vars,globalKnownVars,mark,rowmark,outRow,false,shared); then (); else (); @@ -4648,7 +4647,7 @@ algorithm end for; if isSome(oelsepart) then SOME(elsepart) := oelsepart; - outRow := adjacencyRowWhenEnhanced(elsepart, mark, rowmark, vars, globalKnownVars, lst, outRow, staticStateSelection); + outRow := adjacencyRowWhenEnhanced(elsepart, mark, rowmark, vars, globalKnownVars, lst, outRow, shared); end if; end adjacencyRowWhenEnhanced; @@ -4678,7 +4677,7 @@ protected function adjacencyRowEnhanced1 input array rowmark; input BackendDAE.AdjacencyMatrixElementEnhanced inRow; input Boolean trytosolve; - input Boolean staticStateSelection; + input BackendDAE.Shared shared; output BackendDAE.AdjacencyMatrixElementEnhanced outRow; algorithm outRow := matchcontinue(lst,e1,e2,vars,globalKnownVars,mark,rowmark,inRow) @@ -4701,24 +4700,6 @@ algorithm then adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_UNSOLVABLE())::inRow,trytosolve); */ - // Consider vars with stateSelect.never as solvable. Only for staticStateSelection - case(r::rest,_,_,_,_,_,_,_) - guard(staticStateSelection) - equation - BackendDAE.VAR(varName=cr,varKind=BackendDAE.STATE(),values = SOME(DAE.VAR_ATTR_REAL(stateSelectOption = SOME(DAE.NEVER())))) = BackendVariable.getVarAt(vars, intAbs(r)); - de = Differentiate.differentiateExpSolve(e1, cr, NONE()); - (de,_) = ExpressionSimplify.simplify(de); - (_,crlst) = Expression.traverseExpBottomUp(de, Expression.traversingComponentRefFinder, {}); - solvab = adjacencyRowEnhanced2(cr,de,crlst,vars,globalKnownVars); - b = isSolvable(solvab,true); - de = Differentiate.differentiateExpSolve(e2, cr, NONE()); - (de,_) = ExpressionSimplify.simplify(de); - (_,crlst) = Expression.traverseExpBottomUp(de, Expression.traversingComponentRefFinder, {}); - solvab = adjacencyRowEnhanced2(cr,de,crlst,vars,globalKnownVars); - b2 = isSolvable(solvab,true); - true = (b or b2) and not (b and b2); // either lhs or rhs is solvable - then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_CONSTONE(),{})::inRow,trytosolve,staticStateSelection); case(r::rest,DAE.CALL(path= Absyn.IDENT("der"),expLst={DAE.CREF(componentRef = cr)}),_,_,_,_,_,_) guard intGt(r,0) @@ -4730,7 +4711,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasDerCref(e2,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.CALL(path= Absyn.IDENT("der"),expLst={DAE.CREF(componentRef = cr)}),_,_,_,_,_) equation rabs = intAbs(r); @@ -4741,7 +4722,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasDerCref(e1,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.CREF(componentRef=cr),_,_,_,_,_,_) equation rabs = intAbs(r); @@ -4752,7 +4733,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasCrefNoPreorDer(e2,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.CREF(componentRef=cr),_,_,_,_,_,_) equation rabs = intAbs(r); @@ -4764,7 +4745,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, crarr); false = Expression.expHasCrefNoPreorDer(e2,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.LUNARY(operator=DAE.NOT(_),exp=DAE.CREF(componentRef=cr)),_,_,_,_,_,_) equation rabs = intAbs(r); @@ -4775,7 +4756,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasCrefNoPreorDer(e2,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.UNARY(operator=DAE.UMINUS(_),exp=DAE.CREF(componentRef=cr)),_,_,_,_,_,_) equation rabs = intAbs(r); @@ -4786,7 +4767,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasCrefNoPreorDer(e2,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.UNARY(operator=DAE.UMINUS_ARR(_),exp=DAE.CREF(componentRef=cr)),_,_,_,_,_,_) equation rabs = intAbs(r); @@ -4798,7 +4779,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, crarr); false = Expression.expHasCrefNoPreorDer(e2,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.CREF(componentRef=cr),_,_,_,_,_) equation rabs = intAbs(r); @@ -4809,7 +4790,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasCrefNoPreorDer(e1,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.CREF(componentRef=cr),_,_,_,_,_) equation rabs = intAbs(r); @@ -4821,7 +4802,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, crarr); false = Expression.expHasCrefNoPreorDer(e1,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.LUNARY(operator=DAE.NOT(_),exp=DAE.CREF(componentRef=cr)),_,_,_,_,_) equation rabs = intAbs(r); @@ -4832,7 +4813,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasCrefNoPreorDer(e1,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.UNARY(operator=DAE.UMINUS(_),exp=DAE.CREF(componentRef=cr)),_,_,_,_,_) equation rabs = intAbs(r); @@ -4843,7 +4824,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, cr1); false = Expression.expHasCrefNoPreorDer(e1,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.UNARY(operator=DAE.UMINUS_ARR(_),exp=DAE.CREF(componentRef=cr)),_,_,_,_,_) equation rabs = intAbs(r); @@ -4855,7 +4836,7 @@ algorithm true = ComponentReference.crefEqualNoStringCompare(cr, crarr); false = Expression.expHasCrefNoPreorDer(e1,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.CREF(componentRef=cr),_,_,_,_,_,_) equation rabs = intAbs(r); @@ -4866,7 +4847,7 @@ algorithm true = ComponentReference.crefPrefixOf(cr, cr1); false = Expression.expHasCrefNoPreorDer(e2,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.CREF(componentRef=cr),_,_,_,_,_) equation rabs = intAbs(r); @@ -4877,7 +4858,7 @@ algorithm true = ComponentReference.crefPrefixOf(cr, cr1); false = Expression.expHasCrefNoPreorDer(e1,cr); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.CALL(path=path,expLst=explst,attr=DAE.CALL_ATTR(ty= DAE.T_COMPLEX(complexClassType=ClassInf.RECORD(path1)))),_,_,_,_,_,_) equation true = AbsynUtil.pathEqual(path,path1); @@ -4889,7 +4870,7 @@ algorithm true = expCrefLstHasCref(explst,cr1); false = Expression.expHasCrefNoPreorDer(e2,cr1); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,DAE.CALL(path=path,expLst=explst,attr=DAE.CALL_ATTR(ty= DAE.T_COMPLEX(complexClassType=ClassInf.RECORD(path1)))),_,_,_,_,_) equation true = AbsynUtil.pathEqual(path,path1); @@ -4901,7 +4882,7 @@ algorithm true = expCrefLstHasCref(explst,cr1); false = Expression.expHasCrefNoPreorDer(e1,cr1); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,DAE.TUPLE(PR=explst),DAE.CALL(),_,_,_,_,_) equation rabs = intAbs(r); @@ -4916,7 +4897,7 @@ algorithm true = expCrefLstHasCref(crexplst,cr1); false = Expression.expHasCrefNoPreorDer(e2,cr1); then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_SOLVED(),{})::inRow,trytosolve,shared); case(r::rest,_,_,_,_,_,_,_) // case: state derivative equation @@ -4929,7 +4910,7 @@ algorithm e = Expression.crefExp(cr); ((e,_)) = Expression.replaceExp(Expression.expSub(e1,e2), DAE.CALL(Absyn.IDENT("der"),{e},DAE.callAttrBuiltinReal), Expression.crefExp(cr1)); e_derAlias = Expression.traverseExpDummy(e, replaceDerCall); - (de,solved,derived,cons) = tryToSolveOrDerive(e_derAlias, cr1, vars, NONE(),trytosolve); + (de,solved,derived,cons) = tryToSolveOrDerive(e_derAlias, cr1, vars, SOME(shared.functionTree),trytosolve); if not solved then (de,_) = ExpressionSimplify.simplify(de); (_,crlst) = Expression.traverseExpBottomUp(de, Expression.traversingComponentRefFinder, {}); @@ -4945,7 +4926,7 @@ algorithm end if; end if; then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,solvab,cons)::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,solvab,cons)::inRow,trytosolve,shared); case(r::rest,_,_,_,_,_,_,_) equation rabs = intAbs(r); @@ -4955,7 +4936,7 @@ algorithm BackendDAE.VAR(varName=cr) = BackendVariable.getVarAt(vars, rabs); e = Expression.expSub(e1,e2); e_derAlias = Expression.traverseExpDummy(e, replaceDerCall); - (de,solved,derived,cons) = tryToSolveOrDerive(e_derAlias, cr, vars, NONE(),trytosolve); + (de,solved,derived,cons) = tryToSolveOrDerive(e_derAlias, cr, vars, SOME(shared.functionTree),trytosolve); if not solved then (de,_) = ExpressionSimplify.simplify(de); (_,crlst) = Expression.traverseExpTopDown(de, Expression.traversingComponentRefFinderNoPreDer, {}); @@ -4971,10 +4952,10 @@ algorithm end if; end if; then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,solvab,cons)::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,solvab,cons)::inRow,trytosolve,shared); case(r::rest,_,_,_,_,_,_,_) then - adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_UNSOLVABLE(),{})::inRow,trytosolve,staticStateSelection); + adjacencyRowEnhanced1(rest,e1,e2,vars,globalKnownVars,mark,rowmark,(r,BackendDAE.SOLVABILITY_UNSOLVABLE(),{})::inRow,trytosolve,shared); end matchcontinue; end adjacencyRowEnhanced1; diff --git a/OMCompiler/Compiler/BackEnd/BackendDump.mo b/OMCompiler/Compiler/BackEnd/BackendDump.mo index 07f1160944d..59556715bc6 100644 --- a/OMCompiler/Compiler/BackEnd/BackendDump.mo +++ b/OMCompiler/Compiler/BackEnd/BackendDump.mo @@ -454,6 +454,40 @@ algorithm outTpl := (varNo + 1, buffer); end var1String; +public function varListStringShort + input list inVars; + input String heading; + output String outString; +algorithm + outString := match(inVars, heading) + local + String buffer; + + case (_, "") equation + ((_, buffer)) = List.fold(inVars, varNameString, (1, "")); + then buffer; + + else equation + ((_, buffer)) = List.fold(inVars, varNameString, (1, "")); + buffer = heading + "\n" + UNDERLINE + "\n" + buffer; + then buffer; + end match; +end varListStringShort; + +protected function varNameString + input BackendDAE.Var inVar; + input tuple inTpl; + output tuple outTpl; +protected + Integer varNo; + String buffer; +algorithm + (varNo, buffer) := inTpl; + buffer := buffer + intString(varNo) + ": "; + buffer := buffer + ComponentReference.printComponentRefStr(inVar.varName) + "\n"; + outTpl := (varNo + 1, buffer); +end varNameString; + public function varListStringIndented input list inVars; input String heading; diff --git a/OMCompiler/Compiler/BackEnd/BackendVariable.mo b/OMCompiler/Compiler/BackEnd/BackendVariable.mo index dc60f69ce04..fb9bd66358e 100644 --- a/OMCompiler/Compiler/BackEnd/BackendVariable.mo +++ b/OMCompiler/Compiler/BackEnd/BackendVariable.mo @@ -460,6 +460,17 @@ algorithm end match; end varStateSelect; +public function varStateSelectNever "author: kabdelhak + Returns true, if the state select attribute is DAE.NEVER()" + input BackendDAE.Var inVar; + output Boolean isNever; +algorithm + isNever := match(varStateSelect(inVar)) + case DAE.NEVER() then true; + else false; + end match; +end varStateSelectNever; + public function setVarStateSelect "Sets the state select attribute of a variable." input BackendDAE.Var inVar; input DAE.StateSelect stateSelect; diff --git a/OMCompiler/Compiler/BackEnd/ExpressionSolve.mo b/OMCompiler/Compiler/BackEnd/ExpressionSolve.mo index 08ce572e099..91a5bbf5bab 100644 --- a/OMCompiler/Compiler/BackEnd/ExpressionSolve.mo +++ b/OMCompiler/Compiler/BackEnd/ExpressionSolve.mo @@ -193,6 +193,7 @@ public function solve input DAE.Exp inExp1 "lhs"; input DAE.Exp inExp2 "rhs"; input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'"; + input Option functions = NONE() "need for solve modelica functions"; output DAE.Exp outExp; output list outAsserts; protected @@ -209,7 +210,7 @@ algorithm (outExp,outAsserts,dummy1, dummy2, dummyI) := matchcontinue inExp1 case _ then solveSimple(inExp1, inExp2, inExp3, 0); case _ then solveSimple(inExp2, inExp1, inExp3, 0); - case _ then solveWork(inExp1, inExp2, inExp3, NONE(), NONE(), 0, false, false); + case _ then solveWork(inExp1, inExp2, inExp3, functions, NONE(), 0, false, false); else equation if Flags.isSet(Flags.FAILTRACE) then Error.addInternalError("Failed to solve \"" + ExpressionDump.printExpStr(inExp1) + " = " + ExpressionDump.printExpStr(inExp2) + "\" w.r.t. \"" + ExpressionDump.printExpStr(inExp3) + "\"", sourceInfo()); diff --git a/OMCompiler/Compiler/BackEnd/HpcOmEqSystems.mo b/OMCompiler/Compiler/BackEnd/HpcOmEqSystems.mo index 9c8582ffc6d..855afaced32 100644 --- a/OMCompiler/Compiler/BackEnd/HpcOmEqSystems.mo +++ b/OMCompiler/Compiler/BackEnd/HpcOmEqSystems.mo @@ -505,7 +505,7 @@ algorithm resEqsOut := hs; // some optimization - (eqsNewOut,varsNewOut,resEqsOut) := simplifyNewEquations(eqsNewOut,varsNewOut,resEqsOut,listLength(List.flatten(arrayList(xa_iArr))),2); + (eqsNewOut,varsNewOut,resEqsOut) := simplifyNewEquations(eqsNewOut,varsNewOut,resEqsOut,listLength(List.flatten(arrayList(xa_iArr))),2,ishared); // handle the strongComponent (system of equations) to solve the tearing vars (compsEqSys,resEqsOut,tVarsOut,addEqLst,addVarLst) := buildEqSystemComponent(resEqIdcs0,tVarIdcs0,resEqsOut,tVarsOut,a_iArr,ishared); @@ -559,6 +559,7 @@ protected function simplifyNewEquations input list resEqsIn; input Integer numAuxiliaryVars; // to prevent replacement of coefficients input Integer numIter; + input BackendDAE.Shared shared; output list eqsOut; output list varsOut; output list resEqsOut; @@ -578,7 +579,7 @@ algorithm eqSys := BackendDAEUtil.createEqSystem(varArr, eqArr); (m,mT) := BackendDAEUtil.incidenceMatrix(eqSys,BackendDAE.ABSOLUTE(),NONE()); size := listLength(eqsIn); - (eqIdcs,varIdcs,resEqsOut) := List.fold(List.intRange(size),function simplifyNewEquations1(eqArr=eqArr,varArr=varArr,m=m,mt=mT,numAuxiliaryVars=numAuxiliaryVars),({},{},resEqsIn)); + (eqIdcs,varIdcs,resEqsOut) := List.fold(List.intRange(size),function simplifyNewEquations1(eqArr=eqArr,varArr=varArr,m=m,mt=mT,numAuxiliaryVars=numAuxiliaryVars,shared=shared),({},{},resEqsIn)); numAux := numAuxiliaryVars-listLength(varIdcs); if listEmpty(varIdcs) then numIterNew:=0; else numIterNew := numIter; @@ -588,7 +589,7 @@ algorithm (_,eqIdcs,_) := List.intersection1OnTrue(List.intRange(size),eqIdcs,intEq); eqsOut := BackendEquation.getList(eqIdcs,eqArr); varsOut := List.map1(varIdcs,BackendVariable.getVarAtIndexFirst,varArr); - if numIterNew<>0 then (eqsOut,varsOut,resEqsOut) := simplifyNewEquations(eqsOut,varsOut,resEqsOut,numAux,numIterNew-1); + if numIterNew<>0 then (eqsOut,varsOut,resEqsOut) := simplifyNewEquations(eqsOut,varsOut,resEqsOut,numAux,numIterNew-1,shared); else (eqsOut,varsOut,resEqsOut) := (eqsOut,varsOut,resEqsOut); end if; end simplifyNewEquations; @@ -600,6 +601,7 @@ protected function simplifyNewEquations1 input BackendDAE.IncidenceMatrix m; input BackendDAE.IncidenceMatrix mt; input Integer numAuxiliaryVars; + input BackendDAE.Shared shared; input tuple,list,list> tplIn; //these can be removed afterwards (eqIdcs,varIdcs,_) output tuple,list,list> tplOut; algorithm @@ -628,7 +630,7 @@ algorithm varExp := Expression.crefExp(varCref); rhs := BackendEquation.getEquationRHS(eq); lhs := BackendEquation.getEquationLHS(eq); - (rhs,_) := ExpressionSolve.solve(lhs,rhs,varExp); + (rhs,_) := ExpressionSolve.solve(lhs,rhs,varExp,NONE()); if Expression.isAsubExp(rhs) then rhs := List.fold1(Expression.allTerms(rhs),Expression.makeBinaryExp,DAE.ADD(Expression.typeof(varExp)),DAE.RCONST(0.0)); //in case ({a,0,b}+funcCall(1,2,3,4,5))[2] I need to get 0+funcCAll(1,2,3,4,5)[2] end if; diff --git a/OMCompiler/Compiler/BackEnd/IndexReduction.mo b/OMCompiler/Compiler/BackEnd/IndexReduction.mo index fddb74d15f3..88c2712188b 100644 --- a/OMCompiler/Compiler/BackEnd/IndexReduction.mo +++ b/OMCompiler/Compiler/BackEnd/IndexReduction.mo @@ -390,7 +390,6 @@ algorithm list ilst, unassignedEqns, eqnsLst, discEqns, stateIndxs; list> rest; Boolean b; - list foundStates, artificial, natural; case {} then (inEqnsLstAcc, inStateIndxsAcc, inUnassEqnsAcc, inDiscEqnsAcc); @@ -402,14 +401,6 @@ algorithm stateIndxs := List.fold2(ilst, statesInEquations, (m, statemark, mark), inAssignments1, {}); // print("stateIndxs " + stringDelimitList(List.map(stateIndxs, intString), ", ") + "\n"); b := intGe(listLength(stateIndxs), listLength(unassignedEqns)); - foundStates := List.map1r(stateIndxs,BackendVariable.getVarAt,vars); - - artificial := list(var for var guard BackendVariable.isArtificialState(var) in foundStates); - natural := list(var for var guard BackendVariable.isNaturalState(var) in foundStates); - - BackendDump.dumpVarList(artificial, "artificial states"); - BackendDump.dumpVarList(natural, "natural states"); - singulareSystemError(b, stateIndxs, unassignedEqns, eqnsLst, inSystem, inShared, inAssignments1, inAssignments2, inArg); (outEqnsLst, outStateIndxs, outUnassEqnsAcc, outDiscEqns) := minimalStructurallySingularSystemMSS(rest, inSystem, inShared, inAssignments1, inAssignments2, inArg, statemark, mark+1, m, vars, eqnsLst::inEqnsLstAcc, stateIndxs::inStateIndxsAcc, unassignedEqns::inUnassEqnsAcc, discEqns); then (outEqnsLst, outStateIndxs, outUnassEqnsAcc, outDiscEqns); @@ -1340,7 +1331,7 @@ algorithm BackendDAE.EQSYSTEM(orderedVars=vars) := inSystem; BackendDAE.SHARED(functionTree=funcs) := inShared; // do late Inline also in orgeqnslst - // orgEqnsLst := inlineOrgEqns(orgEqnsLst,(SOME(funcs),{DAE.NORM_INLINE(),DAE.AFTER_INDEX_RED_INLINE()})); + orgEqnsLst := inlineOrgEqns(orgEqnsLst,(SOME(funcs),{DAE.NORM_INLINE(), DAE.AFTER_INDEX_RED_INLINE()})); if Flags.isSet(Flags.BLT_DUMP) then print("########################### STATE SELECTION ###########################\n"); end if; @@ -2038,7 +2029,7 @@ algorithm (outDummyVars,oStateSets) := matchcontinue inSystem local - list dummyVars,vlst; + list dummyVars,stateVars,vlst; array> mapEqnIncRow; array mapIncRowEqn; Integer nv,nv1,ne,ne1,neqnarr; @@ -2073,8 +2064,8 @@ algorithm hovvars = BackendVariable.listVar1(statecandidates); eqns1 = BackendEquation.listEquation(eqnslst); syst = BackendDAEUtil.createEqSystem(hovvars, eqns1); - (me,meT,_,_) = BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(syst,inShared,false,true); - m1 = incidenceMatrixfromEnhanced2(me,hovvars); + (me,meT,_,_) = BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(syst,inShared,false); + m1 = incidenceMatrixfromEnhancedStrict(me,hovvars); mT1 = AdjacencyMatrix.transposeAdjacencyMatrix(m1,nfreeStates); // BackendDump.printEqSystem(syst); hovvars = sortStateCandidatesVars(hovvars,BackendVariable.daeVars(inSystem),SOME(mT1)); @@ -2134,14 +2125,14 @@ algorithm vars = BackendVariable.addVars(BackendVariable.varList(hovvars), vars); syst = BackendDAEUtil.createEqSystem(vars, eqns); // get advanced incidence Matrix - (me,meT,mapEqnIncRow,mapIncRowEqn) = BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(syst,inShared,false,true); + (me,meT,mapEqnIncRow,mapIncRowEqn) = BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(syst,inShared,false); if Flags.isSet(Flags.BLT_DUMP) then BackendDump.dumpAdjacencyMatrixEnhanced(me); print("\n"); BackendDump.dumpAdjacencyMatrixTEnhanced(meT); end if; // get indicenceMatrix from Enhanced - m = incidenceMatrixfromEnhanced2(me,vars); + m = incidenceMatrixfromEnhancedStrict(me,vars); nv = BackendVariable.varsSize(vars); ne = BackendEquation.equationArraySize(eqns); mT = AdjacencyMatrix.transposeAdjacencyMatrix(m,nv); @@ -2151,28 +2142,36 @@ algorithm vec1 = arrayCreate(nv,-1); vec2 = arrayCreate(ne,-1); BackendDAEEXT.getAssignment(vec1,vec2); + // get the matched state candidates -> dummys + // get the unmatched state candidates -> actual states + // force StateSelect.never vars + (dstates, states, vec1, vec2) = forceStateSelectNever(vec1, vec2, vars, eqns, me, inShared, so); if Flags.isSet(Flags.BLT_DUMP) then print("\n"); BackendDump.dumpMatching(vec1); print("\n"); BackendDump.dumpMatching(vec2); end if; - // get the matched state candidates -> dummyVars (dstates,_) = checkAssignment(1,nv,vec1,vars); dummyVars = List.map1r(List.map(dstates,Util.tuple22),BackendVariable.getVarAt,vars); + stateVars = List.map1r(List.map(states,Util.tuple22),BackendVariable.getVarAt,vars); dummyVars = List.select(dummyVars, BackendVariable.isStateVar); - if Flags.isSet(Flags.BLT_DUMP) then - BackendDump.dumpVarList(dummyVars, "Selected dummy states:"); - end if; // get assigned and unassigned equations unassigned = Matching.getUnassigned(ne, vec2, {}); _ = Matching.getAssigned(ne, vec2, {}); if Flags.isSet(Flags.BLT_DUMP) then if listEmpty(unassigned) then print("Perfect Matching, no dynamic index reduction needed! There are no unassigned equations.\n\n"); + if Flags.isSet(Flags.BLT_DUMP) then + BackendDump.dumpVarList(dummyVars, "Selected dummy states:"); + BackendDump.dumpVarList(stateVars, "Selected continuous states:"); + end if; else print("No perfect matching possible, dynamic index reduction needed.\n"); BackendDump.dumpEquationList(BackendEquation.getEquationArraySubsetLst(eqns,unassigned),"Unassigned equations:"); + if Flags.isSet(Flags.BLT_DUMP) then + BackendDump.dumpVarList(dummyVars, "Statically selected dummy states:"); + end if; print("\n"); end if; end if; @@ -2221,6 +2220,121 @@ algorithm end matchcontinue; end selectStatesWork1; +protected function forceStateSelectNever + input array vec_old1; + input array vec_old2; + input BackendDAE.Variables vars; + input BackendDAE.EquationArray eqns; + input BackendDAE.AdjacencyMatrixEnhanced me; + input BackendDAE.Shared inShared; + input BackendDAE.StateOrder so; + output list> dummyStates; + output list> states; + output array vec1 = vec_old1; + output array vec2 = vec_old2; +protected + Integer nv, nv2, ne, never_i, eq_i, old_i; + BackendDAE.Var var; + HashTable3.HashTable ht; + array vec_res1, vec_res2; + list neverVars = {}; + BackendDAE.Variables neverVarsArray; + list neverIdx = {}; + BackendDAE.EqSystem syst2; + BackendDAE.AdjacencyMatrixEnhanced me2; + BackendDAE.IncidenceMatrix m, mT; + list> tplLst; + String msg; +algorithm + nv := BackendVariable.varsSize(vars); + ne := BackendEquation.equationArraySize(eqns); + BackendDAE.STATEORDER(invHashTable = ht) := so; + (dummyStates,states) := checkAssignment(1,nv,vec_old1,vars); + + for state in states loop + var := BackendVariable.getVarAt(vars,Util.tuple22(state)); + if varStateSelectNever(var) and not BaseHashTable.hasKey(BackendVariable.varCref(var), ht) then + neverVars := var::neverVars; + neverIdx := Util.tuple22(state)::neverIdx; + end if; + end for; + + // If there are any unmatched stateSelect.never vars dump them and start algorithm + if not listEmpty(neverVars) then + if Flags.isSet(Flags.BLT_DUMP) then + BackendDump.dumpVarList(neverVars, "StateSelect.never variables that will tried to be forced as dummys"); + else + msg := System.gettext(BackendDump.varListStringShort(neverVars,"") + + "They were forced to be statically selected as dummys, this could lead to errors during simulation, please use -d=bltdump for more information.\n"); + Error.addMessage(Error.STATE_STATESELECT_NEVER_FORCED, {msg}); + end if; + + // If there are unmatched equations try to match full system with casual rules for stateSelect.never vars + if Matching.anyUnassigned(ne, vec2) then + // non strict rules to be able to solve for nonlinear stateSelect.never states + m := incidenceMatrixfromEnhanced(me,vars,so); + mT := AdjacencyMatrix.transposeAdjacencyMatrix(m,nv); + + (vec2,vec1,_) := Matching.ContinueMatching(mT,ne,nv,vec2,vec1); + (dummyStates,states) := checkAssignment(1,nv,vec1,vars); + + // look for leftover unmatched stateSelect.never vars + neverVars := {}; + neverIdx := {}; + for state in states loop + var := BackendVariable.getVarAt(vars,Util.tuple22(state)); + if varStateSelectNever(var) and not BaseHashTable.hasKey(BackendVariable.varCref(var), ht) then + neverVars := var::neverVars; + neverIdx := Util.tuple22(state)::neverIdx; + end if; + end for; + end if; + + // if there are (still) unmatched never vars try to solve them individually and swap out in original matching + if not listEmpty(neverVars) then + neverVarsArray := BackendVariable.listVar1(neverVars); + + nv2 := BackendVariable.varsSize(neverVarsArray); + syst2 := BackendDAEUtil.createEqSystem(neverVarsArray, eqns); + (me2,_,_,_) := BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(syst2,inShared,false); + + // non strict rules to be able to solve for nonlinear stateSelect.never states + m := incidenceMatrixfromEnhancedPartial(me2,vars,neverVarsArray,vec2,so); + + if not AdjacencyMatrix.isEmpty(m) then + mT := AdjacencyMatrix.transposeAdjacencyMatrix(m,nv2); + Matching.matchingExternalsetIncidenceMatrix(ne,nv2,mT); + BackendDAEEXT.matching(ne,nv2,3,-1,1.0,1); + vec_res1 := arrayCreate(nv2,-1); + vec_res2 := arrayCreate(ne,-1); + BackendDAEEXT.getAssignment(vec_res1,vec_res2); + + // swap out new matchings + tplLst := List.zip(neverIdx,List.intRange(listLength(neverIdx))); + for tpl in tplLst loop + (never_i, eq_i) := tpl; + // assigning matched never var to matched eq + vec1[never_i] := vec_res1[eq_i]; + old_i := vec2[vec_res1[eq_i]]; + vec2[vec_res1[eq_i]] := never_i; + // unassigning old matched var if the equation was matched + if not intEq(old_i,-1) then + vec1[old_i] := -1; + end if; + end for; + end if; + end if; + + if Flags.isSet(Flags.BLT_DUMP) then + print("\n###################################\n" + + "INCLUDES FORCED STATESELECT.NEVER()\n" + + "###################################\n"); + end if; + + (dummyStates,states) := checkAssignment(1,nv,vec1,vars); + end if; +end forceStateSelectNever; + protected function selectBlock input list comp; input Integer ne; @@ -3253,29 +3367,29 @@ algorithm end match; end varStateSelectNever; -protected function incidenceMatrixfromEnhanced2 +protected function incidenceMatrixfromEnhancedStrict "author: Frenkel TUD 2012-11 converts an AdjacencyMatrixEnhanced into a IncidenceMatrix" input BackendDAE.AdjacencyMatrixEnhanced me; input BackendDAE.Variables vars; output BackendDAE.IncidenceMatrix m; algorithm - m := Array.map1(me,incidenceMatrixElementfromEnhanced2,vars); -end incidenceMatrixfromEnhanced2; + m := Array.map1(me,incidenceMatrixElementfromEnhancedStrict,vars); +end incidenceMatrixfromEnhancedStrict; -protected function incidenceMatrixElementfromEnhanced2 +protected function incidenceMatrixElementfromEnhancedStrict "author: Frenkel TUD 2012-11 - helper for incidenceMatrixfromEnhanced2" + helper for incidenceMatrixfromEnhancedStrict" input BackendDAE.AdjacencyMatrixElementEnhanced iRow; input BackendDAE.Variables vars; output BackendDAE.IncidenceMatrixElement oRow; algorithm - oRow := List.fold1(iRow, incidenceMatrixElementElementfromEnhanced2, vars, {}); + oRow := List.fold1(iRow, incidenceMatrixElementElementfromEnhancedStrict, vars, {}); oRow := List.map(oRow,intAbs); oRow := listReverse(oRow); -end incidenceMatrixElementfromEnhanced2; +end incidenceMatrixElementfromEnhancedStrict; -protected function incidenceMatrixElementElementfromEnhanced2 +protected function incidenceMatrixElementElementfromEnhancedStrict "author: Frenkel TUD 2012-11 converts an AdjacencyMatrix entry into a IncidenceMatrix entry" input tuple inTpl; @@ -3289,32 +3403,117 @@ algorithm case ((i,BackendDAE.SOLVABILITY_CONSTONE(),_),_,_) then i::iRow; case ((i,BackendDAE.SOLVABILITY_CONST(),_),_,_) then i::iRow; case ((i,BackendDAE.SOLVABILITY_PARAMETER(b=true),_),_,_) then i::iRow; -// case ((i,BackendDAE.SOLVABILITY_PARAMETER(b=false)),_,_) then incidenceMatrixElementElementfromEnhanced2_1(i,vars,iRow); -// case ((i,BackendDAE.SOLVABILITY_LINEAR(b=_)),_,_) then incidenceMatrixElementElementfromEnhanced2_1(i,vars,iRow); -// case ((i,BackendDAE.SOLVABILITY_NONLINEAR()),_,_) then incidenceMatrixElementElementfromEnhanced2_1(i,vars,iRow); -// case ((i,BackendDAE.SOLVABILITY_NONLINEAR()),_,_) then iRow; else iRow; end match; -end incidenceMatrixElementElementfromEnhanced2; - -// protected function incidenceMatrixElementElementfromEnhanced2_1 -// input Integer i; -// input BackendDAE.Variables vars; -// input list iRow; -// output list oRow; -// protected -// BackendDAE.Var v; -// DAE.StateSelect s; -// Integer si; -// Boolean b; -// algorithm -// v := BackendVariable.getVarAt(vars,intAbs(i)); -// s := BackendVariable.varStateSelect(v); -// si := BackendVariable.stateSelectToInteger(s); -// // oRow := List.consOnTrue(intLt(si,0),i,iRow); -// b := BackendVariable.isStateVar(v); -// oRow := List.consOnTrue(intLt(si,0) or not b,i,iRow); -// end incidenceMatrixElementElementfromEnhanced2_1; +end incidenceMatrixElementElementfromEnhancedStrict; + +protected function incidenceMatrixfromEnhanced +"author: kabdelhak FHB 2019-08 + converts an AdjacencyMatrixEnhanced into a IncidenceMatrix" + input BackendDAE.AdjacencyMatrixEnhanced me; + input BackendDAE.Variables vars; + input BackendDAE.StateOrder so; + output BackendDAE.IncidenceMatrix m; +algorithm + m := Array.map1(me,incidenceMatrixElementfromEnhanced,(vars,so)); +end incidenceMatrixfromEnhanced; + +protected function incidenceMatrixElementfromEnhanced +"author: kabdelhak FHB 2019-08 + helper for incidenceMatrixfromEnhanced" + input BackendDAE.AdjacencyMatrixElementEnhanced iRow; + input tuple tpl; + output BackendDAE.IncidenceMatrixElement oRow; +algorithm + oRow := List.fold1(iRow, incidenceMatrixElementElementfromEnhanced, tpl, {}); + oRow := List.map(oRow,intAbs); + oRow := listReverse(oRow); +end incidenceMatrixElementfromEnhanced; + +protected function incidenceMatrixfromEnhancedPartial +"author: kabdelhak FHB 2019-08 + converts an AdjacencyMatrixEnhanced into a IncidenceMatrix. + does not allow unmatching of stateSelect.never variables." + input BackendDAE.AdjacencyMatrixEnhanced me; + input BackendDAE.Variables vars; + input BackendDAE.Variables neverVars; + input array ass; + input BackendDAE.StateOrder so; + output BackendDAE.IncidenceMatrix m; +algorithm + m := Array.map1Ind(me,incidenceMatrixElementfromEnhancedPartial,(vars,neverVars,ass,so)); +end incidenceMatrixfromEnhancedPartial; + +protected function incidenceMatrixElementfromEnhancedPartial +"author: kabdelhak FHB 2019-08 + helper for incidenceMatrixfromEnhanced + does not allow unmatching of stateSelect.never variables." + input BackendDAE.AdjacencyMatrixElementEnhanced iRow; + input Integer index; + input tuple,BackendDAE.StateOrder> varsAssTpl; + output BackendDAE.IncidenceMatrixElement oRow; +protected + BackendDAE.Variables vars; + BackendDAE.Variables neverVars; + array ass; + BackendDAE.StateOrder so; +algorithm + (vars, neverVars, ass, so) := varsAssTpl; + if intEq(ass[index], -1) or (not BackendVariable.varStateSelectNever(BackendVariable.getVarAt(vars,ass[index]))) then + oRow := List.fold1(iRow, incidenceMatrixElementElementfromEnhanced, (neverVars, so), {}); + oRow := List.map(oRow,intAbs); + oRow := listReverse(oRow); + else + oRow := {}; + end if; +end incidenceMatrixElementfromEnhancedPartial; + +protected function incidenceMatrixElementElementfromEnhanced +"author: kabdelhak FHB 2019-08 + converts an AdjacencyMatrix entry into an IncidenceMatrix entry + ToDo: remove nonlinear (?)" + input tuple inTpl; + input tuple tpl; + input list iRow; + output list oRow; +algorithm + oRow := match inTpl + local Integer i; + case (i,BackendDAE.SOLVABILITY_SOLVED(),_) then i::iRow; + case (i,BackendDAE.SOLVABILITY_CONSTONE(),_) then i::iRow; + case (i,BackendDAE.SOLVABILITY_CONST(),_) then i::iRow; + case (i,BackendDAE.SOLVABILITY_PARAMETER(b=true),_) then i::iRow; + case (i,BackendDAE.SOLVABILITY_PARAMETER(b=false),_) then incidenceMatrixElementElementfromEnhanced_1(i,tpl,iRow); + case (i,BackendDAE.SOLVABILITY_LINEAR(b=true),_) then incidenceMatrixElementElementfromEnhanced_1(i,tpl,iRow); + case (i,BackendDAE.SOLVABILITY_NONLINEAR(),_) then incidenceMatrixElementElementfromEnhanced_1(i,tpl,iRow); + else iRow; + end match; +end incidenceMatrixElementElementfromEnhanced; + +protected function incidenceMatrixElementElementfromEnhanced_1 + "author: kabdelhak FHB 2019-08 + Only append if stateSelect.never and the state is not the derivative of another state + - ticket #5459 + Non-strict rules can only be applied on stateSelect.never, since they are forced in by the user. + If they are derivatives of another state, this rule does not apply, since (theoretically) the derivative + of the other state is actually what is simulated and the original state with stateSelect.never is not + considered anyways." + input Integer i; + input tuple tpl; + input list iRow; + output list oRow; +protected + BackendDAE.Variables vars; + BackendDAE.StateOrder so; + HashTable3.HashTable ht; + BackendDAE.Var v; + Boolean b; +algorithm + (vars, so as BackendDAE.STATEORDER(invHashTable = ht)) := tpl; + v := BackendVariable.getVarAt(vars,intAbs(i)); + b := BackendVariable.varStateSelectNever(v) and not BaseHashTable.hasKey(BackendVariable.varCref(v), ht); + oRow := List.consOnTrue(b,i,iRow); +end incidenceMatrixElementElementfromEnhanced_1; protected function checkAssignment "author: Frenkel TUD 2012-05 @@ -3944,7 +4143,7 @@ algorithm eqns := BackendEquation.listEquation(eqnslst); syst := BackendDAEUtil.createEqSystem(vars, eqns); (me, _, mapEqnIncRow, mapIncRowEqn) := BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(syst, shared, false); - m := incidenceMatrixfromEnhanced2(me, vars); + m := incidenceMatrixfromEnhancedStrict(me, vars); // match the equations, umatched are constrained equations nv := BackendVariable.varsSize(vars); ne := BackendEquation.equationArraySize(eqns); diff --git a/OMCompiler/Compiler/BackEnd/Matching.mo b/OMCompiler/Compiler/BackEnd/Matching.mo index a2468760dba..4588dd6535c 100644 --- a/OMCompiler/Compiler/BackEnd/Matching.mo +++ b/OMCompiler/Compiler/BackEnd/Matching.mo @@ -55,8 +55,10 @@ import Differentiate; import DumpGraphML; import ElementSource; import Error; +import Expression; import Flags; import IndexReduction; +import Inline; import List; import MetaModelica.Dangerous; import Util; @@ -82,7 +84,7 @@ algorithm end PerfectMatching; public function RegularMatching " - This function returns at least a partial matching for singular systems. + This function returns at least a partial matching for singular systems, starting from scratch. Unmatched nodes are represented by -1." input BackendDAE.IncidenceMatrix m; input Integer nVars; @@ -90,6 +92,21 @@ public function RegularMatching " output array ass1 "eqn := ass1[var]"; output array ass2 "var := ass2[eqn]"; output Boolean outPerfectMatching=true; +algorithm + ass2 := arrayCreate(nEqns, -1); + ass1 := arrayCreate(nVars, -1); + (ass1, ass2, outPerfectMatching) := ContinueMatching(m, nVars, nEqns, ass1, ass2); +end RegularMatching; + +public function ContinueMatching " + This function returns at least a partial matching for singular systems. + Unmatched nodes are represented by -1." + input BackendDAE.IncidenceMatrix m; + input Integer nVars; + input Integer nEqns; + input output array ass1 "eqn := ass1[var]"; + input output array ass2 "var := ass2[eqn]"; + output Boolean outPerfectMatching=true; protected Integer i, j; array eMark, vMark; @@ -115,7 +132,7 @@ algorithm end if; i := i+1; end while; -end RegularMatching; +end ContinueMatching; public function BBMatching input BackendDAE.EqSystem inSys; @@ -5590,25 +5607,27 @@ algorithm then (ass1,ass2,isyst,ishared,inArg); case ({},false,BackendDAE.EQSYSTEM(m=SOME(m),mT=SOME(mt)),_) - equation + algorithm matchingExternalsetIncidenceMatrix(nv,ne,m); BackendDAEEXT.matching(nv,ne,algIndx,cheapMatching,1.0,clearMatching); BackendDAEEXT.getAssignment(ass1,ass2); - unmatched1 = getUnassigned(ne, ass1, {}); + unmatched1 := getUnassigned(ne, ass1, {}); //BackendDump.dumpEqSystem(isyst, "EQSYS"); if Flags.isSet(Flags.BLT_DUMP) and Flags.isSet(Flags.GRAPHML) then BackendDump.dumpBipartiteGraphEqSystem(isyst, ishared, "BeforMatching_"+intString(arrayLength(m))+"_unmatched "+intString(listLength(unmatched1))); end if; if Flags.isSet(Flags.BLT_DUMP) and listLength(unmatched1) > 0 then print("unmatched equations: "+stringDelimitList(List.map(unmatched1,intString),", ")+"\n\n"); end if; // remove some edges which do not have to be traversed when finding the MSSS - m1 = arrayCopy(m); - m1t = arrayCopy(mt); - (m1,m1t) = removeEdgesForNoDerivativeFunctionInputs(m1,m1t,isyst,ishared); - (m1,m1t) = removeEdgesToDiscreteEquations(m1,m1t,isyst,ishared); + m1 := arrayCopy(m); + m1t := arrayCopy(mt); + (m1,m1t) := removeEdgesForNoDerivativeFunctionInputs(m1,m1t,isyst,ishared); + (m1,m1t) := removeEdgesToDiscreteEquations(m1,m1t,isyst,ishared); - meqns1 = getEqnsforIndexReduction(unmatched1,ne,m1,m1t,ass1,ass2,inArg); + meqns1 := getEqnsforIndexReduction(unmatched1,ne,m1,m1t,ass1,ass2,inArg); if not listEmpty(meqns1) then - sanityCheckArtificialStates(isyst, ishared, meqns1, ass2, ass1, inArg); + (syst, meqns1, ass1_1, ass2_1) := sanityCheckArtificialStates(isyst, ishared, nv, ne, meqns1, ass1, ass2, algIndx, cheapMatching, clearMatching, inArg); + else + (syst, meqns1, ass1_1, ass2_1) := (isyst, meqns1, ass1, ass2); end if; if Flags.isSet(Flags.BLT_DUMP) and listLength(meqns1) > 0 then print("Index Reduction neccessary!\n" @@ -5619,19 +5638,19 @@ algorithm //if listLength(List.flatten(meqns1)) >= 5 then meqs_short = List.firstN(List.flatten(meqns1),5); else meqs_short = List.flatten(meqns1); end if; //BackendDump.dumpBipartiteGraphEqSystem(isyst,ishared,"MSSS_"+stringDelimitList(List.map(meqs_short,intString),"_")); - (ass1_1,ass2_1,syst,shared,arg) = matchingExternal(meqns1,true,algIndx,-1,0,isyst,ishared,nv,ne,ass1,ass2,inMatchingOptions,sssHandler,inArg); + (ass1_1,ass2_1,syst,shared,arg) := matchingExternal(meqns1,true,algIndx,-1,0,syst,ishared,nv,ne,ass1_1,ass2_1,inMatchingOptions,sssHandler,inArg); then (ass1_1,ass2_1,syst,shared,arg); case (_::_,_,_,(BackendDAE.INDEX_REDUCTION(),_)) - equation - memsize = arrayLength(ass1); - (_,_,syst,shared,ass2_1,ass1_1,arg) = sssHandler(meqns,0,isyst,ishared,ass2,ass1,inArg); - ne_1 = BackendDAEUtil.systemSize(syst); - nv_1 = BackendVariable.daenumVariables(syst); - ass1_2 = assignmentsArrayExpand(ass1_1,ne_1,memsize,-1); - ass2_2 = assignmentsArrayExpand(ass2_1,nv_1,memsize,-1); - true = BackendDAEEXT.setAssignment(ne_1,nv_1,ass1_2,ass2_2); - (ass1_3,ass2_3,syst,shared,arg1) = matchingExternal({},false,algIndx,cheapMatching,clearMatching,syst,shared,nv_1,ne_1,ass1_2,ass2_2,inMatchingOptions,sssHandler,arg); + algorithm + memsize := arrayLength(ass1); + (_,_,syst,shared,ass2_1,ass1_1,arg) := sssHandler(meqns,0,isyst,ishared,ass2,ass1,inArg); + ne_1 := BackendDAEUtil.systemSize(syst); + nv_1 := BackendVariable.daenumVariables(syst); + ass1_2 := assignmentsArrayExpand(ass1_1,ne_1,memsize,-1); + ass2_2 := assignmentsArrayExpand(ass2_1,nv_1,memsize,-1); + true := BackendDAEEXT.setAssignment(ne_1,nv_1,ass1_2,ass2_2); + (ass1_3,ass2_3,syst,shared,arg1) := matchingExternal({},false,algIndx,cheapMatching,clearMatching,syst,shared,nv_1,ne_1,ass1_2,ass2_2,inMatchingOptions,sssHandler,arg); then (ass1_3,ass2_3,syst,shared,arg1); @@ -5645,58 +5664,131 @@ algorithm end matchingExternal; protected function sanityCheckArtificialStates - input BackendDAE.EqSystem syst; + "author: kabdelhak 2019-09 FHB + Checks if all equations that have to be differentiated for index reduction can + be differentiated by all artificial states. Reverts all those for which we + cannot derive to be algebraic variables. + Ticket: 5459 + ----------------------------- INFO ----------------------------- + Artificial states are those which do not naturally appear + differentiated in the system of DAEs, but have been forced + to be states with 'StateSelect.always' or 'StateSelect.prefer'. + ----------------------------------------------------------------" + input output BackendDAE.EqSystem syst; input BackendDAE.Shared shared; - input list> eqns; - input array ass1; - input array ass2; + input Integer nv; + input Integer ne; + input output list> eqns; + input output array ass1; + input output array ass2; + input Integer algIndx "Index of the algorithm, see BackendDAEEXT.matching"; + input Integer cheapMatching "Method for cheap Matching"; + input Integer clearMatching; input BackendDAE.StructurallySingularSystemHandlerArg arg; protected - BackendDAE.Variables orderedVars; - BackendDAE.EquationArray orderedEqs; - list> eqns_1, unassignedStates, unassignedEqs; - list flat_unassignedStates, flat_unassignedEqs; - list artificialStates = {}; + list> eqns_1, unassignedStates; + list flat_unassignedStates, flat_eqns, unmatched1; + list artificialStates = {}, undiffable_artificial = {}; list residuals, equations = {}; BackendDAE.Var var; + DAE.ComponentRef cr; DAE.Exp residualExp; + BackendDAE.AdjacencyMatrix m, mt, m1, m1t; + Boolean unique_flag; + String msg; algorithm - for i_l in eqns loop - for i in i_l loop - print(intString(i) + " "); - end for; - print("\n"); - end for; - BackendDAE.EQSYSTEM(orderedVars = orderedVars, orderedEqs = orderedEqs) := syst; - (eqns_1, unassignedStates, unassignedEqs, _) := IndexReduction.minimalStructurallySingularSystem(eqns, syst, shared, ass1, ass2, arg); - print("test 2\n"); - + try + // Get the information about MSSS from index reduction + (eqns_1, unassignedStates, _, _) := IndexReduction.minimalStructurallySingularSystem(eqns, syst, shared, ass2, ass1, arg); + else + if Flags.isSet(Flags.BLT_DUMP) then + singularSystemError(eqns, 0, syst, shared, ass1, ass2, arg); + end if; + fail(); + end try; flat_unassignedStates := List.flatten(unassignedStates); + // Collect artificial states for state in flat_unassignedStates loop - var := BackendVariable.getVarAt(orderedVars, state); + var := BackendVariable.getVarAt(syst.orderedVars, state); if BackendVariable.isArtificialState(var) and not listMember(var, artificialStates) then artificialStates := var :: artificialStates; end if; end for; - BackendDump.dumpVarList(artificialStates, "artificial"); - flat_unassignedEqs := List.flatten(unassignedEqs); - for eqn in flat_unassignedEqs loop - equations := BackendEquation.get(orderedEqs, eqn) :: equations; + flat_eqns := List.flatten(eqns_1); // use eqns_1 (without discrete) + for eqn in flat_eqns loop + try + equations := BackendEquation.get(syst.orderedEqs, eqn) :: equations; + else + // equation can't be found -- array index to normal index? + end try; end for; for var in artificialStates loop + unique_flag := false; for eqn in equations loop try residuals := BackendEquation.equationToScalarResidualForm(eqn, shared.functionTree); for res in residuals loop BackendDAE.RESIDUAL_EQUATION(exp = residualExp) := res; - //(_, _) := Differentiate.differentiateEquationFragile(eqn, BackendVariable.varCref(var), BackendDAE.emptyInputData, BackendDAE.SIMPLE_DIFFERENTIATION(), shared.functionTree); - _ := Differentiate.differentiateExpSolve(residualExp, BackendVariable.varCref(var), SOME(shared.functionTree)); + cr := BackendVariable.varCref(var); + // Only check if the equation actually contains the cref + if Expression.expHasCref(residualExp, cr) then + (residualExp,_,_) := Inline.forceInlineExp(residualExp,(SOME(shared.functionTree),{DAE.NORM_INLINE(),DAE.DEFAULT_INLINE()}),DAE.emptyElementSource); + // Check if the equation can be differentiated by the artificial state + _ := Differentiate.differentiateExpSolve(residualExp, cr, SOME(shared.functionTree)); + // Check if the equation contains smooth(0, cr) + false := Expression.expHasCrefInSmoothZero(residualExp, cr); + end if; end for; else - print("failed for " + BackendDump.varString(var) + ".\n"); + if Flags.isSet(Flags.BLT_DUMP) then + print("### The Equation ### \n" + BackendDump.equationString(eqn) + + "\n\n--- could not be differentiated for artificial variable ---\n " + + ComponentReference.printComponentRefStr(var.varName) + ".\n\n"); + end if; + // revert the varKind from STATE to VARIABLE + if not unique_flag then + undiffable_artificial := BackendVariable.setVarKind(var, BackendDAE.VARIABLE()) :: undiffable_artificial; + unique_flag := true; + end if; end try; end for; end for; + + if not listEmpty(undiffable_artificial) then + syst.orderedVars := BackendVariable.addVars(undiffable_artificial, syst.orderedVars); + (syst, _, _, _, _) := BackendDAEUtil.getIncidenceMatrixScalar(syst, BackendDAE.SOLVABLE(), SOME(shared.functionTree)); + + if isSome(syst.m) and isSome(syst.mT) then + SOME(m) := syst.m; + SOME(mt) := syst.mT; + // It would be great to replace full matching by continue matching, + // to reuse old information, but somehow it does something different. + // (ass1, ass2) := ContinueMatching(m, nv, ne, ass1, ass2); + matchingExternalsetIncidenceMatrix(nv, ne, m); + BackendDAEEXT.matching(nv, ne, algIndx, cheapMatching, 1.0, clearMatching); + BackendDAEEXT.getAssignment(ass1, ass2); + unmatched1 := getUnassigned(ne, ass1, {}); + m1 := arrayCopy(m); + m1t := arrayCopy(mt); + (m1,m1t) := removeEdgesForNoDerivativeFunctionInputs(m1, m1t, syst, shared); + (m1,m1t) := removeEdgesToDiscreteEquations(m1, m1t, syst, shared); + eqns := getEqnsforIndexReduction(unmatched1, ne, m1, m1t, ass1, ass2, arg); + end if; + + if Flags.isSet(Flags.BLT_DUMP) then + print("----------------------------- INFO -----------------------------\n" + + " Artificial states are those which do not naturally appear\n" + + " differentiated in the system of DAEs, but have been forced\n" + + " to be states with 'StateSelect.always' or 'StateSelect.prefer'.\n" + + " The ones mentioned above will be treated as if they had \n" + + "'StateSelect.default'.\n" + + "----------------------------------------------------------------\n\n"); + end if; + + msg := System.gettext(BackendDump.varListStringShort(undiffable_artificial,"They will be treated as if they had stateSelect=StateSelect.default") + + "Please use -d=bltdump for more information.\n"); + Error.addMessage(Error.STATE_STATESELECT_PREFER_REVERT, {msg}); + end if; end sanityCheckArtificialStates; protected function removeEdgesToDiscreteEquations"" @@ -5924,6 +6016,24 @@ algorithm end match; end getUnassigned; +public function anyUnassigned + "author: kabdelhak FHB 2019-08 + returns true if any assignment is -1." + input Integer ne; + input array ass; + output Boolean b; +algorithm + b := match(ne) + local + list unassigned; + case 0 + then false; + case _ guard(intLt(ass[ne],1)) + then true; + else anyUnassigned(ne-1,ass); + end match; +end anyUnassigned; + public function getAssignedArray "author: lochel" input array ass; output array outIsAssigned; diff --git a/OMCompiler/Compiler/BackEnd/SymbolicJacobian.mo b/OMCompiler/Compiler/BackEnd/SymbolicJacobian.mo index f648aae6783..3b7d75129ae 100644 --- a/OMCompiler/Compiler/BackEnd/SymbolicJacobian.mo +++ b/OMCompiler/Compiler/BackEnd/SymbolicJacobian.mo @@ -3095,6 +3095,7 @@ algorithm end if; backendDAE := BackendDAEUtil.transformBackendDAE(backendDAE, SOME((BackendDAE.NO_INDEX_REDUCTION(), BackendDAE.EXACT())), NONE(), NONE()); + BackendDAE.DAE({BackendDAE.EQSYSTEM(orderedVars = dependentVars)}, BackendDAE.SHARED(globalKnownVars = globalKnownVars)) := backendDAE; // prepare creation of symbolic jacobian @@ -3114,6 +3115,7 @@ algorithm outJacobian := BackendDAE.GENERIC_JACOBIAN(symJacBDAE, sparsePattern, sparseColoring); outShared := BackendDAEUtil.setSharedFunctionTree(inShared, funcs); else + if Flags.isSet(Flags.JAC_DUMP) then Error.addInternalError("function getSymbolicJacobian failed", sourceInfo()); end if; diff --git a/OMCompiler/Compiler/FrontEnd/Expression.mo b/OMCompiler/Compiler/FrontEnd/Expression.mo index ba37e233674..89c9ffc1a98 100644 --- a/OMCompiler/Compiler/FrontEnd/Expression.mo +++ b/OMCompiler/Compiler/FrontEnd/Expression.mo @@ -6808,6 +6808,39 @@ algorithm end match; end expHasCrefInIfWork; +public function expHasCrefInSmoothZero + "author: kabdelhak FHB 2019-09 + Returns true if the expression contains a function call of the form + smooth(0, cr). Used to determine if an artificial state should + not be differentiated and reverted to be an algebraic variable + instead. " + input DAE.Exp exp; + input DAE.ComponentRef cr; + output Boolean b; +algorithm + (_, (_, b)) := traverseExpBottomUp(exp, expHasCrefInSmoothZeroWork, (cr, false)); +end expHasCrefInSmoothZero; + +protected function expHasCrefInSmoothZeroWork + "author: kabdelhak FHB 2019-09 + Work function to detect smooth(0, cr). + Q: Should it also return true if the smooth call expression + contains the cref in any way not only for direct identity?" + input output DAE.Exp exp; + input output tuple tpl; +algorithm + tpl := match (exp, tpl) + local + DAE.ComponentRef cr, sCr; + Boolean b; + case (DAE.CALL(path = Absyn.IDENT(name = "smooth"), expLst = {DAE.ICONST(0), DAE.CREF(componentRef = sCr)}), (cr, false)) + algorithm + b := ComponentReference.crefEqual(sCr, cr); + //expHasCref(e, cr); use this instead for full check? + then (cr, b); + else tpl; + end match; +end expHasCrefInSmoothZeroWork; public function traverseCrefsFromExp " Author: Frenkel TUD 2011-05, traverses all ComponentRef from an Expression." diff --git a/OMCompiler/Compiler/Util/Array.mo b/OMCompiler/Compiler/Util/Array.mo index 930e30663a3..d3171718239 100644 --- a/OMCompiler/Compiler/Util/Array.mo +++ b/OMCompiler/Compiler/Util/Array.mo @@ -248,6 +248,41 @@ algorithm end if; end map1; +function map1Ind + "Takes an array, an extra arguments, and a function over the elements of the + array, which is applied to each element. The index is passed as an extra + argument. The updated elements will form a new + array, leaving the original array unchanged." + input array inArray; + input FuncType inFunc; + input ArgT inArg; + output array outArray; + + partial function FuncType + input TI inElement; + input Integer index; + input ArgT inArg; + output TO outElement; + end FuncType; +protected + Integer len = arrayLength(inArray); + TO res; +algorithm + // If the array is empty, use list transformations to fix the types! + if len == 0 then + outArray := listArray({}); + else + // If the array isn't empty, use the first element to create the new array. + res := inFunc(arrayGetNoBoundsChecking(inArray, 1), 1, inArg); + outArray := arrayCreateNoInit(len, res); + arrayUpdate(outArray, 1, res); + + for i in 2:len loop + arrayUpdate(outArray, i, inFunc(arrayGetNoBoundsChecking(inArray, i), i, inArg)); + end for; + end if; +end map1Ind; + function map0 "Applies a non-returning function to all elements in an array." input array inArray; diff --git a/OMCompiler/Compiler/Util/Error.mo b/OMCompiler/Compiler/Util/Error.mo index 824c5d3b9b4..65a49577f7d 100644 --- a/OMCompiler/Compiler/Util/Error.mo +++ b/OMCompiler/Compiler/Util/Error.mo @@ -1008,34 +1008,38 @@ public constant Message MOVING_PARAMETER_BINDING_TO_INITIAL_EQ_SECTION = MESSAGE Util.gettext("Moving binding to initial equation section and setting fixed attribute of %s to false.")); public constant Message MIXED_DETERMINED = MESSAGE(584, SYMBOLIC(), ERROR(), Util.gettext("The initialization problem of given system is mixed-determined. It is under- as well as overdetermined and the mixed-determination-index is too high. [index > %s]\nPlease checkout the option \"--maxMixedDeterminedIndex\" to simulate with a higher threshold or consider changing some initial equations, fixed variables and start values.")); -public constant Message STACK_OVERFLOW_DETAILED = MESSAGE(584, SCRIPTING(), ERROR(), +public constant Message STACK_OVERFLOW_DETAILED = MESSAGE(585, SCRIPTING(), ERROR(), Util.gettext("Stack overflow occurred while evaluating %s:\n%s")); -public constant Message NF_VECTOR_INVALID_DIMENSIONS = MESSAGE(585, TRANSLATION(), ERROR(), +public constant Message NF_VECTOR_INVALID_DIMENSIONS = MESSAGE(586, TRANSLATION(), ERROR(), Util.gettext("Invalid dimensions %s in %s, no more than one dimension may have size > 1.")); -public constant Message NF_ARRAY_TYPE_MISMATCH = MESSAGE(586, TRANSLATION(), ERROR(), +public constant Message NF_ARRAY_TYPE_MISMATCH = MESSAGE(587, TRANSLATION(), ERROR(), Util.gettext("Array types mismatch. Argument %s (%s) has type %s whereas previous arguments have type %s.")); -public constant Message NF_DIFFERENT_NUM_DIM_IN_ARGUMENTS = MESSAGE(587, TRANSLATION(), ERROR(), +public constant Message NF_DIFFERENT_NUM_DIM_IN_ARGUMENTS = MESSAGE(588, TRANSLATION(), ERROR(), Util.gettext("Different number of dimensions (%s) in arguments to %s.")); -public constant Message NF_CAT_WRONG_DIMENSION = MESSAGE(588, TRANSLATION(), ERROR(), +public constant Message NF_CAT_WRONG_DIMENSION = MESSAGE(589, TRANSLATION(), ERROR(), Util.gettext("The first argument of cat characterizes an existing dimension in the other arguments (1..%s), but got dimension %s.")); -public constant Message NF_CAT_FIRST_ARG_EVAL = MESSAGE(589, TRANSLATION(), ERROR(), +public constant Message NF_CAT_FIRST_ARG_EVAL = MESSAGE(590, TRANSLATION(), ERROR(), Util.gettext("The first argument of cat must be possible to evaluate during compile-time. Expression %s has variability %s.")); -public constant Message COMMA_OPERATOR_DIFFERENT_SIZES = MESSAGE(590, TRANSLATION(), ERROR(), +public constant Message COMMA_OPERATOR_DIFFERENT_SIZES = MESSAGE(591, TRANSLATION(), ERROR(), Util.gettext("Arguments of concatenation comma operator have different sizes for the first dimension: %s has dimension %s and %s has dimension %s.")); -public constant Message NON_STATE_STATESELECT_ALWAYS = MESSAGE(591, SYMBOLIC(), WARNING(), +public constant Message NON_STATE_STATESELECT_ALWAYS = MESSAGE(592, SYMBOLIC(), WARNING(), Util.gettext("Variable %s has attribute stateSelect=StateSelect.always, but was selected as a continuous variable.")); -public constant Message STATE_STATESELECT_NEVER = MESSAGE(592, SYMBOLIC(), WARNING(), +public constant Message STATE_STATESELECT_NEVER = MESSAGE(593, SYMBOLIC(), WARNING(), Util.gettext("Variable %s has attribute stateSelect=StateSelect.never, but was selected as a state")); -public constant Message FUNCTION_HIGHER_VARIABILITY_BINDING = MESSAGE(593, TRANSLATION(), WARNING(), +public constant Message FUNCTION_HIGHER_VARIABILITY_BINDING = MESSAGE(594, TRANSLATION(), WARNING(), Util.gettext("Component ‘%s’ of variability %s has binding %s of higher variability %s.")); -public constant Message OCG_MISSING_BRANCH = MESSAGE(594, TRANSLATION(), WARNING(), +public constant Message OCG_MISSING_BRANCH = MESSAGE(595, TRANSLATION(), WARNING(), Util.gettext("Connections.rooted(%s) needs exactly one statement Connections.branch(%s, B.R) involving %s but we found none in the graph. Run with -d=cgraphGraphVizFile to debug")); -public constant Message UNBOUND_PARAMETER_EVALUATE_TRUE = MESSAGE(594, TRANSLATION(), WARNING(), +public constant Message UNBOUND_PARAMETER_EVALUATE_TRUE = MESSAGE(596, TRANSLATION(), WARNING(), Util.gettext("Parameter %s has annotation(Evaluate=true) and no binding.")); -public constant Message FMI_URI_RESOLVE = MESSAGE(595, TRANSLATION(), WARNING(), +public constant Message FMI_URI_RESOLVE = MESSAGE(597, TRANSLATION(), WARNING(), Util.gettext("Could not resolve URI (%s) at compile-time; copying all loaded packages into the FMU")); -public constant Message PATTERN_MIXED_POS_NAMED = MESSAGE(596, TRANSLATION(), WARNING(), +public constant Message PATTERN_MIXED_POS_NAMED = MESSAGE(598, TRANSLATION(), WARNING(), Util.gettext("Call to %s contains mixed positional and named arguments.")); +public constant Message STATE_STATESELECT_NEVER_FORCED = MESSAGE(599, TRANSLATION(), WARNING(), + Util.gettext("Following variables have attribute stateSelect=StateSelect.never, but cant be statically chosen. %s")); +public constant Message STATE_STATESELECT_PREFER_REVERT = MESSAGE(600, TRANSLATION(), WARNING(), + Util.gettext("Some equations could not be differentiated for following variables having attribute stateSelect=StateSelect.prefer. %s")); public constant Message MATCH_SHADOWING = MESSAGE(5001, TRANSLATION(), ERROR(), Util.gettext("Local variable '%s' shadows another variable.")); diff --git a/OMCompiler/Compiler/Util/List.mo b/OMCompiler/Compiler/Util/List.mo index e29bcb59e6a..b73b0e836c2 100644 --- a/OMCompiler/Compiler/Util/List.mo +++ b/OMCompiler/Compiler/Util/List.mo @@ -4502,12 +4502,12 @@ public function flatten of the sublists. O(len(outList)) Example: flatten({{1, 2}, {3, 4, 5}, {6}, {}}) => {1, 2, 3, 4, 5, 6}" input list> inList; - output list outList = listAppend(lst for lst in listReverse(inList)); + output list outList = if listEmpty(inList) then {} elseif intEq(listLength(inList), 1) then first(inList) else listAppend(lst for lst in listReverse(inList)); end flatten; public function flattenReverse input list> inList; - output list outList = listAppend(lst for lst in inList); + output list outList = if listEmpty(inList) then {} elseif intEq(listLength(inList), 1) then first(inList) else listAppend(lst for lst in inList); end flattenReverse; public function thread @@ -4567,6 +4567,27 @@ algorithm outTuples := list((e1, e2) threaded for e1 in inList1, e2 in inList2); end threadTuple; +public function zip + "Takes two lists and returns a list of two-element tuples contaning the + elements in the same order. Fails if the lists are not of the same length. + Example: zip({1, 3}, {2, 4}) => {(1, 2), (3, 4)}" + input list inList1; + input list inList2; + output list> outTuples = {}; +protected + list dummyList = inList2; + T2 t2; +algorithm + if intEq(listLength(inList1),listLength(inList2)) then + for t1 in inList1 loop + t2::dummyList := dummyList; + outTuples := (t1, t2)::outTuples; + end for; + else fail(); + end if; + outTuples := listReverse(outTuples); +end zip; + public function unzip "Takes a list of two-element tuples and splits the tuples into two separate lists. Example: unzip({(1, 2), (3, 4)}) => ({1, 3}, {2, 4})"