Skip to content

Commit

Permalink
Replace shared crefs with zero in resolve loops (#12167)
Browse files Browse the repository at this point in the history 
Fixes #11236

---------

Co-authored-by: kabdelhak <38032125+kabdelhak@users.noreply.github.com>
  • Loading branch information
@phannebohm @kabdelhak
phannebohm and kabdelhak committed Mar 29, 2024
1 parent ec34d31 commit 6c2fe7e
Show file tree
Hide file tree
Showing 6 changed files with 89 additions and 65 deletions.
144 changes: 84 additions & 60 deletions OMCompiler/Compiler/BackEnd/ResolveLoops.mo
View file
Expand Up @@ -45,6 +45,7 @@ import AvlSetInt;
import BackendDAEUtil;
import BackendEquation;
import BackendVariable;
import BackendVarTransform;
import BackendDump;
import ComponentReference;
import Expression;
Expand All @@ -54,8 +55,8 @@ import ExpressionDump;
import Flags;
import HpcOmTaskGraph;
import List;
import Util;
import Tearing;
import Util;

public function resolveLoops "author:Waurich TUD 2013-12
traverses the equations and finds simple equations(i.e. linear functions
Expand Down Expand Up @@ -855,7 +856,7 @@ algorithm
(daeEqsOut,replEqsOut) := matchcontinue(loopsIn,eqCrossLstIn,varCrossLstIn,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn,replEqsIn)
local
Integer pos,crossEq,crossVar,eq1,eq2;
list<Integer> loop1, eqs, vars, crossEqs, crossEqs2, removeCrossEqs, crossVars, replEqs, loopVars, adjVars;
list<Integer> loop1, eqs, vars, crossEqs, crossEqs2, removeCrossEqs, crossVars, replEqs, loopVars, adjVars, m_row;
list<list<Integer>> rest, eqVars;
BackendDAE.Equation resolvedEq;
BackendDAE.EquationArray daeEqs;
Expand All @@ -868,7 +869,7 @@ algorithm
// only eqCrossNodes
//print("only eqCrossNodes\n");
loop1 = List.unique(loop1);
resolvedEq = resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
(resolvedEq, m_row) = resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);

// get the equation that will be replaced and the rest
(crossEqs,eqs,_) = List.intersection1OnTrue(loop1,eqCrossLstIn,intEq); // replace a crossEq in the loop
Expand Down Expand Up @@ -909,6 +910,7 @@ algorithm
// replace Equation
//print("replace equation "+intString(pos)+"\n");
replEqs = pos::replEqsIn;
arrayUpdate(mIn,pos,m_row);
pos = arrayGet(eqMap,pos);
daeEqs = BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);

Expand All @@ -920,7 +922,7 @@ algorithm
// only varCrossNodes
//print("only varCrossNodes\n");
loop1 = List.unique(loop1);
resolvedEq = resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
(resolvedEq, m_row) = resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);

// get the equation that will be replaced and the rest
(replEqs,_,eqs) = List.intersection1OnTrue(replEqsIn,loop1,intEq); // just consider the already replaced equations in this loop
Expand Down Expand Up @@ -962,6 +964,7 @@ algorithm
// replace Equation
//print("replace equation "+intString(pos)+"\n");
replEqs = pos::replEqsIn;
arrayUpdate(mIn,pos,m_row);
pos = arrayGet(eqMap,pos);
daeEqs = BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);

Expand All @@ -973,7 +976,7 @@ algorithm
// single Loop
loop1 = List.unique(loop1);
//print("single loop\n");
resolvedEq = resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
(resolvedEq, m_row) = resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);

// update AdjacencyMatrix
(_,crossEqs,_) = List.intersection1OnTrue(loop1,replEqsIn,intEq); // do not replace an already replaced Eq
Expand All @@ -987,8 +990,10 @@ algorithm
// replace Equation
//print("replace equation "+intString(pos)+"\n");
replEqs = pos::replEqsIn;
arrayUpdate(mIn,pos,m_row);
pos = arrayGet(eqMap,pos);
daeEqs = BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);

