Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 69.0% 863 / 0 / 1251
Functions: -% 0 / 1 / 1
Branches: 61.2% 452 / 0 / 739
OMCompiler/Compiler/BackEnd/SymbolicJacobian.mo
Line Branch Exec Source
1 /*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF AGPL VERSION 3 LICENSE OR
11 * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GNU AGPL
14 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the OSMC (Open Source Modelica Consortium)
17 * Public License (OSMC-PL) are obtained from OSMC, either from the above
18 * address, from the URLs:
19 * http://www.openmodelica.org or
20 * https://github.com/OpenModelica/ or
21 * http://www.ida.liu.se/projects/OpenModelica,
22 * and in the OpenModelica distribution.
23 *
24 * GNU AGPL version 3 is obtained from:
25 * https://www.gnu.org/licenses/licenses.html#GPL
26 *
27 * This program is distributed WITHOUT ANY WARRANTY; without
28 * even the implied warranty of MERCHANTABILITY or FITNESS
29 * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
30 * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL.
31 *
32 * See the full OSMC Public License conditions for more details.
33 *
34 */
35
36 encapsulated package SymbolicJacobian
37 " file: SymbolicJacobian.mo
38 package: SymbolicJacobian
39 description: This package contains stuff that is related to symbolic jacobian or sparsity structure."
40
41
42 public import Absyn;
43 public import BackendDAE;
44 public import DAE;
45 public import FCore;
46 public import FGraph;
47
48 protected
49 import Array;
50 import BackendDAEOptimize;
51 import BackendDAETransform;
52 import BackendDAEUtil;
53 import BackendDump;
54 import BackendEquation;
55 import Coloring;
56 import BackendVariable;
57 import BackendVarTransform;
58 import BaseHashSet;
59 import Ceval;
60 import ClockIndexes;
61 import Config;
62 import ComponentReference;
63 protected import ComponentReferenceBasics;
64 import Debug;
65 import Differentiate;
66 import DynamicOptimization;
67 import ElementSource;
68 import ExecStat.execStat;
69 import ExpandableArray;
70 import Expression;
71 protected import ExpressionBasics;
72 import ExpressionDump;
73 import ExpressionSimplify;
74 import Error;
75 import Flags;
76 import FlagsUtil;
77 import GCExt;
78 import Global;
79 import Graph;
80 import HashSet;
81 import IndexReduction;
82 import List;
83 import StringUtil;
84 import System;
85 import UnorderedMap;
86 import UnorderedSet;
87 import Util;
88 import Values;
89 import ValuesUtil;
90
91 // =============================================================================
92 // section for postOptModule >>symbolicJacobian<<
93 //
94 // Detects the sparse pattern of the ODE system and calculates also the symbolic
95 // Jacobian if flag "--generateDynamicJacobian=symbolic".
96 // =============================================================================
97
98 // From User Documentation for ida v5.4.0 equation (2.5) aka Alpha
99 // is the scalar in the system Jacobian, proportional to the inverse of the step
100 // size used for DAE_Mode symbolic jacobians
101 public constant String DAE_CJ = "$DAE_CJ";
102
103 public function symbolicJacobian "author: lochel
104 Detects the sparse pattern of the ODE system and calculates also the symbolic
105 Jacobian if flag '--generateDynamicJacobian=symbolic'."
106 input BackendDAE.BackendDAE inDAE;
107 output BackendDAE.BackendDAE outDAE;
108 algorithm
109 outDAE := match Flags.getConfigString(Flags.GENERATE_DYNAMIC_JACOBIAN)
110 case "none" then inDAE;
111 1055 case "numeric" then detectSparsePatternODE(inDAE);
112 6 case "symbolic" then generateSymbolicJacobianPast(inDAE);
113 end match;
114 end symbolicJacobian;
115
116 // =============================================================================
117 // section for postOptModule >>calculateStateSetsJacobians<<
118 //
119 // =============================================================================
120
121 public function calculateStateSetsJacobians "author: wbraun
122 Calculates the Jacobian matrix with directional derivative method for dynamic
123 state selection."
124 input BackendDAE.BackendDAE inDAE;
125 output BackendDAE.BackendDAE outDAE;
126 algorithm
127 2138 outDAE := BackendDAEUtil.mapEqSystem(inDAE, calculateEqSystemStateSetsJacobians);
128 end calculateStateSetsJacobians;
129
130 // =============================================================================
131 // section for postOptModule >>calculateStrongComponentJacobians<<
132 //
133 // Module for to calculate strong component Jacobian matrices
134 // =============================================================================
135
136 public function calculateStrongComponentJacobians "author: wbraun
137 Calculates Jacobian matrix with directional derivative method for each SCC."
138 input BackendDAE.BackendDAE inDAE;
139 output BackendDAE.BackendDAE outDAE;
140 algorithm
141 try
142 3935 outDAE := BackendDAEUtil.mapEqSystem(inDAE, calculateEqSystemJacobians);
143 else
144 outDAE := inDAE;
145 end try;
146 end calculateStrongComponentJacobians;
147
148 // =============================================================================
149 // section for postOptModule >>constantLinearSystem<<
150 //
151 // constant Jacobian matrices. Linear system of equations (A x = b) where
152 // A and b are constant.
153 // =============================================================================
154
155 public function constantLinearSystem
156 input BackendDAE.BackendDAE inDAE;
157 output BackendDAE.BackendDAE outDAE;
158 algorithm
159 2815 (outDAE, _) := BackendDAEUtil.mapEqSystemAndFold(inDAE, constantLinearSystem0, (false,1));
160 end constantLinearSystem;
161
162 // =============================================================================
163 // section for postOptModule >>detectSparsePatternODE<<
164 //
165 // Generate sparse pattern
166 // =============================================================================
167 protected function detectSparsePatternODE
168 input BackendDAE.BackendDAE inBackendDAE;
169 output BackendDAE.BackendDAE outBackendDAE;
170 protected
171 BackendDAE.BackendDAE DAE;
172 BackendDAE.EqSystems eqs;
173 BackendDAE.Shared shared;
174 BackendDAE.SparseColoring coloredCols;
175 BackendDAE.SparsePattern sparsePattern;
176 list<BackendDAE.Var> states;
177 BackendDAE.Variables v;
178 constant Boolean debug = false;
179 algorithm
180 // lochel: This module fails for some models (e.g. #3543)
181 try
182 if debug then execStat("detectSparsePatternODE -> start "); end if;
183 1055 BackendDAE.DAE(eqs = eqs) := inBackendDAE;
184
185 // prepare a DAE
186 1055 DAE := BackendDAEUtil.copyBackendDAE(inBackendDAE);
187 if debug then execStat("detectSparsePatternODE -> copy dae "); end if;
188 1055 DAE := BackendDAEOptimize.collapseIndependentContinuousBlocks(DAE);
189 if debug then execStat("detectSparsePatternODE -> collapse blocks "); end if;
190 1055 DAE := BackendDAEUtil.transformBackendDAE(DAE, SOME((BackendDAE.NO_INDEX_REDUCTION(), BackendDAE.EXACT())), NONE(), NONE());
191 if debug then execStat("detectSparsePatternODE -> transform backend dae "); end if;
192
193 // get states for DAE
194
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1055 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1055 times.
 1055 BackendDAE.DAE(eqs = {BackendDAE.EQSYSTEM(orderedVars = v)}, shared=shared) := DAE;
195 1055 states := BackendVariable.getAllStateVarFromVariables(v);
196 if debug then execStat("detectSparsePatternODE -> get all vars "); end if;
197
198 // generate sparse pattern
199 1055 (sparsePattern, coloredCols) := generateSparsePattern(DAE, states, states);
200 if debug then execStat("detectSparsePatternODE -> generateSparsePattern "); end if;
201 1055 shared := addBackendDAESharedJacobianSparsePattern(sparsePattern, coloredCols, BackendDAE.SymbolicJacobianAIndex, shared);
202 if debug then execStat("detectSparsePatternODE -> addBackendDAESharedJacobianSparsePattern "); end if;
203
204 1055 outBackendDAE := BackendDAE.DAE(eqs, shared);
205 else
206 // skip this optimization module
207 ✗ Error.addCompilerWarning("The optimization module detectJacobianSparsePattern failed. This module will be skipped and the transformation process continued.");
208 outBackendDAE := inBackendDAE;
209 end try;
210 end detectSparsePatternODE;
211
212 // =============================================================================
213 // section for postOptModule >>symbolicJacobianDAE<<
214 //
215 // Generate symbolic jacobian for DAEMode
216 // =============================================================================
217
218 public function symbolicJacobianDAE
219 input BackendDAE.BackendDAE inBackendDAE;
220 output BackendDAE.BackendDAE outBackendDAE;
221 protected
222 BackendDAE.BackendDAE DAE;
223 BackendDAE.EqSystems eqs;
224 BackendDAE.Shared shared;
225 BackendDAE.SparseColoring coloredCols;
226 BackendDAE.SparsePattern sparsePattern;
227 BackendDAE.NonlinearPattern nonlinearPattern;
228 list<BackendDAE.Var> inDepVars;
229 list<BackendDAE.Var> depVars;
230 BackendDAE.Variables v, resVars;
231 BackendDAE.Variables emptyVars = BackendVariable.emptyVars();
232 Option<BackendDAE.SymbolicJacobian> symjac;
233 AvlTreePathFunction.Tree funcs;
234 constant Boolean debug = false;
235 algorithm
236 try
237 if debug then execStat(getInstanceName() + "-> start "); end if;
238 9 BackendDAE.DAE(eqs = eqs) := inBackendDAE;
239
240 // prepare a DAE
241 9 DAE := BackendDAEUtil.copyBackendDAE(inBackendDAE);
242 if debug then execStat(getInstanceName() + "-> copy dae "); end if;
243 9 DAE := BackendDAEOptimize.collapseIndependentBlocks(DAE);
244 if debug then execStat(getInstanceName() + "-> collapse blocks "); end if;
245 9 DAE := BackendDAEUtil.transformBackendDAE(DAE, SOME((BackendDAE.NO_INDEX_REDUCTION(), BackendDAE.EXACT())), NONE(), NONE());
246 if debug then execStat(getInstanceName() + "-> transform backend dae "); end if;
247
248 // get states for DAE
249
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
 9 BackendDAE.DAE(eqs = {BackendDAE.EQSYSTEM(orderedVars = v)}, shared=shared) := DAE;
250 9 (_, resVars) := BackendVariable.traverseBackendDAEVars(v, BackendVariable.collectVarKindVarinVariables, (BackendVariable.isDAEmodeResVar, emptyVars));
251 9 depVars := BackendVariable.varList(resVars);
252
253 9 inDepVars := listAppend(shared.daeModeData.stateVars, shared.daeModeData.algStateVars);
254
255 if debug then execStat(getInstanceName() + "-> get all vars "); end if;
256
257
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 9 if Flags.getConfigString(Flags.GENERATE_DYNAMIC_JACOBIAN) == "symbolic" then
258 // generate symbolic jacobian and sparsity pattern
259 ✗ (symjac, funcs, sparsePattern, coloredCols, nonlinearPattern) := generateGenericJacobian(
260 inBackendDAE = DAE,
261 inDiffVars = inDepVars,
262 inStateVars = BackendVariable.emptyVars(),
263 inInputVars = BackendVariable.emptyVars(),
264 inParameterVars = shared.globalKnownVars,
265 inDifferentiatedVars = resVars,
266 inVars = BackendVariable.varList(v),
267 inName = "A",
268 onlySparsePattern = false,
269 daeMode = true);
270 if debug then execStat(getInstanceName() + "-> generateGenericJacobian "); end if;
271
272 ✗ shared.symjacs := List.set(shared.symjacs, BackendDAE.SymbolicJacobianAIndex, (symjac, sparsePattern, coloredCols, nonlinearPattern));
273 ✗ shared.functionTree := funcs;
274
275 if debug then BackendDump.dumpJacobianString(BackendDAE.GENERIC_JACOBIAN(symjac, sparsePattern, coloredCols, nonlinearPattern)); end if;
276 else
277 // only generate sparsity pattern
278 9 (sparsePattern, coloredCols) := generateSparsePattern(DAE, inDepVars, depVars);
279 if debug then execStat(getInstanceName() + "-> generateSparsePattern "); end if;
280 9 shared := addBackendDAESharedJacobianSparsePattern(sparsePattern, coloredCols, BackendDAE.SymbolicJacobianAIndex, shared);
281 if debug then execStat(getInstanceName() + "-> addBackendDAESharedJacobianSparsePattern "); end if;
282 end if;
283
284 9 outBackendDAE := BackendDAE.DAE(eqs, shared);
285 else
286 // skip this optimization module
287 ✗ Error.addCompilerWarning("The optimization module " + getInstanceName() + " failed. This module will be skipped and the transformation process continued.");
288 outBackendDAE := inBackendDAE;
289 end try;
290 end symbolicJacobianDAE;
291
292 // =============================================================================
293 // section for postOptModule >>generateSymbolicJacobianPast<<
294 //
295 // Symbolic Jacobian subsection
296 // =============================================================================
297
298 protected function generateSymbolicJacobianPast
299 input BackendDAE.BackendDAE inBackendDAE;
300 output BackendDAE.BackendDAE outBackendDAE;
301 protected
302 BackendDAE.EqSystems eqs;
303 BackendDAE.Shared shared;
304 Option<BackendDAE.SymbolicJacobian> symJacA;
305 BackendDAE.SparsePattern sparsePattern;
306 BackendDAE.SparseColoring sparseColoring;
307 BackendDAE.NonlinearPattern nonlinearPattern;
308 AvlTreePathFunction.Tree funcs, functionTree;
309 algorithm
310 6 System.realtimeTick(ClockIndexes.RT_CLOCK_EXECSTAT_JACOBIANS);
311 6 BackendDAE.DAE(eqs=eqs,shared=shared) := inBackendDAE;
312 6 (symJacA, funcs, sparsePattern, sparseColoring, nonlinearPattern) := createSymbolicJacobianforStates(inBackendDAE);
313 6 shared := addBackendDAESharedJacobian(symJacA, sparsePattern, sparseColoring, nonlinearPattern, shared);
314 6 functionTree := BackendDAEUtil.getFunctions(shared);
315 6 functionTree := AvlTreePathFunction.join(functionTree, funcs);
316 6 shared := BackendDAEUtil.setSharedFunctionTree(shared, functionTree);
317 6 outBackendDAE := BackendDAE.DAE(eqs,shared);
318 6 System.realtimeTock(ClockIndexes.RT_CLOCK_EXECSTAT_JACOBIANS);
319 end generateSymbolicJacobianPast;
320
321 protected function createSymbolicJacobianforStates "author: wbraun
322 all functionODE equation are differentiated with respect to the states."
323 input BackendDAE.BackendDAE inBackendDAE;
324 output Option<BackendDAE.SymbolicJacobian> outJacobian;
325 output AvlTreePathFunction.Tree outFunctionTree;
326 output BackendDAE.SparsePattern outSparsePattern;
327 output BackendDAE.SparseColoring outSparseColoring;
328 output BackendDAE.NonlinearPattern outNonlinearPattern;
329 protected
330 BackendDAE.BackendDAE backendDAE2;
331 list<BackendDAE.Var> varlst, knvarlst, states, inputvars, paramvars;
332 BackendDAE.Variables v, globalKnownVars;
333 algorithm
334
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
 6 if Flags.isSet(Flags.JAC_DUMP2) then
335 ✗ print("analytical Jacobians -> start generate system for matrix A time : " + realString(clock()) + "\n");
336 end if;
337 6 backendDAE2 := BackendDAEUtil.copyBackendDAE(inBackendDAE);
338 6 backendDAE2 := BackendDAEOptimize.collapseIndependentContinuousBlocks(backendDAE2);
339 6 backendDAE2 := BackendDAEUtil.transformBackendDAE(backendDAE2,SOME((BackendDAE.NO_INDEX_REDUCTION(),BackendDAE.EXACT())),NONE(),NONE());
340
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
 6 BackendDAE.DAE({BackendDAE.EQSYSTEM(orderedVars = v)},BackendDAE.SHARED(globalKnownVars = globalKnownVars)) := backendDAE2;
341
342 // Prepare all needed variables
343 6 varlst := BackendVariable.varList(v);
344 6 knvarlst := BackendVariable.varList(globalKnownVars);
345 6 states := BackendVariable.getAllStateVarFromVariables(v);
346 6 inputvars := List.select(knvarlst,BackendVariable.isInput);
347 6 paramvars := List.select(knvarlst, BackendVariable.isParam);
348
349
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
 6 if Flags.isSet(Flags.JAC_DUMP2) then
350 ✗ print("analytical Jacobians -> prepared vars for symbolic matrix A time: " + realString(clock()) + "\n");
351 end if;
352
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
 6 if Flags.isSet(Flags.JAC_DUMP2) then
353 ✗ BackendDump.bltdump("System to create symbolic jacobian of: ",backendDAE2);
354 end if;
355 6 (outJacobian, outFunctionTree, outSparsePattern, outSparseColoring, outNonlinearPattern) := generateGenericJacobian(backendDAE2,states,BackendVariable.listVar1(states),BackendVariable.listVar1(inputvars),BackendVariable.listVar1(paramvars),BackendVariable.listVar1(states),varlst,"A",false);
356 end createSymbolicJacobianforStates;
357
358 // =============================================================================
359 // section for postOptModule >>generateSymbolicSensitivities<<
360 //
361 // That function generates symbolic sentivities for parameters
362 // by differentiatiating the states with respect to the parameters
363 // =============================================================================
364
365 public function generateSymbolicSensitivities
366 input BackendDAE.BackendDAE inBackendDAE;
367 output BackendDAE.BackendDAE outBackendDAE;
368 protected
369 BackendDAE.EqSystems eqs;
370 BackendDAE.Shared shared;
371 Option<BackendDAE.SymbolicJacobian> symJacS;
372 BackendDAE.SparsePattern sparsePattern;
373 BackendDAE.SparseColoring sparseColoring;
374 BackendDAE.NonlinearPattern nonlinearPattern;
375 AvlTreePathFunction.Tree funcs, functionTree;
376 algorithm
377 ✗ System.realtimeTick(ClockIndexes.RT_CLOCK_EXECSTAT_JACOBIANS);
378 ✗ BackendDAE.DAE(eqs=eqs,shared=shared) := inBackendDAE;
379 ✗ (symJacS, funcs, sparsePattern, sparseColoring, nonlinearPattern) := createSymbolicJacobianforParameters(inBackendDAE);
380 ✗ shared := addBackendDAESharedJacobian(symJacS, sparsePattern, sparseColoring, nonlinearPattern, shared);
381 ✗ functionTree := BackendDAEUtil.getFunctions(shared);
382 ✗ functionTree := AvlTreePathFunction.join(functionTree, funcs);
383 ✗ shared := BackendDAEUtil.setSharedFunctionTree(shared, functionTree);
384 ✗ outBackendDAE := BackendDAE.DAE(eqs,shared);
385 ✗ System.realtimeTock(ClockIndexes.RT_CLOCK_EXECSTAT_JACOBIANS);
386 end generateSymbolicSensitivities;
387
388 protected function createSymbolicJacobianforParameters
389 "author: wbraun
390 all functionODE equation are differentiated with respect to the parameters."
391 input BackendDAE.BackendDAE inBackendDAE;
392 output Option<BackendDAE.SymbolicJacobian> outJacobian;
393 output AvlTreePathFunction.Tree outFunctionTree;
394 output BackendDAE.SparsePattern outSparsePattern;
395 output BackendDAE.SparseColoring outSparseColoring;
396 output BackendDAE.NonlinearPattern outNonlinearPattern;
397 protected
398 BackendDAE.BackendDAE backendDAE2;
399 list<BackendDAE.Var> varlst, knvarlst, states, inputvars, paramvars;
400 BackendDAE.Variables v, globalKnownVars;
401 algorithm
402 ✗ if Flags.isSet(Flags.JAC_DUMP2) then
403 ✗ print("analytical Jacobians -> start generate system for matrix S time : " + realString(clock()) + "\n");
404 end if;
405
406 ✗ backendDAE2 := BackendDAEUtil.copyBackendDAE(inBackendDAE);
407 ✗ backendDAE2 := BackendDAEOptimize.collapseIndependentContinuousBlocks(backendDAE2);
408 ✗ backendDAE2 := BackendDAEUtil.transformBackendDAE(backendDAE2,SOME((BackendDAE.NO_INDEX_REDUCTION(),BackendDAE.EXACT())),NONE(),NONE());
409 ✗ BackendDAE.DAE({BackendDAE.EQSYSTEM(orderedVars = v)},BackendDAE.SHARED(globalKnownVars = globalKnownVars)) := backendDAE2;
410
411 // Prepare all needed variables
412 ✗ varlst := BackendVariable.varList(v);
413 ✗ knvarlst := BackendVariable.varList(globalKnownVars);
414 ✗ states := BackendVariable.getAllStateVarFromVariables(v);
415 ✗ inputvars := List.select(knvarlst,BackendVariable.isInput);
416 ✗ paramvars := List.select(knvarlst, BackendVariable.isParam);
417
418 ✗ if Flags.isSet(Flags.JAC_DUMP2) then
419 ✗ print("analytical Jacobians -> prepared vars for symbolic matrix S time: " + realString(clock()) + "\n");
420 end if;
421 ✗ if Flags.isSet(Flags.JAC_DUMP2) then
422 ✗ BackendDump.bltdump("System to create symbolic jacobian of: ",backendDAE2);
423 end if;
424 ✗ (outJacobian, outFunctionTree, outSparsePattern, outSparseColoring, outNonlinearPattern) := generateGenericJacobian(backendDAE2,paramvars,BackendVariable.listVar1(states),BackendVariable.listVar1(inputvars),BackendVariable.listVar1(states),BackendVariable.listVar1(states),varlst,"S",false);
425 end createSymbolicJacobianforParameters;
426
427 // =============================================================================
428 // section for postOptModule >>generateSymbolicLinearizationPast<<
429 //
430 // =============================================================================
431
432 public function generateSymbolicLinearizationPast
433 input BackendDAE.BackendDAE inBackendDAE;
434 output BackendDAE.BackendDAE outBackendDAE;
435 algorithm
436 outBackendDAE := matchcontinue inBackendDAE
437 local
438 BackendDAE.EqSystems eqs;
439 BackendDAE.Shared shared;
440 BackendDAE.SymbolicJacobians linearModelMatrices;
441 AvlTreePathFunction.Tree funcs, functionTree;
442 case _ algorithm
443
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 51 times.
 51 true := Flags.getConfigBool(Flags.GENERATE_SYMBOLIC_LINEARIZATION);
444 51 BackendDAE.DAE(eqs=eqs,shared=shared) := inBackendDAE;
445 51 (linearModelMatrices, funcs) := createLinearModelMatrices(inBackendDAE, Config.acceptOptimicaGrammar());
446 51 shared := BackendDAEUtil.setSharedSymJacs(shared, linearModelMatrices);
447 51 functionTree := BackendDAEUtil.getFunctions(shared);
448 51 functionTree := AvlTreePathFunction.join(functionTree, funcs);
449 51 shared := BackendDAEUtil.setSharedFunctionTree(shared, functionTree);
450 51 outBackendDAE := BackendDAE.DAE(eqs,shared);
451 then outBackendDAE;
452
453 else inBackendDAE;
454 end matchcontinue;
455 end generateSymbolicLinearizationPast;
456
457 // =============================================================================
458 // section for postOptModule >>inputDerivativesUsed<<
459 //
460 // check for derivatives of inputs
461 // =============================================================================
462
463 public function inputDerivativesUsed "author: Frenkel TUD 2012-10
464 checks if der(input) is used and report a warning/error."
465 input BackendDAE.BackendDAE inDAE;
466 output BackendDAE.BackendDAE outDAE;
467 algorithm
468 1061 (outDAE, _) := BackendDAEUtil.mapEqSystemAndFold(inDAE, inputDerivativesUsedWork, false);
469 end inputDerivativesUsed;
470
471 protected function inputDerivativesUsedWork "author: Frenkel TUD 2012-10"
472 input BackendDAE.EqSystem isyst;
473 input BackendDAE.Shared inShared;
474 input Boolean inChanged;
475 output BackendDAE.EqSystem osyst;
476 output BackendDAE.Shared outShared = inShared "unused";
477 output Boolean outChanged;
478 protected
479 Boolean hasFailed = false;
480 algorithm
481 (osyst, outChanged) := matchcontinue isyst
482 local
483 BackendDAE.EquationArray orderedEqs;
484 list<DAE.Exp> explst;
485 String s;
486 case BackendDAE.EQSYSTEM(orderedEqs=orderedEqs) algorithm
487
1/2
✓ Branch 3 taken 3919 times.
✗ Branch 4 not taken.
 3919 (_, explst as _::_) := BackendDAEUtil.traverseBackendDAEExpsEqns(orderedEqs, traverserinputDerivativesUsed, (BackendVariable.daeGlobalKnownVars(inShared), {}));
488 ✗ s := stringDelimitList(List.map(explst, ExpressionBasics.printExpStr), "\n");
489 ✗ Error.addMessage(Error.DERIVATIVE_INPUT, {s});
490 ✗ hasFailed := true;
491 ✗ then (BackendDAEUtil.setEqSystEqs(isyst, orderedEqs), true);
492
493 else (isyst, inChanged);
494 end matchcontinue;
495
496 // Fail after error is displayed.
497 // We do it this way, because I was to lazy to rewrite all of this function.
498
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3919 times.
 3919 if hasFailed then fail(); end if;
