-
Notifications
You must be signed in to change notification settings - Fork 298
/
NBTearing.mo
524 lines (471 loc) · 21.1 KB
/
NBTearing.mo
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
/*
* This file is part of OpenModelica.
*
* Copyright (c) 1998-2020, Open Source Modelica Consortium (OSMC),
* c/o Linköpings universitet, Department of Computer and Information Science,
* SE-58183 Linköping, Sweden.
*
* All rights reserved.
*
* THIS PROGRAM IS PROVIDED UNDER THE TERMS OF GPL VERSION 3 LICENSE OR
* THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
* ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
* RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
* ACCORDING TO RECIPIENTS CHOICE.
*
* The OpenModelica software and the Open Source Modelica
* Consortium (OSMC) Public License (OSMC-PL) are obtained
* from OSMC, either from the above address,
* from the URLs: http://www.ida.liu.se/projects/OpenModelica or
* http://www.openmodelica.org, and in the OpenModelica distribution.
* GNU version 3 is obtained from: http://www.gnu.org/copyleft/gpl.html.
*
* This program is distributed WITHOUT ANY WARRANTY; without
* even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
* IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL.
*
* See the full OSMC Public License conditions for more details.
*
*/
encapsulated uniontype NBTearing
"file: NBTearing.mo
package: NBTearing
description: This file contains the data-types used for tearing. It is a
uniontype and therefore also contains some structures for tearing.
"
public
import BackendDAE = NBackendDAE;
import Module = NBModule;
import Slice = NBSlice;
import NBVariable.{VariablePointer, VariablePointers};
import NBEquation.{Equation, EquationPointer, EquationPointers};
import StrongComponent = NBStrongComponent;
protected
// selfimport
import Tearing = NBTearing;
// NF imports
import NFFlatten.FunctionTree;
import Variable = NFVariable;
// Backend imports
import Adjacency = NBAdjacency;
import BEquation = NBEquation;
import BJacobian = NBJacobian;
import BVariable = NBVariable;
import Causalize = NBCausalize;
import Differentiate = NBDifferentiate;
import Inline = NBInline;
import Jacobian = NBackendDAE.BackendDAE;
import Matching = NBMatching;
import Sorting = NBSorting;
import System = NBSystem;
//Util imports
import BackendUtil = NBBackendUtil;
import StringUtil;
public
record TEARING_SET
list<Slice<VariablePointer>> iteration_vars "the variables used for iteration";
list<Slice<EquationPointer>> residual_eqns "implicitely solved residual equations";
array<StrongComponent> innerEquations "array of matched equations and variables";
Option<Jacobian> jac "optional jacobian";
end TEARING_SET;
function hash
input Tearing tearing;
output Integer i = -1;
algorithm
// ToDo!
end hash;
function isEqual
input Tearing tearing1;
input Tearing tearing2;
output Boolean b = false;
algorithm
// ToDo!
end isEqual;
function toString
input Tearing set;
input output String str;
algorithm
str := StringUtil.headline_4(str);
str := str + "### Iteration Variables:\n" + Slice.lstToString(set.iteration_vars, BVariable.pointerToString);
str := str + "\n### Residual Equations:\n" + Slice.lstToString(set.residual_eqns, function Equation.pointerToString(str = ""));
str := str + "\n### Inner Equations:\n" + List.toString(arrayList(set.innerEquations), function StrongComponent.toString(index = -1), "", "\t", "\n\t", "");
if Util.isSome(set.jac) then
str := str + "\n" + BJacobian.toString(Util.getOption(set.jac), "NLS");
end if;
end toString;
function main
"Wrapper function for any tearing function. This will be
called during simulation and gets the corresponding subfunction from
Config."
extends Module.wrapper;
input System.SystemType systemType;
protected
constant list<Module.tearingInterface> funcs = getModule();
FunctionTree funcTree;
algorithm
if Flags.isSet(Flags.TEARING_DUMP) then
print(StringUtil.headline_1("[" + System.System.systemTypeString(systemType) + "] Tearing") + "\n");
end if;
bdae := match (systemType, bdae)
local
list<System.System> systems;
VariablePointers variables;
Pointer<Integer> eq_index;
case (NBSystem.SystemType.ODE, BackendDAE.MAIN(ode = systems, funcTree = funcTree, varData = BVariable.VAR_DATA_SIM(variables = variables), eqData = BEquation.EQ_DATA_SIM(uniqueIndex = eq_index)))
algorithm
(systems, funcTree) := tearingTraverser(systems, funcs, funcTree, variables, eq_index, systemType);
bdae.ode := systems;
bdae.funcTree := funcTree;
then bdae;
case (NBSystem.SystemType.INI, BackendDAE.MAIN(init = systems, funcTree = funcTree, varData = BVariable.VAR_DATA_SIM(variables = variables), eqData = BEquation.EQ_DATA_SIM(uniqueIndex = eq_index)))
algorithm
(systems, funcTree) := tearingTraverser(systems, funcs, funcTree, variables, eq_index, systemType);
bdae.init := systems;
if Util.isSome(bdae.init_0) then
(systems, funcTree) := tearingTraverser(Util.getOption(bdae.init_0), funcs, funcTree, variables, eq_index, systemType);
bdae.init_0 := SOME(systems);
end if;
bdae.funcTree := funcTree;
then bdae;
case (NBSystem.SystemType.DAE, BackendDAE.MAIN(dae = SOME(systems), funcTree = funcTree, varData = BVariable.VAR_DATA_SIM(variables = variables), eqData = BEquation.EQ_DATA_SIM(uniqueIndex = eq_index)))
algorithm
(systems, funcTree) := tearingTraverser(systems, funcs, funcTree, variables, eq_index, systemType);
bdae.dae := SOME(systems);
bdae.funcTree := funcTree;
then bdae;
// ToDo: all the other cases: e.g. Jacobian, Hessian
end match;
end main;
function implicit
input output StrongComponent comp "the suspected algebraic loop.";
input output FunctionTree funcTree "Function call bodies";
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()
Adjacency.Matrix dummy = Adjacency.EMPTY(NBAdjacency.MatrixStrictness.FULL);
algorithm
(comp, dummy, funcTree, index) := match comp
// create implicit equations
case StrongComponent.SINGLE_COMPONENT()
then tearingNone(StrongComponent.ALGEBRAIC_LOOP(
idx = index,
strict = singleImplicit(comp.var, comp.eqn),
casual = NONE(),
linear = false,
mixed = false,
status = NBSolve.Status.IMPLICIT), dummy, funcTree, index, VariablePointers.empty(), Pointer.create(0), systemType);
case StrongComponent.MULTI_COMPONENT()
then tearingNone(StrongComponent.ALGEBRAIC_LOOP(
idx = index,
strict = singleImplicit(List.first(comp.vars), comp.eqn), // this is wrong! need to take all vars
casual = NONE(),
linear = false,
mixed = false,
status = NBSolve.Status.IMPLICIT), dummy, funcTree, index, VariablePointers.empty(), Pointer.create(0), systemType);
// do nothing otherwise
else (comp, dummy, funcTree, index);
end match;
end implicit;
function singleImplicit
input VariablePointer var;
input EquationPointer eqn;
output NBTearing tearingSet = Tearing.TEARING_SET(
iteration_vars = {Slice.SLICE(var, {})},
residual_eqns = {Slice.SLICE(eqn, {})},
innerEquations = listArray({}),
jac = NONE());
end singleImplicit;
function getModule
"Returns the module function that was chosen by the user."
output list<Module.tearingInterface> funcs;
protected
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};
/* ... New tearing modules have to be added here */
else fail();
end match;
end getModule;
function getResidualVars
input Tearing tearing;
output list<Pointer<Variable>> residuals;
algorithm
residuals := list(Equation.getResidualVar(Slice.getT(eqn)) for eqn in tearing.residual_eqns);
end getResidualVars;
function setResidualEqns
input output Tearing tearing;
input list<Slice<EquationPointer>> residuals;
algorithm
tearing.residual_eqns := residuals;
end setResidualEqns;
protected
// Traverser function
function tearingTraverser
input list<System.System> systems;
input list<Module.tearingInterface> funcs;
output list<System.System> new_systems = {};
input output FunctionTree funcTree;
input VariablePointers variables;
input Pointer<Integer> eq_index;
input System.SystemType systemType;
protected
array<StrongComponent> strongComponents;
StrongComponent tmp;
Integer idx = 0;
Adjacency.Matrix full "full adjacency matrix containing solvability info";
algorithm
for syst in systems loop
if isSome(syst.strongComponents) and isSome(syst.adjacencyMatrix) then
SOME(strongComponents) := syst.strongComponents;
SOME(full) := syst.adjacencyMatrix;
for i in 1:arrayLength(strongComponents) loop
// each module has a list of functions that need to be applied
tmp := strongComponents[i];
for func in funcs loop
(tmp, full, funcTree, idx) := func(tmp, full, funcTree, idx, variables, eq_index, systemType);
end for;
// only update if it changed
if not referenceEq(tmp, strongComponents[i]) then
arrayUpdate(strongComponents, i, tmp);
end if;
end for;
syst.strongComponents := SOME(strongComponents);
syst.adjacencyMatrix := SOME(full);
end if;
new_systems := syst :: new_systems;
end for;
new_systems := listReverse(new_systems);
end tearingTraverser;
function tearingFinalize extends Module.tearingInterface;
protected
list<StrongComponent> residual_comps;
Option<Jacobian> jacobian;
Tearing strict;
protected
list<Slice<EquationPointer>> tmp;
list<list<Slice<EquationPointer>>> acc = {};
String tag;
algorithm
(comp, 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);
end match;
end tearingFinalize;
function tearingNone extends Module.tearingInterface;
// does nothing but set index and call the jacobian
protected
list<StrongComponent> residual_comps;
Option<Jacobian> jacobian;
Tearing strict;
protected
list<Slice<EquationPointer>> tmp;
list<list<Slice<EquationPointer>>> acc = {};
algorithm
(comp, index) := match comp
case StrongComponent.ALGEBRAIC_LOOP(strict = strict) algorithm
index := index + 1;
comp.idx := index;
for eqn in listReverse(strict.residual_eqns) loop
tmp := Inline.inlineRecordSliceEquation(eqn, variables, eq_index, true);
if listEmpty(tmp) then
acc := {eqn} :: acc;
else
acc := tmp :: acc;
end if;
end for;
// create residual equations
strict.residual_eqns := list(Slice.apply(eqn, function Equation.createResidual(new = true)) for eqn in List.flatten(acc));
comp.strict := strict;
then (comp, index);
else (comp, index);
end match;
end tearingNone;
function tearingMinimal extends Module.tearingInterface;
// only extracts discrete variables to be solved as inner equations and calls jacobian module
protected
list<StrongComponent> residual_comps;
Option<Jacobian> jacobian;
Tearing strict;
list<Pointer<Variable>> vars_lst, cont_vars, disc_vars;
list<Pointer<Equation>> eqns_lst, cont_eqns, disc_eqns;
list<Slice<VariablePointer>> iter_lst;
list<Slice<EquationPointer>> residual_lst;
VariablePointers discreteVars;
EquationPointers eqns;
Adjacency.Matrix adj;
Matching matching;
list<StrongComponent> inner_comps, residual_comps;
algorithm
//print("######## minimal ########\n");
(comp, index) := match comp
case StrongComponent.ALGEBRAIC_LOOP(strict = strict) algorithm
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");
//print(List.toString(disc_eqns, function Equation.pointerToString(str = "")) + "\n");
/*
// 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);
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");
end if;
then (comp, index);
else (comp, index);
end match;
end tearingMinimal;
/*
function tearingMinimal extends Module.tearingInterface;
protected
list<StrongComponent> residual_comps;
Option<Jacobian> jacobian;
StrongComponent new_comp;
list<Slice<EquationPointer>> residuals = {};
list<Pointer<Variable>> cont_lst, disc_lst;
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;
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);
end match;
end tearingMinimal;
function tearingMinimalWork
input String name;
input list<Pointer<Variable>> variables;
input list<Pointer<Equation>> equations;
input Boolean mixed;
output StrongComponent comp;
input output FunctionTree funcTree;
input output Integer index;
protected
list<Pointer<Variable>> cont_lst, disc_lst;
list<Slice<VariablePointer>> iteration_lst;
list<Slice<EquationPointer>> residual_lst;
VariablePointers discreteVars;
EquationPointers eqns;
Adjacency.Matrix adj;
Matching matching;
list<StrongComponent> inner_comps, residual_comps;
Tearing tearingSet;
Option<Jacobian> jacobian;
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;