(daeEqs,replEqs) = resolveLoops_resolveAndReplace(rest,eqCrossLstIn,varCrossLstIn,mIn,mTIn,eqMap,varMap,daeEqs,daeVarsIn,replEqs);
then
(daeEqs,replEqs);
Expand All @@ -998,7 +1003,7 @@ algorithm
//print("both eqCrossNodes and varCrossNodes, loopLength"+intString(listLength(loop1))+"\n");
loop1 := List.unique(loop1);
true := listLength(loop1) == 2;
resolvedEq := resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
(resolvedEq, m_row) := resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
//print("resolved eq to "+BackendDump.equationString(resolvedEq)+"\n");

if eqIsConst(resolvedEq)then
Expand All @@ -1016,6 +1021,7 @@ algorithm
end if;
replEqs := pos::replEqsIn;
//print("contract eqs: "+stringDelimitList(List.map(loop1,intString),",")+" to eq "+intString(pos)+"\n");
arrayUpdate(mIn,pos,m_row);
pos := arrayGet(eqMap,pos);
daeEqs := BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);
else
Expand All @@ -1025,10 +1031,6 @@ algorithm
(daeEqs,replEqs) := resolveLoops_resolveAndReplace(rest,eqCrossLstIn,varCrossLstIn,mIn,mTIn,eqMap,varMap,daeEqs,daeVarsIn,replEqs);
then
(daeEqs,replEqs);
else
equation
print("resolveLoops_resolveAndReplace failed!\n");
then fail();
end matchcontinue;
end resolveLoops_resolveAndReplace;

Expand Down Expand Up @@ -1218,6 +1220,7 @@ protected function resolveClosedLoop "author:Waurich TUD 2014-02
input BackendDAE.EquationArray daeEqsIn;
input BackendDAE.Variables daeVarsIn;
output BackendDAE.Equation eqOut;
output list<Integer> m_row;
protected
Integer startEqIdx,startEqDaeIdx;
list<Integer> loop1, restLoop;
Expand All @@ -1226,69 +1229,90 @@ algorithm
startEqIdx::restLoop := loopIn;
startEqDaeIdx := arrayGet(eqMap,startEqIdx);
loop1 := sortLoop(restLoop,m,mT,{startEqIdx});
//print("solve the loop: "+stringDelimitList(List.map(loop1,intString),",")+"\n");
if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) and listLength(loop1) > 1 then
print("solve the loop: " + List.toString(loop1, intString) + "\n");
end if;
eq := BackendEquation.get(daeEqsIn,startEqDaeIdx);
eqOut := resolveClosedLoop2(eq,loop1,m,mT,eqMap,varMap,daeEqsIn,daeVarsIn);
(eqOut, m_row) := resolveClosedLoop2(eq,loop1,m, arrayGet(m,startEqIdx), eqMap,varMap,daeEqsIn,daeVarsIn);
end resolveClosedLoop;

protected function resolveClosedLoop2 "author:Waurich TUD 2013-12"
input BackendDAE.Equation eqIn;
input output BackendDAE.Equation eq;
input list<Integer> loopIn;
input BackendDAE.AdjacencyMatrix m;
input BackendDAE.AdjacencyMatrixT mT;
input output list<Integer> m_row "row of eq in m";
input array<Integer> eqMap;
input array<Integer> varMap;
input BackendDAE.EquationArray daeEqsIn;
input BackendDAE.Variables daeVarsIn;
output BackendDAE.Equation eqOut;
algorithm
(eqOut) := matchcontinue(eqIn,loopIn,m,mT,eqMap,varMap,daeEqsIn,daeVarsIn)
(eq, m_row) := match loopIn
local
Boolean isPosOnRhs1, isPosOnRhs2, algSign;
Integer eqIdx1, eqIdx2, eqDaeIdx2, varIdx, daeVarIdx;
list<Integer> adjVars, adjVars1 ,adjVars2, eqIdcs,varIdcs,varNodes,restLoop, row;
list<BackendDAE.Equation> eqs;
BackendDAE.Equation eq2,resolvedEq;
BackendDAE.Var var;
DAE.ComponentRef cref, eqCrefs;
case(_,{_},_,_,_,_,_,_)
equation
//print("finished loop\n");
then
eqIn;
case(_,(eqIdx1::restLoop),_,_,_,_,_,_)
equation
// the equation to add
eqIdx2 = listHead(restLoop);
eqDaeIdx2 = arrayGet(eqMap,eqIdx2);
eq2 = BackendEquation.get(daeEqsIn,eqDaeIdx2);

