Skip to content
This repository was archived by the owner on May 18, 2019. It is now read-only.

Commit d55e7e1

Browse files
ptaeuberOpenModelica-Hudson
ptaeuber
authored and
OpenModelica-Hudson
committed
Allow adaptive step size also for local homotopy
Belonging to [master]: - #2014
1 parent e82a713 commit d55e7e1

File tree

6 files changed

+143
-52
lines changed

6 files changed

+143
-52
lines changed

‎Compiler/BackEnd/BackendDAEOptimize.mo‎

Lines changed: 88 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -5706,10 +5706,15 @@ protected
57065706
BackendDAE.StrongComponents comps;
57075707
BackendDAE.EqSystems newEqSystems={};
57085708
array<Integer> ass1, ass2;
5709+
Boolean globalHomotopy = stringEq(Flags.getConfigString(Flags.HOMOTOPY_APPROACH), "global");
57095710
algorithm
57105711
for syst in outDAE.eqs loop
57115712
BackendDAE.MATCHING(ass1=ass1, ass2=ass2, comps=comps) := syst.matching;
5712-
(comps, syst) := traverseStrongComponentsForHomotopyLoop(comps, syst);
5713+
if globalHomotopy then
5714+
(comps, syst) := traverseStrongComponentsForHomotopyLoop(comps, syst);
5715+
else
5716+
(comps, syst) := traverseStrongComponentsAddLambda(comps, syst);
5717+
end if;
57135718
syst.matching := BackendDAE.MATCHING(ass1=ass1, ass2=ass2, comps=comps);
57145719
newEqSystems := syst::newEqSystems;
57155720
end for;
@@ -5750,7 +5755,7 @@ algorithm
57505755
homotopyLoopBeginning = compIndex;
57515756
end if;
57525757
end if;
5753-
then();
5758+
then();
57545759

57555760
case(BackendDAE.EQUATIONSYSTEM(eqns=eqnIndexes))
57565761
equation
@@ -5765,7 +5770,7 @@ algorithm
57655770
else
57665771
homotopyLoopEnd = compIndex;
57675772
end if;
5768-
then();
5773+
then();
57695774

57705775
case(BackendDAE.SINGLEARRAY(eqn=eqnIndex))
57715776
equation
@@ -5778,7 +5783,7 @@ algorithm
57785783
homotopyLoopBeginning = compIndex;
57795784
end if;
57805785
end if;
5781-
then();
5786+
then();
57825787

57835788
case(BackendDAE.SINGLEALGORITHM(eqn=eqnIndex))
57845789
equation
@@ -5791,7 +5796,7 @@ algorithm
57915796
homotopyLoopBeginning = compIndex;
57925797
end if;
57935798
end if;
5794-
then();
5799+
then();
57955800

57965801
case(BackendDAE.SINGLECOMPLEXEQUATION(eqn=eqnIndex))
57975802
equation
@@ -5804,7 +5809,7 @@ algorithm
58045809
homotopyLoopBeginning = compIndex;
58055810
end if;
58065811
end if;
5807-
then();
5812+
then();
58085813

58095814
case(BackendDAE.SINGLEWHENEQUATION(eqn=eqnIndex))
58105815
equation
@@ -5817,7 +5822,7 @@ algorithm
58175822
homotopyLoopBeginning = compIndex;
58185823
end if;
58195824
end if;
5820-
then();
5825+
then();
58215826

58225827
case(BackendDAE.SINGLEIFEQUATION(eqn=eqnIndex))
58235828
equation
@@ -5830,7 +5835,7 @@ algorithm
58305835
homotopyLoopBeginning = compIndex;
58315836
end if;
58325837
end if;
5833-
then();
5838+
then();
58345839

58355840
case(BackendDAE.TORNSYSTEM(strictTearingSet=BackendDAE.TEARINGSET(residualequations=resEqnIndexes, innerEquations=innerEquations)))
58365841
equation
@@ -5851,7 +5856,7 @@ algorithm
58515856
else
58525857
homotopyLoopEnd = compIndex;
58535858
end if;
5854-
then();
5859+
then();
58555860
end match;
58565861
end for;
58575862

