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 / 595
Functions: 0.0% 0 / 0 / 20
Branches: 0.0% 0 / 0 / 309
OMCompiler/SimulationRuntime/c/simulation/solver/kinsolSolver.c
Line Branch Exec Source
1 /*
2 * This file belongs to the OpenModelica Run-Time System
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings
5 * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights
6 * reserved.
7 *
8 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
9 * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY
10 * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
11 * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL
12 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
13 *
14 * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License
15 * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs:
16 * http://www.openmodelica.org or https://github.com/OpenModelica/ or
17 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU
18 * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW
19 * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY
23 * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF
24 * OSMC-PL.
25 *
26 */
27
28 /*! \file kinsolSolver.c
29 */
30
31 #include "kinsolSolver.h"
32
33 #include "nonlinearSystem.h"
34 #include "omc_config.h"
35 #include "omc_math.h"
36 #include "../options.h"
37 #include "../simulation_info_json.h"
38 #include "../jacobian_util.h"
39 #include "sundials_util.h"
40 #include "util/omc_error.h"
41
42 #ifdef WITH_SUNDIALS
43
44 #include "events.h"
45 #include "model_help.h"
46 #include "openmodelica.h"
47 #include "openmodelica_func.h"
48 #include "util/read_matlab4.h"
49 #include "util/varinfo.h"
50
51 #include <math.h>
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include <string.h>
55
56 /* Function prototypes */
57 static int nlsKinsolResiduals(N_Vector x, N_Vector f, void* userData);
58 static int nlsSparseJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
59 void* userData, N_Vector tmp1, N_Vector tmp2);
60 int nlsSparseSymJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
61 void* userData, N_Vector tmp1, N_Vector tmp2);
62 static int nlsDenseJac(long int N, N_Vector vecX, N_Vector vecFX,
63 SUNMatrix Jac, NLS_USERDATA *kinsolUserData,
64 N_Vector tmp1, N_Vector tmp2);
65 static void nlsKinsolJacSumSparse(SUNMatrix A);
66 static void nlsKinsolJacSumDense(SUNMatrix A);
67
68 /**
69 * @brief Set KINSOL configuration.
70 *
71 * @param kinsolData Kinsol data with configuration settings.
72 */
73 ✗ static void nlsKinsolConfigSetup(NLS_KINSOL_DATA *kinsolData) {
74 /* Variables */
75 int flag;
76
77 /* configuration */
78 ✗ flag = KINSetFuncNormTol(kinsolData->kinsolMemory,
79 kinsolData->fnormtol); /* Set function-norm stopping tolerance */
80 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetFuncNormTol");
81 ✗ kinsolData->resetTol = FALSE;
82
83 ✗ flag = KINSetScaledStepTol(kinsolData->kinsolMemory,
84 kinsolData->scsteptol); /* Set scaled-step stopping tolerance */
85 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetScaledStepTol");
86
87 ✗ flag = KINSetNumMaxIters(kinsolData->kinsolMemory,
88 ✗ 100 * kinsolData->size); /* Set max. number of nonlinear iterations */
89 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNumMaxIters");
90
91 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH; /* Newton with globalization strategy to solve nonlinear systems */
92
93 ✗ flag = KINSetNoInitSetup(kinsolData->kinsolMemory, SUNFALSE); /* TODO: This is the default value. Is there a point in calling this function? */
94 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNoInitSetup");
95
96 ✗ kinsolData->retries = 0;
97 ✗ kinsolData->countResCalls = 0;
98 ✗ }
99
100 /**
101 * @brief Initialize KINSOL data.
102 *
103 * Allocate memory for KINSOL data and Jacobian.
104 *
105 * @param kinsolData KINSOL data.
106 */
107 ✗ void initKinsolMemory(NLS_KINSOL_DATA *kinsolData) {
108 int flag;
109 ✗ int size = kinsolData->size;
110 ✗ NONLINEAR_SYSTEM_DATA *nlsData = kinsolData->userData->nlsData;
111 ✗ SPARSE_PATTERN* sparsePattern = nlsData->sparsePattern;
112
113 /* Free KINSOL memory block */
114 ✗ if (kinsolData->kinsolMemory != NULL || kinsolData->J != NULL || kinsolData->scaledJ != NULL) {
115 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
116 "KINSOL: Already allocated kinsol memory. Loosing memory!");
117 }
118
119 /* Create KINSOL memory block. The SUNDIALS context was created by
120 * nlsKinsolAllocate, which has to happen before any SUNDIALS object. */
121 ✗ kinsolData->kinsolMemory = KINCreate(kinsolData->sunctx);
122 ✗ if (kinsolData->kinsolMemory == NULL) {
123 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
124 "KINSOL: In function KINCreate: An error occurred.");
125 }
126
127 ✗ flag = KINSetUserData(kinsolData->kinsolMemory, (void*)kinsolData->userData);
128 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetUserData");
129
130 /* Initialize KINSOL object */
131 ✗ flag = KINInit(kinsolData->kinsolMemory, nlsKinsolResiduals,
132 kinsolData->initialGuess);
133 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINInit");
134
135 /* Create matrix object */
136 ✗ if (kinsolData->linearSolverMethod == NLS_LS_DEFAULT ||
137 kinsolData->linearSolverMethod == NLS_LS_LAPACK) {
138 ✗ kinsolData->J = SUNDenseMatrix(size, size, kinsolData->sunctx);
139 ✗ } else if (kinsolData->linearSolverMethod == NLS_LS_KLU) {
140 ✗ if (!sparsePattern) {
141 ✗ kinsolData->nnz = size*size;
142 } else {
143 ✗ kinsolData->nnz = sparsePattern->nnz;
144 }
145 ✗ kinsolData->J = SUNSparseMatrix(size, size, kinsolData->nnz, SUN_CSC_MAT, kinsolData->sunctx);
146 ✗ kinsolData->scaledJ = SUNSparseMatrix(size, size, kinsolData->nnz, SUN_CSC_MAT, kinsolData->sunctx);
147 }
148
149 /* Create linear solver object */
150 ✗ if (kinsolData->linearSolverMethod == NLS_LS_DEFAULT ||
151 kinsolData->linearSolverMethod == NLS_LS_TOTALPIVOT) {
152 ✗ kinsolData->linSol = SUNLinSol_Dense(kinsolData->y, kinsolData->J, kinsolData->sunctx);
153 ✗ if (kinsolData->linSol == NULL) {
154 ✗ throwStreamPrint(NULL, "KINSOL: In function SUNLinSol_Dense: Input incompatible.");
155 }
156 ✗ } else if (kinsolData->linearSolverMethod == NLS_LS_LAPACK) {
157 ✗ kinsolData->linSol = SUNLinSol_LapackDense(kinsolData->y, kinsolData->J, kinsolData->sunctx);
158 ✗ if (kinsolData->linSol == NULL) {
159 ✗ throwStreamPrint(NULL, "KINSOL: In function SUNLinSol_LapackDense: Input incompatible.");
160 }
161 ✗ } else if (kinsolData->linearSolverMethod == NLS_LS_KLU) {
162 ✗ kinsolData->linSol = SUNLinSol_KLU(kinsolData->y, kinsolData->J, kinsolData->sunctx);
163 ✗ if (kinsolData->linSol == NULL) {
164 ✗ throwStreamPrint(NULL, "KINSOL: In function SUNLinSol_KLU: Input incompatible.");
165 }
166 } else {
167 ✗ throwStreamPrint(NULL, "KINSOL: Unknown linear solver method.");
168 }
169 /* Log used solver */
170 ✗ infoStreamPrint(OMC_LOG_NLS, 0, "KINSOL: Using linear solver method %s", NLS_LS_METHOD_NAME[kinsolData->linearSolverMethod]);
171
172 /* Set linear solver */
173 ✗ flag = KINSetLinearSolver(kinsolData->kinsolMemory, kinsolData->linSol,
174 kinsolData->J);
175 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KINLS_FLAG, "KINSetLinearSolver");
176
177 /* Set Jacobian for non-linear solver */
178 ✗ if (kinsolData->linearSolverMethod == NLS_LS_KLU) {
179 ✗ if (nlsData->analyticalJacobianColumn != NULL && sparsePattern != NULL) {
180 ✗ flag = KINSetJacFn(kinsolData->kinsolMemory, nlsSparseSymJac); /* Use symbolic Jacobian with sparsity pattern*/
181 ✗ } else if (sparsePattern != NULL) {
182 ✗ flag = KINSetJacFn(kinsolData->kinsolMemory, nlsSparseJac); /* Use numeric Jacobian with sparsity pattern */
183 } else {
184 ✗ throwStreamPrint(NULL, "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.");
185 }
186 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KINLS_FLAG, "KINSetJacFn");
187 }
188
189 /* Configuration */
190 ✗ nlsKinsolConfigSetup(kinsolData);
191 ✗ }
192
193 /**
194 * @brief Allocate memory for kinsol solver data and initialize KINSOL solver.
195 *
196 * @param size Size of non-linear problem.
197 * @param userData Pointer to set NLS user data.
198 * @param attemptRetry True if KINSOL should retry with different settings after solution failed.
199 * @param isPatternAvailable True if sparsity pattern of Jacobian is available. Allocate work vectors for KLU in that case.
200 * @return NLS_KINSOL_DATA* Pointer to allocated KINSOL data.
201 */
202 ✗ NLS_KINSOL_DATA* nlsKinsolAllocate(int size, NLS_USERDATA* userData, modelica_boolean attemptRetry, modelica_boolean isPatternAvailable) {
203 /* Allocate system data */
204 ✗ NLS_KINSOL_DATA *kinsolData = (NLS_KINSOL_DATA *)calloc(1, sizeof(NLS_KINSOL_DATA));
205
206 ✗ kinsolData->size = size;
207 ✗ kinsolData->linearSolverMethod = userData->nlsData->nlsLinearSolver;
208 ✗ kinsolData->solved = NLS_FAILED;
209 ✗ kinsolData->userData = userData;
210
211 ✗ if (SUNContext_Create(SUN_COMM_NULL, &kinsolData->sunctx) != SUN_SUCCESS) {
212 ✗ throwStreamPrint(NULL, "KINSOL: In function SUNContext_Create: An error occurred.");
213 }
214 ✗ sundialsSilenceLogger(kinsolData->sunctx);
215
216 /* Set error handler */
217 ✗ if (SUNContext_PushErrHandler(kinsolData->sunctx, kinsolErrorHandlerFunction, kinsolData) != SUN_SUCCESS) {
218 ✗ throwStreamPrint(NULL, "KINSOL: In function SUNContext_PushErrHandler: An error occurred.");
219 }
220
221 ✗ kinsolData->fnormtol = newtonFTol; /* function tolerance */
222 ✗ kinsolData->scsteptol = newtonXTol; /* step tolerance */
223
224 ✗ kinsolData->maxstepfactor = maxStepFactor; /* step tolerance */
225 ✗ kinsolData->nominalJac = 0; /* calculate for scaling the scaled matrix */
226 ✗ kinsolData->attemptRetry = attemptRetry;
227
228 ✗ kinsolData->initialGuess = N_VNew_Serial(size, kinsolData->sunctx);
229 ✗ kinsolData->xScale = N_VNew_Serial(size, kinsolData->sunctx);
230 ✗ kinsolData->fScale = N_VNew_Serial(size, kinsolData->sunctx);
231 ✗ kinsolData->constraints = N_VNew_Serial(size, kinsolData->sunctx);
232 ✗ kinsolData->fRes = N_VNew_Serial(size, kinsolData->sunctx);
233 ✗ kinsolData->fTmp = N_VNew_Serial(size, kinsolData->sunctx);
234
235 ✗ kinsolData->y = N_VNew_Serial(size, kinsolData->sunctx);
236 ✗ kinsolData->J = NULL;
237
238 /* tmp1, tmp2 only needed for numeric Jacobian */
239 ✗ if (userData->nlsData->analyticalJacobianColumn != NULL &&
240 ✗ isPatternAvailable &&
241 ✗ kinsolData->linearSolverMethod == NLS_LS_KLU)
242 {
243 ✗ kinsolData->tmp1 = NULL;
244 ✗ kinsolData->tmp2 = NULL;
245 } else {
246 ✗ kinsolData->tmp1 = N_VNew_Serial(size, kinsolData->sunctx);
247 ✗ kinsolData->tmp2 = N_VNew_Serial(size, kinsolData->sunctx);
248 }
249 /* Scaled Jacobian is allocated with J */
250 ✗ kinsolData->scaledJ = NULL;
251
252 ✗ kinsolData->kinsolMemory = NULL;
253
254 ✗ initKinsolMemory(kinsolData);
255
256 ✗ return kinsolData;
257 }
258
259 /**
260 * @brief Deallocates memory for KINSOL solver.
261 *
262 * Free memory that was allocated with `nlsKinsolAllocate`.
263 *
264 * @param kinsolData Pointer to KINSOL data.
265 */
266 ✗ void nlsKinsolFree(NLS_KINSOL_DATA* kinsolData) {
267 ✗ KINFree((void *)&kinsolData->kinsolMemory);
268
269 ✗ N_VDestroy_Serial(kinsolData->initialGuess);
270 ✗ N_VDestroy_Serial(kinsolData->xScale);
271 ✗ N_VDestroy_Serial(kinsolData->fScale);
272 ✗ N_VDestroy_Serial(kinsolData->constraints);
273 ✗ N_VDestroy_Serial(kinsolData->fRes);
274 ✗ N_VDestroy_Serial(kinsolData->fTmp);
275
276 /* Free linear solver data */
277 ✗ SUNLinSolFree(kinsolData->linSol);
278 ✗ SUNMatDestroy(kinsolData->J);
279 ✗ SUNMatDestroy(kinsolData->scaledJ);
280 ✗ N_VDestroy_Serial(kinsolData->y);
281 ✗ if (kinsolData->tmp1 != NULL) {
282 ✗ N_VDestroy_Serial(kinsolData->tmp1);
283 ✗ N_VDestroy_Serial(kinsolData->tmp2);
284 }
285
286 /* The context has to outlive every SUNDIALS object created with it */
287 ✗ SUNContext_Free(&kinsolData->sunctx);
288
289 ✗ freeNlsUserData(kinsolData->userData);
290 ✗ free(kinsolData);
291
292 ✗ return;
293 }
294
295 /**
296 * @brief Residual function for non-linear problem.
297 *
298 * @param x The current value of the variable vector.
299 * @param f Output vector.
300 * @param userData Pointer to Kinsol user data.
301 * @return int Return 0 on success, return 1 on recoverable error.
302 */
303 ✗ static int nlsKinsolResiduals(N_Vector x, N_Vector f, void* userData) {
304 ✗ double *xdata = NV_DATA_S(x);
305 ✗ double *fdata = NV_DATA_S(f);
306
307 NLS_USERDATA* kinsolUserData = (NLS_USERDATA*)userData;
308 ✗ DATA* data = kinsolUserData->data;
309 ✗ threadData_t* threadData = kinsolUserData->threadData;
310 ✗ NONLINEAR_SYSTEM_DATA* nlsData = kinsolUserData->nlsData;
311 ✗ NLS_KINSOL_DATA* kinsolData = (NLS_KINSOL_DATA*)nlsData->solverData;
312 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=kinsolUserData->solverData};
313 ✗ int iflag = 1 /* recoverable error */;
314
315 /* Update statistics */
316 ✗ kinsolData->countResCalls++;
317
318 #ifndef OMC_EMCC
319 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
320 #endif
321
322 /* call residual function */
323 ✗ nlsData->residualFunc(&resUserData, xdata, fdata, (const int *)&iflag);
324 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { iflag = 0 /* success */; }
325
326 #ifndef OMC_EMCC
327 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
328 #endif
329
330 ✗ return iflag;
331 }
332
333 /**
334 * @brief Calculate dense Jacobian matrix.
335 *
336 * @param N Size of vecX and vecFX.
337 * @param vecX Vector x.
338 * @param vecFX Residual vector f(x).
339 * @param Jac Dense Jacobian matrix J(x).
340 * @param kinsolUserData Pointer to Kinsol user data.
341 * @param tmp1 Unused, only to match interface of KINLsJacFn
342 * @param tmp2 Unused, only to match interface of KINLsJacFn
343 * @return int Return 0 on success, -1 on failure.
344 */
345 ✗ static int nlsDenseJac(long int N,
346 N_Vector vecX,
347 N_Vector vecFX,
348 SUNMatrix Jac,
349 NLS_USERDATA *kinsolUserData,
350 N_Vector tmp1,
351 N_Vector tmp2) {
352 DATA *data = kinsolUserData->data;
353 threadData_t *threadData = kinsolUserData->threadData;
354 ✗ NONLINEAR_SYSTEM_DATA *nlsData = kinsolUserData->nlsData;
355 ✗ NLS_KINSOL_DATA *kinsolData = (NLS_KINSOL_DATA *)nlsData->solverData;
356
357 ✗ if (SUNMatGetID(Jac) != SUNMATRIX_DENSE) {
358 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
359 "KINSOL: nlsDenseJac illegal input Jac. Matrix is not dense!");
360 ✗ return -1;
361 }
362
363 /* prepare variables */
364 ✗ double *x = N_VGetArrayPointer(vecX);
365 ✗ double *fx = N_VGetArrayPointer(vecFX);
366 ✗ double *xScaling = NV_DATA_S(kinsolData->xScale);
367 ✗ double *fRes = NV_DATA_S(kinsolData->fRes);
368 double xsave, xscale, sign;
369 double delta_hh;
370 const double delta_h = sqrt(DBL_EPSILON * 2e1);
371
372 long int i, j;
373
374 /* performance measurement */
375 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
376
377 /* Use forward difference quotient to approximate Jacobian */
378 ✗ for (i = 0; i < N; i++) {
379 ✗ xsave = x[i];
380 ✗ delta_hh = delta_h * (fabs(xsave) + 1.0);
381 ✗ if ((xsave + delta_hh >= nlsData->max[i])) {
382 ✗ delta_hh *= -1.0;
383 }
384 ✗ x[i] += delta_hh;
385
386 /* Evaluate Jacobian function */
387 ✗ nlsKinsolResiduals(vecX, kinsolData->fRes, kinsolUserData);
388
389 /* Calculate scaled difference quotient */
390 ✗ delta_hh = 1.0 / delta_hh;
391
392 ✗ for (j = 0; j < N; j++) {
393 ✗ if (kinsolData->nominalJac) {
394 ✗ SM_ELEMENT_D(Jac, j, i) = (fRes[j] - fx[j]) * delta_hh / xScaling[i];
395 } else {
396 ✗ SM_ELEMENT_D(Jac, j, i) =
397 ✗ (fRes[j] - fx[j]) * delta_hh; /* TODO: Or now Jac(i,j) ??? */
398 }
399 }
400 ✗ x[i] = xsave;
401 }
402
403 /* debug */
404 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC)) {
405 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "KINSOL: Dense matrix.");
406 ✗ SUNDenseMatrix_Print(Jac, stdout); /* TODO: Print in OMC_LOG_NLS_JAC */
407 ✗ nlsKinsolJacSumDense(Jac);
408 ✗ messageClose(OMC_LOG_NLS_JAC);
409 }
410
411 /* performance measurement and statistics */
412 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
413 ✗ nlsData->numberOfJEval++;
414
415 ✗ return 0;
416 }
417
418 /**
419 * @brief Finish sparse matrix by fixing colprts.
420 *
421 * Last value of indexptrs should always be nnz.
422 * Search for empty rows which would mean the matrix is singular.
423 *
424 * @param A CSC matrix
425 */
426 ✗ static void finishSparseColPtr(SUNMatrix A, int nnz) {
427 int i;
428
429 /* TODO: Remove this check for performance reasons? */
430 ✗ if (SM_SPARSETYPE_S(A) != SUN_CSC_MAT) {
431 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
432 "KINSOL: In function finishSparseColPtr: Wrong sparse format of SUNMatrix A.");
433 }
434
435 /* Set last value of indexptrs to nnz */
436 ✗ SM_INDEXPTRS_S(A)[SM_COLUMNS_S(A)] = nnz;
437
438 /* Check for empty rows */
439 ✗ for (i = 1; i < SM_COLUMNS_S(A) + 1; ++i) {
440 ✗ if (SM_INDEXPTRS_S(A)[i] == SM_INDEXPTRS_S(A)[i - 1]) {
441 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0,
442 "KINSOL: Jacobian row %d singular. See OMC_LOG_NLS for "
443 "more information.",
444 i);
445 ✗ SM_INDEXPTRS_S(A)[i] = SM_INDEXPTRS_S(A)[i - 1];
446 }
447 }
448 ✗ }
449
450 // TODO: unify this up to a generic level, such that we can use this from pretty much every solver and
451 // it does not take kinsolData and the matrix only as flat buffer + SPARSE_PATTERN
452
453 /**
454 * @brief Perform derivative test comparing symbolic and numerical Jacobians for KINSOL
455 *
456 * Compares the symbolic Jacobian (sparse CSC format) with a numerically approximated
457 * dense Jacobian, checking for numerical and structural anomalies. The numerical
458 * Jacobian is computed using finite differences via nlsDenseJac.
459 *
460 * @param data Runtime data structure
461 * @param nlsData Nonlinear system data
462 * @param kinsolData KINSOL solver data structure
463 * @param Jsym Symbolic Jacobian in sparse CSC format
464 * @param tol Tolerance, all relative errors above tol are considered anomalies
465 * @param newJac TRUE if called from jacobian evaluation, FALSE if called from solver entry point
466 *
467 * @return int 1 derivative test failed and no error
468 * 0 derivative test successful and no error
469 * -1 internal error
470 */
471 ✗ static int nlsKinsolDenseDerivativeTest(DATA *data, NONLINEAR_SYSTEM_DATA *nlsData,
472 NLS_KINSOL_DATA *kinsolData, SUNMatrix Jsym, SolverCaller caller)
473 {
474 int row, col, nz, numericalErrorCount, structuralErrorCount;
475 ✗ const int size = nlsData->size;
476 int ret = 0;
477
478 modelica_real symValue, numValue, absError, relError;
479 modelica_real maxError = 0.0;
480
481 modelica_boolean errorFound;
482
483 ✗ sunindextype nnz = SUNSparseMatrix_NNZ(Jsym);
484 ✗ sunindextype columns = SUNSparseMatrix_Columns(Jsym);
485 ✗ sunindextype rows = SUNSparseMatrix_Rows(Jsym);
486
487 ✗ sunindextype *colPointers = SM_INDEXPTRS_S(Jsym);
488 ✗ sunindextype *rowIndices = SM_INDEXVALS_S(Jsym);
489 ✗ sunrealtype *symValues = SM_DATA_S(Jsym);
490
491 // allocate temporary memory for dense finite-diff matrix
492 ✗ N_Vector vecX = N_VNew_Serial(size, kinsolData->sunctx);
493 ✗ N_Vector vecFX = N_VNew_Serial(size, kinsolData->sunctx);
494 ✗ N_Vector tmp1 = N_VNew_Serial(size, kinsolData->sunctx);
495 ✗ N_Vector tmp2 = N_VNew_Serial(size, kinsolData->sunctx);
496 ✗ SUNMatrix Jnum = SUNDenseMatrix(size, size, kinsolData->sunctx);
497
498 // set tolerances
499 ✗ modelica_real Atol = omc_flag[FLAG_NLS_JAC_TEST_ATOL] ? atof(omc_flagValue[FLAG_NLS_JAC_TEST_ATOL]) : 100 * DBL_EPSILON;
500 ✗ modelica_real Rtol = omc_flag[FLAG_NLS_JAC_TEST_RTOL] ? atof(omc_flagValue[FLAG_NLS_JAC_TEST_RTOL]) : 1e-4;
501
502 // copy current x into new vector, compute f(x) and corresponding dense finite-diff Jacobian
503 ✗ SUNMatZero(Jnum);
504 ✗ N_VScale(1.0, kinsolData->initialGuess, vecX);
505 ✗ nlsKinsolResiduals(vecX, vecFX, kinsolData->userData);
506 ✗ if (nlsDenseJac(size, vecX, vecFX, Jnum, kinsolData->userData, tmp1, tmp2) != 0)
507 {
508 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Numerical Jacobian computation failed in nlsKinsolDenseDerivativeTest");
509 ret = -1;
510 ✗ SUNMatDestroy(Jnum);
511 ✗ N_VDestroy_Serial(vecX);
512 ✗ N_VDestroy_Serial(vecFX);
513 ✗ N_VDestroy_Serial(tmp1);
514 ✗ N_VDestroy_Serial(tmp2);
515 ✗ return ret;
516 }
517
518 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "%s: Derivative test (atol=%.5e, rtol=%.5e, scaled = %s, Caller: %s):",
519 ✗ SolverCaller_callerString(caller), Atol, Rtol, kinsolData->nominalJac ? "true" : "false", SolverCaller_toString(caller));
520 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Matrix Info");
521 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "NLS index = " OMC_INT_FORMAT, nlsData->equationIndex);
522 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Columns = " OMC_INT_FORMAT, columns);
523 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Rows = " OMC_INT_FORMAT, rows);
524 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "NNZ = " OMC_INT_FORMAT, nnz);
525 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Curr Time = %-11.5e", data->localData[0]->timeValue);
526
527 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
528
529 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Anomalies");
530
531 nz = 0;
532 numericalErrorCount = 0;
533 structuralErrorCount = 0;
534
535 ✗ for (col = 0; col < size; col++)
536 {
537 errorFound = FALSE;
538
539 ✗ for (row = 0; row < size; row++)
540 {
541 ✗ numValue = SM_ELEMENT_D(Jnum, row, col);
542
543 ✗ if (colPointers[col] <= nz && nz < colPointers[col+1] && rowIndices[nz] == row)
544 {
545 // structural non-zero -> compare values
546 ✗ symValue = symValues[nz++];
547 ✗ absError = fabs(symValue - numValue);
548 ✗ relError = (absError < Atol) ? 0.0 : absError / fmax(fabs(numValue), fabs(symValue));
549
550 ✗ if (relError > maxError)
551 {
552 maxError = relError;
553 }
554
555 ✗ if (relError > Rtol)
556 {
557 // tolerance exceeded -> numerical error
558 ✗ if (!errorFound)
559 {
560 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Column / Variable: %i, Name: %s",
561 ✗ col + 1, modelInfoGetEquation(&data->modelData->modelDataXml, nlsData->equationIndex).vars[col]);
562 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6s %-6s %-15s %-15s %-8s",
563 "Type", "Col", "Row", "Symbolic", "Numerical", "RelError");
564 errorFound = TRUE;
565 }
566 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6d %-6d %+15.8e %+15.8e %+13.8e",
567 "Numerical", col + 1, row + 1, symValue, numValue, relError);
568 ✗ numericalErrorCount++;
569 }
570 }
571 ✗ else if (fabs(numValue) > Atol)
572 {
573 // structural error with tolerance exceeded -> non-zero in numerical Jacobian but zero in symbolic
574 ✗ if (!errorFound)
575 {
576 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Column / Variable: %i, Name: %s",
577 ✗ col + 1, modelInfoGetEquation(&data->modelData->modelDataXml, nlsData->equationIndex).vars[col]);
578 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6s %-6s %-15s %-15s %-8s",
579 "Type", "Col", "Row", "Symbolic", "Numerical", "RelError");
580 errorFound = TRUE;
581 }
582 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "%-12s %-6d %-6d %+15.8e %+15.8e %+13.8e",
583 "Structural", col + 1, row + 1, 0.0, numValue, 1.0);
584 ✗ structuralErrorCount++;
585 }
586 }
587
588 ✗ if (errorFound)
589 {
590 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
591 }
592 }
593 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
594
595 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 1, "Summary");
596 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Numerical errors: %d (value mismatch w.r.t. reference)", numericalErrorCount);
597 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Structural errors: %d (non-zero not in sparsity pattern)", structuralErrorCount);
598 ✗ infoStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Max relative error: %.3e", maxError);
599
600 ✗ if (numericalErrorCount + structuralErrorCount > 0)
601 {
602 ✗ warningStreamPrint(OMC_LOG_NLS_DERIVATIVE_TEST, 0, "Derivative test failed (%d numerical, %d structural errors)",
603 numericalErrorCount, structuralErrorCount);
604 ret = 1;
605 }
606 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
607
608 ✗ SUNMatDestroy(Jnum);
609 ✗ N_VDestroy_Serial(vecX);
610 ✗ N_VDestroy_Serial(vecFX);
611 ✗ N_VDestroy_Serial(tmp1);
612 ✗ N_VDestroy_Serial(tmp2);
613
614 ✗ messageClose(OMC_LOG_NLS_DERIVATIVE_TEST);
615
616 ✗ return ret;
617 }
618
619 /**
620 * @brief Computes symbolic Jacobian matrix Jac(vecX)
621 *
622 * @param vecX
623 * @param vecFX just for interface compatibility, will not be used here
624 * @param Jac Allocated Jacobian, contains symbolic Jacobian on exit
625 * @param userData Void pointer to user data of type NLS_USERDATA*.
626 * @param tmp1 Unused, only to match interface of KINLsJacFn
627 * @param tmp2 Unused, only to match interface of KINLsJacFn
628 * @return int
629 */
630 ✗ int nlsSparseSymJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
631 void *userData, N_Vector tmp1, N_Vector tmp2) {
632 /* Variables */
633 NLS_USERDATA* kinsolUserData = (NLS_USERDATA *)userData;;
634 ✗ DATA* data = kinsolUserData->data;
635 ✗ threadData_t* threadData = kinsolUserData->threadData;
636 ✗ NONLINEAR_SYSTEM_DATA* nlsData = kinsolUserData->nlsData;
637 ✗ NLS_KINSOL_DATA* kinsolData = (NLS_KINSOL_DATA *)nlsData->solverData;
638 ✗ JACOBIAN* jacobian = kinsolUserData->analyticJacobian;
639 ✗ assertStreamPrint(threadData, NULL != jacobian, "jacobian is NULL");
640 ✗ const SPARSE_PATTERN* sp = jacobian->sparsePattern;
641 ✗ assertStreamPrint(threadData, NULL != sp, "sp is NULL");
642 ✗ double *xScaling = NV_DATA_S(kinsolData->xScale);
643 long int column, nz;
644
645 ✗ if (SUNMatGetID(Jac) != SUNMATRIX_SPARSE || SM_SPARSETYPE_S(Jac) == SUN_CSR_MAT) {
646 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
647 "KINSOL: nlsSparseJac illegal input Jac. Matrix is not sparse!");
648 ✗ return -1;
649 }
650
651 /* performance measurement */
652 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
653
654 /* call generic sparse Jacobian with CSC buffer "SM_DATA_S(Jac)" */
655 ✗ evalJacobian(data, threadData, jacobian, NULL, SM_DATA_S(Jac), FALSE);
656 ✗ setSundialsSparsePattern(jacobian, Jac);
657
658 /* scaling */
659 ✗ if (kinsolData->nominalJac) {
660 ✗ for (column = 0; column < jacobian->sizeCols; column++) {
661 ✗ for (nz = sp->leadindex[column]; nz < sp->leadindex[column + 1]; nz++) {
662 ✗ SM_DATA_S(Jac)[nz] /= xScaling[column];
663 }
664 }
665 }
666
667 /* Finish sparse matrix and do a cheap check for singularity */
668 ✗ finishSparseColPtr(Jac, sp->nnz);
669
670 /* Debug print */
671 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC)) {
672 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "KINSOL: Sparse Matrix.");
673 ✗ SUNSparseMatrix_Print(Jac, stdout); /* TODO: Print in OMC_LOG_NLS_JAC */
674 ✗ nlsKinsolJacSumSparse(Jac);
675 ✗ messageClose(OMC_LOG_NLS_JAC);
676 }
677
678 ✗ if (omc_useStream[OMC_LOG_NLS_DERIVATIVE_TEST])
679 {
680 ✗ nlsKinsolDenseDerivativeTest(data, nlsData, kinsolData, Jac, KINSOL_JAC_EVAL);
681 }
682
683 ✗ if (omc_useStream[OMC_LOG_NLS_JAC_SUMS])
684 {
685 ✗ nlsJacobianRowColSums(data, nlsData, Jac, KINSOL_JAC_EVAL /* called at evaluation */, kinsolData->nominalJac /* scaled */);
686 }
687
688 /* performance measurement and statistics */
689 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
690 ✗ nlsData->numberOfJEval++;
691
692 ✗ return 0;
693 }
694
695 /**
696 * @brief Colored numeric Jacobian evaluation.
697 *
698 * Finite differences while using coloring of Jacobian.
699 * Jacobian matrix format has to be compressed sparse columns (CSC).
700 *
701 * @param vecX Input vector x.
702 * @param vecFX Vector for residual evaluation: f(x)
703 * @param Jac Jacobian to calculate: J(x)
704 * @param userData Pointer to user data, tpyecasted to `NLS_USERDATA`.
705 * @param tmp1 Work vector.
706 * @param tmp2 Work vector.
707 * @return int Return 0 on success.
708 */
709 ✗ static int nlsSparseJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
710 void *userData, N_Vector tmp1, N_Vector tmp2) {
711 /* Variables */
712 NLS_USERDATA *kinsolUserData;
713 DATA *data;
714 NONLINEAR_SYSTEM_DATA *nlsData;
715 NLS_KINSOL_DATA *kinsolData;
716 SPARSE_PATTERN *sparsePattern;
717
718 ✗ if (SUNMatGetID(Jac) != SUNMATRIX_SPARSE || SM_SPARSETYPE_S(Jac) == SUN_CSR_MAT) {
719 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
720 "KINSOL: nlsSparseJac illegal input Jac. Matrix is not sparse!");
721 ✗ return -1;
722 }
723
724 double *x;
725 double *fx;
726 double *xsave;
727 double *delta_hh;
728 double *xScaling;
729 double *fRes;
730
731 const double delta_h = sqrt(DBL_EPSILON * 2e1);
732
733 long int i, j, ii;
734 int nth;
735
736 /* Access userData and nonlinear system data */
737 kinsolUserData = (NLS_USERDATA *)userData;
738 ✗ data = kinsolUserData->data;
739 ✗ nlsData = kinsolUserData->nlsData;
740 ✗ kinsolData = (NLS_KINSOL_DATA *)nlsData->solverData;
741 ✗ sparsePattern = nlsData->sparsePattern;
742
743 /* Access N_Vector variables */
744 ✗ x = N_VGetArrayPointer(vecX);
745 ✗ fx = N_VGetArrayPointer(vecFX);
746 ✗ xsave = N_VGetArrayPointer(tmp1);
747 ✗ delta_hh = N_VGetArrayPointer(tmp2);
748 ✗ xScaling = NV_DATA_S(kinsolData->xScale);
749 ✗ fRes = NV_DATA_S(kinsolData->fRes);
750
751 nth = 0;
752
753 /* performance measurement */
754 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
755
756 /* reset matrix */
757 ✗ SUNMatZero(Jac);
758
759 /* Approximate Jacobian */
760 ✗ for (i = 0; i < sparsePattern->maxColors; i++) {
761 ✗ for (ii = 0; ii < kinsolData->size; ii++) {
762 ✗ if (sparsePattern->colorCols[ii] - 1 == i) {
763 ✗ xsave[ii] = x[ii];
764 ✗ delta_hh[ii] = delta_h * (fabs(xsave[ii]) + 1.0);
765 ✗ if ((xsave[ii] + delta_hh[ii] >= nlsData->max[ii])) {
766 ✗ delta_hh[ii] *= -1;
767 }
768 ✗ x[ii] += delta_hh[ii];
769
770 /* Calculate scaled difference quotient */
771 ✗ delta_hh[ii] = 1. / delta_hh[ii];
772 }
773 }
774 /* Evaluate residual function */
775 ✗ nlsKinsolResiduals(vecX, kinsolData->fRes, userData);
776
777 /* Save column in Jac and unset seed variables */
778 ✗ for (ii = 0; ii < kinsolData->size; ii++) {
779 ✗ if (sparsePattern->colorCols[ii] - 1 == i) {
780 ✗ nth = sparsePattern->leadindex[ii];
781 ✗ while (nth < sparsePattern->leadindex[ii + 1]) {
782 ✗ j = sparsePattern->index[nth];
783 ✗ if (kinsolData->nominalJac) {
784 ✗ setJacElementSundialsSparse(j, ii, nth, (fRes[j] - fx[j]) * delta_hh[ii] / xScaling[ii], Jac, SM_CONTENT_S(Jac)->M);
785 } else {
786 ✗ setJacElementSundialsSparse(j, ii, nth, (fRes[j] - fx[j]) * delta_hh[ii], Jac, SM_CONTENT_S(Jac)->M);
787 }
788 ✗ nth++;
789 }
790 ✗ x[ii] = xsave[ii];
791 }
792 }
793 }
794 /* Finish sparse matrix */
795 ✗ setSundialsSparseColPtrs(sparsePattern, Jac);
796 ✗ finishSparseColPtr(Jac, sparsePattern->nnz);
797
798 /* Debug print */
799 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC)) {
800 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "KINSOL: Sparse Matrix.");
801 ✗ SUNSparseMatrix_Print(Jac, stdout);
802 ✗ nlsKinsolJacSumSparse(Jac);
803 ✗ messageClose(OMC_LOG_NLS_JAC);
804 }
805 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_DEBUG)) {
806 ✗ sundialsPrintSparseMatrix(Jac, "A", OMC_LOG_JAC);
807 }
808
809 ✗ if (omc_useStream[OMC_LOG_NLS_DERIVATIVE_TEST])
810 {
811 ✗ nlsKinsolDenseDerivativeTest(data, nlsData, kinsolData, Jac, KINSOL_JAC_EVAL);
812 }
813
814 ✗ if (omc_useStream[OMC_LOG_NLS_JAC_SUMS])
815 {
816 ✗ nlsJacobianRowColSums(data, nlsData, Jac, KINSOL_JAC_EVAL, kinsolData->nominalJac /* scaled */);
817 }
818
819 /* performance measurement and statistics */
820 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
821 ✗ nlsData->numberOfJEval++;
822
823 ✗ return 0;
824 }
825
826 /**
827 * @brief Check for zero columns of matrix and print absolute sums.
828 *
829 * Compute absolute sum for each column and print the result.
830 * Report a warning if it is zero, since the matrix is singular in that case.
831 *
832 * @param A Dense matrix stored columnwise
833 */
834 ✗ static void nlsKinsolJacSumDense(SUNMatrix A) {
835 /* Variables */
836 int i, j;
837 double sum;
838
839 ✗ for (i = 0; i < SM_ROWS_D(A); ++i) {
840 sum = 0.0;
841 ✗ for (j = 0; j < SM_COLUMNS_D(A); ++j) {
842 ✗ sum += fabs(SM_ELEMENT_D(A, j, i));
843 }
844
845 ✗ if (sum == 0.0) { /* TODO: Don't check for equality(!), maybe use DBL_EPSILON */
846 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
847 "KINSOL: Column %d of Jacobian is zero. Jacobian is singular.",
848 i);
849 } else {
850 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "Column %d of Jacobian absolute sum = %g",
851 i, sum);
852 }
853 }
854 ✗ }
855
856 /**
857 * @brief Check for zero columns of matrix and print absolute sums.
858 *
859 * Compute absolute sum for each column and print the result.
860 * Report a warning if it is zero, since the matrix is singular in that case.
861 *
862 * @param A CSC matrix
863 */
864 ✗ static void nlsKinsolJacSumSparse(SUNMatrix A) {
865 /* Variables */
866 int i, j;
867 double sum;
868
869 /* Check format of A */
870 ✗ if (SM_SPARSETYPE_S(A) != SUN_CSC_MAT) {
871 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
872 "KINSOL: In function nlsKinsolJacSumSparse: Wrong sparse format "
873 "of SUNMatrix A.");
874 }
875
876 /* Check sums of each column of A */
877 ✗ for (i = 0; i < SM_COLUMNS_S(A); ++i) {
878 sum = 0.0;
879 ✗ for (j = SM_INDEXPTRS_S(A)[i]; j < SM_INDEXPTRS_S(A)[i + 1]; ++j) {
880 ✗ sum += fabs(SM_DATA_S(A)[j]);
881 }
882
883 ✗ if (sum == 0.0) { /* TODO: Don't check for equality(!), maybe use DBL_EPSILON */
884 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
885 "KINSOL: Column %d of Jacobian is zero. Jacobian is singular.",
886 i);
887 } else {
888 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "Column %d of Jacobian absolute sum = %g",
889 i, sum);
890 }
891 }
892 ✗ }
893
894 /**
895 * @brief Set maximum scaled length of Newton step.
896 *
897 * Will be set to the weighted Euclidean l_2 norm of xScale with maxstepfactor
898 * as weights. maxStep = sqrt(sum_{1=0}^{n-1} (xScale[i]*maxstepfactor)^2)
899 *
900 * @param kinsolData
901 * @param maxstepfactor
902 */
903 ✗ static void nlsKinsolSetMaxNewtonStep(NLS_KINSOL_DATA *kinsolData,
904 double maxstepfactor) {
905 /* Variables */
906 int flag;
907
908 ✗ N_VConst(maxstepfactor, kinsolData->fTmp);
909 ✗ kinsolData->mxnstepin = N_VWL2Norm(kinsolData->xScale, kinsolData->fTmp);
910
911 /* Set maximum step size */
912 ✗ flag = KINSetMaxNewtonStep(kinsolData->kinsolMemory, kinsolData->mxnstepin);
913 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetMaxNewtonStep");
914 ✗ }
915
916 /**
917 * @brief Set initial guess for KINSOL
918 *
919 * Depending on mode extrapolate start value or use old value for
920 * initialization.
921 *
922 * @param data
923 * @param kinsolData
924 * @param nlsData
925 * @param mode Has to be `INITIAL_EXTRAPOLATION` for extrapolation or
926 * `INITIAL_OLDVALUES` for using old values.
927 */
928 ✗ static void nlsKinsolResetInitial(DATA *data, NLS_KINSOL_DATA *kinsolData,
929 NONLINEAR_SYSTEM_DATA *nlsData,
930 initialMode mode) {
931 ✗ double *xStart = NV_DATA_S(kinsolData->initialGuess);
932
933 /* Set x vector */
934 ✗ switch (mode) {
935 ✗ case INITIAL_EXTRAPOLATION:
936 ✗ if (data->simulationInfo->discreteCall) {
937 ✗ memcpy(xStart, nlsData->nlsx, nlsData->size * (sizeof(double)));
938 } else {
939 ✗ memcpy(xStart, nlsData->nlsxExtrapolation,
940 ✗ nlsData->size * (sizeof(double)));
941 }
942 break;
943 ✗ case INITIAL_OLDVALUES:
944 ✗ memcpy(xStart, nlsData->nlsxOld, nlsData->size * (sizeof(double)));
945 break;
946 ✗ default:
947 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
948 "KINSOL: Function nlsKinsolResetInitial: Unknown mode %d.",
949 (int)mode);
950 }
951 ✗ }
952
953 /**
954 * @brief Scale x vector.
955 *
956 * Scale with 1.0 for mode `SCALING_ONES`.
957 * Scale with 1/fmax(nominal,|xStart|) for mode `SCALING_NOMINALSTART`.
958 *
959 * @param data unused
960 * @param kinsolData
961 * @param nlsData
962 * @param mode Mode for scaling. Use `SCALING_NOMINALSTART` for nominal
963 * scaling and `SCALING_ONES` for no scaling. Will be
964 * overwritten by simulation flag `FLAG_NO_SCALING`.
965 */
966 ✗ static void nlsKinsolXScaling(DATA *data, NLS_KINSOL_DATA *kinsolData,
967 NONLINEAR_SYSTEM_DATA *nlsData,
968 scalingMode mode) {
969 ✗ double *xStart = NV_DATA_S(kinsolData->initialGuess);
970 ✗ double *xScaling = NV_DATA_S(kinsolData->xScale);
971 int i;
972
973 /* if noScaling flag is used overwrite mode */
974 ✗ if (omc_flag[FLAG_NO_SCALING]) {
975 mode = SCALING_ONES;
976 }
977
978 /* Use nominal value or the actual working point for scaling */
979 ✗ switch (mode) {
980 case SCALING_NOMINALSTART:
981 ✗ for (i = 0; i < nlsData->size; i++) {
982 ✗ xScaling[i] = 1.0 / fmax(nlsData->nominal[i], fabs(xStart[i]));
983 }
984 break;
985 case SCALING_ONES:
986 ✗ for (i = 0; i < nlsData->size; i++) {
987 ✗ xScaling[i] = 1.0;
988 }
989 break;
990 ✗ case SCALING_JACOBIAN:
991 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
992 "KINSOL: Function nlsKinsolXScaling: Invalid mode SCALING_JACOBIAN.");
993 ✗ default:
994 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
995 "KINSOL: Function nlsKinsolXScaling: Unknown mode %d.", (int)mode);
996 }
997 ✗ }
998
999 /**
1000 * @brief Scale f(x) vector.
1001 *
1002 * @param data
1003 * @param kinsolData
1004 * @param nlsData
1005 * @param mode
1006 */
1007 ✗ static void nlsKinsolFScaling(DATA *data, NLS_KINSOL_DATA *kinsolData,
1008 NONLINEAR_SYSTEM_DATA *nlsData,
1009 scalingMode mode) {
1010 ✗ double *fScaling = NV_DATA_S(kinsolData->fScale);
1011 ✗ N_Vector x = kinsolData->initialGuess;
1012
1013 int i, j;
1014 SUNErrCode ret;
1015
1016 /* If noScaling flag is used overwrite mode */
1017 ✗ if (omc_flag[FLAG_NO_SCALING]) {
1018 mode = SCALING_ONES;
1019 }
1020
1021 /* Use nominal value or the actual working point for scaling */
1022 ✗ switch (mode) {
1023 ✗ case SCALING_JACOBIAN:
1024 /* Enable scaled jacobian evaluation */
1025 ✗ kinsolData->nominalJac = 1;
1026
1027 /* Calculate the scaled Jacobian */
1028 ✗ if (nlsData->sparsePattern && kinsolData->linearSolverMethod == NLS_LS_KLU) {
1029 ✗ if (kinsolData->solved != NLS_SOLVED) {
1030 ✗ kinsolData->nominalJac = 0;
1031 ✗ if (nlsData->analyticalJacobianColumn != NULL) {
1032 /* Calculate the sparse Jacobian symbolically */
1033 ✗ nlsSparseSymJac(x, kinsolData->fTmp, kinsolData->J, kinsolData->userData, NULL, NULL);
1034 } else {
1035 /* Update f(x) for the numerical jacobian matrix */
1036 ✗ nlsKinsolResiduals(x, kinsolData->fTmp, kinsolData->userData);
1037 ✗ nlsSparseJac(x, kinsolData->fTmp, kinsolData->J, kinsolData->userData, kinsolData->tmp1, kinsolData->tmp2);
1038 }
1039 }
1040 /* Scale the current Jacobian */
1041 ✗ SUNMatCopy_Sparse(kinsolData->J, kinsolData->scaledJ); /* Copy J into scaledJ */
1042 ✗ ret = _omc_SUNSparseMatrixVecScaling(kinsolData->scaledJ, kinsolData->xScale);
1043 ✗ if (ret != SUN_SUCCESS) {
1044 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "KINSOL: _omc_SUNSparseMatrixVecScaling failed.");
1045 }
1046 } else {
1047 /* Update f(x) for the numerical jacobian matrix */
1048 ✗ nlsKinsolResiduals(x, kinsolData->fTmp, kinsolData->userData);
1049 ✗ nlsDenseJac(nlsData->size, x, kinsolData->fTmp, kinsolData->J,
1050 kinsolData->userData, NULL, NULL);
1051 }
1052
1053 /* Disable scaled Jacobian evaluation */
1054 ✗ kinsolData->nominalJac = 0;
1055
1056 ✗ for (i = 0; i < nlsData->size; i++) {
1057 ✗ fScaling[i] = 1e-12;
1058 }
1059
1060 ✗ switch (SUNMatGetID(kinsolData->J))
1061 {
1062 case SUNMATRIX_SPARSE:
1063 ✗ for (i = 0; i < SM_NNZ_S(kinsolData->scaledJ); ++i) {
1064 ✗ if (fScaling[SM_INDEXVALS_S(kinsolData->scaledJ)[i]] < fabs(SM_DATA_S(kinsolData->scaledJ)[i])) {
1065 ✗ fScaling[SM_INDEXVALS_S(kinsolData->scaledJ)[i]] = fabs(SM_DATA_S(kinsolData->scaledJ)[i]);
1066 }
1067 }
1068 break;
1069 case SUNMATRIX_DENSE:
1070 ✗ for (i = 0; i < nlsData->size; i++) {
1071 ✗ for (j = 0; j < nlsData->size; j++) {
1072 ✗ if (fScaling[i] < fabs(SM_ELEMENT_D(kinsolData->J, j, i))) {
1073 ✗ fScaling[i] = fabs(SM_ELEMENT_D(kinsolData->J, j, i));
1074 }
1075 }
1076 }
1077 break;
1078 ✗ default:
1079 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1080 "KINSOL: Function nlsKinsolFScaling: Unknown matrix type.");
1081 }
1082
1083 /* inverse fScale */
1084 ✗ N_VInv(kinsolData->fScale, kinsolData->fScale);
1085
1086 ✗ break;
1087 case SCALING_ONES:
1088 ✗ for (i = 0; i < nlsData->size; i++) {
1089 ✗ fScaling[i] = 1.0;
1090 }
1091 break;
1092 ✗ case SCALING_NOMINALSTART:
1093 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1094 "KINSOL: Function nlsKinsolFScaling: Invalid mode SCALING_NOMINALSTART.");
1095 ✗ default:
1096 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1097 "KINSOL: Function nlsKinsolFScaling: Unknown mode %d.", (int)mode);
1098 }
1099 ✗ }
1100
1101 /**
1102 * @brief Print KINSOL configuration.
1103 *
1104 * Only prints if stream `LOG_NLS_V` is active.
1105 *
1106 * @param kinsolData
1107 * @param nlsData
1108 */
1109 ✗ static void nlsKinsolConfigPrint(NLS_KINSOL_DATA *kinsolData,
1110 NONLINEAR_SYSTEM_DATA *nlsData) {
1111 int retValue;
1112 double fNorm;
1113 ✗ DATA *data = kinsolData->userData->data;
1114 ✗ int eqSystemNumber = nlsData->equationIndex;
1115 _omc_vector vecStart, vecXScaling, vecFScaling;
1116
1117 ✗ if (!omc_useStream[OMC_LOG_NLS_V]) {
1118 ✗ return;
1119 }
1120
1121 ✗ _omc_initVector(&vecStart, kinsolData->size,
1122 ✗ NV_DATA_S(kinsolData->initialGuess));
1123 _omc_initVector(&vecXScaling, kinsolData->size,
1124 ✗ NV_DATA_S(kinsolData->xScale));
1125 _omc_initVector(&vecFScaling, kinsolData->size,
1126 ✗ NV_DATA_S(kinsolData->fScale));
1127
1128 ✗ if (eqSystemNumber>0) {
1129 ✗ _omc_printVectorWithEquationInfo(
1130 &vecStart, "Initial guess values", OMC_LOG_NLS_V,
1131 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber));
1132
1133 ✗ _omc_printVectorWithEquationInfo(
1134 &vecXScaling, "xScaling", OMC_LOG_NLS_V,
1135 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber));
1136 }
1137
1138 ✗ _omc_printVector(&vecFScaling, "fScaling", OMC_LOG_NLS_V);
1139
1140 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL F tolerance: %g", kinsolData->fnormtol);
1141 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL minimal step size %g",
1142 kinsolData->scsteptol);
1143 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL max iterations %d",
1144 ✗ 20 * kinsolData->size);
1145 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL strategy %d",
1146 kinsolData->kinsolStrategy);
1147 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL current retry %d", kinsolData->retries);
1148 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL max step %g", kinsolData->mxnstepin);
1149 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL linear solver %d",
1150 ✗ kinsolData->linearSolverMethod);
1151 }
1152
1153 /**
1154 * @brief Try to handle errors of KINSol().
1155 *
1156 * @param errorCode Error code from KINSOL.
1157 * @param data Pointer to data struct.
1158 * @param nlsData Non-linear solver data.
1159 * @param kinsolData Kinsol data.
1160 * @return modelica_boolean Return true, if it is possible to retry KINSol().
1161 */
1162 ✗ static modelica_boolean nlsKinsolErrorHandler(int errorCode, DATA *data,
1163 NONLINEAR_SYSTEM_DATA *nlsData,
1164 NLS_KINSOL_DATA *kinsolData) {
1165 int flag; /* KIN_* and KINLS_* codes, which are plain macros */
1166 SUNErrCode sunFlag; /* SUNLinearSolver codes, which are not */
1167 double fNorm;
1168 double *xStart = NV_DATA_S(kinsolData->initialGuess);
1169 double *xScaling = NV_DATA_S(kinsolData->xScale);
1170 long outL;
1171
1172 ✗ flag = KINSetNoInitSetup(kinsolData->kinsolMemory, SUNFALSE);
1173 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNoInitSetup");
1174
1175 ✗ switch (errorCode) {
1176 ✗ case KIN_MEM_NULL:
1177 ✗ throwStreamPrint(NULL, "KINSOL: Memory NULL ERROR %d\n", errorCode);
1178 return FALSE;
1179 break;
1180 ✗ case KIN_ILL_INPUT:
1181 ✗ throwStreamPrint(NULL, "KINSOL: Ill input ERROR %d\n", errorCode);
1182 return FALSE;
1183 break;
1184 ✗ case KIN_NO_MALLOC:
1185 ✗ throwStreamPrint(NULL, "KINSOL: Memory issue ERROR %d\n", errorCode);
1186 return FALSE;
1187 break;
1188 /* Just retry with new initial guess */
1189 ✗ case KIN_MXNEWT_5X_EXCEEDED:
1190 ✗ warningStreamPrint(
1191 OMC_LOG_NLS_V, 0,
1192 "Newton step exceed the maximum step size several times. Try again "
1193 "after increasing maximum step size.\n");
1194 ✗ kinsolData->maxstepfactor *= 1e5;
1195 ✗ nlsKinsolSetMaxNewtonStep(kinsolData, kinsolData->maxstepfactor);
1196 ✗ return TRUE;
1197 break;
1198 /* Just retry without line search */
1199 ✗ case KIN_LINESEARCH_NONCONV:
1200 ✗ warningStreamPrint(
1201 OMC_LOG_NLS_V, 0,
1202 "kinsols line search did not convergence. Try without.\n");
1203 ✗ kinsolData->kinsolStrategy = KIN_NONE;
1204 ✗ kinsolData->retries--;
1205 ✗ return TRUE;
1206 break;
1207 /* Maybe happened because of an out-dated factorization, so just retry */
1208 ✗ case KIN_LSOLVE_FAIL:
1209 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1210 "KINSOL: Matrix need new factorization. Try again.\n");
1211 ✗ if (kinsolData->linearSolverMethod == NLS_LS_KLU &&
1212 ✗ nlsData->sparsePattern) {
1213 /* Complete symbolic and numeric factorizations */
1214 ✗ sunFlag = SUNLinSol_KLUReInit(kinsolData->linSol, kinsolData->J,
1215 ✗ kinsolData->nnz, SUNKLU_REINIT_PARTIAL);
1216 ✗ checkReturnFlag_SUNDIALS(sunFlag, SUNDIALS_SUNLS_FLAG, "SUNLinSol_KLUReInit");
1217 ✗ return TRUE;
1218 }
1219 break;
1220 ✗ case KIN_MAXITER_REACHED:
1221 case KIN_REPTD_SYSFUNC_ERR:
1222 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1223 "KINSOL: Runs into issues retry with different configuration.\n");
1224 ✗ break;
1225 ✗ case KIN_LINIT_FAIL:
1226 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1227 "KINSOL: The linear solver's initialization function failed.\n");
1228 ✗ return errorCode;
1229 ✗ case KIN_LSETUP_FAIL:
1230 /* In case something goes wrong with the symbolic jacobian try the numerical */
1231 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0,
1232 "KINSOL: The kinls setup routine (lsetup) encountered an error. "
1233 "Retry with numerical Jacobian.\n");
1234 /* KLU always has a sparsity pattern (initKinsolMemory), and without an
1235 * analytic Jacobian it is numeric already */
1236 ✗ if (kinsolData->linearSolverMethod == NLS_LS_KLU && nlsData->analyticalJacobianColumn != NULL) {
1237 ✗ flag = KINSetJacFn(kinsolData->kinsolMemory, nlsSparseJac);
1238 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KINLS_FLAG, "KINSetJacFn");
1239 ✗ if (flag < 0) {
1240 return FALSE;
1241 }
1242 }
1243 break;
1244 /* the step got too small but the residual is not (checked by the caller) */
1245 case KIN_STEP_LT_STPTOL:
1246 break;
1247 ✗ case KIN_LINESEARCH_BCFAIL:
1248 ✗ KINGetNumBetaCondFails(kinsolData->kinsolMemory, &outL);
1249 ✗ warningStreamPrint(
1250 OMC_LOG_NLS_V, 0,
1251 "kinsols runs into issues with beta-condition fails: %ld\n", outL);
1252 ✗ break;
1253 ✗ default:
1254 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1255 "kinsol has a serious solving issue ERROR %d\n",
1256 errorCode);
1257 ✗ return FALSE;
1258 break;
1259 }
1260
1261 /* check if the current solution is sufficient anyway (a stalled step was checked already) */
1262 ✗ KINGetFuncNorm(kinsolData->kinsolMemory, &fNorm);
1263 ✗ if (errorCode != KIN_STEP_LT_STPTOL && fNorm < FTOL_WITH_LESS_ACCURACY) {
1264 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL: Move forward with a less accurate solution.");
1265 ✗ KINSetFuncNormTol(kinsolData->kinsolMemory, FTOL_WITH_LESS_ACCURACY);
1266 ✗ KINSetScaledStepTol(kinsolData->kinsolMemory, FTOL_WITH_LESS_ACCURACY);
1267 ✗ kinsolData->resetTol = TRUE;
1268 ✗ return TRUE;
1269 } else {
1270 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL: Current status of fx = %f", fNorm);
1271 }
1272
1273 /* reconfigure kinsol for another try */
1274 ✗ switch (kinsolData->retries) {
1275 ✗ case 0:
1276 /* try without scaling */
1277 ✗ nlsKinsolXScaling(data, kinsolData, nlsData, SCALING_ONES);
1278 ✗ nlsKinsolFScaling(data, kinsolData, nlsData, SCALING_ONES);
1279 ✗ break;
1280 ✗ case 1:
1281 /* try without line-search and oldValues */
1282 ✗ nlsKinsolResetInitial(data, kinsolData, nlsData, INITIAL_OLDVALUES);
1283 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH;
1284 ✗ break;
1285 ✗ case 2:
1286 /* try without line-search and oldValues */
1287 ✗ nlsKinsolResetInitial(data, kinsolData, nlsData, INITIAL_EXTRAPOLATION);
1288 ✗ kinsolData->kinsolStrategy = KIN_NONE;
1289 ✗ break;
1290 ✗ case 3:
1291 /* try with exact newton */
1292 ✗ nlsKinsolXScaling(data, kinsolData, nlsData, SCALING_NOMINALSTART);
1293 ✗ nlsKinsolFScaling(data, kinsolData, nlsData, SCALING_JACOBIAN);
1294 ✗ nlsKinsolResetInitial(data, kinsolData, nlsData, INITIAL_EXTRAPOLATION);
1295 ✗ KINSetMaxSetupCalls(kinsolData->kinsolMemory, 1);
1296 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH;
1297 ✗ break;
1298 ✗ case 4:
1299 /* try with exact newton to with out x scaling values */
1300 ✗ nlsKinsolXScaling(data, kinsolData, nlsData, SCALING_ONES);
1301 ✗ nlsKinsolFScaling(data, kinsolData, nlsData, SCALING_ONES);
1302 ✗ nlsKinsolResetInitial(data, kinsolData, nlsData, INITIAL_OLDVALUES);
1303 ✗ KINSetMaxSetupCalls(kinsolData->kinsolMemory, 1);
1304 ✗ kinsolData->kinsolStrategy = KIN_LINESEARCH;
1305 ✗ break;
1306 default:
1307 /* Too many retries */
1308 return FALSE;
1309 break;
1310 }
1311
1312 return TRUE;
1313 }
1314
1315 /**
1316 * @brief Solve non-linear system with KINSol
1317 *
1318 * @param data Runtime data struct.
1319 * @param threadData Thread data for error handling.
1320 * @param nlsData Pointer to non-linear system data.
1321 * @return NLS_SOLVER_STATUS Return NLS_SOLVED on success and NLS_FAILED otherwise.
1322 */
1323 /**
1324 * @brief Set sign constraints from the min and max attributes of the iteration variables.
1325 *
1326 * Only for variables whose initial guess already fulfills the constraint,
1327 * otherwise KINSol() rejects the initial guess.
1328 *
1329 * @param kinsolData Kinsol data with the initial guess.
1330 * @param nlsData Nonlinear system data with min and max values.
1331 */
1332 ✗ static void nlsKinsolSetConstraints(NLS_KINSOL_DATA *kinsolData, NONLINEAR_SYSTEM_DATA *nlsData) {
1333 int i, flag;
1334 ✗ double *x = NV_DATA_S(kinsolData->initialGuess);
1335 ✗ double *c = NV_DATA_S(kinsolData->constraints);
1336
1337 ✗ if (nlsData->min == NULL || nlsData->max == NULL) {
1338 return;
1339 }
1340
1341 ✗ for (i = 0; i < kinsolData->size; i++) {
1342 ✗ c[i] = 0.0;
1343 /* a variable on the bound would block every step that points outside */
1344 ✗ if (nlsData->min[i] >= 0.0 && x[i] > 0.0) {
1345 ✗ c[i] = nlsData->min[i] > 0.0 ? 2.0 : 1.0; /* x > 0 or x >= 0 */
1346 ✗ } else if (nlsData->max[i] <= 0.0 && x[i] < 0.0) {
1347 ✗ c[i] = nlsData->max[i] < 0.0 ? -2.0 : -1.0; /* x < 0 or x <= 0 */
1348 }
1349 }
1350 ✗ flag = KINSetConstraints(kinsolData->kinsolMemory, kinsolData->constraints);
1351 ✗ checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetConstraints");
1352 }
1353
1354 ✗ NLS_SOLVER_STATUS nlsKinsolSolve(DATA* data, threadData_t* threadData, NONLINEAR_SYSTEM_DATA* nlsData) {
1355
1356 ✗ NLS_KINSOL_DATA *kinsolData = (NLS_KINSOL_DATA *)nlsData->solverData;
1357 ✗ int eqSystemNumber = nlsData->equationIndex;
1358 ✗ int indexes[2] = {1, eqSystemNumber};
1359
1360 int flag;
1361 long nFEval;
1362 modelica_boolean success = FALSE;
1363 modelica_boolean retry = TRUE;
1364 modelica_boolean stalled;
1365 NLS_SOLVER_STATUS solver_status;
1366 ✗ double *xStart = NV_DATA_S(kinsolData->initialGuess);
1367 double fNormValue;
1368
1369 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_NLS_V, omc_dummyFileInfo, 1, indexes,
1370 "Start solving Non-Linear System %d (size %d) at time %g with Kinsol Solver",
1371 ✗ eqSystemNumber, (int) nlsData->size, data->localData[0]->timeValue);
1372
1373 /* Solve nonlinear system with KINSol() */
1374 ✗ kinsolData->retries = 0;
1375 do {
1376 ✗ nlsKinsolResetInitial(data, kinsolData, nlsData, INITIAL_EXTRAPOLATION);
1377
1378 /* Set x scaling */
1379 ✗ nlsKinsolXScaling(data, kinsolData, nlsData, SCALING_NOMINALSTART);
1380
1381 /* Set f scaling */
1382 ✗ nlsKinsolFScaling(data, kinsolData, nlsData, SCALING_JACOBIAN);
1383
1384 /* Set maximum step size */
1385 ✗ nlsKinsolSetMaxNewtonStep(kinsolData, kinsolData->maxstepfactor);
1386
1387 /* Keep the sign of variables with a non-negative min or non-positive max attribute */
1388 ✗ nlsKinsolSetConstraints(kinsolData, nlsData);
1389
1390 /* Dump configuration */
1391 ✗ nlsKinsolConfigPrint(kinsolData, nlsData);
1392
1393 /* TODO: This should be another flag, e.g. LOG_NLS_JAC_UPDATE and not OMC_LOG_NLS_DERIVATIVE_TEST
1394 only in some cases this derivative test makes sense, since the scaled Jacobian is outdated frequently!
1395 in most cases, we use an outdated jacobian here, such that errors explode and it detects wrong Jacobian mismatches
1396 that are due to the dense Jacobian evaluated at the new point x_new.
1397
1398 if (omc_useStream[OMC_LOG_NLS_DERIVATIVE_TEST])
1399 {
1400 nlsKinsolDenseDerivativeTest(data, nlsData, kinsolData, kinsolData->J, KINSOL_ENTRY_POINT);
1401 }
1402 */
1403
1404 ✗ if (omc_useStream[OMC_LOG_NLS_JAC_SUMS])
1405 {
1406 ✗ nlsJacobianRowColSums(data, nlsData, kinsolData->J, KINSOL_ENTRY_POINT /* called at entry point */, kinsolData->nominalJac /* scaled */);
1407 }
1408
1409 ✗ if (omc_useStream[OMC_LOG_NLS_SVD] || omc_useStream[OMC_LOG_NLS_SVD_V])
1410 {
1411 ✗ svd_compute(data, nlsData, SM_DATA_S(kinsolData->J), FALSE /* scaled */, KINSOL_ENTRY_POINT /* called at entry point */);
1412 }
1413
1414 ✗ flag = KINSol(
1415 kinsolData->kinsolMemory, /* KINSol memory block */
1416 kinsolData->initialGuess, /* initial guess on input; solution vector */
1417 kinsolData->kinsolStrategy, /* global strategy choice */
1418 kinsolData->xScale, /* scaling vector, for the variable cc */
1419 kinsolData->fScale); /* scaling vector for function values fval */
1420
1421 ✗ if (flag < 0 && kinsolData->attemptRetry) {
1422 ✗ warningStreamPrint(OMC_LOG_NLS, 0, "KINSol finished with errorCode %d.", flag);
1423 } else {
1424 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "KINSol finished with errorCode %d.", flag);
1425 }
1426 /* a step below the tolerance without any iteration only solves the system if the residual is small */
1427 stalled = FALSE;
1428 ✗ KINGetNumNonlinSolvIters(kinsolData->kinsolMemory, &nFEval);
1429 ✗ if (flag == KIN_STEP_LT_STPTOL && nFEval == 0) {
1430 /* KINGetFuncNorm is not set if no step was taken, evaluate the scaled residual */
1431 ✗ nlsKinsolResiduals(kinsolData->initialGuess, kinsolData->fRes, kinsolData->userData);
1432 ✗ fNormValue = N_VWL2Norm(kinsolData->fRes, kinsolData->fScale);
1433 ✗ stalled = !(fNormValue < FTOL_WITH_LESS_ACCURACY);
1434 ✗ if (stalled) {
1435 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "KINSOL: Step below tolerance but fx = %g is not small.", fNormValue);
1436 }
1437 }
1438
1439 /* Try to handle recoverable errors */
1440 ✗ retry = (flag < 0 || stalled) && kinsolData->attemptRetry && nlsKinsolErrorHandler(flag, data, nlsData, kinsolData);
1441
1442 /* solution found */
1443 ✗ if ((flag == KIN_SUCCESS) || (flag == KIN_INITIAL_GUESS_OK) ||
1444 ✗ (flag == KIN_STEP_LT_STPTOL && !stalled)) {
1445 success = TRUE;
1446 }
1447 ✗ kinsolData->retries++;
1448
1449 /* write statistics */
1450 ✗ KINGetNumNonlinSolvIters(kinsolData->kinsolMemory, &nFEval);
1451 ✗ nlsData->numberOfIterations += nFEval;
1452 ✗ nlsData->numberOfFEval = kinsolData->countResCalls;
1453
1454 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "Next try? success = %d, retry = %d, retries = %d = %s\n",
1455 success, retry, kinsolData->retries,
1456 ✗ !success && !retry && kinsolData->retries < RETRY_MAX ? "true" : "false");
1457 ✗ } while (!success && retry && kinsolData->retries < RETRY_MAX);
1458
1459 /* Check solution status */
1460 ✗ if (success && kinsolData->resetTol) {
1461 ✗ kinsolData->solved = NLS_SOLVED_LESS_ACCURACY;
1462 ✗ } else if (success) {
1463 ✗ kinsolData->solved = NLS_SOLVED;
1464 } else {
1465 ✗ kinsolData->solved = NLS_FAILED;
1466 }
1467
1468 /* Reset solver tolerance */
1469 ✗ if (kinsolData->resetTol) {
1470 ✗ KINSetFuncNormTol(kinsolData->kinsolMemory, kinsolData->fnormtol);
1471 ✗ KINSetScaledStepTol(kinsolData->kinsolMemory, kinsolData->scsteptol);
1472 ✗ kinsolData->resetTol = FALSE;
1473 }
1474
1475 ✗ if (success) {
1476 ✗ memcpy(nlsData->nlsx, xStart, nlsData->size * (sizeof(double)));
1477 }
1478
1479 ✗ messageClose(OMC_LOG_NLS_V);
1480
1481 ✗ return kinsolData->solved;
1482 }
1483
1484 #else /* WITH_SUNDIALS */
1485
1486 void* nlsKinsolAllocate(int size, void* userData, int attemptRetry, modelica_boolean isPatternAvailable) {
1487
1488 throwStreamPrint(NULL, "No sundials/kinsol support activated.");
1489 return 0;
1490 }
1491
1492 int nlsKinsolFree(void* kinsolData) {
1493
1494 throwStreamPrint(NULL, "No sundials/kinsol support activated.");
1495 return 0;
1496 }
1497
1498 int nlsKinsolSolve(void *data, threadData_t *threadData, void* nlsData) {
1499
1500 throwStreamPrint(threadData, "No sundials/kinsol support activated.");
1501 return 0;
1502 }
1503
1504 #endif /* WITH_SUNDIALS */
1505