499 end inputDerivativesUsedWork;
500
501 protected function traverserinputDerivativesUsed "author: Frenkel TUD 2012-10"
502 input DAE.Exp inExp;
503 input tuple<BackendDAE.Variables,list<DAE.Exp>> itpl;
504 output DAE.Exp e;
505 output tuple<BackendDAE.Variables,list<DAE.Exp>> tpl;
506 algorithm
507 111098 (e,tpl) := Expression.traverseExpTopDown(inExp,traverserExpinputDerivativesUsed,itpl);
508 end traverserinputDerivativesUsed;
509
510 protected function traverserExpinputDerivativesUsed
511 input DAE.Exp inExp;
512 input tuple<BackendDAE.Variables,list<DAE.Exp>> tpl;
513 output DAE.Exp outExp;
514 output Boolean cont;
515 output tuple<BackendDAE.Variables,list<DAE.Exp>> outTpl;
516 algorithm
517 (outExp,cont,outTpl) := matchcontinue (inExp,tpl)
518 local
519 BackendDAE.Variables vars;
520 DAE.Exp e;
521 DAE.ComponentRef cr;
522 BackendDAE.Var var;
523 list<DAE.Exp> explst;
524 case (e as DAE.CALL(path=Absyn.IDENT(name = "der"),expLst={DAE.CALL(path=Absyn.IDENT(name = "der"),expLst={DAE.CREF(componentRef=cr)})}),(vars,explst))
525 algorithm
526 ✗ (var,_) := BackendVariable.getVarSingle(cr, vars);
527 ✗ true := BackendVariable.isVarOnTopLevelAndInput(var);
528 ✗ then (e,false,(vars,e::explst));
529 case (e as DAE.CALL(path=Absyn.IDENT(name = "der"),expLst={DAE.CREF(componentRef=cr)}),(vars,explst))
530 algorithm
531 6136 (var,_) := BackendVariable.getVarSingle(cr, vars);
532 ✗ true := BackendVariable.isVarOnTopLevelAndInput(var);
533 ✗ then (e,false,(vars,e::explst));
534 else (inExp,true,tpl);
535 end matchcontinue;
536 end traverserExpinputDerivativesUsed;
537
538 // =============================================================================
539 // solve linear systems with constant jacobian and variable b-Vector
540 //
541 // =============================================================================
542
543 protected function jacobianIsConstant
544 input list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
545 output Boolean isConst;
546 protected
547 list<BackendDAE.Equation> eqs;
548 algorithm
549 589 eqs := List.map(jac, Util.tuple33);
550 589 isConst := not List.any(eqs, variableResidual);
551 end jacobianIsConstant;
552
553 protected function variableResidual
554 input BackendDAE.Equation eq;
555 output Boolean isNotConst;
556 algorithm
557 isNotConst := match eq
558 case BackendDAE.RESIDUAL_EQUATION(exp=DAE.RCONST(_))
559 then false;
560
561 else true;
562 end match;
563 end variableResidual;
564
565 protected function replaceStrongComponent "replaces the indexed component with compsNew and adds compsAdd at the end. the assignments will be updated"
566 input BackendDAE.EqSystem systIn;
567 input Integer idx;
568 input BackendDAE.StrongComponents compsNew;
569 input BackendDAE.StrongComponents compsAdd;
570 output BackendDAE.EqSystem systOut = systIn;
571 protected
572 BackendDAE.Matching matching;
573 array<Integer> ass1, ass2, assAdd;
574 BackendDAE.StrongComponents comps;
575 algorithm
576 ✗ BackendDAE.EQSYSTEM(matching=BackendDAE.MATCHING(ass1=ass1, ass2=ass2, comps=comps)) := systIn;
577 ✗ if not listEmpty(compsAdd) then
578 ✗ assAdd := arrayCreate(listLength(compsAdd), 0);
579 ✗ ass1 := arrayAppend(ass1, assAdd);
580 ✗ ass2 := arrayAppend(ass2, assAdd);
581 ✗ List.map2_0(compsAdd, updateAssignment, ass1, ass2);
582 end if;
583 ✗ List.map2_0(compsNew, updateAssignment, ass1, ass2);
584 ✗ comps := List.replaceAtWithList(compsNew, idx, comps);
585 ✗ systOut.matching := BackendDAE.MATCHING(ass1, ass2, listAppend(comps, compsAdd));
586 ✗ systOut := BackendDAEUtil.setEqSystMatrices(systOut);
587 end replaceStrongComponent;
588
589 protected function updateAssignment
590 input BackendDAE.StrongComponent comp;
591 input array<Integer> ass1;
592 input array<Integer> ass2;
593 algorithm
594 () := matchcontinue comp
595 local
596 Integer eq,var;
597 case BackendDAE.SINGLEEQUATION(eqn=eq,var=var)
598 algorithm
599 ✗ arrayUpdate(ass2,eq,var);
600 ✗ arrayUpdate(ass1,var,eq);
601 then ();
602 else
603 then ();
604 end matchcontinue;
605 end updateAssignment;
606
607 protected function solveConstJacLinearSystem
608 input BackendDAE.EqSystem syst;
609 input BackendDAE.Shared ishared;
610 input list<BackendDAE.Equation> eqn_lst;
611 input list<Integer> eqn_indxs;
612 input list<BackendDAE.Var> var_lst;
613 input list<Integer> var_indxs;
614 input list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
615 input Integer sysIdxIn;
616 input Integer compIdxIn;
617 output list<BackendDAE.Equation> sysEqsOut;
618 output list<BackendDAE.Equation> bEqsOut;
619 output list<BackendDAE.Var> bVarsOut;
620 output array<Integer> orderOut;
621 output Integer sysIdxOut;
622 protected
623 BackendDAE.Variables vars,v;
624 BackendDAE.EquationArray eqns,eqns1;
625 list<DAE.Exp> beqs;
626 list<DAE.ElementSource> sources;
627 BackendDAE.Matching matching;
628 AvlTreePathFunction.Tree funcs;
629 BackendDAE.StateSets stateSets;
630 BackendDAE.BaseClockPartitionKind partitionKind;
631
632 array<array<Real>> A;
633 array<Real> b;
634 Integer row,n;
635 array<Integer> order;
636 algorithm
637 BackendDAE.EQSYSTEM(orderedVars=vars,orderedEqs=eqns,matching=matching,stateSets=stateSets,partitionKind=partitionKind) := syst;
638 ✗ BackendDAE.SHARED(functionTree=funcs) := ishared;
639 ✗ eqns1 := BackendEquation.listEquation(eqn_lst);
640 ✗ v := BackendVariable.listVar1(var_lst);
641 ✗ n := listLength(var_lst);
642 ✗ (beqs,sources) := BackendDAEUtil.getEqnSysRhs(eqns1,v,SOME(funcs));
643 ✗ beqs := listReverse(beqs);
644 //print("bside: \n"+ExpressionDump.printExpListStr(beqs)+"\n");
645 ✗ A := evaluateConstantJacobianArray(listLength(var_lst),jac);
646 //print("JacVals\n"+stringDelimitList(List.map(jacVals,rListStr),"\n")+"\n\n");
647
648 ✗ b := arrayCreate(n*n,0.0); // i.e. a matrix for the b-vars to get their coefficients independently [(b1,0,0);(0,b2,0),(0,0,b3)]
649 ✗ order := arrayCreate(n,0);
650 ✗ for row in 1:n loop
651 ✗ arrayUpdate(b,(row-1)*n+row,1.0);
652 end for;
653 //print("b\n"+stringDelimitList(List.mapArray(b,realString),", ")+"\n\n");
654 //print("A\n"+stringDelimitList(List.mapArray(A,realString),", ")+"\n\n");
655 ✗ gauss(A,b,1,n,List.intRange(n),order);
656 //print("the order: "+stringDelimitList(List.mapArray(order,intString),",")+"\n");
657
658 ✗ (bVarsOut,bEqsOut) := createBVecVars(sysIdxIn,compIdxIn,n,DAE.T_REAL_DEFAULT,beqs);
659 ✗ sysEqsOut := createSysEquations(A,b,n,order,var_lst,bVarsOut);
660 ✗ for a in A loop
661 ✗ GCExt.free(a);
662 end for;
663 ✗ GCExt.free(A);
664 ✗ GCExt.free(b);
665 ✗ sysIdxOut := sysIdxIn+1;
666 orderOut := order;
667 end solveConstJacLinearSystem;
668
669 protected function createSysEquations "creates new equations for a linear system with constant Jacobian matrix.
670 author: Waurich TUD 2015-03"
671 input array<array<Real>> A;
672 input array<Real> b;
673 input Integer n;
674 input array<Integer> order;
675 input list<BackendDAE.Var> xVars;
676 input list<BackendDAE.Var> bVars;
677 output list<BackendDAE.Equation> sysEqs = {};
678 protected
679 Integer i;
680 Integer row;
681 DAE.Exp lhs, rhs;
682 list<DAE.Exp> coeffExps, xExps, bExps, xProds, bProds;
683 list<Real> coeffs;
684 BackendDAE.Equation eq;
685 algorithm
686 ✗ xExps := List.map(xVars, BackendVariable.varExp2);
687 ✗ bExps := List.map(bVars, BackendVariable.varExp2);
688 ✗ for i in 1:n loop
689 ✗ row := arrayGet(order,i);
690 ✗ coeffs := arrayList(A[row]);
691 ✗ coeffExps := List.map(coeffs,Expression.makeRealExp);
692 ✗ xProds := List.threadMap1(coeffExps,xExps,makeBinaryExp,DAE.MUL(DAE.T_REAL_DEFAULT));
693 ✗ lhs := List.fold1(xProds,Expression.makeBinaryExp,DAE.ADD(DAE.T_REAL_DEFAULT),DAE.RCONST(0.0));
694 ✗ (lhs,_) := ExpressionSimplify.simplify(lhs);
695 ✗ coeffs := Array.getRange((row-1)*n+1,(row*n),b);
696 ✗ coeffExps := List.map(coeffs,Expression.makeRealExp);
697 ✗ bProds := List.threadMap1(coeffExps,bExps,makeBinaryExp,DAE.MUL(DAE.T_REAL_DEFAULT));
698 ✗ rhs := List.fold1(bProds,Expression.makeBinaryExp,DAE.ADD(DAE.T_REAL_DEFAULT),DAE.RCONST(0.0));
699 ✗ (rhs,_) := ExpressionSimplify.simplify(rhs);
700 ✗ eq := BackendDAE.EQUATION(lhs,rhs,DAE.emptyElementSource,BackendDAE.EQ_ATTR_DEFAULT_DYNAMIC);
701 sysEqs := eq::sysEqs;
702 end for;
703 end createSysEquations;
704
705 public function makeBinaryExp
706 input DAE.Exp inLhs;
707 input DAE.Exp inRhs;
708 input DAE.Operator inOp;
709 output DAE.Exp outExp;
710 algorithm
711 ✗ outExp := DAE.BINARY(inLhs, inOp, inRhs);
712 end makeBinaryExp;
713
714 protected function createBVecVars "creates variables for the b-Vector of a linear system with constant Jacobian
715 author:Waurich TUD 2015-03"
716 input Integer sysIdx;
717 input Integer compIdx;
718 input Integer size;
719 input DAE.Type typ;
720 input list<DAE.Exp> bExps;
721 output list<BackendDAE.Var> varLst = {};
722 output list<BackendDAE.Equation> eqLst = {};
723 protected
724 String ident;
725 Integer i;
726 DAE.ComponentRef cref;
727 BackendDAE.Var var;
728 BackendDAE.Equation beq;
729 algorithm
730 ✗ for i in 1:size loop
731 ✗ ident := "$sys"+intString(sysIdx)+"_"+intString(compIdx)+"_b"+intString(i);
732 ✗ cref := ComponentReferenceBasics.makeCrefIdent(ident,typ,{});
733 ✗ var := BackendVariable.makeVar(cref);
734 varLst := var::varLst;
735 ✗ beq := BackendDAE.EQUATION(listGet(bExps,i),Expression.crefExp(cref),DAE.emptyElementSource,BackendDAE.EQ_ATTR_DEFAULT_DYNAMIC);
736 eqLst := beq::eqLst;
737 end for;
738 end createBVecVars;
739
740 protected function gauss
741 input array<array<Real>> A;
742 input array<Real> b;
743 input Integer indxIn;
744 input Integer n;
745 input list<Integer> rangeIn;
746 input array<Integer> permutation;
747 protected
748 Integer pivotIdx,pos, ir, ic;// ir=rowIdx, ic=columnIdx, p_ir=permuted row idx
749 Real pivot, entry, b_entry, first;
750 list<Integer> range;
751 algorithm
752 () := matchcontinue permutation
753 case _
754 algorithm
755 ✗ true := intLe(indxIn,n);
756 ✗ (pivotIdx,pivot) := getPivotElement(A,rangeIn,indxIn,n);
757 //print("pivot: "+intString(pivotIdx)+" has value: "+realString(pivot)+"\n");
758 ✗ arrayUpdate(permutation,indxIn,pivotIdx);
759 ✗ range := List.deleteMemberOnTrue(pivotIdx,rangeIn,intEq);
760
761 // the pivot row in the A-matrix divided by the pivot element
762 ✗ for ic in indxIn:n loop
763 ✗ entry := arrayGet(A[pivotIdx],ic);
764 ✗ entry := realDiv(entry,pivot); //divide column entry with pivot element
765 //print(" pos "+intString(pos)+" entry "+realString(arrayGet(A,pos))+"\n");
766 ✗ arrayUpdate(A[pivotIdx],ic,entry);
767 end for;
768 // the complete pivot row of the b-vector divided by the pivot element
769 ✗ for ic in 1:n loop
770 ✗ pos := (pivotIdx-1)*n+ic;
771 ✗ b_entry := arrayGet(b,pos);
772 ✗ b_entry := realDiv(b_entry,pivot);
773 ✗ arrayUpdate(b,pos,b_entry);
774 end for;
775
776 // the remaining rows
777 ✗ for ir in range loop
778 ✗ first := arrayGet(A[ir],indxIn); //the first row element, that is going to be zero
779 //print("first "+realString(first)+"\n");
780 ✗ for ic in indxIn:n loop
781 ✗ pos := (ir-1)*n+ic;
782 ✗ entry := arrayGet(A[ir],ic); // the current entry
783 ✗ pivot := arrayGet(A[pivotIdx],ic); // the element from the column in the pivot row
784 //print("pivot "+realString(pivot)+"\n");
785 //print("ir "+intString(ir)+" pos "+intString(pos)+" entry0 "+realString(entry)+" entry1 "+realString(realSub(entry,realDiv(first,pivot)))+"\n");
786 ✗ entry := realSub(entry,realMul(first,pivot));
787 ✗ arrayUpdate(A[ir],ic,entry);
788 ✗ b_entry := arrayGet(b,pos);
789 ✗ pivot := arrayGet(b,(pivotIdx-1)*n+ic);
790 ✗ b_entry := b_entry - realMul(first,pivot);
791 ✗ arrayUpdate(b,pos,b_entry);
792 end for;
793 end for;
794 //print("A\n"+stringDelimitList(List.mapArray(A, realString),", ")+"\n\n");
795 //print("b\n"+stringDelimitList(List.mapArray(b, realString),", ")+"\n\n");
796
797 //print("new permutation: "+stringDelimitList(List.mapArray(permutation, intString),",")+"\n");
798 //print("JACB "+intString(indxIn)+" \n"+stringDelimitList(List.mapArray(jacB, rListStr),"\n ")+"\n\n");
799 ✗ gauss(A,b,indxIn+1,n,range,permutation);
800 then();
801 else ();
802 end matchcontinue;
803 end gauss;
804
805 protected function getPivotElement "gets the highest element in the startIdx'th to n'th rows and the startidx'th column"
806 input array<array<Real>> A;
807 input list<Integer> rangeIn;
808 input Integer startIdx;
809 input Integer n;
810 output Integer pos = 0;
811 output Real value = 0.0;
812 protected
813 Integer i;
814 Real entry;
815 algorithm
816 ✗ for i in rangeIn loop
817 ✗ entry := arrayGet(A[i],startIdx);
818 //print("i "+intString(i)+" pi "+intString(p_i)+" entry "+realString(entry)+"\n");
819 ✗ if realAbs(entry) > value then
820 value := entry;
821 pos := i;
822 end if;
823 end for;
824 end getPivotElement;
825
826 protected function rListStr
827 input list<Real> l;
828 output String s;
829 algorithm
830 ✗ s := stringDelimitList(List.map(l,realString)," , ");
831 end rListStr;
832
833
834
835 // =============================================================================
836 // unsorted section
837 //
838 // =============================================================================
839
840 protected function constantLinearSystem0
841 input BackendDAE.EqSystem isyst;
842 input BackendDAE.Shared inShared;
843 input tuple<Boolean, Integer> iTpl "<inChanged,sysIdxIn>";
844 output BackendDAE.EqSystem osyst;
845 output BackendDAE.Shared outShared;
846 output tuple<Boolean,Integer> oTpl "<oChanged,sysIdxOut>";
847 protected
848 Boolean changed;
849 Integer sysIdx;
850 BackendDAE.StrongComponents comps;
851 algorithm
852 3309 (changed,sysIdx) := iTpl;
853
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3309 times.
 3309 BackendDAE.EQSYSTEM(matching=BackendDAE.MATCHING(comps=comps)) := isyst;
854 3309 (osyst, outShared, changed, sysIdx) := constantLinearSystem1(isyst, inShared, comps, changed, sysIdx, 1);
855 3309 osyst := constantLinearSystem2(changed, osyst);
856
2/2
✓ Branch 0 taken 3304 times.
✓ Branch 1 taken 5 times.
 6613 oTpl := (changed,sysIdx+1);
857 end constantLinearSystem0;
858
859 protected function constantLinearSystem2
860 input Boolean b;
861 input BackendDAE.EqSystem isyst;
862 output BackendDAE.EqSystem osyst;
863 algorithm
864 osyst := match(b,isyst)
865 local
866 BackendDAE.Variables vars;
867 BackendDAE.EquationArray eqns;
868 BackendDAE.StateSets stateSets;
869 BackendDAE.BaseClockPartitionKind partitionKind;
870
871 case (false,_) then isyst;
872 // case (true,BackendDAE.EQSYSTEM(orderedVars=vars,orderedEqs=eqns,matching=BackendDAE.NO_MATCHING()))
873 case (true,BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns, stateSets=stateSets, partitionKind=partitionKind))
874 algorithm
875 // remove empty entries from vars/eqns
876 5 vars := BackendVariable.listVar1(BackendVariable.varList(vars));
877 5 eqns := BackendEquation.listEquation(BackendEquation.equationList(eqns));
878 5 then
879 BackendDAEUtil.createEqSystem(vars, eqns, stateSets, partitionKind);
880 /* case (true,BackendDAE.EQSYSTEM(orderedVars=vars,orderedEqs=eqns,matching=BackendDAE.MATCHING(ass1=ass1,ass2=ass2,comps=comps)))
881 then
882 updateEquationSystemMatching(vars,eqns,ass1,ass2,comps);
883 */ end match;
884 end constantLinearSystem2;
885
886 protected function constantLinearSystem1
887 input BackendDAE.EqSystem isyst;
888 input BackendDAE.Shared ishared;
889 input BackendDAE.StrongComponents inComps;
890 input Boolean inRunMatching;
891 input Integer sysIdxIn;
892 input Integer compIdxIn;
893 output BackendDAE.EqSystem osyst = isyst;
894 output BackendDAE.Shared oshared = ishared;
895 output Boolean runMatching = inRunMatching;
896 output Integer sysIdxOut = sysIdxIn;
897 protected
898 Integer compIdx = compIdxIn;
899 Boolean b;
900 algorithm
901
2/2
✓ Branch 0 taken 56930 times.
✓ Branch 1 taken 3309 times.
 60239 for comp in inComps loop
902 56930 (osyst, oshared, b, sysIdxOut, compIdx) := constantLinearSystemWork(osyst, oshared, comp, sysIdxOut, compIdx);
903 56930 runMatching := b or runMatching;
904 end for;
905 end constantLinearSystem1;
906
907 protected function constantLinearSystemWork
908 input BackendDAE.EqSystem isyst;
909 input BackendDAE.Shared ishared;
910 input BackendDAE.StrongComponent comp;
911 input Integer sysIdxIn;
912 input Integer compIdxIn;
913 output BackendDAE.EqSystem osyst;
914 output BackendDAE.Shared oshared;
915 output Boolean outRunMatching;
916 output Integer sysIdxOut;
917 output Integer compIdxOut;
918 algorithm
919 (osyst, oshared, outRunMatching, sysIdxOut, compIdxOut):=
920 matchcontinue (isyst, ishared, comp)
921 local
922 BackendDAE.Variables vars;
923 BackendDAE.EquationArray eqns;
924 list<BackendDAE.Equation> eqn_lst;
925 list<BackendDAE.Var> var_lst;
926 list<Integer> eindex,vindx;
927 list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
928 BackendDAE.EqSystem syst;
929 BackendDAE.Shared shared;
930
931 Integer sysIdx;
932 array<Integer> order;
933 list<Integer> bVarIdcs,bEqIdcs;
934 list<BackendDAE.Var> bVars;
935 list<BackendDAE.Equation> bEqs,sysEqs;
936 BackendDAE.StrongComponents bComps,sysComps;
937
938 case (syst, shared, (BackendDAE.EQUATIONSYSTEM( eqns=eindex, vars=vindx, jac=BackendDAE.FULL_JACOBIAN(SOME(jac)),
939 jacType=BackendDAE.JAC_CONSTANT() )))
940 algorithm
941 //the A-matrix and the b-Vector are constant
942 7 eqn_lst := BackendEquation.getList(eindex, syst.orderedEqs);
943 7 var_lst := List.map1r(vindx, BackendVariable.getVarAt, syst.orderedVars);
944 7 (syst,shared) := solveLinearSystem(syst, shared, eqn_lst, eindex, var_lst, vindx, jac);
945 7 then (syst,shared,true,sysIdxIn,compIdxIn+1);
946
947 case ( syst as BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns), shared,
948 BackendDAE.EQUATIONSYSTEM( eqns=eindex, vars=vindx, jac=BackendDAE.FULL_JACOBIAN(SOME(jac)),
949 jacType=BackendDAE.JAC_LINEAR() ) )
950 algorithm
951
2/2
✓ Branch 1 taken 29 times.
✓ Branch 2 taken 589 times.
 618 true := BackendDAEUtil.isSimulationDAE(ishared);
952 //only the A-matrix is constant, apply Gaussian Elimination
953 589 eqn_lst := BackendEquation.getList(eindex, eqns);
954 589 var_lst := List.map1r(vindx, BackendVariable.getVarAt, vars);
955
2/2
✓ Branch 1 taken 540 times.
✓ Branch 2 taken 49 times.
 589 true := jacobianIsConstant(jac);
956
1/2
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
 49 true := Flags.isSet(Flags.CONSTJAC);
957 //true = intEq(compIdxIn,37) and intEq(sysIdxIn,1);
958 //print("ITS CONSTANT\n");
959 //print("THE COMPIDX: "+intString(compIdxIn)+" THE SYSIDX"+intString(sysIdxIn)+"\n");
960 //BackendDump.dumpEqnsSolved2({comp},eqns,vars);
961 ✗ eqn_lst := BackendEquation.getList(eindex,eqns);
962 ✗ var_lst := List.map1r(vindx, BackendVariable.getVarAt, vars);
963 ✗ (sysEqs, bEqs, bVars, order, sysIdx) :=
964 solveConstJacLinearSystem(syst, shared, eqn_lst, eindex, listReverse(var_lst), vindx, jac, sysIdxIn, compIdxIn);
965 //print("the b-vector stuff \n");
966 //BackendDump.printEquationList(bEqs);
967 //BackendDump.printVarList(bVars);
968 //print("the sysEqs stuff \n");
969 //BackendDump.printEquationList(sysEqs);
970 //build comps
971 //print("size"+intString(BackendEquation.equationArraySize(eqns))+"\n");
972 //print("numberOfElement"+intString(BackendEquation.getNumberOfEquations(eqns))+"\n");
973 //print("arrSize"+intString(BackendDAEUtil.equationArraySize2(eqns))+"\n");
974 //print("length"+intString(listLength(BackendEquation.equationList(eqns)))+"\n");
975 ✗ bVarIdcs := List.intRange2(BackendVariable.varsSize(vars)+1, BackendVariable.varsSize(vars)+listLength(bVars));
976 ✗ bEqIdcs := List.intRange2(BackendEquation.getNumberOfEquations(eqns)+1, BackendEquation.getNumberOfEquations(eqns)+listLength(bEqs));
977 ✗ bComps := List.threadMap(bEqIdcs, bVarIdcs, BackendDAEUtil.makeSingleEquationComp);
978 ✗ sysComps := List.threadMap( List.map1(arrayList(order), List.getIndexFirst, eindex), listReverse(vindx),
979 BackendDAEUtil.makeSingleEquationComp );
980 //print("bCOMPS\n");
981 //BackendDump.dumpComponents(bComps);
982 //print("SYSCOMPS\n");
983 //BackendDump.dumpComponents(sysComps);
984 //build system
985 ✗ syst.orderedVars := List.fold(bVars, BackendVariable.addVar, vars);
986 ✗ eqns := BackendEquation.addList(bEqs, eqns);
987 ✗ syst.orderedEqs := List.threadFold(eindex, sysEqs, BackendEquation.setAtIndexFirst, eqns);
988 ✗ syst := BackendDAEUtil.setEqSystMatrices(syst);
989 ✗ syst := replaceStrongComponent(syst,compIdxIn,sysComps,bComps);
990 //print("compIdxIn"+intString(compIdxIn)+"\n");
991 ✗ then (syst, ishared, false, sysIdx, compIdxIn+listLength(sysComps));
992 56923 else (isyst, ishared, false, sysIdxIn, compIdxIn+1);
993 end matchcontinue;
994 end constantLinearSystemWork;
995
996 protected function solveLinearSystem
997 input BackendDAE.EqSystem inSyst;
998 input BackendDAE.Shared ishared;
999 input list<BackendDAE.Equation> eqn_lst;
1000 input list<Integer> eqn_indxs;
1001 input list<BackendDAE.Var> var_lst;
1002 input list<Integer> var_indxs;
1003 input list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
1004 output BackendDAE.EqSystem osyst;
1005 output BackendDAE.Shared oshared;
1006 algorithm
1007 (osyst, oshared) := match (inSyst, ishared)
1008 local
1009 BackendDAE.Variables v;
1010 BackendDAE.EquationArray eqns, eqns1;
1011 list<DAE.Exp> beqs;
1012 list<DAE.ElementSource> sources;
1013 list<Real> rhsVals,solvedVals;
1014 list<list<Real>> jacVals;
1015 Integer linInfo;
1016 list<DAE.ComponentRef> names;
1017 AvlTreePathFunction.Tree funcs;
1018 BackendDAE.Shared shared;
1019 BackendDAE.EqSystem syst;
1020
1021 case (syst as BackendDAE.EQSYSTEM(), BackendDAE.SHARED(functionTree=funcs))
1022 algorithm
1023 7 eqns1 := BackendEquation.listEquation(eqn_lst);
1024 7 v := BackendVariable.listVar1(var_lst);
1025 7 (beqs, sources) := BackendDAEUtil.getEqnSysRhs(eqns1, v, SOME(funcs));
1026 7 beqs := listReverse(beqs);
1027 7 rhsVals := ValuesUtil.valueReals(List.map(beqs, Ceval.cevalSimple));
1028 7 jacVals := evaluateConstantJacobian(listLength(var_lst), jac);
1029 7 (solvedVals, linInfo) := System.dgesv(jacVals, rhsVals);
1030 7 names := List.map(var_lst, BackendVariable.varCref);
1031 7 checkLinearSystem(linInfo, names, jacVals, rhsVals, eqn_lst);
1032 7 sources := List.map1( sources, ElementSource.addSymbolicTransformation,
1033 DAE.LINEAR_SOLVED(names, jacVals, rhsVals, solvedVals) );
1034 7 (v, eqns, shared) := changeConstantLinearSystemVars( var_lst, solvedVals, sources, var_indxs,
1035 syst.orderedVars, syst.orderedEqs, ishared );
1036 7 syst.orderedVars := v;
1037 7 syst.orderedEqs := List.fold(eqn_indxs, BackendEquation.delete, eqns);
1038
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 then
1039 (BackendDAEUtil.setEqSystMatrices(syst), shared);
1040 end match;
1041 end solveLinearSystem;
1042
1043 protected function changeConstantLinearSystemVars
1044 input list<BackendDAE.Var> inVarLst;
1045 input list<Real> inSolvedVals;
1046 input list<DAE.ElementSource> inSources;
1047 input list<Integer> var_indxs;
1048 input BackendDAE.Variables inVars;
1049 input BackendDAE.EquationArray ieqns;
1050 input BackendDAE.Shared ishared;
1051 output BackendDAE.Variables outVars;
1052 output BackendDAE.EquationArray oeqns;
1053 output BackendDAE.Shared oshared;
1054 algorithm
1055 (outVars,oeqns,oshared) := match (inVarLst, inSolvedVals, inSources, var_indxs, inVars, ieqns)
1056 local
1057 BackendDAE.Var v,v1;
1058 list<BackendDAE.Var> varlst;
1059 list<DAE.ElementSource> slst;
1060 BackendDAE.Variables vars,vars1,vars2;
1061 Real r;
1062 list<Real> rlst;
1063 BackendDAE.Shared shared;
1064 BackendDAE.EquationArray eqns;
1065 Integer indx;
1066 list<Integer> vindxs;
1067 DAE.ComponentRef cref;
1068 DAE.Type tp;
1069 DAE.Exp e;
1070 case ({}, {}, {}, _, vars, eqns) then (vars,eqns,ishared);
1071 case ((BackendDAE.VAR(varName=cref,varKind=BackendDAE.STATE(),varType=tp))::varlst, r::rlst, _::slst, _::vindxs, vars, eqns)
1072 algorithm
1073 1 e := Expression.makeCrefExp(cref, tp);
1074 1 e := Expression.expDer(e);
1075 1 eqns := BackendEquation.add(BackendDAE.EQUATION(e, DAE.RCONST(r), DAE.emptyElementSource, BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN), eqns);
1076 1 (vars2,eqns,shared) := changeConstantLinearSystemVars(varlst,rlst,slst,vindxs,vars,eqns,ishared);
1077 then (vars2,eqns,shared);
1078 case (v::varlst, r::rlst, _::slst, indx::vindxs, vars, eqns)
1079 algorithm
1080 32 v1 := BackendVariable.setBindExp(v, SOME(DAE.RCONST(r)));
1081 16 v1 := BackendVariable.setVarStartValue(v1,DAE.RCONST(r));
1082 // ToDo: merge source of var and equation
1083 16 (vars1,_) := BackendVariable.removeVar(indx, vars);
1084 16 shared := BackendVariable.addGlobalKnownVarDAE(v1,ishared);
1085 16 (vars2,eqns,shared) := changeConstantLinearSystemVars(varlst,rlst,slst,vindxs,vars1,eqns,shared);
1086 then (vars2,eqns,shared);
1087 end match;
1088 end changeConstantLinearSystemVars;
1089
1090 public function evaluateConstantJacobian
1091 "Evaluate a constant Jacobian so we can solve a linear system during runtime"
1092 input Integer size;
1093 input list<tuple<Integer,Integer,BackendDAE.Equation>> jac;
1094 output list<list<Real>> vals;
1095 protected
1096 array<array<Real>> valarr;
1097 list<array<Real>> tmp2;
1098 algorithm
1099 269 valarr := evaluateConstantJacobianArray(size, jac);
1100 269 tmp2 := arrayList(valarr);
1101 269 vals := List.map(tmp2,arrayList);
1102 end evaluateConstantJacobian;
1103
1104 protected function evaluateConstantJacobianArray
1105 "Evaluate a constant Jacobian so we can solve a linear system during runtime"
1106 input Integer size;
1107 input list<tuple<Integer,Integer,BackendDAE.Equation>> jac;
1108 output array<array<Real>> valarr;
1109 protected
1110 array<Real> tmp;
1111 list<array<Real>> tmp2;
1112 algorithm
1113 269 tmp := arrayCreate(size,0.0);
1114 269 tmp2 := List.map(List.fill(tmp,size),arrayCopy);
1115 269 valarr := listArray(tmp2);
1116 269 List.map1_0(jac,evaluateConstantJacobian2,valarr);
1117 end evaluateConstantJacobianArray;
1118
1119 protected function evaluateConstantJacobian2
1120 input tuple<Integer,Integer,BackendDAE.Equation> jac;
1121 input array<array<Real>> vals;
1122 algorithm
1123 () := match jac
1124 local
1125 DAE.Exp exp;
1126 Integer i1,i2;
1127 Real r;
1128 case (i1,i2,BackendDAE.RESIDUAL_EQUATION(exp=exp))
1129 algorithm
1130
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4396 times.
 4396 Values.REAL(r) := Ceval.cevalSimple(exp);
1131 4396 arrayUpdate(arrayGet(vals,i1),i2,r);
1132 then ();
1133 end match;
1134 end evaluateConstantJacobian2;
1135
1136 protected function checkLinearSystem
1137 input Integer info;
1138 input list<DAE.ComponentRef> vars;
1139 input list<list<Real>> jac;
1140 input list<Real> rhs;
1141 input list<BackendDAE.Equation> eqnlst;
1142 algorithm
1143 () := matchcontinue info
1144 local
1145 String infoStr,syst,varnames,varname,rhsStr,jacStr,eqnstr;
1146 case 0 then ();
1147 case _
1148 algorithm
1149 ✗ true := info > 0;
1150 ✗ varname := ComponentReferenceBasics.printComponentRefStr(listGet(vars,info));
1151 ✗ infoStr := intString(info);
1152 ✗ varnames := stringDelimitList(List.map(vars,ComponentReferenceBasics.printComponentRefStr)," ;\n ");
1153 ✗ rhsStr := stringDelimitList(List.map(rhs, realString)," ;\n ");
1154 ✗ jacStr := stringDelimitList(List.map1(List.mapList(jac,realString),stringDelimitList," , ")," ;\n ");
1155 ✗ eqnstr := BackendDump.dumpEqnsStr(eqnlst);
1156 ✗ syst := stringAppendList({"\n",eqnstr,"\n[\n ", jacStr, "\n]\n *\n[\n ",varnames,"\n]\n =\n[\n ",rhsStr,"\n]"});
1157 ✗ Error.addMessage(Error.LINEAR_SYSTEM_SINGULAR, {syst,infoStr,varname});
1158 ✗ then fail();
1159 case _
1160 algorithm
1161 ✗ true := info < 0;
1162 ✗ varnames := stringDelimitList(List.map(vars,ComponentReferenceBasics.printComponentRefStr)," ;\n ");
1163 ✗ rhsStr := stringDelimitList(List.map(rhs, realString)," ; ");
1164 ✗ jacStr := stringDelimitList(List.map1(List.mapList(jac,realString),stringDelimitList," , ")," ; ");
1165 ✗ eqnstr := BackendDump.dumpEqnsStr(eqnlst);
1166 ✗ syst := stringAppendList({eqnstr,"\n[", jacStr, "] * [",varnames,"] = [",rhsStr,"]"});
1167 ✗ Error.addMessage(Error.LINEAR_SYSTEM_INVALID, {"LAPACK/dgesv",syst});
1168 ✗ then fail();
1169 end matchcontinue;
1170 end checkLinearSystem;
1171
1172 public function generateSparsePattern "author: wbraun
1173 Function generated for a given set of variables and
1174 equations the sparsity pattern and a coloring of Jacobian matrix A^(NxM).
1175 col: N = size(diffVars)
1176 rows : M = size(diffedVars)
1177 The sparsity pattern is represented basically as a list of lists, every list
1178 represents the non-zero elements of a row.
1179
1180 The coloring is saved as a list of lists, every list contains the
1181 cols with the same color."
1182 input BackendDAE.BackendDAE inBackendDAE;
1183 input list<BackendDAE.Var> inIndependentVars "vars";
1184 input list<BackendDAE.Var> inDependentVars "eqns";
1185 input Boolean nonlinearPattern = false;
1186 input Boolean withColoring = true "false gives one colour per column";
1187 output BackendDAE.SparsePattern outSparsePattern;
1188 output BackendDAE.SparseColoring outColoredCols;
1189 protected
1190 constant Boolean debug = false;
1191 String patternName = if nonlinearPattern then "Nonlinear" else "Sparsity";
1192 algorithm
1193 (outSparsePattern,outColoredCols) := matchcontinue(inBackendDAE,inIndependentVars,inDependentVars)
1194 local
1195 BackendDAE.Shared shared;
1196 BackendDAE.EqSystem syst, syst1;
1197 BackendDAE.StrongComponents comps;
1198 BackendDAE.AdjacencyMatrix adjMatrix, adjMatrixT;
1199 BackendDAE.Matching bdaeMatching;
1200
1201
1202 Integer sizeN, sizeM, adjSize, adjSizeT;
1203 Integer nonZeroElements;
1204 list<Integer> nodesEqnsIndex;
1205 list<list<Integer>> sparsepattern,sparsepatternT;
1206 list<BackendDAE.Var> jacDiffVars, dependentVars, independentVars;
1207 BackendDAE.Variables varswithDiffs;
1208 BackendDAE.EquationArray orderedEqns;
1209 array<Integer> ass1;
1210 array<list<Integer>> coloredArray;
1211
1212 list<DAE.ComponentRef> depCompRefsLst, inDepCompRefsLst;
1213 array<DAE.ComponentRef> depCompRefs, inDepCompRefs;
1214
1215 array<list<Integer>> eqnSparse, varSparse, sparseArray, sparseArrayT;
1216 array<Integer> mark, usedvar;
1217
1218 BackendDAE.SparseColoring coloring;
1219 list<list<DAE.ComponentRef>> translated;
1220 list<tuple<DAE.ComponentRef,list<DAE.ComponentRef>>> sparsetuple, sparsetupleT;
1221
1222 // if there are no independent var, no pattern needed, otherwise there
1223 // is an empty pattern for the dependent variables
1224 case (_,_,{}) then (({},{},({},{}),-1),{});
1225 case(BackendDAE.DAE(eqs = (syst as BackendDAE.EQSYSTEM(matching=bdaeMatching as BackendDAE.MATCHING(comps=comps, ass1=ass1)))::{}),independentVars,dependentVars)
1226 algorithm
1227
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4718 times.
 4718 if Flags.isSet(Flags.DUMP_SPARSE_VERBOSE) then
1228 ✗ print(" start getting " + patternName + " pattern for variables : " + intString(listLength(dependentVars)) + " and the independent vars: " + intString(listLength(independentVars)) +"\n");
1229 end if;
1230 if debug then execStat("generateSparsePattern -> do start "); end if;
1231 // prepare crefs
1232 4718 depCompRefsLst := List.map(dependentVars, BackendVariable.varCref);
1233 4718 depCompRefs := listArray(depCompRefsLst);
1234 sizeM := arrayLength(depCompRefs);
1235
1236 // create jacobian vars
1237 4718 (jacDiffVars,inDepCompRefsLst) := createInDepVars(independentVars);
1238 4718 inDepCompRefs := listArray(inDepCompRefsLst);
1239 sizeN := arrayLength(inDepCompRefs);
1240
1241 // generate adjacency matrix including diff vars
1242 4718 syst1 as BackendDAE.EQSYSTEM(orderedVars=varswithDiffs,orderedEqs=orderedEqns) := BackendDAEUtil.addVarsToEqSystem(syst,jacDiffVars);
1243 4718 (adjMatrix, adjMatrixT) := BackendDAEUtil.adjacencyMatrix(syst1,BackendDAE.SPARSE(),NONE(),BackendDAEUtil.isInitializationDAE(inBackendDAE.shared));
1244 adjSize := arrayLength(adjMatrix) "number of equations";
1245
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4718 times.
 4718 adjSizeT := arrayLength(adjMatrixT) "number of variables";
1246
1247 // Debug dumping
1248
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4718 times.
 4718 if Flags.isSet(Flags.DUMP_SPARSE_VERBOSE) then
