diff --git a/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBCausalize.mo b/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBCausalize.mo index 1a37657973..3af66ad61b 100644 --- a/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBCausalize.mo +++ b/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBCausalize.mo @@ -254,8 +254,8 @@ protected // ################################################# // Phase I: match initial equations <-> unfixable vars // ################################################# - vn := UnorderedMap.subSet(system.unknowns.map, list(BVariable.getVarName(var) for var in unfixable)); - en := UnorderedMap.subSet(system.equations.map, list(Equation.getEqnName(eqn) for eqn in initials)); + vn := UnorderedMap.subMap(system.unknowns.map, list(BVariable.getVarName(var) for var in unfixable)); + en := UnorderedMap.subMap(system.equations.map, list(Equation.getEqnName(eqn) for eqn in initials)); adj_matching := Adjacency.Matrix.fromFull(full, vn, en, system.equations, NBAdjacency.MatrixStrictness.MATCHING); matching := Matching.regular(NBMatching.EMPTY_MATCHING, adj_matching, true, true); @@ -266,7 +266,7 @@ protected vo := vn; eo := en; vn := UnorderedMap.new(ComponentRef.hash, ComponentRef.isEqual); - en := UnorderedMap.subSet(system.equations.map, list(Equation.getEqnName(eqn) for eqn in simulation)); + en := UnorderedMap.subMap(system.equations.map, list(Equation.getEqnName(eqn) for eqn in simulation)); (adj_matching, full) := Adjacency.Matrix.expand(adj_matching, full, vo, vn, eo, en, system.unknowns, system.equations); matching := Matching.regular(matching, adj_matching, true, true); @@ -276,7 +276,7 @@ protected // ################################################# vo := UnorderedMap.merge(vo, vn, sourceInfo()); eo := UnorderedMap.merge(eo, en, sourceInfo()); - vn := UnorderedMap.subSet(system.unknowns.map, list(BVariable.getVarName(var) for var in fixable)); + vn := UnorderedMap.subMap(system.unknowns.map, list(BVariable.getVarName(var) for var in fixable)); en := UnorderedMap.new(ComponentRef.hash, ComponentRef.isEqual); (adj_matching, full) := Adjacency.Matrix.expand(adj_matching, full, vo, vn, eo, en, system.unknowns, system.equations); (matching, adj_matching, full, variables, equations, funcTree, varData, eqData) := Matching.singular(matching, adj_matching, full, system.unknowns, system.equations, funcTree, varData, eqData, system.systemType, false, true, false); diff --git a/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBResolveSingularities.mo b/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBResolveSingularities.mo index f4f92b1869..808906b246 100644 --- a/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBResolveSingularities.mo +++ b/OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBResolveSingularities.mo @@ -375,7 +375,7 @@ public // update adjacency matrices vn := UnorderedMap.new(ComponentRef.hash, ComponentRef.isEqual); - en := UnorderedMap.subSet(equations.map, list(Equation.getEqnName(eqn) for eqn in start_eqns)); + en := UnorderedMap.subMap(equations.map, list(Equation.getEqnName(eqn) for eqn in start_eqns)); (adj, full) := Adjacency.Matrix.expand(adj, full, vo, vn, eo, en, variables, equations); if Flags.isSet(Flags.INITIALIZATION) then diff --git a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo index 6e56fe3912..798bce962d 100644 --- a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo +++ b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBSolve.mo @@ -200,6 +200,8 @@ public list> rest, sliced_eqns = {}; StrongComponent generic_comp; list entwined_slices = {}; + Tearing strict; + list tmp, inner_comps = {}; case StrongComponent.SINGLE_COMPONENT() algorithm (eqn, funcTree, solve_status, implicit_index) := solveSingleStrongComponent(Pointer.access(comp.eqn), Pointer.access(comp.var), funcTree, systemType, implicit_index, slicing_map); @@ -209,8 +211,13 @@ public (eqn, funcTree, solve_status, implicit_index) := solveMultiStrongComponent(Pointer.access(comp.eqn), comp.vars, funcTree, systemType, implicit_index, slicing_map); then ({StrongComponent.MULTI_COMPONENT(comp.vars, Pointer.create(eqn), solve_status)}, solve_status); - case StrongComponent.ALGEBRAIC_LOOP() algorithm - // do we need to do smth here? e.g. solve inner equations? call tearing from here? + case StrongComponent.ALGEBRAIC_LOOP(strict = strict) algorithm + for inner_comp in listReverse(arrayList(strict.innerEquations)) loop + (tmp, funcTree, implicit_index) := solveStrongComponent(inner_comp, funcTree, systemType, implicit_index, slicing_map); + inner_comps := listAppend(tmp, inner_comps); + end for; + strict.innerEquations := listArray(inner_comps); + comp.strict := strict; comp.status := Status.IMPLICIT; then ({comp}, Status.IMPLICIT); @@ -540,8 +547,8 @@ public else // If eqn is non-linear in cref if Flags.isSet(Flags.FAILTRACE) then - Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " cref: " - + ComponentRef.toString(fixed_cref) + " has to be solved implicitely in equation:\n" + Equation.toString(eqn)}); + Error.addCompilerWarning(getInstanceName() + " cref: " + ComponentRef.toString(fixed_cref) + + " has to be solved implicitely in equation:\n" + Equation.toString(eqn)); end if; invertRelation := false; status := Status.IMPLICIT; @@ -583,7 +590,6 @@ public end if; end solveIfBody; -protected function solveSimple input output Equation eqn; input ComponentRef cref; @@ -601,11 +607,13 @@ protected case Equation.WHEN_EQUATION() then (eqn, Status.EXPLICIT, false); // ToDo: more cases + // ToDo: tuples, record elements, array constructors else (eqn, Status.UNPROCESSED, false); end match; end solveSimple; +protected function solveSimpleLhsRhs input Expression lhs; input Expression rhs; diff --git a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBTearing.mo b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBTearing.mo index 7ef57454e5..05ce9198ee 100644 --- a/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBTearing.mo +++ b/OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBTearing.mo @@ -48,12 +48,13 @@ protected import Tearing = NBTearing; // NF imports - import ComponentRef = NFComponentRef; import NFFlatten.FunctionTree; import Variable = NFVariable; + import ComponentRef = NFComponentRef; // Backend imports import Adjacency = NBAdjacency; + import NBAdjacency.Solvability; import BEquation = NBEquation; import BJacobian = NBJacobian; import BVariable = NBVariable; @@ -159,9 +160,8 @@ public input output Integer index "current unique loop index"; input System.SystemType systemType = NBSystem.SystemType.ODE "system type"; protected - // dummy adjacency matrix, dont need it for tearingNone() + // dummy adjacency matrix, dont need it for implicit Adjacency.Matrix dummy = Adjacency.EMPTY(NBAdjacency.MatrixStrictness.FULL); - list funcs = {tearingNone, tearingFinalize}; StrongComponent new_comp; algorithm (comp, dummy, funcTree, index) := match comp @@ -174,10 +174,7 @@ public linear = false, mixed = false, status = NBSolve.Status.IMPLICIT); - for func in funcs loop - (new_comp, dummy, funcTree, index) := func(new_comp, dummy, funcTree, index, VariablePointers.empty(), EquationPointers.empty(), Pointer.create(0), systemType); - end for; - then (new_comp, dummy, funcTree, index); + then finalize(new_comp, dummy, funcTree, index, VariablePointers.empty(), EquationPointers.empty(), Pointer.create(0), systemType); case StrongComponent.MULTI_COMPONENT() algorithm new_comp := StrongComponent.ALGEBRAIC_LOOP( @@ -187,10 +184,7 @@ public linear = false, mixed = false, status = NBSolve.Status.IMPLICIT); - for func in funcs loop - (new_comp, dummy, funcTree, index) := func(new_comp, dummy, funcTree, index, VariablePointers.empty(), EquationPointers.empty(), Pointer.create(0), systemType); - end for; - then (new_comp, dummy, funcTree, index); + then finalize(new_comp, dummy, funcTree, index, VariablePointers.empty(), EquationPointers.empty(), Pointer.create(0), systemType); // do nothing otherwise else (comp, dummy, funcTree, index); @@ -214,10 +208,10 @@ public String flag = Flags.getConfigString(Flags.TEARING_METHOD); algorithm funcs := match flag - case "cellier" then {tearingNone, tearingFinalize}; - case "noTearing" then {tearingNone, tearingFinalize}; - case "omcTearing" then {tearingMinimal, tearingFinalize}; - case "minimalTearing" then {tearingMinimal, tearingFinalize}; + case "cellier" then {initialize, none, finalize}; + case "noTearing" then {initialize, none, finalize}; + case "omcTearing" then {initialize, minimal, finalize}; + case "minimalTearing" then {initialize, minimal, finalize}; /* ... New tearing modules have to be added here */ else fail(); end match; @@ -275,40 +269,39 @@ protected new_systems := listReverse(new_systems); end tearingTraverser; - function tearingFinalize extends Module.tearingInterface; + function initialize extends Module.tearingInterface; protected - list residual_comps; - Option jacobian; Tearing strict; - protected - list> tmp; - list>> acc = {}; - String tag; + list> vars_lst; + list> eqns_lst; + UnorderedSet vars_set "all loop vars, used to determine solvability"; + UnorderedMap v, e "all loop vars and equations map"; algorithm - (comp, index) := match comp + (comp, full, index) := match comp case StrongComponent.ALGEBRAIC_LOOP(strict = strict) algorithm - tag := if comp.linear then "_LS_JAC_" else "_NLS_JAC_"; - // create residual equations - residual_comps := list(StrongComponent.fromSolvedEquationSlice(eqn) for eqn in strict.residual_eqns); - // update jacobian to take slices (just to have correct inner variables and such) - (jacobian, funcTree) := BJacobian.nonlinear( - variables = VariablePointers.fromList(list(Slice.getT(var) for var in strict.iteration_vars)), - equations = EquationPointers.fromList(list(Slice.getT(eqn) for eqn in strict.residual_eqns)), - comps = listArray(residual_comps), - funcTree = funcTree, - name = System.System.systemTypeString(systemType) + tag + intString(index)); - strict.jac := jacobian; - comp.strict := strict; - if Flags.isSet(Flags.TEARING_DUMP) then - print(StrongComponent.toString(comp) + "\n"); - end if; - then (comp, index); - else (comp, index); + index := index + 1; + comp.idx := index; + + // get all loop variables and equations + vars_lst := list(Slice.getT(var) for var in strict.iteration_vars); + eqns_lst := list(Slice.getT(eqn) for eqn in strict.residual_eqns); + + // the set of all loop variables used to determine solvability + vars_set := UnorderedSet.fromList(list(BVariable.getVarName(var) for var in vars_lst), ComponentRef.hash, ComponentRef.isEqual); + + // the sets of discrete variables and discrete equations + v := UnorderedMap.subMap(variables.map, list(BVariable.getVarName(var) for var in vars_lst)); + e := UnorderedMap.subMap(equations.map, list(Equation.getEqnName(eqn) for eqn in eqns_lst)); + + // refining the adjacency matrix by updating solvability information using differentiation + (full, funcTree) := Adjacency.Matrix.refine(full, funcTree, v, e, variables, equations, vars_set); + comp.linear := checkLinearity(full, v, e); + then (comp, full, index); + else (comp, full, index); end match; - end tearingFinalize; + end initialize; - function tearingNone extends Module.tearingInterface; - // does nothing but set index and call the jacobian + function finalize extends Module.tearingInterface; protected list residual_comps; Option jacobian; @@ -316,12 +309,11 @@ protected protected list> tmp; list>> acc = {}; + String tag; algorithm - (comp, index) := match comp + comp := match comp case StrongComponent.ALGEBRAIC_LOOP(strict = strict) algorithm - index := index + 1; - comp.idx := index; - + // inline potential records for eqn in listReverse(strict.residual_eqns) loop tmp := Inline.inlineRecordSliceEquation(eqn, variables, eq_index, true); if listEmpty(tmp) then @@ -333,205 +325,107 @@ protected // create residual equations strict.residual_eqns := list(Slice.apply(eqn, function Equation.createResidual(new = true)) for eqn in List.flatten(acc)); + + tag := if comp.linear then "_LS_JAC_" else "_NLS_JAC_"; + // create residual equations + residual_comps := list(StrongComponent.fromSolvedEquationSlice(eqn) for eqn in strict.residual_eqns); + // update jacobian to take slices (just to have correct inner variables and such) + (jacobian, funcTree) := BJacobian.nonlinear( + variables = VariablePointers.fromList(list(Slice.getT(var) for var in strict.iteration_vars)), + equations = EquationPointers.fromList(list(Slice.getT(eqn) for eqn in strict.residual_eqns)), + comps = listArray(residual_comps), + funcTree = funcTree, + name = System.System.systemTypeString(systemType) + tag + intString(index)); + strict.jac := jacobian; comp.strict := strict; - then (comp, index); - else (comp, index); + if Flags.isSet(Flags.TEARING_DUMP) then + print(StrongComponent.toString(comp) + "\n"); + end if; + then comp; + else comp; end match; - end tearingNone; + end finalize; + + function none extends Module.tearingInterface; + // does nothing + end none; - function tearingMinimal extends Module.tearingInterface; - // only extracts discrete variables to be solved as inner equations and calls jacobian module + function minimal extends Module.tearingInterface; + // only extracts discrete variables to be solved as inner equations protected - list residual_comps; - Option jacobian; Tearing strict; list> vars_lst, cont_vars, disc_vars; list> eqns_lst, cont_eqns, disc_eqns; - list> iter_lst; list> residual_lst; - VariablePointers discreteVars; - EquationPointers eqns; Adjacency.Matrix adj; Matching matching; - list inner_comps, residual_comps; - UnorderedMap v, e; + list inner_comps; + UnorderedMap v, e "discrete variables and equations we have to refine"; algorithm - print("######## minimal ########\n"); - (comp, index) := match comp + comp := match comp case StrongComponent.ALGEBRAIC_LOOP(strict = strict) algorithm - + // split equations and variables for discretes and continuous vars_lst := list(Slice.getT(var) for var in strict.iteration_vars); eqns_lst := list(Slice.getT(eqn) for eqn in strict.residual_eqns); (cont_vars, disc_vars) := List.splitOnTrue(vars_lst, BVariable.isContinuous); (cont_eqns, disc_eqns) := List.splitOnTrue(eqns_lst, Equation.isContinuous); - print(List.toString(disc_vars, function BVariable.pointerToString(), "", "", "\n", "") + "\n"); - print(List.toString(disc_eqns, function Equation.pointerToString(str = ""), "", "", "\n", "") + "\n"); - v := UnorderedMap.subSet(variables.map, list(BVariable.getVarName(var) for var in disc_vars)); - e := UnorderedMap.subSet(equations.map, list(Equation.getEqnName(eqn) for eqn in disc_eqns)); - (full, funcTree) := Adjacency.Matrix.refine(full, funcTree, v, e, variables, equations); - - - /* - - // ToDo: if other tearing modules have been used before - // we should not only look in residuals and iteration arrays. - // have minimal tearing explicitely as starting method? - - // Equation attributes update! -> discrete initial etc are not mutually exclusive - // - - // split variables and equations in discrete and continuous - iter_lst := list(Slice.SLICE(var, {}) for var in cont_vars); - - if listEmpty(disc_vars) and listEmpty(disc_eqns) then - // if there are no discrete variables > don't do anything - residual_lst := strict.residual_eqns; - inner_comps := {}; - else - // fail and report if length is not equal! - if listLength(disc_vars) <> listLength(disc_eqns) then - Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() - + " failed because there is an unequal amount of discrete variables and equations:\n" - + List.toString(disc_vars, BVariable.pointerToString, "discrete variables", "\t", "\n\t", "\n\n") - + List.toString(disc_eqns, function Equation.pointerToString(str = ""), "discrete equations", "\t", "\n\t", "\n\n")}); - fail(); - end if; - // solve the equations for linear occurences of the discrete variables - // it should be: - // solve discrete vars <-> discrete eqs - discreteVars := VariablePointers.fromList(disc_vars); - eqns := EquationPointers.fromList(disc_eqns); - - (matching, inner_comps) := Causalize.simple(discreteVars, eqns, NBAdjacency.MatrixStrictness.LINEAR); - (_, _, _, residual_lst) := Matching.getMatches(matching, NONE(), discreteVars, eqns); + if listLength(disc_vars) <> listLength(disc_eqns) then + Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + + " failed because there is an unequal amount of discrete variables and equations:\n" + + List.toString(disc_vars, BVariable.pointerToString, "discrete variables", "\t", "\n\t", "\n\n") + + List.toString(disc_eqns, function Equation.pointerToString(str = ""), "discrete equations", "\t", "\n\t", "\n\n")}); + fail(); end if; - // create residual equations - strict.residual_eqns := list(Slice.apply(eqn, function Equation.createResidual(new = true)) for eqn in strict.residual_eqns); - residual_comps := list(StrongComponent.fromSolvedEquationSlice(eqn) for eqn in strict.residual_eqns); - - // update jacobian to take slices (just to have correct inner variables and such) - (jacobian, funcTree) := BJacobian.nonlinear( - variables = VariablePointers.fromList(list(Slice.getT(var) for var in comp.strict.iteration_vars)), - equations = EquationPointers.fromList(list(Slice.getT(eqn) for eqn in comp.strict.residual_eqns)), - comps = listArray(residual_comps), - funcTree = funcTree, - name = System.System.systemTypeString(systemType) + "_NLS_JAC_" + intString(index)); - - strict.iteration_vars := iter_lst; - strict.residual_eqns := residual_lst; - strict.innerEquations := listArray(inner_comps); - strict.jac := jacobian; - comp.strict := strict; - */ - if Flags.isSet(Flags.TEARING_DUMP) then - print(StrongComponent.toString(comp) + "\n"); + if not listEmpty(disc_vars) then + comp.mixed := true; + + // the sets of discrete variables and discrete equations + v := UnorderedMap.subMap(variables.map, list(BVariable.getVarName(var) for var in disc_vars)); + e := UnorderedMap.subMap(equations.map, list(Equation.getEqnName(eqn) for eqn in disc_eqns)); + + // match the discretes to create inner components + adj := Adjacency.Matrix.fromFull(full, v, e, equations, NBAdjacency.MatrixStrictness.MATCHING); + matching := Matching.regular(NBMatching.EMPTY_MATCHING, adj, true, true); + adj := Adjacency.Matrix.upgrade(adj, full, v, e, equations, NBAdjacency.MatrixStrictness.SORTING); + inner_comps := Sorting.tarjan(adj, matching, variables, equations); //probably need other variables and equations here? + strict.innerEquations := listArray(inner_comps); + + // create residuals equations and iteration variables + (_, _, _, residual_lst) := Matching.getMatches(matching, Adjacency.Matrix.getMappingOpt(full), variables, equations); + strict.residual_eqns := residual_lst; + strict.iteration_vars := list(Slice.SLICE(var, {}) for var in cont_vars); + comp.strict := strict; end if; - then (comp, index); - else (comp, index); + then comp; + else comp; end match; - end tearingMinimal; + end minimal; -/* - function tearingMinimal extends Module.tearingInterface; + function checkLinearity + input Adjacency.Matrix full; + input UnorderedMap v; + input UnorderedMap e; + output Boolean linear = true; protected - list residual_comps; - Option jacobian; - StrongComponent new_comp; - list> residuals = {}; - list> cont_lst, disc_lst; + list var_lst = UnorderedMap.keyList(v); algorithm - (comp, index) := match comp - case StrongComponent.ALGEBRAIC_LOOP() algorithm - index := index + 1; - comp.idx := index; - - //(cont_lst, disc_lst) := List.splitOnTrue(comp.strict.residual_vars, BVariable.isContinuous); - - // create residual equations - for eqn in listReverse(comp.strict.iteration_vars) loop - residuals := Slice.apply(eqn, function Equation.createResidual(new = true)) :: residuals; + linear := match full + case Adjacency.Matrix.FULL() algorithm + for eqn_idx in UnorderedMap.valueList(e) loop + for var in var_lst loop + if UnorderedSet.contains(var, full.occurences[eqn_idx]) and Solvability.nonlinearOrImplicit(UnorderedMap.getSafe(var, full.solvabilities[eqn_idx], sourceInfo())) then + linear := false; break; + end if; + end for; end for; - comp.strict := setResidualEqns(comp.strict, residuals); - residual_comps := list(StrongComponent.fromSolvedEquationSlice(eqn) for eqn in comp.strict.residual_eqns); - - // update jacobian to take slices (just to have correct inner variables and such) - (jacobian, funcTree) := BJacobian.nonlinear( - variables = VariablePointers.fromList(list(Slice.getT(var) for var in comp.strict.iteration_vars)), - equations = EquationPointers.fromList(list(Slice.getT(eqn) for eqn in comp.strict.residual_eqns)), - comps = listArray(residual_comps), - funcTree = funcTree, - name = System.System.systemTypeString(systemType) + "_NLS_JAC_" + intString(index)); - - new_comp := StrongComponent.addLoopJacobian(comp, jacobian); - if Flags.isSet(Flags.TEARING_DUMP) then - print(StrongComponent.toString(comp) + "\n"); - end if; - then (new_comp, index); - else (comp, index); + then true; + else algorithm + Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " expected type full, got type " + Adjacency.Matrix.strictnessString(Adjacency.Matrix.getStrictness(full)) + "."}); + then fail(); end match; - end tearingMinimal; - - function tearingMinimalWork - input String name; - input list> variables; - input list> equations; - input Boolean mixed; - output StrongComponent comp; - input output FunctionTree funcTree; - input output Integer index; - protected - list> cont_lst, disc_lst; - list> iteration_lst; - list> residual_lst; - VariablePointers discreteVars; - EquationPointers eqns; - Adjacency.Matrix adj; - Matching matching; - list inner_comps, residual_comps; - Tearing tearingSet; - Option jacobian; + end checkLinearity; - algorithm - index := index + 1; - - (cont_lst, disc_lst) := List.splitOnTrue(variables, BVariable.isContinuous); - iteration_lst := list(Slice.SLICE(var, {}) for var in cont_lst); - if listEmpty(disc_lst) then - residual_lst := list(Slice.SLICE(eqn, {}) for eqn in equations); - inner_comps := {}; - else - discreteVars := VariablePointers.fromList(disc_lst); - eqns := EquationPointers.fromList(equations); - - (adj, SOME(funcTree)) := Adjacency.Matrix.create(discreteVars, eqns, NBAdjacency.MatrixStrictness.LINEAR, SOME(funcTree)); - matching := Matching.regular(NBMatching.EMPTY_MATCHING, adj, true, false); - (_, _, _, residual_lst) := Matching.getMatches(matching, NONE(), discreteVars, eqns); - inner_comps := Sorting.tarjan(adj, matching, discreteVars, eqns); - end if; - - // create residual equations - residual_lst := list(Slice.apply(eqn, function Equation.createResidual(new = true)) for eqn in residual_lst); - residual_comps := list(StrongComponent.fromSolvedEquationSlice(eqn) for eqn in residual_lst); - - tearingSet := TEARING_SET(iteration_lst, residual_lst, listArray(inner_comps), NONE()); - comp := StrongComponent.ALGEBRAIC_LOOP(index, tearingSet, NONE(), false, mixed, NBSolve.Status.IMPLICIT); - - // inner equations are part of the jacobian - (jacobian, funcTree) := BJacobian.nonlinear( - variables = VariablePointers.fromList(cont_lst), - equations = EquationPointers.fromList(list(Slice.getT(res) for res in residual_lst)), - comps = listArray(listAppend(inner_comps, residual_comps)), - funcTree = funcTree, - name = name + intString(index) - ); - - comp := StrongComponent.addLoopJacobian(comp, jacobian); - if Flags.isSet(Flags.TEARING_DUMP) and not listEmpty(disc_lst) then - print(StrongComponent.toString(comp) + "\n"); - end if; - end tearingMinimalWork; -*/ annotation(__OpenModelica_Interface="backend"); end NBTearing; diff --git a/OMCompiler/Compiler/NBackEnd/Util/NBAdjacency.mo b/OMCompiler/Compiler/NBackEnd/Util/NBAdjacency.mo index cf446f9146..6d74af4fbf 100644 --- a/OMCompiler/Compiler/NBackEnd/Util/NBAdjacency.mo +++ b/OMCompiler/Compiler/NBackEnd/Util/NBAdjacency.mo @@ -55,6 +55,7 @@ protected import NBDifferentiate.{DifferentiationArguments, DifferentiationType}; import BEquation = NBEquation; import NBEquation.{Equation, EquationAttributes, EquationPointers, Iterator, IfEquationBody, WhenEquationBody, WhenStatement}; + import Solve = NBSolve; import BVariable = NBVariable; import NBVariable.VariablePointers; @@ -701,43 +702,79 @@ public function refine "refines the solvability kind using differentiation - Note: only updates the solvabilites of the variables and equations from the maps" + Note: only updates the solvabilites of the variables and equations from the maps v and e" input output Matrix full; input output FunctionTree funcTree; input UnorderedMap v "variables to refine"; input UnorderedMap e "equations to refine"; input VariablePointers vars "all variables"; input EquationPointers eqns "all equations"; + input UnorderedSet vars_set "context variables to determine solvability"; algorithm full := match full local ComponentRef eqn, var; Integer eqn_idx, var_idx; - DifferentiationArguments diffArgs; + DifferentiationArguments diffArgs = DifferentiationArguments.default(NBDifferentiate.DifferentiationType.SIMPLE, funcTree); Pointer eqn_ptr; Expression exp; + Solve.Status status; + Solvability sol; + UnorderedSet linear_set; + list param_lst, var_lst; case FULL() algorithm for v_tpl in UnorderedMap.toList(v) loop (var, var_idx) := v_tpl; - diffArgs := DifferentiationArguments.default(NBDifferentiate.DifferentiationType.SIMPLE, funcTree); diffArgs.diffCref := var; for e_tpl in UnorderedMap.toList(e) loop (eqn, eqn_idx) := e_tpl; - eqn_ptr := EquationPointers.getEqnAt(eqns, eqn_idx); - // get the residual expression, differentiate and simplify it - exp := Equation.getResidualExp(Pointer.access(eqn_ptr)); - (exp, diffArgs) := Differentiate.differentiateExpressionDump(exp, diffArgs, getInstanceName()); - exp := SimplifyExp.simplifyDump(exp, true, getInstanceName()); - print("the partial derivative by " + ComponentRef.toString(var) + " of equation\n" + Equation.pointerToString(eqn_ptr) + "\nis: " + Expression.toString(exp) + "\n"); + // only do something if there is a value to refine + if UnorderedSet.contains(var, full.occurences[eqn_idx]) then + // only do something if it is not implicit or unsolvable) + sol := UnorderedMap.getSafe(var, full.solvabilities[eqn_idx], sourceInfo()); + if Solvability.rank(sol) < Solvability.rank(Solvability.IMPLICIT()) then + eqn_ptr := EquationPointers.getEqnAt(eqns, eqn_idx); + // booleans or (todo: enumerations) + if Type.isBoolean(Equation.getType(Pointer.access(eqn_ptr))) or Type.isBoolean(ComponentRef.getSubscriptedType(var)) then + // if the equation or cref type is boolean, it can only be solved if its isolated in the LHS or RHS + // Use solveSimple for this and check if status is EXPLICIT + (_, status, _) := Solve.solveSimple(Pointer.access(eqn_ptr), var); + sol := if status == NBSolve.Status.EXPLICIT then Solvability.EXPLICIT_LINEAR(NONE(), NONE()) else Solvability.UNSOLVABLE(); + else + // get the residual expression, differentiate and simplify it + exp := Equation.getResidualExp(Pointer.access(eqn_ptr)); + (exp, diffArgs) := Differentiate.differentiateExpressionDump(exp, diffArgs, getInstanceName()); + exp := SimplifyExp.simplifyDump(exp, true, getInstanceName()); + if Expression.isZero(exp) then + sol := Solvability.UNSOLVABLE(); + elseif Expression.containsCrefSet(exp, vars_set) then + // nonlinear -> unique solution if does not contain the variable itself + sol := Solvability.EXPLICIT_NONLINEAR(Expression.containsCref(exp, var)); + else + // linear -> find all contained crefs and split them by kind. remove constants and save params / variables + linear_set := Expression.extractCrefs(exp); + (_, var_lst) := List.splitOnTrue(UnorderedSet.toList(linear_set), function BVariable.checkCref(func = BVariable.isConst)); + (param_lst, var_lst) := List.splitOnTrue(UnorderedSet.toList(linear_set), function BVariable.checkCref(func = BVariable.isParamOrConst)); + sol := Solvability.EXPLICIT_LINEAR( + pars = if listEmpty(param_lst) then NONE() else SOME(UnorderedSet.fromList(param_lst, ComponentRef.hash, ComponentRef.isEqual)), + vars = if listEmpty(var_lst) then NONE() else SOME(UnorderedSet.fromList(var_lst, ComponentRef.hash, ComponentRef.isEqual))); + end if; + end if; + UnorderedMap.add(var, sol, full.solvabilities[eqn_idx]); + end if; + end if; end for; end for; then full; - else algorithm Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " expected type full, got type " + strictnessString(getStrictness(full)) + "."}); then fail(); end match; + + if Flags.isSet(Flags.BLT_MATRIX_DUMP) then + print(toString(full, "Refined Full") + "\n"); + end if; end refine; function compress @@ -998,8 +1035,9 @@ public output Option mapping; algorithm mapping := match adj + case FULL() then SOME(adj.mapping); case FINAL() then SOME(adj.mapping); - else NONE(); + else NONE(); end match; end getMappingOpt; @@ -1380,7 +1418,7 @@ public end EXPLICIT_NONLINEAR; record EXPLICIT_LINEAR - Boolean param "true if we need to divide by a parameter to solve"; + Option> pars "parameters we need to divide by to solve"; Option> vars "variables we need to divide by to solve"; end EXPLICIT_LINEAR; @@ -1392,7 +1430,7 @@ public case UNSOLVABLE() then "XX"; case IMPLICIT() then "II"; case EXPLICIT_NONLINEAR() then "N" + (if sol.unique then "+" else "-"); - case EXPLICIT_LINEAR() then "L" + (if Util.isSome(sol.vars) then "V" elseif sol.param then "P" else "C"); + case EXPLICIT_LINEAR() then "L" + (if Util.isSome(sol.vars) then "V" elseif Util.isSome(sol.pars) then "P" else "C"); case UNKNOWN() then "||"; else algorithm Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because of unknown solvability kind."}); @@ -1410,7 +1448,7 @@ public case EXPLICIT_NONLINEAR(unique = false) then 5; case EXPLICIT_NONLINEAR() then 4; case EXPLICIT_LINEAR(vars = SOME(_)) then 3; - case EXPLICIT_LINEAR(param = true) then 2; + case EXPLICIT_LINEAR(pars = SOME(_)) then 2; case EXPLICIT_LINEAR() then 1; case UNKNOWN() then 0; else algorithm @@ -1460,7 +1498,7 @@ public case EXPLICIT_NONLINEAR(unique = false) algorithm NM := cref :: NM; then(); case EXPLICIT_NONLINEAR() algorithm NP := cref :: NP; then(); case EXPLICIT_LINEAR(vars = SOME(_)) algorithm LV := cref :: LV; then(); - case EXPLICIT_LINEAR(param = true) algorithm LP := cref :: LP; then(); + case EXPLICIT_LINEAR(pars = SOME(_)) algorithm LP := cref :: LP; then(); case EXPLICIT_LINEAR() algorithm LC := cref :: LC; then(); else algorithm QQ := cref :: QQ; then(); end match; @@ -1494,12 +1532,23 @@ public output Solvability sol; algorithm sol := match st - case MatrixStrictness.LINEAR then EXPLICIT_LINEAR(false, SOME(UnorderedSet.new(ComponentRef.hash, ComponentRef.isEqual))); + case MatrixStrictness.LINEAR then EXPLICIT_LINEAR(NONE(), SOME(UnorderedSet.new(ComponentRef.hash, ComponentRef.isEqual))); case MatrixStrictness.MATCHING then IMPLICIT(); case MatrixStrictness.SORTING then UNSOLVABLE(); else UNKNOWN(); end match; end fromStrictness; + + function nonlinearOrImplicit + input Solvability sol; + output Boolean b; + algorithm + b := match sol + case EXPLICIT_NONLINEAR() then true; + case IMPLICIT() then true; + else false; + end match; + end nonlinearOrImplicit; end Solvability; function collectDependenciesEquation @@ -1539,7 +1588,7 @@ public Dependency.addListFull(eqn.alg.outputs, dep_map, rep_set); // make inputs unsolvable and outputs solvable Solvability.updateList(eqn.alg.inputs, Solvability.UNSOLVABLE(), sol_map); - Solvability.updateList(eqn.alg.outputs, Solvability.EXPLICIT_LINEAR(false, NONE()), sol_map); + Solvability.updateList(eqn.alg.outputs, Solvability.EXPLICIT_LINEAR(NONE(), NONE()), sol_map); then UnorderedSet.fromList(listAppend(eqn.alg.inputs, eqn.alg.outputs), ComponentRef.hash, ComponentRef.isEqual); case Equation.FOR_EQUATION(body = {body}) algorithm @@ -1776,7 +1825,7 @@ public if not UnorderedMap.contains(cref, dep_map) then UnorderedMap.add(cref, Dependency.create(ComponentRef.getSubscriptedType(cref)), dep_map); end if; - Solvability.update(cref, Solvability.EXPLICIT_LINEAR(false, NONE()), sol_map); + Solvability.update(cref, Solvability.EXPLICIT_LINEAR(NONE(), NONE()), sol_map); crefs := {cref}; else var := BVariable.getVarPointer(cref); diff --git a/OMCompiler/Compiler/NFFrontEnd/NFExpression.mo b/OMCompiler/Compiler/NFFrontEnd/NFExpression.mo index e97fe7a31d..15efedf2eb 100644 --- a/OMCompiler/Compiler/NFFrontEnd/NFExpression.mo +++ b/OMCompiler/Compiler/NFFrontEnd/NFExpression.mo @@ -41,14 +41,15 @@ protected import Builtin = NFBuiltin; import BuiltinCall = NFBuiltinCall; + import Ceval = NFCeval; + import ComplexType = NFComplexType; + import ExpandExp = NFExpandExp; import Expression = NFExpression; import Function = NFFunction; import NFPrefixes.{Variability, Purity}; import Prefixes = NFPrefixes; - import Ceval = NFCeval; - import ComplexType = NFComplexType; - import ExpandExp = NFExpandExp; import TypeCheck = NFTypeCheck; + import UnorderedSet; import ValuesUtil; import MetaModelica.Dangerous.*; import RangeIterator = NFRangeIterator; @@ -4237,32 +4238,22 @@ public CREF(cref = cref) := exp; end toCref; - function extract - "author: kabdelhak 2020-06 - Extracts all sub expressions from an expression using a filter function." + function extractCrefs input Expression exp; - input filter func; - output list exp_lst; - partial function filter - input Expression exp; - output Boolean b; - end filter; - protected - // traverse helper function only needed in this function - function traverser - input Expression exp; - input filter func; - input output list exp_lst; - partial function filter - input Expression exp; - output Boolean b; - end filter; - algorithm - exp_lst := if func(exp) then exp :: exp_lst else exp_lst; - end traverser; + output UnorderedSet crefs = fold(exp, extractCref, UnorderedSet.new(ComponentRef.hash, ComponentRef.isEqual)); + end extractCrefs; + + function extractCref + input Expression exp; + input output UnorderedSet crefs; algorithm - exp_lst := fold(exp, function traverser(func = func), {}); - end extract; + crefs := match exp + case CREF() algorithm + UnorderedSet.add(exp.cref, crefs); + then crefs; + else crefs; + end match; + end extractCref; function isIterator input Expression exp; @@ -5578,6 +5569,7 @@ public end isPure; function containsCref + "returns true if the expression contains the cref" input Expression exp; input ComponentRef cref; output Boolean b; @@ -5596,6 +5588,26 @@ public end match; end isCrefEqual; + function containsCrefSet + "returns true if the expression contains any crefs in the set" + input Expression exp; + input UnorderedSet set; + output Boolean b; + algorithm + b := fold(exp, function isCrefEqualSet(set = set), false); + end containsCrefSet; + + function isCrefEqualSet + input Expression exp; + input output Boolean b; + input UnorderedSet set; + algorithm + b := match (b, exp) + case (false, CREF()) then UnorderedSet.contains(exp.cref, set); + else b; + end match; + end isCrefEqualSet; + function filterSplitIndices input output Expression exp; input InstNode node; diff --git a/OMCompiler/Compiler/Util/UnorderedMap.mo b/OMCompiler/Compiler/Util/UnorderedMap.mo index ab875aacd3..17123f2ddc 100644 --- a/OMCompiler/Compiler/Util/UnorderedMap.mo +++ b/OMCompiler/Compiler/Util/UnorderedMap.mo @@ -556,7 +556,7 @@ public end for; end merge; - function subSet + function subMap input UnorderedMap map; input list lst; output UnorderedMap sub_set; @@ -571,7 +571,7 @@ public for k in lst loop add(k, getSafe(k, map, sourceInfo()), sub_set); end for; - end subSet; + end subMap; function all "Returns true if the given function returns true for all values in the map,