// get the vars that are shared of the 2 equations
adjVars1 = arrayGet(m,eqIdx1);
adjVars2 = arrayGet(m,eqIdx2);
(adjVars,adjVars1,_) = List.intersection1OnTrue(adjVars1,adjVars2,intEq);

// just take the first
varIdx = listHead(adjVars);
daeVarIdx = arrayGet(varMap,varIdx);
var = BackendVariable.getVarAt(daeVarsIn,daeVarIdx);
cref = BackendVariable.varCref(var);

// check the algebraic signs
isPosOnRhs1 = CRefIsPosOnRHS(cref,eqIn);
isPosOnRhs2 = CRefIsPosOnRHS(cref,eq2);
algSign = boolOr((not isPosOnRhs1) and isPosOnRhs2,(not isPosOnRhs2) and isPosOnRhs1); // XOR
resolvedEq = sumUp2Equations(algSign,eqIn,eq2);
if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) then
print("From eqs \n"+BackendDump.equationString(eqIn)+"\n"+BackendDump.equationString(eq2)+"\n");
print("resolved the eq \n"+BackendDump.equationString(resolvedEq)+"\n\n");
end if;
resolvedEq = resolveClosedLoop2(resolvedEq,restLoop,m,mT,eqMap,varMap,daeEqsIn,daeVarsIn);
then
resolvedEq;
end matchcontinue;
Boolean algSign;
Integer eqIdx1, eqIdx2;
list<Integer> adjVars, adjVars1, adjVars2, restLoop, posVars, negVars;
list<DAE.ComponentRef> adjCrefs;
BackendDAE.Equation eq2, eq3, resolvedEq;
BackendVarTransform.VariableReplacements replacements;
case {_} then (eq, m_row);
case eqIdx1::eqIdx2::restLoop algorithm
// the equation to add
eq2 := BackendEquation.get(daeEqsIn, arrayGet(eqMap,eqIdx2));

// get the vars that are shared of the 2 equations
adjVars1 := m_row;
adjVars2 := arrayGet(m,eqIdx2);
(adjVars, adjVars1, adjVars2) := List.intersection1OnTrue(adjVars1, adjVars2, intEq);

// split shared vars by sign
(posVars, negVars) := List.splitOnTrue(adjVars, function varSign(varMap = varMap, daeVarsIn = daeVarsIn, eq1 = eq, eq2 = eq2));
algSign := listLength(posVars) > listLength(negVars); // choose set with more canceling vars
adjCrefs := list(crefFromIndex(idx, varMap, daeVarsIn) for idx in (if algSign then posVars else negVars));

// shared crefs are removed from adjacecy
m_row := List.flatten({adjVars1, adjVars2, (if algSign then negVars else posVars)});

// replace `cref` with zero to make the job easier for `simplify`
replacements := BackendVarTransform.emptyReplacementsSized(listLength(adjCrefs));
replacements := BackendVarTransform.addReplacements(replacements, adjCrefs, list(Expression.createZeroExpression(ComponentReference.crefTypeFull(c)) for c in adjCrefs), NONE());
({resolvedEq, eq3}, _) := BackendVarTransform.replaceEquations({eq, eq2}, replacements, NONE());