1249 ✗ BackendDump.printVarList(BackendVariable.varList(varswithDiffs));
1250 ✗ BackendDump.printEquationList(BackendEquation.equationList(orderedEqns));
1251 ✗ BackendDump.dumpAdjacencyMatrix(adjMatrix);
1252 ✗ BackendDump.dumpAdjacencyMatrixT(adjMatrixT);
1253 ✗ BackendDump.dumpFullMatching(bdaeMatching);
1254 end if;
1255
1256 // get indexes of diffed vars (rows)
1257 4718 nodesEqnsIndex := BackendVariable.getVarIndexFromVars(dependentVars,varswithDiffs);
1258 4718 nodesEqnsIndex := List.map1(nodesEqnsIndex, Array.getIndexFirst, ass1);
1259
1260 // debug dump
1261
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4718 times.
 4718 if Flags.isSet(Flags.DUMP_SPARSE_VERBOSE) then
1262 ✗ print("nodesEqnsIndexs: ");
1263 ✗ BackendDump.dumpAdjacencyRow(nodesEqnsIndex);
1264 ✗ print("\n");
1265 ✗ print("analytical Jacobians[" + patternName + "] -> build sparse graph: " + realString(clock()) + "\n");
1266 end if;
1267
1268 // prepare data for getSparsePattern
1269 4718 eqnSparse := arrayCreate(adjSize, {});
1270 4718 varSparse := arrayCreate(adjSizeT, {});
1271 4718 mark := arrayCreate(adjSizeT, 0);
1272 4718 usedvar := arrayCreate(adjSizeT, 0);
1273
1274 // make dependent variables as used if there are some
1275 // otherwise Array.setRange fails start is greater than end
1276
2/2
✓ Branch 0 taken 4692 times.
✓ Branch 1 taken 26 times.
 4718 if (sizeN>0) then
1277 4692 usedvar := Array.setRange(adjSizeT-(sizeN-1), adjSizeT, usedvar, 1);
1278 end if;
1279
1280 if debug then execStat("generateSparsePattern -> start "); end if;
1281 4718 eqnSparse := getSparsePattern(comps, eqnSparse, varSparse, mark, usedvar, 1, adjMatrix, adjMatrixT);
1282 if debug then execStat("generateSparsePattern -> end "); end if;
1283 // debug dump
1284
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4718 times.
 4718 if Flags.isSet(Flags.DUMP_SPARSE_VERBOSE) then
1285 ✗ BackendDump.dumpSparsePatternArray(eqnSparse);
1286 ✗ print("analytical Jacobians[" + patternName + "] -> prepared arrayList for transpose list: " + realString(clock()) + "\n");
1287 end if;
1288
1289 // select nodesEqnsIndex and map index to incoming vars
1290 4718 sparseArray := Array.select(eqnSparse, nodesEqnsIndex);
1291 4718 sparsepattern := arrayList(sparseArray);
1292
1293 4718 sparsepattern := List.map1List(sparsepattern, intSub, adjSizeT-sizeN);
1294 4718 sparseArray := listArray(sparsepattern);
1295
1296 if debug then execStat("generateSparsePattern -> postProcess "); end if;
1297
1298 // transpose the column-based pattern to row-based pattern
1299 4718 sparseArrayT := arrayCreate(sizeN,{});
1300 4718 sparseArrayT := transposeSparsePattern(sparsepattern, sparseArrayT, 1);
1301 4718 sparsepatternT := arrayList(sparseArrayT);
1302 4718 nonZeroElements := List.lengthListElements(sparsepattern);
1303 if debug then execStat("generateSparsePattern -> transpose done "); end if;
1304
1305
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4718 times.
 4718 if Flags.isSet(Flags.DUMP_SPARSE_VERBOSE) then
1306 // dump statistics
1307 ✗ dumpSparsePatternStatistics(nonZeroElements,sparsepatternT);
1308 ✗ BackendDump.dumpSparsePattern(sparsepattern);
1309 ✗ BackendDump.dumpSparsePattern(sparsepatternT);
1310 //execStat("generateSparsePattern -> nonZeroElements: " + intString(nonZeroElements) + " " ,ClockIndexes.RT_CLOCK_EXECSTAT_BACKEND_MODULES);
1311 end if;
1312
1313 // translated to DAE.ComRefs
1314
1/2
✓ Branch 0 taken 4718 times.
✗ Branch 1 not taken.
 4718 if listEmpty(sparsepattern) then
1315 sparsetuple := {};
1316 sparsetupleT := {};
1317 else
1318
8/8
✓ Branch 0 taken 23654 times.
✓ Branch 1 taken 4718 times.
✓ Branch 2 taken 23654 times.
✓ Branch 3 taken 4718 times.
✓ Branch 4 taken 79934 times.
✓ Branch 5 taken 23654 times.
✓ Branch 6 taken 79934 times.
✓ Branch 7 taken 23654 times.
 108306 translated := list(list(arrayGet(inDepCompRefs, i) for i in lst) for lst in sparsepattern);
1319
7/8
✓ Branch 0 taken 23654 times.
✓ Branch 1 taken 4718 times.
✓ Branch 3 taken 23654 times.
✓ Branch 4 taken 4718 times.
✓ Branch 5 taken 23654 times.
✓ Branch 6 taken 4718 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4718 times.
 56744 sparsetuple := list((cr,t) threaded for cr in depCompRefs, t in translated);
1320
8/8
✓ Branch 0 taken 24987 times.
✓ Branch 1 taken 4718 times.
✓ Branch 2 taken 24987 times.
✓ Branch 3 taken 4718 times.
✓ Branch 4 taken 79934 times.
✓ Branch 5 taken 24987 times.
✓ Branch 6 taken 79934 times.
✓ Branch 7 taken 24987 times.
 109639 translated := list(list(arrayGet(depCompRefs, i) for i in lst) for lst in sparsepatternT);
1321
7/8
✓ Branch 0 taken 24987 times.
✓ Branch 1 taken 4718 times.
✓ Branch 3 taken 24987 times.
✓ Branch 4 taken 4718 times.
✓ Branch 5 taken 24987 times.
✓ Branch 6 taken 4718 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4718 times.
 59410 sparsetupleT := list((cr,t) threaded for cr in inDepCompRefs, t in translated);
1322 end if;
1323
1324 if debug then execStat("generateSparsePattern -> coloring start "); end if;
1325
3/4
✓ Branch 0 taken 3058 times.
✓ Branch 1 taken 1660 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3058 times.
 4718 if nonlinearPattern or not withColoring or Flags.isSet(Flags.DISABLE_COLORING) then
1326 //without coloring
1327
2/2
✓ Branch 0 taken 8766 times.
✓ Branch 1 taken 1660 times.
 10426 coloring := list({arrayGet(inDepCompRefs, i)} for i in 1:sizeN);
1328 else
1329 // get coloring based on sparse pattern
1330 3058 coloredArray := Coloring.createColoring(sparseArray, sparseArrayT, sizeN, sizeM);
1331
8/8
✓ Branch 0 taken 6963 times.
✓ Branch 1 taken 3058 times.
✓ Branch 3 taken 6963 times.
✓ Branch 4 taken 3058 times.
✓ Branch 5 taken 16221 times.
✓ Branch 6 taken 6963 times.
✓ Branch 7 taken 16221 times.
✓ Branch 8 taken 6963 times.
 29300 coloring := list(list(arrayGet(inDepCompRefs, i) for i in lst) for lst in coloredArray);
1332 end if;
1333 if debug then execStat("generateSparsePattern -> coloring done "); end if;
1334
1335
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4718 times.
 4718 if Flags.isSet(Flags.DUMP_SPARSE_VERBOSE) then
1336 ✗ print("analytical Jacobians[" + patternName + "] -> ready! " + realString(clock()) + "\n");
1337 end if;
1338
1339 4718 outSparsePattern := (sparsetupleT, sparsetuple, (inDepCompRefsLst, depCompRefsLst), nonZeroElements);
1340
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4718 times.
 4718 if Flags.isSet(Flags.DUMP_SPARSE) then
