Bug Summary

File:OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSystem.c
Warning:line 985, column 5
Value stored to 'constraintsSatisfied' is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2014, 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 THE BSD NEW LICENSE OR THE
11 * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
14 * 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: http://www.openmodelica.org or
19 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica
20 * distribution. GNU version 3 is obtained from:
21 * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from:
22 * http://www.opensource.org/licenses/BSD-3-Clause.
23 *
24 * This program is distributed WITHOUT ANY WARRANTY; without even the implied
25 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS
26 * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE
27 * CONDITIONS OF OSMC-PL.
28 *
29 */
30
31/*! \file nonlinearSystem.c
32 */
33
34#include <math.h>
35#include <string.h>
36
37#include "../../util/simulation_options.h"
38#include "../../util/omc_error.h"
39#include "../../util/omc_file.h"
40#include "nonlinearSystem.h"
41#include "nonlinearValuesList.h"
42#if !defined(OMC_MINIMAL_RUNTIME1)
43#include "kinsolSolver.h"
44#include "nonlinearSolverHybrd.h"
45#include "nonlinearSolverNewton.h"
46#include "newtonIteration.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/* for try and catch simulationJumpBuffer */
55#include "../../meta/meta_modelica.h"
56
57int check_nonlinear_solution(DATA *data, int printFailingSystems, int sysNumber);
58
59extern int init_lambda_steps;
60
61struct dataSolver
62{
63 void* ordinaryData;
64 void* initHomotopyData;
65};
66
67struct dataMixedSolver
68{
69 void* newtonHomotopyData;
70 void* hybridData;
71};
72
73#if !defined(OMC_MINIMAL_RUNTIME1)
74#include "../../util/write_csv.h"
75/*! \fn int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData)
76 *
77 * This function initializes csv files for analysis propose.
78 *
79 * \param [ref] [data]
80 * \param [ref] [systemData]
81 */
82int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData)
83{
84 struct csvStats* stats = (struct csvStats*) malloc(sizeof(struct csvStats));
85 char buffer[100];
86 sprintf(buffer, "%s_NLS%dStatsCall.csv", data->modelData->modelFilePrefix, (int)systemData->equationIndex);
87 stats->callStats = omc_write_csv_init(buffer, ',', '"');
88
89 sprintf(buffer, "%s_NLS%dStatsIter.csv", data->modelData->modelFilePrefix, (int)systemData->equationIndex);
90 stats->iterStats = omc_write_csv_init(buffer, ',', '"');
91
92 systemData->csvData = stats;
93
94 return 0;
95}
96#else
97int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData)
98{
99 fprintf(stderrstderr, "initializeNLScsvData not implemented for OMC_MINIMAL_RUNTIME");
100 abort();
101}
102#endif
103
104#if !defined(OMC_MINIMAL_RUNTIME1)
105/*! \fn int print_csvLineCallStatsHeader(OMC_WRITE_CSV* csvData)
106 *
107 * This function initializes csv files for analysis propose.
108 *
109 * \param [ref] [data]
110 * \param [ref] [systemData]
111 */
112int print_csvLineCallStatsHeader(OMC_WRITE_CSV* csvData)
113{
114 char buffer[1024];
115 buffer[0] = 0;
116
117 /* number of call */
118 sprintf(buffer,"numberOfCall");
119 omc_write_csv(csvData, buffer);
120 fputc(csvData->seperator,csvData->handle);
121
122 /* simulation time */
123 sprintf(buffer,"simulationTime");
124 omc_write_csv(csvData, buffer);
125 fputc(csvData->seperator,csvData->handle);
126
127 /* solving iterations */
128 sprintf(buffer,"iterations");
129 omc_write_csv(csvData, buffer);
130 fputc(csvData->seperator,csvData->handle);
131
132 /* solving fCalls */
133 sprintf(buffer,"numberOfFunctionCall");
134 omc_write_csv(csvData, buffer);
135 fputc(csvData->seperator,csvData->handle);
136
137 /* solving Time */
138 sprintf(buffer,"solvingTime");
139 omc_write_csv(csvData, buffer);
140 fputc(csvData->seperator,csvData->handle);
141
142 /* solving Time */
143 sprintf(buffer,"solvedSystem");
144 omc_write_csv(csvData, buffer);
145
146 /* finish line */
147 fputc('\n',csvData->handle);
148
149 return 0;
150}
151
152/*! \fn int print_csvLineCallStats(OMC_WRITE_CSV* csvData)
153 *
154 * This function initializes csv files for analysis propose.
155 *
156 * \param [ref] [csvData]
157 * \param [in] [number of calls]
158 * \param [in] [simulation time]
159 * \param [in] [iterations]
160 * \param [in] [number of function call]
161 * \param [in] [solving time]
162 * \param [in] [solved system]
163 */
164int print_csvLineCallStats(OMC_WRITE_CSV* csvData, int num, double time,
165 int iterations, int fCalls, double solvingTime,
166 int solved)
167{
168 char buffer[1024];
169
170 /* number of call */
171 sprintf(buffer, "%d", num);
172 omc_write_csv(csvData, buffer);
173 fputc(csvData->seperator,csvData->handle);
174
175 /* simulation time */
176 sprintf(buffer, "%g", time);
177 omc_write_csv(csvData, buffer);
178 fputc(csvData->seperator,csvData->handle);
179
180 /* solving iterations */
181 sprintf(buffer, "%d", iterations);
182 omc_write_csv(csvData, buffer);
183 fputc(csvData->seperator,csvData->handle);
184
185 /* solving fCalls */
186 sprintf(buffer, "%d", fCalls);
187 omc_write_csv(csvData, buffer);
188 fputc(csvData->seperator,csvData->handle);
189
190 /* solving Time */
191 sprintf(buffer, "%f", solvingTime);
192 omc_write_csv(csvData, buffer);
193 fputc(csvData->seperator,csvData->handle);
194
195 /* solved system */
196 sprintf(buffer, "%s", solved?"TRUE":"FALSE");
197 omc_write_csv(csvData, buffer);
198
199 /* finish line */
200 fputc('\n',csvData->handle);
201
202 return 0;
203}
204
205/*! \fn int print_csvLineIterStatsHeader(OMC_WRITE_CSV* csvData)
206 *
207 * This function initializes csv files for analysis propose.
208 *
209 * \param [ref] [data]
210 * \param [ref] [systemData]
211 */
212int print_csvLineIterStatsHeader(DATA* data, NONLINEAR_SYSTEM_DATA* systemData, OMC_WRITE_CSV* csvData)
213{
214 char buffer[1024];
215 int j;
216 int size = modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).numVar;
217
218 /* number of call */
219 sprintf(buffer,"numberOfCall");
220 omc_write_csv(csvData, buffer);
221 fputc(csvData->seperator,csvData->handle);
222
223 /* solving iterations */
224 sprintf(buffer,"iteration");
225 omc_write_csv(csvData, buffer);
226 fputc(csvData->seperator,csvData->handle);
227
228 /* variables x */
229 for(j=0; j<size; ++j) {
230 sprintf(buffer, "%s", modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[j]);
231 omc_write_csv(csvData, buffer);
232 fputc(csvData->seperator,csvData->handle);
233 }
234
235 /* residuals */
236 for(j=0; j<size; ++j) {
237 sprintf(buffer, "r%d", j+1);
238 omc_write_csv(csvData, buffer);
239 fputc(csvData->seperator,csvData->handle);
240 }
241
242 /* delta x */
243 sprintf(buffer,"delta_x");
244 omc_write_csv(csvData, buffer);
245 fputc(csvData->seperator,csvData->handle);
246
247 /* delta x scaled */
248 sprintf(buffer,"delta_x_scaled");
249 omc_write_csv(csvData, buffer);
250 fputc(csvData->seperator,csvData->handle);
251
252 /* error in f */
253 sprintf(buffer,"error_f");
254 omc_write_csv(csvData, buffer);
255 fputc(csvData->seperator,csvData->handle);
256
257 /* error in f scaled */
258 sprintf(buffer,"error_f_scaled");
259 omc_write_csv(csvData, buffer);
260 fputc(csvData->seperator,csvData->handle);
261
262 /* damping lambda */
263 sprintf(buffer,"lambda");
264 omc_write_csv(csvData, buffer);
265
266 /* finish line */
267 fputc('\n',csvData->handle);
268
269 return 0;
270}
271
272/*! \fn int print_csvLineIterStatsHeader(OMC_WRITE_CSV* csvData)
273 *
274 * This function initializes csv files for analysis propose.
275 *
276 * \param [ref] [csvData]
277 * \param [in] [size, num, ...]
278 */
279int print_csvLineIterStats(void* voidCsvData, int size, int num,
280 int iteration, double* x, double* f, double error_f,
281 double error_fs, double delta_x, double delta_xs,
282 double lambda)
283{
284 OMC_WRITE_CSV* csvData = voidCsvData;
285 char buffer[1024];
286 int j;
287
288 /* number of call */
289 sprintf(buffer, "%d", num);
290 omc_write_csv(csvData, buffer);
291 fputc(csvData->seperator,csvData->handle);
292
293 /* simulation time */
294 sprintf(buffer, "%d", iteration);
295 omc_write_csv(csvData, buffer);
296 fputc(csvData->seperator,csvData->handle);
297
298 /* x */
299 for(j=0; j<size; ++j) {
300 sprintf(buffer, "%g", x[j]);
301 omc_write_csv(csvData, buffer);
302 fputc(csvData->seperator,csvData->handle);
303 }
304
305 /* r */
306 for(j=0; j<size; ++j) {
307 sprintf(buffer, "%g", f[j]);
308 omc_write_csv(csvData, buffer);
309 fputc(csvData->seperator,csvData->handle);
310 }
311
312 /* error_f */
313 sprintf(buffer, "%g", error_f);
314 omc_write_csv(csvData, buffer);
315 fputc(csvData->seperator,csvData->handle);
316
317 /* error_f */
318 sprintf(buffer, "%g", error_fs);
319 omc_write_csv(csvData, buffer);
320 fputc(csvData->seperator,csvData->handle);
321
322 /* delta_x */
323 sprintf(buffer, "%g", delta_x);
324 omc_write_csv(csvData, buffer);
325 fputc(csvData->seperator,csvData->handle);
326
327 /* delta_xs */
328 sprintf(buffer, "%g", delta_xs);
329 omc_write_csv(csvData, buffer);
330 fputc(csvData->seperator,csvData->handle);
331
332 /* lambda */
333 sprintf(buffer, "%g", lambda);
334 omc_write_csv(csvData, buffer);
335
336 /* finish line */
337 fputc('\n',csvData->handle);
338
339 return 0;
340}
341#endif
342
343/*! \fn int initializeNonlinearSystems(DATA *data)
344 *
345 * This function allocates memory for all nonlinear systems.
346 *
347 * \param [ref] [data]
348 */
349int initializeNonlinearSystems(DATA *data, threadData_t *threadData)
350{
351 TRACE_PUSH
352 int i,j;
353 int size, nnz;
354 NONLINEAR_SYSTEM_DATA *nonlinsys = data->simulationInfo->nonlinearSystemData;
355 struct dataSolver *solverData;
356 struct dataMixedSolver *mixedSolverData;
357
358 infoStreamPrint(LOG_NLS, 1, "initialize non-linear system solvers");
359 infoStreamPrint(LOG_NLS, 0, "%ld non-linear systems", data->modelData->nNonLinearSystems);
360 if (data->simulationInfo->nlsLinearSolver == NLS_LS_DEFAULT) {
361#if !defined(OMC_MINIMAL_RUNTIME1)
362 if (data->simulationInfo->nlsMethod == NLS_KINSOL) {
363 data->simulationInfo->nlsLinearSolver = NLS_LS_KLU;
364 } else {
365 data->simulationInfo->nlsLinearSolver = NLS_LS_LAPACK;
366 }
367#else
368 data->simulationInfo->nlsLinearSolver = NLS_LS_LAPACK;
369#endif
370 }
371
372 for(i=0; i<data->modelData->nNonLinearSystems; ++i)
373 {
374 size = nonlinsys[i].size;
375 nonlinsys[i].numberOfFEval = 0;
376 nonlinsys[i].numberOfIterations = 0;
377
378 /* check if residual function pointer are valid */
379 assertStreamPrint(threadData, ((0 != nonlinsys[i].residualFunc)) || ((nonlinsys[i].strictTearingFunctionCall != NULL) ? (0 != nonlinsys[i].strictTearingFunctionCall) : 0), "residual function pointer is invalid" )if (!(((0 != nonlinsys[i].residualFunc)) || ((nonlinsys[i].strictTearingFunctionCall
!= ((void*)0)) ? (0 != nonlinsys[i].strictTearingFunctionCall
) : 0))) {throwStreamPrint((threadData), "residual function pointer is invalid"
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "simulation/solver/nonlinearSystem.c", 379
, __extension__ __PRETTY_FUNCTION__); }));}
;
380
381 /* check if analytical jacobian is created */
382 if(nonlinsys[i].jacobianIndex != -1)
383 {
384 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[nonlinsys[i].jacobianIndex]);
385 assertStreamPrint(threadData, 0 != nonlinsys[i].analyticalJacobianColumn, "jacobian function pointer is invalid" )if (!(0 != nonlinsys[i].analyticalJacobianColumn)) {throwStreamPrint
((threadData), "jacobian function pointer is invalid"); ((void
) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else __assert_fail
("0", "simulation/solver/nonlinearSystem.c", 385, __extension__
__PRETTY_FUNCTION__); }));}
;
386 if(nonlinsys[i].initialAnalyticalJacobian(data, threadData, jacobian))
387 {
388 nonlinsys[i].jacobianIndex = -1;
389 }
390 }
391
392#if !defined(OMC_MINIMAL_RUNTIME1)
393 if (nonlinsys[i].isPatternAvailable)
394 {
395 nnz = nonlinsys[i].sparsePattern->numberOfNoneZeros;
396
397 if(nnz/(double)(size*size)<=nonlinearSparseSolverMaxDensity && size >= nonlinearSparseSolverMinSize)
398 {
399 data->simulationInfo->nlsMethod = NLS_KINSOL;
400 infoStreamPrint(LOG_STDOUT, 0, "Using sparse solver kinsol for nonlinear system %d,\nbecause density of %.2f remains under threshold of %.2f and size of %d exceeds threshold of %d.\nThe maximum density and the minimal system size for using sparse solvers can be specified\nusing the runtime flags '<-nlsMaxDensity=value>' and '<-nlsMinSize=value>'.", i, nnz/(double)(size*size), nonlinearSparseSolverMaxDensity, size, nonlinearSparseSolverMinSize);
401 }
402 }
403#endif
404
405 /* allocate system data */
406 nonlinsys[i].nlsx = (double*) malloc(size*sizeof(double));
407 nonlinsys[i].nlsxExtrapolation = (double*) malloc(size*sizeof(double));
408 nonlinsys[i].nlsxOld = (double*) malloc(size*sizeof(double));
409 nonlinsys[i].resValues = (double*) malloc(size*sizeof(double));
410
411 /* allocate value list*/
412 nonlinsys[i].oldValueList = (void*) allocValueList(1);
413
414 nonlinsys[i].lastTimeSolved = 0.0;
415
416 nonlinsys[i].nominal = (double*) malloc(size*sizeof(double));
417 nonlinsys[i].min = (double*) malloc(size*sizeof(double));
418 nonlinsys[i].max = (double*) malloc(size*sizeof(double));
419 nonlinsys[i].initializeStaticNLSData(data, threadData, &nonlinsys[i]);
420
421#if !defined(OMC_MINIMAL_RUNTIME1)
422 /* csv data call stats*/
423 if (data->simulationInfo->nlsCsvInfomation)
424 {
425 if (initializeNLScsvData(data, &nonlinsys[i]))
426 {
427 throwStreamPrint(threadData, "csvData initialization failed");
428 }
429 else
430 {
431 print_csvLineCallStatsHeader(((struct csvStats*) nonlinsys[i].csvData)->callStats);
432 print_csvLineIterStatsHeader(data, &nonlinsys[i], ((struct csvStats*) nonlinsys[i].csvData)->iterStats);
433 }
434 }
435#endif
436
437 switch(data->simulationInfo->nlsMethod)
438 {
439#if !defined(OMC_MINIMAL_RUNTIME1)
440 case NLS_HYBRID:
441 solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
442 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
443 allocateHybrdData(size-1, &(solverData->ordinaryData));
444 allocateHomotopyData(size-1, &(solverData->initHomotopyData));
445 } else {
446 allocateHybrdData(size, &(solverData->ordinaryData));
447 }
448 nonlinsys[i].solverData = (void*) solverData;
449 break;
450 case NLS_KINSOL:
451 solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
452 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
453 // Try without homotopy not supported for kinsol
454 // nlsKinsolAllocate(size-1, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
455 // solverData->ordinaryData = nonlinsys[i].solverData;
456 allocateHomotopyData(size-1, &(solverData->initHomotopyData));
457 } else {
458 nlsKinsolAllocate(size, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
459 solverData->ordinaryData = nonlinsys[i].solverData;
460 }
461 nonlinsys[i].solverData = (void*) solverData;
462 break;
463 case NLS_NEWTON:
464 solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
465 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
466 allocateNewtonData(size-1, &(solverData->ordinaryData));
467 allocateHomotopyData(size-1, &(solverData->initHomotopyData));
468 } else {
469 allocateNewtonData(size, &(solverData->ordinaryData));
470 }
471 nonlinsys[i].solverData = (void*) solverData;
472 break;
473#endif
474 case NLS_HOMOTOPY:
475 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
476 allocateHomotopyData(size-1, &nonlinsys[i].solverData);
477 } else {
478 allocateHomotopyData(size, &nonlinsys[i].solverData);
479 }
480 break;
481#if !defined(OMC_MINIMAL_RUNTIME1)
482 case NLS_MIXED:
483 mixedSolverData = (struct dataMixedSolver*) malloc(sizeof(struct dataMixedSolver));
484 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
485 allocateHomotopyData(size-1, &(mixedSolverData->newtonHomotopyData));
486 allocateHybrdData(size-1, &(mixedSolverData->hybridData));
487 } else {
488 allocateHomotopyData(size, &(mixedSolverData->newtonHomotopyData));
489 allocateHybrdData(size, &(mixedSolverData->hybridData));
490 }
491 nonlinsys[i].solverData = (void*) mixedSolverData;
492 break;
493#endif
494 default:
495 throwStreamPrint(threadData, "unrecognized nonlinear solver");
496 }
497 }
498 messageClose(LOG_NLS);
499
500 TRACE_POP
501 return 0;
502}
503
504/*! \fn int updateStaticDataOfNonlinearSystems(DATA *data)
505 *
506 * This function allocates memory for all nonlinear systems.
507 *
508 * \param [ref] [data]
509 */
510int updateStaticDataOfNonlinearSystems(DATA *data, threadData_t *threadData)
511{
512 TRACE_PUSH
513 int i;
514 int size;
515 NONLINEAR_SYSTEM_DATA *nonlinsys = data->simulationInfo->nonlinearSystemData;
516
517 infoStreamPrint(LOG_NLS, 1, "update static data of non-linear system solvers");
518
519 for(i=0; i<data->modelData->nNonLinearSystems; ++i)
520 {
521 nonlinsys[i].initializeStaticNLSData(data, threadData, &nonlinsys[i]);
522 }
523
524 messageClose(LOG_NLS);
525
526 TRACE_POP
527 return 0;
528}
529
530/*! \fn freeNonlinearSystems
531 *
532 * This function frees memory of nonlinear systems.
533 *
534 * \param [ref] [data]
535 */
536int freeNonlinearSystems(DATA *data, threadData_t *threadData)
537{
538 TRACE_PUSH
539 int i;
540 NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
541 struct csvStats* stats;
542
543 infoStreamPrint(LOG_NLS, 1, "free non-linear system solvers");
544
545 for(i=0; i<data->modelData->nNonLinearSystems; ++i)
546 {
547 free(nonlinsys[i].nlsx);
548 free(nonlinsys[i].nlsxExtrapolation);
549 free(nonlinsys[i].nlsxOld);
550 free(nonlinsys[i].resValues);
551 free(nonlinsys[i].nominal);
552 free(nonlinsys[i].min);
553 free(nonlinsys[i].max);
554 freeValueList(nonlinsys[i].oldValueList, 1);
555
556
557#if !defined(OMC_MINIMAL_RUNTIME1)
558 if (data->simulationInfo->nlsCsvInfomation)
559 {
560 stats = nonlinsys[i].csvData;
561 omc_write_csv_free(stats->callStats);
562 omc_write_csv_free(stats->iterStats);
563 }
564#endif
565 /* free solver data */
566 switch(data->simulationInfo->nlsMethod)
567 {
568#if !defined(OMC_MINIMAL_RUNTIME1)
569 case NLS_HYBRID:
570 freeHybrdData(&((struct dataSolver*) nonlinsys[i].solverData)->ordinaryData);
571 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
572 freeHomotopyData(&((struct dataSolver*) nonlinsys[i].solverData)->initHomotopyData);
573 }
574 free(nonlinsys[i].solverData);
575 break;
576 case NLS_KINSOL:
577 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
578 freeHomotopyData(&((struct dataSolver*) nonlinsys[i].solverData)->initHomotopyData);
579 } else {
580 nlsKinsolFree(&((struct dataSolver*) nonlinsys[i].solverData)->ordinaryData);
581 }
582 free(nonlinsys[i].solverData);
583 break;
584 case NLS_NEWTON:
585 freeNewtonData(&((struct dataSolver*) nonlinsys[i].solverData)->ordinaryData);
586 if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
587 freeHomotopyData(&((struct dataSolver*) nonlinsys[i].solverData)->initHomotopyData);
588 }
589 free(nonlinsys[i].solverData);
590 break;
591#endif
592 case NLS_HOMOTOPY:
593 freeHomotopyData(&nonlinsys[i].solverData);
594 free(nonlinsys[i].solverData);
595 break;
596#if !defined(OMC_MINIMAL_RUNTIME1)
597 case NLS_MIXED:
598 freeHomotopyData(&((struct dataMixedSolver*) nonlinsys[i].solverData)->newtonHomotopyData);
599 freeHybrdData(&((struct dataMixedSolver*) nonlinsys[i].solverData)->hybridData);
600 free(nonlinsys[i].solverData);
601 break;
602#endif
603 default:
604 throwStreamPrint(threadData, "unrecognized nonlinear solver");
605 }
606 }
607
608 messageClose(LOG_NLS);
609
610 TRACE_POP
611 return 0;
612}
613
614/*! \fn printNonLinearSystemSolvingStatistics
615 *
616 * This function prints memory for all non-linear systems.
617 *
618 * \param [ref] [data]
619 * [in] [sysNumber] index of corresponding non-linear system
620 */
621void printNonLinearSystemSolvingStatistics(DATA *data, int sysNumber, int logLevel)
622{
623 NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
624 infoStreamPrint(logLevel, 1, "Non-linear system %d of size %d solver statistics:", (int)nonlinsys[sysNumber].equationIndex, (int)nonlinsys[sysNumber].size);
625 infoStreamPrint(logLevel, 0, " number of calls : %ld", nonlinsys[sysNumber].numberOfCall);
626 infoStreamPrint(logLevel, 0, " number of iterations : %ld", nonlinsys[sysNumber].numberOfIterations);
627 infoStreamPrint(logLevel, 0, " number of function evaluations : %ld", nonlinsys[sysNumber].numberOfFEval);
628 infoStreamPrint(logLevel, 0, " number of jacobian evaluations : %ld", nonlinsys[sysNumber].numberOfJEval);
629 infoStreamPrint(logLevel, 0, " time of jacobian evaluations : %f", nonlinsys[sysNumber].jacobianTime);
630 infoStreamPrint(logLevel, 0, " average time per call : %f", nonlinsys[sysNumber].totalTime/nonlinsys[sysNumber].numberOfCall);
631 infoStreamPrint(logLevel, 0, " total time : %f", nonlinsys[sysNumber].totalTime);
632 messageClose(logLevel);
633}
634
635/*! \fn printNonLinearInitialInfo
636 *
637 * This function prints information of an non-linear systems before an solving step.
638 *
639 * \param [in] [logName] log level in general LOG_NLS
640 * [ref] [data]
641 * [ref] [nonlinsys] index of corresponding non-linear system
642 */
643void printNonLinearInitialInfo(int logName, DATA* data, NONLINEAR_SYSTEM_DATA *nonlinsys)
644{
645 long i;
646
647 if (!ACTIVE_STREAM(logName)(useStream[logName])) return;
648 infoStreamPrint(logName, 1, "initial variable values:");
649
650 for(i=0; i<nonlinsys->size; i++)
651 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g", i+1,
652 modelInfoGetEquation(&data->modelData->modelDataXml,nonlinsys->equationIndex).vars[i],
653 nonlinsys->nlsx[i], nonlinsys->nominal[i]);
654 messageClose(logName);
655}
656
657/*! \fn printNonLinearFinishInfo
658 *
659 * This function prints information of an non-linear systems after a solving step.
660 *
661 * \param [in] [logName] log level in general LOG_NLS
662 * [ref] [data]
663 * [ref] [nonlinsys] index of corresponding non-linear system
664 */
665void printNonLinearFinishInfo(int logName, DATA* data, NONLINEAR_SYSTEM_DATA *nonlinsys)
666{
667 long i;
668
669 if (!ACTIVE_STREAM(logName)(useStream[logName])) return;
670
671 infoStreamPrint(logName, 1, "Solution status: %s", nonlinsys->solved?"SOLVED":"FAILED");
672 infoStreamPrint(logName, 0, " number of iterations : %ld", nonlinsys->numberOfIterations);
673 infoStreamPrint(logName, 0, " number of function evaluations : %ld", nonlinsys->numberOfFEval);
674 infoStreamPrint(logName, 0, " number of jacobian evaluations : %ld", nonlinsys->numberOfJEval);
675 infoStreamPrint(logName, 0, "solution values:");
676 for(i=0; i<nonlinsys->size; i++)
677 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g", i+1,
678 modelInfoGetEquation(&data->modelData->modelDataXml,nonlinsys->equationIndex).vars[i],
679 nonlinsys->nlsx[i]);
680
681 messageClose(logName);
682}
683
684/*! \fn getInitialGuess
685 *
686 * This function writes initial guess to nonlinsys->nlsx and nonlinsys->nlsOld.
687 *
688 * \param [in] [nonlinsys]
689 * \param [in] [time] time for extrapolation
690 *
691 */
692int getInitialGuess(NONLINEAR_SYSTEM_DATA *nonlinsys, double time)
693{
694 /* value extrapolation */
695 printValuesListTimes((VALUES_LIST*)nonlinsys->oldValueList);
696 /* if list is empty use current start values */
697 if (listLen(((VALUES_LIST*)nonlinsys->oldValueList)->valueList)==0)
698 {
699 /* use old value if no values are stored in the list */
700 memcpy(nonlinsys->nlsx, nonlinsys->nlsxOld, nonlinsys->size*(sizeof(double)));
701 }
702 else
703 {
704 /* get extrapolated values */
705 getValues((VALUES_LIST*)nonlinsys->oldValueList, time, nonlinsys->nlsxExtrapolation, nonlinsys->nlsxOld);
706 memcpy(nonlinsys->nlsx, nonlinsys->nlsxOld, nonlinsys->size*(sizeof(double)));
707 }
708
709 return 0;
710}
711
712/*! \fn updateInitialGuessDB
713 *
714 * This function writes new values to solution list.
715 *
716 * \param [in] [nonlinsys]
717 * \param [in] [time] time for extrapolation
718 * \param [in] [context] current context of evaluation
719 *
720 */
721int updateInitialGuessDB(NONLINEAR_SYSTEM_DATA *nonlinsys, double time, int context)
722{
723 /* write solution to oldValue list for extrapolation */
724 if (nonlinsys->solved == 1)
725 {
726 /* do not use solution of jacobian for next extrapolation */
727 if (context < 4)
728 {
729 addListElement((VALUES_LIST*)nonlinsys->oldValueList,
730 createValueElement(nonlinsys->size, time, nonlinsys->nlsx));
731 }
732 }
733 else if (nonlinsys->solved == 2)
734 {
735 if (listLen(((VALUES_LIST*)nonlinsys->oldValueList)->valueList)>0)
736 {
737 cleanValueList((VALUES_LIST*)nonlinsys->oldValueList, NULL((void*)0));
738 }
739 /* do not use solution of jacobian for next extrapolation */
740 if (context < 4)
741 {
742 addListElement((VALUES_LIST*)nonlinsys->oldValueList,
743 createValueElement(nonlinsys->size, time, nonlinsys->nlsx));
744 }
745 }
746 messageClose(LOG_NLS_EXTRAPOLATE);
747 return 0;
748}
749
750/*! \fn updateInnerEquation
751 *
752 * This function updates inner equation with the current x.
753 *
754 * \param [ref] [data]
755 * [in] [sysNumber] index of corresponding non-linear system
756 */
757int updateInnerEquation(void **dataIn, int sysNumber, int discrete)
758{
759 DATA *data = (DATA*) ((void**)dataIn[0]);
760 threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]);
761
762 NONLINEAR_SYSTEM_DATA* nonlinsys = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
763 int success = 0;
764 int constraintViolated = 0;
765
766 /* solve non continuous at discrete points*/
767 if(discrete)
768 {
769 data->simulationInfo->solveContinuous = 0;
770 }
771
772 /* try */
773#ifndef OMC_EMCC
774 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
775#endif
776
777 /* call residual function */
778 if (nonlinsys->strictTearingFunctionCall != NULL((void*)0))
779 constraintViolated = nonlinsys->residualFuncConstraints((void*) dataIn, nonlinsys->nlsx, nonlinsys->resValues, (int*)&nonlinsys->size);
780 else
781 nonlinsys->residualFunc((void*) dataIn, nonlinsys->nlsx, nonlinsys->resValues, (int*)&nonlinsys->size);
782
783 /* replace extrapolated values by current x for discrete step */
784 memcpy(nonlinsys->nlsxExtrapolation, nonlinsys->nlsx, nonlinsys->size*(sizeof(double)));
785
786 if (!constraintViolated)
787 success = 1;
788 /*catch */
789#ifndef OMC_EMCC
790 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
791#endif
792
793 if (!success && !constraintViolated)
794 {
795 warningStreamPrint(LOG_STDOUT, 0, "Non-Linear Solver try to handle a problem with a called assert.");
796 }
797
798 if(discrete)
799 {
800 data->simulationInfo->solveContinuous = 1;
801 }
802
803 return success;
804}
805
806
807int solveNLS(DATA *data, threadData_t *threadData, int sysNumber)
808{
809 int success = 0, constraintsSatisfied = 1;
810 NONLINEAR_SYSTEM_DATA* nonlinsys = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
811 int casualTearingSet = nonlinsys->strictTearingFunctionCall != NULL((void*)0);
812 struct dataSolver *solverData;
813 struct dataMixedSolver *mixedSolverData;
814
815 /* use the selected solver for solving nonlinear system */
816 switch(data->simulationInfo->nlsMethod)
817 {
818#if !defined(OMC_MINIMAL_RUNTIME1)
819 case NLS_HYBRID:
820 solverData = nonlinsys->solverData;
821 nonlinsys->solverData = solverData->ordinaryData;
822 success = solveHybrd(data, threadData, sysNumber);
823 nonlinsys->solverData = solverData;
824 break;
825 case NLS_KINSOL:
826 solverData = nonlinsys->solverData;
827 nonlinsys->solverData = solverData->ordinaryData;
828 success = nlsKinsolSolve(data, threadData, sysNumber);
829 nonlinsys->solverData = solverData;
830 break;
831 case NLS_NEWTON:
832 solverData = nonlinsys->solverData;
833 nonlinsys->solverData = solverData->ordinaryData;
834 success = solveNewton(data, threadData, sysNumber);
835 /* check if solution process was successful, if not use alternative tearing set if available (dynamic tearing)*/
836 if (!success && casualTearingSet){
837 debugString(LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
838 success = nonlinsys->strictTearingFunctionCall(data, threadData);
839 if (success) success=2;
840 }
841 nonlinsys->solverData = solverData;
842 break;
843#endif
844 case NLS_HOMOTOPY:
845 success = solveHomotopy(data, threadData, sysNumber);
846 break;
847#if !defined(OMC_MINIMAL_RUNTIME1)
848 case NLS_MIXED:
849 mixedSolverData = nonlinsys->solverData;
850 nonlinsys->solverData = mixedSolverData->newtonHomotopyData;
851
852 /* try to solve the system if it is the strict set, only try to solve the casual set if the constraints are satisfied */
853 if ((!casualTearingSet) || constraintsSatisfied)
854 success = solveHomotopy(data, threadData, sysNumber);
855
856 /* check if solution process was successful, if not use alternative tearing set if available (dynamic tearing)*/
857 if (!success && casualTearingSet){
858 debugString(LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
859 success = nonlinsys->strictTearingFunctionCall(data, threadData);
860 if (success){
861 success=2;
862 }
863 }
864
865 if (!success) {
866 nonlinsys->solverData = mixedSolverData->hybridData;
867 success = solveHybrd(data, threadData, sysNumber);
868 }
869
870 /* update iteration variables of nonlinsys->nlsx */
871 if (success){
872 nonlinsys->getIterationVars(data, nonlinsys->nlsx);
873 }
874
875 nonlinsys->solverData = mixedSolverData;
876 break;
877#endif
878 default:
879 throwStreamPrint(threadData, "unrecognized nonlinear solver");
880 }
881
882 return success;
883}
884
885/*! \fn solve system with homotopy solver
886 *
887 * \param [in] [data]
888 * \param [in] [threadData]
889 * \param [in] [sysNumber] index of corresponding non-linear system
890 *
891 * \author ptaeuber
892 */
893int solveWithInitHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
894{
895 int success = 0;
896 NONLINEAR_SYSTEM_DATA* nonlinsys = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
897 struct dataSolver *solverData;
898 struct dataMixedSolver *mixedSolverData;
899
900 /* use the homotopy solver for solving the initial system */
901 switch(data->simulationInfo->nlsMethod)
902 {
903#if !defined(OMC_MINIMAL_RUNTIME1)
904 case NLS_HYBRID:
905 case NLS_KINSOL:
906 case NLS_NEWTON:
907 solverData = nonlinsys->solverData;
908 nonlinsys->solverData = solverData->initHomotopyData;
909 success = solveHomotopy(data, threadData, sysNumber);
910 nonlinsys->solverData = solverData;
911 break;
912#endif
913 case NLS_HOMOTOPY:
914 success = solveHomotopy(data, threadData, sysNumber);
915 break;
916#if !defined(OMC_MINIMAL_RUNTIME1)
917 case NLS_MIXED:
918 mixedSolverData = nonlinsys->solverData;
919 nonlinsys->solverData = mixedSolverData->newtonHomotopyData;
920 success = solveHomotopy(data, threadData, sysNumber);
921 nonlinsys->solverData = mixedSolverData;
922 break;
923#endif
924 default:
925 throwStreamPrint(threadData, "unrecognized nonlinear solver");
926 }
927
928 return success;
929}
930
931/*! \fn solve non-linear systems
932 *
933 * \param [in] [data]
934 * \param [in] [threadData]
935 * \param [in] [sysNumber] index of corresponding non-linear system
936 *
937 * \author ptaeuber
938 */
939int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
940{
941 void *dataAndThreadData[2] = {data, threadData};
942 int success = 0, saveJumpState, constraintsSatisfied = 1;
943 NONLINEAR_SYSTEM_DATA* nonlinsys = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
944 int casualTearingSet = nonlinsys->strictTearingFunctionCall != NULL((void*)0);
945 int step;
946 int equidistantHomotopy = 0;
947 int solveWithHomotopySolver = 0;
948 int homotopyDeactivated = 0;
949 int j;
950 int nlsLs;
951 int kinsol = 0;
952 struct dataSolver *solverData;
953 struct dataMixedSolver *mixedSolverData;
954 char buffer[4096];
955 FILE *pFile = NULL((void*)0);
956
957#if !defined(OMC_MINIMAL_RUNTIME1)
958 kinsol = (data->simulationInfo->nlsMethod == NLS_KINSOL);
959#endif
960
961 data->simulationInfo->currentNonlinearSystemIndex = sysNumber;
962
963 /* enable to avoid division by zero */
964 data->simulationInfo->noThrowDivZero = 1;
965 ((DATA*)data)->simulationInfo->solveContinuous = 1;
966
967 /* performance measurement */
968 rt_ext_tp_tick(&nonlinsys->totalTimeClock);
969
970 infoStreamPrint(LOG_NLS_EXTRAPOLATE, 1, "Nonlinear system %ld dump LOG_NLS_EXTRAPOLATE", nonlinsys->equationIndex);
971 /* grab the initial guess */
972 /* if last solving is too long ago use just old values */
973 if (fabs(data->localData[0]->timeValue - nonlinsys->lastTimeSolved) < 5*data->simulationInfo->stepSize || casualTearingSet)
974 {
975 getInitialGuess(nonlinsys, data->localData[0]->timeValue);
976 }
977 else
978 {
979 nonlinsys->getIterationVars(data, nonlinsys->nlsx);
980 memcpy(nonlinsys->nlsx, nonlinsys->nlsxOld, nonlinsys->size*(sizeof(double)));
981 }
982 /* update non continuous */
983 if (data->simulationInfo->discreteCall)
984 {
985 constraintsSatisfied = updateInnerEquation(dataAndThreadData, sysNumber, 1);
Value stored to 'constraintsSatisfied' is never read
986 }
987
988 /* print debug initial information */
989 infoStreamPrint(LOG_NLS, 1, "############ Solve nonlinear system %ld at time %g ############", nonlinsys->equationIndex, data->localData[0]->timeValue);
990 printNonLinearInitialInfo(LOG_NLS, data, nonlinsys);
991
992 /* try */
993#ifndef OMC_EMCC
994 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
995#endif
996
997 /* handle asserts */
998 saveJumpState = threadData->currentErrorStage;
999 threadData->currentErrorStage = ERROR_NONLINEARSOLVER;
1000
1001 equidistantHomotopy = data->simulationInfo->initial
1002 && nonlinsys->homotopySupport
1003 && (data->callback->useHomotopy == 0 && init_lambda_steps > 1);
1004
1005 solveWithHomotopySolver = data->simulationInfo->initial
1006 && nonlinsys->homotopySupport
1007 && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3);
1008
1009 homotopyDeactivated = !data->simulationInfo->initial // Not an initialization system
1010 || !nonlinsys->homotopySupport // There is no homotopy in this component
1011 || (data->callback->useHomotopy == 1) // Equidistant homotopy is performed globally in symbolic_initialization()
1012 || (data->callback->useHomotopy == 0 // Equidistant local homotopy is selected but homotopy is deactivated ...
1013 && init_lambda_steps <= 1); // ... by the number of steps
1014
1015 nonlinsys->solved = 0;
1016 nonlinsys->initHomotopy = 0;
1017
1018 /* If homotopy is deactivated in this place or flag homotopyOnFirstTry is not set,
1019 solve the system with the selected solver */
1020 if (homotopyDeactivated || !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY]) {
1021 if (solveWithHomotopySolver && kinsol) {
1022 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "Automatically set -homotopyOnFirstTry, because trying without homotopy first is not supported for the local global approach in combination with KINSOL.");
1023 } else {
1024 if (!homotopyDeactivated && !omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
1025 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "Try to solve nonlinear initial system %d without homotopy first.", sysNumber);
1026
1027 data->simulationInfo->lambda = 1.0;
1028 /* SOLVE! */
1029 nonlinsys->solved = solveNLS(data, threadData, sysNumber);
1030 }
1031 }
1032
1033 /* The following cases are only valid for initial systems with homotopy */
1034 /* **********************************************************************/
1035
1036 /* If the adaptive local/global homotopy approach is activated and trying without homotopy failed or is not wanted,
1037 use the HOMOTOPY SOLVER */
1038 if (solveWithHomotopySolver && !nonlinsys->solved) {
1039 if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY] && !kinsol)
1040 warningStreamPrint(LOG_ASSERT, 0, "Failed to solve the initial system %d without homotopy method.", sysNumber);
1041 data->simulationInfo->lambda = 0.0;
1042 if (data->callback->useHomotopy == 3) {
1043 // First solve the lambda0-system separately
1044 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "Local homotopy with adaptive step size started for nonlinear system %d.", sysNumber);
1045 infoStreamPrint(LOG_INIT_HOMOTOPY, 1, "homotopy process\n---------------------------");
1046 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "solve lambda0-system");
1047 nonlinsys->homotopySupport = 0;
1048 if (!kinsol) {
1049 nonlinsys->solved = solveNLS(data, threadData, sysNumber);
1050 } else {
1051 nlsLs = data->simulationInfo->nlsLinearSolver;
1052 data->simulationInfo->nlsLinearSolver = NLS_LS_LAPACK;
1053 nonlinsys->solved = solveWithInitHomotopy(data, threadData, sysNumber);
1054 data->simulationInfo->nlsLinearSolver = nlsLs;
1055 }
1056 nonlinsys->homotopySupport = 1;
1057 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "solving lambda0-system done with%s success\n---------------------------", nonlinsys->solved ? "" : " no");
1058 messageClose(LOG_INIT_HOMOTOPY);
1059 }
1060 /* SOLVE! */
1061 if (data->callback->useHomotopy == 2 || nonlinsys->solved) {
1062 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "run along the homotopy path and solve the actual system");
1063 nonlinsys->initHomotopy = 1;
1064 nonlinsys->solved = solveWithInitHomotopy(data, threadData, sysNumber);
1065 }
1066 }
1067
1068 /* If equidistant local homotopy is activated and trying without homotopy failed or is not wanted,
1069 use EQUIDISTANT LOCAL HOMOTOPY */
1070 if (equidistantHomotopy && !nonlinsys->solved) {
1071 if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY])
1072 warningStreamPrint(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);
1073 else
1074 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "Local homotopy with equidistant step size started for nonlinear system %d.", sysNumber);
1075#if !defined(OMC_NO_FILESYSTEM)
1076 const char sep[] = ",";
1077 if(ACTIVE_STREAM(LOG_INIT_HOMOTOPY)(useStream[LOG_INIT_HOMOTOPY]))
1078 {
1079 sprintf(buffer, "%s_nonlinsys%d_equidistant_local_homotopy.csv", data->modelData->modelFilePrefix, sysNumber);
1080 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "The homotopy path of system %d will be exported to %s.", sysNumber, buffer);
1081 pFile = omc_fopen(buffer, "wt");
1082 fprintf(pFile, "\"sep=%s\"\n%s", sep, "\"lambda\"");
1083 for(j=0; j<nonlinsys->size; ++j)
1084 fprintf(pFile, "%s\"%s\"", sep, modelInfoGetEquation(&data->modelData->modelDataXml, nonlinsys->equationIndex).vars[j]);
1085 fprintf(pFile, "\n");
1086 }
1087#endif
1088
1089 for(step=0; step<init_lambda_steps; ++step)
1090 {
1091 data->simulationInfo->lambda = ((double)step)/(init_lambda_steps-1);
1092
1093 if (data->simulationInfo->lambda > 1.0) {
1094 data->simulationInfo->lambda = 1.0;
1095 }
1096
1097 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "[system %d] homotopy parameter lambda = %g", sysNumber, data->simulationInfo->lambda);
1098 /* SOLVE! */
1099 nonlinsys->solved = solveNLS(data, threadData, sysNumber);
1100 if (!nonlinsys->solved) break;
1101
1102#if !defined(OMC_NO_FILESYSTEM)
1103 if(ACTIVE_STREAM(LOG_INIT_HOMOTOPY)(useStream[LOG_INIT_HOMOTOPY]))
1104 {
1105 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "[system %d] homotopy parameter lambda = %g done\n---------------------------", sysNumber, data->simulationInfo->lambda);
1106 fprintf(pFile, "%.16g", data->simulationInfo->lambda);
1107 for(j=0; j<nonlinsys->size; ++j)
1108 fprintf(pFile, "%s%.16g", sep, nonlinsys->nlsx[j]);
1109 fprintf(pFile, "\n");
1110 }
1111#endif
1112 }
1113
1114#if !defined(OMC_NO_FILESYSTEM)
1115 if(ACTIVE_STREAM(LOG_INIT_HOMOTOPY)(useStream[LOG_INIT_HOMOTOPY]))
1116 {
1117 fclose(pFile);
1118 }
1119#endif
1120 data->simulationInfo->homotopySteps += init_lambda_steps;
1121 }
1122
1123 /* handle asserts */
1124 threadData->currentErrorStage = saveJumpState;
1125
1126 /*catch */
1127#ifndef OMC_EMCC
1128 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1129#endif
1130
1131 /* update value list database */
1132 updateInitialGuessDB(nonlinsys, data->localData[0]->timeValue, data->simulationInfo->currentContext);
1133 if (nonlinsys->solved == 1)
1134 {
1135 nonlinsys->lastTimeSolved = data->localData[0]->timeValue;
1136 }
1137 printNonLinearFinishInfo(LOG_NLS, data, nonlinsys);
1138 messageClose(LOG_NLS);
1139
1140
1141 /* enable to avoid division by zero */
1142 data->simulationInfo->noThrowDivZero = 0;
1143 ((DATA*)data)->simulationInfo->solveContinuous = 0;
1144
1145 /* performance measurement and statistics */
1146 nonlinsys->totalTime += rt_ext_tp_tock(&(nonlinsys->totalTimeClock));
1147 nonlinsys->numberOfCall++;
1148
1149 /* write csv file for debugging */
1150#if !defined(OMC_MINIMAL_RUNTIME1)
1151 if (data->simulationInfo->nlsCsvInfomation)
1152 {
1153 print_csvLineCallStats(((struct csvStats*) nonlinsys->csvData)->callStats,
1154 nonlinsys->numberOfCall,
1155 data->localData[0]->timeValue,
1156 nonlinsys->numberOfIterations,
1157 nonlinsys->numberOfFEval,
1158 nonlinsys->totalTime,
1159 nonlinsys->solved
1160 );
1161 }
1162#endif
1163 return check_nonlinear_solution(data, 1, sysNumber);
1164}
1165
1166/*! \fn check_nonlinear_solutions
1167 *
1168 * This function check whether some of non-linear systems
1169 * are failed to solve. If one is failed it returns 1 otherwise 0.
1170 *
1171 * \param [in] [data]
1172 * \param [in] [printFailingSystems]
1173 * \param [out] [returnValue] It returns >0 if fail otherwise 0.
1174 *
1175 * \author wbraun
1176 */
1177int check_nonlinear_solutions(DATA *data, int printFailingSystems)
1178{
1179 long i;
1180
1181 for(i=0; i<data->modelData->nNonLinearSystems; ++i) {
1182 if(check_nonlinear_solution(data, printFailingSystems, i))
1183 return 1;
1184 }
1185
1186 return 0;
1187}
1188
1189/*! \fn check_nonlinear_solution
1190 *
1191 * This function checks if a non-linear system
1192 * is solved. Returns a warning and 1 in case it's not
1193 * solved otherwise 0.
1194 *
1195 * \param [in] [data]
1196 * \param [in] [printFailingSystems]
1197 * \param [in] [sysNumber] index of corresponding non-linear System
1198 * \param [out] [returnValue] Returns 1 if fail otherwise 0.
1199 *
1200 * \author wbraun
1201 */
1202int check_nonlinear_solution(DATA *data, int printFailingSystems, int sysNumber)
1203{
1204 NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
1205 long j;
1206 int i = sysNumber;
1207
1208 if(nonlinsys[i].solved == 0)
1209 {
1210 int index = nonlinsys[i].equationIndex, indexes[2] = {1,index};
1211 if (!printFailingSystems) return 1;
1212 warningStreamPrintWithEquationIndexes(LOG_NLS, 0, indexes, "nonlinear system %d fails: at t=%g", index, data->localData[0]->timeValue);
1213 if(data->simulationInfo->initial)
1214 {
1215 warningStreamPrint(LOG_INIT, 1, "proper start-values for some of the following iteration variables might help");
1216 }
1217 for(j=0; j<modelInfoGetEquation(&data->modelData->modelDataXml, (nonlinsys[i]).equationIndex).numVar; ++j) {
1218 int done=0;
1219 long k;
1220 const MODEL_DATA *mData = data->modelData;
1221 for(k=0; k<mData->nVariablesReal && !done; ++k)
1222 {
1223 if (!strcmp(mData->realVarsData[k].info.name, modelInfoGetEquation(&data->modelData->modelDataXml, (nonlinsys[i]).equationIndex).vars[j]))
1224 {
1225 done = 1;
1226 warningStreamPrint(LOG_INIT, 0, "[%ld] Real %s(start=%g, nominal=%g)", j+1,
1227 mData->realVarsData[k].info.name,
1228 mData->realVarsData[k].attribute.start,
1229 mData->realVarsData[k].attribute.nominal);
1230 }
1231 }
1232 if (!done)
1233 {
1234 warningStreamPrint(LOG_INIT, 0, "[%ld] Real %s(start=?, nominal=?)", j+1, modelInfoGetEquation(&data->modelData->modelDataXml, (nonlinsys[i]).equationIndex).vars[j]);
1235 }
1236 }
1237 messageCloseWarning(LOG_INIT);
1238 return 1;
1239 }
1240 if(nonlinsys[i].solved == 2)
1241 {
1242 nonlinsys[i].solved = 1;
1243 return 2;
1244 }
1245
1246
1247 return 0;
1248}
1249
1250/*! \fn cleanUpOldValueListAfterEvent
1251 *
1252 * This function clean old value list up to parameter time for all
1253 * non-linear systems.
1254 *
1255 * \param [in] [data]
1256 * \param [in] [time]
1257 *
1258 * \author wbraun
1259 */
1260void cleanUpOldValueListAfterEvent(DATA *data, double time)
1261{
1262 long i;
1263 NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
1264
1265 for(i=0; i<data->modelData->nNonLinearSystems; ++i) {
1266 cleanValueListbyTime(nonlinsys[i].oldValueList, time);
1267 }
1268}
1269