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 / 716
Functions: 0.0% 0 / 0 / 27
Branches: 0.0% 0 / 0 / 362
OMCompiler/SimulationRuntime/c/simulation/solver/kinsol_b.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 // this a quick copy / rewrite of the Kinsol NLS interface
29 // it is only meant for experimental builds in order to test scalings
30
31 /*! \file kinsol_B.c
32 */
33
34
35 #include "kinsol_b.h"
36
37 #include "nonlinearSystem.h"
38 #include "omc_config.h"
39 #include "omc_math.h"
40 #include "../options.h"
41 #include "../simulation_info_json.h"
42 #include "../jacobian_util.h"
43 #include "sundials_util.h"
44 #include "util/omc_error.h"
45
46 #ifdef WITH_SUNDIALS
47
48 #include "events.h"
49 #include "model_help.h"
50 #include "openmodelica.h"
51 #include "../results/simulation_result_rust.h"
52 #include "openmodelica_func.h"
53 #include "util/read_matlab4.h"
54 #include "util/varinfo.h"
55
56 #include <math.h>
57 #include <stdio.h>
58 #include <stdlib.h>
59 #include <string.h>
60
61 /* Function prototypes */
62 static int B_nlsKinsolResiduals(N_Vector x, N_Vector f, void* userData);
63 static int B_nlsSparseJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
64 void* userData, N_Vector tmp1, N_Vector tmp2);
65 static int B_nlsSparseSymJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
66 void* userData, N_Vector tmp1, N_Vector tmp2);
67 static int B_nlsDenseJac(long int N, N_Vector vecX, N_Vector vecFX,
68 SUNMatrix Jac, NLS_USERDATA *kinsolUserData,
69 N_Vector tmp1, N_Vector tmp2);
70 static void B_nlsKinsolJacSumSparse(SUNMatrix A);
71 static void B_nlsKinsolJacSumDense(SUNMatrix A);
72
73 static void B_print_jac(NONLINEAR_SYSTEM_DATA* nlsData, SUNMatrix J, const char* name) {
74 int i, j, size, nz, nnz, col, row;
75
76 SPARSE_PATTERN *sp = nlsData->sparsePattern;
77 size = nlsData->size;
78 if (SUNMatGetID(J) == SUNMATRIX_DENSE) {
79 for (col = 0; col < size; col++) {
80 for (row = 0; row < size; row++) {
81 infoStreamPrint(OMC_LOG_STDOUT, 0, "%s(row = %d, col = %d) = %.3e", name, row, col, SM_ELEMENT_D(J, row, col));
82 }
83 }
84 }
85 else if (SUNMatGetID(J) == SUNMATRIX_SPARSE && SM_SPARSETYPE_S(J) == SUN_CSC_MAT) {
86 for (col = 0; col < size; col++) {
87 for (nz = sp->leadindex[col]; nz < sp->leadindex[col + 1]; nz++) {
88 row = sp->index[nz];
89 infoStreamPrint(OMC_LOG_STDOUT, 0, "%s(row = %d, col = %d) = %.3e", name, row, col, SM_DATA_S(J)[nz]);
90 }
91 }
92 }
93 }
94
95 // debug print
96 static void B_print_X(B_NLS_KINSOL_DATA* kinsolData, N_Vector x, const char* name) {
97 int i, j, size, nz, nnz, col, row;
98 sunrealtype *values;
99
100 size = kinsolData->size;
101 values = N_VGetArrayPointer(x);
102 for (i = 0; i < size; i++) {
103 infoStreamPrint(OMC_LOG_STDOUT, 0, "%s(%d) = %.3e", name, i, values[i]);
104 }
105 }
106
107 ✗ static void nlsKinsolInplaceScaleJac(NONLINEAR_SYSTEM_DATA *nlsData, B_NLS_KINSOL_DATA *kinsolData, SUNMatrix Jac) {
108 /* scaling pointers */
109 int i, j, size, nz, col, row;
110 ✗ SPARSE_PATTERN *sp = nlsData->sparsePattern;
111
112 ✗ double *x_scaling = N_VGetArrayPointer(kinsolData->xScale);
113 ✗ double *f_scaling = N_VGetArrayPointer(kinsolData->fScale);
114
115 ✗ size = kinsolData->size;
116
117 ✗ if (SUNMatGetID(Jac) == SUNMATRIX_DENSE) {
118 ✗ for (col = 0; col < size; col++) {
119 ✗ for (row = 0; row < size; row++) {
120 ✗ SM_ELEMENT_D(Jac, row, col) *= f_scaling[row] / x_scaling[col];
121 }
122 }
123 }
124 ✗ else if (SUNMatGetID(Jac) == SUNMATRIX_SPARSE && SM_SPARSETYPE_S(Jac) == SUN_CSC_MAT) {
125 ✗ for (col = 0; col < size; col++) {
126 ✗ for (nz = sp->leadindex[col]; nz < sp->leadindex[col + 1]; nz++) {
127 ✗ row = sp->index[nz];
128 ✗ SM_DATA_S(Jac)[nz] *= f_scaling[row] / x_scaling[col];
129 }
130 }
131 }
132 else {
133 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "kinsol-experimental: Matrix not supported in nlsKinsolInplaceScaleJac.");
134 }
135 ✗ }
136
137 ✗ static void nlsKinsolInplaceUnscaleJac(NONLINEAR_SYSTEM_DATA *nlsData, B_NLS_KINSOL_DATA *kinsolData, SUNMatrix Jac) {
138 /* scaling pointers */
139 int i, j, size, nz, col, row;
140 ✗ SPARSE_PATTERN *sp = nlsData->sparsePattern;
141
142 ✗ double *x_scaling = N_VGetArrayPointer(kinsolData->xScale);
143 ✗ double *f_scaling = N_VGetArrayPointer(kinsolData->fScale);
144
145 ✗ size = kinsolData->size;
146
147 ✗ if (SUNMatGetID(Jac) == SUNMATRIX_DENSE) {
148 ✗ for (col = 0; col < size; col++) {
149 ✗ for (row = 0; row < size; row++) {
150 ✗ SM_ELEMENT_D(Jac, row, col) *= x_scaling[col] / f_scaling[row];
151 }
152 }
153 }
154 ✗ else if (SUNMatGetID(Jac) == SUNMATRIX_SPARSE && SM_SPARSETYPE_S(Jac) == SUN_CSC_MAT) {
155 ✗ for (col = 0; col < size; col++) {
156 ✗ for (nz = sp->leadindex[col]; nz < sp->leadindex[col + 1]; nz++) {
157 ✗ row = sp->index[nz];
158 ✗ SM_DATA_S(Jac)[nz] *= x_scaling[col] / f_scaling[row];
159 }
160 }
161 }
162 else {
163 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "kinsol-experimental: Matrix not supported in nlsKinsolInplaceScaleJac.");
164 }
165 ✗ }
166
167 ✗ static void nlsKinsolInplaceScaleX(B_NLS_KINSOL_DATA *kinsolData, N_Vector x) {
168 ✗ int i, size = kinsolData->size;
169 ✗ double *x_data = N_VGetArrayPointer(x);
170 ✗ double *x_scaling = N_VGetArrayPointer(kinsolData->xScale);
171
172 ✗ for (i = 0; i < size; i++) {
173 ✗ x_data[i] *= x_scaling[i];
174 }
175 ✗ }
176
177 ✗ static void nlsKinsolInplaceScaleF(B_NLS_KINSOL_DATA *kinsolData, N_Vector f) {
178 ✗ int i, size = kinsolData->size;
179 ✗ double *f_data = N_VGetArrayPointer(f);
180 ✗ double *f_scaling = N_VGetArrayPointer(kinsolData->fScale);
181
182 ✗ for (i = 0; i < size; i++) {
183 ✗ f_data[i] *= f_scaling[i];
184 }
185 ✗ }
186
187 ✗ static void nlsKinsolInplaceUnscaleX(B_NLS_KINSOL_DATA *kinsolData, N_Vector x) {
188 ✗ int i, size = kinsolData->size;
189 ✗ double *x_data = N_VGetArrayPointer(x);
190 ✗ double *x_scaling = N_VGetArrayPointer(kinsolData->xScale);
191
192 ✗ for (i = 0; i < size; i++) {
193 ✗ x_data[i] /= x_scaling[i];
194 }
195 ✗ }
196
197 /**
198 * @brief Set KINSOL configuration.
199 *
200 * @param kinsolData Kinsol data with configuration settings.
201 */
202 ✗ static void B_nlsKinsolConfigSetup(B_NLS_KINSOL_DATA *kinsolData) {
203 /* Variables */
204 int flag;
205
206 /* configuration */
207 ✗ flag = KINSetFuncNormTol(kinsolData->kinsolMemory,
208 kinsolData->fnormtol); /* Set function-norm stopping tolerance */
209 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetFuncNormTol");
210 ✗ kinsolData->resetTol = FALSE;
211
212 ✗ flag = KINSetScaledStepTol(kinsolData->kinsolMemory,
213 kinsolData->scsteptol); /* Set scaled-step stopping tolerance */
214 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetScaledStepTol");
215
216 ✗ flag = KINSetNumMaxIters(kinsolData->kinsolMemory,
217 ✗ 100 * kinsolData->size); /* Set max. number of nonlinear iterations */
218 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNumMaxIters");
219
220 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH; /* Newton with globalization strategy to solve nonlinear systems */
221
222 ✗ flag = KINSetNoInitSetup(kinsolData->kinsolMemory, SUNFALSE); /* TODO: This is the default value. Is there a point in calling this function? */
223 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNoInitSetup");
224
225 ✗ kinsolData->retries = 0;
226 ✗ kinsolData->countResCalls = 0;
227 ✗ }
228
229 /**
230 * @brief Error handler function given to the SUNContext.
231 *
232 * @param line Line in the SUNDIALS source where the error was raised.
233 * @param func Name of the SUNDIALS function in which the error occurred.
234 * @param file SUNDIALS source file where the error was raised.
235 * @param msg Error message.
236 * @param err_code SUNDIALS error code.
237 * @param err_user_data Pointer to user data given with SUNContext_PushErrHandler.
238 * @param sunctx SUNDIALS context, unused.
239 */
240 ✗ static void B_kinsolErrorHandlerFunction(int line, const char *func, const char *file,
241 const char *msg, SUNErrCode err_code,
242 void *err_user_data, SUNContext sunctx) {
243 /* Variables */
244 B_NLS_KINSOL_DATA* kinsolData;
245 DATA* data;
246 NONLINEAR_SYSTEM_DATA* nlsData;
247 long eqSystemNumber = -1;
248
249 (void)(sunctx); /* Disables compiler warning */
250
251 ✗ if (err_user_data != NULL) {
252 kinsolData = (B_NLS_KINSOL_DATA *)err_user_data;
253 ✗ data = kinsolData->userData->data;
254 ✗ nlsData = kinsolData->userData->nlsData;
255 ✗ if (nlsData) {
256 ✗ eqSystemNumber = nlsData->equationIndex;
257 }
258 }
259
260 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS)) {
261 ✗ if (err_user_data != NULL && eqSystemNumber > 0) {
262 ✗ warningStreamPrint(
263 OMC_LOG_NLS, 1, "kinsol failed for system %d",
264 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber).id);
265 } else {
266 ✗ warningStreamPrint(
267 OMC_LOG_NLS, 1, "kinsol failed");
268 }
269
270 ✗ warningStreamPrint(OMC_LOG_NLS, 0,
271 "[function] %s | [at] %s:%d | [error_code] %d",
272 func, file, line, err_code);
273 /* Package level codes (KIN_* and friends) are not SUNErrCodes, so
274 * SUNGetErrMsg() only makes sense when SUNDIALS did not supply a message. */
275 ✗ warningStreamPrint(OMC_LOG_NLS, 0, "%s", msg ? msg : SUNGetErrMsg(err_code));
276
277 ✗ messageCloseWarning(OMC_LOG_NLS);
278 }
279 ✗ }
280
281 /**
282 * @brief Initialize KINSOL data.
283 *
284 * Allocate memory for KINSOL data and Jacobian.
285 *
286 * @param kinsolData KINSOL data.
287 */
288 ✗ static void initKinsolMemory(B_NLS_KINSOL_DATA *kinsolData) {
289 int flag;
290 ✗ int size = kinsolData->size;
291 ✗ NONLINEAR_SYSTEM_DATA *nlsData = kinsolData->userData->nlsData;
292 ✗ SPARSE_PATTERN* sparsePattern = nlsData->sparsePattern;
293
294 /* Free KINSOL memory block */
295 ✗ if (kinsolData->kinsolMemory != NULL || kinsolData->J != NULL || kinsolData->scaledJ != NULL) {
296 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
297 "experimental-kinsol: Already allocated kinsol memory. Loosing memory!");
298 }
299
300 /* Create KINSOL memory block */
301 /* The SUNDIALS context was created by B_nlsKinsolAllocate, which has to happen
302 * before any SUNDIALS object. */
303 ✗ kinsolData->kinsolMemory = KINCreate(kinsolData->sunctx);
304 ✗ if (kinsolData->kinsolMemory == NULL) {
305 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
306 "experimental-kinsol: In function KINCreate: An error occurred.");
307 }
308
309 /* Set error handler and print level */
310 ✗ flag = KINSetUserData(kinsolData->kinsolMemory, (void*)kinsolData->userData);
311 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetUserData");
312
313 /* Initialize KINSOL object */
314 ✗ flag = KINInit(kinsolData->kinsolMemory, B_nlsKinsolResiduals,
315 kinsolData->initialGuess);
316 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINInit");
317
318 /* Create matrix object */
319 ✗ if (kinsolData->linearSolverMethod == NLS_LS_DEFAULT ||
320 kinsolData->linearSolverMethod == NLS_LS_LAPACK) {
321 ✗ kinsolData->J = SUNDenseMatrix(size, size, kinsolData->sunctx);
322 ✗ } else if (kinsolData->linearSolverMethod == NLS_LS_KLU) {
323 ✗ if (!sparsePattern) {
324 ✗ kinsolData->nnz = size*size;
325 } else {
326 ✗ kinsolData->nnz = sparsePattern->nnz;
327 }
328 ✗ kinsolData->J = SUNSparseMatrix(size, size, kinsolData->nnz, SUN_CSC_MAT, kinsolData->sunctx);
329 ✗ kinsolData->scaledJ = SUNSparseMatrix(size, size, kinsolData->nnz, SUN_CSC_MAT, kinsolData->sunctx);
330 }
331
332 /* Create linear solver object */
333 ✗ if (kinsolData->linearSolverMethod == NLS_LS_DEFAULT ||
334 kinsolData->linearSolverMethod == NLS_LS_TOTALPIVOT) {
335 ✗ kinsolData->linSol = SUNLinSol_Dense(kinsolData->y, kinsolData->J, kinsolData->sunctx);
336 ✗ if (kinsolData->linSol == NULL) {
337 ✗ throwStreamPrint(NULL, "experimental-kinsol: In function SUNLinSol_Dense: Input incompatible.");
338 }
339 ✗ } else if (kinsolData->linearSolverMethod == NLS_LS_LAPACK) {
340 ✗ kinsolData->linSol = SUNLinSol_LapackDense(kinsolData->y, kinsolData->J, kinsolData->sunctx);
341 ✗ if (kinsolData->linSol == NULL) {
342 ✗ throwStreamPrint(NULL, "experimental-kinsol: In function SUNLinSol_LapackDense: Input incompatible.");
343 }
344 ✗ } else if (kinsolData->linearSolverMethod == NLS_LS_KLU) {
345 ✗ kinsolData->linSol = SUNLinSol_KLU(kinsolData->y, kinsolData->J, kinsolData->sunctx);
346 ✗ if (kinsolData->linSol == NULL) {
347 ✗ throwStreamPrint(NULL, "experimental-kinsol: In function SUNLinSol_KLU: Input incompatible.");
348 }
349 } else {
350 ✗ throwStreamPrint(NULL, "experimental-kinsol: Unknown linear solver method.");
351 }
352 /* Log used solver */
353 ✗ infoStreamPrint(OMC_LOG_NLS, 0, "experimental-kinsol: Using linear solver method %s", NLS_LS_METHOD_NAME[kinsolData->linearSolverMethod]);
354
355 /* Set linear solver */
356 ✗ flag = KINSetLinearSolver(kinsolData->kinsolMemory, kinsolData->linSol,
357 kinsolData->J);
358 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KINLS_FLAG, "KINSetLinearSolver");
359
360 /* Set Jacobian for non-linear solver */
361 ✗ if (kinsolData->linearSolverMethod == NLS_LS_KLU) {
362 ✗ if (nlsData->analyticalJacobianColumn != NULL && sparsePattern != NULL) {
363 ✗ flag = KINSetJacFn(kinsolData->kinsolMemory, B_nlsSparseSymJac); /* Use symbolic Jacobian with sparsity pattern*/
364 ✗ } else if (sparsePattern != NULL) {
365 ✗ flag = KINSetJacFn(kinsolData->kinsolMemory, B_nlsSparseJac); /* Use numeric Jacobian with sparsity pattern */
366 } else {
367 ✗ throwStreamPrint(NULL, "experimental-kinsol: In function initKinsolMemory: Sparse linear solver KLU needs sparse Jacobian, but no sparsity pattern is available. Use a dense non-linear solver instead of KINSOL.");
368 }
369 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KINLS_FLAG, "KINSetJacFn");
370 }
371
372 /* Configuration */
373 ✗ B_nlsKinsolConfigSetup(kinsolData);
374 ✗ }
375
376 /**
377 * @brief Allocate memory for kinsol solver data and initialize KINSOL solver.
378 *
379 * @param size Size of non-linear problem.
380 * @param userData Pointer to set NLS user data.
381 * @param attemptRetry True if KINSOL should retry with different settings after solution failed.
382 * @param isPatternAvailable True if sparsity pattern of Jacobian is available. Allocate work vectors for KLU in that case.
383 * @return B_NLS_KINSOL_DATA* Pointer to allocated KINSOL data.
384 */
385 ✗ B_NLS_KINSOL_DATA* B_nlsKinsolAllocate(int size, NLS_USERDATA* userData, modelica_boolean attemptRetry, modelica_boolean isPatternAvailable) {
386 /* Allocate system data */
387 ✗ B_NLS_KINSOL_DATA *kinsolData = (B_NLS_KINSOL_DATA *)calloc(1, sizeof(B_NLS_KINSOL_DATA));
388 int i;
389 double *ones_x, *ones_f;
390
391 ✗ kinsolData->size = size;
392 ✗ kinsolData->linearSolverMethod = userData->nlsData->nlsLinearSolver;
393 ✗ kinsolData->solved = NLS_FAILED;
394 ✗ kinsolData->userData = userData;
395
396 ✗ if (SUNContext_Create(SUN_COMM_NULL, &kinsolData->sunctx) != SUN_SUCCESS) {
397 ✗ throwStreamPrint(NULL, "experimental-kinsol: In function SUNContext_Create: An error occurred.");
398 }
399 ✗ sundialsSilenceLogger(kinsolData->sunctx);
400
401 /* Set error handler */
402 ✗ if (SUNContext_PushErrHandler(kinsolData->sunctx, B_kinsolErrorHandlerFunction, kinsolData) != SUN_SUCCESS) {
403 ✗ throwStreamPrint(NULL, "experimental-kinsol: In function SUNContext_PushErrHandler: An error occurred.");
404 }
405
406 ✗ kinsolData->fnormtol = newtonFTol; /* function tolerance */
407 ✗ kinsolData->scsteptol = newtonXTol; /* step tolerance */
408
409 ✗ kinsolData->maxstepfactor = maxStepFactor; /* step tolerance */
410 ✗ kinsolData->useScaling = FALSE; /* calculate for scaling the scaled matrix */
411 ✗ kinsolData->attemptRetry = attemptRetry;
412
413 ✗ kinsolData->initialGuess = N_VNew_Serial(size, kinsolData->sunctx);
414 ✗ kinsolData->xScale = N_VNew_Serial(size, kinsolData->sunctx);
415 ✗ kinsolData->fScale = N_VNew_Serial(size, kinsolData->sunctx);
416 ✗ kinsolData->ONES_xScale = N_VNew_Serial(size, kinsolData->sunctx);
417 ✗ kinsolData->ONES_fScale = N_VNew_Serial(size, kinsolData->sunctx);
418 ✗ ones_x = N_VGetArrayPointer(kinsolData->ONES_xScale);
419 ✗ ones_f = N_VGetArrayPointer(kinsolData->ONES_fScale);
420
421 ✗ for (i = 0; i < size; i++) {
422 ✗ ones_x[i] = 1.0;
423 ✗ ones_f[i] = 1.0;
424 }
425
426 ✗ kinsolData->fRes = N_VNew_Serial(size, kinsolData->sunctx);
427 ✗ kinsolData->fTmp = N_VNew_Serial(size, kinsolData->sunctx);
428
429 ✗ kinsolData->y = N_VNew_Serial(size, kinsolData->sunctx);
430 ✗ kinsolData->J = NULL;
431
432 /* tmp1, tmp2 only needed for numeric Jacobian */
433 ✗ if (userData->nlsData->analyticalJacobianColumn != NULL &&
434 ✗ isPatternAvailable &&
435 ✗ kinsolData->linearSolverMethod == NLS_LS_KLU)
436 {
437 ✗ kinsolData->tmp1 = NULL;
438 ✗ kinsolData->tmp2 = NULL;
439 } else {
440 ✗ kinsolData->tmp1 = N_VNew_Serial(size, kinsolData->sunctx);
441 ✗ kinsolData->tmp2 = N_VNew_Serial(size, kinsolData->sunctx);
442 }
443 /* Scaled Jacobian is allocated with J */
444 ✗ kinsolData->scaledJ = NULL;
445
446 ✗ kinsolData->kinsolMemory = NULL;
447
448 ✗ initKinsolMemory(kinsolData);
449
450 ✗ return kinsolData;
451 }
452
453 /**
454 * @brief Deallocates memory for KINSOL solver.
455 *
456 * Free memory that was allocated with `nlsKinsolAllocate`.
457 *
458 * @param kinsolData Pointer to KINSOL data.
459 */
460 ✗ void B_nlsKinsolFree(B_NLS_KINSOL_DATA* kinsolData) {
461 ✗ KINFree((void *)&kinsolData->kinsolMemory);
462
463 ✗ N_VDestroy_Serial(kinsolData->initialGuess);
464 ✗ N_VDestroy_Serial(kinsolData->xScale);
465 ✗ N_VDestroy_Serial(kinsolData->fScale);
466 ✗ N_VDestroy_Serial(kinsolData->ONES_xScale);
467 ✗ N_VDestroy_Serial(kinsolData->ONES_fScale);
468 ✗ N_VDestroy_Serial(kinsolData->fRes);
469 ✗ N_VDestroy_Serial(kinsolData->fTmp);
470
471 /* Free linear solver data */
472 ✗ SUNLinSolFree(kinsolData->linSol);
473 ✗ SUNMatDestroy(kinsolData->J);
474 ✗ SUNMatDestroy(kinsolData->scaledJ);
475 ✗ N_VDestroy_Serial(kinsolData->y);
476 ✗ if (kinsolData->tmp1 != NULL) {
477 ✗ N_VDestroy_Serial(kinsolData->tmp1);
478 ✗ N_VDestroy_Serial(kinsolData->tmp2);
479 }
480
481 /* The context has to outlive every SUNDIALS object created with it */
482 ✗ SUNContext_Free(&kinsolData->sunctx);
483
484 ✗ freeNlsUserData(kinsolData->userData);
485 ✗ free(kinsolData);
486
487 ✗ return;
488 }
489
490 /**
491 * @brief Residual function for non-linear problem.
492 *
493 * @param x The current value of the variable vector.
494 * @param f Output vector.
495 * @param userData Pointer to Kinsol user data.
496 * @return int Return 0 on success, return 1 on recoverable error.
497 */
498 ✗ static int B_nlsKinsolResiduals(N_Vector x, N_Vector f, void* userData) {
499
500 ✗ double *xdata = NV_DATA_S(x);
501 ✗ double *fdata = NV_DATA_S(f);
502
503 NLS_USERDATA* kinsolUserData = (NLS_USERDATA*)userData;
504 ✗ DATA* data = kinsolUserData->data;
505 ✗ threadData_t* threadData = kinsolUserData->threadData;
506 ✗ NONLINEAR_SYSTEM_DATA* nlsData = kinsolUserData->nlsData;
507 ✗ B_NLS_KINSOL_DATA* kinsolData = (B_NLS_KINSOL_DATA*)nlsData->solverData;
508 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=kinsolUserData->solverData};
509 ✗ int iflag = 1 /* recoverable error */;
510
511 /* Update statistics */
512 ✗ kinsolData->countResCalls++;
513
514 #ifndef OMC_EMCC
515 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
516 #endif
517
518 ✗ if (kinsolData->useScaling) {
519 ✗ nlsKinsolInplaceUnscaleX(kinsolData, x);
520 }
521
522 /* call residual function */
523 ✗ nlsData->residualFunc(&resUserData, xdata, fdata, (const int *)&iflag);
524
525 ✗ if (kinsolData->useScaling) {
526 ✗ nlsKinsolInplaceScaleX(kinsolData, x);
527 ✗ nlsKinsolInplaceScaleF(kinsolData, f);
528 };
529
530 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { iflag = 0 /* success */; }
531
532 #ifndef OMC_EMCC
533 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
534 #endif
535
536 ✗ return iflag;
537 }
538
539 /**
540 * @brief Calculate dense Jacobian matrix.
541 *
542 * @param N Size of vecX and vecFX.
543 * @param vecX Vector x.
544 * @param vecFX Residual vector f(x).
545 * @param Jac Dense Jacobian matrix J(x).
546 * @param kinsolUserData Pointer to Kinsol user data.
547 * @param tmp1 Unused, only to match interface of KINLsJacFn
548 * @param tmp2 Unused, only to match interface of KINLsJacFn
549 * @return int Return 0 on success, -1 on failure.
550 */
551 ✗ static int B_nlsDenseJac(long int N,
552 N_Vector vecX,
553 N_Vector vecFX,
554 SUNMatrix Jac,
555 NLS_USERDATA *kinsolUserData,
556 N_Vector tmp1,
557 N_Vector tmp2) {
558 DATA *data = kinsolUserData->data;
559 threadData_t *threadData = kinsolUserData->threadData;
560 ✗ NONLINEAR_SYSTEM_DATA *nlsData = kinsolUserData->nlsData;
561 ✗ B_NLS_KINSOL_DATA *kinsolData = (B_NLS_KINSOL_DATA *)nlsData->solverData;
562
563 ✗ if (SUNMatGetID(Jac) != SUNMATRIX_DENSE) {
564 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
565 "experimental-kinsol: B_nlsDenseJac illegal input Jac. Matrix is not dense!");
566 ✗ return -1;
567 }
568
569 /* prepare variables */
570 ✗ double *x = N_VGetArrayPointer(vecX);
571 ✗ double *fx = N_VGetArrayPointer(vecFX);
572 ✗ double *fRes = NV_DATA_S(kinsolData->fRes);
573 double xsave, xscale, sign;
574 double delta_hh;
575 const double delta_h = sqrt(DBL_EPSILON * 2e1);
576
577 long int col, row;
578
579 ✗ modelica_boolean stored_nominal_jac = kinsolData->useScaling;
580 ✗ if (kinsolData->useScaling) {
581 ✗ nlsKinsolInplaceUnscaleX(kinsolData, vecX);
582 ✗ kinsolData->useScaling = FALSE;
583 }
584
585 ✗ SUNMatZero_Dense(Jac);
586 ✗ B_nlsKinsolResiduals(vecX, vecFX, kinsolUserData);
587
588 /* performance measurement */
589 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
590
591 /* Use forward difference quotient to approximate Jacobian */
592 ✗ for (col = 0; col < N; col++) {
593 ✗ xsave = x[col];
594 ✗ delta_hh = delta_h * (fabs(xsave) + 1.0);
595 ✗ if ((xsave + delta_hh >= nlsData->max[col])) {
596 ✗ delta_hh *= -1.0;
597 }
598 ✗ x[col] += delta_hh;
599
600 /* Evaluate Jacobian function */
601 ✗ B_nlsKinsolResiduals(vecX, kinsolData->fRes, kinsolUserData);
602
603 /* Calculate scaled difference quotient */
604 ✗ delta_hh = 1.0 / delta_hh;
605
606 ✗ for (row = 0; row < N; row++) {
607 ✗ SM_ELEMENT_D(Jac, row, col) = (fRes[row] - fx[row]) * delta_hh;
608 }
609 ✗ x[col] = xsave;
610 }
611
612 ✗ kinsolData->useScaling = stored_nominal_jac;
613 ✗ if (kinsolData->useScaling) {
614 ✗ nlsKinsolInplaceScaleX(kinsolData, vecX);
615 ✗ nlsKinsolInplaceScaleJac(nlsData, kinsolData, Jac);
616 };
617
618 /* debug */
619 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC)) {
620 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "experimental-kinsol: Dense matrix (scaled = %s).", kinsolData->useScaling ? "true" : "false");
621 ✗ SUNDenseMatrix_Print(Jac, stdout); /* TODO: Print in OMC_LOG_NLS_JAC */
622 ✗ B_nlsKinsolJacSumDense(Jac);
623 ✗ messageClose(OMC_LOG_NLS_JAC);
624 }
625
626 /* performance measurement and statistics */
627 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
628 ✗ nlsData->numberOfJEval++;
629
630 ✗ return 0;
631 }
632
633 /**
634 * @brief Finish sparse matrix by fixing colprts.
635 *
636 * Last value of indexptrs should always be nnz.
637 * Search for empty rows which would mean the matrix is singular.
638 *
639 * @param A CSC matrix
640 */
641 ✗ static void finishSparseColPtr(SUNMatrix A, int nnz) {
642 int i;
643
644 /* TODO: Remove this check for performance reasons? */
645 ✗ if (SM_SPARSETYPE_S(A) != SUN_CSC_MAT) {
646 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
647 "experimental-kinsol: In function finishSparseColPtr: Wrong sparse format of SUNMatrix A.");
648 }
649
650 /* Set last value of indexptrs to nnz */
651 ✗ SM_INDEXPTRS_S(A)[SM_COLUMNS_S(A)] = nnz;
652
653 /* Check for empty rows */
654 ✗ for (i = 1; i < SM_COLUMNS_S(A) + 1; ++i) {
655 ✗ if (SM_INDEXPTRS_S(A)[i] == SM_INDEXPTRS_S(A)[i - 1]) {
656 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0,
657 "experimental-kinsol: Jacobian row %d singular. See OMC_LOG_NLS for "
658 "more information.",
659 i);
660 ✗ SM_INDEXPTRS_S(A)[i] = SM_INDEXPTRS_S(A)[i - 1];
661 }
662 }
663 ✗ }
664
665 /**
666 * @brief Perform derivative test comparing symbolic and numerical Jacobians for KINSOL
667 *
668 * Compares the symbolic Jacobian (sparse CSC format) with a numerically approximated
669 * dense Jacobian, checking for numerical and structural anomalies. The numerical
670 * Jacobian is computed using finite differences via B_nlsDenseJac.
671 *
672 * @param data Runtime data structure
673 * @param nlsData Nonlinear system data
674 * @param kinsolData KINSOL solver data structure
675 * @param Jsym Symbolic Jacobian in sparse CSC format
676 * @param tol Tolerance, all relative errors above tol are considered anomalies
677 * @param newJac TRUE if called from jacobian evaluation, FALSE if called from solver entry point
678 *
679 * @return int 1 derivative test failed and no error
680 * 0 derivative test successful and no error
681 * -1 internal error
682 */
683 ✗ static int nlsKinsolDenseDerivativeTest(DATA *data, NONLINEAR_SYSTEM_DATA *nlsData, B_NLS_KINSOL_DATA *kinsolData,
684 SUNMatrix Jsym, SolverCaller caller)
685 {
686 int row, col, nz, numericalErrorCount, structuralErrorCount;
687 ✗ const int size = nlsData->size;
688 int ret = 0;
689
690 modelica_real symValue, numValue, absError, relError;
691 modelica_real maxError = 0.0;
692
693 modelica_boolean errorFound;
694
695 ✗ sunindextype nnz = SUNSparseMatrix_NNZ(Jsym);
696 ✗ sunindextype columns = SUNSparseMatrix_Columns(Jsym);
697 ✗ sunindextype rows = SUNSparseMatrix_Rows(Jsym);
698
699 ✗ sunindextype *colPointers = SM_INDEXPTRS_S(Jsym);
700 ✗ sunindextype *rowIndices = SM_INDEXVALS_S(Jsym);
701 ✗ sunrealtype *symValues = SM_DATA_S(Jsym);
702
703 // allocate temporary memory for dense finite-diff matrix
704 ✗ N_Vector vecX = N_VNew_Serial(size, kinsolData->sunctx);
705 ✗ N_Vector vecFX = N_VNew_Serial(size, kinsolData->sunctx);
706 ✗ N_Vector tmp1 = N_VNew_Serial(size, kinsolData->sunctx);
707 ✗ N_Vector tmp2 = N_VNew_Serial(size, kinsolData->sunctx);
708 ✗ SUNMatrix Jnum = SUNDenseMatrix(size, size, kinsolData->sunctx);
709
710 // set tolerances
711 ✗ modelica_real Atol = omc_flag[FLAG_NLS_JAC_TEST_ATOL] ? atof(omc_flagValue[FLAG_NLS_JAC_TEST_ATOL]) : 100 * DBL_EPSILON;
712 ✗ modelica_real Rtol = omc_flag[FLAG_NLS_JAC_TEST_RTOL] ? atof(omc_flagValue[FLAG_NLS_JAC_TEST_RTOL]) : 1e-4;
713
714 if (kinsolData->useScaling) {
715 errorFound = FALSE;
716 }
717
718 // copy current x into new vector, compute f(x) and corresponding dense finite-diff Jacobian
719 ✗ SUNMatZero(Jnum);
720 ✗ N_VScale(1.0, kinsolData->initialGuess, vecX);
721 ✗ B_nlsKinsolResiduals(vecX, vecFX, kinsolData->userData);
722 ✗ if (B_nlsDenseJac(size, vecX, vecFX, Jnum, kinsolData->userData, tmp1, tmp2) != 0)
723 {
724 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Numerical Jacobian computation failed in nlsKinsolDenseDerivativeTest");
725 ret = -1;
726 ✗ SUNMatDestroy(Jnum);
727 ✗ N_VDestroy_Serial(vecX);
728 ✗ N_VDestroy_Serial(vecFX);
729 ✗ N_VDestroy_Serial(tmp1);
730 ✗ N_VDestroy_Serial(tmp2);
731 ✗ return ret;
732 }
733
734 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "%s: Derivative test (atol=%.5e, rtol=%.5e, scaled = %s, Caller: %s):",
735 ✗ SolverCaller_callerString(caller), Atol, Rtol, kinsolData->useScaling ? "true" : "false", SolverCaller_toString(caller));
736 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Matrix Info");
737 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "NLS index = " OMC_INT_FORMAT, nlsData->equationIndex);
738 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Columns = " OMC_INT_FORMAT, columns);
739 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Rows = " OMC_INT_FORMAT, rows);
740 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "NNZ = " OMC_INT_FORMAT, nnz);
741 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Curr Time = %-11.5e", data->localData[0]->timeValue);
742
743 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
744
745 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Anomalies");
746
747 nz = 0;
748 numericalErrorCount = 0;
749 structuralErrorCount = 0;
750
751 ✗ for (col = 0; col < size; col++)
752 {
753 errorFound = FALSE;
754
755 ✗ for (row = 0; row < size; row++)
756 {
757 ✗ numValue = SM_ELEMENT_D(Jnum, row, col);
758
759 ✗ if (colPointers[col] <= nz && nz < colPointers[col+1] && rowIndices[nz] == row)
760 {
761 // structural non-zero -> compare values
762 ✗ symValue = symValues[nz++];
763 ✗ absError = fabs(symValue - numValue);
764 ✗ relError = (absError < Atol) ? 0.0 : absError / fmax(fabs(numValue), fabs(symValue));
765
766 ✗ if (relError > maxError)
767 {
768 maxError = relError;
769 }
770
771 ✗ if (relError > Rtol)
772 {
773 // tolerance exceeded -> numerical error
774 ✗ if (!errorFound)
775 {
776 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Column / Variable: %i, Name: %s",
777 ✗ col + 1, modelInfoGetEquation(&data->modelData->modelDataXml, nlsData->equationIndex).vars[col]);
778 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6s %-6s %-15s %-15s %-8s",
779 "Type", "Col", "Row", "Symbolic", "Numerical", "RelError");
780 errorFound = TRUE;
781 }
782 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6d %-6d %+15.8e %+15.8e %+13.8e",
783 "Numerical", col + 1, row + 1, symValue, numValue, relError);
784 ✗ numericalErrorCount++;
785 }
786 }
787 ✗ else if (fabs(numValue) > Atol)
788 {
789 // structural error with tolerance exceeded -> non-zero in numerical Jacobian but zero in symbolic
790 ✗ if (!errorFound)
791 {
792 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Column / Variable: %i, Name: %s",
793 ✗ col + 1, modelInfoGetEquation(&data->modelData->modelDataXml, nlsData->equationIndex).vars[col]);
794 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6s %-6s %-15s %-15s %-8s",
795 "Type", "Col", "Row", "Symbolic", "Numerical", "RelError");
796 errorFound = TRUE;
797 }
798 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6d %-6d %+15.8e %+15.8e %+13.8e",
799 "Structural", col + 1, row + 1, 0.0, numValue, 1.0);
800 ✗ structuralErrorCount++;
801 }
802 }
803
804 ✗ if (errorFound)
805 {
806 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
807 }
808 }
809 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
810
811 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Summary");
812 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Numerical errors: %d (value mismatch w.r.t. reference)", numericalErrorCount);
813 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Structural errors: %d (non-zero not in sparsity pattern)", structuralErrorCount);
814 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Max relative error: %.3e", maxError);
815
816 ✗ if (numericalErrorCount + structuralErrorCount > 0)
817 {
818 ✗ warningStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Derivative test failed (%d numerical, %d structural errors)",
819 numericalErrorCount, structuralErrorCount);
820 ret = 1;
821 }
822 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
823
824 ✗ SUNMatDestroy(Jnum);
825 ✗ N_VDestroy_Serial(vecX);
826 ✗ N_VDestroy_Serial(vecFX);
827 ✗ N_VDestroy_Serial(tmp1);
828 ✗ N_VDestroy_Serial(tmp2);
829
830 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
831
832 ✗ return ret;
833 }
834
835 /**
836 * @brief Computes symbolic Jacobian matrix Jac(vecX)
837 *
838 * @param vecX
839 * @param vecFX just for interface compatibility, will not be used here
840 * @param Jac Allocated Jacobian, contains symbolic Jacobian on exit
841 * @param userData Void pointer to user data of type NLS_USERDATA*.
842 * @param tmp1 Unused, only to match interface of KINLsJacFn
843 * @param tmp2 Unused, only to match interface of KINLsJacFn
844 * @return int
845 */
846 ✗ static int B_nlsSparseSymJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
847 void *userData, N_Vector tmp1, N_Vector tmp2) {
848 /* Variables */
849 NLS_USERDATA* kinsolUserData = (NLS_USERDATA *)userData;;
850 ✗ DATA* data = kinsolUserData->data;
851 ✗ threadData_t* threadData = kinsolUserData->threadData;
852 ✗ NONLINEAR_SYSTEM_DATA* nlsData = kinsolUserData->nlsData;
853 ✗ B_NLS_KINSOL_DATA* kinsolData = (B_NLS_KINSOL_DATA *)nlsData->solverData;
854 ✗ JACOBIAN* jacobian = kinsolUserData->analyticJacobian;
855 ✗ assertStreamPrint(threadData, NULL != jacobian, "jacobian is NULL");
856 ✗ const SPARSE_PATTERN* sp = jacobian->sparsePattern;
857 ✗ assertStreamPrint(threadData, NULL != sp, "sp is NULL");
858 long int column, nz;
859
860 ✗ if (SUNMatGetID(Jac) != SUNMATRIX_SPARSE || SM_SPARSETYPE_S(Jac) == SUN_CSR_MAT) {
861 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
862 "experimental-kinsol: B_nlsSparseSymJac illegal input Jac. Matrix is not sparse!");
863 ✗ return -1;
864 }
865
866 /* performance measurement */
867 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
868
869 ✗ if (kinsolData->useScaling) {
870 ✗ nlsKinsolInplaceUnscaleX(kinsolData, vecX);
871 }
872
873 /* call generic sparse Jacobian with CSC buffer "SM_DATA_S(Jac)" */
874 ✗ evalJacobian(data, threadData, jacobian, NULL, SM_DATA_S(Jac), FALSE);
875 ✗ setSundialsSparsePattern(jacobian, Jac);
876
877 /* Finish sparse matrix and do a cheap check for singularity */
878 ✗ finishSparseColPtr(Jac, sp->nnz);
879
880 ✗ if (kinsolData->useScaling) {
881 ✗ nlsKinsolInplaceScaleX(kinsolData, vecX);
882 ✗ nlsKinsolInplaceScaleJac(nlsData, kinsolData, Jac);
883 };
884
885 /* Debug print */
886 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC)) {
887 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "experimental-kinsol: Sparse Matrix.");
888 ✗ SUNSparseMatrix_Print(Jac, stdout); /* TODO: Print in OMC_LOG_NLS_JAC */
889 ✗ B_nlsKinsolJacSumSparse(Jac);
890 ✗ messageClose(OMC_LOG_NLS_JAC);
891 }
892
893 ✗ if (omc_useStream[OMC_LOG_NLS_DERIVATIVE_TEST])
894 {
895 ✗ nlsKinsolDenseDerivativeTest(data, nlsData, kinsolData, Jac, KINSOL_B_JAC_EVAL);
896 }
897
898 ✗ if (omc_useStream[OMC_LOG_NLS_JAC_SUMS])
899 {
900 ✗ nlsJacobianRowColSums(data, nlsData, Jac, KINSOL_B_JAC_EVAL /* called at evaluation */, kinsolData->useScaling /* scaled */);
901 }
902
903 /* performance measurement and statistics */
904 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
905 ✗ nlsData->numberOfJEval++;
906
907 ✗ return 0;
908 }
909
910 /**
911 * @brief Colored numeric Jacobian evaluation.
912 *
913 * Finite differences while using coloring of Jacobian.
914 * Jacobian matrix format has to be compressed sparse columns (CSC).
915 *
916 * @param vecX Input vector x.
917 * @param vecFX Vector for residual evaluation: f(x)
918 * @param Jac Jacobian to calculate: J(x)
919 * @param userData Pointer to user data, tpyecasted to `NLS_USERDATA`.
920 * @param tmp1 Work vector.
921 * @param tmp2 Work vector.
922 * @return int Return 0 on success.
923 */
924 ✗ static int B_nlsSparseJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
925 void *userData, N_Vector tmp1, N_Vector tmp2) {
926 /* Variables */
927 NLS_USERDATA *kinsolUserData;
928 DATA *data;
929 NONLINEAR_SYSTEM_DATA *nlsData;
930 B_NLS_KINSOL_DATA *kinsolData;
931 SPARSE_PATTERN *sparsePattern;
932
933 ✗ if (SUNMatGetID(Jac) != SUNMATRIX_SPARSE || SM_SPARSETYPE_S(Jac) == SUN_CSR_MAT) {
934 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
935 "experimental-kinsol: B_nlsSparseJac illegal input Jac. Matrix is not sparse!");
936 ✗ return -1;
937 }
938
939 double *x;
940 double *fx;
941 double *xsave;
942 double *delta_hh;
943 double *xScaling;
944 double *fRes;
945
946 const double delta_h = sqrt(DBL_EPSILON * 2e1);
947
948 modelica_real result;
949 long int i, j, ii;
950 int nth;
951
952 modelica_boolean stored_nominal_jac;
953
954 /* Access userData and nonlinear system data */
955 kinsolUserData = (NLS_USERDATA *)userData;
956 ✗ data = kinsolUserData->data;
957 ✗ nlsData = kinsolUserData->nlsData;
958 ✗ kinsolData = (B_NLS_KINSOL_DATA *)nlsData->solverData;
959 ✗ sparsePattern = nlsData->sparsePattern;
960
961 /* Access N_Vector variables */
962 ✗ x = N_VGetArrayPointer(vecX);
963 ✗ fx = N_VGetArrayPointer(vecFX);
964 ✗ xsave = N_VGetArrayPointer(tmp1);
965 ✗ delta_hh = N_VGetArrayPointer(tmp2);
966 ✗ fRes = N_VGetArrayPointer(kinsolData->fRes);
967
968 nth = 0;
969
970 /* performance measurement */
971 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
972
973 /* reset matrix */
974 ✗ SUNMatZero(Jac);
975
976 ✗ stored_nominal_jac = kinsolData->useScaling;
977 ✗ if (kinsolData->useScaling) {
978 ✗ nlsKinsolInplaceUnscaleX(kinsolData, vecX);
979 ✗ kinsolData->useScaling = FALSE;
980 }
981
982 ✗ B_nlsKinsolResiduals(vecX, vecFX, userData);
983
984 /* Approximate Jacobian */
985 ✗ for (i = 0; i < sparsePattern->maxColors; i++) {
986 ✗ for (ii = 0; ii < kinsolData->size; ii++) {
987 ✗ if (sparsePattern->colorCols[ii] - 1 == i) {
988 ✗ xsave[ii] = x[ii];
989 ✗ delta_hh[ii] = delta_h * (fabs(xsave[ii]) + 1.0);
990 ✗ if ((xsave[ii] + delta_hh[ii] >= nlsData->max[ii])) {
991 ✗ delta_hh[ii] *= -1;
992 }
993 ✗ x[ii] += delta_hh[ii];
994
995 /* Calculate scaled difference quotient */
996 ✗ delta_hh[ii] = 1. / delta_hh[ii];
997 }
998 }
999 /* Evaluate residual function */
1000 ✗ B_nlsKinsolResiduals(vecX, kinsolData->fRes, userData);
1001
1002 /* Save column in Jac and unset seed variables */
1003 ✗ for (ii = 0; ii < kinsolData->size; ii++) {
1004 ✗ if (sparsePattern->colorCols[ii] - 1 == i) {
1005 ✗ nth = sparsePattern->leadindex[ii];
1006 ✗ while (nth < sparsePattern->leadindex[ii + 1]) {
1007 ✗ j = sparsePattern->index[nth];
1008
1009 // TODO: investigate NaN values stemming from residual functions
1010 // Hypothesis: lambda = 0 system forces variables to be 0, while for lambda = eps, we divide by them?!
1011 ✗ result = (fRes[j] - fx[j]) * delta_hh[ii];
1012
1013 // (IN)SANITY CHECK
1014 ✗ if (isnan(result) || isinf(result)) {
1015 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0,
1016 "WARNING: NaN (%d) or Inf (%d) detected at col %ld row %ld: fRes=%g, fx=%g, delta_hh=%g, x=%g, xsave=%g\n"
1017 "ACTION: setting Jacobian entry := 0.0 and trying to recover...",
1018 ✗ isnan(result), isinf(result), ii, j, fRes[j], fx[j], delta_hh[ii], x[ii], xsave[ii]);
1019 result = 0.0;
1020 }
1021 ✗ setJacElementSundialsSparse(j, ii, nth, result, Jac, SM_CONTENT_S(Jac)->M);
1022 ✗ nth++;
1023 }
1024 ✗ x[ii] = xsave[ii];
1025 }
1026 }
1027 }
1028 /* Finish sparse matrix */
1029 ✗ setSundialsSparseColPtrs(sparsePattern, Jac);
1030 ✗ finishSparseColPtr(Jac, sparsePattern->nnz);
1031
1032 ✗ kinsolData->useScaling = stored_nominal_jac;
1033 ✗ if (kinsolData->useScaling) {
1034 ✗ nlsKinsolInplaceScaleX(kinsolData, vecX);
1035 ✗ nlsKinsolInplaceScaleF(kinsolData, vecFX);
1036 ✗ nlsKinsolInplaceScaleJac(nlsData, kinsolData, Jac);
1037 };
1038
1039 /* Debug print */
1040 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC)) {
1041 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "experimental-kinsol: Sparse Matrix.");
1042 ✗ SUNSparseMatrix_Print(Jac, stdout);
1043 ✗ B_nlsKinsolJacSumSparse(Jac);
1044 ✗ messageClose(OMC_LOG_NLS_JAC);
1045 }
1046
1047 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_DEBUG)) {
1048 ✗ sundialsPrintSparseMatrix(Jac, "A", OMC_LOG_JAC);
1049 }
1050
1051 ✗ if (omc_useStream[OMC_LOG_NLS_DERIVATIVE_TEST])
1052 {
1053 ✗ nlsKinsolDenseDerivativeTest(data, nlsData, kinsolData, Jac, KINSOL_B_JAC_EVAL);
1054 }
1055
1056 ✗ if (omc_useStream[OMC_LOG_NLS_JAC_SUMS])
1057 {
1058 ✗ nlsJacobianRowColSums(data, nlsData, Jac, KINSOL_B_JAC_EVAL /* called at evaluation */, kinsolData->useScaling /* scaled */);
1059 }
1060
1061 /* performance measurement and statistics */
1062 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
1063 ✗ nlsData->numberOfJEval++;
1064
1065 ✗ return 0;
1066 }
1067
1068 /**
1069 * @brief Check for zero columns of matrix and print absolute sums.
1070 *
1071 * Compute absolute sum for each column and print the result.
1072 * Report a warning if it is zero, since the matrix is singular in that case.
1073 *
1074 * @param A Dense matrix stored columnwise
1075 */
1076 ✗ static void B_nlsKinsolJacSumDense(SUNMatrix A) {
1077 /* Variables */
1078 int i, j;
1079 double sum;
1080
1081 ✗ for (i = 0; i < SM_ROWS_D(A); ++i) {
1082 sum = 0.0;
1083 ✗ for (j = 0; j < SM_COLUMNS_D(A); ++j) {
1084 ✗ sum += fabs(SM_ELEMENT_D(A, j, i));
1085 }
1086
1087 ✗ if (sum == 0.0) { /* TODO: Don't check for equality(!), maybe use DBL_EPSILON */
1088 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1089 "experimental-kinsol: Column %d of Jacobian is zero. Jacobian is singular.",
1090 i);
1091 } else {
1092 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "Column %d of Jacobian absolute sum = %g",
1093 i, sum);
1094 }
1095 }
1096 ✗ }
1097
1098 /**
1099 * @brief Check for zero columns of matrix and print absolute sums.
1100 *
1101 * Compute absolute sum for each column and print the result.
1102 * Report a warning if it is zero, since the matrix is singular in that case.
1103 *
1104 * @param A CSC matrix
1105 */
1106 ✗ static void B_nlsKinsolJacSumSparse(SUNMatrix A) {
1107 /* Variables */
1108 int i, j;
1109 double sum;
1110
1111 /* Check format of A */
1112 ✗ if (SM_SPARSETYPE_S(A) != SUN_CSC_MAT) {
1113 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1114 "experimental-kinsol: In function B_nlsKinsolJacSumSparse: Wrong sparse format "
1115 "of SUNMatrix A.");
1116 }
1117
1118 /* Check sums of each column of A */
1119 ✗ for (i = 0; i < SM_COLUMNS_S(A); ++i) {
1120 sum = 0.0;
1121 ✗ for (j = SM_INDEXPTRS_S(A)[i]; j < SM_INDEXPTRS_S(A)[i + 1]; ++j) {
1122 ✗ sum += fabs(SM_DATA_S(A)[j]);
1123 }
1124
1125 ✗ if (sum == 0.0) { /* TODO: Don't check for equality(!), maybe use DBL_EPSILON */
1126 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1127 "experimental-kinsol: Column %d of Jacobian is zero. Jacobian is singular.",
1128 i);
1129 } else {
1130 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "Column %d of Jacobian absolute sum = %g",
1131 i, sum);
1132 }
1133 }
1134 ✗ }
1135
1136 /**
1137 * @brief Set maximum scaled length of Newton step.
1138 *
1139 * Will be set to the weighted Euclidean l_2 norm of xScale with maxstepfactor
1140 * as weights. maxStep = sqrt(sum_{1=0}^{n-1} (xScale[i]*maxstepfactor)^2)
1141 *
1142 * @param kinsolData
1143 * @param maxstepfactor
1144 */
1145 ✗ static void B_nlsKinsolSetMaxNewtonStep(B_NLS_KINSOL_DATA *kinsolData,
1146 double maxstepfactor) {
1147 /* Variables */
1148 int flag;
1149
1150 ✗ N_VConst(maxstepfactor, kinsolData->fTmp);
1151 ✗ kinsolData->mxnstepin = N_VWL2Norm(kinsolData->xScale, kinsolData->fTmp); // TODO: ?
1152
1153 /* Set maximum step size */
1154 ✗ flag = KINSetMaxNewtonStep(kinsolData->kinsolMemory, kinsolData->mxnstepin);
1155 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetMaxNewtonStep");
1156 ✗ }
1157
1158 /**
1159 * @brief Set initial guess for KINSOL and unscale previous Jacobian
1160 *
1161 * Depending on mode extrapolate start value or use old value for
1162 * initialization.
1163 *
1164 * @param data
1165 * @param kinsolData
1166 * @param nlsData
1167 * @param mode Has to be `INITIAL_EXTRAPOLATION` for extrapolation or
1168 * `INITIAL_OLDVALUES` for using old values.
1169 */
1170 ✗ static void B_nlsKinsolResetInitialUnscaled(DATA *data, B_NLS_KINSOL_DATA *kinsolData,
1171 NONLINEAR_SYSTEM_DATA *nlsData,
1172 B_initialMode mode) {
1173 ✗ double *xStart = NV_DATA_S(kinsolData->initialGuess);
1174
1175 /* Set x vector */
1176 ✗ switch (mode) {
1177 ✗ case B_INITIAL_EXTRAPOLATION:
1178 ✗ if (data->simulationInfo->discreteCall) {
1179 ✗ memcpy(xStart, nlsData->nlsx, nlsData->size * (sizeof(double)));
1180 } else {
1181 ✗ memcpy(xStart, nlsData->nlsxExtrapolation,
1182 ✗ nlsData->size * (sizeof(double)));
1183 }
1184 break;
1185 ✗ case B_INITIAL_OLDVALUES:
1186 ✗ memcpy(xStart, nlsData->nlsxOld, nlsData->size * (sizeof(double)));
1187 break;
1188 ✗ default:
1189 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1190 "experimental-kinsol: Function B_nlsKinsolResetInitialUnscaled: Unknown mode %d.",
1191 (int)mode);
1192 }
1193 ✗ }
1194
1195 /**
1196 * @brief Scale x vector.
1197 *
1198 * Scale with 1.0 for mode `SCALING_ONES`.
1199 * Scale with 1/fmax(nominal,|xStart|) for mode `SCALING_NOMINALSTART`.
1200 *
1201 * @param data unused
1202 * @param kinsolData
1203 * @param nlsData
1204 * @param mode Mode for scaling. Use `SCALING_NOMINALSTART` for nominal
1205 * scaling and `SCALING_ONES` for no scaling. Will be
1206 * overwritten by simulation flag `FLAG_NO_SCALING`.
1207 */
1208 ✗ static void B_nlsKinsolXScaling(DATA *data, B_NLS_KINSOL_DATA *kinsolData,
1209 NONLINEAR_SYSTEM_DATA *nlsData,
1210 B_scalingMode mode) {
1211 ✗ double *xStart = NV_DATA_S(kinsolData->initialGuess);
1212 ✗ double *xScaling = NV_DATA_S(kinsolData->xScale);
1213 int i;
1214
1215 /* if noScaling flag is used overwrite mode */
1216 ✗ if (omc_flag[FLAG_NO_SCALING]) {
1217 mode = B_SCALING_ONES;
1218 }
1219
1220 /* Use nominal value or the actual working point for scaling */
1221 ✗ switch (mode) {
1222 case B_SCALING_NOMINALSTART:
1223 ✗ for (i = 0; i < nlsData->size; i++) {
1224 ✗ xScaling[i] = 1.0 / fmax(nlsData->nominal[i], fabs(xStart[i]));
1225 }
1226 break;
1227 case B_SCALING_ONES:
1228 ✗ for (i = 0; i < nlsData->size; i++) {
1229 ✗ xScaling[i] = 1.0;
1230 }
1231 break;
1232 ✗ case B_SCALING_JACOBIAN:
1233 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1234 "experimental-kinsol: Function B_nlsKinsolXScaling: Invalid mode SCALING_JACOBIAN.");
1235 ✗ default:
1236 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1237 "experimental-kinsol: Function B_nlsKinsolXScaling: Unknown mode %d.", (int)mode);
1238 }
1239 ✗ }
1240
1241 /**
1242 * @brief Scale f(x) vector.
1243 *
1244 * @param data
1245 * @param kinsolData
1246 * @param nlsData
1247 * @param mode
1248 */
1249 ✗ static void B_nlsKinsolFScaling(DATA *data, B_NLS_KINSOL_DATA *kinsolData,
1250 NONLINEAR_SYSTEM_DATA *nlsData,
1251 B_scalingMode mode) {
1252 ✗ double *fScaling = NV_DATA_S(kinsolData->fScale);
1253 ✗ N_Vector x = kinsolData->initialGuess;
1254
1255 int i, j;
1256 SUNErrCode ret;
1257
1258 /* If noScaling flag is used overwrite mode */
1259 ✗ if (omc_flag[FLAG_NO_SCALING]) {
1260 mode = B_SCALING_ONES;
1261 }
1262
1263 /* Disable scaled jacobian evaluation */
1264 ✗ kinsolData->useScaling = FALSE;
1265
1266 /* Use nominal value or the actual working point for scaling */
1267 ✗ switch (mode) {
1268 ✗ case B_SCALING_JACOBIAN:
1269
1270 /* Calculate the scaled Jacobian */
1271 ✗ if (nlsData->sparsePattern && kinsolData->linearSolverMethod == NLS_LS_KLU) {
1272 ✗ if (kinsolData->solved != NLS_SOLVED) {
1273 ✗ if (nlsData->analyticalJacobianColumn != NULL) {
1274 /* Calculate the sparse Jacobian symbolically */
1275 ✗ B_nlsSparseSymJac(x, kinsolData->fTmp, kinsolData->J, kinsolData->userData, NULL, NULL);
1276 } else {
1277 /* Update f(x) for the numerical jacobian matrix */
1278 ✗ B_nlsKinsolResiduals(x, kinsolData->fTmp, kinsolData->userData);
1279 ✗ B_nlsSparseJac(x, kinsolData->fTmp, kinsolData->J, kinsolData->userData, kinsolData->tmp1, kinsolData->tmp2);
1280 }
1281 }
1282 /* Scale the current Jacobian */
1283 ✗ SUNMatCopy_Sparse(kinsolData->J, kinsolData->scaledJ); /* Copy J into scaledJ */
1284 ✗ ret = _omc_SUNSparseMatrixVecScaling(kinsolData->scaledJ, kinsolData->xScale);
1285 ✗ if (ret != SUN_SUCCESS) {
1286 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "experimental-kinsol: _omc_SUNSparseMatrixVecScaling failed.");
1287 }
1288 } else {
1289 /* Update f(x) for the numerical jacobian matrix */
1290 ✗ B_nlsKinsolResiduals(x, kinsolData->fTmp, kinsolData->userData);
1291 ✗ B_nlsDenseJac(nlsData->size, x, kinsolData->fTmp, kinsolData->J,
1292 kinsolData->userData, NULL, NULL);
1293 }
1294
1295 ✗ for (i = 0; i < nlsData->size; i++) {
1296 ✗ fScaling[i] = 1e-12;
1297 }
1298
1299 ✗ switch (SUNMatGetID(kinsolData->J))
1300 {
1301 case SUNMATRIX_SPARSE:
1302 ✗ for (i = 0; i < SM_NNZ_S(kinsolData->scaledJ); ++i) {
1303 ✗ if (fScaling[SM_INDEXVALS_S(kinsolData->scaledJ)[i]] < fabs(SM_DATA_S(kinsolData->scaledJ)[i])) {
1304 ✗ fScaling[SM_INDEXVALS_S(kinsolData->scaledJ)[i]] = fabs(SM_DATA_S(kinsolData->scaledJ)[i]);
1305 }
1306 }
1307 break;
1308 case SUNMATRIX_DENSE:
1309 ✗ for (i = 0; i < nlsData->size; i++) {
1310 ✗ for (j = 0; j < nlsData->size; j++) {
1311 ✗ if (fScaling[i] < fabs(SM_ELEMENT_D(kinsolData->J, j, i))) {
1312 ✗ fScaling[i] = fabs(SM_ELEMENT_D(kinsolData->J, j, i));
1313 }
1314 }
1315 }
1316 break;
1317 ✗ default:
1318 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1319 "KINSOL: Function B_nlsKinsolFScaling: Unknown matrix type.");
1320 }
1321
1322 /* inverse fScale */
1323 ✗ N_VInv(kinsolData->fScale, kinsolData->fScale);
1324
1325 ✗ break;
1326 case B_SCALING_ONES:
1327 ✗ for (i = 0; i < nlsData->size; i++) {
1328 ✗ fScaling[i] = 1.0;
1329 }
1330 break;
1331 ✗ case B_SCALING_NOMINALSTART:
1332 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1333 "experimental-kinsol: Function B_nlsKinsolFScaling: Invalid mode SCALING_NOMINALSTART.");
1334 ✗ default:
1335 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1336 "experimental-kinsol: Function B_nlsKinsolFScaling: Unknown mode %d.", (int)mode);
1337 }
1338 ✗ }
1339
1340 /**
1341 * @brief Print KINSOL configuration.
1342 *
1343 * Only prints if stream `LOG_NLS_V` is active.
1344 *
1345 * @param kinsolData
1346 * @param nlsData
1347 */
1348 ✗ static void B_nlsKinsolConfigPrint(B_NLS_KINSOL_DATA *kinsolData,
1349 NONLINEAR_SYSTEM_DATA *nlsData) {
1350 int retValue;
1351 double fNorm;
1352 ✗ DATA *data = kinsolData->userData->data;
1353 ✗ int eqSystemNumber = nlsData->equationIndex;
1354 _omc_vector vecStart, vecXScaling, vecFScaling;
1355
1356 ✗ if (!omc_useStream[OMC_LOG_NLS_V]) {
1357 ✗ return;
1358 }
1359
1360 ✗ _omc_initVector(&vecStart, kinsolData->size,
1361 ✗ NV_DATA_S(kinsolData->initialGuess));
1362 _omc_initVector(&vecXScaling, kinsolData->size,
1363 ✗ NV_DATA_S(kinsolData->xScale));
1364 _omc_initVector(&vecFScaling, kinsolData->size,
1365 ✗ NV_DATA_S(kinsolData->fScale));
1366
1367 ✗ if (eqSystemNumber>0) {
1368 ✗ _omc_printVectorWithEquationInfo(
1369 &vecStart, "Initial guess values", OMC_LOG_NLS_V,
1370 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber));
1371
1372 ✗ _omc_printVectorWithEquationInfo(
1373 &vecXScaling, "xScaling", OMC_LOG_NLS_V,
1374 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber));
1375 }
1376
1377 ✗ _omc_printVector(&vecFScaling, "fScaling", OMC_LOG_NLS_V);
1378
1379 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol F tolerance: %g", kinsolData->fnormtol);
1380 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol minimal step size %g",
1381 kinsolData->scsteptol);
1382 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol max iterations %d",
1383 ✗ 20 * kinsolData->size);
1384 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol strategy %d",
1385 kinsolData->kinsolStrategy);
1386 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol current retry %d", kinsolData->retries);
1387 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol max step %g", kinsolData->mxnstepin);
1388 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol linear solver %d",
1389 ✗ kinsolData->linearSolverMethod);
1390 }
1391
1392 /**
1393 * @brief Try to handle errors of KINSol().
1394 *
1395 * @param errorCode Error code from KINSOL.
1396 * @param data Pointer to data struct.
1397 * @param nlsData Non-linear solver data.
1398 * @param kinsolData Kinsol data.
1399 * @return modelica_boolean Return true, if it is possible to retry KINSol().
1400 */
1401 ✗ static modelica_boolean nlsKinsolErrorHandler(int errorCode, DATA *data,
1402 NONLINEAR_SYSTEM_DATA *nlsData,
1403 B_NLS_KINSOL_DATA *kinsolData) {
1404 int flag; /* KIN_* and KINLS_* codes, which are plain macros */
1405 SUNErrCode sunFlag; /* SUNLinearSolver codes, which are not */
1406 double fNorm;
1407 double *xStart = NV_DATA_S(kinsolData->initialGuess);
1408 long outL;
1409
1410 ✗ flag = KINSetNoInitSetup(kinsolData->kinsolMemory, SUNFALSE);
1411 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNoInitSetup");
1412
1413 ✗ switch (errorCode) {
1414 ✗ case KIN_MEM_NULL:
1415 ✗ throwStreamPrint(NULL, "experimental-kinsol: Memory NULL ERROR %d\n", errorCode);
1416 return FALSE;
1417 break;
1418 ✗ case KIN_ILL_INPUT:
1419 ✗ throwStreamPrint(NULL, "experimental-kinsol: Ill input ERROR %d\n", errorCode);
1420 return FALSE;
1421 break;
1422 ✗ case KIN_NO_MALLOC:
1423 ✗ throwStreamPrint(NULL, "experimental-kinsol: Memory issue ERROR %d\n", errorCode);
1424 return FALSE;
1425 break;
1426 /* Just retry with new initial guess */
1427 ✗ case KIN_MXNEWT_5X_EXCEEDED:
1428 ✗ warningStreamPrint(
1429 OMC_LOG_NLS_V, 0,
1430 "Newton step exceed the maximum step size several times. Try again "
1431 "after increasing maximum step size.\n");
1432 ✗ kinsolData->maxstepfactor *= 1e5;
1433 ✗ B_nlsKinsolSetMaxNewtonStep(kinsolData, kinsolData->maxstepfactor);
1434 ✗ return TRUE;
1435 break;
1436 /* Just retry without line search */
1437 ✗ case KIN_LINESEARCH_NONCONV:
1438 ✗ warningStreamPrint(
1439 OMC_LOG_NLS_V, 0,
1440 "kinsols line search did not convergence. Try without.\n");
1441 ✗ kinsolData->kinsolStrategy = KIN_NONE;
1442 ✗ kinsolData->retries--;
1443 ✗ return TRUE;
1444 break;
1445 /* Maybe happened because of an out-dated factorization, so just retry */
1446 ✗ case KIN_LSOLVE_FAIL:
1447 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1448 "experimental-kinsol: Matrix need new factorization. Try again.\n");
1449 ✗ if (kinsolData->linearSolverMethod == NLS_LS_KLU &&
1450 ✗ nlsData->sparsePattern) {
1451 /* Complete symbolic and numeric factorizations */
1452 ✗ sunFlag = SUNLinSol_KLUReInit(kinsolData->linSol, kinsolData->J,
1453 ✗ kinsolData->nnz, SUNKLU_REINIT_PARTIAL);
1454 ✗ checkReturnFlag_SUNDIALS(sunFlag, SUNDIALS_SUNLS_FLAG, "SUNLinSol_KLUReInit");
1455 ✗ return TRUE;
1456 }
1457 break;
1458 ✗ case KIN_MAXITER_REACHED:
1459 case KIN_REPTD_SYSFUNC_ERR:
1460 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1461 "experimental-kinsol: Runs into issues retry with different configuration.\n");
1462 ✗ break;
1463 ✗ case KIN_LINIT_FAIL:
1464 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1465 "experimental-kinsol: The linear solver's initialization function failed.\n");
1466 ✗ return errorCode;
1467 ✗ case KIN_LSETUP_FAIL:
1468 /* In case something goes wrong with the symbolic jacobian try the numerical */
1469 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1470 "experimental-kinsol: The kinls setup routine (lsetup) encountered an error. "
1471 "Retry with numerical Jacobian.\n");
1472 ✗ if (kinsolData->linearSolverMethod == NLS_LS_KLU) {
1473 ✗ if (nlsData->sparsePattern && nlsData->analyticalJacobianColumn != NULL) {
1474 ✗ flag = KINSetJacFn(kinsolData->kinsolMemory, B_nlsSparseJac);
1475 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KINLS_FLAG, "KINSetJacFn");
1476 ✗ if (flag < 0) {
1477 return FALSE;
1478 }
1479 } else {
1480 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "experimental-kinsol: Trying to switch to numeric Jacobian for sparse solver KLU, but no sparsity pattern is available.");
1481 ✗ return FALSE;
1482 }
1483 }
1484 break;
1485 ✗ case KIN_LINESEARCH_BCFAIL:
1486 ✗ KINGetNumBetaCondFails(kinsolData->kinsolMemory, &outL);
1487 ✗ warningStreamPrint(
1488 OMC_LOG_NLS_V, 0,
1489 "kinsols runs into issues with beta-condition fails: %ld\n", outL);
1490 ✗ break;
1491 ✗ default:
1492 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1493 "kinsol has a serious solving issue ERROR %d\n",
1494 errorCode);
1495 ✗ return FALSE;
1496 break;
1497 }
1498
1499 // TODO: configure the retry strategies properly!!
1500 // currently this is does not make sense with the new scaling
1501
1502 /* check if the current solution is sufficient anyway */
1503 ✗ KINGetFuncNorm(kinsolData->kinsolMemory, &fNorm);
1504 ✗ if (fNorm < B_FTOL_WITH_LESS_ACCURACY) {
1505 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol: Move forward with a less accurate solution.");
1506 ✗ KINSetFuncNormTol(kinsolData->kinsolMemory, B_FTOL_WITH_LESS_ACCURACY);
1507 ✗ KINSetScaledStepTol(kinsolData->kinsolMemory, B_FTOL_WITH_LESS_ACCURACY);
1508 ✗ kinsolData->resetTol = TRUE;
1509 ✗ return TRUE;
1510 } else {
1511 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "experimental-kinsol: Current status of fx = %f", fNorm);
1512 }
1513
1514 /* reconfigure kinsol for another try */
1515 ✗ switch (kinsolData->retries) {
1516 ✗ case 0:
1517 /* try without scaling */
1518 ✗ B_nlsKinsolXScaling(data, kinsolData, nlsData, B_SCALING_ONES);
1519 ✗ B_nlsKinsolFScaling(data, kinsolData, nlsData, B_SCALING_ONES);
1520 ✗ break;
1521 ✗ case 1:
1522 /* try without line-search and oldValues */
1523 ✗ B_nlsKinsolResetInitialUnscaled(data, kinsolData, nlsData, B_INITIAL_OLDVALUES);
1524 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH;
1525 ✗ break;
1526 ✗ case 2:
1527 /* try without line-search and oldValues */
1528 ✗ B_nlsKinsolResetInitialUnscaled(data, kinsolData, nlsData, B_INITIAL_EXTRAPOLATION);
1529 ✗ kinsolData->kinsolStrategy = KIN_NONE;
1530 ✗ break;
1531 ✗ case 3:
1532 /* try with exact newton */
1533 ✗ B_nlsKinsolXScaling(data, kinsolData, nlsData, B_SCALING_NOMINALSTART);
1534 ✗ B_nlsKinsolFScaling(data, kinsolData, nlsData, B_SCALING_JACOBIAN);
1535 ✗ B_nlsKinsolResetInitialUnscaled(data, kinsolData, nlsData, B_INITIAL_EXTRAPOLATION);
1536 ✗ KINSetMaxSetupCalls(kinsolData->kinsolMemory, 1);
1537 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH;
1538 ✗ break;
1539 ✗ case 4:
1540 /* try with exact newton to with out x scaling values */
1541 ✗ B_nlsKinsolXScaling(data, kinsolData, nlsData, B_SCALING_ONES);
1542 ✗ B_nlsKinsolFScaling(data, kinsolData, nlsData, B_SCALING_ONES);
1543 ✗ B_nlsKinsolResetInitialUnscaled(data, kinsolData, nlsData, B_INITIAL_OLDVALUES);
1544 ✗ KINSetMaxSetupCalls(kinsolData->kinsolMemory, 1);
1545 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH;
1546 ✗ break;
1547 default:
1548 /* Too many retries */
1549 return FALSE;
1550 break;
1551 }
1552
1553 return TRUE;
1554 }
1555
1556 /**
1557 * @brief Function dedicated for flag '--saveInitialGuess_system=/path/to/file.mat,nls_index' that computes only the
1558 * torn part of a given nonlinear system with index nls_index and writes it to a specified path. Will abort the simulation
1559 * if used incorrectly or if the file has been written.
1560 *
1561 * @param data Runtime data struct.
1562 * @param threadData Thread data for error handling.
1563 * @param nlsData Pointer to non-linear system data.
1564 * @return NLS_SOLVER_STATUS Return NLS_SOLVED on success and NLS_FAILED otherwise.
1565 */
1566 ✗ static void B_save_initial_guess_system(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA *nlsData)
1567 {
1568 char buf[512];
1569 ✗ strncpy(buf, omc_flagValue[FLAG_SAVE_INITIAL_GUESS_SYSTEM], sizeof(buf));
1570 ✗ buf[sizeof(buf)-1] = '\0';
1571
1572 ✗ char *comma = strchr(buf, ',');
1573 ✗ if (!comma)
1574 {
1575 ✗ throwStreamPrint(threadData, "Error: Invalid format for --saveInitialGuess_system flag - Expected: '--saveInitialGuess_system=/path/to/file.mat,nls_index')");
1576 }
1577 ✗ *comma = '\0';
1578
1579 const char *path = buf;
1580 ✗ int nls_idx = atoi(comma + 1);
1581 ✗ if (nls_idx < 0 && strcmp(comma + 1, "0") != 0)
1582 {
1583 ✗ throwStreamPrint(threadData, "Error: Invalid format for --saveInitialGuess_system flag - Expected: '--saveInitialGuess_system=/path/to/file.mat,nls_index')");
1584 }
1585
1586 ✗ B_NLS_KINSOL_DATA *kinsolData = (B_NLS_KINSOL_DATA *) nlsData->solverData;
1587 ✗ if (nlsData->equationIndex == nls_idx)
1588 {
1589 /* compute the residuals (will compute the torn part of the system and write them to localData[0]->realVars)
1590 * this function should have been called before as long as scaling is used (as we require the residual for the scalings),
1591 * but if no scaling is involved we must compute the residual here
1592 * TODO: replace this with call to list of inner equations */
1593 ✗ B_nlsKinsolResiduals(kinsolData->initialGuess, kinsolData->fTmp, kinsolData->userData);
1594
1595 /* write out current localData[0] to .mat file */
1596 simulation_result file_result;
1597 ✗ const char *outputFormat = data->simulationInfo->outputFormat;
1598
1599 ✗ file_result.filename = path;
1600 ✗ file_result.numpoints = 2;
1601 ✗ file_result.cpuTime = 0;
1602
1603 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Trying to write write initial guess for NLS system with index %d to file %s.\n", nls_idx, path);
1604
1605 ✗ data->simulationInfo->outputFormat = "mat";
1606 ✗ rust_result_init(&file_result, data, threadData);
1607 ✗ rust_result_writeParameterData(&file_result, data, threadData);
1608 ✗ rust_result_emit(&file_result, data, threadData);
1609 ✗ rust_result_free(&file_result, data, threadData);
1610 ✗ data->simulationInfo->outputFormat = outputFormat;
1611
1612 /* exit as we do not want to compute any of the following variables; the
1613 * throw unwinds into the solver, which keeps its message to OMC_LOG_NLS */
1614 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Success: Initial guess has been written to disk (path = %s). The program will terminate now.", path);
1615 ✗ throwStreamPrint(threadData, "Initial guess for non-linear system %d written to disk.", nls_idx);
1616 }
1617 else
1618 {
1619 /* not the NLS index that had been specified, continue with the systems prior to the user specified one */
1620 ✗ return;
1621 }
1622 }
1623
1624 ✗ static void B_check_stop_at_system(threadData_t *threadData, NONLINEAR_SYSTEM_DATA *nlsData)
1625 {
1626 ✗ int nls_idx = atoi(omc_flagValue[FLAG_STOP_AT_SYSTEM]);
1627 ✗ if (nls_idx == nlsData->equationIndex)
1628 {
1629 ✗ throwStreamPrint(threadData, "Success: Finished solving specified NLS system with index %d. The program will terminate now.\n", nls_idx);
1630 }
1631 ✗ }
1632
1633 /**
1634 * @brief Solve non-linear system with KINSol
1635 *
1636 * @param data Runtime data struct.
1637 * @param threadData Thread data for error handling.
1638 * @param nlsData Pointer to non-linear system data.
1639 * @return NLS_SOLVER_STATUS Return NLS_SOLVED on success and NLS_FAILED otherwise.
1640 */
1641 ✗ NLS_SOLVER_STATUS B_nlsKinsolSolve(DATA* data, threadData_t* threadData, NONLINEAR_SYSTEM_DATA* nlsData) {
1642
1643 ✗ B_NLS_KINSOL_DATA *kinsolData = (B_NLS_KINSOL_DATA *)nlsData->solverData;
1644 ✗ int eqSystemNumber = nlsData->equationIndex;
1645 ✗ int indexes[2] = {1, eqSystemNumber};
1646
1647 int flag;
1648 long nFEval;
1649 modelica_boolean success = FALSE;
1650 modelica_boolean retry = TRUE;
1651 NLS_SOLVER_STATUS solver_status;
1652 ✗ double *xStart = NV_DATA_S(kinsolData->initialGuess);
1653 double fNormValue;
1654
1655 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_NLS_V, omc_dummyFileInfo, 1, indexes,
1656 "Start solving Non-Linear System %d (size %d) at time %g with Kinsol Solver",
1657 ✗ eqSystemNumber, (int) nlsData->size, data->localData[0]->timeValue);
1658
1659 /* Solve nonlinear system with KINSol() */
1660 ✗ kinsolData->retries = 0;
1661 do {
1662 // FIXME: This entire interface needs a complete redesign from the ground up.
1663 // With the new scaling logic, everything becomes tangled and the control flow increasingly unclear.
1664
1665 ✗ kinsolData->useScaling = FALSE;
1666
1667 // set x := unscaled x (solution of prev / initial guess)
1668 ✗ B_nlsKinsolResetInitialUnscaled(data, kinsolData, nlsData, B_INITIAL_EXTRAPOLATION);
1669
1670 // create new x scaling based on x and x_nominal
1671 ✗ B_nlsKinsolXScaling(data, kinsolData, nlsData, B_SCALING_NOMINALSTART);
1672
1673 // create new f scaling based on Jacobian in memory (from prev solve) or if none compute new scaling
1674 ✗ B_nlsKinsolFScaling(data, kinsolData, nlsData, B_SCALING_JACOBIAN);
1675
1676 /* Set maximum step size */
1677 ✗ B_nlsKinsolSetMaxNewtonStep(kinsolData, kinsolData->maxstepfactor);
1678
1679 /* Dump configuration */
1680 ✗ B_nlsKinsolConfigPrint(kinsolData, nlsData);
1681
1682 ✗ kinsolData->useScaling = TRUE;
1683 ✗ nlsKinsolInplaceScaleX(kinsolData, kinsolData->initialGuess);
1684 ✗ nlsKinsolInplaceScaleJac(nlsData, kinsolData, kinsolData->J);
1685
1686 /* TODO: This should be another flag, e.g. LOG_NLS_JAC_UPDATE and not OMC_LOG_NLS_DERIVATIVE_TEST
1687 only in some cases this derivative test makes sense, since the scaled Jacobian is outdated frequently!
1688 in many cases, we use an outdated jacobian here, such that errors explode and it detects wrong Jacobian mismatches
1689 that are due to the dense Jacobian evaluated at the new point x_new.
1690
1691 if (omc_useStream[OMC_LOG_NLS_DERIVATIVE_TEST])
1692 {
1693 nlsKinsolDenseDerivativeTest(data, nlsData, kinsolData, kinsolData->J, KINSOL_B_ENTRY_POINT);
1694 }
1695 */
1696
1697 ✗ if (omc_useStream[OMC_LOG_NLS_JAC_SUMS])
1698 {
1699 ✗ nlsJacobianRowColSums(data, nlsData, kinsolData->J, KINSOL_B_ENTRY_POINT /* called at entry point */, kinsolData->useScaling /* scaled */);
1700 }
1701
1702 ✗ if (omc_useStream[OMC_LOG_NLS_SVD] || omc_useStream[OMC_LOG_NLS_SVD_V])
1703 {
1704 ✗ svd_compute(data, nlsData, SM_DATA_S(kinsolData->J), kinsolData->useScaling, KINSOL_B_ENTRY_POINT /* called at entry point */);
1705 }
1706
1707 ✗ if (omc_flag[FLAG_SAVE_INITIAL_GUESS_SYSTEM]) {
1708 ✗ B_save_initial_guess_system(data, threadData, nlsData);
1709 }
1710
1711 ✗ flag = KINSol(
1712 kinsolData->kinsolMemory, /* KINSol memory block */
1713 kinsolData->initialGuess, /* initial guess on input; solution vector */
1714 kinsolData->kinsolStrategy, /* global strategy choice */
1715 kinsolData->ONES_xScale, /* (1, ..., 1)^T */
1716 kinsolData->ONES_fScale); /* (1, ..., 1)^T */
1717
1718 ✗ if (flag < 0 && kinsolData->attemptRetry) {
1719 ✗ warningStreamPrint(OMC_LOG_NLS, 0, "KINSol finished with errorCode %d.", flag);
1720 } else {
1721 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSol finished with errorCode %d.", flag);
1722 }
1723 /* Try to handle recoverable errors */
1724 ✗ retry = flag < 0 && kinsolData->attemptRetry && nlsKinsolErrorHandler(flag, data, nlsData, kinsolData);
1725
1726 /* solution found */
1727 ✗ if ((flag == KIN_SUCCESS) || (flag == KIN_INITIAL_GUESS_OK) ||
1728 (flag == KIN_STEP_LT_STPTOL)) {
1729 success = TRUE;
1730 }
1731 ✗ kinsolData->retries++;
1732
1733 /* write statistics */
1734 ✗ KINGetNumNonlinSolvIters(kinsolData->kinsolMemory, &nFEval);
1735 ✗ nlsData->numberOfIterations += nFEval;
1736 ✗ nlsData->numberOfFEval = kinsolData->countResCalls;
1737
1738 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "Next try? success = %d, retry = %d, retries = %d = %s\n",
1739 success, retry, kinsolData->retries,
1740 ✗ !success && !retry && kinsolData->retries < B_RETRY_MAX ? "true" : "false");
1741 ✗ } while (!success && retry && kinsolData->retries < B_RETRY_MAX);
1742
1743 /* Check solution status */
1744 ✗ if (success && kinsolData->resetTol) {
1745 ✗ kinsolData->solved = NLS_SOLVED_LESS_ACCURACY;
1746 ✗ } else if (success) {
1747 ✗ kinsolData->solved = NLS_SOLVED;
1748 } else {
1749 ✗ kinsolData->solved = NLS_FAILED;
1750 }
1751
1752 /* Reset solver tolerance */
1753 ✗ if (kinsolData->resetTol) {
1754 ✗ KINSetFuncNormTol(kinsolData->kinsolMemory, kinsolData->fnormtol);
1755 ✗ KINSetScaledStepTol(kinsolData->kinsolMemory, kinsolData->scsteptol);
1756 ✗ kinsolData->resetTol = FALSE;
1757 }
1758
1759 ✗ if (success) {
1760 ✗ if (kinsolData->useScaling) {
1761 ✗ nlsKinsolInplaceUnscaleX(kinsolData, kinsolData->initialGuess);
1762
1763 /* repeated solve; we must unscale the previous Jacobian, since we will compute the new fScalings
1764 in the next step from the Jacobian that is already in memory */
1765 ✗ nlsKinsolInplaceUnscaleJac(nlsData, kinsolData, kinsolData->J);
1766 ✗ kinsolData->useScaling = FALSE;
1767 }
1768 ✗ memcpy(nlsData->nlsx, xStart, nlsData->size * (sizeof(double)));
1769 }
1770
1771 ✗ messageClose(OMC_LOG_NLS_V);
1772
1773 ✗ if (omc_flag[FLAG_STOP_AT_SYSTEM]) {
1774 ✗ B_check_stop_at_system(threadData, nlsData);
1775 }
1776
1777 ✗ return kinsolData->solved;
1778 }
1779
1780 #else /* WITH_SUNDIALS */
1781
1782 void* B_nlsKinsolAllocate(int size, void* userData, int attemptRetry, modelica_boolean isPatternAvailable) {
1783
1784 throwStreamPrint(NULL, "No sundials/kinsol support activated.");
1785 return 0;
1786 }
1787
1788 int B_nlsKinsolFree(void* kinsolData) {
1789
1790 throwStreamPrint(NULL, "No sundials/kinsol support activated.");
1791 return 0;
1792 }
1793
1794 int B_nlsKinsolSolve(void *data, threadData_t *threadData, void* nlsData) {
1795
1796 throwStreamPrint(threadData, "No sundials/kinsol support activated.");
1797 return 0;
1798 }
1799
1800 #endif /* WITH_SUNDIALS */
1801