1341 ✗ BackendDump.dumpSparsityPattern(outSparsePattern, " --- " + patternName + " Pattern ---");
1342 ✗ BackendDump.dumpSparseColoring(coloring, " --- " + patternName + " Coloring ---");
1343 end if;
1344 if debug then execStat("generateSparsePattern -> final end "); end if;
1345 then (outSparsePattern, coloring);
1346 else
1347 algorithm
1348 ✗ Error.addInternalError("function generateSparsePattern failed", sourceInfo());
1349 ✗ then fail();
1350 end matchcontinue;
1351 end generateSparsePattern;
1352
1353 protected function dumpSparsePatternStatistics
1354 input Integer nonZeroElements;
1355 input list<list<Integer>> sparsepatternT;
1356 protected
1357 Integer maxDegree;
1358 algorithm
1359 ✗ (_, maxDegree) := List.mapFold(sparsepatternT, findDegrees, 0);
1360 ✗ print("analytical Jacobians[SPARSE] -> got sparse pattern nonZeroElements: "+ String(nonZeroElements) + " maxNodeDegree: " + String(maxDegree) + " time : " + String(clock()) + "\n");
1361 end dumpSparsePatternStatistics;
1362
1363 protected function findDegrees<T>
1364 input list<T> inList;
1365 input Integer inValue;
1366 output Integer outDegree;
1367 output Integer outMaxDegree;
1368 algorithm
1369 ✗ outDegree := listLength(inList);
1370 outMaxDegree := intMax(inValue, outDegree);
1371 end findDegrees;
1372
1373 protected function getSparsePattern
1374 input BackendDAE.StrongComponents inComponents;
1375 input array<list<Integer>> ineqnSparse; //
1376 input array<list<Integer>> invarSparse; //
1377 input array<Integer> inMark; //
1378 input array<Integer> inUsed; //
1379 input Integer inmarkValue;
1380 input BackendDAE.AdjacencyMatrix inMatrix;
1381 input BackendDAE.AdjacencyMatrix inMatrixT;
1382 output array<list<Integer>> outSparsePattern;
1383 algorithm
1384 outSparsePattern := match (inComponents, ineqnSparse)
1385 local
1386 list<Integer> vars, vars1, eqns, eqns1;
1387 list<Integer> inputVars;
1388 list<list<Integer>> inputVarsLst;
1389 list<Integer> solvedVars;
1390 array<list<Integer>> result;
1391 Integer var, eqn;
1392 BackendDAE.StrongComponents rest;
1393 BackendDAE.StrongComponent comp;
1394 BackendDAE.InnerEquations innerEquations;
1395 case ({}, result) then result;
1396
1397 case(BackendDAE.SINGLEEQUATION(eqn=eqn,var=var)::rest, result)
1398 algorithm
1399 98835 inputVars := arrayGet(inMatrix, eqn);
1400 98835 inputVars := List.removeOnTrue(var, intEq, inputVars);
1401
1402 98835 getSparsePattern2(inputVars, {var}, {eqn}, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1403
1404 98835 result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1405 then result;
1406 case(BackendDAE.SINGLEARRAY(eqn=eqn,vars=solvedVars)::rest, result)
1407 algorithm
1408 216 inputVars := arrayGet(inMatrix, eqn);
1409
6/6
✓ Branch 1 taken 1866 times.
✓ Branch 2 taken 146 times.
✓ Branch 3 taken 2012 times.
✓ Branch 4 taken 216 times.
✓ Branch 5 taken 146 times.
✓ Branch 6 taken 216 times.
 2228 inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1410
1411 216 getSparsePattern2(inputVars, solvedVars, {eqn}, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1412
1413 216 result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1414 then result;
1415 case(BackendDAE.SINGLEIFEQUATION(eqn=eqn,vars=solvedVars)::rest, result)
1416 algorithm
1417 ✗ inputVars := arrayGet(inMatrixT, eqn);
1418 ✗ inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1419
1420 ✗ getSparsePattern2(inputVars, solvedVars, {eqn}, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1421
1422 ✗ result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1423 then result;
1424 case(BackendDAE.SINGLEALGORITHM(eqn=eqn,vars=solvedVars)::rest, result)
1425 algorithm
1426 126 inputVars := arrayGet(inMatrix, eqn);
1427
6/6
✓ Branch 1 taken 715 times.
✓ Branch 2 taken 273 times.
✓ Branch 3 taken 988 times.
✓ Branch 4 taken 126 times.
✓ Branch 5 taken 273 times.
✓ Branch 6 taken 126 times.
 1114 inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1428
1429 126 getSparsePattern2(inputVars, solvedVars, {eqn}, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1430
1431 126 result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1432 then result;
1433 case(BackendDAE.SINGLECOMPLEXEQUATION(eqn=eqn,vars=solvedVars)::rest, result)
1434 algorithm
1435 874 inputVars := arrayGet(inMatrix, eqn);
1436
6/6
✓ Branch 1 taken 8862 times.
✓ Branch 2 taken 1946 times.
✓ Branch 3 taken 10808 times.
✓ Branch 4 taken 874 times.
✓ Branch 5 taken 1946 times.
✓ Branch 6 taken 874 times.
 11682 inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1437
1438 874 getSparsePattern2(inputVars, solvedVars, {eqn}, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1439
1440 874 result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1441 then result;
1442 case(BackendDAE.SINGLEWHENEQUATION(eqn=eqn,vars=solvedVars)::rest, result)
1443 algorithm
1444 277 inputVars := arrayGet(inMatrix, eqn);
1445
6/6
✓ Branch 1 taken 277 times.
✓ Branch 2 taken 488 times.
✓ Branch 3 taken 765 times.
✓ Branch 4 taken 277 times.
✓ Branch 5 taken 488 times.
✓ Branch 6 taken 277 times.
 1042 inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1446
1447 277 getSparsePattern2(inputVars, solvedVars, {eqn}, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1448
1449 277 result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1450 then result;
1451 case(BackendDAE.SINGLEIFEQUATION(eqn=eqn,vars=solvedVars)::rest, result)
1452 algorithm
1453 ✗ inputVars := arrayGet(inMatrix, eqn);
1454 ✗ inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1455
1456 ✗ getSparsePattern2(inputVars, solvedVars, {eqn}, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1457
1458 ✗ result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1459 then result;
1460 case(BackendDAE.EQUATIONSYSTEM(eqns=eqns,vars=solvedVars)::rest, result)
1461 algorithm
1462 898 inputVarsLst := List.map1(eqns, Array.getIndexFirst, inMatrix);
1463 898 inputVars := List.flatten(inputVarsLst);
1464
6/6
✓ Branch 1 taken 34528 times.
✓ Branch 2 taken 24148 times.
✓ Branch 3 taken 58676 times.
✓ Branch 4 taken 898 times.
✓ Branch 5 taken 24148 times.
✓ Branch 6 taken 898 times.
 59574 inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1465
1466 898 getSparsePattern2(inputVars, solvedVars, eqns, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1467
1468 898 result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1469 then result;
1470 case(BackendDAE.TORNSYSTEM(BackendDAE.TEARINGSET(residualequations=eqns,tearingvars=vars,innerEquations=innerEquations))::rest, result)
1471 algorithm
1472 9 (eqns1,inputVarsLst,_) := List.map_3(innerEquations, BackendDAEUtil.getEqnAndVarsFromInnerEquation);
1473 9 vars1 := List.flatten(inputVarsLst);
1474 9 eqns1 := listAppend(eqns, eqns1);
1475 9 solvedVars := listAppend(vars, vars1);
1476
1477 9 inputVarsLst := List.map1(eqns1, Array.getIndexFirst, inMatrix);
1478 9 inputVars := List.flatten(inputVarsLst);
1479
6/6
✓ Branch 1 taken 78 times.
✓ Branch 2 taken 114 times.
✓ Branch 3 taken 192 times.
✓ Branch 4 taken 9 times.
✓ Branch 5 taken 114 times.
✓ Branch 6 taken 9 times.
 201 inputVars := list(v for v guard not listMember(v, solvedVars) in inputVars);
1480
1481 9 getSparsePattern2(inputVars, solvedVars, eqns1, ineqnSparse, invarSparse, inMark, inUsed, inmarkValue);
1482
1483 9 result := getSparsePattern(rest, result, invarSparse, inMark, inUsed, inmarkValue+1, inMatrix, inMatrixT);
1484 then result;
1485 else
1486 algorithm
1487 ✗ comp::_ := inComponents;
1488 ✗ BackendDump.dumpComponent(comp);
1489 ✗ Error.addInternalError("function getSparsePattern failed", sourceInfo());
1490 ✗ then fail();
1491 end match;
1492 end getSparsePattern;
1493
1494 protected function getSparsePattern2
1495 input list<Integer> inInputVars;
1496 input list<Integer> inSolvedVars;
1497 input list<Integer> inEqns;
1498 input array<list<Integer>> ineqnSparse;
1499 input array<list<Integer>> invarSparse;
1500 input array<Integer> inMark;
1501 input array<Integer> inUsed;
1502 input Integer inmarkValue;
1503 protected
1504 list<Integer> localList;
1505 algorithm
1506 101235 localList := getSparsePatternHelp(inInputVars, invarSparse, inMark, inUsed, inmarkValue);
1507 101235 List.map2_0(inSolvedVars, Array.updateIndexFirst, localList, invarSparse);
1508 101235 List.map2_0(inEqns, Array.updateIndexFirst, localList, ineqnSparse);
1509 end getSparsePattern2;
1510
1511 protected function getSparsePatternHelp
1512 input list<Integer> inInputVars;
1513 input array<list<Integer>> invarSparse;
1514 input array<Integer> inMark;
1515 input array<Integer> inUsed;
1516 input Integer inmarkValue;
1517 output list<Integer> outLocalList = {};
1518 protected
1519 Integer arrayElement;
1520 list<Integer> varSparse;
1521 algorithm
1522
2/2
✓ Branch 0 taken 224949 times.
✓ Branch 1 taken 101235 times.
 326184 for var in inInputVars loop
1523 224949 arrayElement := arrayGet(inUsed, var);
1524
2/2
✓ Branch 0 taken 67725 times.
✓ Branch 1 taken 157224 times.
 224949 if intEq(1, arrayElement) then
1525 67725 arrayElement := arrayGet(inMark, var);
1526
2/2
✓ Branch 0 taken 66051 times.
✓ Branch 1 taken 1674 times.
 67725 if not intEq(inmarkValue, arrayElement) then
1527 66051 arrayUpdate(inMark, var, inmarkValue);
1528 outLocalList := var::outLocalList;
1529 end if;
1530 end if;
1531
1532 224949 varSparse := arrayGet(invarSparse, var);
1533
2/2
✓ Branch 0 taken 658601 times.
✓ Branch 1 taken 224949 times.
 883550 for v in varSparse loop
1534 658601 arrayElement := arrayGet(inMark, v);
1535
2/2
✓ Branch 0 taken 287924 times.
✓ Branch 1 taken 370677 times.
 658601 if not intEq(inmarkValue, arrayElement) then
1536 287924 arrayUpdate(inMark, v, inmarkValue);
1537 outLocalList := v::outLocalList;
1538 end if;
1539 end for;
1540 end for;
1541 end getSparsePatternHelp;
1542
1543 public function transposeSparsePattern
1544 input list<list<Integer>> inSparsePattern;
1545 input array<list<Integer>> inAccumList;
1546 input Integer inValue;
1547 output array<list<Integer>> outSparsePattern = inAccumList;
1548 protected
1549 Integer value = inValue;
1550 list<Integer> tmplist;
1551 algorithm
1552
2/2
✓ Branch 0 taken 23654 times.
✓ Branch 1 taken 4718 times.
 28372 for oneList in inSparsePattern loop
1553
2/2
✓ Branch 0 taken 79934 times.
✓ Branch 1 taken 23654 times.
 103588 for oneElem in oneList loop
1554 79934 tmplist := arrayGet(outSparsePattern,oneElem);
1555 MetaModelica.Dangerous.arrayUpdateNoBoundsChecking(outSparsePattern, oneElem, value::tmplist);
1556 end for;
1557 23654 value := value + 1;
1558 end for;
1559 end transposeSparsePattern;
1560
1561 public function transposeSparsePatternTuple
1562 input list<tuple<Integer, list<Integer>>> inSparsePattern;
1563 input array<tuple<Integer,list<Integer>>> inAccumList;
1564 output array<tuple<Integer,list<Integer>>> outSparsePattern = inAccumList;
1565 protected
1566 Integer value;
1567 list<Integer> tmplist;
1568 list<Integer> oneList;
1569 tuple<Integer,list<Integer>> tmpTuple;
1570 Integer i;
1571 algorithm
1572 ✗ for oneListTuple in inSparsePattern loop
1573 ✗ (value, oneList) := oneListTuple;
1574 ✗ for oneElem in oneList loop
1575 ✗ tmpTuple := arrayGet(outSparsePattern,oneElem+1);
1576 ✗ (_, tmplist) := tmpTuple;
1577 tmplist := value::tmplist;
1578 ✗ tmpTuple := (oneElem, tmplist);
1579 MetaModelica.Dangerous.arrayUpdateNoBoundsChecking(outSparsePattern, oneElem+1, tmpTuple);
1580 end for;
1581 end for;
1582 // sort all transposed lists
1583 ✗ for i in 1:listLength(inSparsePattern) loop
1584 ✗ tmpTuple := arrayGet(outSparsePattern,i);
1585 ✗ (value, tmplist) := tmpTuple;
1586 ✗ tmplist := List.heapSortIntList(tmplist);
1587 ✗ tmpTuple := (value, tmplist);
1588 MetaModelica.Dangerous.arrayUpdateNoBoundsChecking(outSparsePattern, i, tmpTuple);
1589 end for;
1590 end transposeSparsePatternTuple;
1591
1592 protected function createInDepVars
1593 "This function creates variables for the dependecy
1594 analysis, this needs to cosider different behavoir
1595 clock stated and continuous states.
1596 continuous states: der(x) > dependent and x > independent
1597 clocked states: previous(x) > independent and x > dependent
1598 "
1599 input list<BackendDAE.Var> independentVars;
1600 input Boolean createpDerStates = true;
1601 output list<BackendDAE.Var> outVars = {};
1602 output list<DAE.ComponentRef> outCrefs = {};
1603 protected
1604 BackendDAE.Var var;
1605 algorithm
1606
2/2
✓ Branch 0 taken 30998 times.
✓ Branch 1 taken 6538 times.
 37536 for v in independentVars loop
1607
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 30986 times.
 30998 if BackendVariable.isClockedStateVar(v) then
1608 12 var := BackendVariable.createClockedState(v);
1609 outVars := var::outVars;
1610 12 outCrefs := var.varName::outCrefs;
1611 elseif createpDerStates then
1612 24975 outVars := BackendVariable.createpDerVar(v)::outVars;
1613 24975 outCrefs := v.varName::outCrefs;
1614 else
1615 outVars := v::outVars;
1616 6011 outCrefs := v.varName::outCrefs;
1617 end if;
1618 end for;
1619 6538 outVars := listReverse(outVars);
1620 6538 outCrefs := listReverse(outCrefs);
1621 end createInDepVars;
1622
1623 public function createFMIModelDerivatives
1624 "This function genererate the stucture output and the
1625 partial derivatives for FMI, which are basically the jacobian matrices.
1626 author: wbraun"
1627 input BackendDAE.BackendDAE inBackendDAE;
1628 output BackendDAE.SymbolicJacobians outJacobianMatrices = {};
1629 output AvlTreePathFunction.Tree outFunctionTree;
1630 protected
1631 BackendDAE.BackendDAE backendDAE,emptyBDAE;
1632 BackendDAE.EqSystem eqSyst;
1633 Option<BackendDAE.SymbolicJacobian> outJacobian;
1634
1635 list<BackendDAE.Var> varlst, knvarlst, states, inputvars, outputvars, paramvars, indepVars, depVars;
1636
1637 BackendDAE.Variables v,globalKnownVars,statesarr,inputvarsarr,paramvarsarr,depVarsArr;
1638
1639 BackendDAE.SparsePattern sparsePattern;
1640 BackendDAE.SparseColoring sparseColoring;
1641 BackendDAE.NonlinearPattern nonlinearPattern;
1642
1643 AvlTreePathFunction.Tree functionTree;
1644
1645 BackendDAE.ExtraInfo ei;
1646 FCore.Cache cache;
1647 FCore.Graph graph;
1648 algorithm
1649 // Dependency analysis only, nothing to differentiate, and one partition: take
1650 // the pattern from the system as it stands rather than causalizing a collapsed
1651 // copy of it. Several partitions are clocked ones, which sample each other's
1652 // variables, so those still have to be collapsed to be seen across.
1653
4/4
✓ Branch 1 taken 67 times.
✓ Branch 2 taken 13 times.
✓ Branch 4 taken 41 times.
✓ Branch 5 taken 26 times.
 80 if Flags.isSet(Flags.DIS_SYMJAC_FMI20) and listLength(inBackendDAE.eqs) == 1 then
1654 41 (sparsePattern, sparseColoring) := fmiDerSparsePattern(inBackendDAE);
1655 82 outJacobianMatrices := {(
1656 SOME((BackendDAE.DAE({BackendDAEUtil.createEqSystem(BackendVariable.emptyVars(), BackendEquation.emptyEqns())},
1657 BackendDAEUtil.createEmptyShared(BackendDAE.JACOBIAN(), inBackendDAE.shared.info,
1658 inBackendDAE.shared.cache, inBackendDAE.shared.graph)),
1659 "FMIDER", {}, {}, {}, {})),
1660 sparsePattern, sparseColoring, BackendDAE.emptyNonlinearPattern)};
1661 41 outFunctionTree := inBackendDAE.shared.functionTree;
1662 41 return;
1663 end if;
1664 try
1665 // for now perform on collapsed system
1666 39 backendDAE := BackendDAEUtil.copyBackendDAE(inBackendDAE);
1667 39 backendDAE := BackendDAEOptimize.collapseIndependentBlocks(backendDAE);
1668 39 backendDAE := BackendDAEUtil.transformBackendDAE(backendDAE,SOME((BackendDAE.NO_INDEX_REDUCTION(),BackendDAE.EXACT())),NONE(),NONE());
1669
1670 // get all variables
1671
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 39 times.
 39 eqSyst::{} := backendDAE.eqs;
1672 39 v := eqSyst.orderedVars;
1673 39 globalKnownVars := backendDAE.shared.globalKnownVars;
1674
1675 // prepare all needed variables
1676 39 varlst := BackendVariable.varList(v);
1677 39 knvarlst := BackendVariable.varList(globalKnownVars);
1678
1679
1/2
✓ Branch 1 taken 39 times.
✗ Branch 2 not taken.
 39 states := if Config.languageStandardAtLeast(Config.LanguageStandard._3_3) then
1680 BackendVariable.getAllClockedStatesFromVariables(v) else {};
1681
1682 39 states := listAppend(BackendVariable.getAllStateVarFromVariables(v), states);
1683
1684 39 inputvars := List.select(knvarlst,BackendVariable.isVarOnTopLevelAndInput);
1685 39 outputvars := List.select(varlst, BackendVariable.isVarOnTopLevelAndOutput);
1686
1687 // independent varibales states + inputs
1688 39 indepVars := listAppend(states, inputvars);
1689
1690 // dependent varibales der(states) + outputs
1691 39 depVars := listAppend(states, outputvars);
1692
1693 // Generate sparse pattern for matrices states
1694 // prepare more needed variables
1695
2/2
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 13 times.
 39 if Flags.isSet(Flags.DIS_SYMJAC_FMI20) then
1696 // empty BackendDAE in case derivates should not calclulated
1697 26 cache := backendDAE.shared.cache;
1698 26 graph := backendDAE.shared.graph;
1699 26 ei := backendDAE.shared.info;
1700 52 emptyBDAE := BackendDAE.DAE({BackendDAEUtil.createEqSystem(BackendVariable.emptyVars(), BackendEquation.emptyEqns())}, BackendDAEUtil.createEmptyShared(BackendDAE.JACOBIAN(), ei, cache, graph));
1701
1702 26 (sparsePattern, sparseColoring) := generateSparsePattern(backendDAE, indepVars, depVars, withColoring = false);
1703
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 26 times.
 26 if Flags.isSet(Flags.JAC_DUMP2) then
1704 ✗ BackendDump.dumpSparsityPattern(sparsePattern, "FMI sparsity");
1705 end if;
1706 26 outJacobianMatrices := (SOME((emptyBDAE,"FMIDER",{},{},{}, {})), sparsePattern, sparseColoring, BackendDAE.emptyNonlinearPattern)::outJacobianMatrices;
1707 26 outFunctionTree := inBackendDAE.shared.functionTree;
1708 else
1709 // prepare more needed variables
1710 13 paramvars := List.select(knvarlst, BackendVariable.isParam);
1711 13 statesarr := BackendVariable.listVar1(states);
1712 13 inputvarsarr := BackendVariable.listVar1(inputvars);
1713 13 paramvarsarr := BackendVariable.listVar1(paramvars);
1714 13 depVarsArr := BackendVariable.listVar1(depVars);
1715
1716 13 (outJacobian, outFunctionTree, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE,indepVars,statesarr,inputvarsarr,paramvarsarr,depVarsArr,varlst,"FMIDER", Flags.isSet(Flags.DIS_SYMJAC_FMI20));
1717
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 13 times.
 13 if Flags.isSet(Flags.JAC_DUMP2) then
1718 ✗ BackendDump.dumpSparsityPattern(sparsePattern, "FMI sparsity");
1719 end if;
1720 13 outJacobianMatrices := (outJacobian, sparsePattern, sparseColoring, nonlinearPattern)::outJacobianMatrices;
1721 13 outFunctionTree := AvlTreePathFunction.join(inBackendDAE.shared.functionTree, outFunctionTree);
1722 end if;
1723 else
1724 ✗ Error.addInternalError("function createFMIModelDerivatives failed", sourceInfo());
1725 outJacobianMatrices := {};
1726 ✗ outFunctionTree := inBackendDAE.shared.functionTree;
1727 end try;
1728 end createFMIModelDerivatives;
1729
1730 protected function fmiDerSparsePattern
1731 "The FMIDER dependency pattern of a DAE that is a single partition, taken as it
1732 stands: collapsing it is a no-op merge that drops the matching only for
1733 transformBackendDAE to compute it again."
1734 input BackendDAE.BackendDAE inDAE;
1735 output BackendDAE.SparsePattern outSparsePattern;
1736 output BackendDAE.SparseColoring outColoring;
1737 protected
1738 // generateSparsePattern adds the seed variables to the system it is given.
1739 BackendDAE.BackendDAE dae = BackendDAEUtil.copyBackendDAE(inDAE);
1740 BackendDAE.EqSystem syst = listHead(dae.eqs);
1741 list<BackendDAE.Var> states, inputvars, outputvars;
1742 algorithm
1743
1/2
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
 41 states := if Config.languageStandardAtLeast(Config.LanguageStandard._3_3) then
1744 BackendVariable.getAllClockedStatesFromVariables(syst.orderedVars) else {};
1745 41 states := listAppend(BackendVariable.getAllStateVarFromVariables(syst.orderedVars), states);
1746 41 outputvars := List.select(BackendVariable.varList(syst.orderedVars), BackendVariable.isVarOnTopLevelAndOutput);
1747 41 inputvars := List.select(BackendVariable.varList(dae.shared.globalKnownVars), BackendVariable.isVarOnTopLevelAndInput);
1748
1749 41 (outSparsePattern, outColoring) := generateSparsePattern(dae, listAppend(states, inputvars),
1750 listAppend(states, outputvars), withColoring = false);
1751 end fmiDerSparsePattern;
1752
1753 public function createFMIModelDerivativesForInitialization
1754 "This function genererate the stucture output and the
1755 partial derivatives for FMI, which are basically the jacobian matrices."
1756 input BackendDAE.BackendDAE initDAE;
1757 input BackendDAE.BackendDAE simDAE;
1758 input list<BackendDAE.Var> depVars;
1759 input list<BackendDAE.Var> indepVars;
1760 input BackendDAE.Variables orderedVars;
1761 input BackendDAE.SparsePattern sparsePattern_;
1762 input BackendDAE.SparseColoring sparseColoring_;
1763 output BackendDAE.SymbolicJacobians outJacobianMatrices = {};
1764 output AvlTreePathFunction.Tree outFunctionTree "may contain functions created by the differentiation, e.g. partial derivatives";
1765 protected
1766 BackendDAE.BackendDAE backendDAE_1, emptyBDAE;
1767 BackendDAE.EqSystem currentSystem;
1768 Option<BackendDAE.SymbolicJacobian> outJacobian;
1769 list<BackendDAE.Var> varlst, knvarlst, states, clockedStates, inputvars, paramvars;
1770 BackendDAE.Variables statesarr, inputvarsarr, paramvarsarr, depVarsArr;
1771 BackendDAE.ExtraInfo ei;
1772 FCore.Cache cache;
1773 FCore.Graph graph;
1774 BackendDAE.EquationArray newOrderedEquationArray;
1775 BackendDAE.Shared shared;
1776 DAE.Exp lhs, rhs;
1777 BackendDAE.Equation eqn;
1778 DAE.ComponentRef cr, rhsCr;
1779 UnorderedSet<DAE.ComponentRef> crefsVarsToRemove, protectedCrefs;
1780 BackendDAE.Variables newVars;
1781 algorithm
1782 // Generate empty jacobian martices
1783
2/2
✓ Branch 1 taken 57 times.
✓ Branch 2 taken 12 times.
 69 if Flags.isSet(Flags.DIS_SYMJAC_FMI20) then
1784 57 cache := initDAE.shared.cache;
1785 57 graph := initDAE.shared.graph;
1786 57 ei := initDAE.shared.info;
1787 114 emptyBDAE := BackendDAE.DAE({BackendDAEUtil.createEqSystem(BackendVariable.emptyVars(), BackendEquation.emptyEqns())}, BackendDAEUtil.createEmptyShared(BackendDAE.JACOBIAN(), ei, cache, graph));
1788 114 outJacobianMatrices := (SOME((emptyBDAE,"FMIDERINIT",{},{},{}, {})), BackendDAE.emptySparsePattern, {}, BackendDAE.emptyNonlinearPattern)::outJacobianMatrices;
1789 57 outFunctionTree := initDAE.shared.functionTree;
1790 57 return;
1791 end if;
1792 try
1793
1794 12 backendDAE_1 := BackendDAEUtil.copyBackendDAE(initDAE);
1795 12 backendDAE_1 := BackendDAEOptimize.collapseIndependentBlocks(backendDAE_1);
1796
1797 //BackendDump.printBackendDAE(backendDAE_1);
1798 //BackendDump.dumpVariables(simDAE.shared.globalKnownVars, "check global vars");
1799
1800 /* add the calculated parameter equations here which does not have constant binding
1801 parameter Real x = 10;
1802 Real m = x; */
1803
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
 12 BackendDAE.DAE(currentSystem::{}, shared) := backendDAE_1;
1804 12 protectedCrefs := UnorderedSet.new(ComponentReferenceBasics.hashComponentRef, ComponentReferenceBasics.crefEqual);
1805
2/2
✓ Branch 0 taken 30 times.
✓ Branch 1 taken 12 times.
 42 for var in depVars loop
1806 30 UnorderedSet.add(var.varName, protectedCrefs);
1807
3/4
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 29 times.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
 30 if BackendVariable.isParam(var) and not BackendVariable.varHasConstantBindExp(var) then
1808 //print("\n PARAM_CHECK: " + ComponentReferenceBasics.printComponentRefStr(var.varName));
1809 1 lhs := BackendVariable.varExp(var);
1810 1 rhs := BackendVariable.varBindExpStartValueNoFail(var) "bindings are optional";
1811 1 eqn := BackendDAE.EQUATION(lhs, rhs, DAE.emptyElementSource, BackendDAE.EQ_ATTR_DEFAULT_BINDING);
1812 //BackendDump.printEquation(eqn);
1813 1 BackendEquation.add(eqn, currentSystem.orderedEqs);
1814
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
 1 if not BackendVariable.containsCref(var.varName, currentSystem.orderedVars) then
1815 1 currentSystem := BackendVariable.addVarDAE(BackendVariable.makeVar(var.varName), currentSystem);
1816 end if;
1817 end if;
1818 end for;
1819
1820 // Remove initialization start-value helper variables from the system used for
1821 // symbolic differentiation. Differentiate treats $START.* crefs as constants,
1822 // so keeping these variables would create derivative variables without
1823 // remaining equations after simplification.
1824 12 newOrderedEquationArray := BackendEquation.emptyEqns();
1825 12 crefsVarsToRemove := UnorderedSet.new(ComponentReferenceBasics.hashComponentRef, ComponentReferenceBasics.crefEqual);
1826
2/2
✓ Branch 1 taken 53 times.
✓ Branch 2 taken 12 times.
 77 for eq in BackendEquation.equationList(currentSystem.orderedEqs) loop
1827
1/2
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
 53 if not BackendEquation.isAlgorithm(eq) then
1828 53 lhs := BackendEquation.getEquationLHS(eq);
1829 53 rhs := BackendEquation.getEquationRHS(eq);
1830
2/2
✓ Branch 1 taken 47 times.
✓ Branch 2 taken 6 times.
 53 if Expression.isExpCref(lhs) then
1831 47 cr := Expression.expCref(lhs);
1832 // remove lhs equation of type $Start.a = ... as it does not contribute to the jacobian and create a variable a with constant binding which is not wanted
1833
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 44 times.
 47 if ComponentReference.isStartCref(cr) then
1834 3 UnorderedSet.add(cr, crefsVarsToRemove);
1835 elseif Expression.isExpCref(rhs) and not UnorderedSet.contains(cr, protectedCrefs) then
1836 8 rhsCr := Expression.expCref(rhs);
1837 // remove equation of form a = $START.a as it does not contribute to the jacobian and create a variable a with constant binding which is not wanted
1838
3/4
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
 8 if ComponentReference.isStartCref(rhsCr) and ComponentReferenceBasics.crefEqual(ComponentReference.popCref(rhsCr), cr) then
1839 6 UnorderedSet.add(cr, crefsVarsToRemove);
1840 else
1841 2 BackendEquation.add(eq, newOrderedEquationArray);
1842 end if;
1843 else
1844 36 BackendEquation.add(eq, newOrderedEquationArray);
1845 end if;
1846 else
1847 6 BackendEquation.add(eq, newOrderedEquationArray);
1848 end if;
1849 else
1850 ✗ BackendEquation.add(eq, newOrderedEquationArray);
1851 end if;
1852 end for;
1853
1854 12 newVars := BackendVariable.emptyVars();
1855
2/2
✓ Branch 1 taken 53 times.
✓ Branch 2 taken 12 times.
 77 for var in BackendVariable.varList(currentSystem.orderedVars) loop
1856
2/2
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 9 times.
 53 if not UnorderedSet.contains(var.varName, crefsVarsToRemove) then
1857 // make depVars crefs as unreplaceable as it might be removed by removeSimpleEquation and Optimization fails for jacobians
1858
2/2
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 14 times.
 44 if UnorderedSet.contains(var.varName, protectedCrefs) then
1859 30 var := BackendVariable.setVarUnreplaceable(var, true);
1860 end if;
1861 44 newVars := BackendVariable.addVar(var, newVars);
1862 end if;
1863 end for;
1864
1865 12 currentSystem := BackendDAEUtil.setEqSystEqs(currentSystem, newOrderedEquationArray);
1866 12 currentSystem := BackendDAEUtil.setEqSystVars(currentSystem, newVars);
1867
1868 // put the shared globalknown Vars
1869 // for var in BackendVariable.varList(simDAE.shared.globalKnownVars) loop
1870 // if not BackendVariable.containsCref(var.varName, currentSystem.orderedVars) then
1871 // shared := BackendVariable.addGlobalKnownVarDAE(var, shared);
1872 // end if;
1873 // end for;
1874
1875
1876 12 backendDAE_1 := BackendDAE.DAE({currentSystem}, shared);
1877 12 backendDAE_1 := BackendDAEOptimize.collapseIndependentBlocks(backendDAE_1);
1878 12 backendDAE_1 := BackendDAEUtil.transformBackendDAE(backendDAE_1, SOME((BackendDAE.NO_INDEX_REDUCTION(),BackendDAE.EXACT())),NONE(),NONE());
1879
1880 //BackendDump.printBackendDAE(backendDAE_1);
1881
1882 // Only the state variables are read from the simulation DAE, so it needs
1883 // neither a copy nor a collapse. The finders cons and collapsing folds the
1884 // systems in reverse, so reading them forwards keeps the old order.
1885 states := {};
1886 clockedStates := {};
1887
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 12 times.
 24 for syst in simDAE.eqs loop
1888 12 states := List.append_reverse(BackendVariable.getAllStateVarFromVariables(syst.orderedVars), states);
1889
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 if Config.languageStandardAtLeast(Config.LanguageStandard._3_3) then
1890 12 clockedStates := List.append_reverse(BackendVariable.getAllClockedStatesFromVariables(syst.orderedVars), clockedStates);
1891 end if;
1892 end for;
1893 12 states := listAppend(listReverse(states), listReverse(clockedStates));
1894
1895 // prepare all needed variables from initialization DAE
1896 12 varlst := BackendVariable.varList(currentSystem.orderedVars);
1897 12 knvarlst := BackendVariable.varList(simDAE.shared.globalKnownVars);
1898 //BackendDump.dumpVarList(knvarlst, "shared simulation DAE");
1899 12 inputvars := List.select(knvarlst, BackendVariable.isVarOnTopLevelAndInput);
1900
1901 // prepare more needed variables
1902 12 paramvars := List.select(knvarlst, BackendVariable.isParam);
1903 12 statesarr := BackendVariable.listVar1(states);
1904 12 inputvarsarr := BackendVariable.listVar1(inputvars);
1905 12 paramvarsarr := BackendVariable.listVar1(paramvars);
1906 12 depVarsArr := BackendVariable.listVar1(depVars);
1907
1908 //(outJacobian, outFunctionTree, _, _) := generateGenericJacobian(backendDAE_1, indepVars, BackendVariable.emptyVars(), BackendVariable.emptyVars(), BackendVariable.emptyVars(), depVarsArr, depVars, "FMIDERINIT", Flags.isSet(Flags.DIS_SYMJAC_FMI20));
1909 12 (outJacobian, outFunctionTree, _, _) := generateGenericJacobian(backendDAE_1, indepVars, statesarr, inputvarsarr, paramvarsarr, depVarsArr, varlst, "FMIDERINIT", false);
1910
1911
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 if Flags.isSet(Flags.JAC_DUMP2) then
1912 ✗ BackendDump.dumpSparsityPattern(sparsePattern_, "FMI sparsity");
1913 end if;
1914 // kabdelhak: maybe also pass nonlinearity pattern to add it here
1915 12 outJacobianMatrices := (outJacobian, sparsePattern_, sparseColoring_, BackendDAE.emptyNonlinearPattern)::outJacobianMatrices;
1916 12 outFunctionTree := AvlTreePathFunction.join(initDAE.shared.functionTree, outFunctionTree);
1917 else
1918 ✗ Error.addInternalError("function createFMIModelDerivativesForInitialization failed", sourceInfo());
1919 outJacobianMatrices := {};
1920 ✗ outFunctionTree := initDAE.shared.functionTree;
1921 end try;
1922 end createFMIModelDerivativesForInitialization;
1923
1924 protected function createLinearModelMatrices "This function creates the linear model matrices column-wise
1925 author: wbraun"
1926 input BackendDAE.BackendDAE inBackendDAE;
1927 input Boolean useOptimica;
1928 output BackendDAE.SymbolicJacobians outJacobianMatrices;
1929 output AvlTreePathFunction.Tree outFunctionTree;
1930
1931 algorithm
1932 (outJacobianMatrices, outFunctionTree) :=
1933 match (inBackendDAE, useOptimica)
1934 local
1935 BackendDAE.BackendDAE backendDAE,backendDAE2;
1936
1937 list<BackendDAE.Var> varlst, knvarlst, states, inputvars, inputvars2, outputvars, paramvars, states_inputs, conVarsList, fconVarsList, object;
1938
1939 BackendDAE.Variables v,globalKnownVars,statesarr,inputvarsarr,paramvarsarr,outputvarsarr, optimizer_vars, conVars;
1940
1941 BackendDAE.SymbolicJacobians linearModelMatrices;
1942 Option<BackendDAE.SymbolicJacobian> linearModelMatrix;
1943
1944 BackendDAE.SparsePattern sparsePattern;
1945 BackendDAE.SparseColoring sparseColoring;
1946 BackendDAE.NonlinearPattern nonlinearPattern;
1947
1948 AvlTreePathFunction.Tree funcs, functionTree;
1949
1950
1951 case (backendDAE, false)
1952 algorithm
1953 17 backendDAE2 := BackendDAEUtil.copyBackendDAE(backendDAE);
1954 17 backendDAE2 := BackendDAEOptimize.collapseIndependentBlocks(backendDAE2);
1955 17 backendDAE2 := BackendDAEUtil.transformBackendDAE(backendDAE2,SOME((BackendDAE.NO_INDEX_REDUCTION(),BackendDAE.EXACT())),NONE(),NONE());
1956
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 17 times.
 17 BackendDAE.DAE({BackendDAE.EQSYSTEM(orderedVars = v)}, BackendDAE.SHARED(globalKnownVars = globalKnownVars)) := backendDAE2;
1957
1958 // Prepare all needed variables
1959 17 varlst := BackendVariable.varList(v);
1960 17 knvarlst := BackendVariable.varList(globalKnownVars);
1961 17 states := BackendVariable.getAllStateVarFromVariables(v);
1962 17 inputvars := List.select(knvarlst,BackendVariable.isInput);
1963 17 paramvars := List.select(knvarlst, BackendVariable.isParam);
1964 17 inputvars2 := List.select(knvarlst,BackendVariable.isVarOnTopLevelAndInput);
1965 17 outputvars := List.select(varlst, BackendVariable.isVarOnTopLevelAndOutput);
1966
1967 17 statesarr := BackendVariable.listVar1(states);
1968 17 inputvarsarr := BackendVariable.listVar1(inputvars);
1969 17 paramvarsarr := BackendVariable.listVar1(paramvars);
1970 17 outputvarsarr := BackendVariable.listVar1(outputvars);
1971
1972 // Differentiate the System w.r.t states for matrices A
1973 17 (linearModelMatrix, functionTree, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2,states,statesarr,inputvarsarr,paramvarsarr,statesarr,varlst,"A",false);
1974 17 backendDAE2 := BackendDAEUtil.setFunctionTree(backendDAE2, functionTree);
1975 17 linearModelMatrices := {(linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern)};
1976
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 17 times.
 17 if Flags.isSet(Flags.JAC_DUMP2) then
1977 ✗ print("analytical Jacobians -> generated system for matrix A time: " + realString(clock()) + "\n");
1978 end if;
1979
1980 // Differentiate the System w.r.t inputs for matrices B
1981 17 (linearModelMatrix, funcs, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2,inputvars2,statesarr,inputvarsarr,paramvarsarr,statesarr,varlst,"B",false);
1982 17 functionTree := AvlTreePathFunction.join(functionTree, funcs);
1983 17 backendDAE2 := BackendDAEUtil.setFunctionTree(backendDAE2, functionTree);
1984 17 linearModelMatrices := (linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern) :: linearModelMatrices;
1985
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 17 times.
 17 if Flags.isSet(Flags.JAC_DUMP2) then
1986 ✗ print("analytical Jacobians -> generated system for matrix B time: " + realString(clock()) + "\n");
1987 end if;
1988
1989 // Differentiate the System w.r.t states for matrices C
1990 17 (linearModelMatrix, funcs, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2,states,statesarr,inputvarsarr,paramvarsarr,outputvarsarr,varlst,"C",false);
1991 17 functionTree := AvlTreePathFunction.join(functionTree, funcs);
1992 17 backendDAE2 := BackendDAEUtil.setFunctionTree(backendDAE2, functionTree);
1993 17 linearModelMatrices := (linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern) :: linearModelMatrices;
1994
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 17 times.
 17 if Flags.isSet(Flags.JAC_DUMP2) then
1995 ✗ print("analytical Jacobians -> generated system for matrix C time: " + realString(clock()) + "\n");
1996 end if;
1997
1998 // Differentiate the System w.r.t inputs for matrices D
1999 17 (linearModelMatrix, funcs, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2,inputvars2,statesarr,inputvarsarr,paramvarsarr,outputvarsarr,varlst,"D",false);
2000 17 functionTree := AvlTreePathFunction.join(functionTree, funcs);
2001 17 linearModelMatrices := (linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern) :: linearModelMatrices;
2002
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 17 times.
 17 if Flags.isSet(Flags.JAC_DUMP2) then
2003 ✗ print("analytical Jacobians -> generated system for matrix D time: " + realString(clock()) + "\n");
2004 end if;
2005
2006 17 then
2007 (listReverse(linearModelMatrices), functionTree);
2008
2009 case (backendDAE, true) // created linear model (matrices) for optimization
2010 algorithm
2011 // A := der(x)
2012 // B := {der(x), con(x), L(x)}
2013 // C := {der(x), con(x), L(x), M(x)}
2014 // D := {}
2015
2016 34 backendDAE2 := BackendDAEUtil.copyBackendDAE(backendDAE);
2017 34 backendDAE2 := BackendDAEOptimize.collapseIndependentBlocks(backendDAE2);
2018 34 backendDAE2 := BackendDAEUtil.transformBackendDAE(backendDAE2,SOME((BackendDAE.NO_INDEX_REDUCTION(),BackendDAE.EXACT())),NONE(),NONE());
2019
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34 times.
 34 BackendDAE.DAE({BackendDAE.EQSYSTEM(orderedVars = v)}, BackendDAE.SHARED(globalKnownVars = globalKnownVars)) := backendDAE2;
2020
2021 // Prepare all needed variables
2022 34 varlst := BackendVariable.varList(v);
2023 34 knvarlst := BackendVariable.varList(globalKnownVars);
2024 34 states := BackendVariable.getAllStateVarFromVariables(v);
2025 34 inputvars := List.select(knvarlst,BackendVariable.isInput);
2026 34 paramvars := List.select(knvarlst, BackendVariable.isParam);
2027 34 inputvars2 := List.select(knvarlst,BackendVariable.isVarOnTopLevelAndInputNoDerInput); // without der(u)
2028 34 outputvars := List.select(varlst, BackendVariable.isVarOnTopLevelAndOutput);
2029 34 conVarsList := List.select(varlst, BackendVariable.isRealOptimizeConstraintsVars);
2030 34 fconVarsList := List.select(varlst, BackendVariable.isRealOptimizeFinalConstraintsVars); // ToDo: FinalCon
2031
2032 34 states_inputs := listAppend(states, inputvars2);
2033 34 statesarr := BackendVariable.listVar1(states);
2034 34 inputvarsarr := BackendVariable.listVar1(inputvars);
2035 34 paramvarsarr := BackendVariable.listVar1(paramvars);
2036 34 outputvarsarr := BackendVariable.listVar1(outputvars);
2037 34 conVars := BackendVariable.listVar1(conVarsList);
2038
2039 //BackendDump.printVariables(conVars);
2040 //BackendDump.printVariables(object);
2041 //print(intString(BackendVariable.varsSize(object)));
2042 //object = BackendVariable.listVar1(object);
2043
2044 // Differentiate the System w.r.t states for matrices A
2045 34 (linearModelMatrix, functionTree, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2,states,statesarr,inputvarsarr,paramvarsarr,statesarr,varlst,"A",false);
2046
2047 34 backendDAE2 := BackendDAEUtil.setFunctionTree(backendDAE2, functionTree);
2048 34 linearModelMatrices := {(linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern)};
2049
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if Flags.isSet(Flags.JAC_DUMP2) then
2050 ✗ print("analytical Jacobians -> generated system for matrix A time: " + realString(clock()) + "\n");
2051 end if;
2052
2053 // Differentiate the System w.r.t states&inputs for matrices B
2054
2055 34 optimizer_vars := BackendVariable.addVariables(statesarr, BackendVariable.copyVariables(conVars));
2056 34 object := DynamicOptimization.checkObjectIsSet(outputvarsarr, BackendDAE.optimizationLagrangeTermName);
2057 34 optimizer_vars := BackendVariable.addVars(object, optimizer_vars);
2058 //BackendDump.printVariables(optimizer_vars);
2059 34 (linearModelMatrix, funcs, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2,states_inputs,statesarr,inputvarsarr,paramvarsarr,optimizer_vars,varlst,"B",false);
2060 34 functionTree := AvlTreePathFunction.join(functionTree, funcs);
2061 34 backendDAE2 := BackendDAEUtil.setFunctionTree(backendDAE2, functionTree);
2062 34 linearModelMatrices := (linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern) :: linearModelMatrices;
2063
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if Flags.isSet(Flags.JAC_DUMP2) then
2064 ✗ print("analytical Jacobians -> generated system for matrix B time: " + realString(clock()) + "\n");
2065 end if;
2066
2067 // Differentiate the System w.r.t states for matrices C
2068 34 object := DynamicOptimization.checkObjectIsSet(outputvarsarr, BackendDAE.optimizationMayerTermName);
2069 34 optimizer_vars := BackendVariable.addVars(object, optimizer_vars);
2070 //BackendDump.printVariables(optimizer_vars);
2071 34 (linearModelMatrix, funcs, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2,states_inputs,statesarr,inputvarsarr,paramvarsarr,optimizer_vars,varlst,"C",false);
2072 34 functionTree := AvlTreePathFunction.join(functionTree, funcs);
2073 34 backendDAE2 := BackendDAEUtil.setFunctionTree(backendDAE2, functionTree);
2074 34 linearModelMatrices := (linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern) :: linearModelMatrices;
2075
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if Flags.isSet(Flags.JAC_DUMP2) then
2076 ✗ print("analytical Jacobians -> generated system for matrix C time: " + realString(clock()) + "\n");
2077 end if;
2078
2079 // Differentiate the System w.r.t inputs for matrices D
2080 34 optimizer_vars := BackendVariable.emptyVars();
2081 34 optimizer_vars := BackendVariable.listVar1(fconVarsList);
2082
2083 34 (linearModelMatrix, funcs, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE2, states_inputs, statesarr, inputvarsarr, paramvarsarr, optimizer_vars, varlst, "D", false);
2084 34 functionTree := AvlTreePathFunction.join(functionTree, funcs);
2085 34 linearModelMatrices := (linearModelMatrix,sparsePattern,sparseColoring, nonlinearPattern) :: linearModelMatrices;
2086
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if Flags.isSet(Flags.JAC_DUMP2) then
2087 ✗ print("analytical Jacobians -> generated system for matrix D time: " + realString(clock()) + "\n");
2088 end if;
2089
2090 34 then
2091 (listReverse(linearModelMatrices), functionTree);
2092 else
2093 algorithm
2094 ✗ Error.addInternalError("Generation of LinearModel Matrices failed.", sourceInfo());
2095 ✗ then
2096 fail();
2097 end match;
2098 end createLinearModelMatrices;
2099
2100 protected function generateGenericJacobian "author: wbraun"
2101 input BackendDAE.BackendDAE inBackendDAE;
2102 input list<BackendDAE.Var> inDiffVars "independent vars";
2103 input BackendDAE.Variables inStateVars;
2104 input BackendDAE.Variables inInputVars;
2105 input BackendDAE.Variables inParameterVars "globalKnownVars";
2106 input BackendDAE.Variables inDifferentiatedVars "resVars";
2107 input list<BackendDAE.Var> inVars "dependent vars = resVars + other vars";
2108 input String inName;
2109 input Boolean onlySparsePattern;
2110 input Boolean daeMode = false;
2111 output Option<BackendDAE.SymbolicJacobian> outJacobian;
2112 output AvlTreePathFunction.Tree outFunctionTree;
2113 output BackendDAE.SparsePattern outSparsePattern = BackendDAE.emptySparsePattern;
2114 output BackendDAE.SparseColoring outSparseColoring = {};
2115 output BackendDAE.NonlinearPattern nonlinearPattern;
2116 protected
2117 BackendDAE.SymbolicJacobian symbolicJacobian;
2118 BackendDAE.Shared shared = inBackendDAE.shared;
2119 BackendDAE.BackendDAE jacDAE;
2120 list<BackendDAE.Var> jacDiffedVars;
2121 algorithm
2122 try
2123 2528 outFunctionTree := shared.functionTree;
2124
2/2
✓ Branch 0 taken 1820 times.
✓ Branch 1 taken 708 times.
 2528 if not onlySparsePattern then
2125 1820 (symbolicJacobian, outFunctionTree) := createJacobian(inBackendDAE,inDiffVars, inStateVars, inInputVars, inParameterVars, inDifferentiatedVars, inVars, inName, daeMode);
2126
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1819 times.
 1820 true := checkForNonLinearStrongComponents(symbolicJacobian);
2127 outJacobian := SOME(symbolicJacobian);
2128 // nonlinear pattern is the same as the sparse pattern of the jacobian
2129 1819 (jacDAE, _, _, _, _, _) := symbolicJacobian;
2130 1819 jacDiffedVars := getJacobianResiduals(jacDAE);
2131 // copy the jacobian DAE to avoid wrong variables being added
2132 1819 (nonlinearPattern, _) := generateSparsePattern(BackendDAEUtil.copyBackendDAE(jacDAE), inDiffVars, jacDiffedVars, true);
2133 1819 nonlinearPattern := stripPartialDerNonlinearPattern(nonlinearPattern);
2134 else
2135 outJacobian := NONE();
2136 // no jacobian -> no nonlinear pattern
2137 nonlinearPattern := BackendDAE.emptyNonlinearPattern;
2138 end if;
2139 // generate sparse pattern
2140
3/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2515 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 12 times.
 2527 if (not stringEq(inName, "FMIDERINIT")) then
2141 2515 (outSparsePattern,outSparseColoring) := generateSparsePattern(inBackendDAE, inDiffVars, BackendVariable.varList(inDifferentiatedVars));
2142 end if;
2143 else
2144 1 fail();
2145 end try;
2146 end generateGenericJacobian;
2147
2148 protected function createJacobian "author: wbraun"
2149 input BackendDAE.BackendDAE inBackendDAE;
2150 input list<BackendDAE.Var> inDiffVars "independent vars";
2151 input BackendDAE.Variables inStateVars;
2152 input BackendDAE.Variables inInputVars;
2153 input BackendDAE.Variables inParameterVars "globalKnownVars";
2154 input BackendDAE.Variables inDifferentiatedVars "resVars";
2155 input list<BackendDAE.Var> inVars "dependent vars = resVars + other vars";
2156 input String inName;
2157 input Boolean daeMode;
2158 output BackendDAE.SymbolicJacobian outJacobian;
2159 output AvlTreePathFunction.Tree outFunctionTree;
2160 algorithm
2161 (outJacobian, outFunctionTree) :=
2162 matchcontinue inName
2163 local
2164 BackendDAE.BackendDAE backendDAE, reducedDAE;
2165
2166 list<DAE.ComponentRef> comref_vars, comref_differentiatedVars, dependencies;
2167
2168 BackendDAE.Shared shared;
2169 BackendDAE.Variables globalKnownVars;
2170 list<BackendDAE.Var> diffedVars "resVars", seedlst, indepVars;
2171
2172 AvlTreePathFunction.Tree funcs;
2173
2174 case _
2175 algorithm
2176 1820 diffedVars := BackendVariable.varList(inDifferentiatedVars);
2177 1820 comref_differentiatedVars := List.map(diffedVars, BackendVariable.varCref);
2178
2179 1820 reducedDAE := BackendDAEUtil.reduceEqSystemsInDAE(inBackendDAE, diffedVars, true, not Flags.getConfigBool(Flags.CAUSALIZE_DAE_MODE));
2180
2181 1820 indepVars := createInDepVars(inDiffVars, false);
2182 1820 comref_vars := List.map(inDiffVars, BackendVariable.varCref);
2183 1820 seedlst := List.map1(comref_vars, createSeedVars, inName);
2184
2185
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1820 times.
 1820 if Flags.isSet(Flags.JAC_DUMP) then
2186 ✗ print("Create symbolic Jacobians from:\n");
2187 ✗ print(BackendDump.varListString(indepVars, "Independent Variables"));
2188 ✗ print(BackendDump.varListString(diffedVars, "Dependent Variables"));
2189 ✗ print("Basic equation system:\n");
2190 ✗ print(BackendDump.equationListString(BackendEquation.equationSystemsEqnsLst(reducedDAE.eqs), "differentiated equations"));
2191 ✗ print(BackendDump.varListString(BackendVariable.equationSystemsVarsLst(reducedDAE.eqs), "related variables"));
2192 ✗ print(BackendDump.varListString(BackendVariable.varList(reducedDAE.shared.globalKnownVars), "known variables"));
2193 end if;
2194
2195 // Differentiate the eqns system in reducedDAE w.r.t. independents
2196 1820 (backendDAE as BackendDAE.DAE(), funcs) := generateSymbolicJacobian(reducedDAE, indepVars, inDifferentiatedVars, BackendVariable.listVar1(seedlst), inStateVars, inInputVars, inParameterVars, inName, daeMode);
2197
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1820 times.
 1820 if Flags.isSet(Flags.JAC_DUMP2) then
2198 ✗ print("analytical Jacobians -> generated equations for Jacobian " + inName + " time: " + realString(clock()) + "\n");
2199 end if;
2200
2201 // Add the function tree to the jacobian backendDAE
2202 1820 backendDAE := BackendDAEUtil.setFunctionTree(backendDAE, funcs);
2203
2204 1820 backendDAE := optimizeJacobianMatrix(backendDAE,comref_differentiatedVars,comref_vars);
2205
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1820 times.
 1820 if Flags.isSet(Flags.JAC_DUMP2) then
2206 ✗ print("analytical Jacobians -> generated Jacobian DAE time: " + realString(clock()) + "\n");
2207 end if;
2208 1820 dependencies := calcJacobianDependencies((backendDAE, "", {}, {}, {}, {}));
2209
2210 1820 then
2211 ((backendDAE, inName, inDiffVars, diffedVars, inVars, dependencies), funcs);
2212 else
2213 algorithm
2214 ✗ Error.addInternalError("function createJacobian failed", sourceInfo());
2215 ✗ then
2216 fail();
2217 end matchcontinue;
2218 end createJacobian;
2219
2220 protected function optimizeJacobianMatrix "author: wbraun"
2221 input BackendDAE.BackendDAE inBackendDAE;
2222 input list<DAE.ComponentRef> inComRef1 "eqnvars";
2223 input list<DAE.ComponentRef> inComRef2 "vars to differentiate";
2224 output BackendDAE.BackendDAE outJacobian;
2225 protected
2226 array<Integer> ea = listArray({});
2227 BackendDAE.Matching eMatching = BackendDAE.MATCHING(ea, ea, {});
2228 algorithm
2229 outJacobian :=
2230 matchcontinue (inBackendDAE,inComRef1,inComRef2)
2231 local
2232 BackendDAE.BackendDAE backendDAE, backendDAE2;
2233 BackendDAE.EqSystem syst;
2234 BackendDAE.Shared shared;
2235 Boolean b = false;
2236 list<String> strPostOptModules;
2237
2238 case (BackendDAE.DAE(syst::{}, shared), {}, _)
2239 algorithm
2240 56 syst.orderedVars := BackendVariable.listVar({});
2241 56 syst.matching := eMatching;
2242 56 then BackendDAE.DAE(syst::{}, shared);
2243 case (BackendDAE.DAE(syst::{}, shared), _, {})
2244 algorithm
2245 22 syst.orderedVars := BackendVariable.listVar({});
2246 22 syst.matching := eMatching;
2247 22 then BackendDAE.DAE(syst::{}, shared);
2248 case (backendDAE, _, _)
2249 algorithm
2250
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1742 times.
 1742 if Flags.isSet(Flags.JAC_DUMP2) then
2251 ✗ print("analytical Jacobians -> optimize jacobians time: " + realString(clock()) + "\n");
2252 end if;
2253
2254
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1742 times.
 1742 if Flags.isSet(Flags.JAC_DUMP) then
2255 ✗ BackendDump.bltdump("Symbolic Jacobian",backendDAE);
2256 else
2257 1742 b := FlagsUtil.disableDebug(Flags.EXEC_STAT);
2258 end if;
2259
2260 strPostOptModules := {"wrapFunctionCalls",
2261 "inlineArrayEqn",
2262 "constantLinearSystem",
2263 "solveSimpleEquations",
2264 "tearingSystem",
2265 "calculateStrongComponentJacobians",
2266 "removeConstants",
2267 "simplifyTimeIndepFuncCalls"};
2268
2269 // Add removeSimpleEquation to remove constant(= independent of seed) equations.
2270
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1742 times.
 1742 if Flags.isSet(Flags.SPLIT_CONSTANT_PARTS_SYMJAC) then
2271 /* ToDo: removeSimpleEquation can't handle all sorts of equations inside the
2272 * jacobian BackendDAE. E.g. for equations lile
2273 * $cse14 := $DER$$PModelica$PMedia$PWater$PIF97_Utilities$PwaterBaseProp_ph(p[10], h[10], 0, 0, 1.0, 0.0);
2274 * from SteamPipe from ScalableTestsuite.
2275 * Add a new module which finds constant (= independent of seed) equations
2276 * and moves them to a different system.
2277 */
2278 ✗ strPostOptModules := List.insert(strPostOptModules, 4, "removeSimpleEquations");
2279 end if;
2280
2281 1742 backendDAE2 := BackendDAEUtil.getSolvedSystemforJacobians(backendDAE,
2282 {"removeEqualRHS",
2283 "removeSimpleEquations",
2284 "evalFunc"},
2285 NONE(),
2286 NONE(),
2287 strPostOptModules);
2288
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1742 times.
 1742 if Flags.isSet(Flags.JAC_DUMP) then
2289 ✗ BackendDump.bltdump("Symbolic Jacobian",backendDAE2);
2290 else
2291 1742 FlagsUtil.set(Flags.EXEC_STAT, b);
2292 end if;
2293 then backendDAE2;
2294 else
2295 algorithm
2296 ✗ Error.addInternalError("function optimizeJacobianMatrix failed", sourceInfo());
2297 ✗ then fail();
2298 end matchcontinue;
2299 end optimizeJacobianMatrix;
2300
2301 protected function generateSymbolicJacobian "author: lochel"
2302 input BackendDAE.BackendDAE inBackendDAE "reducedDAE (variables and equations needed to calculate resVars)";
2303 input list<BackendDAE.Var> inVars "independent vars";
2304 input BackendDAE.Variables inDiffedVars "resVars";
2305 input BackendDAE.Variables inSeedVars;
2306 input BackendDAE.Variables inStateVars;
2307 input BackendDAE.Variables inInputVars;
2308 input BackendDAE.Variables inParamVars "globalKnownVars";
2309 input String inMatrixName;
2310 input Boolean daeMode;
2311 output BackendDAE.BackendDAE outJacobian;
2312 output AvlTreePathFunction.Tree outFunctions;
2313 algorithm
2314 (outJacobian,outFunctions) := matchcontinue(inBackendDAE, inVars, inDiffedVars, inMatrixName)
2315 local
2316 AvlTreePathFunction.Tree functions;
2317 list<DAE.ComponentRef> comref_diffvars;
2318 DAE.ComponentRef x;
2319 String dummyVarName;
2320
2321 BackendDAE.Variables diffVarsArr;
2322 BackendDAE.Variables diffedVars "resVars";
2323 BackendDAE.BackendDAE jacobian;
2324
2325 // BackendDAE
2326 BackendDAE.Variables orderedVars, jacOrderedVars; // ordered Variables, only states and alg. vars
2327 BackendDAE.Variables globalKnownVars, jacKnownVars; // Known variables, i.e. constants and parameters
2328 BackendDAE.EquationArray orderedEqs, jacOrderedEqs; // ordered Equations
2329 // Removed equations a=b
2330 // end BackendDAE
2331
2332 list<BackendDAE.Var> diffVars "independent vars", derivedVariables;
2333 list<BackendDAE.Equation> eqns, derivedEquations;
2334
2335
2336
2337 FCore.Cache cache;
2338 FCore.Graph graph;
2339 BackendDAE.Shared shared;
2340
2341 String matrixName;
2342 array<Integer> ass2;
2343
2344 BackendDAE.DifferentiateInputData diffData;
2345
2346 BackendDAE.ExtraInfo ei;
2347 Integer size;
2348
2349 case(BackendDAE.DAE(shared=BackendDAE.SHARED(cache=cache, graph=graph, info=ei, functionTree=functions)), {}, _, _) algorithm
2350 74 jacobian := BackendDAE.DAE( {BackendDAEUtil.createEqSystem(BackendVariable.emptyVars(), BackendEquation.emptyEqns())},
2351 BackendDAEUtil.createEmptyShared(BackendDAE.JACOBIAN(), ei, cache, graph));
2352 37 then (jacobian, functions);
2353
2354 case(BackendDAE.DAE( BackendDAE.EQSYSTEM(orderedVars=orderedVars, orderedEqs=orderedEqs, matching=BackendDAE.MATCHING(ass2=ass2))::{},
2355 BackendDAE.SHARED(globalKnownVars=globalKnownVars, cache=cache,graph=graph, functionTree=functions, info=ei) ), diffVars, diffedVars, matrixName) algorithm
2356 // Generate tmp variables
2357 1783 dummyVarName := ("dummyVar" + matrixName);
2358 1783 x := DAE.CREF_IDENT(dummyVarName,DAE.T_REAL_DEFAULT,{});
2359
2360 // differentiate the equation system
2361
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1783 times.
 1783 if Flags.isSet(Flags.JAC_DUMP2) then
2362 ✗ print("*** analytical Jacobians -> derived all algorithms time: " + realString(clock()) + "\n");
2363 end if;
2364 1783 diffVarsArr := BackendVariable.listVar1(diffVars);
2365 1783 comref_diffvars := List.map(diffVars, BackendVariable.varCref);
2366 diffData := BackendDAE.emptyInputData;
2367 1783 diffData.independenentVars := SOME(diffVarsArr);
2368 1783 diffData.dependenentVars := SOME(diffedVars);
2369 1783 diffData.knownVars := SOME(globalKnownVars);
2370 1783 diffData.allVars := SOME(orderedVars);
2371 1783 diffData.diffCrefs := comref_diffvars;
2372 1783 diffData.matrixName := SOME(matrixName);
2373 1783 eqns := BackendEquation.equationList(orderedEqs);
2374
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1783 times.
 1783 if Flags.isSet(Flags.JAC_DUMP2) then
2375 ✗ print("*** analytical Jacobians -> before derive all equation: " + realString(clock()) + "\n");
2376 end if;
2377 1783 (derivedEquations, functions) := deriveAll(eqns, arrayList(ass2), x, diffData, functions, daeMode);
2378
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1783 times.
 1783 if Flags.isSet(Flags.JAC_DUMP2) then
2379 ✗ print("*** analytical Jacobians -> after derive all equation: " + realString(clock()) + "\n");
2380 end if;
2381 // replace all der(x), since ExpressionSolve can't handle der(x) proper
2382 1783 derivedEquations := BackendEquation.replaceDerOpInEquationList(derivedEquations);
2383
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1783 times.
 1783 if Flags.isSet(Flags.JAC_DUMP2) then
2384 ✗ print("*** analytical Jacobians -> created all derived equation time: " + realString(clock()) + "\n");
2385 end if;
2386
2387 // create BackendDAE.DAE with differentiated vars and equations
2388
2389 // all variables for new equation system
2390 // d(ordered vars)/d(dummyVar)
2391 1783 diffVars := BackendVariable.varList(orderedVars);
2392 1783 derivedVariables := createAllDiffedVars(diffVars, x, diffedVars, matrixName);
2393
2394 1783 jacOrderedVars := BackendVariable.listVar1(derivedVariables);
2395 // known vars: all variable from original system + seed
2396 1783 size := BackendVariable.varsSize(orderedVars) +
2397 BackendVariable.varsSize(globalKnownVars) +
2398 BackendVariable.varsSize(inSeedVars);
2399 1783 jacKnownVars := BackendVariable.emptyVarsSized(size);
2400 1783 jacKnownVars := BackendVariable.addVariables(inSeedVars, jacKnownVars);
2401 1783 (jacKnownVars,_) := BackendVariable.traverseBackendDAEVarsWithUpdate(jacKnownVars, BackendVariable.setVarDirectionTpl, (DAE.INPUT()));
2402 1783 jacKnownVars := BackendVariable.addVariables(orderedVars, jacKnownVars);
2403 1783 jacKnownVars := BackendVariable.addVariables(globalKnownVars, jacKnownVars);
2404 1783 jacOrderedEqs := BackendEquation.listEquation(derivedEquations);
2405
2406
2407 1783 shared := BackendDAEUtil.createEmptyShared(BackendDAE.JACOBIAN(), ei, cache, graph);
2408
2409 3566 jacobian := BackendDAE.DAE( BackendDAEUtil.createEqSystem(jacOrderedVars, jacOrderedEqs)::{},
2410 BackendDAEUtil.setSharedGlobalKnownVars(shared, jacKnownVars) );
2411 1783 then (jacobian, functions);
2412
2413 else
2414 algorithm
2415 ✗ Error.addInternalError(getInstanceName() + " failed", sourceInfo());
2416 ✗ then fail();
2417 end matchcontinue;
2418 end generateSymbolicJacobian;
2419
2420 public function createSeedVars "author: wbraun"
2421 input DAE.ComponentRef indiffVar;
2422 input String inMatrixName;
2423 output BackendDAE.Var outSeedVar;
2424 protected
2425 DAE.ComponentRef derivedCref;
2426 algorithm
2427 6011 derivedCref := Differentiate.createSeedCrefName(indiffVar, inMatrixName);
2428 6011 outSeedVar := BackendDAE.VAR(derivedCref, BackendDAE.STATE_DER(), DAE.INPUT(), DAE.NON_PARALLEL(), ComponentReference.crefLastType(derivedCref), NONE(), NONE(), {}, DAE.emptyElementSource, NONE(), NONE(), NONE(), NONE(),DAE.NON_CONNECTOR(), DAE.NOT_INNER_OUTER(), true, false, false);
2429 end createSeedVars;
2430
2431 protected function createAllDiffedVars "author: wbraun"
2432 input list<BackendDAE.Var> inVars;
2433 input DAE.ComponentRef inCref;
2434 input BackendDAE.Variables inAllVars;
2435 input String inMatrixName;
2436 output list<BackendDAE.Var> outVars;
2437 algorithm
2438 try
2439 1783 outVars := createAllDiffedVarsWork(inVars, inCref, inAllVars, 0, inMatrixName, {});
2440 else
2441 ✗ Error.addMessage(Error.INTERNAL_ERROR, {"SymbolicJacobian.createAllDiffedVars failed"});
2442 ✗ fail();
2443 end try;
2444 end createAllDiffedVars;
2445
2446 protected function createAllDiffedVarsWork "author: wbraun,hkiel"
2447 input list<BackendDAE.Var> inVars;
2448 input DAE.ComponentRef inCref;
2449 input BackendDAE.Variables inAllVars;
2450 input Integer inIndex;
2451 input String inMatrixName;
2452 input list<BackendDAE.Var> iVars;
2453 output list<BackendDAE.Var> outVars;
2454 algorithm
2455 outVars := match(inVars, inCref, inIndex)
2456 local
2457 BackendDAE.Var v, r1;
2458 DAE.ComponentRef currVar, cref, derivedCref;
2459 list<BackendDAE.Var> restVar;
2460 Integer index;
2461
2462 case({}, _, _)
2463 1783 then listReverse(iVars);
2464
2465 case((v as BackendDAE.VAR(varName=currVar,varKind=BackendDAE.STATE()))::restVar, cref, index) algorithm
2466 try
2467 452 BackendVariable.getVarSingle(currVar, inAllVars);
2468 334 currVar := ComponentReference.crefPrefixDer(currVar);
2469 334 derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName);
2470 334 r1 := BackendVariable.copyVarNewName(derivedCref, v);
2471 334 r1 := BackendVariable.setVarKind(r1, BackendDAE.STATE_DER());
2472 334 r1.unreplaceable := true;
2473 index := index + 1;
2474 else
2475 118 currVar := ComponentReference.crefPrefixDer(currVar);
2476 118 derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName);
2477 118 r1 := BackendVariable.copyVarNewName(derivedCref, v);
2478 118 r1 := BackendVariable.setVarKind(r1, BackendDAE.STATE_DER());
2479 end try;
2480 452 then
2481 createAllDiffedVarsWork(restVar, cref, inAllVars, index, inMatrixName, r1::iVars);
2482
2483 case((v as BackendDAE.VAR(varName=currVar))::restVar, cref, index) algorithm
2484 try
2485 18054 BackendVariable.getVarSingle(currVar, inAllVars);
2486 5423 derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName);
2487 5423 r1 := BackendVariable.copyVarNewName(derivedCref, v);
2488 5423 r1 := BackendVariable.setVarKind(r1, BackendDAE.VARIABLE());
2489 5423 r1.unreplaceable := true;
2490 index := index + 1;
2491 else
2492 12631 derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName);
2493 12631 r1 := BackendVariable.copyVarNewName(derivedCref, v);
2494 12631 r1 := BackendVariable.setVarKind(r1, BackendDAE.VARIABLE());
2495 end try;
2496 18054 then
2497 createAllDiffedVarsWork(restVar, cref, inAllVars, index, inMatrixName, r1::iVars);
2498
2499 end match;
2500 end createAllDiffedVarsWork;
2501
2502 protected function deriveAll
2503 input list<BackendDAE.Equation> inEquations;
2504 input list<Integer> ass2;
2505 input DAE.ComponentRef inDiffCref;
2506 input BackendDAE.DifferentiateInputData inDiffData;
2507 input AvlTreePathFunction.Tree inFunctions;
2508 input Boolean daeMode;
2509 output list<BackendDAE.Equation> outDerivedEquations = {};
2510 output AvlTreePathFunction.Tree outFunctions = inFunctions;
2511 protected
2512 BackendDAE.Variables allVars;
2513 BackendDAE.Equation currDerivedEquation;
2514 list<BackendDAE.Equation> tmpEquations;
2515 algorithm
2516 try
2517
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1783 times.
✓ Branch 2 taken 1783 times.
✗ Branch 3 not taken.
 1783 BackendDAE.DIFFINPUTDATA(allVars=SOME(allVars)) := inDiffData;
2518
2/2
✓ Branch 0 taken 18135 times.
✓ Branch 1 taken 1783 times.
 19918 for currEquation in inEquations loop
2519
1/2
✓ Branch 0 taken 18135 times.
✗ Branch 1 not taken.
 36270 (currDerivedEquation, outFunctions) := Differentiate.differentiateEquation(currEquation, inDiffCref, inDiffData, BackendDAE.GENERIC_GRADIENT(daeMode), outFunctions);
2520 18135 tmpEquations := BackendEquation.scalarComplexEquations(currDerivedEquation, outFunctions);
2521 18135 outDerivedEquations := listAppend(tmpEquations, outDerivedEquations);
2522 end for;
2523
2524 1783 outDerivedEquations := listReverse(outDerivedEquations);
2525
2526 else
2527 ✗ Error.addMessage(Error.INTERNAL_ERROR, {"SymbolicJacobian.deriveAll failed"});
2528 ✗ fail();
2529 end try;
2530 end deriveAll;
2531
2532 public function getJacobianMatrixbyName
2533 input BackendDAE.SymbolicJacobians injacobianMatrices;
2534 input String inJacobianName;
2535 output Option<tuple<Option<BackendDAE.SymbolicJacobian>, BackendDAE.SparsePattern, BackendDAE.SparseColoring, BackendDAE.NonlinearPattern>> outMatrix;
2536 algorithm
2537 outMatrix := match injacobianMatrices
2538 local
2539 tuple<Option<BackendDAE.SymbolicJacobian>, BackendDAE.SparsePattern, BackendDAE.SparseColoring, BackendDAE.NonlinearPattern> matrix;
2540 BackendDAE.SymbolicJacobians rest;
2541 String name;
2542
2543 case (matrix as (SOME((_,name,_,_,_,_)), _, _, _))::_ guard
2544 stringEq(name, inJacobianName)
2545 then SOME(matrix);
2546
2547 case _::rest
2548 ✗ then getJacobianMatrixbyName(rest, inJacobianName);
2549
2550 else NONE();
2551 end match;
2552 end getJacobianMatrixbyName;
2553
2554 public function updateJacobianDependencies
2555 input output BackendDAE.Jacobian jacobian;
2556 algorithm
2557 jacobian := match jacobian
2558 local
2559 BackendDAE.Jacobian jac;
2560 BackendDAE.SymbolicJacobian symJac;
2561 BackendDAE.EqSystem syst;
2562 BackendDAE.Shared shared;
2563 String name;
2564 list<BackendDAE.Var> diffVars;
2565 list<BackendDAE.Var> diffedVars;
2566 list<BackendDAE.Var> allDiffedVars;
2567 list<DAE.ComponentRef> dependencies;
2568 case jac as BackendDAE.GENERIC_JACOBIAN()
2569 algorithm
2570 ✗ SOME(symJac as (BackendDAE.DAE({syst}, shared),name,diffVars,diffedVars,allDiffedVars,dependencies)) := jac.jacobian;
2571 ✗ dependencies := calcJacobianDependencies(symJac);
2572 ✗ jac.jacobian := SOME((BackendDAE.DAE({syst}, shared),name,diffVars,diffedVars,allDiffedVars,dependencies));
2573 then jac;
2574 else jacobian;
2575 end match;
2576 end updateJacobianDependencies;
2577
2578 public function calcJacobianDependencies
2579 input BackendDAE.SymbolicJacobian jacobian;
2580 output list<DAE.ComponentRef> dependencies;
2581 protected
2582 BackendDAE.EqSystems systems;
2583 BackendDAE.Shared shared;
2584 BackendDAE.EqSystem syst;
2585 algorithm
2586 1820 (BackendDAE.DAE(systems, shared), _, _, _, _, _) := jacobian;
2587 1820 syst := listHead(systems); // Only the first system contains directional derivative,
2588 // the others contain optional constant equations
2589 1820 dependencies := BackendEquation.getCrefsFromEquations(syst.orderedEqs, syst.orderedVars, shared.globalKnownVars);
2590 end calcJacobianDependencies;
2591
2592 public function getJacobianDependencies
2593 input BackendDAE.Jacobian jacobian;
2594 output list<DAE.ComponentRef> dependencies;
2595 algorithm
2596 dependencies := match jacobian
2597 case BackendDAE.GENERIC_JACOBIAN(jacobian=SOME((_, _, _, _, _, dependencies)))
2598 then dependencies;
2599
2600 case BackendDAE.GENERIC_JACOBIAN(jacobian=NONE())
2601 then {};
2602
2603 else algorithm
2604 ✗ Error.addInternalError("function getJacobianDependencies failed", sourceInfo());
2605 ✗ then fail();
2606
2607 end match;
2608 end getJacobianDependencies;
2609
2610 // =============================================================================
2611 // Module for to calculate strong component Jacobains
2612 //
2613 // =============================================================================
2614
2615 protected function calculateEqSystemJacobians
2616 input BackendDAE.EqSystem inSyst;
2617 input BackendDAE.Shared inShared;
2618 output BackendDAE.EqSystem outSyst;
2619 output BackendDAE.Shared outShared;
2620 algorithm
2621 (outSyst, outShared) := match (inSyst, inShared)
2622 local
2623 BackendDAE.EqSystem syst;
2624 BackendDAE.Shared shared;
2625 array<Integer> ass1;
2626 array<Integer> ass2;
2627 BackendDAE.StrongComponents comps;
2628 BackendDAE.Variables vars;
2629 BackendDAE.EquationArray eqns;
2630
2631 case (syst as BackendDAE.EQSYSTEM( orderedVars=vars, orderedEqs=eqns,
2632 matching=BackendDAE.MATCHING(ass1,ass2,comps) ), shared)
2633 algorithm
2634 53242 (comps, shared) := calculateJacobiansComponents(comps, vars, eqns, shared);
2635
1/2
✓ Branch 1 taken 53242 times.
✗ Branch 2 not taken.
 106484 syst.matching := BackendDAE.MATCHING(ass1, ass2, comps);
2636
1/2
✓ Branch 0 taken 53242 times.
✗ Branch 1 not taken.
 53242 then (syst, shared);
2637 end match;
2638 end calculateEqSystemJacobians;
2639
2640 protected function calculateJacobiansComponents
2641 input BackendDAE.StrongComponents inComps;
2642 input BackendDAE.Variables inVars;
2643 input BackendDAE.EquationArray inEqns;
2644 input BackendDAE.Shared inShared;
2645 output BackendDAE.StrongComponents outComps;
2646 output BackendDAE.Shared outShared = inShared;
2647 algorithm
2648
4/4
✓ Branch 0 taken 160546 times.
✓ Branch 1 taken 53242 times.
✓ Branch 2 taken 160546 times.
✓ Branch 3 taken 53242 times.
 213788 outComps := list(match component
2649 local
2650 BackendDAE.StrongComponent comp;
2651 case comp algorithm
2652 160546 (comp, outShared) := calculateJacobianComponent(comp, inVars, inEqns, outShared);
2653 then comp;
2654 end match for component in inComps);
2655 end calculateJacobiansComponents;
2656
2657 public function prepareTornStrongComponentData
2658 input BackendDAE.Variables inVars;
2659 input BackendDAE.EquationArray inEqns;
2660 input list<Integer> inIterationvarsInts;
2661 input list<Integer> inResidualequations;
2662 input BackendDAE.InnerEquations innerEquations;
2663 input AvlTreePathFunction.Tree funcTree;
2664 input String name;
2665 output BackendDAE.Variables outDiffVars;
2666 output BackendDAE.Variables outResidualVars;
2667 output BackendDAE.Variables outOtherVars;
2668 output BackendDAE.EquationArray outResidualEqns;
2669 output BackendDAE.EquationArray outOtherEqns;
2670 protected
2671 list<BackendDAE.Var> iterationvars, resVarsLst, ovarsLst;
2672 list<BackendDAE.Equation> reqns, otherEqnsLst;
2673 list<list<Integer>> otherVarsIntsLst;
2674 list<Integer> otherEqnsInts, otherVarsInts;
2675 algorithm
2676 try
2677 // get iteration vars
2678
4/4
✓ Branch 0 taken 6905 times.
✓ Branch 1 taken 1823 times.
✓ Branch 2 taken 6905 times.
✓ Branch 3 taken 1823 times.
 8728 iterationvars := list(BackendVariable.transformXToXd(BackendVariable.getVarAt(inVars, e)) for e in inIterationvarsInts);
2679 1823 outDiffVars := BackendVariable.listVar1(iterationvars);
2680
2681 // debug
2682
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1823 times.
 1823 if Flags.isSet(Flags.DEBUG_ALGLOOP_JACOBIAN) then
2683 ✗ print("*** got iteration variables at time: " + realString(clock()) + "\n");
2684 ✗ BackendDump.printVarList(iterationvars);
2685 end if;
2686
2687 // get residual eqns
2688 1823 reqns := BackendEquation.getList(inResidualequations, inEqns);
2689 1823 reqns := BackendEquation.replaceDerOpInEquationList(reqns);
2690 1823 outResidualEqns := BackendEquation.listEquation(reqns);
2691
2692 // create residual equations
2693 1823 (_, reqns) := BackendEquation.traverseEquationArray(outResidualEqns, BackendEquation.traverseEquationToScalarResidualForm, (funcTree, {}));
2694 1821 reqns := listReverse(reqns);
2695 1821 (reqns, resVarsLst) := BackendEquation.convertResidualsIntoSolvedEquations(reqns, "$res_" + name + "_", 1);
2696 1811 outResidualVars := BackendVariable.listVar1(resVarsLst);
2697 1811 outResidualEqns := BackendEquation.listEquation(reqns);
2698
2699 // debug
2700
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1811 times.
 1811 if Flags.isSet(Flags.DEBUG_ALGLOOP_JACOBIAN) then
2701 ✗ print("*** got residual equation and created corresponding variables at time: " + realString(clock()) + "\n");
2702 ✗ print("Equations:\n");
2703 ✗ BackendDump.printEquationList(reqns);
2704 end if;
2705
2706 // get other eqns
2707 1811 (otherEqnsInts,otherVarsIntsLst,_) := List.map_3(innerEquations, BackendDAEUtil.getEqnAndVarsFromInnerEquation);
2708 1811 otherEqnsLst := BackendEquation.getList(otherEqnsInts, inEqns);
2709 1811 otherEqnsLst := BackendEquation.replaceDerOpInEquationList(otherEqnsLst);
2710 1811 outOtherEqns := BackendEquation.listEquation(otherEqnsLst);
2711
2712 // get other vars
2713 1811 otherVarsInts := List.flatten(otherVarsIntsLst);
2714
4/4
✓ Branch 0 taken 29347 times.
✓ Branch 1 taken 1811 times.
✓ Branch 2 taken 29347 times.
✓ Branch 3 taken 1811 times.
 31158 ovarsLst := list(BackendVariable.transformXToXd(BackendVariable.getVarAt(inVars, e)) for e in otherVarsInts);
2715 1811 outOtherVars := BackendVariable.listVar1(ovarsLst);
2716
2717 // debug
2718
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1811 times.
 1811 if Flags.isSet(Flags.DEBUG_ALGLOOP_JACOBIAN) then
2719 ✗ print("*** got residual equation and created corresponding variables at time: " + realString(clock()) + "\n");
2720 ✗ print("other Equations:\n");
2721 ✗ BackendDump.printEquationList(otherEqnsLst);
2722 ✗ print("other Variables:\n");
2723 ✗ BackendDump.printVarList(ovarsLst);
2724 end if;
2725 else
2726 12 fail();
2727 end try;
2728 end prepareTornStrongComponentData;
2729
2730 protected function checkForSymbolicJacobian
2731 input list<BackendDAE.Equation> inResidualEqns;
2732 input list<BackendDAE.Equation> inOtherEqns;
2733 input String name;
2734 output Boolean out;
2735 protected
2736 Boolean b1, b2;
2737 algorithm
2738
1/2
✓ Branch 1 taken 1257 times.
✗ Branch 2 not taken.
 1257 if not Flags.isSet(Flags.FORCE_NLS_ANALYTIC_JACOBIAN) then
2739 try // this might fail because of algorithms TODO: fix it!
2740 1257 (b1, _) := BackendEquation.traverseExpsOfEquationList_WithStop(inResidualEqns, traverserhasEqnNonDiffParts, ({}, true, false));
2741 1257 (b2, _) := BackendEquation.traverseExpsOfEquationList_WithStop(inOtherEqns, traverserhasEqnNonDiffParts, ({}, true, false));
2742
2/2
✓ Branch 0 taken 689 times.
✓ Branch 1 taken 547 times.
 1236 if not (b1 and b2) then
2743
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 689 times.
 689 if Flags.isSet(Flags.FAILTRACE) then
2744 ✗ Debug.traceln("Skip symbolic jacobian for non-linear system " + name + "\n");
2745 end if;
2746 out := false;
2747 else
2748 out := true;
2749 end if;
2750 else
2751 out := false;
2752 end try;
2753 else
2754 out := true;
2755 end if;
2756 end checkForSymbolicJacobian;
2757
2758 protected function calculateTearingSetJacobian
2759 input BackendDAE.Variables inVars;
2760 input BackendDAE.EquationArray inEqns;
2761 input BackendDAE.TearingSet inTearingSet;
2762 input BackendDAE.Shared inShared;
2763 input Boolean isLinear;
2764 output BackendDAE.Jacobian outJacobian;
2765 output BackendDAE.Shared outShared;
2766 protected
2767 String name, prename;
2768 Boolean debug = false, onlySparsePattern=false;
2769
2770 BackendDAE.Variables diffVars, oVars, resVars;
2771 BackendDAE.EquationArray resEqns, oEqns;
2772 algorithm
2773 try
2774 // check non-linear flag
2775
3/4
✓ Branch 0 taken 854 times.
✓ Branch 1 taken 969 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 854 times.
 1823 if not isLinear and not Flags.isSet(Flags.NLS_ANALYTIC_JACOBIAN) then
2776 onlySparsePattern := true;
2777 end if;
2778 // generate jacobian name
2779
2/2
✓ Branch 0 taken 854 times.
✓ Branch 1 taken 969 times.
 1823 if isLinear then
2780 prename := "LS";
2781 else
2782 prename := "NLS";
2783 end if;
2784 1823 name := prename + "Jac" + intString(System.tmpTickIndex(Global.backendDAE_jacobianSeq));
2785
2786 if debug then
2787 print("*** "+ prename + "-JAC *** start creating Jacobian for a torn system " + name + " of size " + intString(listLength(inTearingSet.tearingvars)) + " time: " + realString(clock()) + "\n");
2788 end if;
2789
2790 1823 (diffVars, resVars, oVars, resEqns, oEqns) := prepareTornStrongComponentData(inVars, inEqns, inTearingSet.tearingvars, inTearingSet.residualequations, inTearingSet.innerEquations, inShared.functionTree, name);
2791
2792 if debug then
2793 print("*** "+ prename + "-JAC *** prepared all data for differentiation at time: " + realString(clock()) + "\n");
2794 end if;
2795
2796 //check if we are able to calc symbolic jacobian
2797
4/4
✓ Branch 0 taken 842 times.
✓ Branch 1 taken 969 times.
✓ Branch 5 taken 573 times.
✓ Branch 6 taken 269 times.
 1811 if not (isLinear or checkForSymbolicJacobian(BackendEquation.equationList(resEqns), BackendEquation.equationList(oEqns), name)) then
2798 onlySparsePattern := true;
2799 end if;
2800
2801 // generate generic jacobian backend dae
2802 1811 (outJacobian, outShared) := getSymbolicJacobian(diffVars, resEqns, resVars, oEqns, oVars, inShared, inVars, name, onlySparsePattern);
2803 else
2804 12 fail();
2805 end try;
2806 end calculateTearingSetJacobian;
2807
2808 protected function calculateJacobianComponent
2809 "Calculates jacobian matrix for strong components of torn systems and non-linear systems."
2810 input BackendDAE.StrongComponent inComp;
2811 input BackendDAE.Variables inVars;
2812 input BackendDAE.EquationArray inEqns;
2813 input BackendDAE.Shared inShared;
2814 output BackendDAE.StrongComponent outComp;
2815 output BackendDAE.Shared outShared;
2816 algorithm
2817 (outComp, outShared) := matchcontinue inComp
2818 local
2819 BackendDAE.StrongComponent comp;
2820 BackendDAE.Shared shared;
2821 list<Integer> iterationvarsInts;
2822 list<Integer> residualequations;
2823
2824
2825 list<BackendDAE.Var> iterationvars, resVarsLst;
2826 BackendDAE.Variables diffVars, ovars, resVars;
2827 list<BackendDAE.Equation> reqns;
2828 BackendDAE.EquationArray eqns, oeqns;
2829
2830 BackendDAE.Jacobian jacobian,jacobianCausal;
2831
2832 String name;
2833 Boolean mixedSystem, linear;
2834
2835 Boolean onlySparsePattern = true;
2836 BackendDAE.TearingSet strictTearingset, casualTearingSet;
2837 Option<BackendDAE.TearingSet> optCasualTearingSet;
2838
2839 // generate symbolic jacobian for a torn system
2840 case BackendDAE.TORNSYSTEM(strictTearingset, optCasualTearingSet, linear, mixedSystem)
2841 algorithm
2842 // generate generic jacobian backend dae
2843 1820 (jacobian, shared) := calculateTearingSetJacobian(inVars, inEqns, strictTearingset, inShared, linear);
2844
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1808 times.
 1808 strictTearingset.jac := jacobian;
2845
2846
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1808 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 1805 times.
 1808 if isSome(optCasualTearingSet) then
2847 3 casualTearingSet := Util.getOption(optCasualTearingSet);
2848 3 (jacobianCausal, shared) := calculateTearingSetJacobian(inVars, inEqns, casualTearingSet, shared, linear);
2849 3 casualTearingSet.jac := jacobianCausal;
2850 optCasualTearingSet := SOME(casualTearingSet);
2851 end if;
2852
4/4
✓ Branch 0 taken 1601 times.
✓ Branch 1 taken 207 times.
✓ Branch 2 taken 839 times.
✓ Branch 3 taken 969 times.
 4248 then (BackendDAE.TORNSYSTEM(strictTearingset, optCasualTearingSet, linear, mixedSystem), shared);
2853
2854 // do not touch constant systems for now
2855 case comp as BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_CONSTANT()) then (comp, inShared);
2856
2857 // Convert linear system to a torn system with symbolica jacobian, when flag is enabled
2858 case BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_LINEAR(), eqns=residualequations, vars=iterationvarsInts, mixedSystem=mixedSystem)
2859 guard(Flags.isSet(Flags.LS_ANALYTIC_JACOBIAN))
2860 algorithm
2861 ✗ strictTearingset := BackendDAE.TEARINGSET(iterationvarsInts, residualequations, {}, BackendDAE.EMPTY_JACOBIAN());
2862 ✗ (jacobian, shared) := calculateTearingSetJacobian(inVars, inEqns, strictTearingset, inShared, true);
2863 ✗ strictTearingset.jac := jacobian;
2864 ✗ then (BackendDAE.TORNSYSTEM(strictTearingset, NONE(), true, mixedSystem), shared);
2865
2866 // Do not touch linear system
2867 case comp as BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_LINEAR()) then (comp, inShared);
2868
2869 case BackendDAE.EQUATIONSYSTEM(eqns=residualequations, vars=iterationvarsInts, mixedSystem=mixedSystem)
2870 algorithm
2871 //generate jacobian name
2872 415 name := "NLSJac" + intString(System.tmpTickIndex(Global.backendDAE_jacobianSeq));
2873
2874 // get iteration vars
2875 415 iterationvars := List.map1r(iterationvarsInts, BackendVariable.getVarAt, inVars);
2876 415 iterationvars := List.map(iterationvars, BackendVariable.transformXToXd);
2877 415 iterationvars := listReverse(iterationvars);
2878 415 diffVars := BackendVariable.listVar1(iterationvars);
2879
2880 // get residual eqns
2881 415 reqns := BackendEquation.getList(residualequations, inEqns);
2882 415 reqns := BackendEquation.replaceDerOpInEquationList(reqns);
2883
2884 //check if we are able to calc symbolic jacobian
2885
3/4
✓ Branch 1 taken 278 times.
✓ Branch 2 taken 137 times.
✓ Branch 4 taken 278 times.
✗ Branch 5 not taken.
 415 if checkForSymbolicJacobian(reqns, {}, name) and Flags.isSet(Flags.NLS_ANALYTIC_JACOBIAN) then
2886 onlySparsePattern := false;
2887 end if;
2888
2889 415 eqns := BackendEquation.listEquation(reqns);
2890 // create residual equations
2891 415 (_, reqns) := BackendEquation.traverseEquationArray(eqns, BackendEquation.traverseEquationToScalarResidualForm, (inShared.functionTree, {}));
2892 413 reqns := listReverse(reqns);
2893 413 (reqns, resVarsLst) := BackendEquation.convertResidualsIntoSolvedEquations(reqns, "$res_" + name + "_", 1);
2894 413 resVars := BackendVariable.listVar1(resVarsLst);
2895 413 eqns := BackendEquation.listEquation(reqns);
2896
2897 // other eqns and vars are empty
2898 413 oeqns := BackendEquation.listEquation({});
2899 413 ovars := BackendVariable.emptyVars();
2900
2901 // generate generic jacobian backend dae
2902 413 (jacobian, shared) := getSymbolicJacobian(diffVars, eqns, resVars, oeqns, ovars, inShared, inVars, name, onlySparsePattern);
2903
1/2
✓ Branch 0 taken 413 times.
✗ Branch 1 not taken.
 826 then (BackendDAE.EQUATIONSYSTEM(residualequations, iterationvarsInts, jacobian, BackendDAE.JAC_GENERIC(), mixedSystem), shared);
2904
2905 case comp then (comp, inShared);
2906 end matchcontinue;
2907
2908 // Check if all nonlinear iteration variables have start values
2909
2/2
✓ Branch 1 taken 102540 times.
✓ Branch 2 taken 58006 times.
 160546 if BackendDAEUtil.isInitializationDAE(inShared) then
2910 try
2911 102540 checkNonLinDependecies(outComp,inEqns);
2912 else
2913 ✗ Error.addInternalError("function calculateJacobianComponent failed to check all non-linear iteration variables for start values.", sourceInfo());
2914 end try;
2915 end if;
2916 end calculateJacobianComponent;
2917
2918 protected function checkNonLinDependecies
2919 "Check if all non-linear iteartion variables of given non-linear equation
2920 system have a start value and throw warning if not. Only start values for
2921 those have an influence on solver iteration."
2922 input BackendDAE.StrongComponent inComp;
2923 input BackendDAE.EquationArray inEqns;
2924 protected
2925 String name, msg;
2926 Boolean existNonLin;
2927 algorithm
2928
2/2
✓ Branch 1 taken 1226 times.
✓ Branch 2 taken 101314 times.
 102540 if Flags.isSet(Flags.INITIALIZATION) then
2929 // Dump full information.
2930 () := match inComp
2931 local
2932 BackendDAE.Jacobian jac;
2933 list<Integer> resIndices, eqnIndices = {};
2934 BackendDAE.InnerEquations innerEquations;
2935 Boolean linear;
2936 // Case non-linear torn equation system
2937 case BackendDAE.TORNSYSTEM(strictTearingSet=BackendDAE.TEARINGSET(jac=jac, residualequations=resIndices, innerEquations=innerEquations), linear=false)
2938 algorithm
2939
3/5
✓ Branch 0 taken 504 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 504 times.
✓ Branch 4 taken 6 times.
 510 for eq in innerEquations loop
2940 eqnIndices := match eq
2941 local
2942 Integer idx;
2943 case BackendDAE.INNEREQUATION(eqn = idx) then idx::eqnIndices;
2944 case BackendDAE.INNEREQUATIONCONSTRAINTS(eqn = idx) then idx::eqnIndices;
2945 else eqnIndices;
2946 end match;
2947 end for;
2948 6 eqnIndices := listAppend(resIndices,eqnIndices);
2949 6 printNonLinIterVarsAndEqs(jac,eqnIndices,inEqns);
2950 then();
2951
2952 // Case non-linear non-torn equation system
2953 case BackendDAE.EQUATIONSYSTEM(eqns=eqnIndices, jac=jac, jacType=BackendDAE.JAC_NONLINEAR())
2954 algorithm
2955 ✗ printNonLinIterVarsAndEqs(jac,eqnIndices,inEqns);
2956 then();
2957
2958 // ToDo: Check if jacType=BackendDAE.JAC_GENERIC is needed
2959 //case BackendDAE.EQUATIONSYSTEM(jac=jac, jacType=BackendDAE.JAC_GENERIC())
2960 else();
2961 end match;
2962 else
2963 // Only error message.
2964 (existNonLin, name) := match inComp
2965 local
2966 BackendDAE.Jacobian jac;
2967 Boolean linear;
2968 // Case non-linear teared equation system
2969 case BackendDAE.TORNSYSTEM(strictTearingSet=BackendDAE.TEARINGSET(jac=jac), linear=false)
2970 415 then existNonLinIterVars(jac);
2971
2972 // Case non-linear non-teared equation system
2973 case BackendDAE.EQUATIONSYSTEM(jac=jac, jacType=BackendDAE.JAC_NONLINEAR())
2974 ✗ then existNonLinIterVars(jac);
2975
2976 // ToDo: Check if jacType=BackendDAE.JAC_GENERIC is needed
2977 //case BackendDAE.EQUATIONSYSTEM(jac=jac, jacType=BackendDAE.JAC_GENERIC())
2978 else (false,"");
2979 end match;
2980
2/2
✓ Branch 0 taken 386 times.
✓ Branch 1 taken 29 times.
 415 if existNonLin then
2981 msg := "For more information set -d=initialization. In OMEdit Tools->Options->Simulation->Show additional information from the initialization process, in OMNotebook call setCommandLineOptions(\"-d=initialization\")";
2982 29 Error.addMessage(Error.INITIALIZATION_ITERATION_VARIABLES, {name, msg});
2983 end if;
2984 end if;
2985 end checkNonLinDependecies;
2986
2987 protected function existNonLinIterVars
2988 "Helper function for checkNonLinDependecies. Returns true if any non-linear
2989 iteration variables without start value are contained in given jacobian."
2990 input BackendDAE.Jacobian jacobian_in;
2991 output Boolean existNonLin;
2992 output String jacName;
2993 algorithm
2994 (existNonLin, jacName) := match jacobian_in
2995 local
2996 list<BackendDAE.Var> diffVars;
2997 list<DAE.ComponentRef> dependentVarsCref;
2998 DAE.ComponentRef varCref;
2999 BackendDAE.Var var;
3000 String name;
3001 Boolean exist=false;
3002 case BackendDAE.GENERIC_JACOBIAN(SOME((_,name,diffVars,_,_,dependentVarsCref))) algorithm
3003 // Search for non-linear variables without start value
3004
2/2
✓ Branch 0 taken 370 times.
✓ Branch 1 taken 122 times.
 492 for varCref in dependentVarsCref loop
3005
2/2
✓ Branch 0 taken 11935 times.
✓ Branch 1 taken 341 times.
 12276 for var in diffVars loop
3006
2/2
✓ Branch 1 taken 274 times.
✓ Branch 2 taken 11661 times.
 11935 if ComponentReferenceBasics.crefEqual(varCref, var.varName) then
3007
2/2
✓ Branch 1 taken 245 times.
✓ Branch 2 taken 29 times.
 274 if (not BackendVariable.varHasStartValue(var)) then
3008 exist:= true;
3009 break;
3010 end if;
3011 end if;
3012 end for;
3013
2/2
✓ Branch 0 taken 341 times.
✓ Branch 1 taken 29 times.
 370 if exist then
3014 break;
3015 end if;
3016 end for;
3017 then (exist, name);
3018
3019 // ToDo
3020 // case BackendDAE.FULL_JACOBIAN() algorithm
3021 else (false, "");
3022 end match;
3023 end existNonLinIterVars;
3024
3025 protected function printNonLinIterVarsAndEqs
3026 "Helper function for checkNonLinDependecies. Prints relevant information regarding
3027 start attributes of non linear iteration variables."
3028 input BackendDAE.Jacobian jacobian;
3029 input list<Integer> eqnIndices;
3030 input BackendDAE.EquationArray inEqns;
3031 algorithm
3032 () := match jacobian
3033 local
3034 list<BackendDAE.Var> diffVars, allDiffedVars, nonLin = {}, nonLinStart = {}, lin = {};
3035 list<DAE.ComponentRef> dependentVarsCref;
3036 DAE.ComponentRef varCref;
3037 BackendDAE.Var var;
3038 String name;
3039 case BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((BackendDAE.DAE({_}, _),name,diffVars,_,allDiffedVars,dependentVarsCref)))
3040 algorithm
3041 // Get non-linear variables without start value
3042
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 1 time.
 8 for varCref in dependentVarsCref loop
