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 / 647
Functions: 0.0% 0 / 0 / 39
Branches: 0.0% 0 / 0 / 478
OMCompiler/SimulationRuntime/c/simulation/solver/gbode_internal_nls.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 #include <klu.h>
29
30 #include "gbode_main.h"
31 #include "gbode_util.h"
32 #include "gbode_sparse.h"
33 #include "gbode_internal_nls.h"
34
35 #include "../options.h"
36 #include "../arrayIndex.h"
37
38 // TODO: Calibrate safety factor for internal tolerances
39
40 #define DBL_ABSORPTION (10 * DBL_EPSILON)
41 #define MAX(a,b) (((a)>(b))?(a):(b))
42
43 /* some constants for less verbose BLAS calls */
44 static const double DBL_ZERO = 0.0;
45 static const double DBL_ONE = 1.0;
46 static const double DBL_MINUS_ONE = -1.0;
47 static const int INT_ONE = 1;
48 static const int INT_TWO = 2;
49 static const char CHAR_TRANS = 'T';
50 static const char CHAR_NO_TRANS = 'N';
51
52 /* y := a * x + y */
53 extern void daxpy_(const int *n,
54 const double *alpha,
55 const double *x, const int *incX,
56 double *y, const int *incY);
57
58 /* y := alpha * A * x + beta * y */
59 extern void dgemv_(const char *trans,
60 const int *m,
61 const int *n,
62 const double *alpha, const double *A, const int *ldA,
63 const double *x, const int *incX,
64 const double *beta, double *y, const int *incY
65 );
66
67 /* C := alpha * A * B + beta * C */
68 extern void dgemm_(const char *transA,
69 const char *transB,
70 const int *m,
71 const int *n,
72 const int *k,
73 const double *alpha, const double *A, const int *ldA,
74 const double *B, const int *ldB,
75 const double *beta, double *C, const int *ldC
76 );
77
78 /* x := alpha * x */
79 extern void dscal_(const int *n,
80 const double *alpha,
81 double *x, const int *incX
82 );
83
84 /* y[incY * i] := x[incX * i] for i = 0 ... n - 1 */
85 extern void dcopy_(const int *n,
86 const double *x, const int *incX,
87 double *y, const int *incY
88 );
89
90 typedef struct KLUInternals
91 {
92 klu_common common; // shared KLU configuration for all transformed systems
93 klu_symbolic *symbolic; // shared symbolic factorization of the NLS sparse pattern
94 klu_numeric **num_real; // real numerical factorizations
95 klu_numeric **num_cmplx; // complex numerical factorizations
96 int n_real;
97 int n_cmplx;
98 } KLUInternals;
99
100 typedef struct GB_INTERNAL_NLS_DATA
101 {
102 NLS_USERDATA *nls_user_data; // pointer to data, gbode data, etc.
103 KLUInternals klu; // KLU data
104 double *jacobian_callback; // buffer for continuous ODE Jacobian (size = nnz(J_f))
105 int *ode_to_nls; // mapping ODE Jacobian nnz -> NLS Jacobian nnz
106 int *nls_diag_indices; // all diagonal nz indices of NLS Jacobian (size = cols)
107 double *scal; // scaling vector for termination of Newton loop
108 double *etas; // Newton contraction factors for each NLS stage (size == number of stages)
109 double eta_inital_damping; // Initial damping factor eta_new = eta_old^eta_initial_damping
110 double integrator_tol; // Integrator / user provided tolerance
111 double fnewt; // Newton tolerance: if eta * norm(dx) <= fnewt -> convergence
112 double theta_keep; // if norm(dx_k) / norm(dx_{k-1}) = theta_{k} < theta_keep -> keep old jacobian_callback
113 modelica_boolean call_jac; // call jacobian in the next call to NLS solve
114 double theta_divergence; // if norm(dx_k) / norm(dx_{k-1}) = theta_{k} > theta_divergence (<= 1.0) -> divergence of Newton
115 int max_newton_it; // maximum number of Newton iterations
116 int size; // size of the system
117 BUTCHER_TABLEAU *tabl; // butcher tableau of the method
118 modelica_boolean use_t_transform; // use T transform to solve the system (false for (E)SDIRK, true for FIRK)
119 double **real_nls_jacs; // real NLS jacobians
120 double **cmplx_nls_jacs; // complex NLS jacobians (packed as real, imag - memory layout is struct{double real, double imag}[])
121 double **real_nls_res; // real NLS residuum
122 double **cmplx_nls_res; // complex NLS residuum (packed as real, imag - memory layout is struct{double real, double imag}[])
123 double *Z; // update variables Z = Y(t_ij) - Y0 for T-transformation (coupled space)
124 double *W; // update variables W = (T^{-1} otimes I) Z for T-transformation (decoupled space)
125 double *work; // some work memory for the T transformation (size: transform->size * x.size) or other stuff, at least 32 * N_STATES bytes
126
127 // stuff for multirate
128 modelica_boolean multirate; // multirate or singlerate system?
129 modelica_boolean new_fast_states; // if the selection changed and we need to update sparse pattern and symbolic factorization - set from NLS routine
130 } GB_INTERNAL_NLS_DATA;
131
132 ✗ static inline SPARSE_PATTERN *getODEPattern(DATA *data, DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls)
133 {
134 ✗ return nls->multirate ? gbData->gbfData->sparsePattern_ODE : getJacobianCscPattern(getSymbolicOdeJacobian(data));
135 }
136
137 static inline SPARSE_PATTERN *getNLSPattern(DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls)
138 {
139 ✗ return nls->multirate ? gbData->gbfData->sparsePattern_NLS : gbData->sparsePattern_NLS;
140 }
141
142 ✗ static void gbInternal_evalJacobianMR(DATA* data,
143 threadData_t *threadData,
144 DATA_GBODE *gbData,
145 JACOBIAN* fullJac,
146 GB_INTERNAL_NLS_DATA *nls,
147 double* smallJac)
148 {
149 ✗ const SPARSE_PATTERN* smallSp = gbData->gbfData->sparsePattern_ODE;
150
151 ✗ int* fast_idx = gbData->fastStatesIdx;
152 ✗ unsigned int size_fast = gbData->nFastStates;
153
154 ✗ fullJac->evalSelection = NULL;
155 ✗ memset(fullJac->seedVars, 0, fullJac->sizeCols * sizeof(modelica_real));
156
157 int color, col, nz;
158
159 ✗ for (color = 0; color < smallSp->maxColors; color++)
160 {
161 ✗ for (col = 0; col < size_fast; col++)
162 {
163 ✗ unsigned int big_col = fast_idx[col];
164
165 ✗ if (smallSp->colorCols[col] - 1 == color)
166 {
167 ✗ fullJac->seedVars[big_col] = 1.0;
168 }
169 }
170
171 ✗ fullJac->evalColumn(data, threadData, fullJac, NULL);
172
173 ✗ for (col = 0; col < size_fast; col++)
174 {
175 ✗ unsigned int big_col = fast_idx[col];
176
177 ✗ if (smallSp->colorCols[col] - 1 == color)
178 {
179 ✗ for (nz = smallSp->leadindex[col]; nz < smallSp->leadindex[col + 1]; nz++)
180 {
181 ✗ unsigned int small_row = smallSp->index[nz];
182 ✗ unsigned int full_row = fast_idx[small_row];
183
184 ✗ smallJac[nz] = fullJac->resultVars[full_row];
185 }
186
187 ✗ fullJac->seedVars[big_col] = 0.0;
188 }
189 }
190 }
191
192 ✗ fullJac->evalSelection = NULL;
193 ✗ }
194
195 ✗ static void gbInternal_evalNumericalJacobian(DATA *data,
196 threadData_t *threadData,
197 DATA_GBODE *gbData,
198 GB_INTERNAL_NLS_DATA *nls)
199 {
200 ✗ const double delta_h = numericalDifferentiationDeltaXsolver;
201
202 ✗ const SPARSE_PATTERN *sparsity = getODEPattern(data, gbData, nls);
203 EVAL_SELECTION *selection = NULL;
204 int *state_map = NULL;
205 ✗ int size = gbData->nStates;
206 int full_size = gbData->nStates;
207
208 ✗ if (nls->multirate)
209 {
210 ✗ state_map = gbData->fastStatesIdx;
211 ✗ size = gbData->nFastStates;
212 ✗ selection = gbData->gbfData->evalSelectionFast;
213 }
214
215 ✗ double *x = data->localData[0]->realVars;
216 ✗ double *der_x = &data->localData[0]->realVars[full_size];
217
218 // work 0 ... full_size - 1 is used already by the backup data
219 ✗ double *der_x_ref = &nls->work[full_size];
220 ✗ double *x_save = &nls->work[2 * full_size];
221 ✗ double *delta_hh = &nls->work[3 * full_size];
222
223 ✗ const double *nominals = gbData->nominals;
224 ✗ const double *maxs = gbData->maxs;
225
226 memcpy(der_x_ref, der_x, full_size * sizeof(double));
227
228 ✗ for (unsigned int color = 0; color < sparsity->maxColors; color++)
229 {
230 // careful perturbation of the variables (a la DASSL interface)
231 ✗ for (unsigned int col = 0; col < size; col++)
232 {
233 ✗ unsigned int big_col = state_map ? state_map[col] : col;
234
235 ✗ if (sparsity->colorCols[col] - 1 == color)
236 {
237 // we follow the procedure of the DASSL interface for the selection of perturbation h_i
238
239 // h * f(x)_i
240 ✗ double delta_hhh = delta_h * der_x_ref[big_col];
241
242 // scal_raw = ATOL * NOMINAL + RTOL * abs(x_i), we use the real (un-transformed) integrator tolerances though
243 ✗ double raw_weight = nls->integrator_tol * nominals[big_col] + nls->integrator_tol * fabs(x[big_col]);
244
245 // choose h_i := h * max(abs(x_i), h * f(x)_i, ATOL * NOMINAL + RTOL * abs(x_i), 1e-3)
246 ✗ delta_hh[big_col] = delta_h * fmax(fmax(fmax(fabs(x[big_col]), 1e-3), fabs(delta_hhh)), fabs(raw_weight));
247 ✗ delta_hh[big_col] = x[big_col] + delta_hh[big_col] - x[big_col];
248
249 ✗ if (x[big_col] + delta_hh[big_col] >= maxs[big_col])
250 {
251 ✗ delta_hh[big_col] *= -1;
252 }
253
254 ✗ x_save[big_col] = x[big_col];
255 ✗ x[big_col] += delta_hh[big_col];
256 ✗ delta_hh[big_col] = 1.0 / delta_hh[big_col];
257 }
258 }
259
260 // eval f(x + h)
261 ✗ gbode_fODE(data, threadData, NULL, selection);
262
263 // do forward finite differencing (f(x + h) - f(x)) / h and reset states
264 ✗ for (unsigned int col = 0; col < size; col++)
265 {
266 ✗ unsigned int big_col = state_map ? state_map[col] : col;
267
268 ✗ if (sparsity->colorCols[col] - 1 == color)
269 {
270 ✗ for (unsigned int nz = sparsity->leadindex[col]; nz < sparsity->leadindex[col + 1]; nz++)
271 {
272 ✗ unsigned int small_row = sparsity->index[nz];
273 ✗ unsigned int big_row = state_map ? state_map[small_row] : small_row;
274
275 ✗ nls->jacobian_callback[nz] = (der_x[big_row] - der_x_ref[big_row]) * delta_hh[big_col];
276 }
277
278 ✗ x[big_col] = x_save[big_col];
279 }
280 }
281 }
282 ✗ }
283
284 /**
285 * @brief Calculate ODE Jacobian (numerical or analytic) depending on availability of
286 * analytic Jacobian.
287 *
288 * Fills the nls->jacobian_callback callback buffer with the ODE Jacobian.
289 * It requires the sparsity pattern and a work array of size 3 * #states in nls.
290 * This is just a wrapper around evalJacobian to also support numerical ODE Jacobians.
291 *
292 * @param data DATA object
293 * @param threadData Thread data
294 * @param gbData GBODE solver data
295 * @param nls Internal strategy data (nls->jacobian_callback will be filled)
296 */
297 ✗ static int gbInternal_evalJacobian(DATA *data, threadData_t *threadData, DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls)
298 {
299 int ret = -1;
300
301 /* try */
302 #if !defined(OMC_EMCC)
303 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
304 #endif
305
306 ✗ rt_tick(SIM_TIMER_JACOBIAN);
307 /* The multi-rate path drives a sub-set of the columns with its own coloring, so it
308 * always needs the forward Jacobian. All other paths use the Jacobian selected by
309 * the `-jacobian` flag, which may be evaluated forward, adjoint or bidirectionally.
310 * Its currently disabled so its fine but its the safe option */
311 ✗ JACOBIAN* jacobian_ODE = nls->multirate
312 ✗ ? &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A])
313 ✗ : getSymbolicOdeJacobian(data);
314
315 ✗ if (nls->multirate && jacobian_ODE->availability == JACOBIAN_AVAILABLE)
316 {
317 ✗ gbInternal_evalJacobianMR(data, threadData, gbData, jacobian_ODE, nls, nls->jacobian_callback);
318 }
319 ✗ else if (!nls->multirate && jacobian_ODE->availability == JACOBIAN_AVAILABLE)
320 {
321 ✗ evalJacobian(data, threadData, jacobian_ODE, NULL, nls->jacobian_callback, FALSE);
322 }
323 else
324 {
325 ✗ gbInternal_evalNumericalJacobian(data, threadData, gbData, nls);
326 }
327
328 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { ret = 0; }
329
330 #if !defined(OMC_EMCC)
331 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
332 #endif
333
334 ✗ if (nls->multirate)
335 {
336 ✗ gbData->gbfData->stats.nCallsJacobian++;
337 }
338 else
339 {
340 ✗ gbData->stats.nCallsJacobian++;
341 }
342
343 ✗ rt_accumulate(SIM_TIMER_JACOBIAN);
344
345 ✗ return ret;
346 }
347
348 /**
349 * @brief Assemble (E)SDIRK stage Jacobian for the nonlinear system.
350 *
351 * Scales the ODE Jacobian by `h * gamma`, maps it into the NLS Jacobian buffer,
352 * and subtracts the identity on the diagonal: J = -I + h*a_ii * dfdx.
353 *
354 * @param[in] data Runtime data (unused)
355 * @param[in] threadData Thread data
356 * @param[in] gbData GBODE integrator data
357 * @param[in] nls Internal NLS data with index mapping
358 * @param[in] jac_ode ODE Jacobian (structures)
359 * @param[in] jac_buf_ode ODE Jacobian values (already filled)
360 * @param[out] jac_buf_nls Output buffer for NLS Jacobian
361 * @return 0 on success
362 */
363 ✗ static int jacobian_DIRK_assemble(DATA *data,
364 threadData_t *threadData,
365 DATA_GBODE* gbData,
366 GB_INTERNAL_NLS_DATA *nls,
367 SPARSE_PATTERN *ode_jac_sp,
368 double *jac_buf_ode,
369 double *jac_buf_nls)
370 {
371 ✗ memset(jac_buf_nls, 0, getNLSPattern(gbData, nls)->nnz * sizeof(double));
372
373 ✗ DATA_GBODEF *gbfData = gbData->gbfData;
374
375 ✗ const double fac = (nls->multirate ? gbfData->stepSize * gbfData->tableau->A[gbfData->act_stage * gbfData->tableau->nStages + gbfData->act_stage]
376 ✗ : gbData->stepSize * gbData->tableau->A[gbData->act_stage * gbData->tableau->nStages + gbData->act_stage]);
377
378 ✗ for (int nz = 0; nz < ode_jac_sp->nnz; nz++)
379 {
380 ✗ int idx = nls->ode_to_nls[nz];
381 ✗ jac_buf_nls[idx] = fac * jac_buf_ode[nz];
382 }
383
384 ✗ for (int d = 0; d < nls->size; d++)
385 {
386 ✗ jac_buf_nls[nls->nls_diag_indices[d]] -= 1.0;
387 }
388
389 ✗ return 0;
390 }
391
392 /**
393 * @brief Assemble Jacobian for the nonlinear system.
394 *
395 * Jacobian has form gamma / h * I - J_f
396 *
397 * @param[in] data Runtime data (unused)
398 * @param[in] threadData Thread data
399 * @param[in] gbData GBODE integrator data
400 * @param[in] nls Internal NLS data with index mapping
401 * @param[in] jac_ode ODE Jacobian (structures)
402 * @param[in] jac_buf_ode ODE Jacobian values
403 * @param[out] jac_buf_nls Output buffer for NLS Jacobian
404 * @return 0 on success
405 */
406 ✗ static int jacobian_real_assemble(DATA *data,
407 threadData_t *threadData,
408 DATA_GBODE* gbData,
409 GB_INTERNAL_NLS_DATA *nls,
410 double gamma,
411 SPARSE_PATTERN *ode_jac_sp,
412 double *jac_buf_ode,
413 double *jac_buf_nls)
414 {
415 ✗ memset(jac_buf_nls, 0, getNLSPattern(gbData, nls)->nnz * sizeof(double));
416
417 ✗ const double inv_step = 1.0 / (nls->multirate ? gbData->gbfData->stepSize : gbData->stepSize);
418 ✗ const double weight = inv_step * gamma;
419
420 ✗ for (int nz = 0; nz < ode_jac_sp->nnz; nz++)
421 {
422 ✗ int idx = nls->ode_to_nls[nz];
423 ✗ jac_buf_nls[idx] = -jac_buf_ode[nz];
424 }
425
426 ✗ for (int d = 0; d < nls->size; d++)
427 {
428 ✗ jac_buf_nls[nls->nls_diag_indices[d]] += weight;
429 }
430
431 ✗ return 0;
432 }
433
434 /**
435 * @brief Assemble complex Jacobian for the nonlinear system.
436 *
437 * Jacobian has form (alpha + i * beta) / h * I - J_f
438 *
439 * @param[in] data Runtime data (unused)
440 * @param[in] threadData Thread data
441 * @param[in] gbData GBODE integrator data (step size, tableau, stage)
442 * @param[in] nls Internal NLS data with index mapping
443 * @param[in] alpha Real part of complex weight
444 * @param[in] beta Imaginary part of complex weight
445 * @param[in] jac_ode ODE Jacobian (structures)
446 * @param[in] jac_buf_ode ODE Jacobian values
447 * @param[out] jac_buf_nls Output buffer for NLS Jacobian
448 * @return 0 on success
449 */
450 ✗ static int jacobian_cmplx_assemble(DATA *data,
451 threadData_t *threadData,
452 DATA_GBODE* gbData,
453 GB_INTERNAL_NLS_DATA *nls,
454 double alpha,
455 double beta,
456 SPARSE_PATTERN *ode_jac_sp,
457 double *jac_buf_ode,
458 double *jac_buf_nls)
459 {
460 ✗ memset(jac_buf_nls, 0, 2 * getNLSPattern(gbData, nls)->nnz * sizeof(double));
461
462 ✗ const double inv_step = 1.0 / (nls->multirate ? gbData->gbfData->stepSize : gbData->stepSize);
463 ✗ const double weight_real = inv_step * alpha;
464 ✗ const double weight_imag = inv_step * beta;
465
466 ✗ for (int nz = 0; nz < ode_jac_sp->nnz; nz++)
467 {
468 ✗ int idx = nls->ode_to_nls[nz];
469 ✗ jac_buf_nls[2 * idx] = -jac_buf_ode[nz];
470 }
471
472 ✗ for (int d = 0; d < nls->size; d++)
473 {
474 ✗ jac_buf_nls[2 * nls->nls_diag_indices[d]] += weight_real;
475 ✗ jac_buf_nls[2 * nls->nls_diag_indices[d] + 1] += weight_imag;
476 }
477
478 ✗ return 0;
479 }
480
481 ✗ static void gbInternal_KLU_initialize(KLUInternals *klu, int n_real, int n_cmplx)
482 {
483 ✗ assertStreamPrint(NULL, n_real >= 0 && n_cmplx >= 0, "Invalid number of KLU systems: %d real and %d complex.", n_real, n_cmplx);
484 ✗ klu_defaults(&klu->common);
485 ✗ klu->symbolic = NULL;
486 ✗ klu->n_real = n_real;
487 ✗ klu->n_cmplx = n_cmplx;
488 ✗ klu->num_real = n_real ? (klu_numeric **) calloc(n_real, sizeof(klu_numeric *)) : NULL;
489 ✗ klu->num_cmplx = n_cmplx ? (klu_numeric **) calloc(n_cmplx, sizeof(klu_numeric *)) : NULL;
490 ✗ }
491
492 ✗ static void gbInternal_KLU_freeNumerics(KLUInternals *klu)
493 {
494 ✗ for (int sys = 0; sys < klu->n_real; sys++)
495 {
496 ✗ if (klu->num_real[sys])
497 {
498 ✗ klu_free_numeric(&klu->num_real[sys], &klu->common);
499 }
500 }
501 ✗ for (int sys = 0; sys < klu->n_cmplx; sys++)
502 {
503 ✗ if (klu->num_cmplx[sys])
504 {
505 ✗ klu_free_numeric(&klu->num_cmplx[sys], &klu->common);
506 }
507 }
508 ✗ }
509
510 /** @brief Run one shared symbolic analysis for all KLU systems. */
511 ✗ static int gbInternal_KLU_analyze(KLUInternals *klu, int size, int *Ap, int *Ai)
512 {
513 ✗ klu_defaults(&klu->common);
514 ✗ klu->symbolic = klu_analyze(size, Ap, Ai, &klu->common);
515 ✗ if (klu->common.status < 0)
516 {
517 ✗ throwStreamPrint(NULL, "Error in gbInternal_KLU_analyze. Symbolic analysis with KLU failed.");
518 }
519 ✗ return klu->common.status;
520 }
521
522 ✗ static int gbInternal_KLU_reanalyze(KLUInternals *klu, int size, int *Ap, int *Ai)
523 {
524 ✗ gbInternal_KLU_freeNumerics(klu);
525 ✗ if (klu->symbolic)
526 {
527 ✗ klu_free_symbolic(&klu->symbolic, &klu->common);
528 }
529 ✗ return gbInternal_KLU_analyze(klu, size, Ap, Ai);
530 }
531
532 ✗ static void gbInternal_KLU_free(KLUInternals *klu)
533 {
534 ✗ gbInternal_KLU_freeNumerics(klu);
535 ✗ if (klu->symbolic)
536 {
537 ✗ klu_free_symbolic(&klu->symbolic, &klu->common);
538 }
539 ✗ free(klu->num_real);
540 ✗ free(klu->num_cmplx);
541 ✗ }
542
543 /** @brief Perform or update real-valued KLU numeric factorization. */
544 ✗ static int gbInternal_dKLU_factorize(KLUInternals *klu, int system, int *Ap, int *Ai, double *values)
545 {
546 ✗ assertStreamPrint(NULL, system >= 0 && system < klu->n_real, "Invalid real KLU system index %d.", system);
547 ✗ klu_numeric **numeric = &klu->num_real[system];
548 ✗ if (*numeric)
549 {
550 ✗ klu_refactor(Ap, Ai, values, klu->symbolic, *numeric, &klu->common);
551 }
552 else
553 {
554 ✗ *numeric = klu_factor(Ap, Ai, values, klu->symbolic, &klu->common);
555 }
556 ✗ return klu->common.status;
557 }
558
559 /** @brief Solve a real linear system using KLU. */
560 ✗ static int gbInternal_dKLU_solve(KLUInternals *klu, int system, int size, double *rhs)
561 {
562 ✗ assertStreamPrint(NULL, system >= 0 && system < klu->n_real, "Invalid real KLU system index %d.", system);
563 ✗ return klu_solve(klu->symbolic, klu->num_real[system], size, 1, rhs, &klu->common);
564 }
565
566 /** @brief Perform or update complex-valued KLU numeric factorization (values packed as struct {double real, double imag}). */
567 ✗ static int gbInternal_zKLU_factorize(KLUInternals *klu, int system, int *Ap, int *Ai, double *values)
568 {
569 ✗ assertStreamPrint(NULL, system >= 0 && system < klu->n_cmplx, "Invalid complex KLU system index %d.", system);
570 ✗ klu_numeric **numeric = &klu->num_cmplx[system];
571 ✗ if (*numeric)
572 {
573 ✗ klu_z_refactor(Ap, Ai, values, klu->symbolic, *numeric, &klu->common);
574 }
575 else
576 {
577 ✗ *numeric = klu_z_factor(Ap, Ai, values, klu->symbolic, &klu->common);
578 }
579 ✗ return klu->common.status;
580 }
581
582 /** @brief Solve a complex linear system using KLU (values packed as struct {double real, double imag}). */
583 ✗ static int gbInternal_zKLU_solve(KLUInternals *klu, int system, int size, double *rhs)
584 {
585 ✗ assertStreamPrint(NULL, system >= 0 && system < klu->n_cmplx, "Invalid complex KLU system index %d.", system);
586 ✗ return klu_z_solve(klu->symbolic, klu->num_cmplx[system], size, 1, rhs, &klu->common);
587 }
588
589 /** @brief Create scalings for scaled 2-norms: used for Newton convergence and integration acceptance criteria. */
590 ✗ static void createGbScales(GB_INTERNAL_NLS_DATA *nls, DATA_GBODE *gbData, double *y1, double *y2)
591 {
592 ✗ const double *nominals = gbData->nominals;
593
594 ✗ if (!nls->multirate)
595 {
596 ✗ for (int i = 0; i < nls->size; i++)
597 {
598 ✗ nls->scal[i] = 1.0 / (nls->integrator_tol * nominals[i] + fmax(fabs(y1[i]), fabs(y2[i])) * nls->integrator_tol);
599 }
600 }
601 else
602 {
603 ✗ for (int i = 0; i < nls->size; i++)
604 {
605 ✗ const size_t fast_idx = (size_t) gbData->fastStatesIdx[i];
606 ✗ nls->scal[i] = 1.0 / (nls->integrator_tol * nominals[fast_idx] + fmax(fabs(y1[i]), fabs(y2[i])) * nls->integrator_tol);
607 }
608 }
609 ✗ }
610
611 /** @brief Compute scaled norm of possible vector stack vec = (v1, v2, ..., v_{stacksize}). */
612 ✗ static double gbScalesNorm(GB_INTERNAL_NLS_DATA *nls, double *vec, int stack_size)
613 {
614 double sum = 0.0;
615 ✗ for (int j = 0; j < stack_size; j++)
616 {
617 ✗ double *vec_stride = &vec[j * nls->size];
618 ✗ for (int i = 0; i < nls->size; i++)
619 {
620 ✗ double tmp = vec_stride[i] * nls->scal[i];
621 ✗ sum += tmp * tmp;
622 }
623 }
624
625 ✗ return sqrt(sum / ((double)nls->size * (double)stack_size));
626 }
627
628 /** @brief Compute scaled norm of possible vector stack vec = (x + v1, x + v2, ..., x + v_{stacksize}). */
629 ✗ static double gbScalesNormXPlusZ(GB_INTERNAL_NLS_DATA *nls, double *x, double *z, int stack_size)
630 {
631 double sum = 0.0;
632 ✗ for (int j = 0; j < stack_size; j++)
633 {
634 ✗ double *vec_stride = &z[j * nls->size];
635 ✗ for (int i = 0; i < nls->size; i++)
636 {
637 ✗ double tmp = (x[i] + vec_stride[i]) * nls->scal[i];
638 ✗ sum += tmp * tmp;
639 }
640 }
641
642 ✗ return sqrt(sum / ((double)nls->size * (double)stack_size));
643 }
644
645 ✗ static int gbInternalEvaluateSimplifiedJacobian(DATA *data,
646 threadData_t *threadData,
647 DATA_GBODE* gbData,
648 GB_INTERNAL_NLS_DATA *nls,
649 modelica_boolean *jac_called,
650 modelica_boolean isDIRK)
651 {
652 int ret;
653 double time_backup;
654 ✗ int full_size = gbData->nStates;
655 ✗ SOLVERSTATS *stats = (nls->multirate ? &gbData->gbfData->stats : &gbData->stats);
656 ✗ double *y_eval = (nls->multirate ? gbData->gbfData->yOld : gbData->yOld);
657
658 // we need to backup potential interpolation data (e.g. in SDIRK realVars already contains Y_full(t0 + h * c1))
659 // backup the time and all the states and apply them after the simplified Jacobian computation
660 ✗ if (isDIRK)
661 {
662 ✗ time_backup = data->localData[0]->timeValue;
663 ✗ if (nls->multirate) memcpy(nls->work, data->localData[0]->realVars, full_size * sizeof(double));
664 }
665
666 // set values for known last point y_eval (simplified Newton)
667 ✗ memcpy(data->localData[0]->realVars, y_eval, full_size * sizeof(double));
668 ✗ data->localData[0]->timeValue = (nls->multirate ? gbData->gbfData->time : gbData->time);
669
670 // callback ODE + callback Jacobian of ODE -> nls_jacobian buffer
671 ✗ ret = gbode_fODE(data, threadData, &(stats->nCallsODE), nls->multirate ? gbData->gbfData->evalSelectionFast : NULL); // TODO: is this correct?
672 ✗ if (ret < 0) return ret;
673
674 ✗ ret = gbInternal_evalJacobian(data, threadData, gbData, nls);
675 ✗ if (ret < 0) return ret;
676
677 ✗ *jac_called = TRUE;
678
679 // reapply the interpolated data
680 ✗ if (isDIRK)
681 {
682 ✗ data->localData[0]->timeValue = time_backup;
683 ✗ if (nls->multirate) memcpy(data->localData[0]->realVars, nls->work, full_size * sizeof(double));
684 }
685
686 return 0;
687 }
688
689 /** @brief Solve one stage of a DIRK method with the internal solve routine. */
690 ✗ static NLS_SOLVER_STATUS gbInternalSolveNls_DIRK(DATA *data,
691 threadData_t *threadData,
692 NONLINEAR_SYSTEM_DATA* nonlinsys,
693 DATA_GBODE* gbData,
694 GB_INTERNAL_NLS_DATA *nls)
695 {
696 ✗ int size = nls->size;
697 ✗ int stage = (nls->multirate ? gbData->gbfData->act_stage : gbData->act_stage);
698 ✗ double *x = nonlinsys->nlsx;
699 ✗ double *x_start = nonlinsys->nlsxOld; // currently the extrapolated (e.g. dense output / hermite guess)
700 ✗ double *res = nonlinsys->resValues;
701 ✗ BUTCHER_TABLEAU *tabl = (nls->multirate ? gbData->gbfData->tableau : gbData->tableau);
702
703 ✗ double stepSize = (nls->multirate ? gbData->gbfData->stepSize : gbData->stepSize);
704 ✗ double lastStepSize = (nls->multirate ? gbData->gbfData->lastStepSize : gbData->lastStepSize);
705
706 // return code from KLU or RHS / Jac calls
707 int ret;
708
709 ✗ createGbScales(nls, gbData, x, x_start);
710 double *scal = nls->scal;
711
712 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=(nls->multirate ? (void *) gbData->gbfData : (void *) gbData)};
713 ✗ SPARSE_PATTERN *ode_pattern = getODEPattern(data, gbData, nls);
714 SPARSE_PATTERN *nls_pattern = getNLSPattern(gbData, nls);
715
716 ✗ const int flag = 1;
717 ✗ modelica_boolean jac_called = FALSE;
718
719 ✗ modelica_boolean is_esdirk = (tabl->A[0] == 0.0);
720 ✗ modelica_boolean sdirk_first_stage = (stage == 0 && !is_esdirk);
721 ✗ modelica_boolean esdirk_first_stage = (stage == 1 && is_esdirk);
722
723 ✗ if (sdirk_first_stage || esdirk_first_stage)
724 {
725 ✗ if (nls->call_jac || gbData->eventHappened)
726 {
727 ✗ ret = gbInternalEvaluateSimplifiedJacobian(data, threadData, gbData, nls, &jac_called, TRUE);
728 ✗ if (ret < 0) return NLS_FAILED;
729 }
730
731 ✗ if (jac_called || stepSize != lastStepSize)
732 {
733 /* fill NLS Jacobian as h * gamma * J_f - I, where J_f is old or newly computed ODE Jacobian nls->jacobian_callback */
734 ✗ jacobian_DIRK_assemble(data, threadData, gbData, nls, ode_pattern, nls->jacobian_callback, nls->real_nls_jacs[0]);
735
736 /* perform factorization */
737 ✗ ret = gbInternal_dKLU_factorize(&nls->klu, 0, (int *) nls_pattern->leadindex, (int *) nls_pattern->index, nls->real_nls_jacs[0]);
738 ✗ if (ret < 0) return NLS_FAILED;
739 }
740 }
741
742 /* invalidate eta of stage if an event happened */
743 ✗ if (gbData->eventHappened) nls->etas[stage] = DBL_MAX;
744
745 ✗ memcpy(x, x_start, size * sizeof(double));
746
747 // norms, convergence rate
748 double nrm_x = 0;
749 double nrm_delta = 0;
750 double nrm_delta_prev = 0;
751 double theta = 0;
752
753 // Newton iteration count - we start with newt_it = 1, because we need this for the step size selection and conditions below
754 ✗ for (int newt_it = 1 ;; newt_it++)
755 ✗ {
756 ✗ nonlinsys->residualFunc(&resUserData, x, res, &flag);
757
758 ✗ ret = gbInternal_dKLU_solve(&nls->klu, 0, size, res);
759 ✗ if (ret < 0) return NLS_FAILED;
760 ✗ daxpy_(&size, &DBL_MINUS_ONE, res, &INT_ONE, x, &INT_ONE);
761
762 ✗ nrm_delta_prev = fmax(DBL_EPSILON, nrm_delta);
763 ✗ nrm_delta = gbScalesNorm(nls, res, 1);
764
765 // handle absorption effects
766 ✗ nrm_x = gbScalesNorm(nls, x, 1);
767 ✗ modelica_boolean absorption = (nrm_delta <= DBL_ABSORPTION * nrm_x);
768
769 ✗ if (newt_it > 1)
770 {
771 ✗ theta = nrm_delta / nrm_delta_prev;
772
773 // Newton failed -> divergence
774 ✗ if (theta >= nls->theta_divergence && !absorption)
775 {
776 break;
777 }
778
779 ✗ nls->etas[stage] = theta / (1 - theta);
780 }
781 else
782 {
783 ✗ nls->etas[stage] = pow(fmax(nls->etas[stage], DBL_EPSILON), nls->eta_inital_damping);
784 }
785
786 ✗ if (!isfinite(nls->etas[stage]) || !isfinite(nrm_delta))
787 {
788 // Inf or NaN detected
789 // Either RHS or Jacobian or solution of the system contained a Inf or NaN
790 return NLS_FAILED;
791 }
792
793 // Newton converged
794 ✗ if (nls->etas[stage] * nrm_delta < nls->fnewt || absorption)
795 {
796 ✗ if (theta < nls->theta_keep)
797 {
798 ✗ nls->call_jac = FALSE;
799 }
800 else
801 {
802 ✗ nls->call_jac = TRUE;
803 }
804
805 ✗ return NLS_SOLVED;
806 }
807
808 // Newton failed -> iteration limit exceeded or too slow convergence
809 ✗ if (newt_it == nls->max_newton_it || (pow(theta, nls->max_newton_it - newt_it) / (1 - theta) * nrm_delta > nls->fnewt))
810 {
811 break;
812 }
813 }
814
815 ✗ nls->call_jac = TRUE;
816 ✗ return NLS_FAILED;
817 }
818
819 /** @brief Compute (T otimes I) * v for block vectors (applies T to block_count blocks of size block_size). */
820 ✗ static inline void dense_kron_id_vec(int block_count,
821 int block_size,
822 const double *T,
823 const double *v,
824 double *out)
825 {
826 ✗ dgemm_(
827 &CHAR_NO_TRANS, &CHAR_NO_TRANS,
828 &block_size, &block_count, &block_count,
829 &DBL_ONE,
830 v, &block_size,
831 T, &block_count,
832 &DBL_ZERO,
833 out, &block_size
834 );
835 ✗ }
836
837 /**
838 * @brief Multiply a stacked vector by a scaled block-lower triangular or block-diagonal matrix
839 * @par Runtime: O(m * n) for block diagonal -- O(m^2 * n) for fully dense block-lower-triangular
840 *
841 * out += factor * ((Lambda + L) otimes I_n) * v.
842 *
843 * Vector layout is fixed by the transform:
844 * - first nRealBlocks scalar real rows,
845 * - then nComplexBlocks consecutive 2x2 real blocks.
846 *
847 * Lambda is the block diagonal part. Each diagonal block is either:
848 * - 1x1 real row: out_i += factor * gamma[i] * v_i
849 *
850 * - 2x2 complex block: out_i += factor * alpha[j] * v_i - factor * beta[j] * v_{i+1}
851 * out_i+1 += factor * beta[j] * v_i + factor * alpha[j] * v_{i+1}
852 *
853 * where i and j are given from the eigenvalue indices (realEigenvalueIndex, complexEigenvalueIndex).
854 *
855 * L contains only the strict lower triangular couplings outside these diagonal blocks.
856 * Hence L never stores the lower entry inside a complex 2x2 block; that entry is beta[j].
857 * Rows with hasL[row] == FALSE are skipped completely.
858 */
859 ✗ static void scaled_transform_matvec(T_TRANSFORM *transform,
860 int block_size,
861 const double factor,
862 const double *v,
863 double *out)
864 {
865 // Diagonal / Block-Diagonal part:
866
867 // 1x1 real blocks: out_i += factor * gamma_i * v_i
868 ✗ for (int real_row = 0; real_row < transform->nRealBlocks; real_row++)
869 {
870 ✗ int sys = transform->realEigenvalueIndex[real_row];
871 ✗ double a = factor * transform->gamma[sys];
872 ✗ daxpy_(&block_size, &a, &v[real_row * block_size], &INT_ONE, &out[real_row * block_size], &INT_ONE);
873 }
874
875 ✗ int offset = transform->nRealBlocks * block_size;
876
877 // 2x2 blocks: out[j] += factor * (alpha_j * v[j] - beta_j * v[j+1])
878 // out[j+1] += factor * (beta_j * v[j] + alpha_j * v[j+1])
879 ✗ for (int cmplx_block = 0; cmplx_block < transform->nComplexBlocks; cmplx_block++)
880 {
881 ✗ int sys = transform->complexEigenpairIndex[cmplx_block];
882 ✗ double a = factor * transform->alpha[sys];
883 ✗ double b = factor * transform->beta[sys];
884 ✗ double mb = -b;
885
886 ✗ const double *v0 = &v[offset];
887 ✗ const double *v1 = &v[offset + block_size];
888
889 ✗ double *out0 = &out[offset];
890 ✗ double *out1 = &out[offset + block_size];
891
892 // out0 = a*v0 - b*v1
893 ✗ daxpy_(&block_size, &a, v0, &INT_ONE, out0, &INT_ONE); // out0 += a*v0
894 ✗ daxpy_(&block_size, &mb, v1, &INT_ONE, out0, &INT_ONE); // out0 -= b*v1
895
896 // out1 = a*v1 + b*v0
897 ✗ daxpy_(&block_size, &a, v1, &INT_ONE, out1, &INT_ONE); // out1 += a*v1
898 ✗ daxpy_(&block_size, &b, v0, &INT_ONE, out1, &INT_ONE); // out1 += b*v0
899
900 ✗ offset += 2 * block_size;
901 }
902
903 // Strictly lower triangular part:
904
905 // L rows for real blocks: out[row] += factor * L[row,col] * v[col], col < row
906 ✗ for (int row = 1; row < transform->nRealBlocks; row++)
907 {
908 ✗ if (!transform->hasL[row]) continue;
909
910 ✗ double *out_row = &out[row * block_size];
911
912 ✗ for (int col = 0; col < row; col++)
913 {
914 ✗ double a = factor * transform->L[GBODE_L_INDEX(row, col)];
915 ✗ if (a != 0.0) daxpy_(&block_size, &a, &v[col * block_size], &INT_ONE, out_row, &INT_ONE);
916 }
917 }
918
919 // L rows for complex blocks: out[i] += factor * L[i,col] * v[col], col < i
920 // out[i+1] += factor * L[i+1,col] * v[col], col < i
921 int cmplx_row = transform->nRealBlocks;
922 ✗ for (int cmplx_block = 0; cmplx_block < transform->nComplexBlocks; cmplx_block++)
923 {
924 int row0 = cmplx_row;
925 ✗ int row1 = cmplx_row + 1;
926
927 ✗ double *out0 = &out[row0 * block_size];
928 ✗ double *out1 = &out[row1 * block_size];
929
930 ✗ if (transform->hasL[row0])
931 {
932 ✗ for (int col = 0; col < cmplx_row; col++)
933 {
934 ✗ double a = factor * transform->L[GBODE_L_INDEX(row0, col)];
935 ✗ if (a != 0.0) daxpy_(&block_size, &a, &v[col * block_size], &INT_ONE, out0, &INT_ONE);
936 }
937 }
938
939 ✗ if (transform->hasL[row1])
940 {
941 ✗ for (int col = 0; col < cmplx_row; col++)
942 {
943 ✗ double a = factor * transform->L[GBODE_L_INDEX(row1, col)];
944 ✗ if (a != 0.0) daxpy_(&block_size, &a, &v[col * block_size], &INT_ONE, out1, &INT_ONE);
945 }
946 }
947
948 ✗ cmplx_row += 2;
949 }
950 ✗ }
951
952 #define GB_INTERNAL_LEFT_BOUNDARY -1
953
954 /**
955 * @brief Set states and time for a T-transform stage evaluation
956 * @par Runtime: O(nStates)
957 *
958 * In multirate mode, slow states are restored from the slow-state cache and
959 * fast states are overwritten from the input vector using the fast state index
960 * mapping. In single-rate mode, the full state vector is copied from the input.
961 *
962 * The simulation time is set to the corresponding stage time using the tableau
963 * coefficient. If stage < 0, the left interval boundary time is used.
964 *
965 * @param[in,out] data Data
966 * @param[in] gbData GBODE data
967 * @param[in] nls Internal nonlinear system data
968 * @param[in] y_fast Input state vector (size nls->size)
969 * @param[in] stage Stage index of the tableau, or negative for left boundary (use macro GB_INTERNAL_LEFT_BOUNDARY)
970 */
971 ✗ static inline void gbInternal_T_Transform_set_states(DATA *data, DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls, const double *y_fast, int stage)
972 {
973 ✗ if (nls->multirate)
974 {
975 ✗ DATA_GBODEF *gbfData = gbData->gbfData;
976
977 ✗ if (stage >= 0)
978 {
979 ✗ slowStateCache_overwrite_stage(gbData, gbfData->slowStateCache, stage, data->localData[0]->realVars);
980 }
981 else
982 {
983 ✗ slowStateCache_overwrite_left(gbData, gbfData->slowStateCache, data->localData[0]->realVars);
984 }
985
986 ✗ for (int fast_idx = 0; fast_idx < nls->size; fast_idx++)
987 {
988 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
989 ✗ data->localData[0]->realVars[full_idx] = y_fast[fast_idx];
990 }
991
992 ✗ data->localData[0]->timeValue = (stage >= 0 ? gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize : gbfData->time);
993 }
994 else
995 {
996 ✗ memcpy(data->localData[0]->realVars, y_fast, nls->size * sizeof(double));
997 ✗ data->localData[0]->timeValue = (stage >= 0 ? gbData->time + gbData->tableau->c[stage] * gbData->stepSize : gbData->time);
998 }
999 ✗ }
1000
1001 ✗ static inline void gbInternal_T_Transform_copy_full_to_fast(DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls, const double *src_full, double *dest_fast)
1002 {
1003 ✗ if (nls->multirate)
1004 {
1005 ✗ for (int fast_idx = 0; fast_idx < nls->size; fast_idx++)
1006 {
1007 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
1008 ✗ dest_fast[fast_idx] = src_full[full_idx];
1009 }
1010 }
1011 else
1012 {
1013 ✗ memcpy(dest_fast, src_full, nls->size * sizeof(double));
1014 }
1015 ✗ }
1016
1017 ✗ static inline void gbInternal_T_Transform_copy_fast_to_full(DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls, const double *src_fast, double *dest_full)
1018 {
1019 ✗ if (nls->multirate)
1020 {
1021 ✗ for (int fast_idx = 0; fast_idx < nls->size; fast_idx++)
1022 {
1023 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
1024 ✗ dest_full[full_idx] = src_fast[fast_idx];
1025 }
1026 }
1027 else
1028 {
1029 ✗ memcpy(dest_full, src_fast, nls->size * sizeof(double));
1030 }
1031 ✗ }
1032
1033 ✗ static inline void gbInternal_T_Transform_fast_to_full_axpy(DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls, const double alpha, const double *x_fast, double *y_full)
1034 {
1035 ✗ if (nls->multirate)
1036 {
1037 ✗ for (int fast_idx = 0; fast_idx < nls->size; fast_idx++)
1038 {
1039 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
1040 ✗ y_full[full_idx] += alpha * x_fast[fast_idx];
1041 }
1042 }
1043 else
1044 {
1045 ✗ daxpy_(&nls->size, &alpha, x_fast, &INT_ONE, y_full, &INT_ONE);
1046 }
1047 ✗ }
1048
1049 ✗ static inline void gbInternal_T_Transform_full_to_fast_axpy(DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls, const double alpha, const double *x_full, double *y_fast)
1050 {
1051 ✗ if (nls->multirate)
1052 {
1053 ✗ for (int fast_idx = 0; fast_idx < nls->size; fast_idx++)
1054 {
1055 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
1056 ✗ y_fast[fast_idx] += alpha * x_full[full_idx];
1057 }
1058 }
1059 else
1060 {
1061 ✗ daxpy_(&nls->size, &alpha, x_full, &INT_ONE, y_fast, &INT_ONE);
1062 }
1063 ✗ }
1064
1065 static inline void gbInternal_T_Transform_copy_full_to_full_axpy(DATA_GBODE *gbData, GB_INTERNAL_NLS_DATA *nls, const double alpha, const double *x_full, double *y_full)
1066 {
1067 if (nls->multirate)
1068 {
1069 for (int fast_idx = 0; fast_idx < nls->size; fast_idx++)
1070 {
1071 int full_idx = gbData->fastStatesIdx[fast_idx];
1072 y_full[full_idx] += alpha * x_full[full_idx];
1073 }
1074 }
1075 else
1076 {
1077 daxpy_(&nls->size, &alpha, x_full, &INT_ONE, y_full, &INT_ONE);
1078 }
1079 }
1080
1081 /** @brief Solve entire NLS of FIRK with possibly singular Runge-Kutta matrix
1082 * via the T-transformation (decoupled space).
1083 *
1084 * After convergence the solutions are written into nonlinsys->x and the stage updates are
1085 * written into gbData->k or gbfData->kCurrPacked.
1086 */
1087 ✗ static NLS_SOLVER_STATUS gbInternalSolveNls_T_Transform(DATA *data,
1088 threadData_t *threadData,
1089 NONLINEAR_SYSTEM_DATA* nonlinsys,
1090 DATA_GBODE* gbData,
1091 GB_INTERNAL_NLS_DATA *nls)
1092 {
1093 ✗ int size = nls->size;
1094 ✗ int w_size = nls->size * nls->tabl->t_transform->size;
1095 ✗ double stepSize = (nls->multirate ? gbData->gbfData->stepSize : gbData->stepSize);
1096 ✗ double lastStepSize = (nls->multirate ? gbData->gbfData->lastStepSize : gbData->lastStepSize);
1097 ✗ double invh = 1.0 / stepSize;
1098 ✗ double minvh = -invh;
1099 ✗ double *x = nonlinsys->nlsx;
1100 ✗ double *x_start = nonlinsys->nlsxOld;
1101 ✗ double *flat_res = nonlinsys->resValues;
1102
1103 ✗ double *yOld = (nls->multirate ? gbData->gbfData->yOldPacked : gbData->yOld);
1104 ✗ double *kPacked = (nls->multirate ? gbData->gbfData->kCurrPacked : gbData->k);
1105
1106 ✗ EVAL_SELECTION *selection = (nls->multirate ? gbData->gbfData->evalSelectionFast : NULL);
1107
1108 ✗ SOLVERSTATS *stats = (nls->multirate ? &gbData->gbfData->stats : &gbData->stats);
1109
1110 // return code from KLU or RHS / Jac calls
1111 int ret;
1112
1113 ✗ createGbScales(nls, gbData, x, x_start);
1114 double *scal = nls->scal;
1115
1116 ✗ SPARSE_PATTERN *ode_pattern = getODEPattern(data, gbData, nls);
1117 SPARSE_PATTERN *nls_pattern = getNLSPattern(gbData, nls);
1118 ✗ T_TRANSFORM *transform = nls->tabl->t_transform;
1119
1120 modelica_boolean jac_called = FALSE;
1121
1122 ✗ if (nls->call_jac || transform->firstRowZero || gbData->eventHappened)
1123 {
1124 /* set values for known last point (simplified Newton) */
1125 ✗ gbInternal_T_Transform_set_states(data, gbData, nls, yOld, GB_INTERNAL_LEFT_BOUNDARY);
1126
1127 /* callback ODE + callback Jacobian of ODE -> nls_jacobian buffer
1128 => TODO: if method has property a_{s,:} = b, i.e. FSAL or kLeft is available, we could recycle the k_s from before here
1129 or does fODE set algebraic variables that may be needed for the Jacobian? */
1130 ✗ ret = gbode_fODE(data, threadData, &stats->nCallsODE, selection);
1131 ✗ if (ret < 0) return NLS_FAILED;
1132
1133 ✗ if (transform->firstRowZero)
1134 {
1135 // save explicit stage (e.g. Lobatto IIIA)
1136 ✗ gbInternal_T_Transform_copy_full_to_fast(gbData, nls, &data->localData[0]->realVars[gbData->nStates], kPacked);
1137 }
1138
1139 ✗ if (nls->call_jac || gbData->eventHappened)
1140 {
1141 ✗ ret = gbInternal_evalJacobian(data, threadData, gbData, nls);
1142 ✗ if (ret < 0) return ret;
1143
1144 jac_called = TRUE;
1145 }
1146 }
1147
1148 ✗ if (jac_called || stepSize != lastStepSize)
1149 {
1150 ✗ for (int sys_real = 0; sys_real < transform->nRealEigenvalues; sys_real++)
1151 {
1152 /* create Jacobian real: gamma/h * I - J_f */
1153 ✗ jacobian_real_assemble(data, threadData, gbData, nls, transform->gamma[sys_real],
1154 ✗ ode_pattern, nls->jacobian_callback, nls->real_nls_jacs[sys_real]);
1155 ✗ ret = gbInternal_dKLU_factorize(&nls->klu, sys_real,
1156 ✗ (int *) nls_pattern->leadindex,
1157 ✗ (int *) nls_pattern->index,
1158 ✗ nls->real_nls_jacs[sys_real]);
1159 ✗ if (ret < 0) return NLS_FAILED;
1160 }
1161 ✗ for (int sys_cmplx = 0; sys_cmplx < transform->nComplexEigenpairs; sys_cmplx++)
1162 {
1163 /* create Jacobian complex: (alpha + i * beta)/h * I - J_f */
1164 ✗ jacobian_cmplx_assemble(data, threadData, gbData, nls, transform->alpha[sys_cmplx], transform->beta[sys_cmplx],
1165 ✗ ode_pattern, nls->jacobian_callback, nls->cmplx_nls_jacs[sys_cmplx]);
1166 ✗ ret = gbInternal_zKLU_factorize(&nls->klu, sys_cmplx,
1167 ✗ (int *) nls_pattern->leadindex,
1168 ✗ (int *) nls_pattern->index,
1169 ✗ nls->cmplx_nls_jacs[sys_cmplx]);
1170 ✗ if (ret < 0) return NLS_FAILED;
1171 }
1172 }
1173
1174 // we solve for Z = X(t_ij) - yOld or W = (T^{-1} otimes I) * Z, then get K back via K = 1/h * A^{-1} * Z
1175
1176 // set guess Z[j] = X_start[j] - yOld
1177 ✗ for (int j = 0; j < transform->size; j++)
1178 {
1179 ✗ memcpy(&nls->Z[j * size], &x_start[j * size], size * sizeof(double));
1180 ✗ daxpy_(&size, &DBL_MINUS_ONE, yOld, &INT_ONE, &nls->Z[j * size], &INT_ONE);
1181 }
1182
1183 // W = (T^{-1} otimes I) * Z
1184 ✗ dense_kron_id_vec(transform->size, size, transform->T_inv, nls->Z, nls->W);
1185
1186 // norms, convergence rate
1187 double nrm_x = 0;
1188 double nrm_delta = 0;
1189 double nrm_delta_prev = 0;
1190 double theta = 0;
1191
1192 // invalidate eta if an event happened
1193 ✗ if (gbData->eventHappened) *nls->etas = DBL_MAX;
1194
1195 // Newton iteration count - we start with newt_it = 1, because we need this for the step size selection and conditions below
1196 ✗ for (int newt_it = 1 ;; newt_it++)
1197 ✗ {
1198 // compute residuals: rhs := -1 / h (Lambda otimes I) * W + (T^{-1} * I) * F((T otimes I) * W) + Phi (Phi = T^{-1} A_part^{-1} * K_1 if (K_1 explicit else 0))
1199
1200 // work[j] = F((T otimes I) * W)[j]
1201 ✗ for (int j = 0; j < transform->size; j++)
1202 {
1203 ✗ gbInternal_T_Transform_set_states(data, gbData, nls, yOld, j + (int)transform->firstRowZero);
1204 ✗ gbInternal_T_Transform_fast_to_full_axpy(gbData, nls, DBL_ONE, &nls->Z[j * size], data->localData[0]->realVars);
1205 ✗ ret = gbode_fODE(data, threadData, &stats->nCallsODE, selection);
1206 ✗ if (ret < 0) return NLS_FAILED;
1207
1208 ✗ gbInternal_T_Transform_copy_full_to_fast(gbData, nls, &data->localData[0]->realVars[gbData->nStates], &nls->work[j * size]);
1209 }
1210
1211 // rhs[j] = (T^{-1} otimes I) * F((T otimes I) * W)
1212 ✗ dense_kron_id_vec(transform->size, size, transform->T_inv, nls->work, flat_res);
1213
1214 // rhs[j] += -1 / h * ((Lambda + L) otimes I) * W
1215 ✗ scaled_transform_matvec(transform, size, minvh, nls->W, flat_res);
1216
1217 // add Phi = T^{-1} * A_part^{-1} * a{r, 1} * K_1 if first stage is explicit, else skip (we computed nls->phi = T^{-1} * A_part^{-1} * a{r, 1})
1218 // where r are all rows that belong to A_part
1219 ✗ if (transform->firstRowZero)
1220 {
1221 ✗ for (int j = 0; j < transform->size; j++)
1222 {
1223 ✗ daxpy_(&size, &transform->phi[j], kPacked,
1224 ✗ &INT_ONE, &flat_res[j * size], &INT_ONE);
1225 }
1226 }
1227
1228 ✗ for (int real_row = 0; real_row < transform->nRealBlocks; real_row++)
1229 {
1230 ✗ double *res_row = &flat_res[real_row * size];
1231
1232 ✗ if (transform->hasL[real_row])
1233 {
1234 ✗ for (int col = 0; col < real_row; col++)
1235 {
1236 ✗ double a = -invh * transform->L[GBODE_L_INDEX(real_row, col)];
1237 ✗ if (a != 0.0) daxpy_(&size, &a, &flat_res[col * size], &INT_ONE, res_row, &INT_ONE);
1238 }
1239 }
1240
1241 ✗ int sys = transform->realEigenvalueIndex[real_row];
1242 ✗ ret = gbInternal_dKLU_solve(&nls->klu, sys, size, res_row);
1243 ✗ if (ret < 0) return NLS_FAILED;
1244 }
1245
1246 int cmplx_row = transform->nRealBlocks;
1247 ✗ for (int cmplx_block = 0; cmplx_block < transform->nComplexBlocks; cmplx_block++)
1248 {
1249 int row0 = cmplx_row;
1250 ✗ int row1 = cmplx_row + 1;
1251 ✗ double *res0 = &flat_res[row0 * size];
1252 ✗ double *res1 = &flat_res[row1 * size];
1253
1254 ✗ if (transform->hasL[row0])
1255 {
1256 ✗ for (int col = 0; col < cmplx_row; col++)
1257 {
1258 ✗ double a = -invh * transform->L[GBODE_L_INDEX(row0, col)];
1259 ✗ if (a != 0.0) daxpy_(&size, &a, &flat_res[col * size], &INT_ONE, res0, &INT_ONE);
1260 }
1261 }
1262
1263 ✗ if (transform->hasL[row1])
1264 {
1265 ✗ for (int col = 0; col < cmplx_row; col++)
1266 {
1267 ✗ double a = -invh * transform->L[GBODE_L_INDEX(row1, col)];
1268 ✗ if (a != 0.0) daxpy_(&size, &a, &flat_res[col * size], &INT_ONE, res1, &INT_ONE);
1269 }
1270 }
1271
1272 ✗ int sys = transform->complexEigenpairIndex[cmplx_block];
1273 ✗ dcopy_(&size, res0, &INT_ONE, &nls->cmplx_nls_res[sys][0], &INT_TWO); // .real
1274 ✗ dcopy_(&size, res1, &INT_ONE, &nls->cmplx_nls_res[sys][1], &INT_TWO); // .imag
1275
1276 ✗ ret = gbInternal_zKLU_solve(&nls->klu, sys, size, &nls->cmplx_nls_res[sys][0]);
1277 ✗ if (ret < 0) return NLS_FAILED;
1278
1279 ✗ dcopy_(&size, &nls->cmplx_nls_res[sys][0], &INT_TWO, res0, &INT_ONE); // r1
1280 ✗ dcopy_(&size, &nls->cmplx_nls_res[sys][1], &INT_TWO, res1, &INT_ONE); // r2
1281
1282 ✗ cmplx_row += 2;
1283 }
1284
1285 // Newton step (we must do W += dW)
1286 ✗ daxpy_(&w_size, &DBL_ONE, flat_res, &INT_ONE, nls->W, &INT_ONE);
1287
1288 // Z = (T otimes I) * W
1289 ✗ dense_kron_id_vec(transform->size, size, transform->T, nls->W, nls->Z);
1290
1291 ✗ nrm_delta_prev = fmax(DBL_EPSILON, nrm_delta);
1292 ✗ nrm_delta = gbScalesNorm(nls, flat_res, nls->tabl->t_transform->size);
1293
1294 // handle absorption effects
1295 ✗ nrm_x = gbScalesNormXPlusZ(nls, yOld, nls->Z, transform->size);
1296 ✗ modelica_boolean absorption = (nrm_delta <= DBL_ABSORPTION * nrm_x);
1297
1298 ✗ if (newt_it > 1)
1299 {
1300 ✗ theta = nrm_delta / nrm_delta_prev;
1301
1302 // Newton failed -> divergence
1303 ✗ if (theta >= nls->theta_divergence && !absorption)
1304 {
1305 ✗ nls->call_jac = TRUE;
1306 ✗ return NLS_FAILED;
1307 }
1308
1309 ✗ *nls->etas = theta / (1 - theta);
1310 }
1311 else
1312 {
1313 ✗ *nls->etas = pow(fmax(*nls->etas, DBL_EPSILON), nls->eta_inital_damping);
1314 }
1315
1316 ✗ if (!isfinite(*nls->etas) || !isfinite(nrm_delta))
1317 {
1318 // Inf or NaN detected
1319 // Either RHS or Jacobian or solution of the system contained a Inf or NaN
1320 return NLS_FAILED;
1321 }
1322
1323 // Newton converged
1324 ✗ if (*nls->etas * nrm_delta < nls->fnewt || absorption)
1325 {
1326 ✗ if (theta < nls->theta_keep)
1327 {
1328 ✗ nls->call_jac = FALSE;
1329 }
1330 else
1331 {
1332 ✗ nls->call_jac = TRUE;
1333 }
1334
1335 // set solution X[j] = X_0 + Z[j]
1336 ✗ int offset = transform->firstRowZero ? size : 0;
1337 ✗ for (int j = 0; j < transform->size; j++)
1338 {
1339 ✗ memcpy(&x[j * size + offset], &nls->Z[j * size], size * sizeof(double));
1340 ✗ daxpy_(&size, &DBL_ONE, yOld, &INT_ONE, &x[j * size + offset], &INT_ONE);
1341 }
1342
1343 // recompute weights K from Z via K = 1 / h * (A_part^{-1} otimes I) * Z + rho * k_1 if k_1 explicit else 0 (rho := -A_part^{-1} * A_{r, 1})
1344 // where r are all rows that belong to A_part
1345 ✗ dense_kron_id_vec(transform->size, size, transform->A_part_inv, nls->Z, &kPacked[offset]);
1346 ✗ dscal_(&w_size, &invh, &kPacked[offset], &INT_ONE);
1347
1348 ✗ if (transform->firstRowZero)
1349 {
1350 // add k[j] += rho[j] * k_1 or k[j] -= (-A_part^{-1} * A_{r, 1} * k_1)[j] * k_1
1351 ✗ for (int j = 0; j < transform->size; j++)
1352 {
1353 ✗ daxpy_(&size, &transform->rho[j], kPacked,
1354 ✗ &INT_ONE, &kPacked[offset + j * size], &INT_ONE);
1355 }
1356 }
1357
1358 // for explicit first stage: copy x0 into x (e.g. Lobatto IIIA)
1359 ✗ if (transform->firstRowZero)
1360 {
1361 ✗ memcpy(x, yOld, size * sizeof(double));
1362 }
1363
1364 // for explicit last stage: compute final K_s and X_s
1365 ✗ if (transform->lastColumnZero)
1366 {
1367 ✗ int s_minus1 = nls->tabl->nStages - 1;
1368
1369 // X_s = x0 + h * sum{j=1}^{s-1} A_{s, j} * k_j as A_{s, s} == 0!
1370 ✗ memcpy(&x[size * s_minus1], yOld, size * sizeof(double));
1371 ✗ dgemm_(&CHAR_NO_TRANS, &CHAR_TRANS,
1372 &size, &INT_ONE, &s_minus1,
1373 &stepSize,
1374 kPacked, &size,
1375 ✗ &nls->tabl->A[s_minus1 * nls->tabl->nStages], &INT_ONE,
1376 &DBL_ONE,
1377 ✗ &x[size * s_minus1], &size);
1378
1379 ✗ gbInternal_T_Transform_set_states(data, gbData, nls, &x[size * s_minus1], s_minus1);
1380 ✗ ret = gbode_fODE(data, threadData, &stats->nCallsODE, selection);
1381 ✗ if (ret < 0) return NLS_FAILED;
1382
1383 // k_s = f(t + h * c_s, x0 + h * sum{j=1}^{s-1} A_{s, j} * k_j)
1384 ✗ gbInternal_T_Transform_copy_full_to_fast(gbData, nls, &data->localData[0]->realVars[gbData->nStates], &kPacked[size * s_minus1]);
1385 }
1386
1387 ✗ return NLS_SOLVED;
1388 }
1389
1390 // Newton failed -> iteration limit exceeded or too slow convergence
1391 ✗ if (newt_it == nls->max_newton_it || (pow(theta, nls->max_newton_it - newt_it) / (1 - theta) * nrm_delta > nls->fnewt))
1392 {
1393 ✗ nls->call_jac = TRUE;
1394 ✗ return NLS_FAILED;
1395 }
1396 }
1397 }
1398
1399 /* Allocate the internal memory + do symbolic analysis. */
1400 ✗ void *gbInternalNlsAllocate(int size,
1401 NLS_USERDATA* userData,
1402 modelica_boolean attemptRetry,
1403 modelica_boolean isFast)
1404 {
1405 ✗ DATA_GBODE *gbData = isFast ? NULL : (DATA_GBODE *) userData->solverData;
1406 ✗ DATA_GBODEF *gbfData = isFast ? (DATA_GBODEF *) userData->solverData : NULL;
1407 ✗ BUTCHER_TABLEAU *tabl = isFast ? gbfData->tableau : gbData->tableau;
1408 ✗ T_TRANSFORM *transform = tabl->t_transform;
1409 ✗ JACOBIAN* jacobian_ODE = getSymbolicOdeJacobian(userData->data);
1410
1411 ✗ GB_INTERNAL_NLS_DATA *nls = (GB_INTERNAL_NLS_DATA *) malloc(sizeof(GB_INTERNAL_NLS_DATA));
1412 ✗ gbInternal_KLU_initialize(&nls->klu, transform ? transform->nRealEigenvalues : 1, transform ? transform->nComplexEigenpairs : 0);
1413
1414 // multirate stuff
1415 ✗ nls->multirate = isFast;
1416 ✗ nls->new_fast_states = FALSE;
1417
1418 // to have sufficiently large buffers, we need an overestimate for the system of the form struct(I + J) * v = b
1419 unsigned int nls_nnz_estimate = 0;
1420
1421 ✗ nls->nls_user_data = userData;
1422 ✗ nls->size = jacobian_ODE->sizeRows;
1423 ✗ nls->jacobian_callback = (double *) malloc(jacobian_ODE->sparsePattern->nnz * sizeof(double));
1424 ✗ nls->ode_to_nls = (int *) malloc(jacobian_ODE->sparsePattern->nnz * sizeof(int));
1425 ✗ nls->nls_diag_indices = (int *) malloc(jacobian_ODE->sizeRows * sizeof(int));
1426
1427 ✗ nls->tabl = tabl;
1428 ✗ nls->use_t_transform = (transform != NULL);
1429
1430 ✗ SPARSE_PATTERN *odePattern = isFast ? gbfData->sparsePattern_ODE : getJacobianCscPattern(jacobian_ODE);
1431 ✗ SPARSE_PATTERN *nlsPattern = isFast ? gbfData->sparsePattern_NLS : gbData->sparsePattern_NLS;
1432 ✗ assertStreamPrint(NULL, odePattern != NULL && nlsPattern != NULL, "GBODE internal NLS requires sparse patterns.");
1433 ✗ nls_nnz_estimate = nlsPattern->nnz;
1434
1435 ✗ if (!nls->multirate)
1436 {
1437 // create ODE Jac -> NLS Jacobian mapping
1438 ✗ gbodeMapSparsePattern(odePattern, nlsPattern, jacobian_ODE->sizeRows, nls->ode_to_nls, nls->nls_diag_indices);
1439 }
1440
1441 ✗ nls->scal = (double *) malloc(jacobian_ODE->sizeRows * sizeof(double));
1442 ✗ nls->etas = (double *) malloc(tabl->nStages * sizeof(double));
1443
1444 ✗ for (int i = 0; i < tabl->nStages; i++)
1445 {
1446 ✗ nls->etas[i] = DBL_MAX;
1447 }
1448
1449 ✗ nls->integrator_tol = userData->data->simulationInfo->tolerance;
1450
1451 /* Internal Newton convergence criteria is written in terms of the raw (unscaled), i.e.
1452 user-provided tolerance: || v ||_scal <= alpha, where || v ||_scal = sqrt(1/N sum_i=1^N (v_i / (ATOL_i + |y_i| * RTOL_i))^2),
1453 ATOL and RTOL are unscaled tolerances. */
1454 const double alpha_default = 3e-2;
1455 const double alpha_maximal = 5e-2;
1456 const double safety_newt = 0.1;
1457 double target_alpha = alpha_default;
1458
1459 ✗ if (!tabl->richardson && tabl->error_order < tabl->order_b && tabl->order_b - tabl->error_order != 1)
1460 {
1461 ✗ const double order_quot = ((double)tabl->error_order + 1.0) / ((double)tabl->order_b + 1.0);
1462 ✗ target_alpha = pow(safety_newt, 1.0 / order_quot);
1463 }
1464 ✗ nls->fnewt = fmax(DBL_ABSORPTION / nls->integrator_tol, fmin(alpha_maximal, target_alpha));
1465
1466 // damping power for eta
1467 ✗ if (omc_flag[FLAG_SR_NLS_INTERNAL_DAMPING_FAC])
1468 {
1469 ✗ double eta_damping = atof(omc_flagValue[FLAG_SR_NLS_INTERNAL_DAMPING_FAC]);
1470
1471 ✗ if (eta_damping > 1.0 || eta_damping < 0.0)
1472 {
1473 ✗ throwStreamPrint(NULL, "Invalid value %1.6e for flag '-gbnls_internal_damping'. Value must be less or equal to 1 and greater or equal to 0.", eta_damping);
1474 }
1475 else
1476 {
1477 ✗ nls->eta_inital_damping = eta_damping;
1478 }
1479 }
1480 else
1481 {
1482 ✗ nls->eta_inital_damping = 0.8;
1483 }
1484
1485 // add a history of thetas_last + #newt iterations to detect nearly linear systems, similar to err controller
1486 ✗ if (omc_flag[FLAG_SR_NLS_INTERNAL_JACKEEP])
1487 {
1488 ✗ double keep_flag_value = atof(omc_flagValue[FLAG_SR_NLS_INTERNAL_JACKEEP]);
1489 ✗ if (keep_flag_value >= 1.0)
1490 {
1491 ✗ throwStreamPrint(NULL, "Invalid value %1.6e for flag '-gbnls_internal_jackeep'. Value must be strictly less than 1.", keep_flag_value);
1492 }
1493 else
1494 {
1495 ✗ nls->theta_keep = keep_flag_value;
1496 }
1497 }
1498 else
1499 {
1500 // heuristic that takes sparsity into account
1501 ✗ if (nls->size > 8)
1502 {
1503 ✗ nls->theta_keep = pow(10.0, -3.0 + 1.75 * log(1.0 + (double)odePattern->maxColors) / log(1.0 + (double)nls->size));
1504 }
1505 else
1506 {
1507 // for very small systems we can compute the Jacobian frequently
1508 ✗ nls->theta_keep = 1e-3;
1509 }
1510 }
1511
1512 ✗ nls->call_jac = TRUE;
1513 ✗ nls->theta_divergence = 0.99;
1514 ✗ nls->max_newton_it = !transform ? 5 : 4 + 2 * transform->size; // = 5 for each (E)SDIRK stage and e.g. 10 for full RadauIIA 3-step
1515
1516 ✗ if (!transform)
1517 {
1518 ✗ nls->real_nls_jacs = (double **) malloc(sizeof(double *));
1519 ✗ nls->real_nls_jacs[0] = (double *) malloc(nls_nnz_estimate * sizeof(double));
1520
1521 // auxiliary memory
1522 ✗ nls->work = (double *) malloc(4 * nls->size * sizeof(double));
1523
1524 ✗ if (!nls->multirate)
1525 {
1526 ✗ gbInternal_KLU_analyze(&nls->klu, nls->size, (int *) nlsPattern->leadindex, (int *) nlsPattern->index);
1527 }
1528 }
1529 else
1530 {
1531 ✗ nls->real_nls_jacs = (double **) malloc(transform->nRealEigenvalues * sizeof(double *));
1532 ✗ nls->real_nls_res = (double **) malloc(transform->nRealEigenvalues * sizeof(double *));
1533 ✗ nls->cmplx_nls_jacs = (double **) malloc(transform->nComplexEigenpairs * sizeof(double *));
1534 ✗ nls->cmplx_nls_res = (double **) malloc(transform->nComplexEigenpairs * sizeof(double *));
1535
1536 ✗ for (int sys_real = 0; sys_real < transform->nRealEigenvalues; sys_real++)
1537 {
1538 ✗ nls->real_nls_res[sys_real] = (double *) malloc(nls->size * sizeof(double));
1539 ✗ nls->real_nls_jacs[sys_real] = (double *) malloc(nls_nnz_estimate * sizeof(double));
1540 }
1541 ✗ for (int sys_cmplx = 0; sys_cmplx < transform->nComplexEigenpairs; sys_cmplx++)
1542 {
1543 ✗ nls->cmplx_nls_res[sys_cmplx] = (double *) malloc(2 * nls->size * sizeof(double));
1544 ✗ nls->cmplx_nls_jacs[sys_cmplx] = (double *) malloc(2 * nls_nnz_estimate * sizeof(double));
1545 }
1546
1547 ✗ if (!nls->multirate)
1548 {
1549 ✗ gbInternal_KLU_analyze(&nls->klu, nls->size, (int *) nlsPattern->leadindex, (int *) nlsPattern->index);
1550 }
1551
1552 // iterate
1553 ✗ nls->Z = (double *) malloc(nls->size * transform->size * sizeof(double));
1554 ✗ nls->W = (double *) malloc(nls->size * transform->size * sizeof(double));
1555
1556 // auxiliary memory
1557 ✗ nls->work = (double *) malloc(nls->size * MAX(transform->size, 4) * sizeof(double));
1558 }
1559
1560 ✗ return (void *) nls;
1561 }
1562
1563 /* Free the internal memory. */
1564 ✗ void gbInternalNlsFree(void *nls_ptr)
1565 {
1566 GB_INTERNAL_NLS_DATA *nls = (GB_INTERNAL_NLS_DATA *) nls_ptr;
1567 ✗ gbInternal_KLU_free(&nls->klu);
1568 ✗ free(nls->jacobian_callback);
1569 ✗ free(nls->ode_to_nls);
1570 ✗ free(nls->nls_diag_indices);
1571 ✗ free(nls->scal);
1572 ✗ free(nls->etas);
1573 ✗ free(nls->work);
1574
1575 ✗ if (!nls->tabl->t_transform)
1576 {
1577 ✗ free(nls->real_nls_jacs[0]);
1578 ✗ free(nls->real_nls_jacs);
1579 }
1580 else
1581 {
1582 ✗ for (int sys_real = 0; sys_real < nls->tabl->t_transform->nRealEigenvalues; sys_real++)
1583 {
1584 ✗ free(nls->real_nls_jacs[sys_real]);
1585 ✗ free(nls->real_nls_res[sys_real]);
1586 }
1587 ✗ for (int sys_cmplx = 0; sys_cmplx < nls->tabl->t_transform->nComplexEigenpairs; sys_cmplx++)
1588 {
1589 ✗ free(nls->cmplx_nls_jacs[sys_cmplx]);
1590 ✗ free(nls->cmplx_nls_res[sys_cmplx]);
1591 }
1592
1593 ✗ free(nls->real_nls_jacs);
1594 ✗ free(nls->real_nls_res);
1595 ✗ free(nls->cmplx_nls_jacs);
1596 ✗ free(nls->cmplx_nls_res);
1597
1598 ✗ free(nls->Z);
1599 ✗ free(nls->W);
1600 }
1601
1602 ✗ freeNlsUserData(nls->nls_user_data);
1603 ✗ free(nls);
1604 ✗ }
1605
1606 ✗ void gbInternalScheduleFastStatesUpdate(void *nls_ptr)
1607 {
1608 ✗ assert(((GB_INTERNAL_NLS_DATA *) nls_ptr)->multirate);
1609 ✗ ((GB_INTERNAL_NLS_DATA *) nls_ptr)->new_fast_states = TRUE;
1610 ✗ }
1611
1612 ✗ modelica_boolean updateFastStates(DATA *data,
1613 threadData_t *threadData,
1614 NONLINEAR_SYSTEM_DATA* nonlinsys,
1615 DATA_GBODE* gbData,
1616 GB_INTERNAL_NLS_DATA *nls)
1617 {
1618 ✗ DATA_GBODEF *gbfData = gbData->gbfData;
1619 ✗ SPARSE_PATTERN *odePattern = gbfData->sparsePattern_ODE;
1620 ✗ SPARSE_PATTERN *nlsPattern = gbfData->sparsePattern_NLS;
1621
1622 // update size
1623 ✗ nls->size = gbData->nFastStates;
1624 ✗ nls->call_jac = TRUE;
1625
1626 ✗ for (int stage = 0; stage < gbfData->tableau->nStages; stage++)
1627 {
1628 ✗ nls->etas[stage] = DBL_MAX;
1629 }
1630
1631 // update mappings: ode_to_nls and nls_diag_indices
1632 ✗ gbodeMapSparsePattern(odePattern, nlsPattern, nls->size, nls->ode_to_nls, nls->nls_diag_indices);
1633
1634 // all transformed systems have the same sparsity pattern and share one symbolic analysis
1635 ✗ gbInternal_KLU_reanalyze(&nls->klu, nls->size, (int *) nlsPattern->leadindex, (int *) nlsPattern->index);
1636
1637 ✗ return TRUE;
1638 }
1639
1640 /* Entry point for `internal` solve routine: DIRK or FIRK */
1641 ✗ NLS_SOLVER_STATUS gbInternalSolveNls(DATA *data,
1642 threadData_t *threadData,
1643 NONLINEAR_SYSTEM_DATA* nonlinsys,
1644 DATA_GBODE* gbData,
1645 void *nls_ptr)
1646 {
1647 GB_INTERNAL_NLS_DATA *nls = (GB_INTERNAL_NLS_DATA *) nls_ptr;
1648
1649 ✗ if (nls->new_fast_states)
1650 {
1651 // update sparse pattern struct(I + J) and create a new symbolic factorization
1652 ✗ modelica_boolean success = updateFastStates(data, threadData, nonlinsys, gbData, nls);
1653 ✗ if (!success) return NLS_FAILED;
1654 ✗ nls->new_fast_states = FALSE;
1655 }
1656
1657 ✗ if (nls->use_t_transform)
1658 {
1659 ✗ return gbInternalSolveNls_T_Transform(data, threadData, nonlinsys, gbData, nls);
1660 }
1661 else
1662 {
1663 ✗ return gbInternalSolveNls_DIRK(data, threadData, nonlinsys, gbData, nls);
1664 }
1665 }
1666
1667 /**
1668 * @brief Contractive error estimate for stiff problems. (stiffness filter)
1669 *
1670 * Construct an embedded method of order `nStages` for a given collocation method
1671 * with at least one real eigenvalue and uncollocated point 0.0. We exclude
1672 * complex eigenvalues as this work would be even more expensive then.
1673 *
1674 * This estimate is A-stable and of one order higher than the naive embedded method, which
1675 * is crucial for stiff problems.
1676 *
1677 * See notes on struct CONTRACTIVE_DEFECT for more context.
1678 */
1679 ✗ void gbInternalContractiveDefect(DATA *data,
1680 threadData_t *threadData,
1681 NONLINEAR_SYSTEM_DATA *nonlinsys,
1682 DATA_GBODE *gbData,
1683 CONTRACTIVE_DEFECT *contractive,
1684 double *err)
1685 {
1686 ✗ GB_INTERNAL_NLS_DATA *nls = (GB_INTERNAL_NLS_DATA *) (((struct dataSolver *)nonlinsys->solverData)->ordinaryData);
1687 ✗ BUTCHER_TABLEAU *tabl = nls->tabl;
1688 ✗ SOLVERSTATS *stats = (nls->multirate ? &gbData->gbfData->stats : &gbData->stats);
1689
1690 ✗ int nStates = gbData->nStates;
1691 ✗ int size = nls->size;
1692 ✗ int nStages = tabl->nStages;
1693
1694 ✗ double *yOld = (nls->multirate ? gbData->gbfData->yOldPacked : gbData->yOld);
1695 ✗ double *kPacked = (nls->multirate ? gbData->gbfData->kCurrPacked : gbData->k);
1696
1697 // ERR := -d(0)^T * A * k
1698 ✗ dgemm_(&CHAR_NO_TRANS, &CHAR_NO_TRANS,
1699 &size,
1700 &INT_ONE,
1701 &nStages,
1702 &DBL_MINUS_ONE, kPacked, &size,
1703 ✗ contractive->dT_A, &nStages,
1704 &DBL_ZERO, err, &size);
1705
1706 ✗ modelica_boolean sr_valid = (!nls->multirate && !gbData->didFastStep && gbData->time != data->simulationInfo->startTime && !gbData->eventHappened && gbData->extrapolationBaseTime != INFINITY);
1707 ✗ modelica_boolean mr_valid = (nls->multirate && gbData->didFastStep && gbData->gbfData->extrapolationValid);
1708
1709 ✗ if (tabl->isKRightAvailable && (sr_valid || mr_valid))
1710 ✗ {
1711 // f(t_n, y(t_n)) == k_right == f(t_n-1 + h_n-1, y(t_n-1 + h_n-1)) of the previous step up to NLS precision
1712 // similar to the FSAL property of ESDIRK or ERK, but the method by itself does not have c_0 = 0 as a node!
1713 // Note that the order of this k_right is the order p of the method (e.g. p = 2s-1 for Radau IIA), since kRight is a collocated node
1714 // of the previous interval. Therefore, the computed defect will still be O(h^s) even though f(x0, y0) = k_right + O(h^p)
1715 ✗ double *k0Packed = (nls->multirate ? &gbData->gbfData->kLast[(tabl->nStages - 1) * gbData->nFastStates]
1716 ✗ : &gbData->kLast[(tabl->nStages - 1) * gbData->nStates]);
1717
1718 // ERR := f(t_n, y(t_n)) - d(0)^T * A * k
1719 ✗ daxpy_(&nls->size, &DBL_ONE, k0Packed, &INT_ONE, err, &INT_ONE);
1720 }
1721 else
1722 {
1723 // fresh computation of f(t_n, y(t_n)) as previous step is not valid or method does not collocate node 1
1724 ✗ gbInternal_T_Transform_set_states(data, gbData, nls, yOld, GB_INTERNAL_LEFT_BOUNDARY);
1725 ✗ gbode_fODE(data, threadData, &stats->nCallsODE, (nls->multirate ? gbData->gbfData->evalSelectionFast : NULL));
1726
1727 // ERR := f(t_n, y(t_n)) - d(0)^T * A * k
1728 ✗ gbInternal_T_Transform_full_to_fast_axpy(gbData, nls, DBL_ONE, &data->localData[0]->realVars[nStates], err);
1729 }
1730
1731 // ERR := (gamma / h * I - J)^{-1} * yt = (gamma / h * I - J)^{-1} * (f(t_n, y(t_n)) - d(0)^T * A * k) (exact error measure)
1732 ✗ gbInternal_dKLU_solve(&nls->klu, 0, size, err);
1733 ✗ }
1734
1735 ✗ static void gbInternalContractiveFilterPacked(GB_INTERNAL_NLS_DATA *nls,
1736 DATA_GBODE *gbData,
1737 double *err)
1738 {
1739 ✗ int size = nls->size;
1740 ✗ double filter_scale = 1.0;
1741
1742 ✗ if (nls->use_t_transform)
1743 {
1744 ✗ double stepSize = nls->multirate ? gbData->gbfData->stepSize : gbData->stepSize;
1745 ✗ filter_scale = nls->tabl->t_transform->gamma[0] / stepSize;
1746 }
1747
1748 // DIRK systems use h*gamma*J - I, so the solve already applies the filter up to sign.
1749 // FIRK/T systems use gamma/h*I - J = gamma/h * (I - h/gamma*J), so scale by gamma/h after the solve.
1750 ✗ gbInternal_dKLU_solve(&nls->klu, 0, size, err);
1751 ✗ if (filter_scale != 1.0) dscal_(&size, &filter_scale, err, &INT_ONE);
1752 ✗ }
1753
1754 ✗ void gbInternalContractiveFilterError(NONLINEAR_SYSTEM_DATA *nonlinsys,
1755 DATA_GBODE *gbData,
1756 double *err)
1757 {
1758 ✗ GB_INTERNAL_NLS_DATA *nls = (GB_INTERNAL_NLS_DATA *) (((struct dataSolver *)nonlinsys->solverData)->ordinaryData);
1759
1760 ✗ if (!nls->multirate)
1761 {
1762 ✗ gbInternalContractiveFilterPacked(nls, gbData, err);
1763 }
1764 else
1765 {
1766 ✗ double *work = nls->work;
1767
1768 // work := fast(err)
1769 ✗ gbInternal_T_Transform_copy_full_to_fast(gbData, nls, err, work);
1770 ✗ gbInternalContractiveFilterPacked(nls, gbData, work);
1771
1772 // err := full(work)
1773 ✗ gbInternal_T_Transform_copy_fast_to_full(gbData, nls, work, err);
1774 }
1775 ✗ }
1776
1777 // returns a work pointer of at least 32 * N_STATES bytes == 4 * N_STATES * sizeof(double)
1778 ✗ double *gbInternalGetWorkPointer(void *nls_ptr)
1779 {
1780 ✗ return ((GB_INTERNAL_NLS_DATA *) nls_ptr)->work;
1781 }
1782
1783 /**
1784 * @brief Perform intrastep stage-value-prediction (linear combination of known stage values k).
1785 *
1786 * Calculates y_predictor := y0 + h * A_predictor[1] * k[1] + A_predictor[2] * k[2] + ... + A_predictor[s-1] * k[s-1]),
1787 * where values of A_predictor are choosen such that order and stability properties are nice.
1788 *
1789 * @note This function is located in gbode_internal_nls.c as we get symbol conflicts with
1790 * SUNDIALS BLAS symbols that use long int as index type and are included transitively.
1791 */
1792 ✗ void gbInternalLinearCombinationSVP(STAGE_VALUE_PREDICTORS *svp,
1793 int active_stage,
1794 int nStates,
1795 double stepSize,
1796 const double *K,
1797 const double *y0,
1798 double *ypred)
1799 {
1800 ✗ memcpy(ypred, y0, nStates * sizeof(double));
1801
1802 ✗ dgemv_(
1803 &CHAR_NO_TRANS,
1804 &nStates,
1805 &active_stage,
1806 &stepSize, K, &nStates,
1807 ✗ &svp->A_predictor[active_stage * svp->nStages], &INT_ONE,
1808 &DBL_ONE, ypred, &INT_ONE
1809 );
1810 ✗ }
1811