resolvedEq := sumUp2Equations(algSign,resolvedEq,eq3);
if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) then
print("From eqs \n"+BackendDump.equationString(eq)+"\n"+BackendDump.equationString(eq2)+"\n");
print("resolved the eq \n"+BackendDump.equationString(resolvedEq)+"\n\n");
end if;
then resolveClosedLoop2(resolvedEq,eqIdx2::restLoop,m, m_row, eqMap,varMap,daeEqsIn,daeVarsIn);
end match;
end resolveClosedLoop2;

protected function crefFromIndex
input Integer varIdx;
input array<Integer> varMap;
input BackendDAE.Variables daeVarsIn;
output DAE.ComponentRef cref;
protected
Integer daeVarIdx;
BackendDAE.Var var;
algorithm
daeVarIdx := arrayGet(varMap, varIdx);
var := BackendVariable.getVarAt(daeVarsIn, daeVarIdx);
cref := BackendVariable.varCref(var);
end crefFromIndex;

protected function varSign
input Integer index;
input array<Integer> varMap;
input BackendDAE.Variables daeVarsIn;
input BackendDAE.Equation eq1;
input BackendDAE.Equation eq2;
output Boolean algSign;
protected
DAE.ComponentRef cref = crefFromIndex(index, varMap, daeVarsIn);
algorithm
algSign := CRefIsPosOnRHS(cref, eq1) <> CRefIsPosOnRHS(cref, eq2) "check the algebraic signs"; // XOR
end varSign;