3043
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 7 times.
 28 for var in diffVars loop
3044
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 18 times.
 21 if ComponentReferenceBasics.crefEqual(varCref, var.varName) then
3045
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 if (not BackendVariable.varHasStartValue(var)) then
3046 nonLin := var::nonLin;
3047 else
3048 nonLinStart := var::nonLinStart;
3049 end if;
3050 end if;
3051 end for;
3052 end for;
3053
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if not listEmpty(nonLin) then
3054 1 BackendDump.dumpVarList(nonLin, "Nonlinear iteration variables with default zero start attribute in " + name + ".");
3055 end if;
3056
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if not listEmpty(nonLinStart) then
3057 ✗ BackendDump.dumpVarList(nonLinStart, "Nonlinear iteration variables with predefined start attribute in " + name + ".");
3058 end if;
3059
3060 // Get linear variables with start value, but ignore discrete vars
3061 // kabdelhak: i don't get this, how are these the linear ones? these are the inner variables
3062
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1 time.
 15 for var in allDiffedVars loop
3063
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 14 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 14 if (BackendVariable.varHasStartValue(var) and not BackendVariable.isVarDiscrete(var) ) then
3064 lin := var::lin;
3065 end if;
3066 end for;
3067
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if not listEmpty(lin) then
3068 ✗ BackendDump.dumpVarList(lin, "Linear iteration variables with predefined start attributes that are unrelevant in " + name + ".");
3069 end if;
3070
3071
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 1 if not (listEmpty(nonLin) and listEmpty(nonLinStart) and listEmpty(lin)) then
3072
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
 2 print("Info: Only non-linear iteration variables in non-linear eqation systems require start values."
3073 + " All other start values have no influence on convergence and are ignored."
3074 + (if Flags.isSet(Flags.DUMP_LOOPS) then "\n\n"
3075 else " Use \"-d=dumpLoops\" to show all loops. In OMEdit Tools->Options->Simulation->Additional Translation Flags,"
3076 + " in OMNotebook call setCommandLineOptions(\"-d=dumpLoops\")\n\n"));
3077 end if;
3078 then();
3079
3080 else();
3081 end match;
3082 // ToDo
3083 // BackendDAE.FULL_JACOBIAN()
3084 end printNonLinIterVarsAndEqs;
3085
3086 public function getNonLinearVariables
3087 "Returns all nonlinear variables for the jacobian."
3088 input BackendDAE.Jacobian jacobian;
3089 output list<BackendDAE.Var> nonLin = {};
3090 algorithm
3091 nonLin := match jacobian
3092 local
3093 list<BackendDAE.Var> diffVars;
3094 list<DAE.ComponentRef> dependentVarsCref;
3095
3096 case BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((_, _,diffVars, _, _, dependentVarsCref)))
3097 algorithm
3098 // nonlinear variables are those appearing in the jacobian
3099
2/2
✓ Branch 0 taken 621 times.
✓ Branch 1 taken 226 times.
 847 for varCref in dependentVarsCref loop
