Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 4.0% 24 / 0 / 607
Functions: 16.7% 4 / 0 / 24
Branches: 1.7% 6 / 0 / 355
OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSystem.c
Line Branch Exec Source
1 /*
2 * This file belongs to the OpenModelica Run-Time System
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings
5 * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights
6 * reserved.
7 *
8 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
9 * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY
10 * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
11 * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL
12 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
13 *
14 * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License
15 * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs:
16 * http://www.openmodelica.org or https://github.com/OpenModelica/ or
17 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU
18 * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW
19 * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY
23 * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF
24 * OSMC-PL.
25 *
26 */
27
28 /*! \file nonlinearSystem.c
29 */
30
31 #include <math.h>
32 #include <string.h>
33
34 #include "../jacobian_util.h"
35 #include "../../util/simulation_options.h"
36 #include "../../util/omc_error.h"
37 #include "../../util/omc_file.h"
38 #include "nonlinearSystem.h"
39 #include "nonlinearValuesList.h"
40 #if !defined(OMC_MINIMAL_RUNTIME)
41 #include "kinsolSolver.h"
42 #include "kinsol_b.h"
43 #include "nonlinearSolverHybrd.h"
44 #include "nonlinearSolverNewton.h"
45 #include "newtonIteration.h"
46 #include "newton_diagnostics.h"
47 #endif
48 #include "nonlinearSolverHomotopy.h"
49 #include "../options.h"
50 #include "../simulation_info_json.h"
51 #include "../simulation_runtime.h"
52 #include "model_help.h"
53
54 int check_nonlinear_solution(DATA *data, int printFailingSystems, int sysNumber);
55
56 extern int init_lambda_steps;
57
58 struct dataMixedSolver
59 {
60 void* newtonHomotopyData;
61 void* hybridData;
62 };
63
64 #if !defined(OMC_MINIMAL_RUNTIME)
65 #include "../../util/write_csv.h"
66
67 /*! \fn int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData)
68 *
69 * This function initializes csv files for analysis propose.
70 *
71 * \param [ref] [data]
72 * \param [ref] [systemData]
73 */
74 ✗ int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData)
75 {
76 ✗ struct csvStats* stats = (struct csvStats*) malloc(sizeof(struct csvStats));
77 char buffer[100];
78 ✗ sprintf(buffer, "%s_NLS%dStatsCall.csv", data->modelData->modelFilePrefix, (int)systemData->equationIndex);
79 ✗ stats->callStats = omc_write_csv_init(buffer, ',', '"');
80
81 ✗ sprintf(buffer, "%s_NLS%dStatsIter.csv", data->modelData->modelFilePrefix, (int)systemData->equationIndex);
82 ✗ stats->iterStats = omc_write_csv_init(buffer, ',', '"');
83
84 ✗ systemData->csvData = stats;
85
86 ✗ return 0;
87 }
88 #else
89 ✗ int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData)
90 {
91 ✗ fprintf(stderr, "initializeNLScsvData not implemented for OMC_MINIMAL_RUNTIME");
92 ✗ abort();
93 }
94 #endif
95
96 #if !defined(OMC_MINIMAL_RUNTIME)
97 /*! \fn int print_csvLineCallStatsHeader(OMC_WRITE_CSV* csvData)
98 *
99 * This function initializes csv files for analysis propose.
100 *
101 * \param [ref] [data]
102 * \param [ref] [systemData]
103 */
104 ✗ int print_csvLineCallStatsHeader(OMC_WRITE_CSV* csvData)
105 {
106 char buffer[1024];
107 ✗ buffer[0] = 0;
108
109 /* number of call */
110 sprintf(buffer,"numberOfCall");
111 ✗ omc_write_csv(csvData, buffer);
112 ✗ fputc(csvData->seperator,csvData->handle);
113
114 /* simulation time */
115 sprintf(buffer,"simulationTime");
116 ✗ omc_write_csv(csvData, buffer);
117 ✗ fputc(csvData->seperator,csvData->handle);
118
119 /* solving iterations */
120 sprintf(buffer,"iterations");
121 ✗ omc_write_csv(csvData, buffer);
122 ✗ fputc(csvData->seperator,csvData->handle);
123
124 /* solving fCalls */
125 sprintf(buffer,"numberOfFunctionCall");
126 ✗ omc_write_csv(csvData, buffer);
127 ✗ fputc(csvData->seperator,csvData->handle);
128
129 /* solving Time */
130 sprintf(buffer,"solvingTime");
131 ✗ omc_write_csv(csvData, buffer);
132 ✗ fputc(csvData->seperator,csvData->handle);
133
134 /* solving Time */
135 sprintf(buffer,"solvedSystem");
136 ✗ omc_write_csv(csvData, buffer);
137
138 /* finish line */
139 ✗ fputc('\n',csvData->handle);
140
141 ✗ return 0;
142 }
143
144 /*! \fn int print_csvLineCallStats(OMC_WRITE_CSV* csvData)
145 *
146 * This function initializes csv files for analysis propose.
147 *
148 * \param [ref] [csvData]
149 * \param [in] [number of calls]
150 * \param [in] [simulation time]
151 * \param [in] [iterations]
152 * \param [in] [number of function call]
153 * \param [in] [solving time]
154 * \param [in] [solved system]
155 */
156 ✗ int print_csvLineCallStats(OMC_WRITE_CSV* csvData, int num, double time,
157 int iterations, int fCalls, double solvingTime,
158 NLS_SOLVER_STATUS solved)
159 {
160 char buffer[1024];
161
162 /* number of call */
163 sprintf(buffer, "%d", num);
164 ✗ omc_write_csv(csvData, buffer);
165 ✗ fputc(csvData->seperator,csvData->handle);
166
167 /* simulation time */
168 sprintf(buffer, "%g", time);
169 ✗ omc_write_csv(csvData, buffer);
170 ✗ fputc(csvData->seperator,csvData->handle);
171
172 /* solving iterations */
173 sprintf(buffer, "%d", iterations);
174 ✗ omc_write_csv(csvData, buffer);
175 ✗ fputc(csvData->seperator,csvData->handle);
176
177 /* solving fCalls */
178 sprintf(buffer, "%d", fCalls);
179 ✗ omc_write_csv(csvData, buffer);
180 ✗ fputc(csvData->seperator,csvData->handle);
181
182 /* solving Time */
183 sprintf(buffer, "%f", solvingTime);
184 ✗ omc_write_csv(csvData, buffer);
185 ✗ fputc(csvData->seperator,csvData->handle);
186
187 /* solved system */
188 ✗ sprintf(buffer, "%s", (solved == NLS_SOLVED || solved == NLS_SOLVED_LESS_ACCURACY)?"TRUE":"FALSE");
189 ✗ omc_write_csv(csvData, buffer);
190
191 /* finish line */
192 ✗ fputc('\n',csvData->handle);
193
194 ✗ return 0;
195 }
196
197 /*! \fn int print_csvLineIterStatsHeader(OMC_WRITE_CSV* csvData)
198 *
199 * This function initializes csv files for analysis propose.
200 *
201 * \param [ref] [data]
202 * \param [ref] [systemData]
203 */
204 ✗ int print_csvLineIterStatsHeader(DATA* data, NONLINEAR_SYSTEM_DATA* systemData, OMC_WRITE_CSV* csvData)
205 {
206 char buffer[1024];
207 int j;
208 ✗ int size = modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).numVar;
209
210 /* number of call */
211 sprintf(buffer,"numberOfCall");
212 ✗ omc_write_csv(csvData, buffer);
213 ✗ fputc(csvData->seperator,csvData->handle);
214
215 /* solving iterations */
216 sprintf(buffer,"iteration");
217 ✗ omc_write_csv(csvData, buffer);
218 ✗ fputc(csvData->seperator,csvData->handle);
219
220 /* variables x */
221 ✗ for(j=0; j<size; ++j) {
222 ✗ sprintf(buffer, "%s", modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[j]);
223 ✗ omc_write_csv(csvData, buffer);
224 ✗ fputc(csvData->seperator,csvData->handle);
225 }
226
227 /* residuals */
228 ✗ for(j=0; j<size; ++j) {
229 ✗ sprintf(buffer, "r%d", j+1);
230 ✗ omc_write_csv(csvData, buffer);
231 ✗ fputc(csvData->seperator,csvData->handle);
232 }
233
234 /* delta x */
235 sprintf(buffer,"delta_x");
236 ✗ omc_write_csv(csvData, buffer);
237 ✗ fputc(csvData->seperator,csvData->handle);
238
239 /* delta x scaled */
240 sprintf(buffer,"delta_x_scaled");
241 ✗ omc_write_csv(csvData, buffer);
242 ✗ fputc(csvData->seperator,csvData->handle);
243
244 /* error in f */
245 sprintf(buffer,"error_f");
246 ✗ omc_write_csv(csvData, buffer);
247 ✗ fputc(csvData->seperator,csvData->handle);
248
249 /* error in f scaled */
250 sprintf(buffer,"error_f_scaled");
251 ✗ omc_write_csv(csvData, buffer);
252 ✗ fputc(csvData->seperator,csvData->handle);
253
254 /* damping lambda */
255 sprintf(buffer,"lambda");
256 ✗ omc_write_csv(csvData, buffer);
257
258 /* finish line */
259 ✗ fputc('\n',csvData->handle);
260
261 ✗ return 0;
262 }
263
264 /*! \fn int print_csvLineIterStatsHeader(OMC_WRITE_CSV* csvData)
265 *
266 * This function initializes csv files for analysis propose.
267 *
268 * \param [ref] [csvData]
269 * \param [in] [size, num, ...]
270 */
271 ✗ int print_csvLineIterStats(void* voidCsvData, int size, int num,
272 int iteration, double* x, double* f, double error_f,
273 double error_fs, double delta_x, double delta_xs,
274 double lambda)
275 {
276 OMC_WRITE_CSV* csvData = voidCsvData;
277 char buffer[1024];
278 int j;
279
280 /* number of call */
281 sprintf(buffer, "%d", num);
282 ✗ omc_write_csv(csvData, buffer);
283 ✗ fputc(csvData->seperator,csvData->handle);
284
285 /* simulation time */
286 sprintf(buffer, "%d", iteration);
287 ✗ omc_write_csv(csvData, buffer);
288 ✗ fputc(csvData->seperator,csvData->handle);
289
290 /* x */
291 ✗ for(j=0; j<size; ++j) {
292 ✗ sprintf(buffer, "%g", x[j]);
293 ✗ omc_write_csv(csvData, buffer);
294 ✗ fputc(csvData->seperator,csvData->handle);
295 }
296
297 /* r */
298 ✗ for(j=0; j<size; ++j) {
299 ✗ sprintf(buffer, "%g", f[j]);
300 ✗ omc_write_csv(csvData, buffer);
301 ✗ fputc(csvData->seperator,csvData->handle);
302 }
303
304 /* error_f */
305 sprintf(buffer, "%g", error_f);
306 ✗ omc_write_csv(csvData, buffer);
307 ✗ fputc(csvData->seperator,csvData->handle);
308
309 /* error_f */
310 sprintf(buffer, "%g", error_fs);
311 ✗ omc_write_csv(csvData, buffer);
312 ✗ fputc(csvData->seperator,csvData->handle);
313
314 /* delta_x */
315 sprintf(buffer, "%g", delta_x);
316 ✗ omc_write_csv(csvData, buffer);
317 ✗ fputc(csvData->seperator,csvData->handle);
318
319 /* delta_xs */
320 sprintf(buffer, "%g", delta_xs);
321 ✗ omc_write_csv(csvData, buffer);
322 ✗ fputc(csvData->seperator,csvData->handle);
323
324 /* lambda */
325 sprintf(buffer, "%g", lambda);
326 ✗ omc_write_csv(csvData, buffer);
327
328 /* finish line */
329 ✗ fputc('\n',csvData->handle);
330
331 ✗ return 0;
332 }
333 #endif
334
335 /**
336 * @brief Allocate and initialize NLS user data.
337 *
338 * NLS user data is passed to 3rdParty non-linear solvers, to be used
339 * in e.g. residual and Jacobian functions provided by C runtime.
340 *
341 * @param data Pointer to data.
342 * @param threadData Pointer to thread data.
343 * @param sysNumber Index of non-linear system.
344 * A non-negative index indicates that nlsData is a pointer to
345 * data->simulationInfo->nonlinearSystemData[sysNumber]
346 * @param nlsData Pointer to non-linear system data corresponding to sysNumber.
347 * @param analyticJacobian Pointer to analytic Jacobian. Can be NULL.
348 * @return NLS_USERDATA* Newly allocated struct with NLS user data.
349 */
350 ✗ NLS_USERDATA* initNlsUserData(DATA* data, threadData_t* threadData, int sysNumber, NONLINEAR_SYSTEM_DATA* nlsData, JACOBIAN* analyticJacobian) {
351 ✗ NLS_USERDATA* userData = (NLS_USERDATA*) malloc(sizeof(NLS_USERDATA));
352 ✗ assertStreamPrint(threadData, userData != NULL, "setNlsUserData failed: userData is NULL");
353
354 ✗ userData->data = data;
355 ✗ userData->threadData = threadData;
356 ✗ userData->sysNumber = sysNumber;
357 ✗ userData->nlsData = nlsData;
358 ✗ userData->analyticJacobian = analyticJacobian;
359 ✗ userData->solverData = NULL;
360
361 ✗ return userData;
362 }
363
364 /**
365 * @brief Free NLS user data struct.
366 *
367 * Only frees memory of NLS data struct, not of its members.
368 *
369 * @param userData Pointer to NLS user data.
370 */
371 ✗ void freeNlsUserData(NLS_USERDATA* userData) {
372 ✗ free(userData);
373 ✗ }
374
375 /* The adaptive homotopy methods solve for lambda alongside the unknowns, so
376 * they take one unknown off the system and add it back as a column. */
377 static inline modelica_boolean adaptiveHomotopy(DATA *data, NONLINEAR_SYSTEM_DATA *nonlinsys)
378 {
379 ✗ return nonlinsys->homotopySupport
380 ✗ && (data->callback->homotopyMethod == GLOBAL_ADAPTIVE_HOMOTOPY
381 ✗ || data->callback->homotopyMethod == LOCAL_ADAPTIVE_HOMOTOPY);
382 }
383
384 /**
385 * @brief Initialize internal structure of non-linear system.
386 *
387 * @param data Runtime data struct.
388 * @param threadData Thread data for error handling.
389 * @param nonlinsys Pointer to non-linear system.
390 * @param sysNum Number of non-linear system.
391 */
392 ✗ void initializeNonlinearSystemData(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA *nonlinsys, int sysNum) {
393 modelica_integer size;
394 struct dataSolver *solverData;
395 struct dataMixedSolver *mixedSolverData;
396 JACOBIAN* jacobian;
397
398 ✗ size = nonlinsys->size;
399 ✗ nonlinsys->numberOfFEval = 0;
400 ✗ nonlinsys->numberOfIterations = 0;
401
402 /* check if residual function pointer are valid */
403 ✗ assertStreamPrint(threadData, (nonlinsys->residualFunc != NULL) || (nonlinsys->strictTearingFunctionCall != NULL), "residual function pointer is invalid");
404
405 /* check if analytical jacobian is created */
406 ✗ if(nonlinsys->jacobianIndex != -1)
407 {
408 ✗ jacobian = &(data->simulationInfo->analyticJacobians[nonlinsys->jacobianIndex]);
409 ✗ assertStreamPrint(threadData, 0 != nonlinsys->analyticalJacobianColumn, "jacobian function pointer is invalid" );
410 ✗ if(nonlinsys->initialAnalyticalJacobian(data, threadData, jacobian))
411 {
412 ✗ nonlinsys->jacobianIndex = -1;
413 /* simulation_data.h documents analyticalJacobianColumn==NULL as THE signal
414 * that no analytic Jacobian is available; clear it here too, not just
415 * jacobianIndex, so every consumer agrees (e.g. kinsolSolver.c's
416 * initKinsolMemory only checks analyticalJacobianColumn, not jacobianIndex,
417 * to decide whether to wire up the symbolic-Jacobian KINSOL callback -- left
418 * dangling non-NULL here, it would still select that callback, which then
419 * fails an assertion pulling a JACOBIAN* through the now -1 jacobianIndex). */
420 ✗ nonlinsys->analyticalJacobianColumn = NULL;
421 jacobian = NULL;
422 }
423 /* evalJacobian() fills sizeRows*sizeCols entries, the solvers hand it a
424 * rows x (rows+1) buffer. A wider Jacobian would run over it, and its
425 * columns would not be the iteration variables either. */
426 ✗ else if(jacobian->sizeRows != (adaptiveHomotopy(data, nonlinsys) ? size - 1 : size) || jacobian->sizeCols != size)
427 {
428 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Analytic Jacobian of non-linear system %d is %ux%u, but the system has " OMC_INT_FORMAT " iteration variables. "
429 "This indicates that something went wrong during Jacobian generation. "
430 "Using a numeric Jacobian instead.",
431 sysNum, jacobian->sizeRows, jacobian->sizeCols, size);
432 ✗ nonlinsys->jacobianIndex = -1;
433 ✗ nonlinsys->analyticalJacobianColumn = NULL;
434 jacobian = NULL;
435 }
436 } else {
437 jacobian = NULL;
438 }
439
440 /* allocate system data; zeroed because the "last solving is too long ago" branch
441 * of solve_nonlinear_system leaves nlsxExtrapolation -- which solveHomotopy and
442 * solveHybrd start a non-discrete call from -- unwritten. */
443 ✗ nonlinsys->nlsx = (double*) calloc(size, sizeof(double));
444 ✗ nonlinsys->nlsxExtrapolation = (double*) calloc(size, sizeof(double));
445 ✗ nonlinsys->nlsxOld = (double*) calloc(size, sizeof(double));
446 ✗ nonlinsys->resValues = (double*) calloc(size, sizeof(double));
447
448 /* allocate value list*/
449 ✗ nonlinsys->oldValueList = allocValueList(1, nonlinsys->size);
450
451 ✗ nonlinsys->lastTimeSolved = 0.0;
452
453 /* Allocate nomianl, min and max */
454 ✗ nonlinsys->nominal = (double*) malloc(size*sizeof(double));
455 ✗ nonlinsys->min = (double*) malloc(size*sizeof(double));
456 ✗ nonlinsys->max = (double*) malloc(size*sizeof(double));
457 /* Init sparsitiy pattern */
458 ✗ nonlinsys->initializeStaticNLSData(data, threadData, nonlinsys, 1 /* true */, 1 /* true */);
459
460 ✗ if(nonlinsys->sparsePattern) {
461 /* only test for singularity if sparsity pattern is supposed to be there */
462 ✗ modelica_boolean useSparsityPattern = sparsitySanityCheck(nonlinsys->sparsePattern, nonlinsys->size, OMC_LOG_NLS);
463 ✗ if (!useSparsityPattern) {
464 // free sparsity pattern and don't use scaling
465 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Sparsity pattern for non-linear system %d is not regular. "
466 "This indicates that something went wrong during sparsity pattern generation. "
467 "Removing sparsity pattern and disabling NLS scaling.", sysNum);
468 /* DEBUG */
469 //printSparseStructure(nonlinsys->sparsePattern, nonlinsys->size, nonlinsys->size, OMC_LOG_NLS, "NLS sparse pattern");
470 ✗ freeSparsePattern(nonlinsys->sparsePattern);
471 ✗ nonlinsys->sparsePattern = NULL;
472 ✗ omc_flag[FLAG_NO_SCALING] = TRUE;
473 }
474 }
475
476 #if !defined(OMC_MINIMAL_RUNTIME)
477 /* csv data call stats*/
478 ✗ if (data->simulationInfo->nlsCsvInfomation)
479 {
480 ✗ if (initializeNLScsvData(data, nonlinsys))
481 {
482 ✗ throwStreamPrint(threadData, "csvData initialization failed");
483 }
484 else
485 {
486 ✗ print_csvLineCallStatsHeader(((struct csvStats*) nonlinsys->csvData)->callStats);
487 ✗ print_csvLineIterStatsHeader(data, nonlinsys, ((struct csvStats*) nonlinsys->csvData)->iterStats);
488 }
489 }
490 #endif
491
492 ✗ nonlinsys->nlsMethod = data->simulationInfo->nlsMethod;
493 ✗ nonlinsys->nlsLinearSolver = data->simulationInfo->nlsLinearSolver;
494 #if !defined(OMC_MINIMAL_RUNTIME)
495 ✗ if (nonlinsys->matrixFormat == OMC_MATRIX_SPARSE && nonlinsys->sparsePattern
496 ✗ && !(data->simulationInfo->nlsMethod == NLS_KINSOL || data->simulationInfo->nlsMethod == NLS_KINSOL_B))
497 {
498 ✗ nonlinsys->nlsMethod = NLS_KINSOL;
499 ✗ nonlinsys->nlsLinearSolver = NLS_LS_KLU;
500 }
501 #endif
502
503 /* Set NLS user data */
504 ✗ NLS_USERDATA* nlsUserData = initNlsUserData(data, threadData, sysNum, nonlinsys, jacobian);
505
506 // FIXME: add generation of scalar system (total system size = 1) to codegen
507 #if !defined(OMC_MINIMAL_RUNTIME)
508 /* check for trivial sparsity pattern (is not always generated) */
509 ✗ if (nonlinsys->size == 1 && !nonlinsys->sparsePattern && (nonlinsys->nlsMethod == NLS_KINSOL || nonlinsys->nlsMethod == NLS_KINSOL_B)) {
510 ✗ nonlinsys->nlsMethod = NLS_HYBRID;
511 ✗ nonlinsys->nlsLinearSolver = NLS_LS_DEFAULT;
512 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Sparsity pattern for non-linear system %d with size 1x1 (nnz = 1) does not exist. "
513 "Can not use the set sparse NLS solver - changing the method to Hybrid.", sysNum);
514 }
515 #endif
516
517 /* allocate stuff depending on the chosen method */
518 ✗ switch(nonlinsys->nlsMethod)
519 {
520 #if !defined(OMC_MINIMAL_RUNTIME)
521 ✗ case NLS_HYBRID:
522 ✗ solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
523 if (adaptiveHomotopy(data, nonlinsys)) {
524 ✗ solverData->ordinaryData = allocateHybrdData(size-1, nlsUserData);
525 ✗ nlsUserData = initNlsUserData(data, threadData, sysNum, nonlinsys, jacobian); /* Seperate userData for homotopy solver */
526 ✗ solverData->initHomotopyData = (void*) allocateHomotopyData(size-1, nlsUserData);
527 } else {
528 ✗ solverData->ordinaryData = allocateHybrdData(size, nlsUserData);
529 }
530 ✗ nonlinsys->solverData = (void*) solverData;
531 ✗ break;
532 ✗ case NLS_KINSOL:
533 ✗ solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
534 if (adaptiveHomotopy(data, nonlinsys)) {
535 ✗ solverData->initHomotopyData = (void*) allocateHomotopyData(size-1, nlsUserData);
536 } else {
537 ✗ nonlinsys->solverData = (void*) nlsKinsolAllocate(size, nlsUserData, TRUE, !!nonlinsys->sparsePattern);
538 ✗ solverData->ordinaryData = nonlinsys->solverData;
539 }
540 ✗ nonlinsys->solverData = (void*) solverData;
541 ✗ break;
542 ✗ case NLS_KINSOL_B:
543 ✗ solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
544 if (adaptiveHomotopy(data, nonlinsys)) {
545 ✗ solverData->initHomotopyData = (void*) allocateHomotopyData(size-1, nlsUserData);
546 } else {
547 ✗ nonlinsys->solverData = (void*) B_nlsKinsolAllocate(size, nlsUserData, TRUE, !!nonlinsys->sparsePattern);
548 ✗ solverData->ordinaryData = nonlinsys->solverData;
549 }
550 ✗ nonlinsys->solverData = (void*) solverData;
551 ✗ break;
552 ✗ case NLS_NEWTON:
553 ✗ solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
554 if (adaptiveHomotopy(data, nonlinsys)) {
555 ✗ solverData->ordinaryData = (void*) allocateNewtonData(size-1, nlsUserData);
556 ✗ nlsUserData = initNlsUserData(data, threadData, sysNum, nonlinsys, jacobian); /* Seperate userData for homotopy solver */
557 ✗ solverData->initHomotopyData = (void*) allocateHomotopyData(size-1, nlsUserData);
558 } else {
559 ✗ solverData->ordinaryData = (void*) allocateNewtonData(size, nlsUserData);
560 }
561 ✗ nonlinsys->solverData = (void*) solverData;
562 ✗ break;
563 ✗ case NLS_MIXED:
564 ✗ mixedSolverData = (struct dataMixedSolver*) malloc(sizeof(struct dataMixedSolver));
565 if (adaptiveHomotopy(data, nonlinsys)) {
566 ✗ mixedSolverData->newtonHomotopyData = (void*) allocateHomotopyData(size-1, nlsUserData);
567 ✗ nlsUserData = initNlsUserData(data, threadData, sysNum, nonlinsys, jacobian); /* Seperate userData for hybrid solver */
568 ✗ mixedSolverData->hybridData = (void*) allocateHybrdData(size-1, nlsUserData);
569 } else {
570 ✗ mixedSolverData->newtonHomotopyData = (void*) allocateHomotopyData(size, nlsUserData);
571 ✗ nlsUserData = initNlsUserData(data, threadData, sysNum, nonlinsys, jacobian); /* Seperate userData for hybrid solver */
572 ✗ mixedSolverData->hybridData = (void*) allocateHybrdData(size, nlsUserData);
573 }
574 ✗ nonlinsys->solverData = (void*) mixedSolverData;
575 ✗ break;
576 #endif
577 case NLS_HOMOTOPY:
578 if (adaptiveHomotopy(data, nonlinsys)) {
579 ✗ nonlinsys->solverData = (void*) allocateHomotopyData(size-1, nlsUserData);
580 } else {
581 ✗ nonlinsys->solverData = (void*) allocateHomotopyData(size, nlsUserData);
582 }
583 break;
584 ✗ default:
585 ✗ throwStreamPrint(threadData, "unrecognized nonlinear solver");
586 }
587
588 ✗ return;
589 }
590
591 /**
592 * @brief Initialize all non-linear systems.
593 *
594 * Loop over all non-linear systems and call initializeNonlinearSystemData.
595 * Memory for data->simulationInfo->nonlinearSystemData is already allocated.
596 *
597 * @param data Runtime data struct.
598 * @param threadData Thread data for error handling.
599 * @return int Return 0 on success.
600 */
601 1 int initializeNonlinearSystems(DATA *data, threadData_t *threadData)
602 {
603 int i;
604 1 NONLINEAR_SYSTEM_DATA *nonlinsys = data->simulationInfo->nonlinearSystemData;
605
606 1 infoStreamPrint(OMC_LOG_NLS, 1, "initialize non-linear system solvers");
607 1 infoStreamPrint(OMC_LOG_NLS, 0, "%ld non-linear systems", data->modelData->nNonLinearSystems);
608
609 /* set the default nls linear solver depending on the default nls method */
610
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (data->simulationInfo->nlsLinearSolver == NLS_LS_DEFAULT) {
611 #if !defined(OMC_MINIMAL_RUNTIME)
612 /* kinsol works best with KLU,
613 they are both sparse so it makes sense to use them together */
614
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (data->simulationInfo->nlsMethod == NLS_KINSOL || data->simulationInfo->nlsMethod == NLS_KINSOL_B) {
615 ✗ data->simulationInfo->nlsLinearSolver = NLS_LS_KLU;
616 } else {
617 1 data->simulationInfo->nlsLinearSolver = NLS_LS_LAPACK;
618 }
619 #else
620 ✗ data->simulationInfo->nlsLinearSolver = NLS_LS_LAPACK;
621 #endif
622 }
623
624
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for (i=0; i<data->modelData->nNonLinearSystems; ++i) {
625 ✗ initializeNonlinearSystemData(data, threadData, &nonlinsys[i], i);
626 }
627
628 1 messageClose(OMC_LOG_NLS);
629
630 1 return 0;
631 }
632
633 /**
634 * @brief Initialize min, max, nominal for non-linear systems.
635 *
636 * This function allocates memory for sparsity pattern and
637 * initialized nominal, min, max and spsarsity pattern.
638 *
639 * @param data Pointer to data.
640 * @param threadData Thread data for error handling.
641 * @return int Return 0.
642 */
643 1 int updateStaticDataOfNonlinearSystems(DATA *data, threadData_t *threadData)
644 {
645 int i;
646 int size;
647 1 NONLINEAR_SYSTEM_DATA *nonlinsys = data->simulationInfo->nonlinearSystemData;
648
649 1 infoStreamPrint(OMC_LOG_NLS, 1, "update static data of non-linear system solvers");
650
651
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(i=0; i<data->modelData->nNonLinearSystems; ++i)
652 {
653 ✗ nonlinsys[i].initializeStaticNLSData(data, threadData, &nonlinsys[i], FALSE, FALSE);
654 }
655
656 1 messageClose(OMC_LOG_NLS);
657
658 1 return 0;
659 }
660
661 /**
662 * @brief Free non-linear system data.
663 *
664 * Freeing non-linear solver data as well.
665 *
666 * @param data Pointer to data struct.
667 * @param threadData Pointer to thread data. Used for error handling.
668 * @param nonlinsys Pointer to non-linear system data.
669 */
670 ✗ void freeNonlinearSyst(DATA* data, threadData_t* threadData, NONLINEAR_SYSTEM_DATA* nonlinsys)
671 {
672 struct csvStats* stats;
673 NLS_USERDATA* userData;
674
675 ✗ free(nonlinsys->nlsx);
676 ✗ free(nonlinsys->nlsxExtrapolation);
677 ✗ free(nonlinsys->nlsxOld);
678 ✗ free(nonlinsys->resValues);
679 ✗ free(nonlinsys->nominal);
680 ✗ free(nonlinsys->min);
681 ✗ free(nonlinsys->max);
682 ✗ free(nonlinsys->eqn_simcode_indices);
683 ✗ nonlinsys->eqn_simcode_indices = NULL;
684 ✗ nonlinsys->freeStaticNLSData(data, threadData, nonlinsys);
685
686 ✗ freeValueList(nonlinsys->oldValueList, 1);
687 ✗ freeNonlinearPattern(nonlinsys->nonlinearPattern);
688 ✗ nonlinsys->nonlinearPattern = NULL;
689
690 /* Free CSV data */
691 #if !defined(OMC_MINIMAL_RUNTIME)
692 ✗ if (data->simulationInfo->nlsCsvInfomation)
693 {
694 ✗ stats = nonlinsys->csvData;
695 ✗ omc_write_csv_free(stats->callStats);
696 ✗ omc_write_csv_free(stats->iterStats);
697 ✗ free(nonlinsys->csvData);
698 // TODO: Make a function freeNLScsvData
699 }
700 #endif
701
702 /* free solver data */
703 // TODO: Make this simpler. Don't cast nonlinsys->solverData back and forth.
704 // Always have a structure accepting two solver methods (primary and backup) and make the missing one NULL.
705 // Also just set nonlinsys->solverData->ordinaryData to NULL, if no hybrid solver is available
706 // and make freeHybrdData(NULL) work.
707 ✗ switch(nonlinsys->nlsMethod)
708 {
709 #if !defined(OMC_MINIMAL_RUNTIME)
710 ✗ case NLS_HYBRID:
711 ✗ freeHybrdData(((struct dataSolver*) nonlinsys->solverData)->ordinaryData);
712 if (adaptiveHomotopy(data, nonlinsys)) {
713 ✗ freeHomotopyData(((struct dataSolver*) nonlinsys->solverData)->initHomotopyData);
714 }
715 ✗ free(nonlinsys->solverData);
716 ✗ break;
717 case NLS_KINSOL:
718 if (adaptiveHomotopy(data, nonlinsys)) {
719 ✗ freeHomotopyData(((struct dataSolver*) nonlinsys->solverData)->initHomotopyData);
720 } else {
721 ✗ nlsKinsolFree(((struct dataSolver*) nonlinsys->solverData)->ordinaryData);
722 }
723 ✗ free(nonlinsys->solverData);
724 ✗ break;
725 case NLS_KINSOL_B:
726 if (adaptiveHomotopy(data, nonlinsys)) {
727 ✗ freeHomotopyData(((struct dataSolver*) nonlinsys->solverData)->initHomotopyData);
728 } else {
729 ✗ B_nlsKinsolFree(((struct dataSolver*) nonlinsys->solverData)->ordinaryData);
730 }
731 ✗ free(nonlinsys->solverData);
732 ✗ break;
733 ✗ case NLS_NEWTON:
734 ✗ freeNewtonData(((struct dataSolver*) nonlinsys->solverData)->ordinaryData);
735 if (adaptiveHomotopy(data, nonlinsys)) {
736 ✗ freeHomotopyData(((struct dataSolver*) nonlinsys->solverData)->initHomotopyData);
737 }
738 ✗ free(nonlinsys->solverData);
739 ✗ break;
740 #endif
741 ✗ case NLS_HOMOTOPY:
742 ✗ freeHomotopyData(nonlinsys->solverData);
743 ✗ break;
744 #if !defined(OMC_MINIMAL_RUNTIME)
745 ✗ case NLS_MIXED:
746 ✗ freeHomotopyData(((struct dataMixedSolver*) nonlinsys->solverData)->newtonHomotopyData);
747 ✗ freeHybrdData(((struct dataMixedSolver*) nonlinsys->solverData)->hybridData);
748 ✗ free(nonlinsys->solverData);
749 ✗ break;
750 #endif
751 ✗ default:
752 ✗ throwStreamPrint(threadData, "freeNonlinearSyst: Unrecognized non-linear solver method");
753 }
754
755 ✗ return;
756 }
757
758 /**
759 * @brief Free memory for all non-linear systems in simulationInfo.
760 *
761 * Free all non-linear systems from array data->simulationInfo->nonlinearSystemData.
762 *
763 * @param data Pointer to data struct.
764 * @param threadData Pointer to thread data. Used for error handling.
765 */
766 1 void freeNonlinearSystems(DATA *data, threadData_t *threadData)
767 {
768 int i;
769 1 NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
770
771 1 infoStreamPrint(OMC_LOG_NLS, 1, "free non-linear system solvers");
772
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(i=0; i<data->modelData->nNonLinearSystems; ++i)
773 {
774 ✗ freeNonlinearSyst(data, threadData, &nonlinsys[i]);
775 }
776
777 1 messageClose(OMC_LOG_NLS);
778
779 1 return;
780 }
781
782 /**
783 * @brief Print non-linear system statistics to stream.
784 *
785 * @param nonlinsys Non-linear system data to print.
786 * @param stream Log stream to use for logging.
787 */
788 ✗ void printNonLinearSystemSolvingStatistics(NONLINEAR_SYSTEM_DATA* nonlinsys, enum OMC_LOG_STREAM stream)
789 {
790 ✗ if (!OMC_ACTIVE_STREAM(stream)) return;
791 ✗ infoStreamPrint(stream, 1, "Non-linear system %d of size %d solver statistics:", (int)nonlinsys->equationIndex, (int)nonlinsys->size);
792 ✗ infoStreamPrint(stream, 0, " number of calls : %ld", nonlinsys->numberOfCall);
793 ✗ infoStreamPrint(stream, 0, " number of iterations : %ld", nonlinsys->numberOfIterations);
794 ✗ infoStreamPrint(stream, 0, " number of function evaluations : %ld", nonlinsys->numberOfFEval);
795 ✗ infoStreamPrint(stream, 0, " number of jacobian evaluations : %ld", nonlinsys->numberOfJEval);
796 ✗ infoStreamPrint(stream, 0, " time of jacobian evaluations : %f", nonlinsys->jacobianTime);
797 ✗ infoStreamPrint(stream, 0, " average time per call : %f", nonlinsys->totalTime/nonlinsys->numberOfCall);
798 ✗ infoStreamPrint(stream, 0, " total time : %f", nonlinsys->totalTime);
799 ✗ messageClose(stream);
800 }
801
802 /*! \fn printNonLinearInitialInfo
803 *
804 * This function prints information of an non-linear systems before an solving step.
805 *
806 * \param [in] [logName] log level in general OMC_LOG_NLS
807 * [ref] [data]
808 * [ref] [nonlinsys] index of corresponding non-linear system
809 */
810 ✗ void printNonLinearInitialInfo(int logName, DATA* data, NONLINEAR_SYSTEM_DATA *nonlinsys)
811 {
812 long i;
813
814 ✗ if (!OMC_ACTIVE_STREAM(logName)) return;
815 ✗ infoStreamPrint(logName, 1, "initial variable values:");
816
817 ✗ for(i=0; i<nonlinsys->size; i++)
818 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g", i+1,
819 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,nonlinsys->equationIndex).vars[i],
820 ✗ nonlinsys->nlsx[i], nonlinsys->nominal[i]);
821 ✗ messageClose(logName);
822 }
823
824 /*! \fn printNonLinearFinishInfo
825 *
826 * This function prints information of an non-linear systems after a solving step.
827 *
828 * \param [in] [logName] log level in general OMC_LOG_NLS
829 * [ref] [data]
830 * [ref] [nonlinsys] index of corresponding non-linear system
831 */
832 ✗ void printNonLinearFinishInfo(int logName, DATA* data, NONLINEAR_SYSTEM_DATA *nonlinsys)
833 {
834 long i;
835
836 ✗ if (!OMC_ACTIVE_STREAM(logName)) return;
837
838 ✗ switch (nonlinsys->solved)
839 {
840 ✗ case NLS_SOLVED:
841 ✗ infoStreamPrint(logName, 1, "Solution status: SOLVED");
842 ✗ break;
843 ✗ case NLS_SOLVED_LESS_ACCURACY:
844 ✗ infoStreamPrint(logName, 1, "Solution status: SOLVED with less accuracy");
845 ✗ break;
846 ✗ case NLS_FAILED:
847 ✗ infoStreamPrint(logName, 1, "Solution status: FAILED");
848 ✗ break;
849 ✗ default:
850 ✗ throwStreamPrint(NULL, "Unhandled case in printNonLinearFinishInfo");
851 break;
852 }
853 ✗ infoStreamPrint(logName, 0, " number of iterations : %ld", nonlinsys->numberOfIterations);
854 ✗ infoStreamPrint(logName, 0, " number of function evaluations : %ld", nonlinsys->numberOfFEval);
855 ✗ infoStreamPrint(logName, 0, " number of jacobian evaluations : %ld", nonlinsys->numberOfJEval);
856 ✗ infoStreamPrint(logName, 0, "solution values:");
857 ✗ for(i=0; i<nonlinsys->size; i++)
858 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g", i+1,
859 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,nonlinsys->equationIndex).vars[i],
860 ✗ nonlinsys->nlsx[i]);
861
862 ✗ messageClose(logName);
863 }
864
865 /*! \fn getInitialGuess
866 *
867 * This function writes initial guess to nonlinsys->nlsx and nonlinsys->nlsOld.
868 *
869 * \param [in] [nonlinsys]
870 * \param [in] [time] time for extrapolation
871 *
872 */
873 ✗ int getInitialGuess(NONLINEAR_SYSTEM_DATA *nonlinsys, double time)
874 {
875 /* value extrapolation */
876 ✗ printValuesListTimes(nonlinsys->oldValueList->valueList);
877 /* if list is empty use current start values */
878 ✗ if (listLen(nonlinsys->oldValueList->valueList)==0)
879 {
880 /* use old value if no values are stored in the list; the extrapolation too,
881 which `solveHomotopy` starts a non-discrete call from and which is a zeroed
882 allocation until then */
883 ✗ memcpy(nonlinsys->nlsx, nonlinsys->nlsxOld, nonlinsys->size*(sizeof(double)));
884 ✗ memcpy(nonlinsys->nlsxExtrapolation, nonlinsys->nlsxOld, nonlinsys->size*(sizeof(double)));
885 }
886 else
887 {
888 /* get extrapolated values */
889 ✗ getValues(nonlinsys->oldValueList->valueList, time, nonlinsys->nlsxExtrapolation, nonlinsys->nlsxOld);
890 ✗ memcpy(nonlinsys->nlsx, nonlinsys->nlsxOld, nonlinsys->size*(sizeof(double)));
891 }
892
893 ✗ return 0;
894 }
895
896 /*! \fn updateInitialGuessDB
897 *
898 * This function writes new values to solution list.
899 *
900 * \param [in] [nonlinsys]
901 * \param [in] [time] time for extrapolation
902 * \param [in] [context] current context of evaluation
903 *
904 */
905 ✗ int updateInitialGuessDB(NONLINEAR_SYSTEM_DATA *nonlinsys, double time, EVAL_CONTEXT context)
906 {
907 /* Variables */
908 VALUE* tmpNode;
909
910 /* write solution to oldValue list for extrapolation */
911 ✗ if (nonlinsys->solved == NLS_SOLVED)
912 {
913 /* do not use solution of jacobian for next extrapolation */
914 ✗ if (context == CONTEXT_ODE || context == CONTEXT_ALGEBRAIC || context == CONTEXT_EVENTS)
915 {
916 ✗ tmpNode = createValueElement(nonlinsys->size, time, nonlinsys->nlsx);
917 ✗ addListElement(nonlinsys->oldValueList->valueList,
918 tmpNode);
919 ✗ freeValue(tmpNode);
920 }
921 }
922 ✗ else if (nonlinsys->solved == NLS_SOLVED_LESS_ACCURACY)
923 {
924 ✗ if (listLen((nonlinsys->oldValueList)->valueList)>0)
925 {
926 ✗ cleanValueList(nonlinsys->oldValueList->valueList, NULL);
927 }
928 /* do not use solution of jacobian for next extrapolation */
929 ✗ if (context == CONTEXT_ODE || context == CONTEXT_ALGEBRAIC || context == CONTEXT_EVENTS)
930 {
931 ✗ tmpNode = createValueElement(nonlinsys->size, time, nonlinsys->nlsx);
932 ✗ addListElement(nonlinsys->oldValueList->valueList,
933 tmpNode);
934 ✗ freeValue(tmpNode);
935 }
936 }
937 ✗ return 0;
938 }
939
940 /*! \fn updateInnerEquation
941 *
942 * This function updates inner equation with the current x.
943 *
944 * \param [ref] [data]
945 * [in] [sysNumber] index of corresponding non-linear system
946 */
947 ✗ int updateInnerEquation(RESIDUAL_USERDATA* resUserData, int sysNumber, int discrete)
948 {
949 ✗ DATA *data = resUserData->data;
950 ✗ threadData_t *threadData = resUserData->threadData;
951
952 ✗ NONLINEAR_SYSTEM_DATA* nonlinsys = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
953 int success = 0;
954 int constraintViolated = 0;
955
956 /* solve non continuous at discrete points*/
957 ✗ if(discrete)
958 {
959 ✗ data->simulationInfo->solveContinuous = 0;
960 }
961
962 /* try */
963 #ifndef OMC_EMCC
964 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
965 #endif
966
967 /* call residual function */
968 ✗ if (nonlinsys->strictTearingFunctionCall != NULL)
969 ✗ constraintViolated = nonlinsys->residualFuncConstraints(resUserData, nonlinsys->nlsx, nonlinsys->resValues, (int*)&nonlinsys->size);
970 else
971 ✗ nonlinsys->residualFunc(resUserData, nonlinsys->nlsx, nonlinsys->resValues, (int*)&nonlinsys->size);
972
973 /* replace extrapolated values by current x for discrete step */
974 ✗ memcpy(nonlinsys->nlsxExtrapolation, nonlinsys->nlsx, nonlinsys->size*(sizeof(double)));
975
976 ✗ if (!constraintViolated && !OMC_ERROR_RAISED())
977 success = 1;
978 /*catch */
979 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
980 #ifndef OMC_EMCC
981 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
982 #endif
983
984 ✗ if (!success && !constraintViolated)
985 {
986 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Non-Linear Solver try to handle a problem with a called assert.");
987 }
988
989 ✗ if(discrete)
990 {
991 ✗ data->simulationInfo->solveContinuous = 1;
992 }
993
994 ✗ return success;
995 }
996
997 /**
998 * @brief Solve given non-linear system.
999 *
1000 * @param data Runtime data struct.
1001 * @param threadData Thread data for error handling.
1002 * @param nonlinsys Pointer to non-linear system.
1003 * @return NLS_SOLVER_STATUS Return NLS_SOLVED on success,
1004 * NLS_SOLVED_LESS_ACCURACY if a less accurate solution was found and
1005 * NLS_FAILED otherwise.
1006 */
1007 ✗ NLS_SOLVER_STATUS solveNLS(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nonlinsys)
1008 {
1009 ✗ NLS_SOLVER_STATUS solver_status = NLS_FAILED;
1010 ✗ int casualTearingSet = nonlinsys->strictTearingFunctionCall != NULL;
1011 struct dataSolver *solverData;
1012 struct dataMixedSolver *mixedSolverData;
1013
1014 /* use the selected solver for solving nonlinear system */
1015 ✗ switch(nonlinsys->nlsMethod)
1016 {
1017 #if !defined(OMC_MINIMAL_RUNTIME)
1018 ✗ case NLS_HYBRID:
1019 ✗ solverData = nonlinsys->solverData;
1020 ✗ nonlinsys->solverData = solverData->ordinaryData;
1021 /* try */
1022 #ifndef OMC_EMCC
1023 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1024 #endif
1025 ✗ solver_status = solveHybrd(data, threadData, nonlinsys);
1026 /*catch */
1027 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); solver_status = NLS_FAILED; }
1028 #ifndef OMC_EMCC
1029 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1030 #endif
1031 ✗ nonlinsys->solverData = solverData;
1032 ✗ break;
1033 ✗ case NLS_KINSOL:
1034 ✗ solverData = nonlinsys->solverData;
1035 ✗ nonlinsys->solverData = solverData->ordinaryData;
1036 /* try */
1037 #ifndef OMC_EMCC
1038 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1039 #endif
1040 ✗ solver_status = nlsKinsolSolve(data, threadData, nonlinsys);
1041 /*catch */
1042 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); solver_status = NLS_FAILED; }
1043 #ifndef OMC_EMCC
1044 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1045 #endif
1046 ✗ nonlinsys->solverData = solverData;
1047 ✗ break;
1048 ✗ case NLS_KINSOL_B:
1049 ✗ solverData = nonlinsys->solverData;
1050 ✗ nonlinsys->solverData = solverData->ordinaryData;
1051 /* try */
1052 #ifndef OMC_EMCC
1053 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1054 #endif
1055 ✗ solver_status = B_nlsKinsolSolve(data, threadData, nonlinsys);
1056 /*catch */
1057 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); solver_status = NLS_FAILED; }
1058 #ifndef OMC_EMCC
1059 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1060 #endif
1061 ✗ nonlinsys->solverData = solverData;
1062 ✗ break;
1063 ✗ case NLS_NEWTON:
1064 ✗ solverData = nonlinsys->solverData;
1065 ✗ nonlinsys->solverData = solverData->ordinaryData;
1066 /* try */
1067 #ifndef OMC_EMCC
1068 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1069 #endif
1070 ✗ solver_status = solveNewton(data, threadData, nonlinsys);
1071 /*catch */
1072 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); solver_status = NLS_FAILED; }
1073 #ifndef OMC_EMCC
1074 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1075 #endif
1076 /* check if solution process was successful, if not use alternative tearing set if available (dynamic tearing)*/
1077 ✗ if (solver_status != NLS_SOLVED && casualTearingSet){
1078 debugString(OMC_LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
1079 ✗ if(nonlinsys->strictTearingFunctionCall(data, threadData)) {
1080 solver_status = NLS_SOLVED;
1081 } else {
1082 solver_status = NLS_FAILED;
1083 }
1084 }
1085 ✗ nonlinsys->solverData = solverData;
1086 ✗ break;
1087 #endif
1088 ✗ case NLS_HOMOTOPY:
1089 ✗ solver_status = solveHomotopy(data, threadData, nonlinsys);
1090 break;
1091 #if !defined(OMC_MINIMAL_RUNTIME)
1092 ✗ case NLS_MIXED:
1093 ✗ mixedSolverData = nonlinsys->solverData;
1094 ✗ nonlinsys->solverData = mixedSolverData->newtonHomotopyData;
1095
1096 /* try */
1097 #ifndef OMC_EMCC
1098 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1099 #endif
1100 ✗ solver_status = solveHomotopy(data, threadData, nonlinsys);
1101
1102 /* check if solution process was successful, if not use alternative tearing set if available (dynamic tearing)*/
1103 ✗ if (solver_status != NLS_SOLVED && casualTearingSet){
1104 debugString(OMC_LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
1105 ✗ if(nonlinsys->strictTearingFunctionCall(data, threadData)) {
1106 solver_status = NLS_SOLVED;
1107 } else {
1108 solver_status = NLS_FAILED;
1109 }
1110 }
1111
1112 ✗ if (solver_status != NLS_SOLVED ) {
1113 ✗ nonlinsys->solverData = mixedSolverData->hybridData;
1114 ✗ solver_status = solveHybrd(data, threadData, nonlinsys);
1115 }
1116
1117 /* update iteration variables of nonlinsys->nlsx */
1118 ✗ if (solver_status == NLS_SOLVED){
1119 ✗ nonlinsys->getIterationVars(data, nonlinsys->nlsx);
1120 }
1121
1122 /*catch */
1123 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); solver_status = NLS_FAILED; }
1124 #ifndef OMC_EMCC
1125 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1126 #endif
1127 ✗ nonlinsys->solverData = mixedSolverData;
1128 ✗ break;
1129 #endif
1130 ✗ default:
1131 ✗ throwStreamPrint(threadData, "unrecognized nonlinear solver");
1132 }
1133
1134 ✗ return solver_status;
1135 }
1136
1137 /**
1138 * @brief Solve initial non-linear system with homotopy solver
1139 *
1140 * @param data Runtime data struct.
1141 * @param threadData Thread data for error handling.
1142 * @param nonlinsys Pointer to non-linear system data.
1143 * @return NLS_SOLVER_STATUS Return solver status of homotopy solver.
1144 */
1145 ✗ NLS_SOLVER_STATUS solveWithInitHomotopy(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nonlinsys)
1146 {
1147 NLS_SOLVER_STATUS success = NLS_FAILED;
1148 struct dataSolver *solverData;
1149 struct dataMixedSolver *mixedSolverData;
1150
1151 /* use the homotopy solver for solving the initial system */
1152 ✗ switch(nonlinsys->nlsMethod)
1153 {
1154 #if !defined(OMC_MINIMAL_RUNTIME)
1155 ✗ case NLS_HYBRID:
1156 case NLS_KINSOL:
1157 case NLS_KINSOL_B:
1158 case NLS_NEWTON:
1159 ✗ solverData = nonlinsys->solverData;
1160 ✗ nonlinsys->solverData = solverData->initHomotopyData;
1161 ✗ success = solveHomotopy(data, threadData, nonlinsys);
1162 ✗ nonlinsys->solverData = solverData;
1163 ✗ break;
1164 #endif
1165 ✗ case NLS_HOMOTOPY:
1166 ✗ success = solveHomotopy(data, threadData, nonlinsys);
1167 ✗ break;
1168 #if !defined(OMC_MINIMAL_RUNTIME)
1169 ✗ case NLS_MIXED:
1170 ✗ mixedSolverData = nonlinsys->solverData;
1171 ✗ nonlinsys->solverData = mixedSolverData->newtonHomotopyData;
1172 ✗ success = solveHomotopy(data, threadData, nonlinsys);
1173 ✗ nonlinsys->solverData = mixedSolverData;
1174 ✗ break;
1175 #endif
1176 ✗ default:
1177 ✗ throwStreamPrint(threadData, "unrecognized nonlinear solver");
1178 }
1179
1180 ✗ return success;
1181 }
1182
1183 /*! \fn Solve all non-linear systems in data->simulationInfo->nonlinearSystemData.
1184 *
1185 * @param data Runtime data struct.
1186 * @param threadData Thread data for error handling.
1187 * @param sysNumber Index of corresponding non-linear system
1188 *
1189 * @author ptaeuber
1190 */
1191 ✗ int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
1192 {
1193 ✗ assertStreamPrint(NULL, NULL != threadData, "threadData is NULL. Something went horribly wrong!");
1194 ✗ assertStreamPrint(threadData, NULL != data, "data is NULL. Something went horribly wrong!");
1195
1196 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
1197 int saveJumpState, constraintsSatisfied;
1198 ✗ NONLINEAR_SYSTEM_DATA* nonlinsys = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
1199 ✗ int casualTearingSet = nonlinsys->strictTearingFunctionCall != NULL;
1200 int step;
1201 int equidistantHomotopy = 0;
1202 int solveWithHomotopySolver = 0;
1203 int homotopyDeactivated = 0;
1204 int j;
1205 int nlsLs;
1206 ✗ modelica_boolean kinsol = FALSE;
1207 int res;
1208 struct dataSolver *solverData;
1209 struct dataMixedSolver *mixedSolverData;
1210 char buffer[4096];
1211 ✗ FILE *pFile = NULL;
1212 ✗ double originalLambda = data->simulationInfo->lambda;
1213
1214 ✗ if (!nonlinsys->logActive) {
1215 ✗ deactivateLogging();
1216 }
1217
1218 #if !defined(OMC_MINIMAL_RUNTIME)
1219 ✗ kinsol = (nonlinsys->nlsMethod == NLS_KINSOL) || (nonlinsys->nlsMethod == NLS_KINSOL_B);
1220 #endif
1221
1222 /* enable to avoid division by zero */
1223 ✗ data->simulationInfo->noThrowDivZero = 1;
1224 ✗ data->simulationInfo->solveContinuous = 1;
1225
1226 /* performance measurement */
1227 ✗ rt_ext_tp_tick(&nonlinsys->totalTimeClock);
1228
1229 ✗ infoStreamPrint(OMC_LOG_NLS_EXTRAPOLATE, 1, "Nonlinear system " OMC_INT_FORMAT " dump OMC_LOG_NLS_EXTRAPOLATE", nonlinsys->equationIndex);
1230 /* grab the initial guess */
1231 /* if last solving is too long ago use just old values */
1232 ✗ if (fabs(data->localData[0]->timeValue - nonlinsys->lastTimeSolved) < 5*data->simulationInfo->stepSize || casualTearingSet)
1233 {
1234 ✗ getInitialGuess(nonlinsys, data->localData[0]->timeValue);
1235 }
1236 else
1237 {
1238 ✗ nonlinsys->getIterationVars(data, nonlinsys->nlsx);
1239 ✗ memcpy(nonlinsys->nlsx, nonlinsys->nlsxOld, nonlinsys->size*(sizeof(double)));
1240 }
1241 /* update non continuous */
1242 ✗ if (data->simulationInfo->discreteCall)
1243 {
1244 // TODO: constraintsSatisfied never used
1245 ✗ constraintsSatisfied = updateInnerEquation(&resUserData, sysNumber, 1);
1246 }
1247
1248 /* print debug initial information */
1249 ✗ infoStreamPrint(OMC_LOG_NLS, 1, "############ Solve nonlinear system " OMC_INT_FORMAT " at time %g ############", nonlinsys->equationIndex, data->localData[0]->timeValue);
1250 ✗ printNonLinearInitialInfo(OMC_LOG_NLS, data, nonlinsys);
1251
1252 #if !defined(OMC_MINIMAL_RUNTIME)
1253
1254 /* try */
1255 #ifndef OMC_EMCC
1256 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1257 #endif
1258
1259 /* Improve start values with newton diagnostics method */
1260 ✗ if(omc_useStream[OMC_LOG_NLS_NEWTON_DIAGNOSTICS] && data->simulationInfo->initial) {
1261 ✗ EQUATION_INFO eqInfo = modelInfoGetEquation(&data->modelData->modelDataXml, nonlinsys->equationIndex);
1262 ✗ if (eqInfo.section == EQUATION_SECTION_INIT_LAMBDA0 || (eqInfo.section == EQUATION_SECTION_INITIAL && data->callback->functionInitialEquations_lambda0 == NULL)) {
1263 ✗ newtonDiagnostics(data, threadData, sysNumber);
1264 }
1265 }
1266
1267 /*catch */
1268 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
1269 #ifndef OMC_EMCC
1270 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1271 #endif
1272
1273 #endif
1274
1275
1276 /* try */
1277 #ifndef OMC_EMCC
1278 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1279 #endif
1280
1281 // TODO: refactor this logic
1282
1283 /* handle asserts */
1284 ✗ saveJumpState = threadData->currentErrorStage;
1285 ✗ threadData->currentErrorStage = ERROR_NONLINEARSOLVER;
1286
1287 ✗ equidistantHomotopy = data->simulationInfo->initial
1288 ✗ && nonlinsys->homotopySupport
1289 ✗ && (data->callback->homotopyMethod == LOCAL_EQUIDISTANT_HOMOTOPY && init_lambda_steps >= 1);
1290
1291 solveWithHomotopySolver = data->simulationInfo->initial
1292 ✗ && nonlinsys->homotopySupport
1293 ✗ && (data->callback->homotopyMethod == GLOBAL_ADAPTIVE_HOMOTOPY || data->callback->homotopyMethod == LOCAL_ADAPTIVE_HOMOTOPY );
1294
1295 homotopyDeactivated = (!data->simulationInfo->initial // Not an initialization system
1296 ✗ || !nonlinsys->homotopySupport // There is no homotopy in this component
1297 ✗ || (data->callback->homotopyMethod == GLOBAL_EQUIDISTANT_HOMOTOPY) // Equidistant homotopy is performed globally in symbolic_initialization()
1298 ✗ || (data->callback->homotopyMethod == LOCAL_EQUIDISTANT_HOMOTOPY // Equidistant local homotopy is selected but homotopy is deactivated ...
1299 ✗ && init_lambda_steps <= 0))
1300 ✗ || data->callback->homotopyMethod == NO_HOMOTOPY ; // No Homotopy present
1301
1302 ✗ nonlinsys->solved = NLS_FAILED;
1303 ✗ nonlinsys->initHomotopy = 0;
1304
1305 /* If homotopy is deactivated in this place or flag homotopyOnFirstTry is not set,
1306 solve the system with the selected solver */
1307 ✗ if (homotopyDeactivated || !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY]) {
1308 ✗ if (solveWithHomotopySolver && kinsol) {
1309 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the local global approach in combination with KINSOL.");
1310 } else {
1311 ✗ if (!homotopyDeactivated && !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
1312 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "Try to solve nonlinear initial system %d without homotopy first.", sysNumber);
1313
1314 /* SOLVE! */
1315 ✗ nonlinsys->solved = solveNLS(data, threadData, nonlinsys);
1316 }
1317 }
1318
1319 /* The following cases are only valid for initial systems with homotopy */
1320 /* **********************************************************************/
1321
1322 /* If the adaptive local/global homotopy approach is activated and trying without homotopy failed or is not wanted,
1323 use the HOMOTOPY SOLVER */
1324 ✗ if (solveWithHomotopySolver && nonlinsys->solved != NLS_SOLVED) {
1325 ✗ if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY] && !kinsol)
1326 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Failed to solve the initial system %d without homotopy method.", sysNumber);
1327 ✗ data->simulationInfo->lambda = 0.0;
1328 ✗ if (data->callback->homotopyMethod == LOCAL_ADAPTIVE_HOMOTOPY ) {
1329 // First solve the lambda0-system separately
1330 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "Local homotopy with adaptive step size started for nonlinear system %d.", sysNumber);
1331 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 1, "homotopy process\n---------------------------");
1332 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "solve lambda0-system");
1333 ✗ nonlinsys->homotopySupport = 0;
1334 ✗ if (!kinsol) {
1335 ✗ nonlinsys->solved = solveNLS(data, threadData, nonlinsys);
1336 } else {
1337 ✗ nlsLs = data->simulationInfo->nlsLinearSolver;
1338 ✗ data->simulationInfo->nlsLinearSolver = NLS_LS_LAPACK;
1339 ✗ nonlinsys->solved = solveWithInitHomotopy(data, threadData, nonlinsys);
1340 ✗ data->simulationInfo->nlsLinearSolver = nlsLs;
1341 }
1342 ✗ nonlinsys->homotopySupport = 1;
1343 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "solving lambda0-system done with%s success\n---------------------------", nonlinsys->solved==NLS_SOLVED ? "" : " no");
1344 ✗ messageClose(OMC_LOG_INIT_HOMOTOPY);
1345 }
1346 /* SOLVE! */
1347 ✗ if (data->callback->homotopyMethod == GLOBAL_ADAPTIVE_HOMOTOPY || nonlinsys->solved == NLS_SOLVED) {
1348 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "run along the homotopy path and solve the actual system");
1349 ✗ nonlinsys->initHomotopy = 1;
1350 ✗ nonlinsys->solved = solveWithInitHomotopy(data, threadData, nonlinsys);
1351 }
1352 }
1353
1354 /* If equidistant local homotopy is activated and trying without homotopy failed or is not wanted,
1355 use EQUIDISTANT LOCAL HOMOTOPY */
1356 ✗ if (equidistantHomotopy && nonlinsys->solved != NLS_SOLVED) {
1357 ✗ if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
1358 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Failed to solve the initial system %d without homotopy method. The local homotopy method with equidistant step size is used now.", sysNumber);
1359 else
1360 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "Local homotopy with equidistant step size started for nonlinear system %d.", sysNumber);
1361 #if !defined(OMC_NO_FILESYSTEM)
1362 ✗ const char sep[] = ",";
1363 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_INIT_HOMOTOPY))
1364 {
1365 ✗ sprintf(buffer, "%s_nonlinsys%d_equidistant_local_homotopy.csv", data->modelData->modelFilePrefix, sysNumber);
1366 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "The homotopy path of system %d will be exported to %s.", sysNumber, buffer);
1367 ✗ pFile = omc_fopen(buffer, "wt");
1368 ✗ if (pFile) {
1369 fprintf(pFile, "\"sep=%s\"\n%s", sep, "\"lambda\"");
1370 ✗ for(j=0; j<nonlinsys->size; ++j)
1371 ✗ fprintf(pFile, "%s\"%s\"", sep, modelInfoGetEquation(&data->modelData->modelDataXml, nonlinsys->equationIndex).vars[j]);
1372 fprintf(pFile, "\n");
1373 }
1374 }
1375 #endif
1376
1377 ✗ for(step=0; step<=init_lambda_steps; ++step)
1378 {
1379 ✗ data->simulationInfo->lambda = ((double)step)/(init_lambda_steps);
1380
1381 ✗ if (data->simulationInfo->lambda > 1.0) {
1382 ✗ data->simulationInfo->lambda = 1.0;
1383 }
1384
1385 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "[system %d] homotopy parameter lambda = %g", sysNumber, data->simulationInfo->lambda);
1386 /* SOLVE! */
1387 ✗ nonlinsys->solved = solveNLS(data, threadData, nonlinsys);
1388 ✗ if (nonlinsys->solved != NLS_SOLVED) break;
1389
1390 #if !defined(OMC_NO_FILESYSTEM)
1391 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_INIT_HOMOTOPY))
1392 {
1393 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "[system %d] homotopy parameter lambda = %g done\n---------------------------", sysNumber, data->simulationInfo->lambda);
1394 ✗ if (pFile) {
1395 ✗ fprintf(pFile, "%.16g", data->simulationInfo->lambda);
1396 ✗ for(j=0; j<nonlinsys->size; ++j)
1397 ✗ fprintf(pFile, "%s%.16g", sep, nonlinsys->nlsx[j]);
1398 fprintf(pFile, "\n");
1399 }
1400 }
1401 #endif
1402 }
1403
1404 #if !defined(OMC_NO_FILESYSTEM)
1405 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_INIT_HOMOTOPY))
1406 {
1407 ✗ if (pFile) fclose(pFile);
1408 }
1409 #endif
1410 ✗ data->simulationInfo->homotopySteps += init_lambda_steps;
1411 }
1412
1413 /* handle asserts */
1414 ✗ threadData->currentErrorStage = saveJumpState;
1415
1416 /*catch */
1417 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
1418 #ifndef OMC_EMCC
1419 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1420 #endif
1421
1422 ✗ messageClose(OMC_LOG_NLS_EXTRAPOLATE);
1423 /* update value list database */
1424 ✗ updateInitialGuessDB(nonlinsys, data->localData[0]->timeValue, data->simulationInfo->currentContext);
1425 ✗ if (nonlinsys->solved == NLS_SOLVED)
1426 {
1427 ✗ nonlinsys->lastTimeSolved = data->localData[0]->timeValue;
1428 }
1429 ✗ printNonLinearFinishInfo(OMC_LOG_NLS, data, nonlinsys);
1430 ✗ messageClose(OMC_LOG_NLS);
1431
1432
1433 /* enable to avoid division by zero */
1434 ✗ data->simulationInfo->noThrowDivZero = 0;
1435 ✗ data->simulationInfo->solveContinuous = 0;
1436
1437 /* performance measurement and statistics */
1438 ✗ nonlinsys->totalTime += rt_ext_tp_tock(&(nonlinsys->totalTimeClock));
1439 ✗ nonlinsys->numberOfCall++;
1440
1441 /* write csv file for debugging */
1442 #if !defined(OMC_MINIMAL_RUNTIME)
1443 ✗ if (data->simulationInfo->nlsCsvInfomation)
1444 {
1445 ✗ print_csvLineCallStats(((struct csvStats*) nonlinsys->csvData)->callStats,
1446 nonlinsys->numberOfCall,
1447 ✗ data->localData[0]->timeValue,
1448 ✗ nonlinsys->numberOfIterations,
1449 ✗ nonlinsys->numberOfFEval,
1450 nonlinsys->totalTime,
1451 nonlinsys->solved
1452 );
1453 }
1454 #endif
1455 ✗ res = check_nonlinear_solution(data, 1, sysNumber);
1456 ✗ data->simulationInfo->lambda = originalLambda;
1457
1458 ✗ if (!nonlinsys->logActive) {
1459 ✗ reactivateLogging();
1460 }
1461
1462 ✗ return res;
1463 }
1464
1465 /*! \fn check_nonlinear_solutions
1466 *
1467 * This function check whether some of non-linear systems
1468 * are failed to solve. If one is failed it returns 1 otherwise 0.
1469 *
1470 * \param [in] [data]
1471 * \param [in] [printFailingSystems]
1472 * \param [out] [returnValue] It returns >0 if fail otherwise 0.
1473 *
1474 * \author wbraun
1475 */
1476 2 int check_nonlinear_solutions(DATA *data, int printFailingSystems)
1477 {
1478 long i;
1479
1480
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 for(i=0; i<data->modelData->nNonLinearSystems; ++i) {
1481 ✗ if(check_nonlinear_solution(data, printFailingSystems, i))
1482 return 1;
1483 }
1484
1485 return 0;
1486 }
1487
1488 /*! \fn check_nonlinear_solution
1489 *
1490 * This function checks if a non-linear system
1491 * is solved. Returns a warning and 1 in case it's not
1492 * solved otherwise 0.
1493 *
1494 * \param [in] [data]
1495 * \param [in] [printFailingSystems]
1496 * \param [in] [sysNumber] index of corresponding non-linear System
1497 * \param [out] [returnValue] Returns 1 if fail otherwise 0.
1498 *
1499 * \author wbraun
1500 */
1501 ✗ int check_nonlinear_solution(DATA *data, int printFailingSystems, int sysNumber)
1502 {
1503 ✗ NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
1504 long j;
1505 int i = sysNumber;
1506
1507 const size_t buff_size = 2048;
1508 char *start_buffer;
1509 char *nominal_buffer;
1510
1511 ✗ if(nonlinsys[i].solved == NLS_FAILED)
1512 {
1513 ✗ int index = nonlinsys[i].equationIndex, indexes[2] = {1,index};
1514 ✗ if (!printFailingSystems) return 1;
1515 ✗ warningStreamPrintWithEquationIndexes(OMC_LOG_NLS, omc_dummyFileInfo, 0, indexes, "nonlinear system %d fails: at t=%g", index, data->localData[0]->timeValue);
1516 ✗ if(data->simulationInfo->initial)
1517 {
1518 ✗ warningStreamPrint(OMC_LOG_INIT, 1, "The system might not be able to initialize because the iteration variables listed below have no suitable start values. The model was probably developed with another tool that selects different iteration (tearing) variables, so its start values may apply to other variables than the ones OpenModelica iterates on here. Try providing start values for the iteration variables below, or a different tearing method (e.g. --tearingMethod=omcTearing).");
1519 }
1520
1521 ✗ start_buffer = (char*) malloc(buff_size * sizeof(char));
1522 ✗ assertStreamPrint(NULL, start_buffer != NULL, "Out of memory.");
1523 ✗ nominal_buffer = (char*) malloc(buff_size * sizeof(char));
1524 ✗ assertStreamPrint(NULL, nominal_buffer != NULL, "Out of memory.");
1525 ✗ for(j=0; j<modelInfoGetEquation(&data->modelData->modelDataXml, (nonlinsys[i]).equationIndex).numVar; ++j)
1526 {
1527 int done=0;
1528 long k;
1529 ✗ const MODEL_DATA *mData = data->modelData;
1530 ✗ for(k=0; k<mData->nVariablesRealArray && !done; ++k)
1531 {
1532 ✗ if (!strcmp(mData->realVarsData[k].info.name, modelInfoGetEquation(&data->modelData->modelDataXml, (nonlinsys[i]).equationIndex).vars[j]))
1533 {
1534 done = 1;
1535 ✗ real_vector_to_string(&mData->realVarsData[k].attribute.start, mData->realVarsData[k].dimension.numberOfDimensions == 0, start_buffer, buff_size);
1536 ✗ real_vector_to_string(&mData->realVarsData[k].attribute.nominal, mData->realVarsData[k].dimension.numberOfDimensions == 0, nominal_buffer, buff_size);
1537 ✗ warningStreamPrint(OMC_LOG_INIT, 0, "[%ld] Real %s(start=%s, nominal=%s)",
1538 j+1,
1539 ✗ mData->realVarsData[k].info.name,
1540 start_buffer,
1541 nominal_buffer);
1542 }
1543 }
1544 ✗ if (!done)
1545 {
1546 ✗ warningStreamPrint(OMC_LOG_INIT, 0, "[%ld] Real %s(start=?, nominal=?)",
1547 j+1,
1548 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,
1549 ✗ (nonlinsys[i]).equationIndex).vars[j]);
1550 }
1551 }
1552 ✗ free(start_buffer);
1553 ✗ free(nominal_buffer);
1554
1555 ✗ if(data->simulationInfo->initial)
1556 {
1557 ✗ messageCloseWarning(OMC_LOG_INIT);
1558 }
1559 ✗ return 1;
1560 }
1561 ✗ if(nonlinsys[i].solved == NLS_SOLVED_LESS_ACCURACY)
1562 {
1563 ✗ nonlinsys[i].solved = NLS_SOLVED;
1564 ✗ return 2;
1565 }
1566
1567
1568 return 0;
1569 }
1570
1571 /*! \fn cleanUpOldValueListAfterEvent
1572 *
1573 * This function clean old value list up to parameter time for all
1574 * non-linear systems.
1575 *
1576 * \param [in] [data]
1577 * \param [in] [time]
1578 *
1579 * \author wbraun
1580 */
1581 ✗ void cleanUpOldValueListAfterEvent(DATA *data, double time)
1582 {
1583 long i;
1584 ✗ NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
1585
1586 ✗ for(i=0; i<data->modelData->nNonLinearSystems; ++i) {
1587 ✗ cleanValueListbyTime(nonlinsys[i].oldValueList->valueList, time);
1588 }
1589 ✗ }
1590