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 / 141
Functions: 0.0% 0 / 1 / 13
Branches: 0.0% 0 / 0 / 114
OMCompiler/SimulationRuntime/c/simulation/solver/gbode_err.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 "gbode_err.h"
29 #include "gbode_internal_nls.h"
30
31 #include <float.h>
32 #include <math.h>
33
34 /* some constants for less verbose BLAS calls */
35 static const double DBL_ZERO = 0.0;
36 static const double DBL_ONE = 1.0;
37 static const double DBL_MINUS_ONE = -1.0;
38 static const int INT_ONE = 1;
39 static const char CHAR_NO_TRANS = 'N';
40
41 /* y := a * x + y */
42 extern void daxpy_(const int *n,
43 const double *alpha,
44 const double *x, const int *incX,
45 double *y, const int *incY);
46
47 /* y := alpha * A * x + beta * y */
48 extern void dgemv_(const char *trans,
49 const int *m,
50 const int *n,
51 const double *alpha, const double *A, const int *ldA,
52 const double *x, const int *incX,
53 const double *beta, double *y, const int *incY
54 );
55
56 /* x := alpha * x */
57 extern void dscal_(const int *n,
58 const double *alpha,
59 double *x, const int *incX);
60
61 static inline void setErrorEstimatorOrder(GB_ERROR_CONTEXT *context, int order)
62 {
63 ✗ if (context->isFast)
64 {
65 ✗ context->gbfData->currentErrorOrder = order;
66 }
67 else
68 {
69 ✗ context->gbData->currentErrorOrder = order;
70 }
71 }
72
73 static inline int evaluateError(GB_ERROR_CONTEXT *context, const GB_ERROR_ESTIMATOR *estimator)
74 {
75 ✗ if (estimator == NULL || estimator->type == GB_ERROR_UNKNOWN || estimator->evaluate == NULL)
76 {
77 return GB_ERROR_ESTIMATOR_FAILED;
78 }
79
80 ✗ return estimator->evaluate(context, estimator);
81 }
82
83 ✗ int gbEstimateError(GB_ERROR_CONTEXT *context, const GB_ERROR_ESTIMATOR *estimator)
84 {
85 int order = evaluateError(context, estimator);
86 ✗ if (order < 0)
87 {
88 ✗ return order;
89 }
90
91 setErrorEstimatorOrder(context, order);
92 return order;
93 }
94
95 ✗ double gbScaledErrorTolerance(double tol, int methodOrder, int estimatorOrder, modelica_boolean richardson)
96 {
97 ✗ if (richardson || estimatorOrder >= methodOrder)
98 {
99 return tol;
100 }
101
102 ✗ const double order_quot = ((double) estimatorOrder + 1.0) / ((double) methodOrder + 1.0);
103 ✗ const double rtol_pred = GB_TOLERANCE_SCALING_SAFETY * pow(tol, order_quot);
104
105 ✗ return fmax(tol, rtol_pred);
106 }
107
108 ✗ static void embeddedErrorEstimate_gb(BUTCHER_TABLEAU *tableau, const double *weights, const double *K, double stepSize, int nStates, double *err)
109 {
110 double factors[MAX_GBODE_STAGES];
111 ✗ int nStages = tableau->nStages;
112
113 ✗ for (int stage = 0; stage < nStages; stage++)
114 {
115 ✗ factors[stage] = stepSize * (tableau->b[stage] - weights[stage]);
116 }
117
118 /* err := stepSize * (K otimes I) * (b - weights) */
119 ✗ dgemv_(&CHAR_NO_TRANS,
120 &nStates,
121 &nStages,
122 &DBL_ONE, K, &nStates,
123 factors, &INT_ONE,
124 &DBL_ZERO, err, &INT_ONE);
125 ✗ }
126
127 static inline void absErrorEstimate_gb(int nStates, double *errest)
128 {
129 ✗ for (int i = 0; i < nStates; i++)
130 {
131 ✗ errest[i] = fabs(errest[i]);
132 }
133 }
134
135 ✗ static void embeddedErrorEstimate_gbf(BUTCHER_TABLEAU *tableau, const double *weights, DATA_GBODE *gbData, DATA_GBODEF *gbfData, const double *K, double *err)
136 {
137 ✗ int nStates = gbData->nStates;
138 ✗ int nFastStates = gbData->nFastStates;
139 ✗ int nStages = tableau->nStages;
140 double factors[MAX_GBODE_STAGES];
141
142 ✗ for (int stage = 0; stage < nStages; stage++)
143 {
144 ✗ factors[stage] = gbfData->stepSize * (tableau->b[stage] - weights[stage]);
145 }
146
147 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
148 {
149 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
150 ✗ err[full_idx] = 0.0;
151 ✗ for (int stage = 0; stage < nStages; stage++)
152 {
153 ✗ err[full_idx] += factors[stage] * K[stage * nStates + full_idx];
154 }
155 }
156 ✗ }
157
158 static void absErrorEstimate_gbf(DATA_GBODE *gbData, DATA_GBODEF *gbfData)
159 {
160 ✗ for (int fast_idx = 0; fast_idx < gbData->nFastStates; fast_idx++)
161 {
162 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
163 ✗ gbfData->errest[full_idx] = fabs(gbfData->errest[full_idx]);
164 }
165 }
166
167 ✗ int gbEmbeddedErrorEstimator(GB_ERROR_CONTEXT *context, const GB_ERROR_ESTIMATOR *estimator)
168 {
169 ✗ const double *weights = (const double *) estimator->data;
170
171 ✗ if (weights == NULL)
172 {
173 return GB_ERROR_ESTIMATOR_FAILED;
174 }
175
176 ✗ if (context->isFast)
177 {
178 ✗ embeddedErrorEstimate_gbf(context->gbfData->tableau, weights, context->gbData, context->gbfData, context->gbfData->k, context->gbfData->errest);
179 absErrorEstimate_gbf(context->gbData, context->gbfData);
180 }
181 else
182 {
183 ✗ DATA_GBODE *gbData = context->gbData;
184 ✗ embeddedErrorEstimate_gb(gbData->tableau, weights, gbData->k, gbData->stepSize, gbData->nStates, gbData->errest);
185 ✗ absErrorEstimate_gb(gbData->nStates, gbData->errest);
186 }
187
188 ✗ return estimator->order;
189 }
190
191 ✗ int gbContractiveDefectErrorEstimator(GB_ERROR_CONTEXT *context, const GB_ERROR_ESTIMATOR *estimator)
192 {
193 ✗ CONTRACTIVE_DEFECT *contractive = (CONTRACTIVE_DEFECT *) estimator->data;
194
195 ✗ if (contractive == NULL)
196 {
197 return GB_ERROR_ESTIMATOR_FAILED;
198 }
199
200 ✗ if (context->isFast)
201 {
202 ✗ DATA_GBODE *gbData = context->gbData;
203 ✗ DATA_GBODEF *gbfData = context->gbfData;
204
205 ✗ if (gbfData->tableau->t_transform == NULL || gbfData->nlsSolverMethod != GB_NLS_INTERNAL)
206 {
207 return GB_ERROR_ESTIMATOR_FAILED;
208 }
209
210 ✗ double *work = gbInternalGetWorkPointer(((struct dataSolver *)gbfData->nlsData->solverData)->ordinaryData);
211 ✗ gbInternalContractiveDefect(context->data, context->threadData, gbfData->nlsData, gbData, contractive, work);
212 ✗ for (int fast_idx = 0; fast_idx < gbData->nFastStates; fast_idx++)
213 {
214 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
215 ✗ gbfData->errest[full_idx] = fabs(work[fast_idx]);
216 }
217 }
218 else
219 {
220 ✗ DATA_GBODE *gbData = context->gbData;
221
222 ✗ if (gbData->tableau->t_transform == NULL || gbData->nlsSolverMethod != GB_NLS_INTERNAL)
223 {
224 return GB_ERROR_ESTIMATOR_FAILED;
225 }
226
227 ✗ gbInternalContractiveDefect(context->data, context->threadData, gbData->nlsData, gbData, contractive, gbData->errest);
228 ✗ absErrorEstimate_gb(gbData->nStates, gbData->errest);
229 }
230
231 ✗ return estimator->order;
232 }
233
234 ✗ int gbContractiveFilterErrorEstimator(GB_ERROR_CONTEXT *context, const GB_ERROR_ESTIMATOR *estimator)
235 {
236 ✗ if (context->isFast)
237 {
238 ✗ const double *weights = (const double *) context->gbfData->tableau->error.embedded.data;
239 ✗ if (weights == NULL || context->gbfData->nlsSolverMethod != GB_NLS_INTERNAL)
240 {
241 return GB_ERROR_ESTIMATOR_FAILED;
242 }
243 ✗ embeddedErrorEstimate_gbf(context->gbfData->tableau, weights, context->gbData, context->gbfData, context->gbfData->k, context->gbfData->errest);
244 ✗ gbInternalContractiveFilterError(context->gbfData->nlsData, context->gbData, context->gbfData->errest);
245 ✗ absErrorEstimate_gbf(context->gbData, context->gbfData);
246 }
247 else
248 {
249 ✗ DATA_GBODE *gbData = context->gbData;
250 ✗ const double *weights = (const double *) gbData->tableau->error.embedded.data;
251 ✗ if (weights == NULL || gbData->nlsSolverMethod != GB_NLS_INTERNAL)
252 {
253 return GB_ERROR_ESTIMATOR_FAILED;
254 }
255 ✗ embeddedErrorEstimate_gb(gbData->tableau, weights, gbData->k, gbData->stepSize, gbData->nStates, gbData->errest);
256 ✗ gbInternalContractiveFilterError(gbData->nlsData, gbData, gbData->errest);
257 ✗ absErrorEstimate_gb(gbData->nStates, gbData->errest);
258 }
259
260 ✗ return estimator->order;
261 }
262
263 ✗ static inline modelica_boolean twoStepScaleMu(double tol,
264 int methodOrder,
265 int estimatorOrder,
266 modelica_boolean richardson,
267 double *mu)
268 {
269 ✗ const double scaled_tol = gbScaledErrorTolerance(tol, methodOrder, estimatorOrder, richardson);
270 ✗ const double order_quot = ((double) estimatorOrder + 1.0) / ((double) methodOrder + 1.0);
271
272 ✗ *mu *= scaled_tol / pow(tol, order_quot);
273
274 // These cases should not occur for a well-conditioned two-step estimator:
275 // - num(r) -> 0, then mu(r) -> inf: the raw estimator loses its leading term and becomes locally one order higher than designed
276 // - den(r) -> 0, then mu(r) -> 0: the coefficient part already has a pole
277 // in both cases the estimator is invalid / inadequate
278
279 ✗ if (!isfinite(*mu))
280 {
281 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) warningStreamPrint(OMC_LOG_GBODE, 0, "Two-step estimator mu(r) is not finite - falling back.");
282 ✗ return FALSE;
283 }
284 ✗ else if (fabs(*mu) < 1e-6)
285 {
286 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) warningStreamPrint(OMC_LOG_GBODE, 0, "Two-step estimator mu(r) is below 1e-6 - falling back.");
287 ✗ return FALSE;
288 }
289 ✗ else if (fabs(*mu) > 1e6)
290 {
291 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) warningStreamPrint(OMC_LOG_GBODE, 0, "Two-step estimator mu(r) is above 1e6 - falling back.");
292 ✗ return FALSE;
293 }
294
295 return TRUE;
296 }
297
298 ✗ static int twoStepEstimate_gb(TWO_STEP_ESTIMATOR *two_step, DATA_GBODE *gbData, double tol, int estimatorOrder)
299 {
300 ✗ BUTCHER_TABLEAU *tableau = gbData->tableau;
301 double d_old[MAX_GBODE_FIRK_STAGES];
302 double g_new[MAX_GBODE_FIRK_STAGES];
303 double mu;
304 double minus_mu;
305 ✗ int nStages = tableau->nStages;
306 ✗ int nStates = gbData->nStates;
307
308 ✗ if (tableau->nStages > MAX_GBODE_FIRK_STAGES || gbData->lastStepSize <= 0.0 || gbData->extrapolationBaseTime == INFINITY || gbData->eventHappened || gbData->didFastStep)
309 {
310 return GB_ERROR_ESTIMATOR_FAILED;
311 }
312
313 ✗ double r = gbData->stepSize / gbData->lastStepSize;
314 ✗ two_step->weights(r, d_old, g_new, &mu);
315 ✗ if (!twoStepScaleMu(tol, tableau->order_b, estimatorOrder, tableau->richardson, &mu))
316 {
317 return GB_ERROR_ESTIMATOR_FAILED;
318 }
319
320 ✗ for (int stage = 0; stage < nStages; stage++)
321 {
322 ✗ d_old[stage] *= gbData->lastStepSize;
323 ✗ g_new[stage] *= gbData->stepSize;
324 }
325
326 ✗ dgemv_(&CHAR_NO_TRANS,
327 &nStates,
328 &nStages,
329 ✗ &DBL_ONE, gbData->kLast, &nStates,
330 d_old, &INT_ONE,
331 &DBL_ZERO, gbData->errest, &INT_ONE);
332
333 ✗ dgemv_(&CHAR_NO_TRANS,
334 &nStates,
335 &nStages,
336 ✗ &DBL_ONE, gbData->k, &nStates,
337 g_new, &INT_ONE,
338 &DBL_ONE, gbData->errest, &INT_ONE);
339
340 ✗ daxpy_(&nStates, &DBL_ONE, gbData->yOld, &INT_ONE, gbData->errest, &INT_ONE);
341
342 /* errest := abs(mu * (y - y_emb)) */
343 ✗ minus_mu = DBL_MINUS_ONE * mu;
344 ✗ dscal_(&nStates, &minus_mu, gbData->errest, &INT_ONE);
345 ✗ daxpy_(&nStates, &mu, gbData->y, &INT_ONE, gbData->errest, &INT_ONE);
346 ✗ absErrorEstimate_gb(nStates, gbData->errest);
347
348 return 0;
349 }
350
351 ✗ static int twoStepEstimate_gbf(TWO_STEP_ESTIMATOR *two_step, DATA_GBODE *gbData, DATA_GBODEF *gbfData, double tol, int estimatorOrder)
352 {
353 ✗ BUTCHER_TABLEAU *tableau = gbfData->tableau;
354 double d_old[MAX_GBODE_FIRK_STAGES];
355 double g_new[MAX_GBODE_FIRK_STAGES];
356 double mu;
357 ✗ int nStates = gbData->nStates;
358 ✗ int nFastStates = gbData->nFastStates;
359 ✗ int nStages = tableau->nStages;
360
361 ✗ if (tableau->nStages > MAX_GBODE_FIRK_STAGES || gbfData->lastStepSize <= 0.0 || !gbfData->extrapolationValid)
362 {
363 return GB_ERROR_ESTIMATOR_FAILED;
364 }
365
366 ✗ double r = gbfData->stepSize / gbfData->lastStepSize;
367 ✗ two_step->weights(r, d_old, g_new, &mu);
368 ✗ if (!twoStepScaleMu(tol, tableau->order_b, estimatorOrder, tableau->richardson, &mu))
369 {
370 return GB_ERROR_ESTIMATOR_FAILED;
371 }
372
373 ✗ for (int stage = 0; stage < nStages; stage++)
374 {
375 ✗ d_old[stage] *= gbfData->lastStepSize;
376 ✗ g_new[stage] *= gbfData->stepSize;
377 }
378
379 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
380 {
381 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
382 ✗ double y_emb = gbfData->yOld[full_idx];
383 ✗ for (int stage = 0; stage < nStages; stage++)
384 {
385 ✗ double k_new = gbfData->nlsSolverMethod == GB_NLS_INTERNAL
386 ✗ ? gbfData->kCurrPacked[stage * nFastStates + fast_idx]
387 ✗ : gbfData->k[stage * nStates + full_idx];
388 ✗ y_emb += d_old[stage] * gbfData->kLast[stage * nFastStates + fast_idx] + g_new[stage] * k_new;
389 }
390 ✗ gbfData->errest[full_idx] = fabs(mu * (gbfData->y[full_idx] - y_emb));
391 }
392
393 return 0;
394 }
395
396 ✗ int gbTwoStepErrorEstimator(GB_ERROR_CONTEXT *context, const GB_ERROR_ESTIMATOR *estimator)
397 {
398 ✗ TWO_STEP_ESTIMATOR *two_step = (TWO_STEP_ESTIMATOR *) estimator->data;
399
400 ✗ if (two_step == NULL)
401 {
402 return GB_ERROR_ESTIMATOR_FAILED;
403 }
404
405 ✗ const double tol = context->data->simulationInfo->tolerance;
406 ✗ int order = context->isFast ? twoStepEstimate_gbf(two_step, context->gbData, context->gbfData, tol, estimator->order)
407 ✗ : twoStepEstimate_gb(two_step, context->gbData, tol, estimator->order);
408 ✗ if (order >= 0)
409 {
410 ✗ return estimator->order;
411 }
412
413 ✗ return evaluateError(context, two_step->fallback);
414 }
415
416 ✗ int gbRichardsonErrorEstimator(GB_ERROR_CONTEXT *context, const GB_ERROR_ESTIMATOR *estimator)
417 {
418 (void) context; // suppress it.
419 ✗ return estimator->order;
420 }
421