3100
2/2
✓ Branch 0 taken 5966 times.
✓ Branch 1 taken 234 times.
 6200 for var in diffVars loop
3101
2/2
✓ Branch 1 taken 387 times.
✓ Branch 2 taken 5579 times.
 5966 if ComponentReferenceBasics.crefEqual(varCref, var.varName) then
3102 387 var.initNonlinear := true;
3103 nonLin := var::nonLin;
3104 387 break;
3105 end if;
3106 end for;
3107 end for;
3108 then nonLin;
3109
3110 else {};
3111 end match;
3112 // ToDo
3113 // BackendDAE.FULL_JACOBIAN()
3114 end getNonLinearVariables;
3115
3116 protected function traverserhasEqnNonDiffParts
3117 "function breaks differentiation for
3118 currently not working parts of functions"
3119 input DAE.Exp inExp;
3120 input tuple<list<DAE.Exp>, Boolean, Boolean> inTpl;
3121 output DAE.Exp outExp;
3122 output Boolean cont;
3123 output tuple<list<DAE.Exp>, Boolean, Boolean> outTpl = inTpl;
3124 protected
3125 list<DAE.Exp> expList;
3126 algorithm
3127 11571 (outExp, (expList, cont, _)) := Expression.traverseExpTopDown(inExp, hasEqnNonDiffParts, inTpl);
3128
1/4
✓ Branch 1 taken 11571 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 11571 if Flags.isSet(Flags.DUMP_EXCLUDED_EXP) and not cont then
3129 ✗ print("Traverser for catching functions, that should not be differentiated\n");
3130 ✗ print(stringDelimitList(List.map(expList, ExpressionBasics.printExpStr), "\n"));
3131 ✗ print("\n\n");
3132 end if;
3133 end traverserhasEqnNonDiffParts;
3134
3135 protected function hasEqnNonDiffParts
3136 "function breaks differentiation for
3137 currently not working parts of functions"
3138 input DAE.Exp inExp;
3139 input tuple<list<DAE.Exp>, Boolean, Boolean> inTpl;
3140 output DAE.Exp outExp;
3141 output Boolean cont;
3142 output tuple<list<DAE.Exp>, Boolean, Boolean> outTpl;
3143 algorithm
3144 (outExp, cont, outTpl) := match(inExp, inTpl)
3145 local
3146 list<DAE.Exp> expLst;
3147 Boolean b, insideCall;
3148
3149 ✗ case (DAE.CALL(path=Absyn.IDENT("delay")), (expLst, _, insideCall)) then (inExp, false, (inExp::expLst, false, insideCall));
3150
3151 // For now exclude all not built in calls
3152
1/2
✓ Branch 0 taken 807 times.
✗ Branch 1 not taken.
 1614 case (DAE.CALL(attr=DAE.CALL_ATTR(builtin=false)), (expLst, _, insideCall)) then (inExp, false, (inExp::expLst, false, insideCall));
3153
3154 /*
3155 case (_, (expLst, _, true)) guard(Expression.isRecord(inExp)) then (inExp, false, (inExp::expLst, false, true));
3156 case (_, (expLst, _, true)) guard(Expression.isMatrix(inExp)) then (inExp, false, (inExp::expLst, false, true));
3157 case (DAE.CALL(attr=DAE.CALL_ATTR(ty = ty, builtin=false)), (expLst, b, insideCall))
3158 algorithm
3159 true = isRecordInvoled(ty);
3160 then (inExp, false, (inExp::expLst, false, insideCall));
3161 case (DAE.CALL(expLst=expLst1,attr=DAE.CALL_ATTR(builtin=false)), (expLst, b, insideCall))
3162 algorithm
3163 (_, (_, false, _)) = Expression.traverseExpListTopDown(expLst1, hasEqnNonDiffParts, (expLst, b, true));
3164 then (inExp, false, (inExp::expLst, false, insideCall));
3165 */
3166
3167 case (outExp, (_, b, _)) then (outExp, b, inTpl);
3168 end match;
3169 end hasEqnNonDiffParts;
3170
3171 protected function isRecordInvoled
3172 input DAE.Type inType;
3173 output Boolean out;
3174 algorithm
3175 out := match inType
3176 local
3177 DAE.Type ty;
3178 list<DAE.Type> types;
3179 case DAE.T_COMPLEX() then true;
3180 ✗ case DAE.T_ARRAY(ty=ty) then isRecordInvoled(ty);
3181 ✗ case DAE.T_FUNCTION(funcResultType=ty) then isRecordInvoled(ty);
3182 case DAE.T_TUPLE(types)
3183 ✗ then List.any(types, isRecordInvoled);
3184 else false;
3185 end match;
3186 end isRecordInvoled;
3187
3188 public function getSymbolicJacobian "author: wbraun
3189 This function creates a symbolic Jacobian column for non-linear systems and
3190 tearing systems."
3191 input BackendDAE.Variables inDiffVars;
3192 input BackendDAE.EquationArray inResEquations;
3193 input BackendDAE.Variables inResVars;
3194 input BackendDAE.EquationArray inotherEquations;
3195 input BackendDAE.Variables inotherVars;
3196 input BackendDAE.Shared inShared;
3197 input BackendDAE.Variables inAllVars;
3198 input String inName;
3199 input Boolean inOnlySparsePattern;
3200 output BackendDAE.Jacobian outJacobian;
3201 output BackendDAE.Shared outShared;
3202 protected
3203 BackendDAE.BackendDAE backendDAE;
3204 BackendDAE.EquationArray eqns;
3205 BackendDAE.ExtraInfo einfo;
3206 BackendDAE.Shared shared;
3207 BackendDAE.SparseColoring sparseColoring;
3208 BackendDAE.SparsePattern sparsePattern;
3209 BackendDAE.NonlinearPattern nonlinearPattern;
3210 BackendDAE.Variables dependentVars, globalKnownVars;
3211 AvlTreePathFunction.Tree funcs;
3212 FCore.Cache cache;
3213 FCore.Graph graph;
3214 list<BackendDAE.Var> knvarLst1, knvarLst2, independentVarsLst, dependentVarsLst, otherVarsLst;
3215 list<DAE.ComponentRef> independentComRefs, otherVarsLstComRefs;
3216 Option<BackendDAE.SymbolicJacobian> symJacBDAE;
3217 algorithm
3218 try
3219 2293 globalKnownVars := BackendDAEUtil.getGlobalKnownVarsFromShared(inShared);
3220 2293 funcs := BackendDAEUtil.getFunctions(inShared);
3221 2293 einfo := BackendDAEUtil.getExtraInfo(inShared);
3222
3223
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2293 times.
 2293 if Flags.isSet(Flags.JAC_DUMP2) then
3224 ✗ print("---+++ create analytical jacobian +++---");
3225 ✗ print("\n---+++ independent variables +++---\n");
3226 ✗ BackendDump.printVariables(inDiffVars);
3227 ✗ print("\n---+++ equation system +++---\n");
3228 ✗ BackendDump.printEquationArray(inResEquations);
3229 end if;
3230
3231 2293 independentVarsLst := BackendVariable.varList(inDiffVars);
3232 2293 independentComRefs := List.map(independentVarsLst, BackendVariable.varCref);
3233
3234 2293 otherVarsLst := BackendVariable.varList(inotherVars);
3235 2293 otherVarsLstComRefs := List.map(otherVarsLst, BackendVariable.varCref);
3236
3237
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2293 times.
 2293 if Flags.isSet(Flags.JAC_DUMP2) then
3238 ✗ print("\n---+++ known variables +++---\n");
3239 ✗ BackendDump.printVariables(globalKnownVars);
3240 end if;
3241
3242 // dependentVarsLst = listReverse(dependentVarsLst);
3243 2293 dependentVars := BackendVariable.mergeVariables(inResVars, inotherVars);
3244 2293 eqns := BackendEquation.merge(inResEquations, inotherEquations);
3245
3246
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2293 times.
 2293 if Flags.isSet(Flags.JAC_DUMP2) then
3247 ✗ print("\n---+++ created backend system +++---\n");
3248 ✗ print("\n---+++ vars +++---\n");
3249 ✗ BackendDump.printVariables(dependentVars);
3250 ✗ print("\n---+++ equations +++---\n");
3251 ✗ BackendDump.printEquationArray(eqns);
3252 end if;
3253
3254 // create known variables
3255 2293 knvarLst1 := BackendEquation.equationsVars(eqns, globalKnownVars);
3256 //knvarLst2 := BackendEquation.equationsVars(eqns, inAllVars);
3257 knvarLst2 := {};
3258 // Create a list of known variables true *only* for this shared system
3259 2293 globalKnownVars := BackendVariable.listVar2(knvarLst1,knvarLst2);
3260 // Remove inputs for the jacobian
3261 2293 globalKnownVars := BackendVariable.removeCrefs(independentComRefs, globalKnownVars);
3262 2293 globalKnownVars := BackendVariable.removeCrefs(otherVarsLstComRefs, globalKnownVars);
3263
3264
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2293 times.
 2293 if Flags.isSet(Flags.JAC_DUMP2) then
3265 ✗ print("\n---+++ known variables +++---\n");
3266 ✗ BackendDump.printVariables(globalKnownVars);
3267 end if;
3268
3269 // prepare vars and equations for BackendDAE
3270 2293 cache := FCore.emptyCache();
3271 2293 graph := FGraph.empty();
3272 2293 shared := BackendDAEUtil.createEmptyShared(BackendDAE.ALGEQSYSTEM(), einfo, cache, graph);
3273 2293 shared := BackendDAEUtil.setSharedGlobalKnownVars(shared, globalKnownVars);
3274 2293 shared := BackendDAEUtil.setSharedFunctionTree(shared, funcs);
3275 4586 backendDAE := BackendDAE.DAE({BackendDAEUtil.createEqSystem(dependentVars, eqns)}, shared);
3276
3277
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2293 times.
 2293 if Flags.isSet(Flags.JAC_DUMP2) then
3278 ✗ BackendDump.bltdump("System",backendDAE);
3279 end if;
3280
3281 2293 backendDAE := BackendDAEUtil.transformBackendDAE(backendDAE, SOME((BackendDAE.NO_INDEX_REDUCTION(), BackendDAE.EXACT())), NONE(), NONE());
3282
3283
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2293 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2293 times.
 2293 BackendDAE.DAE({BackendDAE.EQSYSTEM(orderedVars = dependentVars)}, BackendDAE.SHARED(globalKnownVars = globalKnownVars)) := backendDAE;
3284
3285 // prepare creation of symbolic jacobian
3286 // create dependent variables
3287 2293 dependentVarsLst := BackendVariable.varList(dependentVars);
3288
3289 2293 (symJacBDAE, funcs, sparsePattern, sparseColoring, nonlinearPattern) := generateGenericJacobian(backendDAE,
3290 independentVarsLst,
3291 BackendVariable.emptyVars(),
3292 BackendVariable.emptyVars(),
3293 globalKnownVars,
3294 inResVars,
3295 dependentVarsLst,
3296 inName,
3297 inOnlySparsePattern);
3298
3299 2292 outJacobian := BackendDAE.GENERIC_JACOBIAN(symJacBDAE, sparsePattern, sparseColoring, nonlinearPattern);
3300 2292 outShared := BackendDAEUtil.setSharedFunctionTree(inShared, funcs);
3301 else
3302
3303
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if Flags.isSet(Flags.JAC_DUMP) then
3304 ✗ Error.addInternalError("function getSymbolicJacobian failed", sourceInfo());
3305 end if;
3306 outJacobian := BackendDAE.EMPTY_JACOBIAN();
3307 outShared := inShared;
3308 end try;
3309 end getSymbolicJacobian;
3310
3311 public function hasGenericSymbolicJacobian
3312 input BackendDAE.Jacobian inJacobian;
3313 output Boolean out;
3314 algorithm
3315 out := match inJacobian
3316 case BackendDAE.GENERIC_JACOBIAN(jacobian=SOME(_)) then true;
3317 else false;
3318 end match;
3319 end hasGenericSymbolicJacobian;
3320
3321 protected function calculateEqSystemStateSetsJacobians
3322 input BackendDAE.EqSystem inSyst;
3323 input BackendDAE.Shared inShared;
3324 output BackendDAE.EqSystem outSyst;
3325 output BackendDAE.Shared outShared;
3326 algorithm
3327 (outSyst,outShared) := match (inSyst, inShared)
3328 local
3329 BackendDAE.EqSystem syst;
3330 BackendDAE.Shared shared;
3331 BackendDAE.StrongComponents comps;
3332 BackendDAE.Variables vars;
3333 BackendDAE.EquationArray eqns;
3334 BackendDAE.StateSets stateSets;
3335
3336 case (syst as BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs = eqns, stateSets=stateSets), shared)
3337 algorithm
3338 5490 comps := BackendDAEUtil.getStrongComponents(syst);
3339 5490 (stateSets, shared) := calculateStateSetsJacobian(stateSets, vars, eqns, comps, shared);
3340
1/2
✓ Branch 0 taken 5490 times.
✗ Branch 1 not taken.
 5490 syst.stateSets := stateSets;
3341
1/2
✓ Branch 0 taken 5490 times.
✗ Branch 1 not taken.
 5490 then (syst, shared);
3342 end match;
3343 end calculateEqSystemStateSetsJacobians;
3344
3345 protected function calculateStateSetsJacobian
3346 input BackendDAE.StateSets inStateSets;
3347 input BackendDAE.Variables inVars;
3348 input BackendDAE.EquationArray inEqns;
3349 input BackendDAE.StrongComponents inComps;
3350 input BackendDAE.Shared inShared;
3351 output BackendDAE.StateSets outStateSets;
3352 output BackendDAE.Shared outShared = inShared;
3353 algorithm
3354
4/4
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 5490 times.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 5490 times.
 5532 outStateSets := list(match s
3355 local
3356 BackendDAE.StateSet stateSet;
3357 case stateSet
3358 algorithm
3359 42 (stateSet, outShared) := calculateStateSetJacobian(stateSet, inVars, inEqns, inComps, outShared);
3360 then stateSet;
3361 end match for s in inStateSets);
3362 end calculateStateSetsJacobian;
3363
3364 protected function calculateStateSetJacobian
3365 input BackendDAE.StateSet inStateSet;
3366 input BackendDAE.Variables inVars;
3367 input BackendDAE.EquationArray inEqns;
3368 input BackendDAE.StrongComponents inComps;
3369 input BackendDAE.Shared inShared;
3370 output BackendDAE.StateSet outStateSet;
3371 output BackendDAE.Shared outShared;
3372 algorithm
3373 (outStateSet, outShared) := match inStateSet
3374 local
3375 BackendDAE.Shared shared;
3376
3377 Integer index, rang;
3378 list<DAE.ComponentRef> state;
3379 DAE.ComponentRef crA, crJ;
3380 list<BackendDAE.Var> varA, varJ, statescandidates, ovars;
3381
3382 list<DAE.ComponentRef> crstates;
3383 array<Boolean> marked;
3384 HashSet.HashSet hs;
3385
3386 list<BackendDAE.Var> statevars, compvars;
3387 BackendDAE.Variables diffVars, allvars, oVars, resVars;
3388 list<BackendDAE.Equation> eqns, compeqns, ceqns, oeqns;
3389 BackendDAE.EquationArray cEqns, oEqns;
3390
3391 BackendDAE.Jacobian jacobian;
3392
3393 String name;
3394
3395 case BackendDAE.STATESET(index=index, rang=rang, state=state, crA=crA, varA=varA, statescandidates=statescandidates,
3396 ovars=ovars, eqns=eqns, oeqns=oeqns, crJ=crJ, varJ=varJ)
3397 algorithm
3398 // get state names
3399 42 crstates := List.map(statescandidates, BackendVariable.varCref);
3400 42 marked := arrayCreate(BackendVariable.varsSize(inVars), false);
3401 // get Equations for Jac from the strong component
3402 42 marked := List.fold1(crstates, markSetStates, inVars, marked);
3403 42 (compeqns, compvars) := getStateSetCompVarEqns(inComps, marked, inEqns, inVars);
3404 // remove the state set equation
3405 42 compeqns := List.select(compeqns, removeStateSetEqn);
3406 // remove the state candidates to geht the other vars
3407 42 hs := List.fold(crstates, BaseHashSet.add, HashSet.emptyHashSet());
3408 42 compvars := List.select1(compvars, removeStateSetStates, hs);
3409 // match the equations to get the residual equations
3410 42 (ceqns, oeqns) := IndexReduction.splitEqnsinConstraintAndOther(compvars, compeqns, inShared);
3411 // change state vars to ders
3412 42 compvars := List.map(compvars, BackendVariable.transformXToXd);
3413 // replace der in equations
3414 42 ceqns := BackendEquation.replaceDerOpInEquationList(ceqns);
3415 42 oeqns := BackendEquation.replaceDerOpInEquationList(oeqns);
3416 // convert ceqns to res[..] = lhs-rhs
3417 42 ceqns := createResidualSetEquations(ceqns, crJ, 1, intGt(listLength(ceqns), 1));
3418
3419 //add states to allVars
3420 42 allvars := BackendVariable.copyVariables(inVars);
3421 42 statevars := BackendVariable.getAllStateVarFromVariables(allvars);
3422 42 statevars := List.map(statevars, BackendVariable.transformXToXd);
3423 42 allvars := BackendVariable.addVars(statevars, allvars);
3424
3425 // create arrays
3426 42 resVars := BackendVariable.listVar1(varJ);
3427 42 diffVars := BackendVariable.listVar1(statescandidates);
3428 42 oVars := BackendVariable.listVar1(compvars);
3429 42 cEqns := BackendEquation.listEquation(ceqns);
3430 42 oEqns := BackendEquation.listEquation(oeqns);
3431
3432 //generate Jacobian name
3433 42 name := "StateSetJac" + intString(System.tmpTickIndex(Global.backendDAE_jacobianSeq));
3434 // generate generic Jacobian back end dae
3435 42 (jacobian, shared) := getSymbolicJacobian(diffVars, cEqns, resVars, oEqns, oVars, inShared, allvars, name, false);
3436
3437
1/2
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
 42 then (BackendDAE.STATESET(index, rang, state, crA, varA, statescandidates, ovars, eqns, oeqns, crJ, varJ, jacobian), shared);
3438 end match;
3439 end calculateStateSetJacobian;
3440
3441 protected function markSetStates
3442 input DAE.ComponentRef inCr;
3443 input BackendDAE.Variables iVars;
3444 input array<Boolean> iMark;
3445 output array<Boolean> oMark;
3446 protected
3447 Integer index;
3448 algorithm
3449 104 (_, index) := BackendVariable.getVarSingle(inCr, iVars);
3450 104 oMark := arrayUpdate(iMark, index, true);
3451 end markSetStates;
3452
3453 protected function removeStateSetStates
3454 input BackendDAE.Var inVar;
3455 input HashSet.HashSet hs;
3456 output Boolean b;
3457 algorithm
3458 264 b := not BaseHashSet.has(BackendVariable.varCref(inVar), hs);
3459 end removeStateSetStates;
3460
3461 protected function removeStateSetEqn
3462 input BackendDAE.Equation inEqn;
3463 output Boolean b;
3464 algorithm
3465 b := match inEqn
3466 case BackendDAE.ARRAY_EQUATION(source=DAE.SOURCE(info=SOURCEINFO(fileName="stateselection"))) then false;
3467 case BackendDAE.EQUATION(source=DAE.SOURCE(info=SOURCEINFO(fileName="stateselection"))) then false;
3468 else true;
3469 end match;
3470 end removeStateSetEqn;
3471
3472 protected function foundMarked
3473 input list<Integer> ilst;
3474 input array<Boolean> marked;
3475 output Boolean found;
3476 algorithm
3477 found := match ilst
3478 local
3479 Boolean b;
3480 Integer i;
3481 list<Integer> rest;
3482 case {} then false;
3483 case i::rest
3484 algorithm
3485 3220 b := marked[i];
3486
2/2
✓ Branch 0 taken 3178 times.
✓ Branch 1 taken 42 times.
 3220 b := if not b then foundMarked(rest, marked) else b;
3487 then
3488 b;
3489 end match;
3490 end foundMarked;
3491
3492 protected function getStateSetCompVarEqns "author: Frenkel TUD 2013-01
3493 Retrieves the equation and the variable for a state set"
3494 input BackendDAE.StrongComponents inComp;
3495 input array<Boolean> marked;
3496 input BackendDAE.EquationArray inEquationArray;
3497 input BackendDAE.Variables inVariables;
3498 output list<BackendDAE.Equation> outEquations = {};
3499 output list<BackendDAE.Var> outVars = {};
3500 protected
3501 list<Integer> elst, vlst;
3502 list<BackendDAE.Equation> eqnlst;
3503 list<BackendDAE.Var> varlst;
3504 algorithm
3505
2/2
✓ Branch 0 taken 2531 times.
✓ Branch 1 taken 42 times.
 2573 for comp in inComp loop
3506 2531 (elst, vlst) := BackendDAETransform.getEquationAndSolvedVarIndxes(comp);
3507
2/2
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 2489 times.
 2531 if foundMarked(vlst, marked) then
