Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.5% 2 / 0 / 371
Functions: 2.3% 1 / 0 / 44
Branches: 0.5% 1 / 0 / 190
OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.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 gbode_util.c
29 */
30 #include "gbode_util.h"
31
32 #define GBODE_EPSILON DBL_EPSILON
33
34
35 // LA functions
36 /**
37 * @brief Scalar multiplication and vector addition a = b + s*c for selected indices.
38 *
39 * Determines the scalar multiplication of an vector and adds the result
40 * to another vector only for selected indices.
41 *
42 * a = b + s*c
43 *
44 * @return a Output vector
45 * @param b Input vector
46 * @param c Input vector
47 * @param s Scalar value
48 * @param nIdx Length of index vector
49 * @param idx Index vector
50 */
51 ✗ void addSmultVec_gbf(double* a, double* b, double *c, double s, int nIdx, int* idx)
52 {
53 int i, ii;
54
55 ✗ for (ii=0; ii<nIdx; ii++) {
56 ✗ i = idx[ii];
57 ✗ a[i] = b[i] + s*c[i];
58 }
59 ✗ }
60
61 /**
62 * @brief Scalar multiplication and vector addition a = b + s*c.
63 *
64 * Determines the scalar multiplication of an vector and adds the result
65 * to another vector.
66 *
67 * a = b + s*c
68 *
69 * @return a Output vector
70 * @param b Input vector
71 * @param c Input vector
72 * @param s Scalar value
73 * @param n Length of the vectors
74 */
75 ✗ void addSmultVec_gb(double* a, double* b, double *c, double s, int n)
76 {
77 int i;
78
79 ✗ for (i=0; i<n; i++) {
80 ✗ a[i] = b[i] + s*c[i];
81 }
82 ✗ }
83
84 /*
85 * ============================================================================
86 * Interpolation functions
87 * ============================================================================
88 */
89
90 /**
91 * @brief Linear interpolation of specific vector components
92 *
93 * @param ta Time value at the left hand side
94 * @param fa Function values at the left hand side
95 * @param tb Time value at the right hand side
96 * @param fb Function values at the right hand side
97 * @param t Time value at the interpolated time point
98 * @param f Function values at the interpolated time point
99 * @param n Size of vector f or size of index vector if non-NULL.
100 * @param idx Index vector, can be NULL.
101 * Specifies which parts of f should be interpolated.
102 */
103 ✗ void linear_interpolation(double ta, double* fa, double tb, double* fb, double t, double* f, int n, int* idx, int nStates)
104 {
105 double lambda, h0, h1;
106 int i, ii;
107
108 // omit division by zero
109 ✗ if (fabs(tb-ta) <= GBODE_EPSILON) {
110 ✗ copyVector_gbf(f, fb, (idx == NULL) ? nStates : n, idx);
111 ✗ return;
112 }
113
114 ✗ lambda = (t-ta)/(tb-ta);
115 ✗ h0 = 1-lambda;
116 h1 = lambda;
117
118 ✗ if (idx == NULL) {
119 ✗ for (i=0; i<nStates; i++) {
120 ✗ f[i] = h0*fa[i] + h1*fb[i];
121 }
122 } else {
123 ✗ for (ii=0; ii<n; ii++) {
124 ✗ i = idx[ii];
125 ✗ f[i] = h0*fa[i] + h1*fb[i];
126 }
127 }
128 return;
129 }
130
131 /**
132 * @brief Hermite interpolation of specific vector components
133 *
134 * @param ta Time value at the left hand side
135 * @param fa Function values at the left hand side
136 * @param dfa Derivative function values at the left hand side
137 * @param tb Time value at the right hand side
138 * @param fb Function values at the right hand side
139 * @param dfb Derivative function values at the right hand side
140 * @param t Time value at the interpolated time point
141 * @param f Function values at the interpolated time point
142 * @param n Size of vector f or size of index vector if non-NULL.
143 * @param idx Index vector, can be NULL.
144 * Specifies which parts of f should be interpolated.
145 */
146 ✗ void hermite_interpolation(double ta, double* fa, double* dfa, double tb, double* fb, double* dfb, double t, double* f, int n, int* idx, int nStates)
147 {
148 double tt, h00, h01, h10, h11;
149 int i, ii;
150
151 // omit division by zero
152 ✗ if (fabs(tb-ta) <= GBODE_EPSILON) {
153 ✗ copyVector_gbf(f, fb, (idx == NULL) ? nStates : n, idx);
154 ✗ return;
155 }
156
157 ✗ tt = (t-ta)/(tb-ta);
158 ✗ h00 = (1+2*tt)*(1-tt)*(1-tt);
159 ✗ h10 = (tb-ta)*tt*(1-tt)*(1-tt);
160 ✗ h01 = (3-2*tt)*tt*tt;
161 ✗ h11 = (tb-ta)*(tt-1)*tt*tt;
162
163 ✗ if (idx == NULL) {
164 ✗ for (i=0; i<nStates; i++) {
165 ✗ f[i] = h00*fa[i]+h10*dfa[i]+h01*fb[i]+h11*dfb[i];
166 }
167 } else {
168 ✗ for (ii=0; ii<n; ii++) {
169 ✗ i = idx[ii];
170 ✗ f[i] = h00*fa[i]+h10*dfa[i]+h01*fb[i]+h11*dfb[i];
171 }
172 }
173
174 return;
175 }
176
177 /**
178 * @brief Hermite interpolation of specific vector components (only right derivative used)
179 *
180 * @param ta Time value at the left hand side
181 * @param fa Function values at the left hand side
182 * @param tb Time value at the right hand side
183 * @param fb Function values at the right hand side
184 * @param dfb Derivative function values at the right hand side
185 * @param t Time value at the interpolated time point
186 * @param f Function values at the interpolated time point
187 * @param n Size of vector f or size of index vector if non-NULL.
188 * @param idx Index vector, can be NULL.
189 * Specifies which parts of f should be interpolated.
190 */
191 ✗ void hermite_interpolation_b(double ta, double* fa, double tb, double* fb, double* dfb, double t, double* f, int n, int* idx, int nStates)
192 {
193 double tat,tbt,tbta, h00, h01, h11;
194 int i, ii;
195
196 // omit division by zero
197 ✗ if (fabs(tb-ta) <= GBODE_EPSILON) {
198 ✗ copyVector_gbf(f, fb, (idx == NULL) ? nStates : n, idx);
199 ✗ return;
200 }
201
202 ✗ tat = (ta-t);
203 ✗ tbt = (tb-t);
204 tbta = (tb-ta);
205 ✗ h00 = tbt*tbt/(tbta*tbta);
206 ✗ h01 = tat*(tat - tbt)/(tbta*tbta);
207 ✗ h11 = tat*tbt/tbta;
208
209 ✗ if (idx == NULL) {
210 ✗ for (i=0; i<nStates; i++) {
211 ✗ f[i] = h00*fa[i]+h01*fb[i]+h11*dfb[i];
212 }
213 } else {
214 ✗ for (ii=0; ii<n; ii++) {
215 ✗ i = idx[ii];
216 ✗ f[i] = h00*fa[i]+h01*fb[i]+h11*dfb[i];
217 }
218 }
219
220 return;
221 }
222
223 /**
224 * @brief Hermite interpolation of specific vector components (only left derivative used)
225 *
226 * @param ta Time value at the left hand side
227 * @param fa Function values at the left hand side
228 * @param dfa Derivative function values at the left hand side
229 * @param tb Time value at the right hand side
230 * @param fb Function values at the right hand side
231 * @param t Time value at the interpolated time point
232 * @param f Function values at the interpolated time point
233 * @param n Size of vector f or size of index vector if non-NULL.
234 * @param idx Index vector, can be NULL.
235 * Specifies which parts of f should be interpolated.
236 */
237 ✗ void hermite_interpolation_a(double ta, double* fa, double* dfa, double tb, double* fb, double t, double* f, int n, int* idx, int nStates)
238 {
239 double tat,tbt,tbta, h00, h01, h10;
240 int i, ii;
241
242 // omit division by zero
243 ✗ if (fabs(tb-ta) <= GBODE_EPSILON) {
244 ✗ copyVector_gbf(f, fb, (idx == NULL) ? nStates : n, idx);
245 ✗ return;
246 }
247
248 ✗ tat = (ta-t);
249 ✗ tbt = (tb-t);
250 tbta = (tb-ta);
251 ✗ h01 = tat*tat/(tbta*tbta);
252 ✗ h00 = 1 - h01;
253 ✗ h10 = -tat*tbt/tbta;
254
255 ✗ if (idx == NULL) {
256 ✗ for (i=0; i<nStates; i++) {
257 ✗ f[i] = h00*fa[i]+h01*fb[i]+h10*dfa[i];
258 }
259 } else {
260 ✗ for (ii=0; ii<n; ii++) {
261 ✗ i = idx[ii];
262 ✗ f[i] = h00*fa[i]+h01*fb[i]+h10*dfa[i];
263 }
264 }
265
266 return;
267 }
268
269 /**
270 * @brief Hermite interpolation of specific vector components
271 *
272 * @param interpolMethod
273 * @param ta Time value at the left hand side
274 * @param fa Function values at the left hand side
275 * @param dfa Derivative function values at the left hand side.
276 * Can be NULL for linear interpolation.
277 * @param tb Time value at the right hand side
278 * @param fb Function values at the right hand side
279 * @param dfb Derivative function values at the right hand side.
280 * Can be NULL for linear interpolation.
281 * @param t Time value at the interpolated time point
282 * @param f Function values at the interpolated time point
283 * @param n Size of vector f or size of index vector if non-NULL.
284 * @param idx Index vector, can be NULL.
285 * Specifies which parts of f should be interpolated.
286 */
287 ✗ void gb_interpolation(enum GB_INTERPOL_METHOD interpolMethod, double ta, double* fa, double* dfa, double tb, double* fb, double* dfb, double t, double* f,
288 int nIdx, int* idx, int nStates, BUTCHER_TABLEAU* tableau, double* x, double *k)
289 {
290 // handle edge case for tiny interval length
291 ✗ if ((tb == ta) || fabs(tb - ta) < GBODE_EPSILON * (fabs(tb) + fabs(ta)))
292 {
293 ✗ if (idx == NULL)
294 {
295 ✗ memcpy(f, fa, nStates * sizeof(double));
296 }
297 else
298 {
299 ✗ for (int ii = 0; ii < nStates; ii++)
300 {
301 ✗ int i = idx[ii];
302 ✗ f[i] = fa[i];
303 }
304 }
305
306 ✗ return;
307 }
308
309 ✗ switch (interpolMethod)
310 {
311 ✗ case GB_INTERPOL_LIN:
312 ✗ linear_interpolation(ta, fa, tb, fb, t, f, nIdx, idx, nStates);
313 ✗ break;
314 ✗ case GB_DENSE_OUTPUT:
315 case GB_DENSE_OUTPUT_ERRCTRL:
316 ✗ if (tableau->withDenseOutput) {
317 ✗ tableau->dense_output(tableau, fa, x, k, (t - ta)/(tb - ta), (tb - ta), f, nIdx, idx, nStates);
318 ✗ break;
319 }
320 case GB_INTERPOL_HERMITE_a:
321 ✗ hermite_interpolation_a(ta, fa, dfa, tb, fb, t, f, nIdx, idx, nStates);
322 ✗ break;
323 ✗ case GB_INTERPOL_HERMITE_b:
324 ✗ hermite_interpolation_b(ta, fa, tb, fb, dfb, t, f, nIdx, idx, nStates);
325 ✗ break;
326 ✗ case GB_INTERPOL_HERMITE_ERRCTRL:
327 case GB_INTERPOL_HERMITE:
328 ✗ hermite_interpolation(ta, fa, dfa, tb, fb, dfb, t, f, nIdx, idx, nStates);
329 ✗ break;
330 ✗ default:
331 ✗ throwStreamPrint(NULL, "Not handled case in gb_interpolation. Unknown interpolation method %i.", interpolMethod);
332 }
333 }
334
335 /**
336 * @brief Difference between linear and hermite interpolation at intermediate points.
337 *
338 * @param gbData
339 */
340 ✗ double error_interpolation_gb(DATA_GBODE* gbData, int nIdx, int* idx, double tol)
341 {
342 int i, ii;
343 double errint = 0.0, errtol;
344
345 ✗ if (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL || gbData->interpolation == GB_DENSE_OUTPUT) {
346 ✗ gb_interpolation(gbData->interpolation, gbData->timeLeft, gbData->yLeft, gbData->kLeft,
347 gbData->timeRight, gbData->yRight, gbData->kRight,
348 ✗ (gbData->timeLeft + gbData->timeRight)/2, gbData->y1,
349 nIdx, idx, gbData->nStates, gbData->tableau, gbData->x, gbData->k);
350 } else {
351 ✗ hermite_interpolation_a(gbData->timeLeft, gbData->yLeft, gbData->kLeft,
352 gbData->timeRight, gbData->yRight,
353 ✗ (gbData->timeLeft + gbData->timeRight)/2, gbData->y1,
354 nIdx, idx, gbData->nStates);
355 }
356 ✗ hermite_interpolation(gbData->timeLeft, gbData->yLeft, gbData->kLeft,
357 gbData->timeRight, gbData->yRight, gbData->kRight,
358 ✗ (gbData->timeLeft + gbData->timeRight)/2, gbData->y2,
359 nIdx, idx, gbData->nStates);
360 ✗ if (idx == NULL) {
361 ✗ for (i=0; i<nIdx; i++) {
362 ✗ errtol = tol * fmax(fabs(gbData->yLeft[i]), fabs(gbData->yRight[i])) + tol;
363 ✗ gbData->errest[i] = fabs(gbData->y2[i] - gbData->y1[i]) / errtol;
364 ✗ errint = fmax(errint, gbData->errest[i]);
365 }
366 } else {
367 ✗ for (ii=0; ii<nIdx; ii++) {
368 ✗ i = idx[ii];
369 ✗ errtol = tol * fmax(fabs(gbData->yLeft[i]), fabs(gbData->yRight[i])) + tol;
370 ✗ gbData->errest[i] = fabs(gbData->y2[i] - gbData->y1[i]) / errtol;
371 ✗ errint = fmax(errint, gbData->errest[i]);
372 }
373 }
374 ✗ return errint;
375 }
376
377 /**
378 * @brief Extrapolation for fast states.
379 *
380 * Using interpolation method specified in gbData->interpolation.
381 *
382 * @param gbData Generic ODE solver data.
383 * @param nlsxExtrapolation On output contains function values at extrapolation point time.
384 * @param time Extrapolation time.
385 */
386 ✗ void extrapolation_gbf(DATA_GBODE* gbData, double* nlsxExtrapolation, double time)
387 {
388 ✗ DATA_GBODEF* gbfData = gbData->gbfData;
389 ✗ const int nStates = gbData->nStates;
390 ✗ const int nFastStates = gbData->nFastStates;
391
392 ✗ if (fabs(gbfData->tv[1]-gbfData->tv[0]) <= GBODE_EPSILON) {
393 ✗ addSmultVec_gbf(nlsxExtrapolation, gbfData->yv, gbfData->kv, time - gbfData->tv[0], nFastStates, gbData->fastStatesIdx);
394 } else {
395 // this is actually extrapolation...
396 ✗ gb_interpolation(GB_INTERPOL_HERMITE,
397 ✗ gbfData->tv[1], gbfData->yv + nStates, gbfData->kv + nStates,
398 gbfData->tv[0], gbfData->yv, gbfData->kv,
399 time, nlsxExtrapolation,
400 nFastStates, gbData->fastStatesIdx, nStates, gbfData->tableau, gbfData->x, gbfData->k);
401 }
402 ✗ }
403
404 /**
405 * @brief Extrapolation for all states.
406 *
407 * Using interpolation method specified in gbData->interpolation.
408 *
409 * @param gbData Generic ODE solver data.
410 * @param nlsxExtrapolation On output contains function values at extrapolation point time.
411 * @param time Extrapolation time.
412 */
413 ✗ void extrapolation_hermite_gb(double* nlsxExtrapolation, int nStates, double t0, double *x0, double* k0, double t1, double *x1, double* k1, double time)
414 {
415 ✗ gb_interpolation(GB_INTERPOL_HERMITE,
416 t0, x0, k0,
417 t1, x1, k1,
418 time, nlsxExtrapolation,
419 nStates, NULL, nStates, NULL, NULL, NULL);
420 ✗ }
421
422 /**
423 * @brief Extrapolation for all states.
424 *
425 * Using interpolation method specified in gbData->interpolation.
426 *
427 * @param gbData Generic ODE solver data.
428 * @param nlsxExtrapolation On output contains function values at extrapolation point time.
429 * @param time Extrapolation time.
430 */
431 ✗ void extrapolation_gb(DATA_GBODE* gbData, double* nlsxExtrapolation, double time)
432 {
433 ✗ int nStates = gbData->nStates;
434
435 ✗ if (fabs(gbData->tv[1]-gbData->tv[0]) <= GBODE_EPSILON || gbData->multi_rate) {
436 ✗ addSmultVec_gb(nlsxExtrapolation, gbData->yv, gbData->kv, time - gbData->tv[0], nStates);
437 } else {
438 // this is actually extrapolation...
439 ✗ gb_interpolation(GB_INTERPOL_HERMITE,
440 ✗ gbData->tv[1], gbData->yv + nStates, gbData->kv + nStates,
441 gbData->tv[0], gbData->yv, gbData->kv,
442 time, nlsxExtrapolation,
443 nStates, NULL, nStates, gbData->tableau, gbData->x, gbData->k);
444 }
445 ✗ }
446
447 /**
448 * @brief Copy specific vector components given by an index vector
449 *
450 * if indx == NULL, the full vector is copied
451 *
452 * @param a Target vector
453 * @param b Source vector
454 * @param nIndx Size of the index vector
455 * @param indx Index vector
456 */
457 ✗ void copyVector_gbf(double* dest, double* src, int nIndx, int* indx)
458 {
459 ✗ if (indx != NULL) {
460 ✗ for (int i = 0; i < nIndx; i++)
461 ✗ dest[indx[i]] = src[indx[i]];
462 } else {
463 ✗ memcpy(dest, src, nIndx*sizeof(double));
464 }
465 ✗ }
466
467 /**
468 * @brief Projection function
469 *
470 * Collects the values in the vector for given indices (idx)
471 * and copy them in an corresponding vector of size (nIdx).
472 *
473 * @return a Target vector
474 * @param b Source Vector
475 * @param nIndx Length of index vector
476 * @param indx Index vector
477 */
478 ✗ void projVector_gbf(double* a, double* b, int nIndx, int* indx)
479 {
480 ✗ for (int i = 0; i < nIndx; i++)
481 ✗ a[i] = b[indx[i]];
482 ✗ }
483
484 /**
485 * @brief Output debug information of the states and derivatives
486 *
487 * that have been evaluated at the past accepted time points.
488 *
489 * @param stream Prints only, if stream is active
490 * @param x States at the past accepted time points
491 * @param k Derivatives at the past accepted time points
492 * @param t Past accepted time points
493 * @param nStates Number of states
494 * @param size Size of buffer
495 */
496 ✗ void debugRingBufferSteps_gb(enum OMC_LOG_STREAM stream, double* x, double* k, double *t, int nStates, int size)
497 {
498 // If stream is not active do nothing
499 ✗ if (!OMC_ACTIVE_STREAM(stream)) return;
500
501 ✗ infoStreamPrint(stream, 1, "States and derivatives at past accepted time steps:");
502
503 int i;
504
505 ✗ infoStreamPrint(stream, 0, "states:");
506 ✗ for (i = 0; i < size; i++) {
507 ✗ printVector_gb(stream, "x", x + i * nStates, nStates, t[i]);
508 }
509 ✗ infoStreamPrint(stream, 0, "derivatives:");
510 ✗ for (i = 0; i < size; i++) {
511 ✗ printVector_gb(stream, "k", k + i * nStates, nStates, t[i]);
512 }
513 ✗ messageClose(stream);
514 }
515
516 /**
517 * @brief Output debug information of the states and derivatives
518 *
519 * that have been evaluated at the past accepted time points.
520 *
521 * @param stream Prints only, if stream is active
522 * @param x States at the past accepted time points
523 * @param k Derivatives at the past accepted time points
524 * @param t Past accepted time points
525 * @param nStates Number of states
526 * @param size Size of buffer
527 * @param nIndx Size of index vector
528 * @param indx Index vector
529 */
530 ✗ void debugRingBufferSteps_gbf(enum OMC_LOG_STREAM stream, double* x, double* k, double *t, int nStates, int size, int nIndx, int* indx)
531 {
532 // If stream is not active do nothing
533 ✗ if (!OMC_ACTIVE_STREAM(stream)) return;
534
535 ✗ infoStreamPrint(stream, 1, "States and derivatives at past accepted time steps (inner integration):");
536
537 int i;
538
539 ✗ infoStreamPrint(stream, 0, "states:");
540 ✗ for (i = 0; i < size; i++) {
541 ✗ printVector_gbf(stream, "x", x + i * nStates, nStates, t[i], nIndx, indx);
542 }
543 ✗ infoStreamPrint(stream, 0, "derivatives:");
544 ✗ for (i = 0; i < size; i++) {
545 ✗ printVector_gbf(stream, "k", k + i * nStates, nStates, t[i], nIndx, indx);
546 }
547 ✗ messageClose(stream);
548 }
549
550 /**
551 * @brief Output debug information of the states and derivatives
552 *
553 * that have been evaluated at the intermediate points given by the
554 * Butcher tableau.
555 *
556 * @param stream Prints only, if stream is active
557 * @param x States at the intermediate time points
558 * @param k Derivatives at the intermediate time points
559 * @param nStates Number of states
560 * @param tableau Tableau of the Runge Kutta method
561 * @param time Current time of the inegrator (left hand side)
562 * @param stepSize Current step size of the integrator
563 */
564 ✗ void debugRingBuffer_gb(enum OMC_LOG_STREAM stream, double* x, double* k, int nStates, BUTCHER_TABLEAU* tableau, double time, double stepSize)
565 {
566 // If stream is not active do nothing
567 ✗ if (!OMC_ACTIVE_STREAM(stream)) return;
568
569 ✗ int nStages = tableau->nStages, stage_;
570
571 ✗ infoStreamPrint(stream, 0, "states:");
572 ✗ for (int stage_ = 0; stage_ < nStages; stage_++) {
573 ✗ printVector_gb(stream, "x", x + stage_ * nStates, nStates, time + tableau->c[stage_] * stepSize);
574 }
575 ✗ infoStreamPrint(stream, 0, "derivatives:");
576 ✗ for (int stage_ = 0; stage_ < nStages; stage_++) {
577 ✗ printVector_gb(stream, "k", k + stage_ * nStates, nStates, time + tableau->c[stage_] * stepSize);
578 }
579 }
580
581 /**
582 * @brief Output debug information of the states and derivatives
583 *
584 * that have been evaluated at the intermediate points given by the
585 * Butcher tableau.
586 *
587 * @param stream Prints only, if stream is active
588 * @param x States at the intermediate time points
589 * @param k Derivatives at the intermediate time points
590 * @param nStates Number of states
591 * @param tableau Tableau of the Runge Kutta method
592 * @param time Current time of the inegrator (left hand side)
593 * @param stepSize Current step size of the integrator
594 * @param nIndx Size of index vector
595 * @param indx Index vector
596 */
597 ✗ void debugRingBuffer_gbf(enum OMC_LOG_STREAM stream, double* x, double* k, int nStates, BUTCHER_TABLEAU* tableau, double time, double stepSize, int nIndx, int* indx)
598 {
599 // If stream is not active do nothing
600 ✗ if (!OMC_ACTIVE_STREAM(stream)) return;
601
602 ✗ int nStages = tableau->nStages, stage_;
603
604 ✗ infoStreamPrint(stream, 0, "states:");
605 ✗ for (int stage_ = 0; stage_ < nStages; stage_++) {
606 ✗ printVector_gbf(stream, "x", x + stage_ * nStates, nStates, time + tableau->c[stage_] * stepSize, nIndx, indx);
607 }
608 ✗ infoStreamPrint(stream, 0, "derivatives:");
609 ✗ for (int stage_ = 0; stage_ < nStages; stage_++) {
610 ✗ printVector_gbf(stream, "k", k + stage_ * nStates, nStates, time + tableau->c[stage_] * stepSize, nIndx, indx);
611 }
612 }
613
614 /**
615 * @brief Prints a vector to stream.
616 *
617 * If vector is larger than 1000 nothing is printed.
618 *
619 * @param stream Prints only, if stream is active
620 * @param name Specific string to print (usually name of the vector)
621 * @param a Vector to print
622 * @param n Size of the vector
623 * @param time Time value
624 */
625 ✗ void printVector_gb(enum OMC_LOG_STREAM stream, char name[], double* a, int n, double time)
626 {
627 // If stream is not active or size of vector to big do nothing
628 ✗ if (!OMC_ACTIVE_STREAM(stream) || n>1000) return;
629
630 // This only works for number of states less than 10!
631 // For large arrays, this is not a good output format!
632 enum { bufSize = 40960 };
633 char row_to_print[bufSize];
634 unsigned int ct;
635 ✗ ct = snprintf(row_to_print, bufSize, "%s(%8g) =\t", name, time);
636 ✗ for (int i=0;i<n;i++)
637 ✗ ct += snprintf(row_to_print+ct, bufSize-ct, " %16.12g", a[i]);
638 ✗ infoStreamPrint(stream, 0, "%s", row_to_print);
639 }
640
641 /**
642 * @brief Prints an integer vector to stream.
643 *
644 * If vector is larger than 1000 nothing is printed.
645 *
646 * @param name Specific string to print (usually name of the vector)
647 * @param a Integer vector to print
648 * @param n Size of the vector
649 * @param time Time value
650 */
651 ✗ void printIntVector_gb(enum OMC_LOG_STREAM stream, char name[], int* a, int n, double time)
652 {
653 // If stream is not active or size of vector to big do nothing
654 ✗ if (!OMC_ACTIVE_STREAM(stream) || n>1000) return;
655
656 enum { bufSize = 40960 };
657 char row_to_print[bufSize];
658 unsigned int ct;
659 ✗ ct = snprintf(row_to_print, bufSize, "%s(%8g) =\t", name, time);
660 ✗ for (int i=0;i<n;i++)
661 ✗ ct += snprintf(row_to_print+ct, bufSize-ct, " %d", a[i]);
662 ✗ infoStreamPrint(stream, 0, "%s", row_to_print);
663 }
664
665 /**
666 * @brief Prints selected vector components given by an index vector.
667 *
668 * If more than 1000 elements should be printed do nothing.
669 *
670 * @param name Specific string to print (usually name of the vector)
671 * @param a Vector to print
672 * @param n Size of the vector
673 * @param time Time value
674 * @param nIndx Size of index vector
675 * @param indx Index vector
676 */
677 ✗ void printVector_gbf(enum OMC_LOG_STREAM stream, char name[], double* a, int n, double time, int nIndx, int* indx)
678 {
679 // If stream is not active or size of vector to big do nothing
680 ✗ if (!OMC_ACTIVE_STREAM(stream) || nIndx>1000) return;
681
682 // This only works for number of states less than 10!
683 // For large arrays, this is not a good output format!
684 enum { bufSize = 40960 };
685 char row_to_print[bufSize];
686 unsigned int ct;
687 ✗ ct = snprintf(row_to_print, bufSize, "%s(%8g) =\t", name, time);
688 ✗ for (int i=0;i<nIndx;i++)
689 ✗ ct += snprintf(row_to_print+ct, bufSize-ct, " %16.12g", a[indx[i]]);
690 ✗ infoStreamPrint(stream, 0, "%s", row_to_print);
691 }
692
693 /**
694 * @brief Prints sparse structure.
695 *
696 * Use to print e.g. sparse Jacobian matrix.
697 * Only prints if stream is active and sparse pattern is non NULL and of size > 0.
698 *
699 * @param sparsePattern Matrix to print.
700 * @param sizeRows Number of rows of matrix.
701 * @param sizeCols Number of columns of matrix.
702 * @param stream Steam to print to.
703 * @param name Name of matrix.
704 */
705 ✗ void printSparseJacobianLocal(JACOBIAN* jacobian, const char* name)
706 {
707 /* Variables */
708 unsigned int row, col, i;
709 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Sparse structure of %s [size: %zux%zu]", name, jacobian->sizeRows, jacobian->sizeCols);
710 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "%u non-zero elements", jacobian->sparsePattern->nnz);
711 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Values of the transposed matrix (rows: states)");
712
713 printf("\n");
714 i=0;
715 ✗ for (row = 0; row < jacobian->sizeRows; row++) {
716 ✗ for (col = 0; col < jacobian->sizeRows; col++) {
717 ✗ if(jacobian->sparsePattern->index[i] == col) {
718 ✗ printf("%20.16g ", jacobian->resultVars[col]);
719 ✗ ++i;
720 } else {
721 printf("%20.16g ", 0.0);
722 }
723 }
724 printf("\n");
725 }
726 printf("\n");
727 ✗ }
728
729 /**
730 * @brief Write information on the active fast states on file (activity diagram)
731 *
732 * @param gbData Pointer to generic GBODE data struct.
733 * @param event If an event has happened, write zeros else ones
734 * @param time Actual time of reporting
735 * @param rejectedType Type of rejection
736 * 0 <= no rejection
737 * 1 <= error of slow states greater than the tolerance
738 * 2 <= interpolation error is too large
739 * 3 <= rejected because solving the NLS failed
740 * -1 <= step is preliminary accepted but needs refinement
741 */
742 ✗ void dumpFastStates_gb(DATA_GBODE* gbData, modelica_boolean event, double time, int rejectedType)
743 {
744 enum { bufSize = 4096 };
745 char fastStates_row[bufSize];
746 unsigned int ct;
747 ✗ ct = snprintf(fastStates_row, bufSize, "%15.10g %2d %15.10g %15.10g %15.10g", time, rejectedType, gbData->err_slow, gbData->err_int, gbData->err_fast);
748 ✗ for (int i = 0; i < gbData->nStates; i++) {
749 ✗ if (event)
750 ✗ ct += snprintf(fastStates_row+ct, bufSize-ct, " 0");
751 else
752 ✗ ct += snprintf(fastStates_row+ct, bufSize-ct, " 1");
753 }
754 ✗ fprintf(gbData->gbfData->fastStatesDebugFile, "%s\n", fastStates_row);
755 ✗ }
756
757 /**
758 * @brief Write information on the active fast states on file (activity diagram)
759 *
760 * @param gbData Pointer to generic GBODE data struct.
761 * @param time Actual time of reporting
762 * @param rejectedType Type of rejection
763 * 0 <= no rejection
764 * 1 <= error of fast states greater than the tolerance
765 * 2 <= interpolation error is too large
766 * 3 <= rejected because solving the NLS failed
767 * -1 <= step is preliminary accepted but needs refinement
768 */
769 ✗ void dumpFastStates_gbf(DATA_GBODE* gbData, double time, int rejectedType)
770 {
771 enum { bufSize = 40960 };
772 char fastStates_row[bufSize];
773 unsigned int ct;
774 int i, ii;
775 ✗ ct = snprintf(fastStates_row, bufSize, "%15.10g %2d %15.10g %15.10g %15.10g", time, rejectedType, gbData->err_slow, gbData->err_int, gbData->err_fast);
776 ✗ for (i = 0, ii = 0; i < gbData->nStates;) {
777 ✗ if (i == gbData->fastStatesIdx[ii]) {
778 ✗ ct += snprintf(fastStates_row+ct, bufSize-ct, " 1");
779 ✗ i++;
780 ✗ if (ii < gbData->nFastStates-1) ii++;
781 } else {
782 ✗ ct += snprintf(fastStates_row+ct, bufSize-ct, " 0");
783 ✗ i++;
784 }
785 }
786 ✗ fprintf(gbData->gbfData->fastStatesDebugFile, "%s\n", fastStates_row);
787 ✗ }
788
789 /**
790 * @brief Check if indices of fast states changed and update indices.
791 *
792 * @param gbData Pointer to gbode data.
793 * @return modelica_boolean TRUE if at least one fast state changed, FALSE otherwise.
794 */
795 ✗ modelica_boolean checkFastStatesChange(DATA_GBODE* gbData)
796 {
797 ✗ DATA_GBODEF* gbfData = gbData->gbfData;
798 modelica_boolean fastStatesChange = FALSE;
799
800 ✗ gbfData->nFastStates = gbData->nFastStates;
801 ✗ gbfData->fastStatesIdx = gbData->fastStatesIdx;
802
803 // check if number of fast states changed
804 ✗ if (gbfData->nFastStates_old != gbData->nFastStates) {
805 fastStatesChange = TRUE;
806 } else {
807 // look for changes in the ordering
808 // TODO memcmp() faster?
809 ✗ for (int k = 0; k < gbData->nFastStates; k++) {
810 ✗ if (gbfData->fastStates_old[k] != gbData->fastStatesIdx[k]) {
811 fastStatesChange = TRUE;
812 break;
813 }
814 }
815 }
816
817 ✗ if (fastStatesChange) {
818 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) {
819 ✗ printIntVector_gb(OMC_LOG_SOLVER, "old fast States:", gbfData->fastStates_old, gbfData->nFastStates_old, gbfData->time);
820 ✗ printIntVector_gb(OMC_LOG_SOLVER, "new fast States:", gbData->fastStatesIdx, gbData->nFastStates, gbfData->time);
821 }
822
823 // Update indices for the current fast states and corresponding counting
824 ✗ gbfData->nFastStates_old = gbData->nFastStates;
825 // TODO memcpy() faster?
826 ✗ for (int k = 0; k < gbData->nFastStates; k++) {
827 ✗ gbfData->fastStates_old[k] = gbData->fastStatesIdx[k];
828 }
829 }
830 ✗ return fastStatesChange;
831 }
832
833 /**
834 * @brief Log ODE integrator solver stats.
835 *
836 * @param name Name of ODE integrator.
837 * @param timeValue Current time value.
838 * @param integratorTime Time value of integrator.
839 * @param stepSize ODE integrator step size.
840 * @param stats Pointer to stats struct.
841 * @param fastStateUpdates Number of fast state updates.
842 */
843 ✗ void logSolverStats(enum OMC_LOG_STREAM stream, const char* name, double timeValue, double integratorTime, double stepSize, SOLVERSTATS* stats, unsigned int *fastStateUpdates, unsigned int *additionalEvalsFODE)
844 {
845 ✗ if (OMC_ACTIVE_STREAM(stream)) {
846 ✗ infoStreamPrint(stream, 1, "%s call statistics:", name);
847 ✗ infoStreamPrint(stream, 0, "number of steps taken so far: %d", stats->nStepsTaken);
848 ✗ infoStreamPrint(stream, 0, "number of calls of functionODE() : %d", stats->nCallsODE);
849 ✗ infoStreamPrint(stream, 0, "number of calculation of jacobian : %d", stats->nCallsJacobian);
850 ✗ infoStreamPrint(stream, 0, "error test failure : %d", stats->nErrorTestFailures);
851 ✗ infoStreamPrint(stream, 0, "convergence failure : %d", stats->nConvergenceTestFailures);
852 ✗ if (fastStateUpdates != NULL) infoStreamPrint(stream, 0, "number of fast state updates : %d", *fastStateUpdates);
853 ✗ if (additionalEvalsFODE != NULL) infoStreamPrint(stream, 0, "number of additional full calls of functionODE() : %d", *additionalEvalsFODE);
854 ✗ messageClose(stream);
855 }
856 ✗ }
857
858 /**
859 * @brief Info message for GBODE replacement.
860 *
861 * Dumps simulation flags to use to OMC_LOG_STDOUT.
862 *
863 * @param gbMethod GBODE method to use.
864 * @param constant If true use constant step size.
865 */
866 ✗ void replacementString(enum GB_METHOD gbMethod, modelica_boolean constant)
867 {
868 ✗ if (constant) {
869 ✗ infoStreamPrint(OMC_LOG_STDOUT, 1, "Use integration method GBODE with method '%s' and constant step size instead:", GB_METHOD_NAME[gbMethod]);
870 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Choose integration method '%s' in Simulation Setup->General and additional simulation flags '-%s=%s -%s=%s' in Simulation Setup->Simulation Flags.",
871 SOLVER_METHOD_NAME[S_GBODE], FLAG_NAME[FLAG_SR], GB_METHOD_NAME[gbMethod], FLAG_NAME[FLAG_SR_CTRL], GB_CTRL_METHOD_NAME[GB_CTRL_CNST]);
872 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "or");
873 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Simulation flags '-s=%s -%s=%s -%s=%s'.",
874 SOLVER_METHOD_NAME[S_GBODE], FLAG_NAME[FLAG_SR], GB_METHOD_NAME[gbMethod], FLAG_NAME[FLAG_SR_CTRL], GB_CTRL_METHOD_NAME[GB_CTRL_CNST]);
875 } else {
876 ✗ infoStreamPrint(OMC_LOG_STDOUT, 1, "Use integration method GBODE with method '%s' instead:", GB_METHOD_NAME[gbMethod]);
877 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Choose integration method '%s' in Simulation Setup->General and additional simulation flags '-%s=%s' in Simulation Setup->Simulation Flags.",
878 SOLVER_METHOD_NAME[S_GBODE], FLAG_NAME[FLAG_SR], GB_METHOD_NAME[gbMethod]);
879 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "or");
880 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Simulation flags '-s=%s -%s=%s'.",
881 SOLVER_METHOD_NAME[S_GBODE], FLAG_NAME[FLAG_SR], GB_METHOD_NAME[gbMethod]);
882 }
883 ✗ messageClose(OMC_LOG_STDOUT);
884 ✗ }
885
886 /**
887 * @brief Display deprecation warning for integration methods replaced by GBODE.
888 *
889 * Deprecated methods: None
890 *
891 * @param solverMethod Integration method.
892 */
893 1 void deprecationWarningGBODE(enum SOLVER_METHOD method)
894 {
895
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 switch (method) {
896 case S_RUNGEKUTTA:
897 break;
898 default:
899 return;
900 }
901
902 ✗ warningStreamPrint(OMC_LOG_STDOUT, 1, "Integration method '%s' is deprecated and will be removed in a future version of OpenModelica.", SOLVER_METHOD_NAME[method]);
903 switch (method) {
904 case S_RUNGEKUTTA:
905 ✗ replacementString(RK_RUNGEKUTTA, TRUE);
906 break;
907 default:
908 throwStreamPrint(NULL, "Not reachable state");
909 }
910
911 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0 , "See OpenModelica User's Guide section on GBODE for more details: https://www.openmodelica.org/doc/OpenModelicaUsersGuide/latest/solving.html#gbode");
912 ✗ messageCloseWarning(OMC_LOG_STDOUT);
913 ✗ return;
914 }
915
916 ✗ SLOW_STATE_CACHE *slowStateCache_alloc(int n_stages, int n_states, double *c)
917 {
918 ✗ SLOW_STATE_CACHE *cache = (SLOW_STATE_CACHE*) malloc(sizeof(SLOW_STATE_CACHE));
919 ✗ cache->n_stages = n_stages;
920 ✗ cache->n_states = n_states;
921 ✗ cache->states = (double *) malloc((n_stages + 2) * n_states * sizeof(double));
922 ✗ cache->work = (double *) malloc(n_states * sizeof(double));
923 ✗ cache->valid = (modelica_boolean *) calloc(n_stages + 2, sizeof(modelica_boolean));
924 ✗ cache->offset = 0;
925
926 ✗ cache->left_stage = n_stages;
927 ✗ cache->right_stage = n_stages + 1;
928
929 ✗ for (int i = 0; i < n_stages; i++)
930 {
931 ✗ if (c[i] == 0.0) cache->left_stage = i;
932 ✗ if (c[i] == 1.0) cache->right_stage = i;
933 }
934 ✗ return cache;
935 }
936
937 ✗ void slowStateCache_free(SLOW_STATE_CACHE *cache)
938 {
939 ✗ free(cache->states);
940 ✗ free(cache->work);
941 ✗ free(cache->valid);
942 ✗ free(cache);
943 ✗ }
944
945 // physical slot for logical index i
946 static inline int slowStateCache_slot(SLOW_STATE_CACHE *cache, int i)
947 {
948 ✗ return (cache->offset + i) % (cache->n_stages + 2);
949 }
950
951 /* simple interpolation wrapper that uses
952 * - GBODE interpolation method
953 * - big step solution -> interpolate to some provided time value
954 * - writes solution to out buffer (field in cache structure)
955 * - uses indirect indexing only if ratio of slow states to all states is very small,
956 * otherwise interpolates all states from the slow step
957 * => the latter behavior is not guaranteed and only done to reduce work */
958 ✗ static inline void slowStateCache_interpolate_slow_to_fast_node(DATA_GBODE *gbData, double time_value, double *out)
959 {
960 ✗ modelica_boolean use_sparse_slow_interp = ((double) gbData->nSlowStates / (double) gbData->nStates < 0.2);
961 ✗ gb_interpolation(gbData->gbfData->interpolation,
962 gbData->timeLeft, gbData->yLeft, gbData->kLeft,
963 gbData->timeRight, gbData->yRight, gbData->kRight,
964 time_value, out,
965 use_sparse_slow_interp ? gbData->nSlowStates : 0, use_sparse_slow_interp ? gbData->slowStatesIdx : NULL,
966 gbData->nStates, gbData->tableau, gbData->x, gbData->k);
967 ✗ }
968
969 // interpolate stage node if not cached, return pointer to states slot
970 ✗ static inline double *slowStateCache_get_or_compute_stage(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, int stage)
971 {
972 int s = slowStateCache_slot(cache, stage);
973 ✗ if (!cache->valid[s])
974 {
975 ✗ double t_stage = gbData->gbfData->time + gbData->gbfData->tableau->c[stage] * gbData->gbfData->stepSize;
976 ✗ slowStateCache_interpolate_slow_to_fast_node(gbData, t_stage, &cache->states[s * cache->n_states]);
977 ✗ cache->valid[s] = TRUE;
978 }
979 ✗ return &cache->states[s * cache->n_states];
980 }
981
982 // interpolate left boundary if not cached, return pointer to states slot
983 ✗ static inline double *slowStateCache_get_or_compute_left(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache)
984 {
985 ✗ int s = slowStateCache_slot(cache, cache->left_stage);
986 ✗ if (!cache->valid[s])
987 {
988 ✗ slowStateCache_interpolate_slow_to_fast_node(gbData, gbData->gbfData->time, &cache->states[s * cache->n_states]);
989 ✗ cache->valid[s] = TRUE;
990 }
991 ✗ return &cache->states[s * cache->n_states];
992 }
993
994 // interpolate right boundary if not cached, return pointer to states slot
995 ✗ static inline double *slowStateCache_get_or_compute_right(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache)
996 {
997 ✗ int s = slowStateCache_slot(cache, cache->right_stage);
998 ✗ if (!cache->valid[s])
999 {
1000 ✗ slowStateCache_interpolate_slow_to_fast_node(gbData, gbData->gbfData->time + gbData->gbfData->stepSize, &cache->states[s * cache->n_states]);
1001 ✗ cache->valid[s] = TRUE;
1002 }
1003 ✗ return &cache->states[s * cache->n_states];
1004 }
1005
1006 // write all slow states of interp to x, preserve fast states using work buffer
1007 ✗ static inline void slowStateCache_merge(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, double *interp, double *x)
1008 {
1009 ✗ if (((double) gbData->nSlowStates / (double) gbData->nStates < 0.2))
1010 {
1011 // do not use memcpy if its really sparse
1012 ✗ for (int i = 0; i < gbData->nSlowStates; i++)
1013 {
1014 ✗ int full_idx = gbData->slowStatesIdx[i];
1015 ✗ x[full_idx] = interp[full_idx];
1016 }
1017 }
1018 else
1019 {
1020 ✗ for (int i = 0; i < gbData->nFastStates; i++) cache->work[i] = x[gbData->fastStatesIdx[i]];
1021 ✗ memcpy(x, interp, cache->n_states * sizeof(double));
1022 ✗ for (int i = 0; i < gbData->nFastStates; i++) x[gbData->fastStatesIdx[i]] = cache->work[i];
1023 }
1024 ✗ }
1025
1026 // write all slow states of interp to x, fast states are potentially overwritten
1027 ✗ static inline void slowStateCache_overwrite(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, double *interp, double *x)
1028 {
1029 ✗ if (((double) gbData->nSlowStates / (double) gbData->nStates < 0.2))
1030 {
1031 // well we do not fill full, if its very sparse
1032 ✗ for (int i = 0; i < gbData->nSlowStates; i++)
1033 {
1034 ✗ int full_idx = gbData->slowStatesIdx[i];
1035 ✗ x[full_idx] = interp[full_idx];
1036 }
1037 }
1038 else
1039 {
1040 ✗ memcpy(x, interp, cache->n_states * sizeof(double));
1041 }
1042 ✗ }
1043
1044 // non-static, public functions
1045
1046 ✗ void slowStateCache_invalidate(SLOW_STATE_CACHE *cache)
1047 {
1048 ✗ memset(cache->valid, 0, sizeof(modelica_boolean) * (cache->n_stages + 2));
1049 ✗ cache->offset = 0;
1050 ✗ }
1051
1052 ✗ void slowStateCache_invalidate_keep_left(SLOW_STATE_CACHE *cache)
1053 {
1054 ✗ modelica_boolean carry = cache->valid[slowStateCache_slot(cache, cache->left_stage)];
1055 ✗ memset(cache->valid, 0, sizeof(modelica_boolean) * (cache->n_stages + 2));
1056 ✗ cache->valid[slowStateCache_slot(cache, cache->left_stage)] = carry;
1057 ✗ }
1058
1059 ✗ void slowStateCache_rotate(SLOW_STATE_CACHE *cache)
1060 {
1061 ✗ int left = cache->left_stage;
1062 ✗ int right = cache->right_stage;
1063
1064 // save old right valid flag before disabling all
1065 ✗ modelica_boolean carry = cache->valid[slowStateCache_slot(cache, right)];
1066
1067 // rotate offset so logical right maps to logical left:
1068 // new_offset + left ≡ old_offset + right (mod n_stages + 2)
1069 int divisor = cache->n_stages + 2;
1070 ✗ cache->offset = ((cache->offset + right - left + divisor) % divisor);
1071
1072 // disable all, restore carried left
1073 ✗ memset(cache->valid, 0, sizeof(modelica_boolean) * (cache->n_stages + 2));
1074 ✗ cache->valid[slowStateCache_slot(cache, left)] = carry;
1075 ✗ }
1076
1077 ✗ void slowStateCache_overwrite_stage(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, int stage, double *x)
1078 {
1079 ✗ double *interp = slowStateCache_get_or_compute_stage(gbData, cache, stage);
1080 ✗ slowStateCache_overwrite(gbData, cache, interp, x);
1081 ✗ }
1082
1083 ✗ void slowStateCache_overwrite_left(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, double *x)
1084 {
1085 ✗ double *interp = slowStateCache_get_or_compute_left(gbData, cache);
1086 ✗ slowStateCache_overwrite(gbData, cache, interp, x);
1087 ✗ }
1088
1089 ✗ void slowStateCache_overwrite_right(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, double *x)
1090 {
1091 ✗ double *interp = slowStateCache_get_or_compute_right(gbData, cache);
1092 ✗ slowStateCache_overwrite(gbData, cache, interp, x);
1093 ✗ }
1094
1095 ✗ void slowStateCache_merge_stage(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, int stage, double *x)
1096 {
1097 ✗ double *interp = slowStateCache_get_or_compute_stage(gbData, cache, stage);
1098 ✗ slowStateCache_merge(gbData, cache, interp, x);
1099 ✗ }
1100
1101 ✗ void slowStateCache_merge_left(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, double *x)
1102 {
1103 ✗ double *interp = slowStateCache_get_or_compute_left(gbData, cache);
1104 ✗ slowStateCache_merge(gbData, cache, interp, x);
1105 ✗ }
1106
1107 ✗ void slowStateCache_merge_right(DATA_GBODE *gbData, SLOW_STATE_CACHE *cache, double *x)
1108 {
1109 ✗ double *interp = slowStateCache_get_or_compute_right(gbData, cache);
1110 ✗ slowStateCache_merge(gbData, cache, interp, x);
1111 ✗ }
1112