Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 1136
Functions: 0.0% 0 / 0 / 55
Branches: 0.0% 0 / 0 / 669
OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.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 #if !defined(OMC_NUM_LINEAR_SYSYTEMS) || OMC_NUM_NONLINEAR_SYSTEMS>0
29
30 #if defined(__AVR__)
31 #warning "AVR CPUs are not suitable for non-linear solvers"
32 #endif
33
34 /*! \file nonlinearSolverHomotopy.c
35 * \author bbachmann
36 */
37
38 #include <math.h>
39 #include <stdlib.h>
40 #include <string.h> /* memcpy */
41
42 #include "../options.h"
43 #include "../simulation_info_json.h"
44 #include "../jacobian_util.h"
45 #include "../../util/omc_error.h"
46 #include "../../util/omc_file.h"
47 #include "../../util/varinfo.h"
48 #include "model_help.h"
49 #if !defined(OMC_MINIMAL_RUNTIME)
50 #include "../../util/write_csv.h"
51 #endif
52
53 #include "nonlinearSystem.h"
54 #include "nonlinearSolverHomotopy.h"
55 #include "nonlinearSolverHybrd.h"
56
57 #ifdef __cplusplus
58 extern "C" {
59 #endif
60
61 extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
62
63 #ifdef __cplusplus
64 }
65 #endif
66
67 /*! \typedef DATA_HOMOTOPY
68 * define memory structure for nonlinear system solver
69 * \author bbachmann
70 */
71 typedef struct DATA_HOMOTOPY
72 {
73 modelica_boolean initialized;
74
75 size_t n; /* dimension; n == size */
76 size_t m; /* dimension: m == size+1 */
77
78 double xtol_sqrd; /* tolerance for updating solution vector */
79 double ftol_sqrd; /* tolerance for accepting accuracy */
80
81 double error_f_sqrd;
82
83 double* resScaling; /* residual scaling */
84 double* fvecScaled; /* function values scaled */
85 double* hvecScaled; /* function values scaled */
86 double* dxScaled; /* scaled solution vector */
87
88 double* minValue; /* min-attribute of variable, only pointer */
89 double* maxValue; /* max-attribute of variable, only pointer */
90 double* xScaling; /* nominal-attrbute [x.nominal,lambda.nominal] with lambda.nominal=1.0 */
91
92 /* used in wrapper_*/
93 double* f1;
94 double* f2;
95 /* used for steepest descent method */
96 double* gradFx;
97
98 /* return value, if success info == 1 */
99 int info;
100 int numberOfIterations; /* over the whole simulation time */
101 int numberOfFunctionEvaluations; /* over the whole simulation time */
102 int maxNumberOfIterations; /* number of Newton steps */
103
104 /* strict tearing set or casual tearing set */
105 int casualTearingSet;
106
107 /* newton algorithm*/
108 double* x;
109 double* x0;
110 double* xStart;
111 double* x1;
112 double* finit;
113 double* fx0;
114 double* fJac; /* n times n Jacobian matrix with additional scaling row at the end */
115 double* fJacx0;
116
117 /* debug arrays */
118 double* debug_fJac;
119 double* debug_dx;
120
121 /* homotopy parameters */
122 int initHomotopy; /* homotopy method used for the initialization with lambda from the homotopy()-operator */
123 double startDirection;
124 double tau;
125 double* y0;
126 double* y1;
127 double* y2;
128 double* yt;
129 double* dy0;
130 double* dy1;
131 double* dy2;
132 double* hvec;
133 double* hJac;
134 double* hJac2;
135 double* hJacInit;
136 double* ones;
137
138 /* linear system */
139 int* indRow;
140 int* indCol;
141
142 int (*f) (struct DATA_HOMOTOPY*, double*, double*);
143 int (*f_con) (struct DATA_HOMOTOPY*, double*, double*);
144 int (*fJac_f) (struct DATA_HOMOTOPY* solverData, double* x, double* fJac);
145 int (*h_function)(struct DATA_HOMOTOPY*, double*, double*);
146 int (*hJac_dh) (struct DATA_HOMOTOPY*, double*, double*);
147
148 NLS_USERDATA* userData;
149 int eqSystemNumber;
150 double timeValue;
151 int mixedSystem;
152
153 DATA_HYBRD* dataHybrid;
154
155 } DATA_HOMOTOPY;
156
157 /**
158 * @brief Allocate memory for non-linear homotopy solver.
159 *
160 * @param size Size of non-linear system.
161 * @param userData Pointer to set NLS user data.
162 * @return DATA_HOMOTOPY* Pointer to allocated KINSOL data.
163 */
164 ✗ DATA_HOMOTOPY* allocateHomotopyData(size_t size, NLS_USERDATA* userData)
165 {
166 ✗ DATA_HOMOTOPY* homotopyData = (DATA_HOMOTOPY*) malloc(sizeof(DATA_HOMOTOPY));
167 ✗ assertStreamPrint(NULL, NULL != homotopyData, "allocationHomotopyData() failed!");
168
169 /* fJac/fJacx0/debug_fJac receive evalJacobian's dense output, which is
170 * strided by the analytic Jacobian's OWN sizeCols -- this can exceed the
171 * NLS's genuine unknown count `size` for a partial-slice Jacobian that
172 * needs extra addressable (but not genuinely-unknown) seed columns for
173 * symbolic subscript correctness (see NBJacobian.mo's
174 * partialSliceSeedCandidates whole-array fallback). Size those three
175 * buffers off the larger of the two so evalJacobian never writes past the
176 * allocation; every other field below stays sized to the genuine `size`
177 * (the Newton/homotopy iteration itself must never see phantom unknowns). */
178 ✗ size_t jacCols = size + 1;
179 ✗ if (userData != NULL && userData->analyticJacobian != NULL &&
180 ✗ (size_t)userData->analyticJacobian->sizeCols > jacCols) {
181 jacCols = (size_t)userData->analyticJacobian->sizeCols;
182 }
183
184 ✗ homotopyData->initialized = FALSE;
185 ✗ homotopyData->n = size;
186 ✗ homotopyData->m = size + 1;
187 ✗ homotopyData->xtol_sqrd = newtonXTol*newtonXTol;
188 ✗ homotopyData->ftol_sqrd = newtonFTol*newtonFTol;
189
190 ✗ homotopyData->error_f_sqrd = 0;
191
192 ✗ homotopyData->maxNumberOfIterations = size*100;
193 ✗ homotopyData->numberOfIterations = 0;
194 ✗ homotopyData->numberOfFunctionEvaluations = 0;
195
196 ✗ homotopyData->resScaling = (double*) calloc(size,sizeof(double));
197 ✗ homotopyData->fvecScaled = (double*) calloc(size,sizeof(double));
198 ✗ homotopyData->hvecScaled = (double*) calloc(size,sizeof(double));
199 ✗ homotopyData->dxScaled = (double*) calloc(size,sizeof(double));
200
201 /* indexed up to jacobian->sizeCols in getAnalyticalJacobianHomotopy's column-scaling loop */
202 ✗ homotopyData->xScaling = (double*) calloc(jacCols,sizeof(double));
203
204 ✗ homotopyData->f1 = (double*) calloc(size,sizeof(double));
205 ✗ homotopyData->f2 = (double*) calloc(size,sizeof(double));
206 ✗ homotopyData->gradFx = (double*) calloc(size,sizeof(double));
207
208 /* damped newton */
209 ✗ homotopyData->x = (double*) calloc((size+1),sizeof(double));
210 ✗ homotopyData->x0 = (double*) calloc((size+1),sizeof(double));
211 ✗ homotopyData->xStart = (double*) calloc(size,sizeof(double));
212 ✗ homotopyData->x1 = (double*) calloc((size+1),sizeof(double));
213 ✗ homotopyData->finit = (double*) calloc(size,sizeof(double));
214 ✗ homotopyData->fx0 = (double*) calloc(size,sizeof(double));
215 ✗ homotopyData->fJac = (double*) calloc((size*jacCols),sizeof(double));
216 ✗ homotopyData->fJacx0 = (double*) calloc((size*jacCols),sizeof(double));
217
218 /* debug arrays */
219 ✗ homotopyData->debug_dx = (double*) calloc(size,sizeof(double));
220 ✗ homotopyData->debug_fJac = (double*) calloc((size*jacCols),sizeof(double));
221
222 /* homotopy */
223 ✗ homotopyData->y0 = (double*) calloc((size+1),sizeof(double));
224 ✗ homotopyData->y1 = (double*) calloc((size+1),sizeof(double));
225 ✗ homotopyData->y2 = (double*) calloc((size+1),sizeof(double));
226 ✗ homotopyData->yt = (double*) calloc((size+1),sizeof(double));
227 ✗ homotopyData->dy0 = (double*) calloc((size+1),sizeof(double));
228 ✗ homotopyData->dy1 = (double*) calloc((size+homBacktraceStrategy),sizeof(double));
229 ✗ homotopyData->dy2 = (double*) calloc((size+1),sizeof(double));
230 ✗ homotopyData->hvec = (double*) calloc(size,sizeof(double));
231 ✗ homotopyData->hJac = (double*) calloc(size*jacCols,sizeof(double));
232 ✗ homotopyData->hJac2 = (double*) calloc((size+1)*(jacCols+1),sizeof(double));
233 ✗ homotopyData->hJacInit = (double*) calloc(size*jacCols,sizeof(double));
234 ✗ homotopyData->ones = (double*) calloc(size+1,sizeof(double));
235
236 /* linear system */
237 ✗ homotopyData->indRow = (int*) calloc(size+homBacktraceStrategy-1,sizeof(int));
238 ✗ homotopyData->indCol = (int*) calloc(size+homBacktraceStrategy,sizeof(int));
239
240 ✗ homotopyData->userData = userData;
241
242 ✗ homotopyData->dataHybrid = allocateHybrdData(size, userData);
243
244 ✗ return homotopyData;
245 }
246
247 /**
248 * @brief Free homotopy data.
249 *
250 * @param homotopyData Pointer to homotopy data.
251 */
252 ✗ void freeHomotopyData(DATA_HOMOTOPY* homotopyData)
253 {
254 ✗ free(homotopyData->resScaling);
255 ✗ free(homotopyData->fvecScaled);
256 ✗ free(homotopyData->hvecScaled);
257 ✗ free(homotopyData->x);
258 ✗ free(homotopyData->debug_dx);
259 ✗ free(homotopyData->finit);
260 ✗ free(homotopyData->f1);
261 ✗ free(homotopyData->f2);
262 ✗ free(homotopyData->gradFx);
263 ✗ free(homotopyData->fJac);
264 ✗ free(homotopyData->fJacx0);
265 ✗ free(homotopyData->debug_fJac);
266
267 /* damped newton */
268 ✗ free(homotopyData->x0);
269 ✗ free(homotopyData->xStart);
270 ✗ free(homotopyData->x1);
271 ✗ free(homotopyData->dxScaled);
272
273 /* homotopy */
274 ✗ free(homotopyData->fx0);
275 ✗ free(homotopyData->hvec);
276 ✗ free(homotopyData->hJac);
277 ✗ free(homotopyData->hJac2);
278 ✗ free(homotopyData->hJacInit);
279 ✗ free(homotopyData->y0);
280 ✗ free(homotopyData->y1);
281 ✗ free(homotopyData->y2);
282 ✗ free(homotopyData->yt);
283 ✗ free(homotopyData->dy0);
284 ✗ free(homotopyData->dy1);
285 ✗ free(homotopyData->dy2);
286 ✗ free(homotopyData->xScaling);
287 ✗ free(homotopyData->ones);
288
289 /* linear system */
290 ✗ free(homotopyData->indRow);
291 ✗ free(homotopyData->indCol);
292
293 /* Don't free userData here, it's done in freeHybrdData */
294 ✗ freeHybrdData(homotopyData->dataHybrid);
295
296 ✗ free(homotopyData);
297 ✗ return;
298 }
299
300 /* Prototypes for debug functions
301 * \author bbachmann
302 */
303 ✗ void printUnknowns(int logName, DATA_HOMOTOPY *solverData)
304 {
305 long i;
306 ✗ int eqSystemNumber = solverData->eqSystemNumber;
307 ✗ DATA *data = solverData->userData->data;
308
309 ✗ if (!OMC_ACTIVE_STREAM(logName)) return;
310 ✗ infoStreamPrint(logName, 1, "nls status");
311 ✗ infoStreamPrint(logName, 0, "variables");
312
313 ✗ for(i=0; i<solverData->n; i++)
314 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g\t\t min = %16.8g\t\t max = %16.8g", i+1,
315 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
316 ✗ solverData->x[i], solverData->xScaling[i], solverData->minValue[i], solverData->maxValue[i]);
317 ✗ messageClose(logName);
318 }
319
320 ✗ void printNewtonStep(int logName, DATA_HOMOTOPY *solverData)
321 {
322 long i;
323 ✗ int eqSystemNumber = solverData->eqSystemNumber;
324 ✗ DATA *data = solverData->userData->data;
325
326 ✗ if (!OMC_ACTIVE_STREAM(logName)) return;
327 ✗ infoStreamPrint(logName, 1, "newton step");
328 ✗ infoStreamPrint(logName, 0, "variables");
329
330 ✗ for(i=0; i<solverData->n; i++)
331 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t step = %16.8g\t\t old = %16.8g", i+1,
332 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
333 ✗ solverData->x1[i], solverData->dy0[i], solverData->x[i]);
334 ✗ messageClose(logName);
335 }
336
337 ✗ void printHomotopyUnknowns(int logName, DATA_HOMOTOPY *solverData)
338 {
339 long i;
340 ✗ int eqSystemNumber = solverData->eqSystemNumber;
341 ✗ DATA *data = solverData->userData->data;
342
343 ✗ if (!OMC_ACTIVE_STREAM(logName)) return;
344 ✗ infoStreamPrint(logName, 1, "homotopy status");
345 ✗ infoStreamPrint(logName, 0, "variables");
346
347 ✗ for(i=0; i<solverData->n; i++)
348 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g\t\t min = %16.8g\t\t max = %16.8g", i+1,
349 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
350 ✗ solverData->y0[i], solverData->xScaling[i], solverData->minValue[i], solverData->maxValue[i]);
351 ✗ if (solverData->initHomotopy) {
352 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g\t\t min = %16.8g\t\t max = %16.8g", i+1,
353 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
354 ✗ solverData->y0[i], solverData->xScaling[i], solverData->minValue[i], solverData->maxValue[i]);
355 }
356 else {
357 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g", i+1,
358 "LAMBDA",
359 ✗ solverData->y0[solverData->n], solverData->xScaling[solverData->n]);
360 }
361 ✗ messageClose(logName);
362 }
363
364 ✗ void printHomotopyPredictorStep(int logName, DATA_HOMOTOPY *solverData)
365 {
366 long i;
367 ✗ int eqSystemNumber = solverData->eqSystemNumber;
368 ✗ DATA *data = solverData->userData->data;
369
370 ✗ if (!OMC_ACTIVE_STREAM(logName)) return;
371 ✗ infoStreamPrint(logName, 1, "predictor status");
372 ✗ infoStreamPrint(logName, 0, "variables");
373
374 ✗ for(i=0; i<solverData->n; i++)
375 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
376 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
377 ✗ solverData->yt[i], solverData->dy0[i], solverData->y0[i], solverData->tau);
378 ✗ if (solverData->initHomotopy) {
379 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
380 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
381 ✗ solverData->yt[i], solverData->dy0[i], solverData->y0[i], solverData->tau);
382 } else {
383 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
384 "LAMBDA",
385 ✗ solverData->yt[solverData->n], solverData->dy0[i], solverData->y0[i], solverData->tau);
386 }
387 ✗ messageClose(logName);
388 }
389
390 ✗ void printHomotopyCorrectorStep(int logName, DATA_HOMOTOPY *solverData)
391 {
392 long i;
393 ✗ int eqSystemNumber = solverData->eqSystemNumber;
394 ✗ DATA *data = solverData->userData->data;
395
396 ✗ if (!OMC_ACTIVE_STREAM(logName)) return;
397 ✗ infoStreamPrint(logName, 1, "corrector status");
398 ✗ infoStreamPrint(logName, 0, "variables");
399
400 ✗ for(i=0; i<solverData->n; i++)
401 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
402 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
403 ✗ solverData->y1[i], solverData->dy1[i], solverData->yt[i], solverData->tau);
404 ✗ if (solverData->initHomotopy) {
405 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
406 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
407 ✗ solverData->y1[i], solverData->dy1[i], solverData->yt[i], solverData->tau);
408 } else {
409 ✗ infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
410 "LAMBDA",
411 ✗ solverData->y1[solverData->n], solverData->dy1[i], solverData->yt[i], solverData->tau);
412 }
413 ✗ messageClose(logName);
414 }
415
416 ✗ void debugMatrixPermutedDouble(int logName, char* matrixName, double* matrix, int n, int m, int* indRow, int* indCol)
417 {
418 ✗ if(OMC_ACTIVE_STREAM(logName))
419 {
420 int i, j;
421 int sparsity = 0;
422 ✗ char *buffer = (char*)malloc(sizeof(char)*m*20);
423
424 ✗ infoStreamPrint(logName, 1, "%s [%dx%d-dim]", matrixName, n, m);
425 ✗ for(i=0; i<n;i++)
426 {
427 char *p = buffer;
428 ✗ for(j=0; j<m; j++)
429 {
430 if (sparsity)
431 {
432 if (fabs(matrix[indRow[i] + indCol[j]*(m-1)])<1e-12)
433 p += sprintf(p, " 0");
434 else
435 p += sprintf(p, " *");
436 }
437 else
438 {
439 ✗ p += sprintf(p, " %16.8g", matrix[indRow[i] + indCol[j]*(m-1)]);
440 }
441 }
442 ✗ infoStreamPrint(logName, 0, "%s", buffer);
443 }
444 ✗ messageClose(logName);
445 ✗ free(buffer);
446 }
447 ✗ }
448
449 ✗ void debugMatrixDouble(int logName, char* matrixName, double* matrix, int n, int m)
450 {
451 ✗ if(OMC_ACTIVE_STREAM(logName))
452 {
453 int i, j;
454 int sparsity = 0;
455 ✗ char *buffer = (char*)malloc(sizeof(char)*m*20);
456
457 ✗ infoStreamPrint(logName, 1, "%s [%dx%d-dim]", matrixName, n, m);
458 ✗ for(i=0; i<n;i++)
459 {
460 char *p = buffer;
461 ✗ for(j=0; j<m; j++)
462 {
463 if (sparsity)
464 {
465 if (fabs(matrix[i + j*(m-1)])<1e-12)
466 p += sprintf(p, " 0");
467 else
468 p += sprintf(p, " *");
469 }
470 else
471 {
472 ✗ p += sprintf(p, " %16.8g", matrix[i + j*(m-1)]);
473 }
474 }
475 ✗ infoStreamPrint(logName, 0, "%s", buffer);
476 }
477 ✗ messageClose(logName);
478 ✗ free(buffer);
479 }
480 ✗ }
481
482 ✗ void debugVectorDouble(int logName, char* vectorName, double* vector, int n)
483 {
484 ✗ if(OMC_ACTIVE_STREAM(logName))
485 {
486 int i;
487 ✗ char *buffer = (char*)malloc(sizeof(char)*n*20);
488
489 ✗ infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
490 {
491 char *p = buffer;
492 ✗ if (vector[0]<-1e+300)
493 ✗ p += sprintf(p, "-INF");
494 ✗ else if (vector[0]>1e+300)
495 ✗ p += sprintf(p, "+INF");
496 else
497 ✗ p += sprintf(p, "%16.8g", vector[0]);
498 ✗ for(i=1; i<n;i++)
499 {
500 ✗ if (vector[i]<-1e+300)
501 ✗ p += sprintf(p, " -INF");
502 ✗ else if (vector[i]>1e+300)
503 ✗ p += sprintf(p, " +INF");
504 else
505 ✗ p += sprintf(p, " %16.8g", vector[i]);
506 }
507 }
508 ✗ infoStreamPrint(logName, 0, "%s", buffer);
509 ✗ messageClose(logName);
510 ✗ free(buffer);
511 }
512 ✗ }
513
514 ✗ void debugVectorBool(int logName, char* vectorName, modelica_boolean* vector, int n)
515 {
516 ✗ if(OMC_ACTIVE_STREAM(logName))
517 {
518 int i;
519 ✗ char *buffer = (char*)malloc(sizeof(char)*n*20);
520
521 ✗ infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
522 {
523 char *p = buffer;
524 if (vector[0]<-1e+300)
525 p += sprintf(p, "-INF");
526 else if (vector[0]>1e+300)
527 p += sprintf(p, "+INF");
528 else
529 ✗ p += sprintf(p, "%d", vector[0]);
530 ✗ for(i=1; i<n;i++)
531 {
532 if (vector[i]<-1e+300)
533 p += sprintf(p, " -INF");
534 else if (vector[i]>1e+300)
535 p += sprintf(p, " +INF");
536 else
537 ✗ p += sprintf(p, " %d", vector[i]);
538 }
539 }
540 ✗ infoStreamPrint(logName, 0, "%s", buffer);
541 ✗ messageClose(logName);
542 ✗ free(buffer);
543 }
544 ✗ }
545
546 ✗ void debugVectorInt(int logName, char* vectorName, int* vector, int n)
547 {
548 ✗ if(OMC_ACTIVE_STREAM(logName))
549 {
550 int i;
551 ✗ char *buffer = (char*)malloc(sizeof(char)*n*20);
552
553 ✗ infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
554 {
555 char *p = buffer;
556 if (vector[0]<-1e+300)
557 p += sprintf(p, "-INF");
558 else if (vector[0]>1e+300)
559 p += sprintf(p, "+INF");
560 else
561 ✗ p += sprintf(p, "%d", vector[0]);
562 ✗ for(i=1; i<n;i++)
563 {
564 if (vector[i]<-1e+300)
565 p += sprintf(p, " -INF");
566 else if (vector[i]>1e+300)
567 p += sprintf(p, " +INF");
568 else
569 ✗ p += sprintf(p, " %d", vector[i]);
570 }
571 }
572 ✗ infoStreamPrint(logName, 0, "%s", buffer);
573 ✗ messageClose(logName);
574 ✗ free(buffer);
575 }
576 ✗ }
577
578
579 /* Prototypes for linear algebra functions
580 * \author bbachmann
581 */
582
583 ✗ double vec2Norm(int n, double *x)
584 {
585 int i;
586 double norm=0.0;
587 ✗ for (i=0;i<n;i++)
588 ✗ norm+=x[i]*x[i];
589 ✗ return sqrt(norm);
590 }
591
592 ✗ double vec2NormSqrd(int n, double *x)
593 {
594 int i;
595 double norm=0.0;
596 ✗ for (i=0;i<n;i++)
597 ✗ norm+=x[i]*x[i];
598 ✗ return norm;
599 }
600
601 ✗ double vecMaxNorm(int n, double *x)
602 {
603 int i;
604 ✗ double norm=fabs(x[0]);
605 ✗ for (i=1;i<n;i++)
606 ✗ if (fabs(x[i])>norm)
607 norm=fabs(x[i]);
608 ✗ return norm;
609 }
610
611 /**
612 * @brief Sets all infs and nans of a vector to 1.
613 */
614 ✗ void vecMakeFinite(int n, double *a)
615 {
616 ✗ for (int i = 0; i < n; i++)
617 {
618 ✗ if (!isfinite(a[i]))
619 {
620 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Entry of scaling vector is inf or nan. Element will be set to 1.0.");
621 ✗ a[i] = 1;
622 }
623 }
624 ✗ }
625
626 ✗ void vecAdd(int n, double *a, double *b, double *c)
627 {
628 int i;
629 ✗ for (i=0;i<n;i++)
630 ✗ c[i] = a[i] + b[i];
631 ✗ }
632
633 ✗ void vecAddScal(int n, double *a, double *b, double s, double *c)
634 {
635 int i;
636 ✗ for (i=0;i<n;i++)
637 ✗ c[i] = a[i] + s*b[i];
638 ✗ }
639
640 ✗ void vecScalarMult(int n, double *a, double s, double *b)
641 {
642 int i;
643 ✗ for (i=0;i<n;i++)
644 ✗ b[i] = s*a[i];
645 ✗ }
646
647 ✗ void vecLinearComb(int n, double *a, double r, double *b, double s, double *c)
648 {
649 int i;
650 ✗ for (i=0;i<n;i++)
651 ✗ c[i] = r*a[i] + s*b[i];
652 ✗ }
653
654 ✗ void vecCopy(int n, double *a, double *b)
655 {
656 ✗ memcpy(b, a, n*(sizeof(double)));
657 ✗ }
658
659 ✗ void vecCopyBool(int n, modelica_boolean *a, modelica_boolean *b)
660 {
661 ✗ memcpy(b, a, n*(sizeof(modelica_boolean)));
662 ✗ }
663
664 ✗ void vecAddInv(int n, double *a, double *b)
665 {
666 int i;
667 ✗ for (i=0;i<n;i++)
668 ✗ b[i] = -a[i];
669 ✗ }
670
671 ✗ void vecDiff(int n, double *a, double *b, double *c)
672 {
673 int i;
674 ✗ for (i=0;i<n;i++)
675 ✗ c[i] = a[i] - b[i];
676 ✗ }
677
678 ✗ int isNotEqualVectorInt(int n, modelica_boolean *a, modelica_boolean *b)
679 {
680 int i, isNotEqual = 0;
681 ✗ for (i=0;i<n;i++)
682 ✗ isNotEqual += abs(a[i] - b[i]);
683 ✗ return isNotEqual;
684 }
685
686 ✗ void vecMultScaling(int n, double *a, double *b, double *c)
687 {
688 int i;
689 ✗ for (i=0;i<n;i++)
690 ✗ c[i] = (fabs(b[i])>0 ? a[i]*fabs(b[i]):a[i]);
691 ✗ }
692
693 ✗ void vecDivScaling(int n, double *a, double *b, double *c)
694 {
695 int i;
696 ✗ for (i=0;i<n;i++)
697 ✗ c[i] = (fabs(b[i])>0 ? a[i]/fabs(b[i]):a[i]);
698 ✗ }
699
700 ✗ void vecNormalize(int n, double *a, double *b)
701 {
702 int i;
703 ✗ double norm = vec2Norm(n,a);
704 ✗ for (i=0;i<n;i++)
705 ✗ b[i] = (norm>0 ? a[i]/norm:a[i]);
706 ✗ }
707
708 ✗ void vecConst(int n, double value, double *a)
709 {
710 int i;
711 ✗ for (i=0;i<n;i++)
712 ✗ a[i] = value;
713 ✗ }
714
715 ✗ double vecScalarProd(int n, double *a, double *b)
716 {
717 int i;
718 double prod;
719
720 ✗ for (i=0,prod=0;i<n;i++)
721 ✗ prod = prod + a[i]*b[i];
722
723 ✗ return prod;
724 }
725
726 /* Matrix has dimension [n x m], vector [m] */
727 ✗ void matVecMult(int n, int m, double *A, double *b, double *c)
728 {
729 int i, j;
730 ✗ for (i=0;i<n;i++)
731 ✗ c[i] = 0.0;
732 ✗ for (j=0;j<m;j++) {
733 ✗ for (i=0;i<n;i++)
734 ✗ c[i] += A[i+j*(m-1)]*b[j];
735 }
736 ✗ }
737
738 /* Matrix has dimension [n x m], vector [m] */
739 ✗ void matVecMultAbs(int n, int m, double *A, double *b, double *c)
740 {
741 int i, j;
742 ✗ for (i=0;i<n;i++)
743 ✗ c[i] = 0.0;
744 ✗ for (j=0;j<m;j++) {
745 ✗ for (i=0;i<n;i++)
746 ✗ c[i] += fabs(A[i+j*(m-1)]*b[j]);
747 }
748 ✗ }
749
750 /* Matrix has dimension [n x (n+1)] */
751 ✗ void matVecMultBB(int n, double *A, double *b, double *c)
752 {
753 int i, j;
754 ✗ for (i=0;i<n;i++)
755 ✗ c[i] = 0.0;
756 ✗ for (j=0;j<n;j++) {
757 ✗ for (i=0;i<n;i++)
758 ✗ c[i] += A[i+j*n]*b[j];
759 }
760 ✗ }
761
762 /* Matrix has dimension [n x (n+1)] */
763 ✗ void matVecMultAbsBB(int n, double *A, double *b, double *c)
764 {
765 int i, j;
766 ✗ for (i=0;i<n;i++)
767 ✗ c[i] = 0.0;
768 ✗ for (j=0;j<n;j++) {
769 ✗ for (i=0;i<n;i++)
770 ✗ c[i] += fabs(A[i+j*n]*b[j]);
771 }
772 ✗ }
773
774 /* Matrix has dimension [n x (n+1)] */
775 ✗ void matAddBB(int n, double* A, double* B, double* C)
776 {
777 int i, j;
778 ✗ for (j=0;j<n+1;j++) {
779 ✗ for (i=0;i<n;i++)
780 ✗ C[i + j*n] = A[i + j*n] + B[i + j*n];
781 }
782 ✗ }
783
784 /* Matrix has dimension [n x (n+1)] */
785 ✗ void matDiffBB(int n, double* A, double* B, double* C)
786 {
787 int i, j;
788 ✗ for (j=0;j<n;j++) {
789 ✗ for (i=0;i<n;i++)
790 ✗ C[i + j*n] = A[i + j*n] - B[i + j*n];
791 }
792 ✗ }
793
794 /* Matrix has dimension [n x m] */
795 ✗ void scaleMatrixRows(int n, int m, double *A)
796 {
797 const double delta = 0; /* This might be changed to sqrt(DBL_EPSILON) */
798 int i, j;
799 ✗ double* rowsMax = (double*) calloc(n,sizeof(double));
800
801 ✗ for (i=0;i<n;i++)
802 ✗ rowsMax[i] = 0;
803
804 /* find maximum of each row */
805 ✗ for (j=0;j<n;j++) {
806 ✗ for (i=0;i<n;i++) {
807 ✗ if (fabs(A[i+j*(m-1)]) > rowsMax[i]) {
808 ✗ rowsMax[i] = fabs(A[i+j*(m-1)]);
809 }
810 }
811 }
812
813 /* remove zero normailzation */
814 ✗ for (i=0;i<n;i++) {
815 ✗ if (rowsMax[i] <= delta)
816 ✗ rowsMax[i] = 1.0;
817 }
818
819 /* scale matrix */
820 ✗ for (j=0;j<m;j++) {
821 ✗ for (i=0;i<n;i++)
822 ✗ A[i+j*(m-1)] /= rowsMax[i];
823 }
824
825 ✗ free(rowsMax);
826 ✗ }
827
828 /* Build the newton matrix for the corrector step with orthogonal backtrace strategy */
829 ✗ void orthogonalBacktraceMatrix(DATA_HOMOTOPY* solverData, double* hJac, double* hvec, double* v, double* hJac2, int n, int m)
830 {
831 int i, j;
832 ✗ for (j=0; j<m; j++) {
833 ✗ for (i=0; i<n; i++) {
834 ✗ hJac2[i + j*m] = hJac[i + j*(m-1)];
835 }
836 ✗ hJac2[n + j*m] = v[j];
837 }
838 ✗ for (i=0; i<n; i++) {
839 ✗ hJac2[i + m*m] = hvec[i];
840 }
841 ✗ hJac2[n + m*m] = 0;
842 ✗ }
843
844 /*! \fn getAnalyticalJacobian
845 *
846 * function calculates analytical jacobian
847 *
848 * \param [ref] [data]
849 * \param [out] [jac]
850 *
851 * \author wbraun
852 * bbachmann: introduce scaling factor
853 *
854 */
855 ✗ int getAnalyticalJacobianHomotopy(DATA_HOMOTOPY* solverData, double* jac)
856 {
857 int j,k,l,ii;
858 ✗ DATA* data = solverData->userData->data;
859 ✗ threadData_t *threadData = solverData->userData->threadData;
860 ✗ JACOBIAN* jacobian = solverData->userData->analyticJacobian;
861 ✗ const SPARSE_PATTERN* sp = jacobian->sparsePattern;
862
863 /* call generic dense Jacobian */
864 ✗ evalJacobian(data, threadData, jacobian, NULL, jac, TRUE);
865
866 ✗ if (!sp) return 0; /* pattern removed; jac is zeroed, solver will fail numerically */
867
868 /* apply scaling to each column; must use the same row stride evalJacobian
869 * used to fill jac (min(sizeRows, sizeCols), not always sizeCols -- see
870 * jacobian_util.c:evalJacobian). Using sizeCols unconditionally here
871 * misaligns every scaling write whenever sizeCols > sizeRows (a genuinely
872 * rectangular Jacobian, not just NLS's "auxiliary rows beyond sizeCols"
873 * case), corrupting entries evalJacobian never touched while leaving the
874 * real ones unscaled. */
875 {
876 ✗ const int denseRows = jacobian->sizeRows < jacobian->sizeCols ? jacobian->sizeRows : jacobian->sizeCols;
877 ✗ for (j = 0; j < jacobian->sizeCols; j++) {
878 ✗ for (ii = sp->leadindex[j]; ii < sp->leadindex[j+1]; ii++) {
879 ✗ l = sp->index[ii];
880 ✗ if (l >= denseRows) continue; /* skip auxiliary rows */
881 ✗ k = j*denseRows + l;
882 ✗ jac[k] *= solverData->xScaling[j];
883 }
884 }
885 }
886
887 return 0;
888 }
889
890 /*! \fn getNumericalJacobianHomotopy
891 *
892 * function calculates a jacobian matrix by
893 * numerical method finite differences
894 * \author bbachmann
895 *
896 */
897 ✗ static int getNumericalJacobianHomotopy(DATA_HOMOTOPY* solverData, double *x, double *fJac)
898 {
899 const double delta_h = sqrt(DBL_EPSILON*2e1);
900 double delta_hh;
901 double xsave;
902 int i,j,l;
903 int N;
904 double* f1;
905 int (*f) (struct DATA_HOMOTOPY*, double*, double*);
906
907 ✗ if (solverData->initHomotopy) {
908 ✗ N = solverData->n + 1; /* also calculate the lambda column */
909 ✗ f1 = solverData->hvec; /* homotopy function values solverData->hvec must be set outside this function based on x */
910 ✗ f = solverData->h_function;
911 } else {
912 ✗ N = solverData->n; /* calculate jacobian without the lambda column */
913 ✗ f1 = solverData->f1; /* normal function values solverData->f1 must be set outside this function based on x */
914 ✗ f = solverData->casualTearingSet ? solverData->f_con : solverData->f;
915 }
916
917 ✗ for(i = 0; i < N; i++) {
918 ✗ xsave = x[i];
919 ✗ delta_hh = delta_h * (fabs(xsave) + 1.0);
920 ✗ if ((xsave + delta_hh >= solverData->maxValue[i]))
921 ✗ delta_hh *= -1;
922 ✗ x[i] += delta_hh;
923 /* Calculate scaled difference quotient */
924 ✗ delta_hh = 1. / delta_hh * solverData->xScaling[i];
925 ✗ f(solverData, x, solverData->f2);
926
927 ✗ for(j = 0; j < solverData->n; j++) {
928 ✗ l = i * solverData->n + j;
929 ✗ fJac[l] = (solverData->f2[j] - f1[j]) * delta_hh;
930 }
931 ✗ x[i] = xsave;
932 }
933 ✗ return 0;
934 }
935
936 /*! \fn wrapper_fvec for the residual Function
937 * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data)
938 *
939 * \author bbachmann
940 *
941 */
942 ✗ static int wrapper_fvec(DATA_HOMOTOPY* solverData, double* x, double* f)
943 {
944 ✗ DATA* data = solverData->userData->data;
945 ✗ threadData_t* threadData = solverData->userData->threadData;
946 ✗ NONLINEAR_SYSTEM_DATA* nlsData = solverData->userData->nlsData;
947 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
948 ✗ int iflag = 0;
949
950 /* TODO: change input to residualFunc from data to systemData */
951 ✗ nlsData->residualFunc(&resUserData, x, f, &iflag);
952 ✗ solverData->numberOfFunctionEvaluations++;
953
954 ✗ return 0;
955 }
956
957 /*! \fn wrapper_fvec_constraints for the residual Function
958 * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data)
959 *
960 * \author ptaeuber
961 *
962 */
963 ✗ int wrapper_fvec_constraints(DATA_HOMOTOPY* solverData, double* x, double* f)
964 {
965 ✗ DATA* data = solverData->userData->data;
966 ✗ threadData_t* threadData = solverData->userData->threadData;
967 ✗ NONLINEAR_SYSTEM_DATA* nlsData = solverData->userData->nlsData;
968 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
969 ✗ int iflag = 0;
970 int retVal;
971
972 /* TODO: change input to residualFunc from data to systemData */
973 ✗ retVal = nlsData->residualFuncConstraints(&resUserData, x, f, &iflag);
974 ✗ solverData->numberOfFunctionEvaluations++;
975
976 ✗ return retVal;
977 }
978
979 /*! \fn wrapper_fvec_der for the residual Function
980 * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data)
981 *
982 * \author bbachmann
983 *
984 */
985 ✗ static int wrapper_fvec_der(DATA_HOMOTOPY* solverData, double* x, double* fJac)
986 {
987 ✗ NONLINEAR_SYSTEM_DATA* nlsData = solverData->userData->nlsData;
988
989 /* performance measurement */
990 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
991
992 /* calculate jacobian */
993 ✗ if(nlsData->jacobianIndex != -1)
994 {
995 /* !!!!!!!!!!! Be sure that actual x is used !!!!!!!!!!! */
996 ✗ getAnalyticalJacobianHomotopy(solverData, fJac);
997 }
998 else
999 {
1000 ✗ getNumericalJacobianHomotopy(solverData, x, fJac);
1001 }
1002
1003 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC_TEST))
1004 {
1005 ✗ int n = solverData->n;
1006 /* debugMatrixDouble(OMC_LOG_NLS_JAC_TEST,"analytical jacobian:",fJac, n, n+1); */
1007 ✗ getNumericalJacobianHomotopy(solverData, x, solverData->debug_fJac);
1008 /* debugMatrixDouble(OMC_LOG_NLS_JAC_TEST,"numerical jacobian:",solverData->debug_fJac, n, n+1); */
1009 ✗ matDiffBB(n, fJac, solverData->debug_fJac, solverData->debug_fJac);
1010 /* debugMatrixDouble(OMC_LOG_NLS_JAC_TEST,"Difference of jacobians:",solverData->debug_fJac, n, n+1); */
1011 ✗ debugDouble(OMC_LOG_NLS_JAC_TEST,"error between analytical and numerical jacobian = ", vecMaxNorm(n*n, solverData->debug_fJac));
1012 ✗ vecDivScaling(n*(n+1), solverData->debug_fJac , fJac, solverData->debug_fJac);
1013 ✗ debugDouble(OMC_LOG_NLS_JAC_TEST,"relative error between analytical and numerical jacobian = ", vecMaxNorm(n*n, solverData->debug_fJac));
1014 ✗ messageClose(OMC_LOG_NLS_JAC_TEST); // FIXME what does this belong to?
1015 }
1016 /* performance measurement and statistics */
1017 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
1018 ✗ nlsData->numberOfJEval++;
1019
1020 ✗ return 0;
1021 }
1022
1023 /*! \fn wrapper_fvec_homotopy_newton for the residual Function
1024 *
1025 * \author bbachmann
1026 *
1027 */
1028 ✗ static int wrapper_fvec_homotopy_newton(DATA_HOMOTOPY* solverData, double* x, double* h)
1029 {
1030 int i;
1031 ✗ int n = solverData->n;
1032
1033 /* Newton homotopy */
1034 ✗ wrapper_fvec(solverData, x, solverData->f1);
1035 ✗ vecAddScal(solverData->n, solverData->f1, solverData->fx0, - (1-x[n]), h);
1036
1037 ✗ return 0;
1038 }
1039
1040 /*! \fn wrapper_fvec_homotopy_newton_der for the residual Function
1041 *
1042 * \author bbachmann
1043 *
1044 */
1045 ✗ static int wrapper_fvec_homotopy_newton_der(DATA_HOMOTOPY* solverData, double* x, double* hJac)
1046 {
1047 int i, j;
1048 ✗ int n = solverData->n;
1049
1050 /* Newton homotopy */
1051 ✗ wrapper_fvec_der(solverData, x, hJac);
1052
1053 /* add f(x0) as the last column of the Jacobian*/
1054 ✗ vecCopy(n, solverData->fx0, hJac + n*n);
1055
1056 ✗ return 0;
1057 }
1058
1059 /*! \fn wrapper_fvec_homotopy_fixpoint for the residual Function
1060 *
1061 * \author bbachmann
1062 *
1063 */
1064 ✗ static int wrapper_fvec_homotopy_fixpoint(DATA_HOMOTOPY* solverData, double* x, double* h)
1065 {
1066 int i;
1067 ✗ int n = solverData->n;
1068
1069 /* Fixpoint homotopy */
1070 ✗ wrapper_fvec(solverData, x, solverData->f1);
1071 ✗ for (i=0; i<n; i++){
1072 ✗ h[i] = x[n]*solverData->f1[i] + (1-x[n]) * (x[i]-solverData->x0[i]);
1073 }
1074
1075 ✗ return 0;
1076 }
1077
1078 /*! \fn wrapper_fvec_homotopy_fixpoint_der for the residual Function
1079 *
1080 * \author bbachmann
1081 *
1082 */
1083 ✗ static int wrapper_fvec_homotopy_fixpoint_der(DATA_HOMOTOPY* solverData, double* x, double* hJac)
1084 {
1085 int i, j;
1086 ✗ int n = solverData->n;
1087
1088 /* Fixpoint homotopy */
1089 ✗ wrapper_fvec_der(solverData, x, hJac);
1090 ✗ for (i=0; i<n; i++){
1091 ✗ for (j=0; j<n; j++) {
1092 ✗ hJac[i+ j * n] = x[n]*hJac[i+ j * n];
1093 }
1094 ✗ hJac[i+ i * n] = hJac[i+ i * n] + (1-x[n]);
1095 ✗ hJac[i+ n * n] = solverData->f1[i]-(x[i] - solverData->x0[i]);
1096 }
1097 ✗ return 0;
1098 }
1099
1100 /*! \fn getIndicesOfPivotElement for calculating pivot element
1101 *
1102 * \author bbachmann
1103 *
1104 */
1105 ✗ void getIndicesOfPivotElement(int *n, int *m, int *l, double* A, int *indRow, int *indCol, int *pRow, int *pCol, double *absMax)
1106 {
1107 int i, j;
1108
1109 ✗ *absMax = fabs(A[indRow[*l] + indCol[*l]* *n]);
1110 ✗ *pCol = *l;
1111 ✗ *pRow = *l;
1112 ✗ for (i = *l; i < *n; i++) {
1113 ✗ for (j = *l; j < *m; j++) {
1114 ✗ if (fabs(A[indRow[i] + indCol[j]* *n]) > *absMax) {
1115 ✗ *absMax = fabs(A[indRow[i] + indCol[j]* *n]);
1116 ✗ *pCol = j;
1117 ✗ *pRow = i;
1118 }
1119 }
1120 }
1121 ✗ }
1122
1123
1124 /*! \fn solveSystemWithTotalPivotSearch for solution of overdetermined linear system
1125 * used for the homotopy solver, for calculating the direction
1126 * used for the newton solver, for calculating the Newton step
1127 *
1128 * \author bbachmann
1129 *
1130 */
1131 ✗ int solveSystemWithTotalPivotSearch(DATA *data, int n, double* x, double* A, int* indRow, int* indCol, int *pos, int *rank, int casualTearingSet)
1132 {
1133 ✗ int i, k, j, m=n+1, nPivot=n;
1134 int pCol, pRow;
1135 double hValue;
1136 double hInt;
1137 double absMax, detJac;
1138 int returnValue = 0;
1139
1140 ✗ debugMatrixDouble(OMC_LOG_NLS_JAC,"Linear System Matrix [Jac res]:",A, n, m);
1141 ✗ debugVectorDouble(OMC_LOG_NLS_JAC,"vector b:", A+n*n, n);
1142
1143 /* assume full rank of matrix [n x (n+1)] */
1144 ✗ *rank = n;
1145
1146 ✗ for (i=0; i<n; i++) {
1147 ✗ indRow[i] = i;
1148 }
1149 ✗ for (i=0; i<m; i++) {
1150 ✗ indCol[i] = i;
1151 }
1152 ✗ if (*pos>=0) {
1153 ✗ indCol[n] = *pos;
1154 ✗ indCol[*pos] = n;
1155 } else {
1156 ✗ nPivot = n+1;
1157 }
1158
1159 ✗ for (i = 0; i < n; i++) {
1160 ✗ getIndicesOfPivotElement(&n, &nPivot, &i, A, indRow, indCol, &pRow, &pCol, &absMax);
1161 ✗ if (absMax<DBL_EPSILON) {
1162 ✗ *rank = i;
1163 ✗ if (data->simulationInfo->initial) {
1164 ✗ warningStreamPrint(OMC_LOG_NLS_V, 1, "Homotopy solver total pivot: Matrix (nearly) singular at initialization.");
1165 } else {
1166 ✗ warningStreamPrint(OMC_LOG_NLS_V, 1, "Homotopy solver total pivot: Matrix (nearly) singular at time %f.", data->localData[0]->timeValue);
1167 }
1168 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Continuing anyway. For more information please use -lv %s.", OMC_LOG_STREAM_NAME[OMC_LOG_NLS_V]);
1169 ✗ messageCloseWarning(OMC_LOG_NLS_V);
1170 ✗ debugInt(OMC_LOG_NLS_V,"rank = ", *rank);
1171 ✗ debugInt(OMC_LOG_NLS_V,"position = ", *pos);
1172 break;
1173 }
1174 /* swap row indices */
1175 ✗ if (pRow!=i) {
1176 ✗ hInt = indRow[i];
1177 ✗ indRow[i] = indRow[pRow];
1178 ✗ indRow[pRow] = hInt;
1179 }
1180 /* swap column indices */
1181 ✗ if (pCol!=i) {
1182 ✗ hInt = indCol[i];
1183 ✗ indCol[i] = indCol[pCol];
1184 ✗ indCol[pCol] = hInt;
1185 }
1186
1187 /* Gauss elimination of row indRow[i] */
1188 ✗ for (k=i+1; k<n; k++) {
1189 ✗ hValue = -A[indRow[k] + indCol[i]*n]/A[indRow[i] + indCol[i]*n];
1190 ✗ for (j=i+1; j<m; j++) {
1191 ✗ A[indRow[k] + indCol[j]*n] = A[indRow[k] + indCol[j]*n] + hValue*A[indRow[i] + indCol[j]*n];
1192 }
1193 ✗ A[indRow[k] + indCol[i]*n] = 0;
1194 }
1195 }
1196
1197 ✗ for (detJac=1.0,k=0; k<n; k++) detJac *= A[indRow[k] + indCol[k]*n];
1198
1199 ✗ debugMatrixPermutedDouble(OMC_LOG_NLS_JAC,"Linear System Matrix [Jac res] after decomposition",A, n, m, indRow, indCol);
1200 debugDouble(OMC_LOG_NLS_JAC,"Determinant = ", detJac);
1201 ✗ if (isnan(detJac)){
1202 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Jacobian determinant is NaN.");
1203 ✗ return -1;
1204 }
1205 ✗ else if (fabs(detJac) < 1e-9 && casualTearingSet)
1206 {
1207 debugString(OMC_LOG_DT, "The determinant of the casual tearing set is vanishing, let's fail if this is not the solution...");
1208 returnValue = 1;
1209 }
1210
1211 /* Solve even singular matrices !!! */
1212 ✗ for (i=n-1;i>=0; i--) {
1213 ✗ if (i>=*rank) {
1214 /* this criteria should be evaluated and may be improved in future */
1215 ✗ if (fabs(A[indRow[i] + indCol[n]*n])>1e-6) {
1216 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "under-determined linear system not solvable!");
1217 ✗ return -1;
1218 } else {
1219 ✗ x[indCol[i]] = 0.0;
1220 }
1221 } else {
1222 ✗ x[indCol[i]] = -A[indRow[i] + indCol[n]*n];
1223 ✗ for (j=n-1; j>i; j--) {
1224 ✗ x[indCol[i]] = x[indCol[i]] - A[indRow[i] + indCol[j]*n]*x[indCol[j]];
1225 }
1226 ✗ x[indCol[i]]=x[indCol[i]]/A[indRow[i] + indCol[i]*n];
1227 }
1228 }
1229 ✗ x[indCol[n]]=1.0;
1230 ✗ debugVectorInt(OMC_LOG_NLS_V,"indRow:", indRow, n);
1231 ✗ debugVectorInt(OMC_LOG_NLS_V,"indCol:", indCol, n+1);
1232 ✗ debugVectorDouble(OMC_LOG_NLS_V,"vector x (solution):", x, n+1);
1233
1234 /* Return position of largest value (1.0) */
1235 ✗ if (*pos<0) {
1236 ✗ *pos=indCol[n];
1237 debugInt(OMC_LOG_NLS_V,"position of largest value = ", *pos);
1238 }
1239
1240 return returnValue;
1241 }
1242
1243
1244 /*! \fn linearSolverWrapper
1245 */
1246 ✗ int linearSolverWrapper(DATA *data, int n, double* x, double* A, int* indRow, int* indCol, int *pos, int *rank, int method, int casualTearingSet)
1247 {
1248 /* First try to use lapack and if it fails then
1249 * use solveSystemWithTotalPivotSearch */
1250 int returnValue = -1;
1251 int solverinfo;
1252 ✗ int nrhs = 1;
1253 ✗ int lda = n;
1254 int k;
1255 double detJac;
1256
1257 ✗ debugMatrixDouble(OMC_LOG_NLS_JAC,"Linear System Matrix [Jac res]:", A, n, n+1);
1258 ✗ debugVectorDouble(OMC_LOG_NLS_JAC,"vector b:", x, n);
1259
1260 ✗ switch(method){
1261 ✗ case NLS_LS_TOTALPIVOT:
1262
1263 ✗ solverinfo = solveSystemWithTotalPivotSearch(data, n, x, A, indRow, indCol, pos, rank, casualTearingSet);
1264 /* in case of failing */
1265 ✗ if (solverinfo == -1)
1266 {
1267 /* debug information */
1268 debugString(OMC_LOG_NLS_V, "Linear total pivot solver failed!!!");
1269 debugString(OMC_LOG_NLS_V, "******************************************************");
1270 }
1271 ✗ else if (solverinfo == 1)
1272 {
1273 returnValue = 1;
1274 }
1275 else
1276 {
1277 returnValue = 0;
1278 }
1279 break;
1280 ✗ case NLS_LS_LAPACK:
1281 /* Solve system with lapack */
1282 ✗ dgesv_((int*) &n,
1283 (int*) &nrhs,
1284 A,
1285 (int*) &lda,
1286 indRow,
1287 x,
1288 (int*) &n,
1289 &solverinfo);
1290
1291 ✗ for (detJac=1.0, k=0; k<n; k++) detJac *= A[k + k*n];
1292
1293 ✗ debugMatrixDouble(OMC_LOG_NLS_JAC,"Linear system matrix [Jac res] after decomposition:", A, n, n+1);
1294 debugDouble(OMC_LOG_NLS_JAC,"Determinant = ", detJac);
1295
1296 /* in case of failing */
1297 ✗ if (solverinfo != 0)
1298 {
1299 /* debug information */
1300 debugString(OMC_LOG_NLS_V, "Linear lapack solver failed!!!");
1301 debugString(OMC_LOG_NLS_V, "******************************************************");
1302 }
1303 ✗ else if (fabs(detJac) < 1e-9 && casualTearingSet)
1304 {
1305 debugString(OMC_LOG_DT, "The determinant of the casual tearing set is vanishing, let's fail if this is not the solution...");
1306 ✗ returnValue = 1;
1307 }
1308 else
1309 {
1310 ✗ vecScalarMult(n, x, -1, x);
1311 returnValue = 0;
1312 }
1313 break;
1314 ✗ default:
1315 ✗ throwStreamPrint(0, "Non-Linear solver try to run with a unknown linear solver (%d).", method);
1316 }
1317
1318 /* Debugging error of linear system */
1319 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC))
1320 {
1321 ✗ double* res = (double*) calloc(n,sizeof(double));
1322 ✗ debugVectorDouble(OMC_LOG_NLS_JAC,"solution:", x, n);
1323 ✗ matVecMult(n, n, A, x, res);
1324 ✗ debugVectorDouble(OMC_LOG_NLS_JAC,"test solution:", res, n);
1325 ✗ debugDouble(OMC_LOG_NLS_JAC,"error of linear system = ", vec2Norm(n, res));
1326 ✗ free(res);
1327 ✗ messageClose(OMC_LOG_NLS_JAC); // FIXME what does this belong to?
1328 }
1329
1330 ✗ return returnValue;
1331 }
1332
1333
1334 /* Pushes e onto hist; true if e moved away from the last residual and back to
1335 * one from 2-4 iterations ago (a limit cycle). */
1336 ✗ static int newtonLimitCycleStep(double *hist, int *nHist, double e)
1337 {
1338 int k, cycle = 0;
1339 ✗ if (*nHist >= 4 && fabs(e - hist[0]) > 1e-2 * e) {
1340 ✗ for (k = 1; k < 4; k++) {
1341 ✗ cycle |= fabs(e - hist[k]) <= 1e-2 * e;
1342 }
1343 }
1344 ✗ for (k = 3; k > 0; k--) {
1345 ✗ hist[k] = hist[k-1];
1346 }
1347 ✗ hist[0] = e;
1348 ✗ (*nHist)++;
1349 ✗ return cycle;
1350 }
1351
1352 /*! \fn solve system with damped Newton-Raphson
1353 *
1354 * \author bbachmann
1355 *
1356 */
1357 ✗ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
1358 {
1359 ✗ int numberOfIterations = 0 ,i, j, n=solverData->n, m=solverData->m;
1360 ✗ int pos = solverData->n, rank;
1361 double error_f_sqrd, error_f1_sqrd, error_f2_sqrd, error_f_sqrd_scaled, error_f1_sqrd_scaled;
1362 double delta_x_sqrd, delta_x_sqrd_scaled, grad_f, grad_f_scaled;
1363 ✗ int numberOfSmallSteps = 0, smallStepsAtHover = 0, lessAccurate, atCycleBottom;
1364 double error_f_old = 1e100, error_f_old_scaled = 1e100;
1365 ✗ int countNegativeSteps = 0;
1366 ✗ int countCycles = 0, nHist = 0, countStalls = 0, countHovers = 0, countCreeps = 0;
1367 ✗ double error_f_best = 1e100, error_f_hover = 1e100, stepLambda;
1368 double errorHist[4];
1369 double lambda;
1370 double lambda1, lambda2;
1371 double lambdaMin = 1e-4;
1372 double a2, a3, rhs1, rhs2, D;
1373 double alpha = 1e-1;
1374 int firstrun;
1375 int constraintViolated;
1376 ✗ int solverinfo = 0;
1377 int lastWasGood = 0; /* boolean, keeps track of previous x */
1378
1379 ✗ int assert = 1;
1380 ✗ DATA* data = solverData->userData->data;
1381 ✗ threadData_t *threadData = solverData->userData->threadData;
1382 ✗ NONLINEAR_SYSTEM_DATA* nlsData = solverData->userData->nlsData;
1383 ✗ int linearSolverMethod = data->simulationInfo->nlsLinearSolver;
1384
1385 /* debug information */
1386 debugString(OMC_LOG_NLS_V, "******************************************************");
1387 ✗ debugInt(OMC_LOG_NLS_V, "NEWTON SOLVER STARTED! equation number: ",solverData->eqSystemNumber);
1388 ✗ debugInt(OMC_LOG_NLS_V, "maximum number of function evaluation: ", solverData->maxNumberOfIterations);
1389 ✗ printUnknowns(OMC_LOG_NLS_V, solverData);
1390
1391 /* set default solver message */
1392 ✗ solverData->info = 0;
1393
1394 /* calculated error of function values */
1395 ✗ error_f_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1396 ✗ vecDivScaling(solverData->n, solverData->f1, solverData->resScaling, solverData->fvecScaled);
1397 ✗ error_f_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->fvecScaled);
1398
1399 while(1)
1400 {
1401 ✗ numberOfIterations++;
1402 /* debug information */
1403 debugInt(OMC_LOG_NLS_V, "Iteration:", numberOfIterations);
1404
1405 /* solve jacobian and function value (both stored in hJac, last column is fvec), side effects: jacobian matrix is changed */
1406 ✗ if (numberOfIterations>1)
1407 ✗ solverinfo = linearSolverWrapper(data, solverData->n, solverData->dy0, solverData->fJac, solverData->indRow, solverData->indCol, &pos, &rank, linearSolverMethod, solverData->casualTearingSet);
1408
1409 ✗ if (solverinfo == -1)
1410 {
1411 /* report solver abortion */
1412 ✗ solverData->info=-1;
1413 /* debug information */
1414 debugString(OMC_LOG_NLS_V, "NEWTON SOLVER DID ---NOT--- CONVERGE TO A SOLUTION!!!");
1415 debugString(OMC_LOG_NLS_V, "******************************************************");
1416 assert = 0;
1417 break;
1418 }
1419 else
1420 {
1421 /* Scaling back to original variables */
1422 ✗ vecMultScaling(solverData->m, solverData->dy0, solverData->xScaling, solverData->dy0);
1423 /* try full Newton step */
1424 ✗ vecAdd(solverData->n, x, solverData->dy0, solverData->x1);
1425 ✗ printNewtonStep(OMC_LOG_NLS_V, solverData);
1426
1427 /* Damping strategy, performance is very sensitive on the value of lambda */
1428 lambda1 = 1.0;
1429 assert = 1;
1430 firstrun = 1;
1431 ✗ while (assert && (lambda1 > lambdaMin))
1432 {
1433 ✗ if (!firstrun){
1434 ✗ lambda1 *= 0.655;
1435 ✗ vecAddScal(solverData->n, x, solverData->dy0, lambda1, solverData->x1);
1436 assert = 1;
1437 }
1438 #ifndef OMC_EMCC
1439 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1440 #endif
1441 ✗ if (solverData->casualTearingSet){
1442 ✗ constraintViolated = solverData->f_con(solverData, solverData->x1, solverData->f1);
1443 ✗ if (constraintViolated){
1444 lambda1 = lambdaMin-1;
1445 ✗ break;
1446 }
1447 }
1448 else
1449 ✗ solverData->f(solverData, solverData->x1, solverData->f1);
1450
1451 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
1452 #ifndef OMC_EMCC
1453 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1454 #endif
1455 firstrun = 0;
1456 ✗ if (assert) {
1457 debugDouble(OMC_LOG_NLS_V, "Assert of Newton step: lambda1 =", lambda1);
1458 }
1459 }
1460
1461 ✗ if (lambda1 < lambdaMin)
1462 {
1463 ✗ debugDouble(OMC_LOG_NLS_V, "UPS! MUST HANDLE A PROBLEM (Newton method), time : ", solverData->timeValue);
1464 ✗ solverData->info = -1;
1465 ✗ break;
1466 }
1467
1468 /* Damping (see Numerical Recipes) */
1469 /* calculate gradient of quadratic function for damping strategy */
1470 ✗ grad_f = -2.0*error_f_sqrd;
1471 ✗ grad_f_scaled = -2.0*error_f_sqrd_scaled;
1472 ✗ error_f1_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1473 ✗ vecDivScaling(solverData->n, solverData->f1, solverData->resScaling, solverData->f2);
1474 ✗ error_f1_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->f2);
1475 debugDouble(OMC_LOG_NLS_V, "Need to damp, grad_f = ", grad_f);
1476 ✗ debugDouble(OMC_LOG_NLS_V, "Need to damp, error_f = ", sqrt(error_f_sqrd));
1477 debugDouble(OMC_LOG_NLS_V, "Need to damp this!! lambda1 = ", lambda1);
1478 ✗ debugDouble(OMC_LOG_NLS_V, "Need to damp, error_f1 = ", sqrt(error_f1_sqrd));
1479 ✗ debugDouble(OMC_LOG_NLS_V, "Need to damp, forced error = ", error_f_sqrd + alpha*lambda1*grad_f);
1480 ✗ stepLambda = lambda1;
1481 ✗ if ((error_f1_sqrd > error_f_sqrd + alpha*lambda1*grad_f)
1482 ✗ && (error_f1_sqrd_scaled > error_f_sqrd_scaled + alpha*lambda1*grad_f_scaled)
1483 ✗ && (error_f_sqrd > 1e-12) && (error_f_sqrd_scaled > 1e-12))
1484 {
1485 ✗ lambda2 = fmax(-lambda1*lambda1*grad_f/(2*(error_f1_sqrd-error_f_sqrd-lambda1*grad_f)),lambdaMin);
1486 ✗ stepLambda = lambda2;
1487 debugDouble(OMC_LOG_NLS_V, "Need to damp this!! lambda2 = ", lambda2);
1488 ✗ vecAddScal(solverData->n, x, solverData->dy0, lambda2, solverData->x1);
1489 assert= 1;
1490 #ifndef OMC_EMCC
1491 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1492 #endif
1493 ✗ if (solverData->casualTearingSet){
1494 ✗ constraintViolated = solverData->f_con(solverData, solverData->x1, solverData->f1);
1495 ✗ if (constraintViolated){
1496 ✗ solverData->info = -1;
1497 ✗ break;
1498 }
1499 }
1500 else
1501 ✗ solverData->f(solverData, solverData->x1, solverData->f1);
1502
1503 ✗ error_f2_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1504 ✗ debugDouble(OMC_LOG_NLS_V, "Need to damp, error_f2 = ", sqrt(error_f2_sqrd));
1505 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
1506 #ifndef OMC_EMCC
1507 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1508 #endif
1509 ✗ if (assert)
1510 {
1511 ✗ debugDouble(OMC_LOG_NLS_V, "UPS! MUST HANDLE A PROBLEM (Newton method), time : ", solverData->timeValue);
1512 ✗ solverData->info = -1;
1513 ✗ break;
1514 }
1515 ✗ if ((error_f1_sqrd > error_f_sqrd + alpha*lambda2*grad_f) && (error_f_sqrd > 1e-12) && (error_f_sqrd_scaled > 1e-12))
1516 {
1517 ✗ rhs1 = error_f1_sqrd - grad_f*lambda1 - error_f_sqrd;
1518 ✗ rhs2 = error_f2_sqrd - grad_f*lambda2 - error_f_sqrd;
1519 ✗ a3 = (rhs1/(lambda1*lambda1) - rhs2/(lambda2*lambda2))/(lambda1 - lambda2);
1520 ✗ a2 = (-lambda2*rhs1/(lambda1*lambda1) + lambda1*rhs2/(lambda2*lambda2))/(lambda1 - lambda2);
1521 ✗ if (a3==0.0)
1522 ✗ lambda = -grad_f/(2.0*a2);
1523 else
1524 {
1525 ✗ D = a2*a2 - 3.0*a3*grad_f;
1526 ✗ if (D <= 0.0)
1527 ✗ lambda = 0.5*lambda1;
1528 else
1529 ✗ if (a2 <= 0.0)
1530 ✗ lambda = (-a2+sqrt(D))/(3.0*a3);
1531 else
1532 ✗ lambda = -grad_f/(a2+sqrt(D));
1533 }
1534 ✗ lambda = fmax(lambda, lambdaMin);
1535 ✗ stepLambda = lambda;
1536 debugDouble(OMC_LOG_NLS_V, "Need to damp this!! lambda = ", lambda);
1537 ✗ vecAddScal(solverData->n, x, solverData->dy0, lambda, solverData->x1);
1538 assert= 1;
1539 #ifndef OMC_EMCC
1540 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1541 #endif
1542 ✗ if (solverData->casualTearingSet){
1543 ✗ constraintViolated = solverData->f_con(solverData, solverData->x1, solverData->f1);
1544 ✗ if (constraintViolated){
1545 ✗ solverData->info = -1;
1546 ✗ break;
1547 }
1548 }
1549 else
1550 ✗ solverData->f(solverData, solverData->x1, solverData->f1);
1551
1552 ✗ error_f1_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1553 ✗ debugDouble(OMC_LOG_NLS_V, "Need to damp, error_f1 = ", sqrt(error_f1_sqrd));
1554 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
1555 #ifndef OMC_EMCC
1556 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1557 #endif
1558 ✗ if (assert)
1559 {
1560 ✗ debugDouble(OMC_LOG_NLS_V, "UPS! MUST HANDLE A PROBLEM (Newton method), time : ", solverData->timeValue);
1561 ✗ solverData->info = -1;
1562 ✗ break;
1563 }
1564 }
1565 }else{
1566 ✗ lambda = lambda1;
1567 }
1568 }
1569
1570 /* Calculate different error measurements */
1571 ✗ vecDivScaling(solverData->n, solverData->f1, solverData->resScaling, solverData->fvecScaled);
1572 ✗ debugVectorDouble(OMC_LOG_NLS_V, "function values:",solverData->f1, n);
1573 ✗ debugVectorDouble(OMC_LOG_NLS_V, "scaled function values:",solverData->fvecScaled, n);
1574
1575 /* update delta_x_sqrd, error_f_sqrd */
1576 ✗ vecDivScaling(solverData->n, solverData->dy0, solverData->xScaling, solverData->dxScaled);
1577 ✗ delta_x_sqrd = vec2NormSqrd(solverData->n, solverData->dy0);
1578 ✗ delta_x_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->dxScaled);
1579
1580 ✗ error_f_old = error_f_sqrd;
1581 ✗ error_f_old_scaled = error_f_sqrd_scaled;
1582 ✗ error_f_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1583 ✗ error_f_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->fvecScaled);
1584
1585 ✗ countNegativeSteps += (error_f_sqrd > 10*error_f_old);
1586 /* a cycle within the less accuracy band is left to the other exits */
1587 ✗ countCycles = newtonLimitCycleStep(errorHist, &nHist, error_f_sqrd)
1588 ✗ && error_f_sqrd >= solverData->ftol_sqrd*1e6 && error_f_sqrd_scaled >= solverData->ftol_sqrd*1e6 ? countCycles + 1 : 0;
1589 ✗ if (error_f_sqrd < 0.99*error_f_best) {
1590 error_f_best = error_f_sqrd;
1591 countStalls = 0;
1592 countHovers = 0;
1593 ✗ } else if (error_f_sqrd < 10*error_f_hover) {
1594 ✗ countStalls++;
1595 ✗ countHovers++;
1596 } else {
1597 ✗ countStalls++;
1598 countHovers = 0;
1599 }
1600 ✗ if (countHovers == 0 || error_f_sqrd < error_f_hover) {
1601 error_f_hover = error_f_sqrd;
1602 }
1603 ✗ countCreeps = stepLambda < 1e-3 ? countCreeps + 1 : 0;
1604 ✗ lastWasGood = error_f_sqrd >= error_f_old;
1605
1606
1607 /* debug information */
1608 ✗ if (omc_useStream[OMC_LOG_NLS_V]) {
1609 debugString(OMC_LOG_NLS_V, "error measurements:");
1610 ✗ debugDouble(OMC_LOG_NLS_V, "delta_x =", sqrt(delta_x_sqrd));
1611 ✗ debugDouble(OMC_LOG_NLS_V, "delta_x_scaled =", sqrt(delta_x_sqrd_scaled));
1612 ✗ debugDouble(OMC_LOG_NLS_V, "newtonXTol =", sqrt(solverData->xtol_sqrd));
1613 ✗ debugDouble(OMC_LOG_NLS_V, "error_f =", sqrt(error_f_sqrd));
1614 ✗ debugDouble(OMC_LOG_NLS_V, "error_f_scaled =", sqrt(error_f_sqrd_scaled));
1615 ✗ debugDouble(OMC_LOG_NLS_V, "newtonFTol =", sqrt(solverData->ftol_sqrd));
1616 }
1617
1618 #if !defined(OMC_MINIMAL_RUNTIME)
1619 ✗ if (data->simulationInfo->nlsCsvInfomation){
1620 ✗ print_csvLineIterStats(((struct csvStats*) nlsData->csvData)->iterStats,
1621 ✗ nlsData->size,
1622 ✗ nlsData->numberOfCall+1,
1623 numberOfIterations,
1624 solverData->x,
1625 solverData->f1,
1626 delta_x_sqrd,
1627 delta_x_sqrd_scaled,
1628 error_f_sqrd,
1629 error_f_sqrd_scaled,
1630 lambda
1631 );
1632 }
1633 #endif
1634 /* away from any solution: no 1% improvement for long, or only heavily damped steps */
1635 ✗ if (countNegativeSteps > 20 || countCycles > 20 || ((countStalls > 400 || countCreeps > 400) && error_f_sqrd >= solverData->ftol_sqrd*1e6 && error_f_sqrd_scaled >= solverData->ftol_sqrd*1e6))
1636 {
1637 debugInt(OMC_LOG_NLS_V, "UPS! Something happened, NegativeSteps = ", countNegativeSteps);
1638 ✗ solverData->info = -1;
1639 ✗ break;
1640 }
1641
1642 /* solution found */
1643 ✗ if (((error_f_sqrd < solverData->ftol_sqrd) || (error_f_sqrd_scaled < solverData->ftol_sqrd)) && ((delta_x_sqrd_scaled < solverData->xtol_sqrd) || (delta_x_sqrd < solverData->xtol_sqrd)))
1644 {
1645 ✗ solverData->info = 1;
1646
1647 /* reject new x if old x is as good, for stability (see issue #6419) */
1648 ✗ if (lastWasGood)
1649 {
1650 debugString(OMC_LOG_NLS_V, "Note: newton solver rejected last x because previous was as good");
1651 }
1652 else
1653 {
1654 ✗ vecCopy(solverData->n, solverData->x1, x);
1655 }
1656
1657 /* update statistics */
1658 ✗ solverData->numberOfIterations += numberOfIterations;
1659 ✗ solverData->error_f_sqrd = error_f_sqrd;
1660
1661 ✗ break;
1662 }
1663 ✗ else if (solverinfo == 1){
1664 ✗ solverData->info = -1;
1665 debugString(OMC_LOG_DT, "It is not the solution.");
1666 break;
1667 }
1668 /* the residual stopped decreasing at an x that meets the tolerance: further steps are round-off */
1669 ✗ else if (lastWasGood && ((error_f_old < solverData->ftol_sqrd) || (error_f_old_scaled < solverData->ftol_sqrd)))
1670 {
1671 ✗ solverData->info = 1;
1672 debugString(OMC_LOG_NLS_V, "Note: newton solver rejected last x because previous was as good");
1673 ✗ solverData->numberOfIterations += numberOfIterations;
1674 ✗ solverData->error_f_sqrd = error_f_old;
1675 ✗ break;
1676 }
1677
1678 /* check if maximum iteration is reached */
1679 ✗ if (numberOfIterations > solverData->maxNumberOfIterations)
1680 {
1681 ✗ solverData->info = -1;
1682 ✗ if (data->simulationInfo->initial) {
1683 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Homotopy solver Newton iteration: Maximum number of iterations reached at initialization, but no root found.");
1684 } else {
1685 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Homotopy solver Newton iteration: Maximum number of iterations reached at time %f, but no root found.", data->localData[0]->timeValue);
1686 }
1687 /* debug information */
1688 debugString(OMC_LOG_NLS_V, "NEWTON SOLVER DID ---NOT--- CONVERGE TO A SOLUTION!!!");
1689 debugString(OMC_LOG_NLS_V, "******************************************************");
1690
1691 /* update statistics */
1692 ✗ solverData->numberOfIterations += numberOfIterations;
1693 ✗ break;
1694 }
1695
1696 ✗ numberOfSmallSteps += (delta_x_sqrd < solverData->xtol_sqrd*1e4) || (delta_x_sqrd_scaled < solverData->xtol_sqrd*1e4);
1697 ✗ if (countHovers == 0) {
1698 ✗ smallStepsAtHover = numberOfSmallSteps;
1699 }
1700 /* the bottom of a stationary cycle without small steps, which are left to their own exit */
1701 ✗ atCycleBottom = countHovers > 20 && numberOfSmallSteps == smallStepsAtHover
1702 ✗ && error_f_sqrd <= errorHist[1] && error_f_sqrd <= errorHist[2] && error_f_sqrd <= errorHist[3];
1703 /* check changes in unknown vector */
1704 ✗ lessAccurate = (error_f_sqrd < solverData->ftol_sqrd*1e6) || (error_f_sqrd_scaled < solverData->ftol_sqrd*1e6);
1705 /* a stationary residual within the less accuracy band is round-off, like small steps */
1706 ✗ if ((delta_x_sqrd < solverData->xtol_sqrd) || (delta_x_sqrd_scaled < solverData->xtol_sqrd) || (numberOfSmallSteps > 20) || (lessAccurate && atCycleBottom))
1707 {
1708 ✗ if (lessAccurate)
1709 {
1710 ✗ solverData->info = 1;
1711 ✗ if (atCycleBottom) {
1712 ✗ vecCopy(solverData->n, solverData->x1, x);
1713 }
1714
1715 /* debug information */
1716 debugString(OMC_LOG_NLS_V, "NEWTON SOLVER DID CONVERGE TO A SOLUTION WITH LESS ACCURACY!!!");
1717 ✗ printUnknowns(OMC_LOG_NLS_V, solverData);
1718 debugString(OMC_LOG_NLS_V, "******************************************************");
1719 ✗ solverData->error_f_sqrd = 0;
1720
1721 } else
1722 {
1723 ✗ solverData->info = -1;
1724 debugString(OMC_LOG_NLS_V, "Warning: newton solver gets stuck!!!");
1725 /* debug information */
1726 debugString(OMC_LOG_NLS_V, "NEWTON SOLVER DID ---NOT--- CONVERGE TO A SOLUTION!!!");
1727 debugString(OMC_LOG_NLS_V, "******************************************************");
1728 }
1729 /* update statistics */
1730 ✗ solverData->numberOfIterations += numberOfIterations;
1731 ✗ break;
1732 }
1733 assert = 1;
1734 #ifndef OMC_EMCC
1735 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1736 #endif
1737 /* updating x */
1738 ✗ vecCopy(solverData->n, solverData->x1, x);
1739
1740 /* calculate jacobian and function values (both stored in fJac, last column is fvec) */
1741 ✗ solverData->fJac_f(solverData, x, solverData->fJac);
1742 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
1743 #ifndef OMC_EMCC
1744 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1745 #endif
1746 ✗ if (assert)
1747 {
1748 /* report solver abortion */
1749 ✗ solverData->info=-1;
1750 debugString(OMC_LOG_NLS_V,"UPS! assert when calculating Jacobian!!!");
1751 break;
1752 }
1753 ✗ vecCopy(n, solverData->f1, solverData->fJac + n*n);
1754 /* calculate scaling factor of residuals */
1755 ✗ matVecMultAbsBB(solverData->n, solverData->fJac, solverData->ones, solverData->resScaling);
1756 ✗ debugVectorDouble(OMC_LOG_NLS_JAC, "residuum scaling:", solverData->resScaling, solverData->n);
1757 ✗ scaleMatrixRows(solverData->n, solverData->m, solverData->fJac);
1758 ✗ vecCopy(n, solverData->fJac + n*n, solverData->dy0);
1759 }
1760 ✗ return 0;
1761 }
1762
1763 /*! \fn solve system with homotopy method
1764 *
1765 * \author bbachmann
1766 */
1767 ✗ static int homotopyAlgorithm(DATA_HOMOTOPY* solverData, double *x)
1768 {
1769 int i, j;
1770 double error_h, error_h_scaled, delta_x;
1771 double vecScalarProduct;
1772
1773 int pos, rank;
1774 ✗ int iter = 0;
1775 ✗ int maxiter = homMaxNewtonSteps;
1776 ✗ int maxTries = homMaxTries;
1777 ✗ int numSteps = 0;
1778 ✗ int stepAccept = 0;
1779 ✗ int correctorStrategy = homBacktraceStrategy; /* 1: go back to the path by fixing one coordinate, 2: go back to the path in an orthogonal direction to the tangent vector */
1780 ✗ double bend = 0;
1781 ✗ double tau = homTauStart, tauMax = homTauMax, tauMin = homTauMin, hEps = homHEps, adaptBend = homAdaptBend;
1782 ✗ double tauDecreasingFactor = homTauDecreasingFactor, tauDecreasingFactorPredictor = homTauDecreasingFactorPredictor;
1783 ✗ double tauIncreasingFactor = homTauIncreasingFactor, tauIncreasingThreshold = homTauIncreasingThreshold;
1784 double preTau;
1785 ✗ int m = solverData->m;
1786 ✗ int n = solverData->n;
1787 ✗ int initialStep = 1;
1788 ✗ int maxLambdaSteps = homMaxLambdaSteps ? homMaxLambdaSteps : solverData->maxNumberOfIterations;
1789
1790 ✗ int assert = 1;
1791 ✗ DATA* data = solverData->userData->data;
1792 ✗ threadData_t *threadData = solverData->userData->threadData;
1793 ✗ int sysNumber = solverData->userData->sysNumber;
1794
1795 // TODO: Make this print a function!
1796 ✗ FILE *pFile = NULL;
1797 char buffer[4096];
1798
1799 #if !defined(OMC_NO_FILESYSTEM)
1800 ✗ const char sep[] = ",";
1801 ✗ if(solverData->initHomotopy && OMC_ACTIVE_STREAM(OMC_LOG_INIT_HOMOTOPY))
1802 {
1803 ✗ if (omc_flag[FLAG_OUTPUT_PATH]) { /* Add output path to file name */
1804 ✗ sprintf(buffer, "%s/%s_nonlinsys%d_adaptive_%s_homotopy_%s.csv", omc_flagValue[FLAG_OUTPUT_PATH], data->modelData->modelFilePrefix, sysNumber, data->callback->homotopyMethod == GLOBAL_ADAPTIVE_HOMOTOPY ? "global" : "local", solverData->startDirection > 0 ? "pos" : "neg");
1805 }
1806 else
1807 {
1808 ✗ sprintf(buffer, "%s_nonlinsys%d_adaptive_%s_homotopy_%s.csv", data->modelData->modelFilePrefix, sysNumber, data->callback->homotopyMethod == GLOBAL_ADAPTIVE_HOMOTOPY ? "global" : "local", solverData->startDirection > 0 ? "pos" : "neg");
1809 }
1810 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "The homotopy path will be exported to %s.", buffer);
1811 ✗ pFile = omc_fopen(buffer, "wt");
1812 fprintf(pFile, "\"sep=%s\"\n%s", sep, "\"lambda\"");
1813 ✗ for(i=0; i<n; ++i)
1814 ✗ fprintf(pFile, "%s\"%s\"", sep, modelInfoGetEquation(&data->modelData->modelDataXml,solverData->eqSystemNumber).vars[i]);
1815 fprintf(pFile, "\n");
1816 fprintf(pFile, "0.0");
1817 ✗ for(i=0; i<n; ++i)
1818 ✗ fprintf(pFile, "%s%.16g", sep, x[i]);
1819 fprintf(pFile, "\n");
1820 }
1821 #endif
1822
1823 /* Initialize vector dy2 using chosen startDirection */
1824 /* set start vector, lambda = 0.0 */
1825 ✗ vecCopy(solverData->n, x, solverData->y0);
1826 ✗ solverData->y0[solverData->n] = 0.0;
1827
1828 ✗ vecConst(solverData->n, 0.0, solverData->dy2);
1829 ✗ solverData->dy2[solverData->n]= solverData->startDirection;
1830 ✗ printHomotopyUnknowns(OMC_LOG_NLS_V, solverData);
1831 assert = 1;
1832 #ifndef OMC_EMCC
1833 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1834 #endif
1835 ✗ solverData->h_function(solverData, solverData->y0, solverData->hvec);
1836 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
1837 #ifndef OMC_EMCC
1838 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1839 #endif
1840 /* start iteration; stop, if lambda = solverData->y0[solverData->n] == 1 */
1841 ✗ while (solverData->y0[solverData->n]<1)
1842 {
1843 ✗ if (solverData->initHomotopy)
1844 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
1845 else
1846 ✗ infoStreamPrint(OMC_LOG_NLS_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
1847 /* Break loop, iff algorithm gets stuck or lambda accelerates to the wrong direction */
1848 ✗ if (iter>=maxTries)
1849 {
1850 ✗ if (solverData->initHomotopy) {
1851 ✗ if (preTau == tau)
1852 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nNo solution for current step size tau found and tau cannot be decreased any further.\nYou can set the minimum step size tau with:\n\t-homTauMin=<value>\nYou can also try to allow more newton steps in the corrector step with:\n\t-homMaxNewtonSteps=<value>\nor change the tolerance for the solution with:\n\t-homHEps=<value>\nYou can also try to use another backtrace stategy in the corrector step with:\n\t-homBacktraceStrategy=<fix|orthogonal>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
1853 else
1854 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe maximum number of tries for one lambda is reached (%d).\nYou can change the number of tries with:\n\t-homMaxTries=<value>\nYou can also try to allow more newton steps in the corrector step with:\n\t-homMaxNewtonSteps=<value>\nor change the tolerance for the solution with:\n\t-homHEps=<value>\nYou can also try to use another backtrace stategy in the corrector step with:\n\t-homBacktraceStrategy=<fix|orthogonal>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", iter);
1855 }
1856 else
1857 debugInt(OMC_LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge: iter = ", iter);
1858 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
1859 return -1;
1860 }
1861 ✗ if (solverData->y0[solverData->n]<(-1))
1862 {
1863 ✗ if (solverData->initHomotopy)
1864 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nlambda is smaller than -1: lambda=%g\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", solverData->y0[solverData->n]);
1865 else
1866 debugDouble(OMC_LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge: lambda = ", solverData->y0[solverData->n]);
1867 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
1868 return -1;
1869 }
1870 ✗ if (numSteps >= maxLambdaSteps)
1871 {
1872 ✗ if (solverData->initHomotopy)
1873 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe maximum number of lambda steps is reached (%d).\nYou can change the maximum number of lambda steps with:\n\t-homMaxLambdaSteps=<value>\nYou can also try to influence the step size tau with the following flags:\n\t-homTauDecFac=<value>\n\t-homTauDecFacPredictor=<value>\n\t-homTauIncFac=<value>\n\t-homTauIncThreshold=<value>\n\t-homTauMax=<value>\n\t-homTauMin=<value>\n\t-homTauStart=<value>\nor you can also set the threshold for accepting the current bending with:\n\t-homAdaptBend=<value>\nYou can also try to use another backtrace stategy in the corrector step with:\n\t-homBacktraceStrategy=<fix|orthogonal>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", maxLambdaSteps);
1874 else
1875 debugInt(OMC_LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge: numSteps = ", numSteps);
1876 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
1877 return -1;
1878 }
1879
1880 stepAccept = 0;
1881
1882 /****************************************************************************
1883 * Predictor step: Calculation of tangent vector! *
1884 ****************************************************************************/
1885 /* If a step succeeded, calculate the homotopy function and corresponding jacobian */
1886 ✗ if (iter==0)
1887 {
1888 /* Handle asserts of function calls, mainly necessary for fluid stuff */
1889 assert = 1;
1890 #ifndef OMC_EMCC
1891 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1892 #endif
1893 ✗ solverData->hJac_dh(solverData, solverData->y0, solverData->hJac);
1894 ✗ debugMatrixDouble(OMC_LOG_NLS_JAC,"Jacobian hJac:",solverData->hJac, solverData->n, solverData->n+1);
1895 ✗ scaleMatrixRows(solverData->n, solverData->m, solverData->hJac);
1896 ✗ debugMatrixDouble(OMC_LOG_NLS_JAC,"Jacobian hJac after scaling:",solverData->hJac, solverData->n, solverData->n+1);
1897 assert = 0;
1898 ✗ pos = -1; /* stable solution algorithm for solving a generalized over-determined linear system */
1899 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
1900 #ifndef OMC_EMCC
1901 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1902 #endif
1903
1904 ✗ if (assert || (solveSystemWithTotalPivotSearch(data, solverData->n, solverData->dy0, solverData->hJac, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1))
1905 {
1906 /* report solver abortion */
1907 ✗ solverData->info=-1;
1908 /* debug information */
1909 ✗ if (assert) {
1910 ✗ if (solverData->initHomotopy)
1911 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nIt was not possible to calculate the jacobian.\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
1912 else {
1913 debugString(OMC_LOG_NLS_HOMOTOPY, "Assert, when calculating Jacobian!");
1914 debugString(OMC_LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge");
1915 }
1916 } else {
1917 ✗ if (solverData->initHomotopy)
1918 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe system is singular and not solvable.\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
1919 else {
1920 debugString(OMC_LOG_NLS_HOMOTOPY, "System singular and not solvable!");
1921 debugString(OMC_LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge");
1922 }
1923 }
1924 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
1925 /* update statistics */
1926 return -1;
1927 }
1928 /* Scaling back to original variables */
1929 ✗ vecMultScaling(solverData->m, solverData->dy0, solverData->xScaling, solverData->dy0);
1930 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY, "tangent vector with original scaling:", solverData->dy0, solverData->m);
1931 ✗ debugDouble(OMC_LOG_NLS_HOMOTOPY,"length of tangent vector with original scaling: ", vec2Norm(solverData->m, solverData->dy0));
1932 // vecNormalize(solverData->m, solverData->dy0, solverData->dy0);
1933 // debugVectorDouble(OMC_LOG_NLS_HOMOTOPY, "normalized tangent vector:", solverData->dy0, solverData->m);
1934 // debugDouble(OMC_LOG_NLS_HOMOTOPY,"length of normalized tangent vector: ", vec2Norm(solverData->m, solverData->dy0));
1935
1936 /* Correct search direction, depending on the last direction (angle < 90 degree) */
1937 ✗ vecScalarProduct = vecScalarProd(solverData->m,solverData->dy0,solverData->dy2);
1938 debugDouble(OMC_LOG_NLS_HOMOTOPY,"scalar product ", vecScalarProduct);
1939 ✗ if (vecScalarProduct<0 || ((fabs(vecScalarProduct)<DBL_EPSILON) && (solverData->startDirection == -1) && initialStep))
1940 {
1941 debugInt(OMC_LOG_NLS_HOMOTOPY,"initialStep = ", initialStep);
1942 ✗ debugInt(OMC_LOG_NLS_HOMOTOPY,"solverData->startDirection = ", solverData->startDirection);
1943 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY,"step:",solverData->dy0, m);
1944 ✗ vecAddInv(solverData->m, solverData->dy0, solverData->dy0);
1945 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY,"corrected step:",solverData->dy0, m);
1946 }
1947 /* adapt tau, if lambda + tau*delta_lambda > 1 */
1948 ✗ if (fabs(solverData->dy0[solverData->n])>1e-8)
1949 {
1950 ✗ tau = fmin(tau,(1-solverData->y0[solverData->n])/fabs(solverData->dy0[solverData->n]));
1951 }
1952 }
1953
1954 assert = 1;
1955 do {
1956 /* do update and store approximated vector in yt */
1957 ✗ vecAddScal(solverData->m, solverData->y0, solverData->dy0, tau, solverData->y1);
1958
1959 /* update function value */
1960 #ifndef OMC_EMCC
1961 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1962 #endif
1963 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY,"y1 (predictor step):",solverData->y1, m);
1964 ✗ solverData->h_function(solverData, solverData->y1, solverData->hvec);
1965 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY,"hvec (predictor step):",solverData->hvec, n);
1966 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
1967 #ifndef OMC_EMCC
1968 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1969 #endif
1970 ✗ if (assert){
1971 debugString(OMC_LOG_NLS_HOMOTOPY, "Assert, when calculating function value!");
1972 debugString(OMC_LOG_NLS_HOMOTOPY, "--- decreasing step size tau in predictor step!");
1973 debugDouble(OMC_LOG_NLS_HOMOTOPY, "old tau =", tau);
1974 ✗ tau = tau/tauDecreasingFactorPredictor;
1975 debugDouble(OMC_LOG_NLS_HOMOTOPY, "new tau =", tau);
1976 }
1977 ✗ } while (assert && (tau > tauMin));
1978
1979 ✗ if (assert)
1980 {
1981 /* report solver abortion */
1982 ✗ solverData->info=-1;
1983 /* debug information */
1984 ✗ if (solverData->initHomotopy)
1985 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe step size tau cannot be decreased anymore and current tau=%g already failed.\nYou can influence the calculation of tau with the following flags:\n\t-homTauDecFac=<value>\n\t-homTauDecFacPredictor=<value>\n\t-homTauIncFac=<value>\n\t-homTauIncThreshold=<value>\n\t-homTauMax=<value>\n\t-homTauMin=<value>\n\t-homTauStart=<value>\nYou can also set the threshold for accepting the current bending with:\n\t-homAdaptBend=<value>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", tau);
1986 else {
1987 debugString(OMC_LOG_NLS_HOMOTOPY, "Assert, because tau cannot be decreased anymore and current tau already failed!");
1988 debugString(OMC_LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge");
1989 }
1990 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
1991 /* update statistics */
1992 return -1;
1993 }
1994 ✗ vecCopy(solverData->m, solverData->y1, solverData->y2);
1995 ✗ vecCopy(solverData->m, solverData->y1, solverData->yt);
1996 ✗ vecCopy(solverData->n, solverData->hvec, solverData->hvecScaled);
1997
1998 ✗ solverData->tau = tau;
1999 ✗ printHomotopyPredictorStep(OMC_LOG_NLS_HOMOTOPY, solverData);
2000
2001 /****************************************************************************
2002 * Corrector step: Newton iteration! *
2003 ****************************************************************************/
2004 debugString(OMC_LOG_NLS_HOMOTOPY, "Newton iteration for corrector step begins!");
2005
2006 /* If this is the last step, use backtrace strategy with one fixed coordinate and fix lambda */
2007 ✗ if (solverData->yt[solverData->n] == 1)
2008 {
2009 debugString(OMC_LOG_NLS_HOMOTOPY, "Force '-homBacktraceStrategy=fix' and fix lambda, because this is the last step!");
2010 ✗ debugDouble(OMC_LOG_NLS_HOMOTOPY, "Set tolerance homHEps to newtonFTol =", newtonFTol);
2011 correctorStrategy = 1;
2012 ✗ pos = solverData->n;
2013 ✗ hEps = newtonFTol;
2014 }
2015
2016 ✗ if (correctorStrategy==1)
2017 debugString(OMC_LOG_NLS_HOMOTOPY, "Using backtrace strategy with one fixed coordinate! To change this use: '-homBacktraceStrategy=orthogonal'");
2018 else
2019 debugString(OMC_LOG_NLS_HOMOTOPY, "Using backtrace strategy orthogonal to the tangent vector! To change this use: '-homBacktraceStrategy=fix'");
2020
2021
2022 ✗ for(j=0;j<maxiter;j++)
2023 {
2024 ✗ debugInt(OMC_LOG_NLS_HOMOTOPY, "Iteration: ", j+1);
2025 ✗ if (vec2Norm(solverData->n, solverData->hvec)<hEps || vec2Norm(solverData->n, solverData->hvecScaled)<hEps)
2026 {
2027 debugString(OMC_LOG_NLS_HOMOTOPY, "step accepted!");
2028 stepAccept = 1;
2029 break;
2030 }
2031 assert = 1;
2032 #ifndef OMC_EMCC
2033 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
2034 #endif
2035 /* calculate homotopy jacobian */
2036 ✗ solverData->hJac_dh(solverData, solverData->y1, solverData->hJac);
2037 ✗ debugMatrixDouble(OMC_LOG_NLS_JAC,"Jacobian hJac:",solverData->hJac, solverData->n, solverData->n+1);
2038
2039 ✗ if (correctorStrategy==2)
2040 {
2041 /* calculate the newton matrix hJac2 for the orthogonal backtrace strategy */
2042 ✗ orthogonalBacktraceMatrix(solverData, solverData->hJac, solverData->hvec, solverData->dy0, solverData->hJac2, solverData->n, solverData->m);
2043 ✗ debugMatrixDouble(OMC_LOG_NLS_JAC,"Enhanced Jacobian hJac2 (orthogonal backtrace strategy):",solverData->hJac2, solverData->n+1, solverData->m+1);
2044 }
2045
2046 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
2047 #ifndef OMC_EMCC
2048 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
2049 #endif
2050 ✗ if (assert)
2051 {
2052 debugString(OMC_LOG_NLS_HOMOTOPY, "step NOT accepted, because hJac_dh could not be calculated!");
2053 stepAccept = 0;
2054 break;
2055 }
2056 ✗ matVecMultAbs(solverData->n, solverData->m, solverData->hJac, solverData->ones, solverData->resScaling);
2057 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY, "residuum scaling of function h:", solverData->resScaling, solverData->n);
2058
2059 ✗ if (correctorStrategy==1) // fix one coordinate
2060 {
2061 /* copy vector h to column "pos" of the jacobian */
2062 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY, "copy vector hvec to column 'pos' of the jacobian:", solverData->hvec, solverData->n);
2063 ✗ vecCopy(solverData->n, solverData->hvec, solverData->hJac + pos*solverData->n);
2064 ✗ scaleMatrixRows(solverData->n, solverData->m, solverData->hJac);
2065 ✗ if (solveSystemWithTotalPivotSearch(data, solverData->n, solverData->dy1, solverData->hJac, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1)
2066 {
2067 debugString(OMC_LOG_NLS_HOMOTOPY, "step NOT accepted, because solveSystemWithTotalPivotSearch failed!");
2068 stepAccept = 0;
2069 break;
2070 }
2071 ✗ solverData->dy1[pos] = 0.0;
2072 }
2073 else // go back in orthogonal direction to tangent vector
2074 {
2075 ✗ scaleMatrixRows(solverData->n+1, solverData->m+1, solverData->hJac2);
2076 ✗ pos = solverData->n+1;
2077 ✗ if (solveSystemWithTotalPivotSearch(data, solverData->n+1, solverData->dy1, solverData->hJac2, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1)
2078 {
2079 debugString(OMC_LOG_NLS_HOMOTOPY, "step NOT accepted, because solveSystemWithTotalPivotSearch failed!");
2080 stepAccept = 0;
2081 break;
2082 }
2083 }
2084
2085 /* Scaling back to original variables */
2086 ✗ vecMultScaling(solverData->m, solverData->dy1, solverData->xScaling, solverData->dy1);
2087 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY, "solution (original scaling):", solverData->dy1, solverData->m);
2088
2089 ✗ vecAdd(solverData->m, solverData->y1, solverData->dy1, solverData->y2);
2090 ✗ vecCopy(solverData->m, solverData->y2, solverData->y1);
2091 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY, "new y in newton:", solverData->y1, solverData->m);
2092 assert = 1;
2093 #ifndef OMC_EMCC
2094 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
2095 #endif
2096 /* calculate homotopy function */
2097 ✗ solverData->h_function(solverData, solverData->y1, solverData->hvec);
2098 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { assert = 0; }
2099 #ifndef OMC_EMCC
2100 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
2101 #endif
2102 ✗ if (assert)
2103 {
2104 debugString(OMC_LOG_NLS_HOMOTOPY, "step NOT accepted, because h_function could not be calculated!");
2105 stepAccept = 0;
2106 break;
2107 }
2108 /* Calculate different error measurements */
2109 ✗ vecDivScaling(solverData->n, solverData->hvec, solverData->resScaling, solverData->hvecScaled);
2110
2111 ✗ delta_x = vec2Norm(solverData->m, solverData->dy1);
2112 ✗ error_h = vec2Norm(solverData->n, solverData->hvec);
2113 ✗ error_h_scaled = vec2Norm(solverData->n, solverData->hvecScaled);
2114
2115
2116 /* debug information */
2117 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY,"function values:",solverData->hvec, n);
2118 ✗ debugVectorDouble(OMC_LOG_NLS_HOMOTOPY,"scaled function values:",solverData->hvecScaled, n);
2119
2120 debugString(OMC_LOG_NLS_HOMOTOPY, "error measurements:");
2121 debugDouble(OMC_LOG_NLS_HOMOTOPY, "delta_x =", delta_x);
2122 debugDouble(OMC_LOG_NLS_HOMOTOPY, "error_h =", error_h);
2123 debugDouble(OMC_LOG_NLS_HOMOTOPY, "error_h_scaled =", error_h_scaled);
2124 debugDouble(OMC_LOG_NLS_HOMOTOPY, "hEps =", hEps);
2125
2126 }
2127 debugString(OMC_LOG_NLS_HOMOTOPY, "Newton iteration for corrector step finished!");
2128
2129 ✗ if (!assert)
2130 {
2131 ✗ vecDiff(solverData->m, solverData->y1, solverData->yt, solverData->dy1);
2132 ✗ vecDiff(solverData->m, solverData->yt, solverData->y0, solverData->dy2);
2133 ✗ printHomotopyCorrectorStep(OMC_LOG_NLS_HOMOTOPY, solverData);
2134 ✗ bend = vec2Norm(solverData->m,solverData->dy1)/vec2Norm(solverData->m,solverData->dy2);
2135
2136 ✗ debugDouble(OMC_LOG_NLS_HOMOTOPY, "vector length of predictor step =", vec2Norm(solverData->m,solverData->dy2));
2137 ✗ debugDouble(OMC_LOG_NLS_HOMOTOPY, "vector length of corrector step =", vec2Norm(solverData->m,solverData->dy1));
2138 debugDouble(OMC_LOG_NLS_HOMOTOPY, "bend =", bend);
2139 debugDouble(OMC_LOG_NLS_HOMOTOPY, "adaptBend =", adaptBend);
2140 }
2141 ✗ if ((bend > adaptBend) || !stepAccept)
2142 {
2143 ✗ if (bend<DBL_EPSILON)
2144 {
2145 /* debug information */
2146 ✗ if (solverData->initHomotopy)
2147 ✗ warningStreamPrint(OMC_LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe value specifying the bending of the homotopy curve is smaller than DBL_EPSILON (increment zero).\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
2148 else
2149 debugString(OMC_LOG_NLS_HOMOTOPY, "\nINCREMENT ZERO: Homotopy algorithm did not converge\n");
2150 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
2151 /* update statistics */
2152 return -1;
2153 }
2154 debugString(OMC_LOG_NLS_HOMOTOPY, "The relation between the vector length of corrector step and predictor step is too big:");
2155 ✗ debugDouble(OMC_LOG_NLS_HOMOTOPY, "bend/adaptBend =", bend/adaptBend);
2156 debugString(OMC_LOG_NLS_HOMOTOPY, "--- decreasing step size tau in corrector step!");
2157 ✗ preTau = tau;
2158 debugDouble(OMC_LOG_NLS_HOMOTOPY, "old tau =", preTau);
2159 ✗ tau = fmax(tauMin,tau/tauDecreasingFactor);
2160 debugDouble(OMC_LOG_NLS_HOMOTOPY, "new tau =", tau);
2161 ✗ if (tau==preTau)
2162 iter = maxTries;
2163 else
2164 ✗ iter++;
2165 } else
2166 {
2167 ✗ initialStep = 0;
2168 ✗ iter = 0;
2169 ✗ numSteps++;
2170 ✗ if (bend < adaptBend/tauIncreasingThreshold)
2171 {
2172 debugString(OMC_LOG_NLS_HOMOTOPY, "--- increasing step size tau in corrector step!");
2173 debugDouble(OMC_LOG_NLS_HOMOTOPY, "old tau =", tau);
2174 ✗ tau = fmin(tauMax, tau*tauIncreasingFactor);
2175 debugDouble(OMC_LOG_NLS_HOMOTOPY, "new tau =", tau);
2176 }
2177 ✗ vecCopy(solverData->m, solverData->y1, solverData->y0);
2178 ✗ vecCopy(solverData->m, solverData->dy0, solverData->dy2);
2179 debugString(OMC_LOG_NLS_HOMOTOPY, "Successfull homotopy step!\n======================================================");
2180 ✗ printHomotopyUnknowns(OMC_LOG_NLS_HOMOTOPY, solverData);
2181
2182 #if !defined(OMC_NO_FILESYSTEM)
2183 ✗ if(solverData->initHomotopy && OMC_ACTIVE_STREAM(OMC_LOG_INIT_HOMOTOPY))
2184 {
2185 ✗ fprintf(pFile, "%.16g", solverData->y0[n]);
2186 ✗ for(i=0; i<n; ++i)
2187 ✗ fprintf(pFile, "%s%.16g", sep, solverData->y0[i]);
2188 fprintf(pFile, "\n");
2189 }
2190 #endif
2191 }
2192 }
2193 ✗ if (solverData->initHomotopy)
2194 ✗ infoStreamPrint(OMC_LOG_INIT_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
2195 else
2196 ✗ infoStreamPrint(OMC_LOG_NLS_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
2197 /* copy solution back to vector x */
2198 ✗ vecCopy(solverData->n, solverData->y1, x);
2199
2200 debugString(OMC_LOG_NLS_HOMOTOPY, "HOMOTOPY ALGORITHM SUCCEEDED");
2201 ✗ if (solverData->initHomotopy) {
2202 ✗ data->simulationInfo->homotopySteps += numSteps;
2203 debugInt(OMC_LOG_INIT_HOMOTOPY, "Total number of lambda steps for this homotopy loop:", numSteps);
2204 }
2205 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
2206 ✗ solverData->info = 1;
2207
2208 #if !defined(OMC_NO_FILESYSTEM)
2209 ✗ if(solverData->initHomotopy && OMC_ACTIVE_STREAM(OMC_LOG_INIT_HOMOTOPY))
2210 ✗ fclose(pFile);
2211 #endif
2212
2213 return 0;
2214 }
2215
2216 /**
2217 * @brief Solve non-linear system with damped Newton method, combined with homotopy approach.
2218 *
2219 * @param data Pointer to data struct.
2220 * @param threadData Pointer to thread data.
2221 * @param nlsData Non-linear system data.
2222 * @return NLS_SOLVER_STATUS Return NLS_SOLVED on success and NLS_FAILED otherwise.
2223 */
2224 ✗ NLS_SOLVER_STATUS solveHomotopy(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nlsData)
2225 {
2226 ✗ DATA_HOMOTOPY* homotopyData = (DATA_HOMOTOPY*)(nlsData->solverData);
2227 DATA_HYBRD* solverDataHybrid;
2228
2229 /*
2230 * Get non-linear equation system
2231 */
2232 ✗ int eqSystemNumber = nlsData->equationIndex;
2233 ✗ int mixedSystem = nlsData->mixedSystem;
2234
2235 int i, j;
2236 ✗ NLS_SOLVER_STATUS success = NLS_FAILED;
2237 double error_f_sqrd, error_f_sqrd_scaled, error_f1_sqrd;
2238
2239 ✗ int assert = 1;
2240 ✗ int giveUp = 0;
2241 ✗ int alreadyTested = 0;
2242 int pos;
2243 int rank;
2244 ✗ int tries = 0;
2245 ✗ int runHomotopy = 0;
2246 ✗ int skipNewton = 0;
2247 ✗ homotopyData->casualTearingSet = nlsData->strictTearingFunctionCall != NULL;
2248 int constraintViolated;
2249 ✗ homotopyData->initHomotopy = nlsData->initHomotopy;
2250
2251 modelica_boolean* relationsPreBackup;
2252 ✗ relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));
2253
2254 ✗ homotopyData->f = wrapper_fvec;
2255 ✗ homotopyData->f_con = wrapper_fvec_constraints;
2256 ✗ homotopyData->fJac_f = wrapper_fvec_der;
2257
2258 ✗ homotopyData->eqSystemNumber = nlsData->equationIndex;
2259 ✗ homotopyData->mixedSystem = mixedSystem;
2260 ✗ homotopyData->timeValue = data->localData[0]->timeValue;
2261 ✗ homotopyData->minValue = nlsData->min;
2262 ✗ homotopyData->maxValue = nlsData->max;
2263 ✗ homotopyData->info = 0;
2264
2265 ✗ vecConst(homotopyData->m,1.0,homotopyData->ones);
2266
2267 ✗ if (!homotopyData->initHomotopy) {
2268 ✗ int indexes[2] = {1,eqSystemNumber};
2269 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_NLS_V, omc_dummyFileInfo, 1, indexes,
2270 "Start solving Non-Linear System %d (size %d) at time %g with Mixed (Newton/Homotopy) Solver",
2271 ✗ eqSystemNumber, (int) nlsData->size, data->localData[0]->timeValue);
2272 } else {
2273 debugString(OMC_LOG_NLS_V, "------------------------------------------------------");
2274 debugString(OMC_LOG_NLS_V, "SOLVING HOMOTOPY INITIALIZATION PROBLEM WITH THE HOMOTOPY SOLVER");
2275 debugInt(OMC_LOG_NLS_V, "EQUATION NUMBER:", eqSystemNumber);
2276 ✗ debugDouble(OMC_LOG_NLS_V, "TIME:", homotopyData->timeValue);
2277 }
2278
2279 /* set x vector */
2280 ✗ if(data->simulationInfo->discreteCall)
2281 {
2282 ✗ vecCopy(homotopyData->n, nlsData->nlsx, homotopyData->xStart);
2283 ✗ debugVectorDouble(OMC_LOG_NLS_V,"System values", homotopyData->xStart, homotopyData->n);
2284 } else
2285 {
2286 ✗ vecCopy(homotopyData->n, nlsData->nlsxExtrapolation, homotopyData->xStart);
2287 ✗ debugVectorDouble(OMC_LOG_NLS_V,"System extrapolation", homotopyData->xStart, homotopyData->n);
2288 }
2289 ✗ vecCopy(homotopyData->n, homotopyData->xStart, homotopyData->x0);
2290 // Initialize lambda variable
2291 ✗ if (homotopyData->userData->nlsData->homotopySupport && !homotopyData->initHomotopy && homotopyData->userData->nlsData->size > homotopyData->n) {
2292 ✗ homotopyData->x0[homotopyData->n] = 1.0;
2293 ✗ homotopyData->x[homotopyData->n] = 1.0;
2294 ✗ homotopyData->x1[homotopyData->n] = 1.0;
2295 } else {
2296 ✗ homotopyData->x0[homotopyData->n] = 0.0;
2297 ✗ homotopyData->x[homotopyData->n] = 0.0;
2298 ✗ homotopyData->x1[homotopyData->n] = 0.0;
2299 }
2300 /* Use actual working point for scaling */
2301 ✗ for (i=0;i<homotopyData->n;i++){
2302 ✗ homotopyData->xScaling[i] = fmax(nlsData->nominal[i],fabs(homotopyData->x0[i]));
2303 }
2304 ✗ homotopyData->xScaling[homotopyData->n] = 1.0;
2305
2306 ✗ debugVectorDouble(OMC_LOG_NLS_V,"Nominal values", nlsData->nominal, homotopyData->n);
2307 ✗ debugVectorDouble(OMC_LOG_NLS_V,"Scaling values", homotopyData->xScaling, homotopyData->m);
2308
2309
2310 ✗ if (!homotopyData->initHomotopy) {
2311 /* Handle asserts of function calls, mainly necessary for fluid stuff */
2312 assert = 1;
2313 giveUp = 1;
2314 ✗ while (tries<=2)
2315 {
2316 ✗ debugVectorDouble(OMC_LOG_NLS_V,"x0", homotopyData->x0, homotopyData->n);
2317 /* evaluate with discontinuities */
2318 ✗ if(data->simulationInfo->discreteCall)
2319 {
2320 ✗ data->simulationInfo->solveContinuous = 0;
2321 }
2322 /* evaluate with discontinuities */
2323 #ifndef OMC_EMCC
2324 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
2325 #endif
2326 ✗ if (mixedSystem)
2327 ✗ memcpy(relationsPreBackup, data->simulationInfo->relations, sizeof(modelica_boolean)*data->modelData->nRelations);
2328
2329 ✗ if (homotopyData->casualTearingSet){
2330 ✗ constraintViolated = homotopyData->f_con(homotopyData, homotopyData->x0, homotopyData->f1);
2331 ✗ if (constraintViolated){
2332 giveUp = 1;
2333 ✗ break;
2334 }
2335 }
2336 else
2337 ✗ homotopyData->f(homotopyData, homotopyData->x0, homotopyData->f1);
2338
2339 /* A raised residual is not a value: skip everything that reads f1, and
2340 leave `assert` set for the retry. */
2341 ✗ if (!OMC_ERROR_RAISED()) {
2342 /* Try to get out of here!!! */
2343 ✗ error_f_sqrd = vec2NormSqrd(homotopyData->n, homotopyData->f1);
2344 ✗ vecDivScaling(homotopyData->n, homotopyData->f1, homotopyData->resScaling, homotopyData->fvecScaled);
2345 ✗ error_f_sqrd_scaled = vec2NormSqrd(homotopyData->n, homotopyData->fvecScaled);
2346
2347 ✗ if (error_f_sqrd < homotopyData->ftol_sqrd*1e-4 || error_f_sqrd_scaled < homotopyData->ftol_sqrd*1e-4)
2348 {
2349 ✗ success = NLS_SOLVED;
2350 /* take the solution */
2351 ✗ vecCopy(homotopyData->n, homotopyData->x, nlsData->nlsx);
2352 /* reset continous flag */
2353 ✗ data->simulationInfo->solveContinuous = 0;
2354 ✗ assert = 0;
2355 } else {
2356 ✗ homotopyData->fJac_f(homotopyData, homotopyData->x0, homotopyData->fJac);
2357 ✗ vecCopy(homotopyData->n, homotopyData->f1, homotopyData->fJac + homotopyData->n*homotopyData->n);
2358 ✗ vecCopy(homotopyData->n*homotopyData->m, homotopyData->fJac, homotopyData->fJacx0);
2359 ✗ if (mixedSystem)
2360 ✗ memcpy(relationsPreBackup, data->simulationInfo->relations, sizeof(modelica_boolean)*data->modelData->nRelations);
2361 /* calculate scaling factor of residuals */
2362 ✗ matVecMultAbsBB(homotopyData->n, homotopyData->fJac, homotopyData->ones, homotopyData->resScaling);
2363 ✗ vecMakeFinite(homotopyData->n, homotopyData->resScaling);
2364 ✗ debugVectorDouble(OMC_LOG_NLS_JAC, "residuum scaling:", homotopyData->resScaling, homotopyData->n);
2365 ✗ scaleMatrixRows(homotopyData->n, homotopyData->m, homotopyData->fJac);
2366
2367 ✗ pos = homotopyData->n;
2368 ✗ assert = (solveSystemWithTotalPivotSearch(data, homotopyData->n, homotopyData->dy0, homotopyData->fJac, homotopyData->indRow, homotopyData->indCol, &pos, &rank, homotopyData->casualTearingSet) == -1);
2369 }
2370 ✗ if (!assert)
2371 debugString(OMC_LOG_NLS_V, "regular initial point!!!");
2372 }
2373 giveUp = 0;
2374 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); }
2375 #ifndef OMC_EMCC
2376 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
2377 #endif
2378 ✗ if (assert && homotopyData->casualTearingSet)
2379 {
2380 giveUp = 1;
2381 break;
2382 }
2383 ✗ if (assert)
2384 {
2385 ✗ tries += 1;
2386 }
2387 else
2388 break;
2389 /* break symmetry, when varying start values */
2390 /* try to find regular initial point, if necessary */
2391 ✗ if (tries == 1)
2392 {
2393 debugString(OMC_LOG_NLS_V, "assert handling:\t vary initial guess by +1%.");
2394 ✗ for(i = 0; i < homotopyData->n; i++)
2395 ✗ homotopyData->x0[i] = homotopyData->xStart[i] + homotopyData->xScaling[i]*(double)i/homotopyData->n*0.01;
2396 }
2397 ✗ if (tries == 2)
2398 {
2399 debugString(OMC_LOG_NLS_V,"assert handling:\t vary initial guess by +10%.");
2400 ✗ for(i = 0; i < homotopyData->n; i++)
2401 ✗ homotopyData->x0[i] = homotopyData->xStart[i] + homotopyData->xScaling[i]*(double)i/homotopyData->n*0.1;
2402 }
2403 }
2404 ✗ if (success != NLS_SOLVED) {
2405 ✗ data->simulationInfo->solveContinuous = 1;
2406 ✗ vecCopy(homotopyData->n, homotopyData->x0, homotopyData->x);
2407 ✗ vecCopy(homotopyData->n, homotopyData->f1, homotopyData->fx0);
2408 }
2409 }
2410
2411 /* start solving loop */
2412 ✗ while(!giveUp && success != NLS_SOLVED)
2413 {
2414 ✗ giveUp = 1;
2415
2416 ✗ if (!skipNewton && !homotopyData->initHomotopy){
2417
2418 /* set x vector */
2419 ✗ if(data->simulationInfo->discreteCall){
2420 ✗ memcpy(nlsData->nlsx, homotopyData->x, homotopyData->n*(sizeof(double)));
2421 }
2422 else{
2423 ✗ memcpy(nlsData->nlsxExtrapolation, homotopyData->x, homotopyData->n*(sizeof(double)));
2424 }
2425
2426 ✗ newtonAlgorithm(homotopyData, homotopyData->x);
2427
2428 // If this is the casual tearing set (only exists for dynamic tearing), break after first try
2429 ✗ if (homotopyData->info == -1 && homotopyData->casualTearingSet){
2430 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "### No Solution for the casual tearing set at the first try! ###");
2431 break;
2432 }
2433
2434 ✗ if (homotopyData->info == -1){
2435 ✗ solverDataHybrid = (DATA_HYBRD*)(homotopyData->dataHybrid);
2436 ✗ nlsData->solverData = solverDataHybrid;
2437
2438 ✗ homotopyData->info = solveHybrd(data, threadData, nlsData);
2439
2440 ✗ memcpy(homotopyData->x, nlsData->nlsx, homotopyData->n*(sizeof(double)));
2441 ✗ nlsData->solverData = homotopyData;
2442 }
2443 }
2444
2445 /* solution found */
2446 ✗ if(homotopyData->info == 1)
2447 {
2448 ✗ success = NLS_SOLVED;
2449 /* This case may be switched off, because of event chattering!!!*/
2450 ✗ if(mixedSystem && data->simulationInfo->discreteCall && (alreadyTested<1))
2451 {
2452 ✗ debugVectorBool(OMC_LOG_NLS_V,"Relations Pre vector", data->simulationInfo->relationsPre, data->modelData->nRelations);
2453 ✗ debugVectorBool(OMC_LOG_NLS_V,"Relations Backup vector", relationsPreBackup, data->modelData->nRelations);
2454 ✗ data->simulationInfo->solveContinuous = 0;
2455
2456 ✗ if (homotopyData->casualTearingSet){
2457 ✗ constraintViolated = homotopyData->f_con(homotopyData, homotopyData->x, homotopyData->f1);
2458 ✗ if (constraintViolated){
2459 success = NLS_FAILED;
2460 break;
2461 }
2462 }
2463 else
2464 ✗ homotopyData->f(homotopyData, homotopyData->x, homotopyData->f1);
2465
2466 ✗ debugVectorBool(OMC_LOG_NLS_V,"Relations vector", data->simulationInfo->relations, data->modelData->nRelations);
2467 ✗ if (isNotEqualVectorInt(data->modelData->nRelations, data->simulationInfo->relations, relationsPreBackup)>0)
2468 {
2469 /* re-run the solution process, since relations in the system have changed */
2470 ✗ success = NLS_FAILED;
2471 ✗ giveUp = 0;
2472 ✗ runHomotopy = 0;
2473 ✗ alreadyTested = 1;
2474 ✗ vecCopy(homotopyData->n, homotopyData->x0, homotopyData->x);
2475 ✗ vecCopy(homotopyData->n, homotopyData->fx0, homotopyData->f1);
2476 ✗ vecCopy(homotopyData->n*homotopyData->m, homotopyData->fJacx0, homotopyData->fJac);
2477
2478 /* calculate scaling factor of residuals */
2479 ✗ matVecMultAbsBB(homotopyData->n, homotopyData->fJac, homotopyData->ones, homotopyData->resScaling);
2480 ✗ scaleMatrixRows(homotopyData->n, homotopyData->m, homotopyData->fJac);
2481
2482 ✗ pos = homotopyData->n;
2483 ✗ solveSystemWithTotalPivotSearch(data, homotopyData->n, homotopyData->dy0, homotopyData->fJac, homotopyData->indRow, homotopyData->indCol, &pos, &rank, homotopyData->casualTearingSet);
2484 ✗ debugDouble(OMC_LOG_NLS_V,"solve mixed system at time : ", homotopyData->timeValue);
2485 ✗ continue;
2486 }
2487 }
2488 if (success == NLS_SOLVED)
2489 {
2490 /* take the solution */
2491 ✗ vecCopy(homotopyData->n, homotopyData->x, nlsData->nlsx);
2492 /* reset continous flag */
2493 ✗ data->simulationInfo->solveContinuous = 0;
2494 ✗ break;
2495 }
2496 }
2497 ✗ if (success != NLS_SOLVED && runHomotopy>=3) break;
2498 /* Start homotopy search for new start values */
2499 ✗ vecCopy(homotopyData->n, homotopyData->x0, homotopyData->x);
2500 ✗ runHomotopy++;
2501 /* debug output */
2502 debugString(OMC_LOG_NLS_HOMOTOPY, "======================================================");
2503
2504 ✗ if (homotopyData->initHomotopy) {
2505 ✗ if (runHomotopy == 1) {
2506 ✗ homotopyData->h_function = wrapper_fvec;
2507 ✗ homotopyData->hJac_dh = wrapper_fvec_der;
2508 ✗ homotopyData->startDirection = omc_flag[FLAG_HOMOTOPY_NEG_START_DIR] ? -1.0 : 1.0;
2509 debugInt(OMC_LOG_INIT_HOMOTOPY, "Homotopy run: ", runHomotopy);
2510 ✗ debugDouble(OMC_LOG_INIT_HOMOTOPY,"startDirection = ", homotopyData->startDirection);
2511 }
2512
2513 ✗ if (runHomotopy == 2) {
2514 ✗ homotopyData->h_function = wrapper_fvec;
2515 ✗ homotopyData->hJac_dh = wrapper_fvec_der;
2516 ✗ homotopyData->startDirection = omc_flag[FLAG_HOMOTOPY_NEG_START_DIR] ? 1.0 : -1.0;
2517 ✗ infoStreamPrint(OMC_LOG_ASSERT, 0, "The homotopy algorithm is started again with opposing start direction.");
2518 debugInt(OMC_LOG_INIT_HOMOTOPY, "Homotopy run: ", runHomotopy);
2519 ✗ debugDouble(OMC_LOG_INIT_HOMOTOPY,"Try again with startDirection = ", homotopyData->startDirection);
2520 }
2521
2522 ✗ if (runHomotopy == 3) {
2523 success = NLS_FAILED;
2524 break;
2525 }
2526 }
2527 else {
2528 debugInt(OMC_LOG_NLS_HOMOTOPY, "Homotopy run: ", runHomotopy);
2529 ✗ if (runHomotopy == 1)
2530 {
2531 /* store x0 and calculate f(x0) -> newton homotopy, fJac(x0) -> taylor, affin homotopy */
2532 ✗ homotopyData->h_function = wrapper_fvec_homotopy_newton;
2533 ✗ homotopyData->hJac_dh = wrapper_fvec_homotopy_newton_der;
2534 ✗ homotopyData->startDirection = 1.0;
2535 debugDouble(OMC_LOG_NLS_HOMOTOPY,"STARTING NEWTON HOMOTOPY METHOD; startDirection = ", homotopyData->startDirection);
2536 }
2537 ✗ if (runHomotopy == 2)
2538 {
2539 /* store x0 and calculate f(x0) -> newton homotopy, fJac(x0) -> taylor, affin homotopy */
2540 ✗ homotopyData->h_function = wrapper_fvec_homotopy_newton;
2541 ✗ homotopyData->hJac_dh = wrapper_fvec_homotopy_newton_der;
2542 ✗ homotopyData->startDirection = -1.0;
2543 debugDouble(OMC_LOG_NLS_HOMOTOPY,"STARTING NEWTON HOMOTOPY METHOD; startDirection = ", homotopyData->startDirection);
2544 }
2545 ✗ if (runHomotopy == 3)
2546 {
2547 ✗ homotopyData->h_function = wrapper_fvec_homotopy_fixpoint;
2548 ✗ homotopyData->hJac_dh = wrapper_fvec_homotopy_fixpoint_der;
2549 ✗ homotopyData->startDirection = 1.0;
2550 debugDouble(OMC_LOG_NLS_HOMOTOPY,"STARTING FIXPOINT HOMOTOPY METHOD = ", homotopyData->startDirection);
2551 }
2552 }
2553
2554 ✗ homotopyAlgorithm(homotopyData, homotopyData->x);
2555
2556 ✗ if (homotopyData->info<1)
2557 {
2558 skipNewton = 1;
2559 ✗ giveUp = runHomotopy>=3;
2560
2561 ✗ } else if (homotopyData->initHomotopy && homotopyData->info==1) {
2562 /* take the solution */
2563 ✗ vecCopy(homotopyData->n, homotopyData->x, nlsData->nlsx);
2564 ✗ debugVectorDouble(OMC_LOG_NLS_V,"Solution", homotopyData->x, homotopyData->n);
2565 success = NLS_SOLVED;
2566 }
2567
2568 else {
2569 assert = 1;
2570 #ifndef OMC_EMCC
2571 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
2572 #endif
2573 ✗ if (homotopyData->casualTearingSet){
2574 ✗ constraintViolated = homotopyData->f_con(homotopyData, homotopyData->x, homotopyData->f1);
2575 ✗ if (constraintViolated){
2576 success = NLS_FAILED;
2577 ✗ break;
2578 }
2579 }
2580 else
2581 ✗ homotopyData->f(homotopyData, homotopyData->x, homotopyData->f1);
2582
2583 ✗ homotopyData->fJac_f(homotopyData, homotopyData->x, homotopyData->fJac);
2584 ✗ vecCopy(homotopyData->n, homotopyData->f1, homotopyData->fJac + homotopyData->n*homotopyData->n);
2585 /* calculate scaling factor of residuals */
2586 ✗ matVecMultAbsBB(homotopyData->n, homotopyData->fJac, homotopyData->ones, homotopyData->resScaling);
2587 ✗ debugVectorDouble(OMC_LOG_NLS_JAC, "residuum scaling:", homotopyData->resScaling, homotopyData->n);
2588 ✗ scaleMatrixRows(homotopyData->n, homotopyData->m, homotopyData->fJac);
2589
2590 ✗ pos = homotopyData->n;
2591 ✗ assert = (solveSystemWithTotalPivotSearch(data, homotopyData->n, homotopyData->dy0, homotopyData->fJac, homotopyData->indRow, homotopyData->indCol, &pos, &rank, homotopyData->casualTearingSet) == -1);
2592 ✗ if (!assert)
2593 debugString(OMC_LOG_NLS_V, "regular initial point!!!");
2594 /* As above: a raised residual leaves the retry armed. */
2595 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); assert = 1; }
2596 #ifndef OMC_EMCC
2597 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
2598 #endif
2599 ✗ if (assert)
2600 {
2601 giveUp = 1;
2602 } else
2603 {
2604 giveUp = 0;
2605 skipNewton = 0;
2606 }
2607 }
2608 }
2609 ✗ if (success != NLS_SOLVED)
2610 {
2611 debugString(OMC_LOG_NLS_V,"Homotopy solver did not converge!");
2612 }
2613 ✗ free(relationsPreBackup);
2614
2615 ✗ if (!homotopyData->initHomotopy) {
2616 ✗ messageClose(OMC_LOG_NLS_V);
2617 }
2618
2619 /* write statistics */
2620 ✗ nlsData->numberOfFEval = homotopyData->numberOfFunctionEvaluations;
2621 ✗ nlsData->numberOfIterations = homotopyData->numberOfIterations;
2622
2623 ✗ return success;
2624 }
2625
2626 /**
2627 * @brief Return pointer to Jacobian.
2628 *
2629 * @param nlsData Non-linear system data.
2630 * @return double* Jacobian in row-major format.
2631 */
2632 ✗ double* getHomotopyJacobian(NONLINEAR_SYSTEM_DATA* nlsData) {
2633 ✗ DATA_HOMOTOPY* homotopyData = (DATA_HOMOTOPY*)(nlsData->solverData);
2634 ✗ return homotopyData->fJac;
2635 }
2636
2637 #endif
2638