public function sortLoop "author:Waurich TUD 2014-01
sorts the equations in a loop so that they are solved in a row."
input list<Integer> loopIn;
Expand Down Expand Up @@ -1674,7 +1698,7 @@ algorithm
exp2 := sumUp2Expressions(sumUp, exp2, exp4);
exp2 := sumUp2Expressions(false, exp2, exp1);
(exp2, _) := ExpressionSimplify.simplify(exp2);
exp1 := DAE.RCONST(0.0);
exp1 := Expression.createZeroExpression(Expression.typeof(exp2));
eqOut := BackendDAE.EQUATION(exp1, exp2, DAE.emptyElementSource, BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN);
eqOut := simplifyZeroAssignment(eqOut);
end sumUp2Equations;
Expand Down
2 changes: 1 addition & 1 deletion ...es/msl32/Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMC_DOL_MultiPhase.mos
View file
Expand Up @@ -38,7 +38,7 @@ runScript(modelTesting);getErrorString();
// OpenModelicaModelTesting.Kind.VerifiedSimulation
// Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMC_DOL_MultiPhase
// {"aimc3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[1].Phi.re","aimc3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[2].Phi.im","aimc3.stator.strayReluctance.port_p.Phi.im","aimc3.stator.strayReluctance.port_p.Phi.re","aimcM.rotorCage.strayInductor.inductor[4].i","aimcM.rotorCage.strayInductor.inductor[5].i","aimcM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[1].Phi.re","aimcM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[2].Phi.im","aimcM.stator.strayReluctance.port_p.Phi.im","aimcM.stator.strayReluctance.port_p.Phi.re","loadInertia3.phi","loadInertia3.w","loadInertiaM.phi","loadInertiaM.w"}
// Simulation options: startTime = 0.0, stopTime = 1.5, numberOfIntervals = 7500, tolerance = 1e-05, method = 'dassl', fileNamePrefix = 'Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMC_DOL_MultiPhase', options = '', outputFormat = 'mat', variableFilter = 'time|aimc3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimc3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimc3.stator.strayReluctance.port_p.Phi.im|aimc3.stator.strayReluctance.port_p.Phi.re|aimcM.rotorCage.strayInductor.inductor.4..i|aimcM.rotorCage.strayInductor.inductor.5..i|aimcM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimcM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimcM.stator.strayReluctance.port_p.Phi.im|aimcM.stator.strayReluctance.port_p.Phi.re|loadInertia3.phi|loadInertia3.w|loadInertiaM.phi|loadInertiaM.w', cflags = '', simflags = ' -abortSlowSimulation -alarm=360 -emit_protected'
// Simulation options: startTime = 0.0, stopTime = 1.5, numberOfIntervals = 7500, tolerance = 1e-5, method = 'dassl', fileNamePrefix = 'Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMC_DOL_MultiPhase', options = '', outputFormat = 'mat', variableFilter = 'time|aimc3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimc3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimc3.stator.strayReluctance.port_p.Phi.im|aimc3.stator.strayReluctance.port_p.Phi.re|aimcM.rotorCage.strayInductor.inductor.4..i|aimcM.rotorCage.strayInductor.inductor.5..i|aimcM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimcM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimcM.stator.strayReluctance.port_p.Phi.im|aimcM.stator.strayReluctance.port_p.Phi.re|loadInertia3.phi|loadInertia3.w|loadInertiaM.phi|loadInertiaM.w', cflags = '', simflags = ' -abortSlowSimulation -alarm=360 -emit_protected'
// Result file: Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMC_DOL_MultiPhase_res.mat
// Messages: LOG_SUCCESS | info | The initialization finished successfully without homotopy method.
// LOG_SUCCESS | info | The simulation finished successfully.
Expand Down
2 changes: 1 addition & 1 deletion .../msl32/Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMS_Start_MultiPhase.mos
View file
Expand Up @@ -38,7 +38,7 @@ runScript(modelTesting);getErrorString();
// OpenModelicaModelTesting.Kind.VerifiedSimulation
// Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMS_Start_MultiPhase
// {"aims3.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter[1].Phi.re","aims3.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter[2].Phi.im","aims3.rotor.zeroInductor.i0","aims3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[1].Phi.re","aims3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[2].Phi.im","aimsM.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter[1].Phi.re","aimsM.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter[2].Phi.im","aimsM.rotor.zeroInductor.i0","aimsM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[1].Phi.re","aimsM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter[2].Phi.im","loadInertia3.phi","loadInertia3.w","loadInertiaM.phi","loadInertiaM.w"}
// Simulation options: startTime = 0.0, stopTime = 1.5, numberOfIntervals = 1500, tolerance = 1e-05, method = 'dassl', fileNamePrefix = 'Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMS_Start_MultiPhase', options = '', outputFormat = 'mat', variableFilter = 'time|aims3.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aims3.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aims3.rotor.zeroInductor.i0|aims3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aims3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimsM.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimsM.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimsM.rotor.zeroInductor.i0|aimsM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimsM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|loadInertia3.phi|loadInertia3.w|loadInertiaM.phi|loadInertiaM.w', cflags = '', simflags = ' -abortSlowSimulation -alarm=360 -emit_protected'
// Simulation options: startTime = 0.0, stopTime = 1.5, numberOfIntervals = 1500, tolerance = 1e-5, method = 'dassl', fileNamePrefix = 'Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMS_Start_MultiPhase', options = '', outputFormat = 'mat', variableFilter = 'time|aims3.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aims3.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aims3.rotor.zeroInductor.i0|aims3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aims3.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimsM.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimsM.rotor.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|aimsM.rotor.zeroInductor.i0|aimsM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.1..Phi.re|aimsM.stator.electroMagneticConverter.singlePhaseElectroMagneticConverter.2..Phi.im|loadInertia3.phi|loadInertia3.w|loadInertiaM.phi|loadInertiaM.w', cflags = '', simflags = ' -abortSlowSimulation -alarm=360 -emit_protected'
// Result file: Modelica.Magnetic.FundamentalWave.Examples.BasicMachines.AIMS_Start_MultiPhase_res.mat
// Messages: LOG_SUCCESS | info | The initialization finished successfully without homotopy method.
// LOG_SUCCESS | info | The simulation finished successfully.
Expand Down

0 comments on commit 6c2fe7e

Please sign in to comment.