@@ -5897,7 +5902,7 @@ algorithm
58975902
else
58985903
outHomotopyComponents = comp::outHomotopyComponents;
58995904
end if;
5900-
then (listReverse(outPreHomotopyComponents), listReverse(outHomotopyComponents), listReverse(outPostHomotopyComponents));
5905+
then (listReverse(outPreHomotopyComponents), listReverse(outHomotopyComponents), listReverse(outPostHomotopyComponents));
59015906

59025907
case(i::indexes, comp::comps)
59035908
equation
@@ -5908,12 +5913,12 @@ algorithm
59085913
else
59095914
outHomotopyComponents = comp::outHomotopyComponents;
59105915
end if;
5911-
then getHomotopyComponents(indexes, comps, homotopyLoopBeginning, homotopyLoopEnd, outPreHomotopyComponents, outHomotopyComponents, outPostHomotopyComponents);
5916+
then getHomotopyComponents(indexes, comps, homotopyLoopBeginning, homotopyLoopEnd, outPreHomotopyComponents, outHomotopyComponents, outPostHomotopyComponents);
59125917
end match;
59135918
end getHomotopyComponents;
59145919

59155920
protected function createOneHomotopyComponent " creates one BackendDAE.TORNSYSTEM out of all homotopy components
5916-
author: ptaeuber"
5921+
author: ptaeuber 2017"
59175922
input BackendDAE.StrongComponents homotopyComponents;
59185923
input BackendDAE.EqSystem inSystem;
59195924
input Integer lambdaIdx;
@@ -5936,39 +5941,39 @@ algorithm
59365941
case(BackendDAE.SINGLEEQUATION(eqn=eqnIndex, var=varIndex))
59375942
equation
59385943
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars={varIndex});
5939-
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
5944+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
59405945

59415946
case(BackendDAE.EQUATIONSYSTEM(eqns=eqnIndexes, vars=varIndexes, mixedSystem=mixedSystem))
59425947
equation
59435948
if mixedSystem then
59445949
isMixed = true;
59455950
end if;
5946-
then (newInnerEquations, listAppend(newResEquations, eqnIndexes), listAppend(newIterationVars, varIndexes));
5951+
then (newInnerEquations, listAppend(newResEquations, eqnIndexes), listAppend(newIterationVars, varIndexes));
59475952

59485953
case(BackendDAE.SINGLEARRAY(eqn=eqnIndex, vars=varIndexes))
59495954
equation
59505955
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
5951-
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
5956+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
59525957

59535958
case(BackendDAE.SINGLEALGORITHM(eqn=eqnIndex, vars=varIndexes))
59545959
equation
59555960
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
5956-
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
5961+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
59575962

59585963
case(BackendDAE.SINGLECOMPLEXEQUATION(eqn=eqnIndex, vars=varIndexes))
59595964
equation
59605965
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
5961-
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
5966+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
59625967

59635968
case(BackendDAE.SINGLEWHENEQUATION(eqn=eqnIndex, vars=varIndexes))
59645969
equation
59655970
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
5966-
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
5971+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
59675972

59685973
case(BackendDAE.SINGLEIFEQUATION(eqn=eqnIndex, vars=varIndexes))
59695974
equation
59705975
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
5971-
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
5976+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
59725977

59735978
case(BackendDAE.TORNSYSTEM(strictTearingSet=BackendDAE.TEARINGSET(residualequations=resEqnIndexes, tearingvars=tVarIndexes, innerEquations=innerEquations), mixedSystem=mixedSystem))
59745979
algorithm
@@ -5978,12 +5983,75 @@ algorithm
59785983
for innerEquation in innerEquations loop
59795984
newInnerEquations := innerEquation::newInnerEquations;
59805985
end for;
5981-
then (newInnerEquations, listAppend(newResEquations, resEqnIndexes), listAppend(newIterationVars, tVarIndexes));
5986+
then (newInnerEquations, listAppend(newResEquations, resEqnIndexes), listAppend(newIterationVars, tVarIndexes));
59825987
end match;
59835988
end for;
59845989