3508 42 eqnlst := BackendEquation.getList(elst, inEquationArray);
3509 42 varlst := List.map1r(vlst, BackendVariable.getVarAt, inVariables);
3510 42 outEquations := listAppend(eqnlst, outEquations);
3511 42 outVars := listAppend(varlst, outVars);
3512 end if;
3513 end for;
3514 end getStateSetCompVarEqns;
3515
3516 protected function createResidualSetEquations
3517 input list<BackendDAE.Equation> iEqs;
3518 input DAE.ComponentRef crJ;
3519 input Integer index;
3520 input Boolean applySubs;
3521 output list<BackendDAE.Equation> oEqs;
3522 protected
3523 Integer idx = index;
3524 algorithm
3525
4/4
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 42 times.
 84 oEqs := list(match eq
3526 local
3527 DAE.ComponentRef crj;
3528 DAE.Exp res, e1, e2, expJ;
3529 BackendDAE.Equation eqn;
3530 DAE.ElementSource source;
3531 BackendDAE.EquationAttributes eqAttr;
3532 case BackendDAE.EQUATION(exp=e1, scalar=e2, source=source, attr=eqAttr)
3533 algorithm
3534
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
 42 crj := if applySubs then ComponentReference.subscriptCrefWithInt(crJ, idx) else crJ;
3535 42 expJ := Expression.crefExp(crj);
3536 42 res := Expression.expSub(e1, e2);
3537 42 eqn := BackendDAE.EQUATION(expJ, res, source, eqAttr);
3538 42 idx := idx + 1;
3539 then eqn;
3540
3541 case BackendDAE.RESIDUAL_EQUATION(exp=e1, source=source, attr=eqAttr)
3542 algorithm
3543 ✗ expJ := Expression.crefExp(ComponentReference.subscriptCrefWithInt(crJ, idx));
3544 ✗ eqn := BackendDAE.EQUATION(expJ, e1, source, eqAttr);
3545 ✗ idx := idx + 1;
3546 then eqn;
3547
3548 case eqn
3549 algorithm
3550 ✗ Error.addInternalError("function createResidualSetEquations failed for equation: " + BackendDump.equationString(eqn), sourceInfo());
3551 ✗ then
3552 fail();
3553 end match for eq in iEqs);
3554 end createResidualSetEquations;
3555
3556 public function calculateJacobian "This function takes an array of equations and the variables of the equation
3557 and calculates the Jacobian of the equations."
3558 input BackendDAE.Variables inVariables;
3559 input BackendDAE.EquationArray inEquationArray;
3560 input BackendDAE.AdjacencyMatrix inAdjacencyMatrix;
3561 input Boolean differentiateIfExp "If true, allow differentiation of if-expressions";
3562 input BackendDAE.Shared iShared;
3563 output Option<list<tuple<Integer, Integer, BackendDAE.Equation>>> outTplIntegerIntegerEquationLstOption;
3564 output BackendDAE.Shared oShared;
3565 algorithm
3566 (outTplIntegerIntegerEquationLstOption, oShared):=
3567 matchcontinue (inVariables, inEquationArray, inAdjacencyMatrix)
3568 local
3569 list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
3570 BackendDAE.Variables vars;
3571 BackendDAE.EquationArray eqns;
3572 BackendDAE.AdjacencyMatrix m;
3573 BackendDAE.Shared shared;
3574 case (vars, eqns, m)
3575 algorithm
3576 5388 (jac, shared) := calculateJacobianRows(eqns,vars,m,1,1,differentiateIfExp,iShared,BackendDAEUtil.varsInEqn);
3577 3431 then
3578 (SOME(jac),shared);
3579 else (NONE(), iShared); /* no analytic jacobian available */
3580 end matchcontinue;
3581 end calculateJacobian;
3582
3583 protected function calculateJacobianRows "author: PA
3584 This function takes a list of Equations and a set of variables and
3585 calculates the Jacobian expression for each variable over each equations,
3586 returned in a sparse matrix representation.
3587 For example, the equation on index e1: 3ax+5yz+ zz given the
3588 variables {x,y,z} on index x1,y1,z1 gives
3589 {(e1,x1,3a), (e1,y1,5z), (e1,z1,5y+2z)}"
3590 replaceable type Type_a subtypeof Any;
3591 input BackendDAE.EquationArray inEquationArray;
3592 input BackendDAE.Variables vars;
3593 input Type_a m;
3594 input Integer eqn_indx;
3595 input Integer scalar_eqn_indx;
3596 input Boolean differentiateIfExp "If true, allow differentiation of if-expressions";
3597 input BackendDAE.Shared iShared;
3598 input varsInEqnFunc varsInEqn;
3599 output list<tuple<Integer, Integer, BackendDAE.Equation>> outLst = {};
3600 output BackendDAE.Shared oShared = iShared;
3601 partial function varsInEqnFunc
3602 input Type_a m;
3603 input Integer indx;
3604 output list<Integer> outIntegerLst;
3605 end varsInEqnFunc;
3606 protected
3607 Integer size, i, j, n, k;
3608 BackendDAE.Equation eqn;
3609 algorithm
3610 i := eqn_indx;
3611 j := scalar_eqn_indx;
3612 5388 size := 0;
3613 5388 n := ExpandableArray.getLastUsedIndex(inEquationArray);
3614 // print("CalcJac(Eqs:" + intString(n) + ")\n");
3615
1/2
✓ Branch 0 taken 5388 times.
✗ Branch 1 not taken.
 56451 for k in 1:n loop
3616
1/2
✓ Branch 1 taken 53020 times.
✗ Branch 2 not taken.
 53020 if ExpandableArray.occupied(k, inEquationArray) then
3617 53020 eqn := ExpandableArray.get(k, inEquationArray);
3618 53020 (outLst, size, oShared) := calculateJacobianRow(eqn, vars, m, i, j, differentiateIfExp, oShared, varsInEqn, outLst);
3619 51063 i := i+1;
3620 51063 j := j+size;
3621 end if;
3622 end for;
3623 3431 outLst := MetaModelica.Dangerous.listReverseInPlace(outLst);
3624 // print("END_CalcJac(Size:" + intString(listLength(outLst)) + ")\n");
3625 end calculateJacobianRows;
3626
3627 protected function calculateJacobianRow "author: PA
3628 Calculates the Jacobian for one equation. See calculateJacobianRows.
3629 inputs: (Equation,
3630 BackendDAE.Variables,
3631 AdjacencyMatrix,
3632 AdjacencyMatrixT,
3633 int /* eqn index */)
3634 outputs: ((int int Equation) list option)"
3635 replaceable type Type_a subtypeof Any;
3636 input BackendDAE.Equation inEquation;
3637 input BackendDAE.Variables vars;
3638 input Type_a m;
3639 input Integer eqn_indx;
3640 input Integer scalar_eqn_indx;
3641 input Boolean differentiateIfExp "If true, allow differentiation of if-expressions";
3642 input BackendDAE.Shared iShared;
3643 input varsInEqnFunc fvarsInEqn;
3644 input list<tuple<Integer, Integer, BackendDAE.Equation>> iAcc;
3645 output list<tuple<Integer, Integer, BackendDAE.Equation>> outLst;
3646 output Integer size;
3647 output BackendDAE.Shared oShared;
3648 partial function varsInEqnFunc
3649 input Type_a m;
3650 input Integer indx;
3651 output list<Integer> outIntegerLst;
3652 end varsInEqnFunc;
3653 algorithm
3654 (outLst, size, oShared):= match inEquation
3655 local
3656 list<Integer> var_indxs,var_indxs_1,ds;
3657 list<tuple<Integer, Integer, BackendDAE.Equation>> eqns;
3658 DAE.Exp e,e1,e2;
3659 list<DAE.Exp> expl;
3660 list<list<DAE.Subscript>> subslst;
3661 DAE.ElementSource source;
3662 DAE.ComponentRef cr;
3663 String str;
3664 BackendDAE.Shared shared;
3665
3666 // residual equations
3667 case BackendDAE.EQUATION(exp = e1,scalar=e2,source=source)
3668 algorithm
3669
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52610 times.
 52610 var_indxs := fvarsInEqn(m, eqn_indx);
3670 // Remove duplicates and get in correct order: ascending index
3671 52610 var_indxs_1 := List.sort(var_indxs,intGt);
3672 52610 var_indxs_1 := List.sortedUnique(var_indxs_1, intEq);
3673 52610 (eqns, shared) := calculateJacobianRow2(Expression.expSub(e1,e2), vars, scalar_eqn_indx, var_indxs_1,differentiateIfExp,iShared,source,iAcc);
3674 51049 then
3675 (eqns, 1, shared);
3676
3677 // residual equations
3678 case BackendDAE.RESIDUAL_EQUATION(exp=e,source=source)
3679 algorithm
3680 ✗ var_indxs := fvarsInEqn(m, eqn_indx);
3681 // Remove duplicates and get in correct order: ascending index
3682 ✗ var_indxs_1 := List.sort(var_indxs,intGt);
3683 ✗ var_indxs_1 := List.sortedUnique(var_indxs_1, intEq);
3684 ✗ (eqns, shared) := calculateJacobianRow2(e, vars, scalar_eqn_indx, var_indxs_1,differentiateIfExp,iShared,source,iAcc);
3685 ✗ then
3686 (eqns, 1, shared);
3687
3688 // solved equations
3689 case BackendDAE.SOLVED_EQUATION(componentRef=cr,exp=e2,source=source)
3690 algorithm
3691 ✗ e1 := Expression.crefExp(cr);
3692
3693 ✗ var_indxs := fvarsInEqn(m, eqn_indx);
3694 // Remove duplicates and get in correct order: ascending index
3695 ✗ var_indxs_1 := List.sort(var_indxs,intGt);
3696 ✗ var_indxs_1 := List.sortedUnique(var_indxs_1, intEq);
3697 ✗ (eqns, shared) := calculateJacobianRow2(Expression.expSub(e1,e2), vars, scalar_eqn_indx, var_indxs_1,differentiateIfExp,iShared,source,iAcc);
3698 ✗ then
3699 (eqns, 1, shared);
3700
3701 // array equations
3702 case BackendDAE.ARRAY_EQUATION(dimSize=ds,left=e1,right=e2,source=source)
3703 algorithm
3704 23 e := Expression.expSub(e1,e2);
3705 23 (e,_) := Expression.extendArrExp(e,false);
3706 23 subslst := Expression.dimensionSizesSubscripts(ds);
3707 23 subslst := Expression.rangesToSubscripts(subslst);
3708 23 expl := List.map1r(subslst,Expression.applyExpSubscripts,e);
3709
3710
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
 23 var_indxs := fvarsInEqn(m, eqn_indx);
3711 // Remove duplicates and get in correct order: ascending index
3712 23 var_indxs_1 := List.sort(var_indxs,intGt);
3713 23 var_indxs_1 := List.sortedUnique(var_indxs_1, intEq);
3714 23 (eqns, shared) := calculateJacobianRowLst(expl, vars, scalar_eqn_indx, var_indxs_1,differentiateIfExp,iShared,source,iAcc);
3715 14 size := List.fold(ds,intMul,1);
3716 14 then
3717 (eqns, size, shared);
3718
3719 else
3720 algorithm
3721
1/2
✓ Branch 1 taken 387 times.
✗ Branch 2 not taken.
 387 true := Flags.isSet(Flags.FAILTRACE);
3722 ✗ str := BackendDump.dumpEqnsStr({inEquation});
3723 ✗ Debug.traceln("- BackendDAE.calculateJacobianRow failed on " + str);
3724 ✗ then
3725 fail();
3726 end match;
3727 end calculateJacobianRow;
3728
3729 protected function calculateJacobianRowLst "author: Frenkel TUD 2012-06
3730 calls calculateJacobianRow2 for a list of DAE.Exp"
3731 input list<DAE.Exp> inExps;
3732 input BackendDAE.Variables vars;
3733 input Integer eqn_indx;
3734 input list<Integer> inIntegerLst;
3735 input Boolean differentiateIfExp "If true, allow differentiation of if-expressions";
3736 input BackendDAE.Shared iShared;
3737 input DAE.ElementSource source;
3738 input list<tuple<Integer, Integer, BackendDAE.Equation>> iAcc;
3739 output list<tuple<Integer, Integer, BackendDAE.Equation>> outLst = iAcc;
3740 output BackendDAE.Shared oShared = iShared;
3741 protected
3742 Integer eqn_indx_arr = eqn_indx;
3743 algorithm
3744
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 14 times.
 65 for e in inExps loop
3745 51 (outLst, oShared) := calculateJacobianRow2(e,vars,eqn_indx_arr,inIntegerLst,differentiateIfExp,oShared,source,outLst);
3746 42 eqn_indx_arr := eqn_indx_arr + 1;
3747 end for;
3748 end calculateJacobianRowLst;
3749
3750 protected function calculateJacobianRow2 "author: PA
3751 Differentiates expression for each variable cref.
3752 inputs: (DAE.Exp,
3753 BackendDAE.Variables,
3754 int, /* equation index */
3755 int list) /* var indexes */
3756 outputs: ((int int Equation) list option)"
3757 input DAE.Exp inExp;
3758 input BackendDAE.Variables vars;
3759 input Integer eqn_indx;
3760 input list<Integer> inIntegerLst;
3761 input Boolean differentiateIfExp "If true, allow differentiation of if-expressions";
3762 input BackendDAE.Shared iShared;
3763 input DAE.ElementSource source;
3764 input list<tuple<Integer, Integer, BackendDAE.Equation>> iAcc;
3765 output list<tuple<Integer, Integer, BackendDAE.Equation>> outLst = iAcc;
3766 output BackendDAE.Shared oShared = iShared;
3767 protected
3768 DAE.Exp e, e_1, dcrexp;
3769 BackendDAE.Var v;
3770 DAE.ComponentRef cr, dcr;
3771 Integer vindx;
3772 String str;
3773 algorithm
3774 try
3775
2/2
✓ Branch 0 taken 151719 times.
✓ Branch 1 taken 51091 times.
 202810 for vindx in inIntegerLst loop
3776 151719 v := BackendVariable.getVarAt(vars, vindx);
3777 151719 cr := BackendVariable.varCref(v);
3778
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 151719 times.
 151719 if BackendVariable.isStateVar(v) then
3779 ✗ dcr := ComponentReference.crefPrefixDer(cr);
3780 ✗ dcrexp := Expression.crefExp(cr);
3781 ✗ dcrexp := DAE.CALL(Absyn.IDENT("der"), {dcrexp}, DAE.callAttrBuiltinReal);
3782 ✗ (e, _) := Expression.replaceExp(inExp, dcrexp, Expression.crefExp(dcr));
3783 end if;
3784 151719 (e_1, oShared) := Differentiate.differentiateExpCrefFullJacobian(inExp, cr, vars, oShared);
3785 // e_1 already simplified in Differentiate.differentiateExpCrefFullJacobian!
3786
2/2
✓ Branch 1 taken 149670 times.
✓ Branch 2 taken 625 times.
 150295 if not Expression.isZero(e_1) then
3787 149670 outLst := (eqn_indx,vindx,BackendDAE.RESIDUAL_EQUATION(e_1,source,BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN))::outLst;
3788 end if;
3789 end for;
3790 else
3791
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1424 times.
 1424 if Flags.isSet(Flags.FAILTRACE) then
3792 ✗ str := ExpressionBasics.printExpStr(inExp);
3793 ✗ Debug.traceln("- BackendDAE.calculateJacobianRow2 failed on " + str);
3794 end if;
3795 1424 fail();
3796 end try;
3797 end calculateJacobianRow2;
3798
3799 protected function addBackendDAESharedJacobian
3800 input Option<BackendDAE.SymbolicJacobian> inSymJac;
3801 input BackendDAE.SparsePattern inSparsePattern;
3802 input BackendDAE.SparseColoring inSparseColoring;
3803 input BackendDAE.NonlinearPattern inNonlinearPattern;
3804 input BackendDAE.Shared inShared;
3805 output BackendDAE.Shared outShared;
3806 protected
3807 BackendDAE.SymbolicJacobians symjacs;
3808 algorithm
3809 6 symjacs := { (inSymJac, inSparsePattern, inSparseColoring, inNonlinearPattern),
3810 (NONE(), ({}, {}, ({}, {}), -1), {}, ({}, {}, ({}, {}), -1)),
3811 (NONE(), ({}, {}, ({}, {}), -1), {}, ({}, {}, ({}, {}), -1)),
3812 (NONE(), ({}, {}, ({}, {}), -1), {}, ({}, {}, ({}, {}), -1))};
3813 6 outShared := BackendDAEUtil.setSharedSymJacs(inShared, symjacs);
3814 end addBackendDAESharedJacobian;
3815
3816 protected function addBackendDAESharedJacobianSparsePattern
3817 input BackendDAE.SparsePattern inSparsePattern;
3818 input BackendDAE.SparseColoring inSparseColoring;
3819 input Integer inIndex;
3820 input BackendDAE.Shared inShared;
3821 output BackendDAE.Shared outShared;
3822 protected
3823 BackendDAE.SymbolicJacobians symjacs;
3824 Option<BackendDAE.SymbolicJacobian> symJac;
3825 BackendDAE.NonlinearPattern nonlinearPattern = BackendDAE.emptyNonlinearPattern;
3826 algorithm
3827 1064 BackendDAE.SHARED(symjacs=symjacs) := inShared;
3828 1064 (symJac, _, _, _) := listGet(symjacs, inIndex);
3829 1064 symjacs := List.set(symjacs, inIndex, ((symJac, inSparsePattern, inSparseColoring, nonlinearPattern)));
3830 1064 outShared := BackendDAEUtil.setSharedSymJacs(inShared, symjacs);
3831 end addBackendDAESharedJacobianSparsePattern;
3832
3833 public function analyzeJacobian "author: PA
3834 Analyse the Jacobian to find out if the Jacobian of system of equations
3835 can be solved at compile time or runtime or if it is a non-linear system
3836 of equations."
3837 input BackendDAE.Variables vars;
3838 input BackendDAE.EquationArray eqns;
3839 input Option<list<tuple<Integer, Integer, BackendDAE.Equation>>> inTplIntegerIntegerEquationLstOption;
3840 output BackendDAE.JacobianType outJacobianType;
3841 output Boolean jacConstant "true if jac is constant, does not check rhs";
3842 algorithm
3843 (outJacobianType,jacConstant):=
3844 matchcontinue inTplIntegerIntegerEquationLstOption
3845 local
3846 list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
3847 Boolean b;
3848 BackendDAE.JacobianType jactype;
3849 case SOME(jac)
3850 algorithm
3851 //str = BackendDump.dumpJacobianStr(SOME(jac));
3852 //print("analyze Jacobian: \n" + str + "\n");
3853 3431 b := jacobianNonlinear(vars, jac);
3854 // check also if variables occur in if expressions
3855
4/4
✓ Branch 0 taken 2949 times.
✓ Branch 1 taken 482 times.
✓ Branch 5 taken 2940 times.
✓ Branch 6 taken 491 times.
 3431 (_,false) := if not b then BackendDAEUtil.traverseBackendDAEExpsEqnsWithStop(eqns,varsNotInRelations,(vars,true)) else (vars,false);
3856 //print("jac type: JAC_NONLINEAR() \n");
3857 then
3858 (BackendDAE.JAC_NONLINEAR(),false);
3859
3860 case SOME(jac)
3861 algorithm
3862
2/2
✓ Branch 1 taken 2678 times.
✓ Branch 2 taken 262 times.
 2940 true := jacobianConstant(jac);
3863 262 b := rhsConstant(vars,eqns);
3864
2/2
✓ Branch 0 taken 235 times.
✓ Branch 1 taken 27 times.
 262 jactype := if b then BackendDAE.JAC_CONSTANT() else BackendDAE.JAC_LINEAR();
3865 //print("jac type: " + if_(b,"JAC_CONSTANT()","JAC_LINEAR()") + "\n");
3866 then
3867 (jactype,true);
3868
3869 case SOME(_) then (BackendDAE.JAC_LINEAR(),false);
3870 case NONE() then (BackendDAE.JAC_NO_ANALYTIC(),false);
3871 end matchcontinue;
3872 end analyzeJacobian;
3873
3874 protected function jacobianNonlinear "author: PA
3875 Check if Jacobian indicates a non-linear system.
3876 TODO: Algorithms and array equations"
3877 input BackendDAE.Variables vars;
3878 input list<tuple<Integer, Integer, BackendDAE.Equation>> inTplIntegerIntegerEquationLst;
3879 output Boolean isNonLinear = false;
3880 protected
3881 DAE.Exp e1,e2,e;
3882 BackendDAE.Equation eq;
3883 tuple<Integer, Integer, BackendDAE.Equation> tpl;
3884 algorithm
3885
2/2
✓ Branch 0 taken 121654 times.
✓ Branch 1 taken 2949 times.
 124603 for tpl in inTplIntegerIntegerEquationLst loop
3886 121654 (_,_,eq) := tpl;
3887 isNonLinear := match eq
3888 case BackendDAE.EQUATION(exp = e1,scalar = e2)
3889 ✗ then jacobianNonlinearExp(vars, e1) or jacobianNonlinearExp(vars, e2);
3890 case BackendDAE.RESIDUAL_EQUATION(exp = e)
3891 121654 then jacobianNonlinearExp(vars, e);
3892 end match;
3893
2/2
✓ Branch 0 taken 482 times.
✓ Branch 1 taken 121172 times.
 121654 if isNonLinear then
3894 482 return;
3895 end if;
3896 end for;
3897 end jacobianNonlinear;
3898
3899 protected function jacobianNonlinearExp "author: PA
3900 Checks whether the Jacobian indicates a non-linear system.
3901 This is true if the Jacobian contains any of the variables
3902 that is solved for."
3903 input BackendDAE.Variables vars;
3904 input DAE.Exp inExp;
3905 output Boolean outBoolean;
3906 algorithm
3907 121654 (_,(_,outBoolean)) := Expression.traverseExpTopDown(inExp,traverserjacobianNonlinearExp,(vars,false));
3908 end jacobianNonlinearExp;
3909
3910 protected function traverserjacobianNonlinearExp "author: Frenkel TUD 2012-08"
3911 input DAE.Exp inExp;
3912 input tuple<BackendDAE.Variables,Boolean> tpl;
3913 output DAE.Exp outExp;
3914 output Boolean cont;
3915 output tuple<BackendDAE.Variables,Boolean> outTpl;
3916 algorithm
3917 (outExp,cont,outTpl) := matchcontinue (inExp,tpl)
3918 local
3919 BackendDAE.Variables vars;
3920 DAE.Exp e;
3921 DAE.ComponentRef cr;
3922 Boolean b;
3923 case (e as DAE.CREF(componentRef=cr),(vars,_))
3924 algorithm
3925
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 511 times.
 75929 (_::_,_) := BackendVariable.getVar(cr, vars);
3926 511 then (e,false,(vars,true));
3927
3928 case (e as DAE.CALL(path=Absyn.IDENT(name = "der"),expLst={DAE.CREF(componentRef=cr)}),(vars,_))
3929 algorithm
3930 ✗ BackendVariable.getVar(cr, vars);
3931 ✗ then (e,false,(vars,true));
3932
3933 case (e as DAE.CALL(path=Absyn.IDENT(name = "pre")),_)
3934 then (e,false,tpl);
3935
3936 case (e as DAE.CALL(path=Absyn.IDENT(name = "previous")),_)
3937 then (e,false,tpl);
3938
3939 331164 case (e,(_,b)) then (e,not b,tpl);
3940 end matchcontinue;
3941 end traverserjacobianNonlinearExp;
3942
3943 protected function jacobianConstant "author: PA
3944 Checks if Jacobian is constant, i.e. all expressions in each equation are constant."
3945 input list<tuple<Integer, Integer, BackendDAE.Equation>> inTplIntegerIntegerEquationLst;
3946 output Boolean outBoolean=true;
3947 protected
3948 DAE.Exp e1,e2, e;
3949 tuple<Integer, Integer, BackendDAE.Equation> tpl;
3950 BackendDAE.Equation eqn;
3951 algorithm
3952 /* TODO: Algorithms and ArrayEquations */
3953
3954
2/2
✓ Branch 0 taken 11996 times.
✓ Branch 1 taken 262 times.
 12258 for tpl in inTplIntegerIntegerEquationLst loop
3955 11996 eqn := Util.tuple33(tpl);
3956 outBoolean := match eqn
3957 case BackendDAE.EQUATION(exp = e1,scalar = e2)
3958 ✗ then Expression.isConst(e1) and Expression.isConst(e2);
3959 case BackendDAE.RESIDUAL_EQUATION(exp = e)
3960 11996 then Expression.isConst(e);
3961 case BackendDAE.SOLVED_EQUATION(exp = e)
3962 ✗ then Expression.isConst(e);
3963 case BackendDAE.ARRAY_EQUATION(left=e1, right=e2)
3964 ✗ then Expression.isConst(e1) and Expression.isConst(e2);
3965 case BackendDAE.COMPLEX_EQUATION(left=e1, right=e2)
3966 ✗ then Expression.isConst(e1) and Expression.isConst(e2);
3967 else false;
3968 end match;
3969
3970
2/2
✓ Branch 0 taken 2678 times.
✓ Branch 1 taken 9318 times.
 11996 if not outBoolean then
3971 break;
3972 end if;
3973
3974 end for;
3975
3976 end jacobianConstant;
3977
3978 public function isJacobianGeneric
3979 input BackendDAE.Jacobian inJac;
3980 output Boolean result;
3981 algorithm
3982 result := match inJac
3983 case BackendDAE.GENERIC_JACOBIAN() then true;
3984 else false;
3985 end match;
3986 end isJacobianGeneric;
3987
3988 protected function varsNotInRelations
3989 input output DAE.Exp exp;
3990 output Boolean cont;
3991 input output tuple<BackendDAE.Variables,Boolean> tpl;
3992 algorithm
3993 (exp,cont,tpl) := match (exp,tpl)
3994 local
3995 DAE.Exp cond,t,f,e1;
3996 BackendDAE.Variables vars;
3997 Boolean b;
3998 Absyn.Path path;
3999 list<DAE.Exp> expLst;
4000 list<DAE.Subscript> subs;
4001
4002 case (DAE.IFEXP(cond,t,f),(vars,b))
4003 algorithm
4004 // check if vars not in condition
4005
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 546 times.
 546 (_,(_,b)) := Expression.traverseExpTopDown(cond, BackendDAEUtil.getEqnsysRhsExp2, (vars,b));
4006
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 539 times.
 553 (t,(_,b)) := Expression.traverseExpTopDown(t, varsNotInRelations, (vars,b));
4007
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 539 times.
 553 (f,(_,b)) := Expression.traverseExpTopDown(f, varsNotInRelations, (vars,b));
4008
2/2
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 539 times.
 553 then (DAE.IFEXP(cond,t,f),false,(vars,b));
4009
4010 case (DAE.CALL(path=Absyn.IDENT(name = "der")),_)
4011 then (exp,true,tpl);
4012 case (DAE.CALL(path = Absyn.IDENT(name = "pre")),_)
4013 then (exp,false,tpl);
4014 case (DAE.CALL(path = Absyn.IDENT(name = "previous")),_)
4015 then (exp,false,tpl);
4016 case (DAE.CALL(path = Absyn.IDENT(name = "smooth")),_)
4017 then (exp,true,tpl);
4018 case (DAE.CALL(path = Absyn.IDENT(name = "noEvent")),_)
4019 then (exp,true,tpl);
4020 case (DAE.CALL(expLst=expLst),_)
4021 algorithm
4022 // check if vars occurs not in argument list
4023 8 (_,tpl) := Expression.traverseExpListTopDown(expLst, BackendDAEUtil.getEqnsysRhsExp2, tpl);
4024 8 then (exp,false,tpl);
4025 case (DAE.LBINARY(),_)
4026 algorithm
4027 // check if vars not in condition
4028 ✗ (_,tpl) := Expression.traverseExpTopDown(exp, BackendDAEUtil.getEqnsysRhsExp2, tpl);
4029 ✗ then (exp,false,tpl);
4030 case (DAE.LUNARY(),tpl)
4031 algorithm
4032 // check if vars not in condition
4033 ✗ (_,tpl) := Expression.traverseExpTopDown(exp, BackendDAEUtil.getEqnsysRhsExp2, tpl);
4034 ✗ then (exp,false,tpl);
4035 case (DAE.RELATION(),tpl)
4036 algorithm
4037 // check if vars not in condition
4038 ✗ (_,tpl) := Expression.traverseExpTopDown(exp, BackendDAEUtil.getEqnsysRhsExp2, tpl);
4039 ✗ then (exp,false,tpl);
4040 case (DAE.ASUB(exp=e1,sub=subs),_)
4041 algorithm
4042
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
 4 expLst := list(Expression.getSubscriptExp(sub) for sub in subs);
4043 // check if vars not in condition
4044 2 (_,tpl as (_,b)) := Expression.traverseExpTopDown(e1, varsNotInRelations, tpl);
4045
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if b then
4046 2 (_,tpl) := Expression.traverseExpListTopDown(expLst, BackendDAEUtil.getEqnsysRhsExp2, tpl);
4047 end if;
4048 2 then (exp,false,tpl);
4049 case (_,(_,b)) then (exp,b,tpl);
4050 end match;
4051 end varsNotInRelations;
4052
4053 protected function rhsConstant "author: PA
4054 Determines if the right hand sides of an equation system,
4055 represented as a BackendDAE, is constant."
4056 input BackendDAE.Variables vars;
4057 input BackendDAE.EquationArray eqns;
4058 output Boolean outBoolean;
4059 protected
4060 BackendVarTransform.VariableReplacements repl;
4061 algorithm
4062
1/2
✓ Branch 1 taken 262 times.
✗ Branch 2 not taken.
 262 if BackendEquation.equationArraySize(eqns) == 0 then
4063 outBoolean:= true;
4064 else
4065 262 repl := BackendDAEUtil.makeZeroReplacements(vars);
4066 262 (_,outBoolean,_) := BackendEquation.traverseEquationArray_WithStop(eqns,rhsConstant2,(vars,true,repl));
4067 end if;
4068 end rhsConstant;
4069
4070 protected function rhsConstant2 "Helper function to rhsConstant, traverses equation list."
4071 input BackendDAE.Equation inEq;
4072 input tuple<BackendDAE.Variables,Boolean,BackendVarTransform.VariableReplacements> inTpl;
4073 output BackendDAE.Equation outEq;
4074 output Boolean cont;
4075 output tuple<BackendDAE.Variables,Boolean,BackendVarTransform.VariableReplacements> outTpl;
4076 algorithm
4077 (outEq,cont,outTpl) := matchcontinue (inEq,inTpl)
4078 local
4079 DAE.Exp new_exp,rhs_exp,e1,e2,e;
4080 Boolean b,res;
4081 BackendDAE.Equation eqn;
4082 BackendDAE.Variables vars;
4083 BackendVarTransform.VariableReplacements repl;
4084 // check rhs for for EQUATION nodes.
4085 case (eqn as BackendDAE.EQUATION(exp = e1,scalar = e2),(vars,b,repl))
4086 algorithm
4087 339 new_exp := Expression.expSub(e1, e2);
4088 339 rhs_exp := BackendDAEUtil.getEqnsysRhsExp(new_exp, vars,NONE(),SOME(repl));
4089 339 res := Expression.isConst(rhs_exp);
4090
2/2
✓ Branch 0 taken 235 times.
✓ Branch 1 taken 104 times.
 574 then (eqn,res,(vars,b and res,repl));
4091 // check rhs for for ARRAY_EQUATION nodes. check rhs for for RESIDUAL_EQUATION nodes.
4092 case (eqn as BackendDAE.ARRAY_EQUATION(left=e1,right=e2),(vars,b,repl))
4093 algorithm
4094 ✗ new_exp := Expression.expSub(e1, e2);
4095 ✗ rhs_exp := BackendDAEUtil.getEqnsysRhsExp(new_exp, vars,NONE(),SOME(repl));
4096 ✗ res := Expression.isConst(rhs_exp);
4097 ✗ then (eqn,res,(vars,b and res,repl));
4098
4099 case (eqn as BackendDAE.COMPLEX_EQUATION(left=e1,right=e2),(vars,b,repl))
4100 algorithm
4101 ✗ new_exp := Expression.expSub(e1, e2);
4102 ✗ rhs_exp := BackendDAEUtil.getEqnsysRhsExp(new_exp, vars,NONE(),SOME(repl));
4103 ✗ res := Expression.isConst(rhs_exp);
4104 ✗ then (eqn,res,(vars,b and res,repl));
4105
4106 case (eqn as BackendDAE.RESIDUAL_EQUATION(exp = e),(vars,b,repl)) /* check rhs for for RESIDUAL_EQUATION nodes. */
4107 algorithm
4108 ✗ rhs_exp := BackendDAEUtil.getEqnsysRhsExp(e, vars,NONE(),SOME(repl));
4109 ✗ res := Expression.isConst(rhs_exp);
4110 ✗ then (eqn,res,(vars,b and res,repl));
4111
4112 ✗ case (eqn,(vars,_,repl)) then (eqn,false,(vars,false,repl));
4113 end matchcontinue;
4114 end rhsConstant2;
4115
4116 function getJacobianResiduals
4117 input BackendDAE.BackendDAE jacDAE;
4118 output list<BackendDAE.Var> diffedRes;
4119 protected
4120 BackendDAE.EqSystem syst;
4121 algorithm
4122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1819 times.
 1819 syst :: _ := jacDAE.eqs;
