Bug Summary

File:OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c
Warning:line 537, column 7
Value stored to 'success' 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 nonlinear_solver.c
32 */
33
34#ifdef __cplusplus
35extern "C" {
36#endif
37
38#include <math.h>
39#include <stdlib.h>
40#include <string.h> /* memcpy */
41
42#include "../simulation_info_json.h"
43#include "../../util/omc_error.h"
44#include "../../util/varinfo.h"
45#include "model_help.h"
46#include "../../gc/omc_gc.h"
47#include "../../meta/meta_modelica.h"
48
49#include "nonlinearSystem.h"
50#include "nonlinearSolverHybrd.h"
51extern double enorm_(integer *n, double *x);
52
53struct dataAndSys {
54 DATA* data;
55 threadData_t *threadData;
56 int sysNumber;
57};
58
59static void wrapper_fvec_hybrj(const integer* n, const double* x, double* f, double* fjac, const integer* ldjac, const integer* iflag, void* data);
60
61/*! \fn allocate memory for nonlinear system solver hybrd
62 *
63 */
64int allocateHybrdData(int size, void** voiddata)
65{
66 DATA_HYBRD* data = (DATA_HYBRD*) malloc(sizeof(DATA_HYBRD));
67
68 *voiddata = (void*)data;
69 assertStreamPrint(NULL, 0 != data, "allocationHybrdData() failed!")if (!(0 != data)) {throwStreamPrint((((void*)0)), "allocationHybrdData() failed!"
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "simulation/solver/nonlinearSolverHybrd.c"
, 69, __extension__ __PRETTY_FUNCTION__); }));}
;
70
71 data->initialized = 0;
72 data->resScaling = (double*) malloc(size*sizeof(double));
73 data->fvecScaled = (double*) malloc(size*sizeof(double));
74 data->useXScaling = 1;
75 data->xScalefactors = (double*) malloc(size*sizeof(double));
76
77 data->n = size;
78 data->x = (double*) malloc((size+1)*sizeof(double));
79 data->xSave = (double*) malloc((size+1)*sizeof(double));
80 data->xScaled = (double*) malloc((size+1)*sizeof(double));
81 data->fvec = (double*) calloc(size, sizeof(double));
82 data->fvecSave = (double*) calloc(size, sizeof(double));
83 data->xtol = 1e-12;
84 data->maxfev = size*10000;
85 data->ml = size - 1;
86 data->mu = size - 1;
87 data->epsfcn = 1e-12;
88 data->diag = (double*) malloc(size*sizeof(double));
89 data->diagres = (double*) malloc(size*sizeof(double));
90 data->mode = 1;
91 data->factor = 100.0;
92 data->nprint = -1;
93 data->info = 0;
94 data->nfev = 0;
95 data->njev = 0;
96 data->fjac = (double*) calloc((size*(size+1)), sizeof(double));
97 data->fjacobian = (double*) calloc((size*(size+1)), sizeof(double));
98 data->ldfjac = size;
99 data->r__ = (double*) malloc(((size*(size+1))/2)*sizeof(double));
100 data->lr = (size*(size + 1)) / 2;
101 data->qtf = (double*) malloc(size*sizeof(double));
102 data->wa1 = (double*) malloc(size*sizeof(double));
103 data->wa2 = (double*) malloc(size*sizeof(double));
104 data->wa3 = (double*) malloc(size*sizeof(double));
105 data->wa4 = (double*) malloc(size*sizeof(double));
106
107 data->numberOfIterations = 0;
108 data->numberOfFunctionEvaluations = 0;
109
110 assertStreamPrint(NULL, 0 != *voiddata, "allocationHybrdData() voiddata failed!")if (!(0 != *voiddata)) {throwStreamPrint((((void*)0)), "allocationHybrdData() voiddata failed!"
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "simulation/solver/nonlinearSolverHybrd.c"
, 110, __extension__ __PRETTY_FUNCTION__); }));}
;
111 return 0;
112}
113
114/*! \fn free memory for nonlinear solver hybrd
115 *
116 */
117int freeHybrdData(void **voiddata)
118{
119 DATA_HYBRD* data = (DATA_HYBRD*) *voiddata;
120
121 free(data->resScaling);
122 free(data->fvecScaled);
123 free(data->xScalefactors);
124 free(data->x);
125 free(data->xSave);
126 free(data->xScaled);
127 free(data->fvec);
128 free(data->fvecSave);
129 free(data->diag);
130 free(data->diagres);
131 free(data->fjac);
132 free(data->fjacobian);
133 free(data->r__);
134 free(data->qtf);
135 free(data->wa1);
136 free(data->wa2);
137 free(data->wa3);
138 free(data->wa4);
139
140 return 0;
141}
142
143/*! \fn printVector
144 *
145 * \param [in] [vector]
146 * \param [in] [size]
147 * \param [in] [logLevel]
148 * \param [in] [name]
149 *
150 * \author wbraun
151 */
152static void printVector(const double *vector, const integer *size, const int logLevel, const char *name)
153{
154 int i;
155 if (!ACTIVE_STREAM(logLevel)(useStream[logLevel])) return;
156 infoStreamPrint(logLevel, 1, "%s", name);
157 for(i=0; i<*size; i++)
158 infoStreamPrint(logLevel, 0, "[%2d] %20.12g", i, vector[i]);
159 messageClose(logLevel);
160}
161
162/*! \fn printStatus
163 *
164 * \param [in] [solverData]
165 * \param [in] [nfunc_evals]
166 * \param [in] [xerror]
167 * \param [in] [xerror_scaled]
168 * \param [in] [logLevel]
169 *
170 * \author wbraun
171 */
172static void printStatus(DATA *data, DATA_HYBRD *solverData, int eqSystemNumber, const int *nfunc_evals, const double *xerror, const double *xerror_scaled, const int logLevel)
173{
174 long i;
175
176 if (!ACTIVE_STREAM(logLevel)(useStream[logLevel])) return;
177 infoStreamPrint(logLevel, 1, "nls status");
178
179 infoStreamPrint(logLevel, 1, "variables");
180 for(i=0; i<solverData->n; i++)
181 infoStreamPrint(logLevel, 0, "[%ld] %s = %.20e\n - scaling factor internal = %.16e\n"
182 " - scaling factor external = %.16e", i+1,
183 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
184 solverData->x[i], solverData->diag[i], solverData->xScalefactors[i]);
185 messageClose(logLevel);
186
187 infoStreamPrint(logLevel, 1, "functions");
188 for(i=0; i<solverData->n; i++)
189 infoStreamPrint(logLevel, 0, "res[%ld] = %.20e [scaling factor = %.16e]", i+1, solverData->fvec[i], solverData->resScaling[i]);
190 messageClose(logLevel);
191
192 infoStreamPrint(logLevel, 1, "statistics");
193 infoStreamPrint(logLevel, 0, "nfunc = %d\nerror = %.20e\nerror_scaled = %.20e", *nfunc_evals, *xerror, *xerror_scaled);
194 messageClose(logLevel);
195
196 messageClose(logLevel);
197
198}
199
200
201/*! \fn getAnalyticalJacobian
202 *
203 * function calculates a jacobian matrix by
204 * numerical method finite differences
205 *
206 * \param [ref] [data]
207 * \param [out] [jac]
208 *
209 * \author wbraun
210 *
211 */
212static int getNumericalJacobian(struct dataAndSys* dataAndSysNum, double* jac, const double* x, double* f)
213{
214 struct dataAndSys *dataSys = (struct dataAndSys*) dataAndSysNum;
215 NONLINEAR_SYSTEM_DATA* systemData = &(dataSys->data->simulationInfo->nonlinearSystemData[dataSys->sysNumber]);
216 DATA_HYBRD* solverData = (DATA_HYBRD*)(systemData->solverData);
217
218 double delta_h = sqrt(solverData->epsfcn);
219 double delta_hh, delta_hhh, deltaInv;
220 integer iflag = 1;
221 int i, j, l;
222
223 memcpy(solverData->xSave, x, solverData->n*sizeof(double));
224
225 for(i = 0; i < solverData->n ; ++i)
226 {
227 delta_hhh = solverData->epsfcn * f[i];
228 delta_hh = fmax(delta_h * fmax(fabs(x[i]), fabs(delta_hhh)), delta_h);
229 delta_hh = ((f[i] >= 0) ? delta_hh : -delta_hh);
230 delta_hh = x[i] + delta_hh - x[i];
231 deltaInv = 1. / delta_hh;
232 solverData->xSave[i] = x[i] + delta_hh;
233
234 infoStreamPrint(LOG_NLS_JAC, 0, "%d. %s = %f (delta_hh = %f)", i+1, modelInfoGetEquation(&dataSys->data->modelData->modelDataXml,systemData->equationIndex).vars[i], solverData->xSave[i], delta_hh);
235 wrapper_fvec_hybrj(&solverData->n, (const double*) solverData->xSave, solverData->fvecSave, solverData->fjacobian, &solverData->ldfjac, &iflag, dataSys);
236
237 for(j = 0; j < solverData->n; ++j)
238 {
239 l = i*solverData->n+j;
240 solverData->fjacobian[l] = jac[l] = (solverData->fvecSave[j] - f[j]) * deltaInv;
241 }
242 solverData->xSave[i] = x[i];
243 }
244
245 return 0;
246}
247
248/*! \fn getAnalyticalJacobian
249 *
250 * function calculates analytical jacobian
251 *
252 * \param [ref] [data]
253 * \param [out] [jac]
254 *
255 * \author wbraun
256 *
257 */
258static int getAnalyticalJacobian(struct dataAndSys* dataSys, double* jac)
259{
260 int i, j, k, l, ii;
261 DATA *data = (dataSys->data);
262 threadData_t *threadData = dataSys->threadData;
263 NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[dataSys->sysNumber]);
264 DATA_HYBRD* solverData = (DATA_HYBRD*)(systemData->solverData);
265 const int index = systemData->jacobianIndex;
266 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]);
267
268 memset(jac, 0, (solverData->n)*(solverData->n)*sizeof(double));
269 memset(solverData->fjacobian, 0, (solverData->n)*(solverData->n)*sizeof(double));
270
271 if (jacobian->constantEqns != NULL((void*)0)) {
272 jacobian->constantEqns(data, threadData, jacobian, NULL((void*)0));
273 }
274
275 for(i=0; i < jacobian->sparsePattern->maxColors; i++)
276 {
277 /* activate seed variable for the corresponding color */
278 for(ii=0; ii < jacobian->sizeCols; ii++)
279 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
280 jacobian->seedVars[ii] = 1;
281
282 ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, NULL((void*)0));
283
284 for(j=0; j<jacobian->sizeCols; j++)
285 {
286 if(jacobian->seedVars[j] == 1)
287 {
288 ii = jacobian->sparsePattern->leadindex[j];
289 while(ii < jacobian->sparsePattern->leadindex[j+1])
290 {
291 l = jacobian->sparsePattern->index[ii];
292 k = j*jacobian->sizeRows + l;
293 solverData->fjacobian[k] = jac[k] = jacobian->resultVars[l];
294 ii++;
295 };
296 }
297 /* de-activate seed variable for the corresponding color */
298 if(jacobian->sparsePattern->colorCols[j]-1 == i)
299 jacobian->seedVars[j] = 0;
300 }
301
302 }
303
304 return 0;
305}
306
307/*! \fn wrapper function of the residual Function
308 * non-linear solver calls this subroutine fcn(n, x, fvec, iflag, data)
309 *
310 *
311 */
312static void wrapper_fvec_hybrj(const integer* n, const double* x, double* f, double* fjac, const integer* ldjac, const integer* iflag, void* dataAndSysNum)
313{
314 int i,j;
315 struct dataAndSys *dataSys = (struct dataAndSys*) dataAndSysNum;
316 DATA *data = (dataSys->data);
317 void *dataAndThreadData[2] = {data, dataSys->threadData};
318 NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[dataSys->sysNumber]);
319 DATA_HYBRD* solverData = (DATA_HYBRD*)(systemData->solverData);
320 int continuous = data->simulationInfo->solveContinuous;
321
322 switch(*iflag)
323 {
324 case 1:
325 /* re-scaling x vector */
326 if(solverData->useXScaling)
327 for(i=0; i<*n; i++)
328 solverData->xScaled[i] = x[i]*solverData->xScalefactors[i];
329
330 /* debug output */
331 if(ACTIVE_STREAM(LOG_NLS_RES)(useStream[LOG_NLS_RES])) {
332 infoStreamPrint(LOG_NLS_RES, 0, "-- residual function call %d -- scaling = %d", (int)solverData->nfev, solverData->useXScaling);
333 printVector(x, n, LOG_NLS_RES, "x vector (scaled)");
334 printVector(solverData->xScaled, n, LOG_NLS_RES, "x vector");
335 }
336
337 /* call residual function */
338 if(solverData->useXScaling){
339 (systemData->residualFunc)(dataAndThreadData, (const double*) solverData->xScaled, f, (const int*)iflag);
340 } else {
341 (systemData->residualFunc)(dataAndThreadData, x, f, (const int*)iflag);
342 }
343
344 /* debug output */
345 if(ACTIVE_STREAM(LOG_NLS_RES)(useStream[LOG_NLS_RES])) {
346 printVector(f, n, LOG_NLS_RES, "residuals");
347 infoStreamPrint(LOG_NLS_RES, 0, "-- end of residual function call %d --", (int)solverData->nfev);
348 }
349
350 solverData->numberOfFunctionEvaluations++;
351 break;
352 case 2:
353 /* set residual function continuous for jacobian calculation */
354 if(continuous)
355 data->simulationInfo->solveContinuous = 0;
356
357 if(ACTIVE_STREAM(LOG_NLS_RES)(useStream[LOG_NLS_RES]))
358 infoStreamPrint(LOG_NLS_RES, 0, "-- begin calculating jacobian --");
359
360 /* performance measurement */
361 rt_ext_tp_tick(&systemData->jacobianTimeClock);
362
363 /* call apropreated jacobian function */
364 if(systemData->jacobianIndex != -1){
365 integer iflagtmp = 1;
366 wrapper_fvec_hybrj(n, x, f, fjac, ldjac, &iflagtmp, dataSys);
367
368 getAnalyticalJacobian(dataSys, fjac);
369 }
370 else{
371 getNumericalJacobian(dataSys, fjac, x, f);
372 }
373
374 /* debug output */
375 if (ACTIVE_STREAM(LOG_NLS_RES)(useStream[LOG_NLS_RES])) {
376 infoStreamPrint(LOG_NLS_RES, 0, "-- end calculating jacobian --");
377
378 if(ACTIVE_STREAM(LOG_NLS_JAC)(useStream[LOG_NLS_JAC]))
379 {
380 char *buffer = (char*)malloc(sizeof(char)*(*n)*25);
381 infoStreamPrint(LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", (int)*n, (int)*n);
382 for(i=0; i<*n; i++)
383 {
384 buffer[0] = 0;
385 for(j=0; j<*n; j++)
386 sprintf(buffer, "%s%20.12g ", buffer, fjac[i*solverData->n+j]);
387 infoStreamPrint(LOG_NLS_JAC, 0, "%s", buffer);
388 }
389 messageClose(LOG_NLS_JAC);
390 free(buffer);
391 }
392 }
393 /* reset residual function again */
394 if(continuous)
395 data->simulationInfo->solveContinuous = 1;
396
397 /* performance measurement and statistics */
398 systemData->jacobianTime += rt_ext_tp_tock(&(systemData->jacobianTimeClock));
399 systemData->numberOfJEval++;
400
401 break;
402
403 default:
404 throwStreamPrint(NULL((void*)0), "Well, this is embarrasing. The non-linear solver should never call this case.%d", (int)*iflag);
405 break;
406 }
407}
408
409/*! \fn solve non-linear system with hybrd method
410 *
411 * \param [in] [data]
412 * [sysNumber] index of the corresponing non-linear system
413 *
414 * \author wbraun
415 */
416int solveHybrd(DATA *data, threadData_t *threadData, int sysNumber)
417{
418 NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
419 DATA_HYBRD* solverData = (DATA_HYBRD*)systemData->solverData;
420 /*
421 * We are given the number of the non-linear system.
422 * We want to look it up among all equations.
423 */
424 int eqSystemNumber = systemData->equationIndex;
425
426 int i, j;
427 integer iflag = 1;
428 double xerror, xerror_scaled;
429 int success = 0;
430 double local_tol = 1e-12;
431 double initial_factor = solverData->factor;
432 int nfunc_evals = 0;
433 int continuous = 1;
434 int nonContinuousCase = 0;
435
436 int giveUp = 0;
437 int retries = 0;
438 int retries2 = 0;
439 int retries3 = 0;
440 int assertCalled = 0;
441 int assertRetries = 0;
442 int assertMessage = 0;
443
444 modelica_boolean* relationsPreBackup;
445
446 struct dataAndSys dataAndSysNumber = {data, threadData, sysNumber};
447
448 relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));
449
450 solverData->numberOfFunctionEvaluations = 0;
451
452 // Initialize lambda variable
453 if (data->simulationInfo->nonlinearSystemData[sysNumber].homotopySupport) {
454 solverData->x[solverData->n] = 1.0;
455 solverData->xSave[solverData->n] = 1.0;
456 solverData->xScaled[solverData->n] = 1.0;
457 }
458 else {
459 solverData->x[solverData->n] = 0.0;
460 solverData->xSave[solverData->n] = 0.0;
461 solverData->xScaled[solverData->n] = 0.0;
462 }
463
464 /* debug output */
465 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
466 {
467 int indexes[2] = {1,eqSystemNumber};
468 infoStreamPrintWithEquationIndexes(LOG_NLS_V, 1, indexes, "Start solving non-linear system >>%d<< using Hybrd solver at time %g", eqSystemNumber, data->localData[0]->timeValue);
469 for(i=0; i<solverData->n; i++)
470 {
471 infoStreamPrint(LOG_NLS_V, 1, "%d. %s = %f", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->nlsx[i]);
472 infoStreamPrint(LOG_NLS_V, 0, " nominal = %f\nold = %f\nextrapolated = %f",
473 systemData->nominal[i], systemData->nlsxOld[i], systemData->nlsxExtrapolation[i]);
474 messageClose(LOG_NLS_V);
475 }
476 messageClose(LOG_NLS_V);
477 }
478
479 /* set x vector */
480 if(data->simulationInfo->discreteCall)
481 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
482 else
483 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
484
485 for(i=0; i<solverData->n; i++){
486 solverData->xScalefactors[i] = fmax(fabs(solverData->x[i]), systemData->nominal[i]);
487 }
488
489 /* start solving loop */
490 while(!giveUp && !success)
491 {
492 for(i=0; i<solverData->n; i++)
493 solverData->xScalefactors[i] = fmax(fabs(solverData->x[i]), systemData->nominal[i]);
494
495 /* debug output */
496 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
497 printVector(solverData->xScalefactors, &(solverData->n), LOG_NLS_V, "scaling factors x vector");
498 printVector(solverData->x, &(solverData->n), LOG_NLS_V, "Iteration variable values");
499 }
500
501 /* Scaling x vector */
502 if(solverData->useXScaling) {
503 for(i=0; i<solverData->n; i++) {
504 solverData->x[i] = (1.0/solverData->xScalefactors[i]) * solverData->x[i];
505 }
506 }
507
508 /* debug output */
509 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
510 {
511 printVector(solverData->x, &solverData->n, LOG_NLS_V, "Iteration variable values (scaled)");
512 }
513
514 /* set residual function continuous
515 */
516 if(continuous) {
517 ((DATA*)data)->simulationInfo->solveContinuous = 1;
518 } else {
519 ((DATA*)data)->simulationInfo->solveContinuous = 0;
520 }
521
522 giveUp = 1;
523
524 /* try */
525 {
526 int success = 0;
527#ifndef OMC_EMCC
528 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) {
529#endif
530 hybrj_(wrapper_fvec_hybrj, &solverData->n, solverData->x,
531 solverData->fvec, solverData->fjac, &solverData->ldfjac, &solverData->xtol,
532 &solverData->maxfev, solverData->diag, &solverData->mode, &solverData->factor,
533 &solverData->nprint, &solverData->info, &solverData->nfev, &solverData->njev, solverData->r__,
534 &solverData->lr, solverData->qtf, solverData->wa1, solverData->wa2,
535 solverData->wa3, solverData->wa4, (void*) &dataAndSysNumber);
536
537 success = 1;
Value stored to 'success' is never read
538 if(assertCalled)
539 {
540 infoStreamPrint(LOG_NLS_V, 0, "After assertions failed, found a solution for which assertions did not fail.");
541 /* re-scaling x vector */
542 for(i=0; i<solverData->n; i++){
543 if(solverData->useXScaling)
544 systemData->nlsxOld[i] = solverData->x[i]*solverData->xScalefactors[i];
545 else
546 systemData->nlsxOld[i] = solverData->x[i];
547 }
548 }
549 assertRetries = 0;
550 assertCalled = 0;
551 success = 1;
552#ifndef OMC_EMCC
553 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
554#endif
555 /* catch */
556 if (!success)
557 {
558 if (!assertMessage)
559 {
560 if (ACTIVE_WARNING_STREAM(LOG_STDOUT)(showAllWarnings || useStream[LOG_STDOUT]))
561 {
562 if(data->simulationInfo->initial)
563 warningStreamPrint(LOG_STDOUT, 1, "While solving non-linear system an assertion failed during initialization.");
564 else
565 warningStreamPrint(LOG_STDOUT, 1, "While solving non-linear system an assertion failed at time %g.", data->localData[0]->timeValue);
566 warningStreamPrint(LOG_STDOUT, 0, "The non-linear solver tries to solve the problem that could take some time.");
567 warningStreamPrint(LOG_STDOUT, 0, "It could help to provide better start-values for the iteration variables.");
568 if (!ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
569 warningStreamPrint(LOG_STDOUT, 0, "For more information simulate with -lv LOG_NLS_V");
570 messageClose(LOG_STDOUT);
571 }
572 assertMessage = 1;
573 }
574
575 solverData->info = -1;
576 xerror_scaled = 1;
577 xerror = 1;
578 assertCalled = 1;
579 }
580 }
581
582 /* set residual function continuous */
583 if(continuous)
584 {
585 ((DATA*)data)->simulationInfo->solveContinuous = 0;
586 }
587 else
588 {
589 ((DATA*)data)->simulationInfo->solveContinuous = 1;
590 }
591
592 /* re-scaling x vector */
593 if(solverData->useXScaling)
594 for(i=0; i<solverData->n; i++)
595 solverData->x[i] = solverData->x[i]*solverData->xScalefactors[i];
596
597 /* check for proper inputs */
598 if(solverData->info == 0) {
599 printErrorEqSyst(IMPROPER_INPUT, modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber),
600 data->localData[0]->timeValue);
601 }
602
603 if(solverData->info != -1)
604 {
605 /* evaluate with discontinuities */
606 if(data->simulationInfo->discreteCall){
607 int scaling = solverData->useXScaling;
608 int success = 0;
609 if(scaling)
610 solverData->useXScaling = 0;
611
612 ((DATA*)data)->simulationInfo->solveContinuous = 0;
613
614 /* try */
615#ifndef OMC_EMCC
616 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) {
617#endif
618 wrapper_fvec_hybrj(&solverData->n, solverData->x, solverData->fvec, solverData->fjac, &solverData->ldfjac, &iflag, (void*) &dataAndSysNumber);
619 success = 1;
620#ifndef OMC_EMCC
621 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
622#endif
623 /* catch */
624 if (!success)
625 {
626 warningStreamPrint(LOG_STDOUT, 0, "Non-Linear Solver try to handle a problem with a called assert.");
627
628 solverData->info = -1;
629 xerror_scaled = 1;
630 xerror = 1;
631 assertCalled = 1;
632 }
633
634 if(scaling)
635 solverData->useXScaling = 1;
636
637 updateRelationsPre(data);
638 }
639 }
640
641 if(solverData->info != -1)
642 {
643 /* scaling residual vector */
644 {
645 int l=0;
646 for(i=0; i<solverData->n; i++){
647 solverData->resScaling[i] = 1e-16;
648 for(j=0; j<solverData->n; j++){
649 solverData->resScaling[i] = (fabs(solverData->fjacobian[l]) > solverData->resScaling[i])
650 ? fabs(solverData->fjacobian[l]) : solverData->resScaling[i];
651 l++;
652 }
653 solverData->fvecScaled[i] = solverData->fvec[i] * (1 / solverData->resScaling[i]);
654 }
655 /* debug output */
656 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
657 {
658 infoStreamPrint(LOG_NLS_V, 1, "scaling factors for residual vector");
659 for(i=0; i<solverData->n; i++)
660 {
661 infoStreamPrint(LOG_NLS_V, 1, "scaled residual [%d] : %.20e", i, solverData->fvecScaled[i]);
662 infoStreamPrint(LOG_NLS_V, 0, "scaling factor [%d] : %.20e", i, solverData->resScaling[i]);
663 messageClose(LOG_NLS_V);
664 }
665 messageClose(LOG_NLS_V);
666 }
667
668 /* debug output */
669 if(ACTIVE_STREAM(LOG_NLS_JAC)(useStream[LOG_NLS_JAC]))
670 {
671 char *buffer = (char*)malloc(sizeof(char)*solverData->n*15);
672
673 infoStreamPrint(LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", (int)solverData->n, (int)solverData->n);
674 for(i=0; i<solverData->n; i++)
675 {
676 buffer[0] = 0;
677 for(j=0; j<solverData->n; j++)
678 sprintf(buffer, "%s%10g ", buffer, solverData->fjacobian[i*solverData->n+j]);
679 infoStreamPrint(LOG_NLS_JAC, 0, "%s", buffer);
680 }
681 messageClose(LOG_NLS_JAC);
682 free(buffer);
683 }
684
685 /* check for error */
686 xerror_scaled = enorm_(&solverData->n, solverData->fvecScaled);
687 xerror = enorm_(&solverData->n, solverData->fvec);
688 }
689 }
690
691 /* reset non-contunuousCase */
692 if(nonContinuousCase && xerror > local_tol && xerror_scaled > local_tol)
693 {
694 memcpy(data->simulationInfo->relationsPre, relationsPreBackup, sizeof(modelica_boolean)*data->modelData->nRelations);
695 nonContinuousCase = 0;
696 }
697
698 if(solverData->info < 4 && xerror > local_tol && xerror_scaled > local_tol)
699 solverData->info = 4;
700
701 /* solution found */
702 if(solverData->info == 1 || xerror <= local_tol || xerror_scaled <= local_tol)
703 {
704 int scaling;
705
706 success = 1;
707 nfunc_evals += solverData->nfev;
708 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
709 {
710 int indexes[2] = {1,eqSystemNumber};
711 /* output solution */
712 infoStreamPrintWithEquationIndexes(LOG_NLS_V, 1, indexes, "solution for NLS %d at t=%g", eqSystemNumber, data->localData[0]->timeValue);
713 for(i=0; i<solverData->n; ++i)
714 {
715 infoStreamPrint(LOG_NLS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], solverData->x[i]);
716 }
717 messageClose(LOG_NLS_V);
718 }else if (ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])){
719 infoStreamPrint(LOG_NLS_V, 1, "system solved");
720 infoStreamPrint(LOG_NLS_V, 0, "%d retries\n%d restarts", retries, retries2+retries3);
721 messageClose(LOG_NLS_V);
722 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
723 }
724 scaling = solverData->useXScaling;
725 if(scaling)
726 solverData->useXScaling = 0;
727
728 /* take the solution */
729 memcpy(systemData->nlsx, solverData->x, solverData->n*(sizeof(double)));
730
731 /* try */
732 {
733 int success = 0;
734#ifndef OMC_EMCC
735 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) {
736#endif
737 wrapper_fvec_hybrj(&solverData->n, solverData->x, solverData->fvec, solverData->fjac, &solverData->ldfjac, &iflag, (void*) &dataAndSysNumber);
738 success = 1;
739#ifndef OMC_EMCC
740 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
741#endif
742 /* catch */
743 if (!success) {
744 warningStreamPrint(LOG_STDOUT, 0, "Non-Linear Solver try to handle a problem with a called assert.");
745
746 solverData->info = 4;
747 xerror_scaled = 1;
748 xerror = 1;
749 assertCalled = 1;
750 success = 0;
751 giveUp = 0;
752 }
753 }
754 if(scaling)
755 solverData->useXScaling = 1;
756 }
757 else if((solverData->info == 4 || solverData->info == 5) && assertRetries < 1+solverData->n && assertCalled)
758 {
759 /* case only used, when the Modelica code called an assert
760 * then, we try to modify start values to avoid the assert call.*/
761 int i;
762
763 memcpy(solverData->x, systemData->nlsxOld, solverData->n*(sizeof(double)));
764
765 /* set all zero values to nominal values */
766 if(assertRetries < 1)
767 {
768 for(i=0; i<solverData->n; i++)
769 {
770 if(systemData->nlsx[i] == 0)
771 {
772 systemData->nlsx[i] = systemData->nominal[i];
773 solverData->x[i] = systemData->nominal[i];
774 }
775 }
776 }
777 /* change initial guess values one by one */
778 else if(assertRetries < solverData->n+1)
779 {
780 i = assertRetries-1;
781 solverData->x[i] += 0.01*systemData->nominal[i];
782 }
783
784 giveUp = 0;
785 nfunc_evals += solverData->nfev;
786 assertRetries++;
787 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
788 {
789 infoStreamPrint(LOG_NLS_V, 0, " - try to handle a problem with a called assert vary initial value a bit. (Retry: %d)",assertRetries);
790 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
791 }
792 }
793 else if((solverData->info == 4 || solverData->info == 5) && retries < 3)
794 {
795 /* first try to decrease factor */
796
797 /* set x vector */
798 if(data->simulationInfo->discreteCall)
799 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
800 else
801 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
802
803 solverData->factor = solverData->factor / 10.0;
804
805 retries++;
806 giveUp = 0;
807 nfunc_evals += solverData->nfev;
808 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
809 {
810 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t decreasing initial step bound to %f.", solverData->factor);
811 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
812 }
813 }
814 else if((solverData->info == 4 || solverData->info == 5) && retries < 4)
815 {
816 /* try to vary the initial values */
817
818 for(i = 0; i < solverData->n; i++)
819 solverData->x[i] += systemData->nominal[i] * 0.1;
820
821 solverData->factor = initial_factor;
822 retries++;
823 giveUp = 0;
824 nfunc_evals += solverData->nfev;
825
826 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
827 {
828 infoStreamPrint(LOG_NLS_V, 0, "iteration making no progress:\t vary solution point by 1%%.");
829 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
830 }
831 }
832 else if((solverData->info == 4 || solverData->info == 5) && retries < 5)
833 {
834 /* try old values as x-Scaling factors */
835
836 /* set x vector */
837 if(data->simulationInfo->discreteCall)
838 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
839 else
840 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
841
842
843 for(i=0; i<solverData->n; i++)
844 solverData->xScalefactors[i] = fmax(fabs(systemData->nlsxOld[i]), systemData->nominal[i]);
845
846 retries++;
847 giveUp = 0;
848 nfunc_evals += solverData->nfev;
849 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
850 {
851 infoStreamPrint(LOG_NLS_V, 0, "iteration making no progress:\t try old values as scaling factors.");
852 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
853 }
854 }
855 else if((solverData->info == 4 || solverData->info == 5) && retries < 6)
856 {
857 int scaling = 0;
858 /* try to disable x-Scaling */
859
860 /* set x vector */
861 if(data->simulationInfo->discreteCall)
862 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
863 else
864 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
865
866 scaling = solverData->useXScaling;
867 if(scaling)
868 solverData->useXScaling = 0;
869
870 /* reset x-scaling factors */
871 for(i=0; i<solverData->n; i++)
872 solverData->xScalefactors[i] = fmax(fabs(solverData->x[i]), systemData->nominal[i]);
873
874 retries++;
875 giveUp = 0;
876 nfunc_evals += solverData->nfev;
877
878 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
879 {
880 infoStreamPrint(LOG_NLS_V, 0, "iteration making no progress:\t try without scaling at all.");
881 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
882 }
883 }
884 else if((solverData->info == 4 || solverData->info == 5) && retries < 7 && data->simulationInfo->discreteCall)
885 {
886 /* try to solve non-continuous
887 * work-a-round: since other wise some model does
888 * stuck in event iteration. e.g.: Modelica.Mechanics.Rotational.Examples.HeatLosses
889 */
890
891 memcpy(solverData->x, systemData->nlsxOld, solverData->n*(sizeof(double)));
892 retries++;
893
894 /* try to solve a discontinuous system */
895 continuous = 0;
896
897 nonContinuousCase = 1;
898 memcpy(relationsPreBackup, data->simulationInfo->relationsPre, sizeof(modelica_boolean)*data->modelData->nRelations);
899
900 giveUp = 0;
901 nfunc_evals += solverData->nfev;
902 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
903 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t try to solve a discontinuous system.");
904 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
905 }
906 /* Then try with old values (instead of extrapolating )*/
907 } else if((solverData->info == 4 || solverData->info == 5) && retries2 < 1) {
908 int scaling = 0;
909 /* set x vector */
910 memcpy(solverData->x, systemData->nlsxOld, solverData->n*(sizeof(double)));
911
912 scaling = solverData->useXScaling;
913 if(!scaling)
914 solverData->useXScaling = 1;
915
916 continuous = 1;
917 solverData->factor = initial_factor;
918
919 retries = 0;
920 retries2++;
921 giveUp = 0;
922 nfunc_evals += solverData->nfev;
923 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
924 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t use old values instead extrapolated.");
925 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
926 }
927 /* try to vary the initial values */
928 } else if((solverData->info == 4 || solverData->info == 5) && retries2 < 2) {
929 /* set x vector */
930 if(data->simulationInfo->discreteCall)
931 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
932 else
933 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
934 for(i = 0; i < solverData->n; i++) {
935 solverData->x[i] *= 1.01;
936 };
937
938 retries = 0;
939 retries2++;
940 giveUp = 0;
941 nfunc_evals += solverData->nfev;
942 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
943 infoStreamPrint(LOG_NLS_V, 0,
944 " - iteration making no progress:\t vary initial point by adding 1%%.");
945 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
946 }
947 /* try to vary the initial values */
948 } else if((solverData->info == 4 || solverData->info == 5) && retries2 < 3) {
949 /* set x vector */
950 if(data->simulationInfo->discreteCall)
951 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
952 else
953 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
954 for(i = 0; i < solverData->n; i++) {
955 solverData->x[i] *= 0.99;
956 };
957
958 retries = 0;
959 retries2++;
960 giveUp = 0;
961 nfunc_evals += solverData->nfev;
962 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
963 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t vary initial point by -1%%.");
964 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
965 }
966 /* try to vary the initial values */
967 } else if((solverData->info == 4 || solverData->info == 5) && retries2 < 4) {
968 /* set x vector */
969 memcpy(solverData->x, systemData->nominal, solverData->n*(sizeof(double)));
970 retries = 0;
971 retries2++;
972 giveUp = 0;
973 nfunc_evals += solverData->nfev;
974 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
975 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t try scaling factor as initial point.");
976 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
977 }
978 /* try own scaling factors */
979 } else if((solverData->info == 4 || solverData->info == 5) && retries2 < 5 && !assertCalled) {
980 /* set x vector */
981 if(data->simulationInfo->discreteCall)
982 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
983 else
984 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
985
986 for(i = 0; i < solverData->n; i++) {
987 solverData->diag[i] = fabs(solverData->resScaling[i]);
988 if(solverData->diag[i] <= 1e-16)
989 solverData->diag[i] = 1e-16;
990 }
991 retries = 0;
992 retries2++;
993 giveUp = 0;
994 solverData->mode = 2;
995 nfunc_evals += solverData->nfev;
996 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
997 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t try with own scaling factors.");
998 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
999 }
1000 /* try without internal scaling */
1001 } else if((solverData->info == 4 || solverData->info == 5) && retries3 < 1) {
1002 /* set x vector */
1003 if(data->simulationInfo->discreteCall)
1004 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
1005 else
1006 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
1007
1008 for(i = 0; i < solverData->n; i++)
1009 solverData->diag[i] = 1.0;
1010
1011 solverData->useXScaling = 1;
1012 retries = 0;
1013 retries2 = 0;
1014 retries3++;
1015 solverData->mode = 2;
1016 giveUp = 0;
1017 nfunc_evals += solverData->nfev;
1018 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
1019 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t disable solver internal scaling.");
1020 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
1021 }
1022 /* try to reduce the tolerance a bit */
1023 } else if((solverData->info == 4 || solverData->info == 5) && retries3 < 6) {
1024 /* set x vector */
1025 if(data->simulationInfo->discreteCall)
1026 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
1027 else
1028 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
1029
1030 /* reduce tolarance */
1031 local_tol = local_tol*10;
1032
1033 solverData->factor = initial_factor;
1034 solverData->mode = 1;
1035
1036 retries = 0;
1037 retries2 = 0;
1038 retries3++;
1039
1040 giveUp = 0;
1041 nfunc_evals += solverData->nfev;
1042 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
1043 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t reduce the tolerance slightly to %e.", local_tol);
1044 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
1045 }
1046 } else if(solverData->info >= 2 && solverData->info <= 5) {
1047
1048 /* while the initialization it's ok to every time a solution */
1049 if(!data->simulationInfo->initial){
1050 printErrorEqSyst(ERROR_AT_TIME, modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber), data->localData[0]->timeValue);
1051 }
1052 if (ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V])) {
1053 infoStreamPrint(LOG_NLS_V, 0, "### No Solution! ###\n after %d restarts", retries*retries2*retries3);
1054 printStatus(data, solverData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, LOG_NLS_V);
1055 }
1056 /* take the best approximation */
1057 memcpy(systemData->nlsx, solverData->x, solverData->n*(sizeof(double)));
1058 }
1059 }
1060
1061 /* reset some solving data */
1062 solverData->factor = initial_factor;
1063 solverData->mode = 1;
1064
1065 /* write statistics */
1066 systemData->numberOfFEval += solverData->numberOfFunctionEvaluations;
1067 /* iteration in hybrid are equal to the nfev numbers */
1068 systemData->numberOfIterations += nfunc_evals;
1069
1070 free(relationsPreBackup);
1071
1072 return success;
1073}
1074
1075#ifdef __cplusplus
1076}
1077#endif