59855990
outHomotopyComponent := BackendDAE.TORNSYSTEM(BackendDAE.TEARINGSET(listAppend(newIterationVars,{lambdaIdx}), newResEquations, listReverse(newInnerEquations), BackendDAE.EMPTY_JACOBIAN()), NONE(), false, isMixed);
59865991
end createOneHomotopyComponent;
59875992

5993+
protected function traverseStrongComponentsAddLambda " traverses all the strong components and adds lambda as the last variable if the system contains homotopy
5994+
author: ptaeuber 2017"
5995+
input output BackendDAE.StrongComponents comps;
5996+
input output BackendDAE.EqSystem system;
5997+
protected
5998+
BackendDAE.StrongComponents newComps={};
5999+
BackendDAE.Var lambda;
6000+
Integer lambdaIdx;
6001+
Boolean hasAnyHomotopy = false;
6002+
algorithm
6003+
lambdaIdx := BackendVariable.varsSize(system.orderedVars) + 1;
6004+
6005+
for comp in comps loop
6006+
comp := match(comp)
6007+
local
6008+
list<Integer> eqnIndexes, varIndexes, resEqnIndexes, tVarIndexes, innerEqnIndexes;
6009+
Option<BackendDAE.TearingSet> casualTearingSet;
6010+
BackendDAE.InnerEquations innerEquations;
6011+
BackendDAE.Jacobian jac;
6012+
BackendDAE.JacobianType jacType;
6013+
list<BackendDAE.Equation> eqnLst;
6014+
Boolean linear, mixedSystem, hasHomotopy;
6015+
6016+
case(BackendDAE.EQUATIONSYSTEM(eqns=eqnIndexes, vars=varIndexes, jac=jac, jacType=jacType, mixedSystem=mixedSystem))
6017+
equation
6018+
eqnLst = BackendEquation.getList(eqnIndexes, system.orderedEqs);
6019+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquationList(eqnLst, BackendDAEUtil.containsHomotopyCall, false);
6020+
6021+
if hasHomotopy then
6022+
hasAnyHomotopy = true;
6023+
comp = BackendDAE.EQUATIONSYSTEM(eqnIndexes, varIndexes, jac, jacType, mixedSystem);
6024+
end if;
6025+
then comp;
6026+
6027+
case(BackendDAE.TORNSYSTEM(strictTearingSet=BackendDAE.TEARINGSET(residualequations=resEqnIndexes, tearingvars=tVarIndexes, innerEquations=innerEquations, jac=jac), casualTearingSet=casualTearingSet, linear=linear, mixedSystem=mixedSystem))
6028+
equation
6029+
eqnLst = BackendEquation.getList(resEqnIndexes, system.orderedEqs);
6030+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquationList(eqnLst, BackendDAEUtil.containsHomotopyCall, false);
6031+
6032+
if not hasHomotopy then
6033+
(innerEqnIndexes,_,_) = List.map_3(innerEquations, BackendDAEUtil.getEqnAndVarsFromInnerEquation);
6034+
eqnLst = BackendEquation.getList(innerEqnIndexes, system.orderedEqs);
6035+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquationList(eqnLst, BackendDAEUtil.containsHomotopyCall, false);
6036+
end if;
6037+
6038+
if hasHomotopy then
6039+
hasAnyHomotopy = true;
6040+
comp = BackendDAE.TORNSYSTEM(BackendDAE.TEARINGSET(listAppend(tVarIndexes, {lambdaIdx}), resEqnIndexes, innerEquations, jac), casualTearingSet, linear, mixedSystem);
6041+
end if;
6042+
then comp;
6043+
else comp;
6044+
end match;
6045+
newComps := comp::newComps;
6046+
end for;
6047+
6048+
if hasAnyHomotopy then
6049+
// Add homotopy lambda to system
6050+
lambda := BackendDAE.VAR(ComponentReference.makeCrefIdent(BackendDAE.homotopyLambda, DAE.T_REAL_DEFAULT, {}), BackendDAE.VARIABLE(), DAE.BIDIR(), DAE.NON_PARALLEL(), DAE.T_REAL_DEFAULT, NONE(), NONE(), {}, DAE.emptyElementSource, NONE(), NONE(), DAE.BCONST(false), NONE(), DAE.NON_CONNECTOR(), DAE.NOT_INNER_OUTER(), true);
6051+
system.orderedVars := BackendVariable.addVar(lambda, system.orderedVars);
6052+
end if;
6053+
comps := listReverse(newComps);
6054+
end traverseStrongComponentsAddLambda;
6055+
59886056
annotation(__OpenModelica_Interface="backend");
59896057
end BackendDAEOptimize;