4123
6/6
✓ Branch 2 taken 11541 times.
✓ Branch 3 taken 5148 times.
✓ Branch 4 taken 16689 times.
✓ Branch 5 taken 1819 times.
✓ Branch 6 taken 5148 times.
✓ Branch 7 taken 1819 times.
 18508 diffedRes := list(var for var guard(BackendVariable.isRESVar(var)) in BackendVariable.varList(syst.orderedVars));
4124 end getJacobianResiduals;
4125
4126 // =============================================================================
4127 // Function detects non-linear strong component in symbolic jacobians
4128 // - non-linear components should never appear in symbolic jacobian and
4129 // indicate an singular or wrong system
4130 // - this modules stops compiling and outputs an error, otherwise we
4131 // would get error at runtime compiling
4132 // =============================================================================
4133
4134 function checkForNonLinearStrongComponents
4135 "Checks for non-linear algebraic strong compontents and break if some found."
4136 input BackendDAE.SymbolicJacobian symbolicJacobian;
4137 output Boolean result;
4138 protected
4139 BackendDAE.BackendDAE jacBDAE;
4140 String name;
4141 algorithm
4142 1820 (jacBDAE, name, _, _, _, _) := symbolicJacobian;
4143 try
4144 1820 BackendDAEUtil.mapEqSystem(jacBDAE, checkForNonLinearStrongComponents_work);
4145 result := true;
4146 else
4147 1 Error.addMessage(Error.INVALID_NONLINEAR_JACOBIAN_COMPONENT, {name});
4148 result := false;
4149 end try;
4150 end checkForNonLinearStrongComponents;
4151
4152 function checkForNonLinearStrongComponents_work
4153 input output BackendDAE.EqSystem syst;
4154 input output BackendDAE.Shared shared "unused";
4155 protected
4156 BackendDAE.StrongComponents comps;
4157 algorithm
4158 try
4159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1820 times.
 1820 BackendDAE.EQSYSTEM(matching=BackendDAE.MATCHING(comps=comps)) := syst;
4160
2/2
✓ Branch 0 taken 16086 times.
✓ Branch 1 taken 1819 times.
 17905 for comp in comps loop
4161 () := match comp
4162 case BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_NONLINEAR()) algorithm
4163 ✗ if Flags.isSet(Flags.JAC_DUMP) then
4164 ✗ print("[symjacdump] Following strong component represents a nonlinear symbolic jacobian:\n" + BackendDump.printComponent(comp, SOME(syst)) + "\n");
4165 end if;
4166 ✗ then fail();
4167 case BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_NO_ANALYTIC())algorithm
4168 ✗ if Flags.isSet(Flags.JAC_DUMP) then
4169 ✗ print("[symjacdump] Following strong component represents a no symbolic jacobian:\n" + BackendDump.printComponent(comp, SOME(syst)) + "\n");
4170 end if;
4171 ✗ then fail();
4172 case BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_GENERIC())algorithm
4173 ✗ if Flags.isSet(Flags.JAC_DUMP) then
4174 ✗ print("[symjacdump] Following strong component represents a generic jacobian:\n" + BackendDump.printComponent(comp, SOME(syst)) + "\n");
4175 end if;
4176 ✗ then fail();
4177 case BackendDAE.TORNSYSTEM(linear=false)algorithm
4178
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if Flags.isSet(Flags.JAC_DUMP) then
4179 ✗ print("[symjacdump] Following (torn) strong component represents a nonlinear symbolic jacobian:\n" + BackendDump.printComponent(comp, SOME(syst)) + "\n");
4180 end if;
4181 1 then fail();
4182 else ();
4183 end match;
4184 end for;
4185 else
4186 1 fail();
4187 end try;
4188 end checkForNonLinearStrongComponents_work;
4189
4190
4191 public function getFixedStatesForSelfdependentSets
4192 " author: kabdelhak
4193 Returns states to fix for initial problem in the case of selfdependent dynamic state sets"
4194 input BackendDAE.StateSet stateSet;
4195 input list<BackendDAE.Var> unfixedStates;
4196 input Integer toFix;
4197 output list<BackendDAE.Var> statesToFix;
4198 protected
4199 list<tuple<Integer,BackendDAE.Var>> nonlinearCountLst = {};
4200 algorithm
4201 _:= match stateSet.jacobian
4202 local
4203 BackendDAE.SymbolicJacobian sJac;
4204 BackendDAE.BackendDAE dae;
4205 String matrixName;
4206 case BackendDAE.GENERIC_JACOBIAN(jacobian=SOME(sJac)) algorithm
4207 ✗ (dae,matrixName,_, _, _,_) := sJac;
4208 ✗ for var in unfixedStates loop
4209 ✗ nonlinearCountLst := getNonlinearStateCount(var,unfixedStates,dae,matrixName)::nonlinearCountLst;
4210 end for;
4211 then 0;
4212 end match;
4213 ✗ statesToFix := fixedVarsFromNonlinearCount(nonlinearCountLst, toFix);
4214 end getFixedStatesForSelfdependentSets;
4215
4216 protected function getNonlinearStateCount
4217 input BackendDAE.Var state;
4218 input list<BackendDAE.Var> diffVars;
4219 input BackendDAE.BackendDAE dae;
4220 input String matrixName;
4221 output tuple<Integer,BackendDAE.Var> outTpl;
4222 protected
4223 algorithm
4224 outTpl:=match dae
4225 local
4226 BackendDAE.EqSystems systs;
4227 tuple<BackendDAE.Var,list<BackendDAE.Var>,Integer,String> tpl;
4228 BackendDAE.Var outState;
4229 Integer nonlinearCount = 0;
4230 case BackendDAE.DAE(eqs=systs) algorithm
4231 ✗ tpl := (state,diffVars,nonlinearCount,matrixName);
4232 ✗ for syst in systs loop
4233 _:= match syst
4234 local
4235 BackendDAE.EquationArray eqnarray;
4236
4237 case BackendDAE.EQSYSTEM(_,eqnarray,_,_,_,_,_,_) algorithm
4238 ✗ tpl := BackendEquation.traverseEquationArray(eqnarray,getNonlinearStateCount0,tpl);
4239 then 0;
4240 end match;
4241 end for;
4242 ✗ (outState,_,nonlinearCount,_) := tpl;
4243 ✗ then (nonlinearCount,outState);
4244 end match;
4245
4246 end getNonlinearStateCount;
4247
4248 protected function getNonlinearStateCount0
4249 input BackendDAE.Equation inEq;
4250 input tuple<BackendDAE.Var,list<BackendDAE.Var>,Integer,String> inTpl;
4251 output BackendDAE.Equation outEq;
4252 output tuple<BackendDAE.Var,list<BackendDAE.Var>,Integer,String> outTpl;
4253 algorithm
4254 outEq := inEq;
4255 outTpl := match inEq
4256 local
4257 DAE.Exp exp, diffExp;
4258 BackendDAE.Var state;
4259 list<BackendDAE.Var> diffVars;
4260 Integer nonlinearCount;
4261 String matrixName;
4262 DAE.ComponentRef seedVar;
4263 case BackendDAE.EQUATION(scalar=exp) algorithm
4264 ✗ (state,diffVars,nonlinearCount,matrixName) := inTpl;
4265 // Differentiate equation to look for nonlinear dependencies
4266 ✗ seedVar := Differentiate.createSeedCrefName(BackendVariable.varCref(state),matrixName);
4267 ✗ diffExp := Differentiate.differentiateExpSolve(exp,seedVar,NONE());
4268 ✗ for var in diffVars loop
4269 ✗ if not ComponentReferenceBasics.crefEqual(var.varName, state.varName) and Expression.expContains(diffExp,Expression.crefExp(var.varName)) then
4270 // Heuristic to punish vars with a value of zero
4271 ✗ if Expression.isZero(BackendVariable.varStartValue(var)) then
4272 ✗ nonlinearCount := nonlinearCount + 2;
4273 else
4274 ✗ nonlinearCount := nonlinearCount + 1;
4275 end if;
4276 end if;
4277 end for;
4278 ✗ then (state,diffVars,nonlinearCount,matrixName);
4279 end match;
4280 end getNonlinearStateCount0;
4281
4282 protected function fixedVarsFromNonlinearCount
4283 input list<tuple<Integer,BackendDAE.Var>> tplLst;
4284 input Integer toFix;
4285 output list<BackendDAE.Var> fixedVars = {};
4286 protected
4287 list<tuple<Integer,BackendDAE.Var>> sortedTplLst, strippedTplLst;
4288 BackendDAE.Var fixVar;
4289 Integer fixInt;
4290 algorithm
4291 ✗ for tpl in tplLst loop
4292 (fixInt,fixVar) := tpl;
4293 end for;
4294 // Sort by nonlinear count and take first N states to fix
4295 ✗ sortedTplLst := List.sort(tplLst, Util.compareTupleIntGt);
4296 ✗ strippedTplLst := List.firstN(sortedTplLst,toFix);
4297 ✗ for tpl in strippedTplLst loop
4298 ✗ (_,fixVar) := tpl;
4299 ✗ fixVar.values := DAEUtil.setFixedAttr(fixVar.values,SOME(DAE.BCONST(true)));
4300 fixedVars := fixVar::fixedVars;
4301 end for;
4302 end fixedVarsFromNonlinearCount;
4303
4304 protected function stripPartialDerNonlinearPattern
4305 "kabdelhak: this function strips a jacobian residual down to the original
4306 residual variable to get the equation mapping right. Used for nonlinear
4307 pattern analysis."
4308 input output BackendDAE.NonlinearPattern pat;
4309 protected
4310 BackendDAE.NonlinearPatternCrefs pat_cref, pat_crefT;
4311 list<DAE.ComponentRef> v1, v2;
4312 Integer index;
4313 algorithm
4314 1819 (pat_cref, pat_crefT, (v1, v2), index) := pat;
4315
4/4
✓ Branch 0 taken 5207 times.
✓ Branch 1 taken 1819 times.
✓ Branch 2 taken 5207 times.
✓ Branch 3 taken 1819 times.
 7026 pat_cref := list(stripPartialDer(cref_tpl) for cref_tpl in pat_cref);
4316
4/4
✓ Branch 0 taken 5148 times.
✓ Branch 1 taken 1819 times.
✓ Branch 2 taken 5148 times.
✓ Branch 3 taken 1819 times.
 6967 pat_crefT := list(stripPartialDer(cref_tpl) for cref_tpl in pat_crefT);
4317
4/4
✓ Branch 0 taken 5207 times.
✓ Branch 1 taken 1819 times.
✓ Branch 2 taken 5207 times.
✓ Branch 3 taken 1819 times.
 7026 v1 := list(stripPartialDerWork(v) for v in v1);
4318
4/4
✓ Branch 0 taken 5148 times.
✓ Branch 1 taken 1819 times.
✓ Branch 2 taken 5148 times.
✓ Branch 3 taken 1819 times.
 6967 v2 := list(stripPartialDerWork(v) for v in v2);
4319 1819 pat := (pat_cref, pat_crefT, (v1, v2), index);
4320 end stripPartialDerNonlinearPattern;
4321
4322 protected function stripPartialDer
4323 input output BackendDAE.NonlinearPatternCref cref_tpl;
4324 protected
4325 DAE.ComponentRef cref;
4326 list<DAE.ComponentRef> dependencies;
4327 algorithm
4328 10355 (cref, dependencies) := cref_tpl;
4329 10355 (cref, _) := stripPartialDerWork(cref);
4330
4/4
✓ Branch 0 taken 2864 times.
✓ Branch 1 taken 10355 times.
✓ Branch 2 taken 2864 times.
✓ Branch 3 taken 10355 times.
 13219 dependencies := list(stripPartialDerWork(dep) for dep in dependencies);
4331 10355 cref_tpl := (cref, dependencies);
4332 end stripPartialDer;
4333
4334 protected function stripPartialDerWork
4335 input output DAE.ComponentRef cref;
4336 output Boolean strip;
4337 algorithm
4338 (cref, strip) := match cref
4339 local
4340 DAE.ComponentRef cr;
4341
4342 case DAE.CREF_IDENT() guard(StringUtil.startsWith(cref.ident, "$pDER")) then (cref, true);
4343
4344 case DAE.CREF_QUAL() guard(StringUtil.startsWith(cref.ident, "$pDER")) then (cref, true);
4345
4346 case DAE.CREF_QUAL() algorithm
4347 32791 (cr, strip) := stripPartialDerWork(cref.componentRef);
4348
2/2
✓ Branch 0 taken 11774 times.
✓ Branch 1 taken 21017 times.
 32791 if strip then
4349 11774 cr := DAE.CREF_IDENT(cref.ident, cref.identType, cref.subscriptLst);
4350 else
4351 21017 cr := DAE.CREF_QUAL(cref.ident, cref.identType, cref.subscriptLst, cr);
4352 end if;
4353 then (cr, false);
4354
4355 else (cref, false);
4356 end match;
4357 end stripPartialDerWork;
4358
4359 // =============================================================================
4360 // [ASSC] section for analytical to symbolical singularity transformation
4361 //
4362 // Generates linear jacobian
4363 // =============================================================================
4364 public
4365 type LinearJacobianRow = UnorderedMap<Integer, Real>;
4366 type LinearJacobianRhs = array<.DAE.Exp>;
4367 type LinearJacobianInd = array<tuple<Integer, Integer>>;
4368
4369 uniontype LinearJacobian
4370 record LINEAR_REAL_JACOBIAN
4371 array<LinearJacobianRow> rows "all loop variables entries";
4372 LinearJacobianRhs rhs "the expression containing all non loop variable entries";
4373 LinearJacobianInd ind "equation indices <array, scalar>";
4374 array<Boolean> eq_marks "changed equations";
4375 end LINEAR_REAL_JACOBIAN;
4376
4377 public function toString
4378 input SymbolicJacobian.LinearJacobian linJac;
4379 input String heading = "";
4380 output String str;
4381 algorithm
4382 2 str := "######################################################\n" +
4383 " LinearJacobian sparsity pattern: " + heading + "\n" +
4384 "######################################################\n" +
4385 "(scal_idx|arr_idx|changed) [var_index, value] || RHS_EXPRESSION\n";
4386
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 8 for idx in 1:arrayLength(linJac.rows) loop
4387 4 str := str + rowToString(linJac.rows[idx], linJac.rhs[idx], linJac.ind[idx], linJac.eq_marks[idx]);
4388 end for;
4389 2 str := str + "\n";
4390 end toString;
4391
4392 protected function rowToString
4393 input SymbolicJacobian.LinearJacobianRow row;
4394 input DAE.Exp rhs;
4395 input tuple<Integer, Integer> indices;
4396 input Boolean changed;
4397 output String str;
4398 protected
4399 Integer i_arr, i_scal, index;
4400 Real value;
4401 list<tuple<Integer, Real>> row_lst = UnorderedMap.toList(row);
4402 algorithm
4403 4 (i_arr, i_scal) := indices;
4404
2/2
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 1 time.
 7 str := "(" + intString(i_arr) + "|" + intString(i_scal) + "|" + boolString(changed) +"): ";
4405
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
 4 if listEmpty(row_lst) then
4406 1 str := str + "EMPTY ROW ";
4407 else
4408
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
 9 for element in row_lst loop
4409 6 (index, value) := element;
4410 6 str := str + "[" + intString(index) + "|" + realString(value) + "] ";
4411 end for;
4412 end if;
4413 4 str := str + " || RHS: " + ExpressionBasics.printExpStr(ExpressionSimplify.simplify(rhs)) + "\n";
4414 end rowToString;
4415
4416 public function generate
4417 "author: kabdelhak FHB 03-2021
4418 Generates a jacobian from algebraic loop equations which are linear
4419 w.r.t. all loopVars. Fails if these criteria are not met."
4420 input list<tuple<BackendDAE.Equation, tuple<Integer, Integer>>> loopEqs;
4421 input list<tuple<BackendDAE.Var, Integer>> loopVars;
4422 input array<Integer> ass1;
4423 output LinearJacobian linJac;
4424 protected
4425 Integer eqn_index = 1, var_index;
4426 Real constReal;
4427 LinearJacobianRow row;
4428 list<LinearJacobianRow> tmp_mat = {};
4429 list<DAE.Exp> tmp_rhs = {};
4430 list<tuple<Integer, Integer>> tmp_idx = {};
4431 BackendDAE.Equation eqn;
4432 tuple<Integer, Integer> index;
4433 Integer scal_idx;
4434 BackendDAE.Var var;
4435 DAE.Exp res, pDer;
4436 BackendVarTransform.VariableReplacements varRep;
4437
4438 // Helper functions to either have integer or real valued coefficients
4439 evaluateFunc eFunc = if Flags.getConfigBool(Flags.REAL_ASSC) then Expression.getEvaluatedConstReal else intWrapperFunc;
4440
4441 partial function evaluateFunc
4442 input DAE.Exp e;
4443 output Real v;
4444 end evaluateFunc;
4445
4446 function intWrapperFunc extends evaluateFunc;
4447 algorithm
4448 135536 v := intReal(Expression.getEvaluatedConstInteger(e));
4449 end intWrapperFunc;
4450
4451 algorithm
4452 /* Add a replacement rule var->0 for each loopVar, so that the RHS can be determined afterwards */
4453 908 varRep := BackendVarTransform.emptyReplacements();
4454
2/2
✓ Branch 0 taken 8460 times.
✓ Branch 1 taken 908 times.
 9368 for loopVar in loopVars loop
4455 8460 (var, _) := loopVar;
4456 8460 varRep := BackendVarTransform.addReplacement(varRep, BackendVariable.varCref(var), DAE.ICONST(0), NONE());
4457 end for;
4458
4459 /* Loop over all equations and create residual expression. */
4460
2/2
✓ Branch 0 taken 7059 times.
✓ Branch 1 taken 898 times.
 7957 for loopEq in loopEqs loop
4461 7059 row := UnorderedMap.new<Real>(Util.id, intEq);
4462 7059 (eqn, index) := loopEq;
4463 7059 res := BackendEquation.createResidualExp(eqn);
4464 /* Loop over all variables and differentiate residual expression for each. */
4465 try
4466
2/2
✓ Branch 0 taken 137012 times.
✓ Branch 1 taken 2226 times.
 139238 for loopVar in loopVars loop
4467 137012 (var, var_index) := loopVar;
4468 137012 pDer := Differentiate.differentiateExpSolve(res, BackendVariable.varCref(var), NONE());
4469 135536 (pDer, _) := ExpressionSimplify.simplify(pDer);
4470
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 135536 times.
 135536 constReal := eFunc(pDer);
4471
2/2
✓ Branch 0 taken 8080 times.
✓ Branch 1 taken 124109 times.
 132189 if not realEq(constReal, 0.0) then
4472 8080 UnorderedMap.add(var_index, constReal, row);
4473 end if;
4474 end for;
4475 /*
4476 Save the full row.
4477 - row entries
4478 - rhs
4479 - equation index
4480 Perform var replacements, multiply by -1 and simplify for rhs.
4481 NOTE: Multiplication with -1 is not really necessary for the
4482 conversion of analytical to structural singularity, but
4483 would be necessary if used for anything else.
4484 */
4485 2226 res := BackendVarTransform.replaceExp(res, varRep, NONE());
4486 2226 tmp_mat := row :: tmp_mat;
4487 4452 tmp_rhs := ExpressionSimplify.simplify(DAE.BINARY(DAE.ICONST(-1), DAE.MUL(DAE.T_UNKNOWN_DEFAULT), res)) :: tmp_rhs;
4488 tmp_idx := index :: tmp_idx;
4489
4490 /* set var as matched so that it can be chosen as pivot element for gaussian elimination */
4491 (_, scal_idx) := index;
4492 eqn_index := eqn_index + 1;
4493 else
4494 /*
4495 Differentiation not possible or not convertible to a real.
4496 Purposely fails.
4497 */
4498 end try;
4499 end for;
4500 /* convert and store all data */
4501 898 linJac := LINEAR_REAL_JACOBIAN(
4502 rows = listArray(tmp_mat),
4503 rhs = listArray(tmp_rhs),
4504 ind = listArray(tmp_idx),
4505 eq_marks = arrayCreate(listLength(tmp_mat), false)
4506 );
4507 end generate;
4508
4509 public function emptyOrSingle
4510 "author: kabdelhak FHB 03-2021
4511 Returns true if the linear real jacobian is empty or has only one single row."
4512 input LinearJacobian linJac;
4513 output Boolean empty = (arrayLength(linJac.rows) < 2)
4514 and (arrayLength(linJac.rhs) < 2)
4515 and (arrayLength(linJac.ind) < 2)
4516 and (arrayLength(linJac.eq_marks) < 2);
4517 end emptyOrSingle;
4518
4519 public function solve
4520 "author: kabdelhak FHB 03-2021
4521 Performs a gaussian elimination algorithm on the jacobian without reducing the
4522 pivot elements to one to maintain the integer structure. This guarantees that
4523 no numerical errors can occur and analytical singularities will be detected.
4524 Also keeps track of the RHS for later equation replacement.
4525
4526 Performs gaussian elimination for one pivot row and all following rows to reduce.
4527 new_row = old_row * pivot_element - pivot_row * row_element
4528 Example:
4529 pivot idx: 2, because the first is zero
4530 pivot row: | 0 -1 -4 |
4531 row-to change: | -3 2 3 |
4532 new_row: | 3 0 5 |"
4533 input output LinearJacobian linJac;
4534 protected
4535 Integer col_index;
4536 Real piv_value, row_value;
4537 algorithm
4538 /*
4539 Gaussian Algorithm without rearranging rows.
4540 */
4541
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 333 times.
 2741 for i in 1:arrayLength(linJac.rows) loop
4542 try
4543 /*
4544 no pivot element can be chosen?
4545 jump over all manipulations, nothing to do
4546 */
4547 2075 (col_index, piv_value) := getPivot(linJac.rows[i]);
4548
4549 //updatePivotRow(linJac.rows[i], piv_value);
4550 // ToDo: updating the pivot row would also need an update for the rhs!
4551
4552
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2074 times.
✓ Branch 2 taken 1742 times.
✓ Branch 3 taken 332 times.
 17927 for j in i+1:arrayLength(linJac.rows) loop
4553 13779 row_value := UnorderedMap.getOrDefault(col_index, linJac.rows[j], 0.0);
4554
2/2
✓ Branch 0 taken 1266 times.
✓ Branch 1 taken 12513 times.
 13779 if not realEq(row_value, 0.0) then
4555 // set row to processed and perform pivot step
4556 1266 linJac.eq_marks[j] := true;
4557 1266 solveRow(linJac.rows[i], linJac.rows[j], piv_value, row_value);
4558 //perform multiplication inside? use simplification of multiplication afterwards?
4559 1266 linJac.rhs[j] := DAE.BINARY(
4560 DAE.BINARY(linJac.rhs[j], DAE.MUL(DAE.T_REAL_DEFAULT), DAE.RCONST(piv_value)), // row_rhs * piv_elem
4561 DAE.SUB(DAE.T_REAL_DEFAULT), // -
4562 DAE.BINARY(linJac.rhs[i], DAE.MUL(DAE.T_REAL_DEFAULT), DAE.RCONST(row_value)) // piv_rhs * row_elem
4563 );
4564 end if;
4565 end for;
4566 else
4567 /* no pivot element, nothing to do */
4568 end try;
4569 end for;
4570 end solve;
4571
4572 public function solveRow
4573 "author: kabdelhak FHB 03-2021
4574 performs one single row update : new_row = old_row * pivot_element - pivot_row * row_element"
4575 input LinearJacobianRow pivot_row;
4576 input LinearJacobianRow row;
4577 input Real piv_value;
4578 input Real row_value;
4579 protected
4580 Integer idx;
4581 Real val, diag_val;
4582 algorithm
4583 // update all elements that are in the pivot row
4584
2/2
✓ Branch 2 taken 4520 times.
✓ Branch 3 taken 1266 times.
 5786 for idx in UnorderedMap.keyList(pivot_row) loop
4585 () := match (UnorderedMap.get(idx, row), UnorderedMap.get(idx, pivot_row))
4586
4587 // row to be updated has and element at this position
4588 case (SOME(val), SOME(diag_val)) algorithm
4589 2216 val := val * piv_value - diag_val * row_value;
4590
2/2
✓ Branch 0 taken 2077 times.
✓ Branch 1 taken 139 times.
 2216 if realAbs(val) < 1e-12 then
4591 /* delete element if zero */
4592 2077 UnorderedMap.remove(idx, row);
4593 else
4594 139 UnorderedMap.add(idx, val, row);
4595 end if;
4596 then ();
4597
4598 // row to be updated does not have an element at this position
4599 case (NONE(), SOME(diag_val)) algorithm
4600 2304 UnorderedMap.add(idx, -diag_val * row_value, row);
4601 then ();
4602
4603 else algorithm
4604 ✗ Error.terminate(getInstanceName() + " key does not have an element in pivot row.", sourceInfo());
4605 then ();
4606 end match;
4607 end for;
4608
4609 // update all row elements that are not in pivot row
4610
2/2
✓ Branch 2 taken 4890 times.
✓ Branch 3 taken 1266 times.
 6156 for idx in UnorderedMap.keyList(row) loop
4611 () := match (UnorderedMap.get(idx, row), UnorderedMap.get(idx, pivot_row))
4612 case (SOME(val), NONE()) algorithm
4613 2447 val := val * piv_value;
4614 2447 UnorderedMap.add(idx, val, row);
4615 then ();
4616 else ();
4617 end match;
4618 end for;
4619 end solveRow;
4620
4621 public function updatePivotRow
4622 "author: kabdelhak FHB 03-2021
4623 updates the pivot row by dividing everything by its pivot value"
4624 input LinearJacobianRow pivot_row;
4625 input Real piv_value;
4626 protected
4627 Real value;
4628 algorithm
4629 ✗ if not realEq(piv_value, 1.0) then
4630 ✗ for idx in UnorderedMap.keyList(pivot_row) loop
4631 ✗ value := UnorderedMap.getOrFail(idx, pivot_row);
4632 ✗ UnorderedMap.add(idx, value/piv_value, pivot_row);
4633 end for;
4634 end if;
4635 end updatePivotRow;
4636
4637 protected function getPivot
4638 "author: kabdelhak FHB 03-2021
4639 Returns the first element that can be chosen as pivot, fails if none can be chosen."
4640 input LinearJacobianRow pivot_row;
4641 output Integer idx;
4642 output Real value;
4643 algorithm
4644
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 2074 times.
 2075 if Vector.isEmpty(pivot_row.keys) then
4645 /* singular row */
4646 1 fail();
4647 else
4648 2074 idx := UnorderedMap.firstKey(pivot_row);
4649 2074 value := UnorderedMap.getOrFail(idx, pivot_row);
4650 end if;
4651 end getPivot;
4652
4653 public function resolveASSC
4654 "author: kabdelhak FHB 03-2021
4655 Resolves analytical singularities by replacing the equations with
4656 zero rows in the jacobian with new equations. Needs preceeding
4657 solving of the linear real jacobian."
4658 input LinearJacobian linJac;
4659 input output array<Integer> ass1;
4660 input output array<Integer> ass2;
4661 input output BackendDAE.EqSystem syst;
4662 input Boolean init;
4663 protected
4664 Integer i_arr, i_scal;
4665 DAE.Exp lhs, rhs;
4666 BackendDAE.Equation newEqn;
4667 list<Integer> updateList_arr = {};
4668 array<list<Integer>> mapEqnIncRow;
4669 array<Integer> mapIncRowEqn;
4670 BackendDAE.IndexType indexType;
4671 Boolean fullASSC = Flags.getConfigBool(Flags.FULL_ASSC);
4672 algorithm
4673
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 333 times.
✓ Branch 2 taken 333 times.
✗ Branch 3 not taken.
 2741 for r in 1:arrayLength(linJac.rows) loop
4674 /*
4675 check if row has been changed
4676 for now also only resolve singularities and not replace full loop
4677 otherwise it sometimes leads to mixed determined systems
4678 */
4679
5/6
✓ Branch 1 taken 851 times.
✓ Branch 2 taken 1224 times.
✓ Branch 5 taken 850 times.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 850 times.
 2075 if linJac.eq_marks[r] and (UnorderedMap.isEmpty(linJac.rows[r]) or fullASSC) then
4680 1 (i_arr, i_scal) := linJac.ind[r];
4681 /* remove assignments */
4682 1 ass2[ass1[i_scal]] := -1;
4683 1 ass1[i_scal] := -1;
4684
4685 /* replace equation */
4686 1 rhs := ExpressionSimplify.simplify(linJac.rhs[r]);
4687 1 lhs := generateLHSfromList(
4688 row_indices = UnorderedMap.keyArray(linJac.rows[r]),
4689 row_values = UnorderedMap.valueArray(linJac.rows[r]),
4690 vars = syst.orderedVars
4691 );
4692 1 newEqn := BackendEquation.generateEquation(lhs, rhs);
4693
4694 /* dump replacements */
4695
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 1 if Flags.isSet(Flags.DUMP_ASSC) or (Flags.isSet(Flags.BLT_DUMP) and UnorderedMap.isEmpty(linJac.rows[r])) then
4696 1 print("[ASSC] The equation: " + BackendDump.equationString(BackendEquation.get(syst.orderedEqs, i_arr)) + "\n");
4697 1 print("[ASSC] Gets replaced by equation: " + BackendDump.equationString(newEqn) + "\n");
4698 end if;
4699
4700 1 syst.orderedEqs := BackendEquation.setAtIndex(syst.orderedEqs, i_arr, newEqn);
4701 updateList_arr := i_arr :: updateList_arr;
4702 end if;
4703 end for;
4704 /*
4705 update adjacency matrix and transposed adjacency matrix
4706 */
4707
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 332 times.
 333 if not listEmpty(updateList_arr) then
4708 try
4709 /* scalar = true */
4710
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 time.
 1 SOME((mapEqnIncRow, mapIncRowEqn, indexType, true, _)) := syst.mapping;
4711 1 syst := BackendDAEUtil.updateAdjacencyMatrixScalar(syst, indexType, NONE(), updateList_arr, mapEqnIncRow, mapIncRowEqn, false);
4712 else
4713 /*
4714 scalar = false,
4715 should never occur, just to have a fallback option if someone wants to use this algorithm somewhere else
4716 */
4717 ✗ syst := BackendDAEUtil.updateAdjacencyMatrix(syst, BackendDAE.SOLVABLE(), NONE(), updateList_arr, false);
4718 end try;
4719 end if;
4720
4721
3/6
✓ Branch 0 taken 332 times.
✓ Branch 1 taken 1 time.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 333 if not listEmpty(updateList_arr) and not Flags.isSet(Flags.DUMP_ASSC) and Flags.isSet(Flags.BLT_DUMP) then
4722 ✗ print("--- Some equations have been changed, for more information please use -d=dumpASSC.---\n\n");
4723 end if;
4724 end resolveASSC;
4725
4726 protected function generateLHSfromList
4727 "author: kabdelhak FHB 03-2021
4728 Generates the LHS expression from a flattened linear real jacobian row.
4729 Only used for full replacement of causalized loop."
4730 input array<Integer> row_indices;
4731 input array<Real> row_values;
4732 input BackendDAE.Variables vars;
4733 output DAE.Exp lhs;
4734 protected
4735 Integer length = arrayLength(row_indices);
4736 algorithm
4737 // add first expression
4738
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if length == 0 then
4739 lhs := DAE.RCONST(0.0);
4740 else
4741 ✗ lhs := DAE.BINARY(
4742 DAE.RCONST(row_values[1]),
4743 DAE.MUL(DAE.T_REAL_DEFAULT),
4744 BackendVariable.varExp(BackendVariable.getVarAt(vars, row_indices[1]))
4745 );
4746 end if;
4747
4748 // add subsequent expressions
4749
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 for i in 2:arrayLength(row_indices) loop
4750 ✗ lhs := DAE.BINARY(lhs, DAE.ADD(DAE.T_REAL_DEFAULT), DAE.BINARY(
4751 DAE.RCONST(row_values[i]),
4752 DAE.MUL(DAE.T_REAL_DEFAULT),
4753 BackendVariable.varExp(BackendVariable.getVarAt(vars, row_indices[i]))
4754 ));
4755 end for;
4756 end generateLHSfromList;
4757
4758 public function anyChanges
4759 "author: kabdelhak FHB 03-2021
4760 Returns true if any row of the jacobian got changed during gaussian elimination."
4761 input LinearJacobian linJac;
4762 output Boolean changed = false;
4763 algorithm
4764
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 175 times.
 958 for i in 1:arrayLength(linJac.eq_marks) loop
4765
2/2
✓ Branch 1 taken 108 times.
✓ Branch 2 taken 608 times.
 716 if linJac.eq_marks[i] then
4766 changed := true;
4767 108 return;
4768 end if;
4769 end for;
4770 end anyChanges;
4771 end LinearJacobian;
4772
4773 annotation(__OpenModelica_Interface="backend");
4774 end SymbolicJacobian;
4775