‎Compiler/Template/CodegenC.tpl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1080,7 +1080,7 @@ template simulationFile(SimCode simCode, String guid, Boolean isModelExchangeFMU
10801080
<%symbolName(modelNamePrefixStr,"function_storeDelayed")%>,
10811081
<%symbolName(modelNamePrefixStr,"updateBoundVariableAttributes")%>,
10821082
<%symbolName(modelNamePrefixStr,"functionInitialEquations")%>,
1083-
<%if listEmpty(initialEquations_lambda0) then '0' else if not BackendDAEUtil.isInitOptModuleActivated("generateHomotopyComponents") then '1' else '2'%>, /* useHomotopy - 0: no homotopy or local homotopy, 1: global homotopy, 2: new global homotopy approach */
1083+
<%if listEmpty(initialEquations_lambda0) then (if not BackendDAEUtil.isInitOptModuleActivated("generateHomotopyComponents") then '0' else '3') else if not BackendDAEUtil.isInitOptModuleActivated("generateHomotopyComponents") then '1' else '2'%>, /* useHomotopy - 0: no homotopy or local homotopy (equidistant lambda), 1: global homotopy (equidistant lambda), 2: new global homotopy approach, 3: new local homotopy approach */
10841084
<%symbolName(modelNamePrefixStr,"functionInitialEquations_lambda0")%>,
10851085
<%symbolName(modelNamePrefixStr,"functionRemovedInitialEquations")%>,
10861086
<%symbolName(modelNamePrefixStr,"updateBoundParameters")%>,

‎SimulationRuntime/c/openmodelica_func.h‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -145,9 +145,10 @@ int (*functionInitialEquations)(DATA *data, threadData_t*);
145145

146146
/*! \var useHomotopy
147147
*
148-
* 0: local homotopy approach
149-
* 1: global homotopy approach
150-
* 2: new global homotopy approach
148+
* 0: local homotopy approach (equidistant lambda steps)
149+
* 1: global homotopy approach (equidistant lambda steps)
150+
* 2: new global homotopy approach (adaptive lambda steps)
151+
* 3: new local homotopy approach (adaptive lambda steps)
151152
*/
152153
const int useHomotopy;
153154

‎SimulationRuntime/c/simulation/solver/initialization/initialization.c‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -185,7 +185,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
185185
FILE *pFile = NULL;
186186
long i;
187187
MODEL_DATA *mData = data->modelData;
188-
int solveWithGlobalHomotopy = data->callback->useHomotopy == 0 || (data->callback->useHomotopy != 2 && init_lambda_steps < 2) ? 0 : 1;
188+
int solveWithGlobalHomotopy = (data->callback->useHomotopy == 1 && init_lambda_steps > 1) || data->callback->useHomotopy == 2;
189189

190190
#if !defined(OMC_NDELAY_EXPRESSIONS) || OMC_NDELAY_EXPRESSIONS>0
191191
/* initial sample and delay before initial the system */
@@ -197,8 +197,10 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
197197

198198
/* If there is no homotopy in the model or local homotopy is activated
199199
or homotopy is disabled by runtime flag '-ils=<lambda_steps>',
200-
solve WITHOUT HOMOTOPY. */
200+
solve WITHOUT HOMOTOPY or LOCAL HOMOTOPY. */
201201
if (!solveWithGlobalHomotopy){
202+
if (data->callback->useHomotopy == 3 && !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
203+
infoStreamPrint(LOG_INIT, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the adaptive local approach yet.");
202204
data->simulationInfo->lambda = 1.0;
203205
data->callback->functionInitialEquations(data, threadData);
204206

@@ -247,6 +249,7 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
247249
}
248250
#endif
249251

252+
infoStreamPrint(LOG_INIT, 0, "Global homotopy with equidistant step size started.");
250253
infoStreamPrint(LOG_INIT, 1, "homotopy process\n---------------------------");
251254
for(step=0; step<init_lambda_steps; ++step)
252255
{
@@ -288,6 +291,10 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
288291
use NEW GLOBAL HOMOTOPY APPROACH. */
289292
if (data->callback->useHomotopy == 2 && solveWithGlobalHomotopy)
290293
{
294+
if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
295+
infoStreamPrint(LOG_INIT, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the adaptive global approach yet.");
296+
297+
infoStreamPrint(LOG_INIT, 0, "Global homotopy with adaptive step size started.");
291298
infoStreamPrint(LOG_INIT, 1, "homotopy process\n---------------------------");
292299

293300
// Solve lambda0-DAE

‎SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1653,7 +1653,7 @@ static int homotopyAlgorithm(DATA_HOMOTOPY* solverData, double *x)
16531653
if (iter>=maxTries)
16541654
{
16551655
if (solverData->initHomotopy)
1656-
warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe maximum number of tries for one lamda is reached (%d).\nYou can change the number of tries with:\n\t-homMaxTries=<value>\nYou can also try to allow more newton steps in the corrector step with:\n\t-homMaxNewtonSteps=<value>\nor change the tolerance for the solution with:\n\t-homHEps=<value>\nYou can use -lv=LOG_INIT,LOG_NLS_HOMOTOPY to get more information.", iter);
1656+
warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe maximum number of tries for one lambda is reached (%d).\nYou can change the number of tries with:\n\t-homMaxTries=<value>\nYou can also try to allow more newton steps in the corrector step with:\n\t-homMaxNewtonSteps=<value>\nor change the tolerance for the solution with:\n\t-homHEps=<value>\nYou can use -lv=LOG_INIT,LOG_NLS_HOMOTOPY to get more information.", iter);
16571657
else
16581658
debugInt(LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge: iter = ", iter);
16591659
debugString(LOG_NLS_HOMOTOPY, "======================================================");
@@ -1787,7 +1787,9 @@ static int homotopyAlgorithm(DATA_HOMOTOPY* solverData, double *x)
17871787

17881788
solverData->tau = tau;
17891789
printHomotopyPredictorStep(LOG_NLS_HOMOTOPY, solverData);
1790+
17901791
/* Corrector step: Newton iteration! */
1792+
debugString(LOG_NLS_HOMOTOPY, "Newton iteration for corrector step begins!");
17911793
for(j=0;j<maxiter;j++)
17921794
{
17931795
if (vec2Norm(solverData->n, solverData->hvec)<hEps || vec2Norm(solverData->n, solverData->hvecScaled)<hEps)
@@ -1864,6 +1866,7 @@ static int homotopyAlgorithm(DATA_HOMOTOPY* solverData, double *x)
18641866
debugDouble(LOG_NLS_HOMOTOPY, "hEps =", hEps);
18651867

18661868
}
1869+
debugString(LOG_NLS_HOMOTOPY, "Newton iteration for corrector step finished!");
18671870

18681871
if (!assert)
18691872
{
@@ -1985,7 +1988,7 @@ int solveHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
19851988
int numberOfFunctionEvaluationsOld = solverData->numberOfFunctionEvaluations;
19861989
solverData->casualTearingSet = systemData->strictTearingFunctionCall != NULL;
19871990
int constraintViolated;
1988-
solverData->initHomotopy = data->callback->useHomotopy == 2 && systemData->homotopySupport;
1991+
solverData->initHomotopy = (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3) && systemData->homotopySupport;
19891992

19901993
modelica_boolean* relationsPreBackup;
19911994
relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));

0 commit comments

Comments
 (0)