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 / 1037
Functions: 0.0% 0 / 0 / 13
Branches: 0.0% 0 / 0 / 480
OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.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_main.c
29 * Implementation of a generic (implicit and explicit) Runge Kutta solver, which works for any
30 * order and stage based on a provided Butcher tableau. Utilizes the sparsity pattern of the ODE
31 * together with the KINSOL / KLU solver
32 *
33 * \author bbachmann
34 */
35
36 #include <time.h>
37
38 #include "gbode_main.h"
39 #include "gbode_util.h"
40
41 #include "gbode_conf.h"
42 #include "gbode_ctrl.h"
43 #include "gbode_err.h"
44 #include "gbode_events.h"
45 #include "gbode_nls.h"
46 #include "gbode_internal_nls.h"
47 #include "gbode_sparse.h"
48 #include "gbode_step.h"
49 #include "gbode_util.h"
50
51 #include <float.h>
52 #include <math.h>
53 #include <string.h>
54
55 #include "../arrayIndex.h"
56 #include "external_input.h"
57 #include "kinsolSolver.h"
58 #include "kinsol_b.h"
59 #include "newtonIteration.h"
60 #include "nonlinearSystem.h"
61 #include "omc_math.h"
62 #include "../options.h"
63 #include "../results/simulation_result.h"
64 #include "../jacobian_util.h"
65 #include "../../util/omc_error.h"
66 #include "../../util/omc_file.h"
67 #include "../../util/simulation_options.h"
68 #include "epsilon.h"
69
70 extern void communicateStatus(const char *phase, double completionPercent, double currentTime, double currentStepSize);
71
72 // TODO: we should add proper return handling of callbacks: ODE, Jacobian, Zero-Crossings, etc.
73 // TODO: make the interface between fast steps and slow steps more clear: It would be best to have one central function, which
74 // copies the required fields from fast -> slow
75
76 /**
77 * @brief Calculate function values of function ODE f(t,y).
78 *
79 * Assuming the correct values for time value and states are set.
80 *
81 * @param data Runtime data struct.
82 * @param threadData Thread data for error handling.
83 * @param counter Counter for function calls. Incremented by 1.
84 * @param selection Equations to evaluate.
85 */
86 ✗ int gbode_fODE(DATA *data, threadData_t *threadData, unsigned int* counter, EVAL_SELECTION* selection)
87 {
88 int ret = -1;
89 /* try */
90 #if !defined(OMC_EMCC)
91 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
92 #endif
93
94 ✗ if (counter)
95 {
96 ✗ (*counter)++;
97 }
98
99 ✗ externalInputUpdate(data);
100 ✗ data->callback->input_function(data, threadData);
101
102 ✗ data->simulationInfo->evalSelection = selection;
103 ✗ data->callback->functionODE(data, threadData);
104
105 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { ret = 0; }
106
107 #if !defined(OMC_EMCC)
108 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
109 #endif
110
111 ✗ return ret;
112 }
113
114 /**
115 * @brief Get the Jacobian method GBODE can use for the given non-linear solver.
116 *
117 * The adjoint and bidirectional evaluation modes produce a full Jacobian matrix in one
118 * go. GBODE's KINSOL and Newton non-linear solvers instead request single columns of the
119 * ODE Jacobian (see jacobian_SR_column() and friends in gbode_nls.c,
120 * These could be adapted of course aswell to allow adjoints), which is a genuine
121 * solver API restriction. Only the internal non-linear solver evaluates the whole ODE
122 * Jacobian at once and can therefore use all evaluation directions.
123 *
124 * @param threadData Used for error handling.
125 * @param nlsSolverMethod Non-linear solver method used by GBODE.
126 * @return JACOBIAN_METHOD Requested method, downgraded to the default if unusable.
127 */
128 ✗ static JACOBIAN_METHOD getGbodeJacobianMethod(threadData_t* threadData, enum GB_NLS_METHOD nlsSolverMethod)
129 {
130 ✗ JACOBIAN_METHOD jacobianMethod = getRequestedJacobianMethod(threadData);
131
132 /* non-internal non-linear solvers cannot use the adjoint or bidirectional Jacobian evaluation methods,
133 * because they only request single columns of the ODE Jacobian and not the whole matrix at once. */
134 ✗ if ((jacobianMethod == COLOREDSYMJACADJ || jacobianMethod == BICOLOREDSYMJAC)
135 ✗ && nlsSolverMethod != GB_NLS_INTERNAL) {
136 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Jacobian method %s requires the internal non-linear solver of GBODE. "
137 "Use `-gbnls=internal` / `-gbfnls=internal`. "
138 "Switching to the forward symbolic Jacobian.",
139 JACOBIAN_METHOD_NAME[jacobianMethod]);
140 jacobianMethod = JAC_UNKNOWN;
141 }
142
143 ✗ return jacobianMethod;
144 }
145
146 /**
147 * @brief Function allocates memory needed for chosen gbodef method.
148 *
149 * @param data Runtime data struct.
150 * @param threadData Thread data for error handling.
151 * @param solverInfo Information about main solver.
152 * @return int Return 0 on success, -1 on failure.
153 */
154 ✗ int gbodef_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, DATA_GBODE *gbData)
155 {
156 ✗ DATA_GBODEF *gbfData = (DATA_GBODEF *)calloc(1, sizeof(DATA_GBODEF));
157 ✗ gbData->gbfData = gbfData;
158
159 JACOBIAN *jacobian = NULL;
160 int i;
161
162 ✗ gbfData->nStates = gbData->nStates;
163
164 ✗ gbfData->GM_method = getGB_method(FLAG_MR);
165 ✗ gbfData->tableau = initButcherTableau(gbfData->GM_method, FLAG_MR_ERR);
166 ✗ if (gbfData->tableau == NULL) {
167 // ERROR
168 ✗ messageClose(OMC_LOG_STDOUT); // FIXME what does this belong to?
169 ✗ omc_throw_function(threadData);
170 }
171
172 // Get size of non-linear system
173 ✗ analyseButcherTableau(gbfData->tableau, gbData->nStates, &gbfData->nlSystemSize, &gbfData->type);
174
175 ✗ if (gbfData->GM_method == MS_ADAMS_MOULTON) {
176 ✗ gbfData->nlSystemSize = gbData->nStates;
177 ✗ gbfData->step_fun = &(full_implicit_MS_MR);
178 ✗ gbfData->type = MS_TYPE_IMPLICIT;
179 ✗ gbfData->isExplicit = FALSE;
180 }
181
182 ✗ switch (gbfData->type)
183 {
184 ✗ case GM_TYPE_EXPLICIT:
185 ✗ gbfData->isExplicit = TRUE;
186 ✗ gbfData->step_fun = &(expl_diag_impl_RK_MR);
187 ✗ break;
188 ✗ case GM_TYPE_DIRK:
189 ✗ gbfData->isExplicit = FALSE;
190 ✗ gbfData->step_fun = &(expl_diag_impl_RK_MR);
191 ✗ break;
192 ✗ case MS_TYPE_IMPLICIT:
193 ✗ gbfData->isExplicit = FALSE;
194 ✗ gbfData->step_fun = &(full_implicit_MS_MR);
195 ✗ break;
196 ✗ case GM_TYPE_IMPLICIT:
197 ✗ if (getGB_NLS_method(FLAG_MR_NLS) != GB_NLS_INTERNAL)
198 {
199 ✗ throwStreamPrint(NULL, "Unsupported configuration: fully implicit Runge-Kutta multirate integration is only available with -gbnls=internal.");
200 }
201 ✗ gbfData->isExplicit = FALSE;
202 ✗ gbfData->step_fun = &(full_implicit_RK_MR);
203 ✗ break;
204 ✗ default:
205 ✗ throwStreamPrint(NULL, "Not handled case for Runge-Kutta method %i", gbfData->type);
206 }
207
208 ✗ gbfData->nlsSolverMethod = gbfData->isExplicit ? GB_NLS_UNKNOWN : getGB_NLS_method(FLAG_MR_NLS);
209 ✗ finalizeButcherTableauError(gbfData->tableau, gbfData->nlsSolverMethod);
210
211 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Step control factor is set to %g", gbfData->tableau->fac);
212
213 ✗ gbfData->ctrl_method = getControllerMethod(FLAG_MR_CTRL);
214 ✗ if (gbfData->ctrl_method == GB_CTRL_CNST) {
215 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Constant step size not supported for inner integration. Using IController.");
216 ✗ gbfData->ctrl_method = GB_CTRL_I;
217 }
218 ✗ gbfData->currentErrorOrder = gbfData->tableau->error_order;
219
220 // allocate memory for the generic RK method
221 ✗ gbfData->y = malloc(gbData->nStates*sizeof(double));
222 ✗ gbfData->yOld = malloc(gbData->nStates*sizeof(double));
223 ✗ gbfData->yt = malloc(gbData->nStates*sizeof(double));
224 ✗ gbfData->y1 = malloc(gbData->nStates*sizeof(double));
225 ✗ gbfData->f = malloc(gbData->nStates*sizeof(double));
226 ✗ gbfData->yLast = malloc(gbData->nStates*sizeof(double));
227 ✗ gbfData->yOldPacked = malloc(gbData->nStates*sizeof(double));
228 ✗ gbfData->kLast = malloc(gbData->nStates*gbfData->tableau->nStages*sizeof(double));
229 ✗ gbfData->kCurrPacked = malloc(gbData->nStates*gbfData->tableau->nStages*sizeof(double));
230 ✗ gbfData->k = malloc(gbData->nStates*gbfData->tableau->nStages*sizeof(double));
231 ✗ gbfData->x = malloc(gbData->nStates*gbfData->tableau->nStages*sizeof(double));
232 ✗ gbfData->yLeft = malloc(gbData->nStates*sizeof(double));
233 ✗ gbfData->kLeft = malloc(gbData->nStates*sizeof(double));
234 ✗ gbfData->yRight = malloc(gbData->nStates*sizeof(double));
235 ✗ gbfData->kRight = malloc(gbData->nStates*sizeof(double));
236 ✗ gbfData->res_const = malloc(gbData->nStates*sizeof(double));
237 ✗ gbfData->errest = malloc(gbData->nStates*sizeof(double));
238 ✗ gbfData->errtol = malloc(gbData->nStates*sizeof(double));
239 ✗ gbfData->err = malloc(gbData->nStates*sizeof(double));
240 ✗ gbfData->ringBufferSize = 4;
241 ✗ gbfData->errValues = calloc(gbfData->ringBufferSize, sizeof(double));
242 ✗ gbfData->stepSizeValues = malloc(gbfData->ringBufferSize*sizeof(double));
243 ✗ gbfData->tv = malloc(gbfData->ringBufferSize*sizeof(double));
244 ✗ gbfData->yv = malloc(gbData->nStates*gbfData->ringBufferSize*sizeof(double));
245 ✗ gbfData->kv = malloc(gbData->nStates*gbfData->ringBufferSize*sizeof(double));
246 ✗ gbfData->slowStateCache = slowStateCache_alloc(gbfData->tableau->nStages, gbData->nStates, gbfData->tableau->c);
247
248 ✗ gbfData->extrapolationBaseTime = INFINITY;
249 ✗ gbfData->extrapolationStepSize = 0.0;
250 ✗ gbfData->extrapolationValid = FALSE;
251
252 ✗ gbData->nFastStates = 0;
253 ✗ gbData->nSlowStates = gbData->nFastStates;
254 ✗ gbfData->fastStates_old = malloc(gbData->nStates*sizeof(int));
255 ✗ gbfData->nFastStates_old = gbData->nFastStates;
256 ✗ for (int i = 0; i < gbData->nStates; i++) {
257 ✗ gbfData->fastStates_old[i] = i;
258 }
259
260 ✗ printButcherTableau(gbfData->tableau);
261
262 /* get DAG for functionODE */
263 ✗ data->callback->getDAG_ODE(data, threadData);
264 /* allocate selective RHS evaluation */
265 ✗ gbfData->evalSelectionFast = allocEvalSelection(data->modelData->dag);
266
267 /* initialize analytic Jacobian, if available and needed */
268 ✗ if (!gbfData->isExplicit) {
269 // Allocate Jacobian, if !gbfData->isExplcit and gbData->isExplicit
270 // Free is done in gbode_freeData
271 jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
272 ✗ if (gbData->isExplicit) {
273 ✗ JACOBIAN_METHOD jacobianMethod = getGbodeJacobianMethod(threadData, gbfData->nlsSolverMethod);
274 /* GBODE always needs the forward Jacobian A for its evaluation DAG and the
275 * multi-rate path so set requireForwardJacobian=TRUE. */
276 ✗ jacobian = initSymbolicOdeJacobian(data, threadData, &jacobianMethod, TRUE);
277 ✗ if (jacobian->availability != JACOBIAN_AVAILABLE && jacobian->availability != JACOBIAN_ONLY_SPARSITY) {
278 ✗ throwStreamPrint(threadData, "##GBODE## Implicit method requires a sparse pattern for the jacobian but no sparse pattern is generated.");
279 }
280
281 ✗ gbfData->symJacAvailable = jacobian->availability == JACOBIAN_AVAILABLE;
282 // change GBODE specific jacobian method
283 ✗ if (jacobianMethod == SYMJAC) {
284 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Symbolic Jacobians without coloring are currently not supported by GBODE."
285 " Colored symbolical Jacobian will be used.");
286 ✗ } else if(jacobianMethod == NUMJAC || jacobianMethod == COLOREDNUMJAC || jacobianMethod == INTERNALNUMJAC) {
287 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Numerical Jacobians without coloring are currently not supported by GBODE."
288 " Colored numerical Jacobian will be used.");
289 ✗ gbfData->symJacAvailable = FALSE;
290 }
291 } else {
292 ✗ gbfData->symJacAvailable = gbData->symJacAvailable;
293 ✗ jacobian = getSymbolicOdeJacobian(data);
294 }
295 /* The evaluation DAG is generated and consumed for the forward Jacobian A,
296 * even when the selected Jacobian evaluates adjoint directions. */
297 ✗ JACOBIAN* forwardJacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
298 ✗ if (forwardJacobian->availability == JACOBIAN_AVAILABLE) {
299 ✗ data->callback->getDAG_JacA(data, threadData, forwardJacobian);
300 }
301 ✗ if (!forwardJacobian->dag) {
302 ✗ throwStreamPrint(threadData,
303 "Cannot create multirate data structures without a valid Jacobian DAG. Use a symbolic Jacobian "
304 "(--generateDynamicJacobian=symbolic), an explicit integrator, or switch to single-rate integration.");
305 }
306
307 ✗ initializeSparsePattern_GBODEF(data, gbfData);
308
309 /* Initialize data for the nonlinear solver */
310 ✗ gbfData->nlsData = initRK_NLS_DATA_MR(data, threadData, gbfData);
311 ✗ if (!gbfData->nlsData) {
312 return -1;
313 }
314 } else {
315 ✗ gbfData->symJacAvailable = FALSE;
316 ✗ gbfData->nlsSolverMethod = GB_NLS_UNKNOWN;
317 ✗ gbfData->nlsData = NULL;
318 ✗ gbfData->jacobian = NULL;
319 }
320
321 ✗ gbfData->interpolation = getInterpolationMethod(FLAG_MR_INT);
322 ✗ if (!gbfData->tableau->withDenseOutput) {
323 ✗ if (gbfData->interpolation == GB_DENSE_OUTPUT) gbfData->interpolation = GB_INTERPOL_HERMITE;
324 }
325 ✗ switch (gbfData->interpolation)
326 {
327 ✗ case GB_INTERPOL_LIN:
328 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Linear interpolation is used for emitting results");
329 ✗ break;
330 ✗ case GB_INTERPOL_HERMITE:
331 case GB_INTERPOL_HERMITE_a:
332 case GB_INTERPOL_HERMITE_b:
333 case GB_INTERPOL_HERMITE_ERRCTRL:
334 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Hermite interpolation is used for the slow states");
335 ✗ break;
336 ✗ case GB_DENSE_OUTPUT:
337 case GB_DENSE_OUTPUT_ERRCTRL:
338 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Dense output is used for emitting results");
339 ✗ break;
340 ✗ default:
341 ✗ throwStreamPrint(NULL, "Unhandled interpolation case.");
342 }
343
344 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
345 enum { bufSize = 4096 };
346 char filename[bufSize];
347 ✗ snprintf(filename, bufSize, "%s_ActiveStates.txt", data->modelData->modelFilePrefix);
348 ✗ gbfData->fastStatesDebugFile = omc_fopen(filename, "w");
349 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "LOG_GBODE_STATES sets -noEquidistantTimeGrid for emitting results!");
350 ✗ solverInfo->solverNoEquidistantGrid = TRUE;
351 } else {
352 ✗ gbfData->fastStatesDebugFile = NULL;
353 }
354 ✗ i = (int) fmin(fmax(round(gbData->nStates * gbData->percentage), 1), gbData->nStates - 1);
355 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Number of states %d (%d slow states, %d fast states)", gbData->nStates, gbData->nStates-i, i);
356
357 /* reset statistics because it is accumulated in solver_main.c */
358 ✗ resetSolverStats(&gbfData->stats);
359 ✗ gbfData->fastStateUpdateCount = 0;
360 ✗ gbfData->additionalFullODEEvaluations = 0;
361
362 ✗ return 0;
363 }
364
365 /**
366 * @brief Read the states' nominal, min and max attributes into the solver data.
367 *
368 * Expensive scalar queries, so cached. Re-read by updateSolverNominals once
369 * initialization has computed the ones that are parameter expressions.
370 *
371 * @param data Runtime data struct.
372 * @param gbData Runge-Kutta solver data struct.
373 */
374 ✗ void gbode_setVarAttributes(DATA* data, DATA_GBODE* gbData)
375 {
376 ✗ for (int i = 0; i < gbData->nStates; i++) {
377 ✗ gbData->nominals[i] = fmax(fabs(getNominalFromScalarIdx(data->simulationInfo, data->modelData, VAR_KIND_STATE, i)), 1e-32);
378 ✗ gbData->mins[i] = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_STATE, i);
379 ✗ gbData->maxs[i] = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_STATE, i);
380 }
381 ✗ }
382
383 /**
384 * @brief Function allocates memory needed for generic RK method.
385 *
386 * @param data Runtime data struct.
387 * @param threadData Thread data for error handling.
388 * @param solverInfo Information about main solver.
389 * @return int Return 0 on success, -1 on failure.
390 */
391 ✗ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
392 {
393 ✗ DATA_GBODE *gbData = (DATA_GBODE *)calloc(1, sizeof(DATA_GBODE));
394
395 // Set backup in simulationInfo
396 ✗ data->simulationInfo->backupSolverData = (void *)gbData;
397
398 ✗ solverInfo->solverData = (void *)gbData;
399
400 ✗ gbData->nStates = data->modelData->nStates;
401
402 JACOBIAN* jacobian = NULL;
403
404 ✗ gbData->GM_method = getGB_method(FLAG_SR);
405 ✗ gbData->tableau = initButcherTableau(gbData->GM_method, FLAG_SR_ERR);
406 ✗ if (gbData->tableau == NULL) {
407 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "allocateDataGm: Failed to initialize gbode tableau for method %s", GB_METHOD_NAME[gbData->GM_method]);
408 ✗ return -1;
409 }
410
411 // Get size of non-linear system
412 ✗ analyseButcherTableau(gbData->tableau, gbData->nStates, &gbData->nlSystemSize, &gbData->type);
413
414 ✗ switch (gbData->type) {
415 ✗ case GM_TYPE_EXPLICIT:
416 ✗ gbData->isExplicit = TRUE;
417 ✗ gbData->step_fun = &(expl_diag_impl_RK);
418 ✗ break;
419 ✗ case GM_TYPE_DIRK:
420 ✗ gbData->isExplicit = FALSE;
421 ✗ gbData->step_fun = &(expl_diag_impl_RK);
422 ✗ break;
423 ✗ case GM_TYPE_IMPLICIT:
424 ✗ gbData->isExplicit = FALSE;
425 ✗ gbData->step_fun = &(full_implicit_RK);
426 ✗ break;
427 ✗ case MS_TYPE_IMPLICIT:
428 ✗ gbData->isExplicit = FALSE;
429 ✗ gbData->step_fun = &(full_implicit_MS);
430 ✗ break;
431 ✗ default:
432 ✗ throwStreamPrint(NULL, "gbode_allocateData: Unknown type %i", gbData->type);
433 }
434 ✗ if (gbData->GM_method == MS_ADAMS_MOULTON) {
435 ✗ gbData->nlSystemSize = gbData->nStates;
436 ✗ gbData->step_fun = &(full_implicit_MS);
437 ✗ gbData->type = MS_TYPE_IMPLICIT;
438 ✗ gbData->isExplicit = FALSE;
439 }
440
441 ✗ gbData->nlsSolverMethod = gbData->isExplicit ? GB_NLS_UNKNOWN : getGB_NLS_method(FLAG_SR_NLS);
442 ✗ finalizeButcherTableauError(gbData->tableau, gbData->nlsSolverMethod);
443
444 // detect controller method
445 ✗ gbData->ctrl_method = getControllerMethod(FLAG_SR_CTRL);
446 ✗ gbData->currentErrorOrder = gbData->tableau->error_order;
447 ✗ use_fhr = (modelica_boolean) omc_flag[FLAG_SR_CTRL_FHR];
448 ✗ use_filter = getGBCtrlFilterValue();
449
450 /* define maximum step size gbode is allowed to go */
451 ✗ if (omc_flag[FLAG_MAX_STEP_SIZE]) {
452 ✗ gbData->maxStepSize = atof(omc_flagValue[FLAG_MAX_STEP_SIZE]);
453 ✗ if (gbData->maxStepSize < 0 || gbData->maxStepSize > DBL_MAX/2) {
454 ✗ throwStreamPrint(NULL, "maximum step size %g is not allowed", gbData->maxStepSize);
455 } else {
456 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "maximum step size %g", gbData->maxStepSize);
457 }
458 } else {
459 ✗ gbData->maxStepSize = -1;
460 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "maximum step size not set");
461 }
462 /* Initial step size */
463 ✗ if (omc_flag[FLAG_INITIAL_STEP_SIZE]) {
464 ✗ gbData->initialStepSize = atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]);
465 ✗ if (gbData->initialStepSize < GB_MINIMAL_STEP_SIZE || gbData->initialStepSize > DBL_MAX/2) {
466 ✗ throwStreamPrint(NULL, "initial step size %g is not allowed, minimal step size is %g", gbData->initialStepSize, GB_MINIMAL_STEP_SIZE);
467 } else {
468 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "initial step size %g", gbData->initialStepSize);
469 }
470 } else {
471 ✗ gbData->initialStepSize = -1; /* use default */
472 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "initial step size not set");
473 }
474
475 /* if FLAG_NO_RESTART is set, configure gbode */
476 ✗ gbData->noRestart = omc_flag[FLAG_NO_RESTART];
477
478 ✗ gbData->eventTime = DBL_MAX;
479
480 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "gbode performs a restart after an event occurs %s", gbData->noRestart?"NO":"YES");
481
482 ✗ gbData->isFirstStep = TRUE;
483 ✗ gbData->didFastStep = FALSE;
484 ✗ gbData->eventHappened = FALSE;
485
486 /* mark initial extrapolation data as invalid () */
487 ✗ gbData->extrapolationBaseTime = INFINITY;
488 ✗ gbData->extrapolationStepSize = 0.0;
489
490 /* Allocate internal memory */
491 ✗ gbData->y = malloc(sizeof(double) * gbData->nStates);
492 ✗ gbData->yOld = malloc(sizeof(double) * gbData->nStates);
493 ✗ gbData->yLeft = malloc(sizeof(double) * gbData->nStates);
494 ✗ gbData->kLeft = malloc(sizeof(double) * gbData->nStates);
495 ✗ gbData->yRight = malloc(sizeof(double) * gbData->nStates);
496 ✗ gbData->kRight = malloc(sizeof(double) * gbData->nStates);
497 ✗ gbData->kLast = malloc(sizeof(double) * gbData->nStates * gbData->tableau->nStages);
498 ✗ gbData->yLast = malloc(sizeof(double) * gbData->nStates);
499 ✗ gbData->yt = malloc(sizeof(double) * gbData->nStates);
500 ✗ gbData->y1 = malloc(sizeof(double) * gbData->nStates);
501 ✗ gbData->y2 = malloc(sizeof(double) * gbData->nStates);
502 ✗ gbData->f = malloc(sizeof(double) * gbData->nStates);
503 ✗ gbData->k = malloc(sizeof(double) * gbData->nStates * gbData->tableau->nStages);
504 ✗ gbData->x = malloc(sizeof(double) * gbData->nStates * gbData->tableau->nStages);
505 ✗ gbData->res_const = malloc(sizeof(double) * gbData->nStates);
506 ✗ gbData->errest = malloc(sizeof(double) * gbData->nStates);
507 ✗ gbData->errtol = malloc(sizeof(double) * gbData->nStates);
508 ✗ gbData->err = malloc(sizeof(double) * gbData->nStates);
509 ✗ gbData->nominals = malloc(sizeof(double) * gbData->nStates);
510 ✗ gbData->mins = malloc(sizeof(double) * gbData->nStates);
511 ✗ gbData->maxs = malloc(sizeof(double) * gbData->nStates);
512 // ring buffer for different purposes (extrapolation, etc.)
513 ✗ gbData->ringBufferSize = 4;
514 ✗ gbData->errValues = malloc(sizeof(double) * gbData->ringBufferSize);
515 ✗ gbData->stepSizeValues = malloc(sizeof(double) * gbData->ringBufferSize);
516 ✗ gbData->tv = malloc(sizeof(double) * gbData->ringBufferSize);
517 ✗ gbData->yv = malloc(gbData->nStates*sizeof(double) * gbData->ringBufferSize);
518 ✗ gbData->kv = malloc(gbData->nStates*sizeof(double) * gbData->ringBufferSize);
519 ✗ gbData->tr = malloc(sizeof(double) * 2);
520 ✗ gbData->yr = malloc(gbData->nStates*sizeof(double) * 2);
521 ✗ gbData->kr = malloc(gbData->nStates*sizeof(double) * 2);
522
523 ✗ printButcherTableau(gbData->tableau);
524
525 ✗ gbode_setVarAttributes(data, gbData);
526
527 /* initialize analytic Jacobian, if available and needed */
528 ✗ if (!gbData->isExplicit) {
529 ✗ JACOBIAN_METHOD jacobianMethod = getGbodeJacobianMethod(threadData, gbData->nlsSolverMethod);
530 /* GBODE always needs the forward Jacobian A for its evaluation DAG and the
531 * multi-rate path, see gbInternal_evalJacobian() and initRK_NLS_DATA_MR(). */
532 ✗ jacobian = initSymbolicOdeJacobian(data, threadData, &jacobianMethod, TRUE);
533 ✗ if (jacobian->availability != JACOBIAN_AVAILABLE && jacobian->availability != JACOBIAN_ONLY_SPARSITY) {
534 ✗ throwStreamPrint(threadData, "##GBODE## Implicit method requires a sparse pattern for the jacobian but no sparse pattern is generated.");
535 }
536
537 ✗ gbData->symJacAvailable = jacobian->availability == JACOBIAN_AVAILABLE;
538 // change GBODE specific jacobian method
539 ✗ if (jacobianMethod == SYMJAC) {
540 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Symbolic Jacobians without coloring are currently not supported by GBODE."
541 " Colored symbolical Jacobian will be used.");
542 ✗ } else if (jacobianMethod == NUMJAC || jacobianMethod == COLOREDNUMJAC || jacobianMethod == INTERNALNUMJAC) {
543 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Numerical Jacobians without coloring are currently not supported by GBODE."
544 " Colored numerical Jacobian will be used.");
545 ✗ gbData->symJacAvailable = FALSE;
546 }
547
548 ✗ initializeSparsePattern_GBODE(data, gbData);
549
550 /* Initialize data for the nonlinear solver */
551 ✗ gbData->nlsData = initRK_NLS_DATA(data, threadData, gbData);
552 ✗ if (!gbData->nlsData) {
553 ✗ return -1;
554 } else {
555 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "Nominal values of the states:");
556 ✗ for (int i = 0; i < gbData->nStates; i++) {
557 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "%s = %g", data->modelData->realVarsData[i].info.name, gbData->nlsData->nominal[i]);
558 }
559 ✗ messageClose(OMC_LOG_SOLVER);
560 }
561 } else {
562 ✗ gbData->symJacAvailable = FALSE;
563 ✗ gbData->nlsSolverMethod = GB_NLS_UNKNOWN;
564 ✗ gbData->nlsData = NULL;
565 ✗ gbData->jacobian = NULL;
566 }
567
568 ✗ gbData->percentage = getGBRatio();
569 ✗ gbData->multi_rate = gbData->percentage > 0 && gbData->percentage < 1;
570
571 ✗ gbData->fastStatesIdx = malloc(sizeof(int) * gbData->nStates);
572 ✗ gbData->slowStatesIdx = malloc(sizeof(int) * gbData->nStates);
573 ✗ gbData->sortedStatesIdx = malloc(sizeof(int) * gbData->nStates);
574
575 ✗ gbData->nFastStates = 0;
576 ✗ gbData->nSlowStates = gbData->nStates;
577 ✗ for (int i = 0; i < gbData->nStates; i++) {
578 // TODO memcpy() faster?
579 ✗ gbData->fastStatesIdx[i] = i;
580 ✗ gbData->slowStatesIdx[i] = i;
581 ✗ gbData->sortedStatesIdx[i] = i;
582 }
583
584 ✗ if (gbData->multi_rate && omc_flagValue[FLAG_SR_INT]==NULL) {
585 ✗ gbData->interpolation = GB_DENSE_OUTPUT;
586 } else {
587 ✗ gbData->interpolation = getInterpolationMethod(FLAG_SR_INT);
588 }
589
590 ✗ if (!gbData->tableau->withDenseOutput) {
591 ✗ switch (gbData->interpolation) {
592 ✗ case GB_DENSE_OUTPUT: gbData->interpolation = GB_INTERPOL_HERMITE; break;
593 ✗ case GB_DENSE_OUTPUT_ERRCTRL: gbData->interpolation = GB_INTERPOL_HERMITE_ERRCTRL; break;
594 default: break;
595 }
596 }
597
598 enum { bufSize = 1024 };
599 char buffer[bufSize];
600 ✗ if (gbData->multi_rate) {
601 snprintf(buffer, bufSize, "%s", " and slow states interpolation");
602 } else {
603 snprintf(buffer, bufSize, "%s"," ");
604 }
605 ✗ switch (gbData->interpolation)
606 {
607 ✗ case GB_INTERPOL_LIN:
608 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Linear interpolation is used for emitting results%s", buffer);
609 ✗ break;
610 ✗ case GB_INTERPOL_HERMITE_ERRCTRL:
611 case GB_INTERPOL_HERMITE_a:
612 case GB_INTERPOL_HERMITE_b:
613 case GB_INTERPOL_HERMITE:
614 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Hermite interpolation is used for emitting results%s", buffer);
615 ✗ break;
616 ✗ case GB_DENSE_OUTPUT:
617 case GB_DENSE_OUTPUT_ERRCTRL:
618 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Dense output is used for emitting results%s", buffer);
619 ✗ break;
620 ✗ default:
621 ✗ throwStreamPrint(NULL, "Unhandled interpolation case.");
622 }
623 ✗ gbData->err_int = 0; // needed, if GB_INTERPOL_HERMITE_ERRCTRL or GB_DENSE_OUTPUT_ERRCTRL is used
624
625 ✗ if (gbData->multi_rate) {
626 ✗ gbodef_allocateData(data, threadData, solverInfo, gbData);
627 ✗ gbData->tableau->isKRightAvailable = FALSE;
628 } else {
629 ✗ gbData->gbfData = NULL;
630 }
631
632 // Value will be handled in the initial step size determination (-1 and 0 means no failure)
633 ✗ gbData->initialFailures = -1;
634
635 ✗ return 0;
636 }
637
638 /**
639 * @brief Free generic RK data.
640 *
641 * @param data Pointer to generik Runge-Kutta data struct.
642 */
643 ✗ void gbodef_freeData(DATA_GBODEF *gbfData)
644 {
645 ✗ freeEvalSelection(gbfData->evalSelectionFast);
646
647 /* Free non-linear system data */
648 ✗ freeRK_NLS_DATA(gbfData->nlsSolverMethod, gbfData->nlsData);
649
650 /* Free Jacobian */
651 ✗ freeJacobianCopy(gbfData->jacobian);
652
653 /* Free sparsity data. */
654 ✗ freeSparsePattern(gbfData->sparsePattern_ODE);
655 ✗ freeSparsePattern(gbfData->sparsePattern_NLS);
656 ✗ free(gbfData->sparseWork);
657
658 /* Free Butcher tableau */
659 ✗ freeButcherTableau(gbfData->tableau);
660
661 ✗ free(gbfData->y);
662 ✗ free(gbfData->yOld);
663 ✗ free(gbfData->yLeft);
664 ✗ free(gbfData->kLeft);
665 ✗ free(gbfData->yRight);
666 ✗ free(gbfData->kRight);
667 ✗ free(gbfData->yLast);
668 ✗ free(gbfData->yOldPacked);
669 ✗ free(gbfData->kLast);
670 ✗ free(gbfData->kCurrPacked);
671 ✗ free(gbfData->yt);
672 ✗ free(gbfData->y1);
673 ✗ free(gbfData->f);
674 ✗ free(gbfData->k);
675 ✗ free(gbfData->x);
676 ✗ free(gbfData->res_const);
677 ✗ free(gbfData->errest);
678 ✗ free(gbfData->errtol);
679 ✗ free(gbfData->err);
680 ✗ free(gbfData->errValues);
681 ✗ free(gbfData->stepSizeValues);
682 ✗ free(gbfData->tv);
683 ✗ free(gbfData->yv);
684 ✗ free(gbfData->kv);
685 ✗ free(gbfData->fastStates_old);
686
687 ✗ slowStateCache_free(gbfData->slowStateCache);
688
689 ✗ if (gbfData->fastStatesDebugFile)
690 ✗ fclose(gbfData->fastStatesDebugFile);
691
692 ✗ free(gbfData);
693 ✗ }
694
695 /**
696 * @brief Free generic RK data.
697 *
698 * @param gbData Pointer to generik Runge-Kutta data struct.
699 */
700 ✗ void gbode_freeData(DATA* data, DATA_GBODE *gbData)
701 {
702 ✗ freeSymbolicOdeJacobian(data);
703
704 /* Free non-linear system data */
705 ✗ freeRK_NLS_DATA(gbData->nlsSolverMethod, gbData->nlsData);
706
707 /* Free Jacobian */
708 ✗ freeJacobianCopy(gbData->jacobian);
709
710 /* Free sparsity data. */
711 ✗ freeSparsePattern(gbData->sparsePattern_NLS);
712
713 /* Free Butcher tableau */
714 ✗ freeButcherTableau(gbData->tableau);
715
716 ✗ if (gbData->multi_rate)
717 {
718 ✗ gbodef_freeData(gbData->gbfData);
719 gbData->gbfData = NULL;
720 }
721 /* Free multi-rate data */
722 ✗ free(gbData->err);
723 ✗ free(gbData->errValues);
724 ✗ free(gbData->stepSizeValues);
725 ✗ free(gbData->tv);
726 ✗ free(gbData->yv);
727 ✗ free(gbData->kv);
728 ✗ free(gbData->tr);
729 ✗ free(gbData->yr);
730 ✗ free(gbData->kr);
731 ✗ free(gbData->fastStatesIdx);
732 ✗ free(gbData->slowStatesIdx);
733 ✗ free(gbData->sortedStatesIdx);
734
735 /* Free remaining arrays */
736 ✗ free(gbData->y);
737 ✗ free(gbData->yOld);
738 ✗ free(gbData->yLast);
739 ✗ free(gbData->yLeft);
740 ✗ free(gbData->kLeft);
741 ✗ free(gbData->kLast);
742 ✗ free(gbData->yRight);
743 ✗ free(gbData->kRight);
744 ✗ free(gbData->yt);
745 ✗ free(gbData->y1);
746 ✗ free(gbData->y2);
747 ✗ free(gbData->f);
748 ✗ free(gbData->k);
749 ✗ free(gbData->x);
750 ✗ free(gbData->res_const);
751 ✗ free(gbData->errest);
752 ✗ free(gbData->errtol);
753 ✗ free(gbData->nominals);
754 ✗ free(gbData->mins);
755 ✗ free(gbData->maxs);
756
757 ✗ free(gbData);
758
759 ✗ return;
760 }
761
762 /**
763 * @brief Calculate initial step size.
764 *
765 * Called at the beginning of simulation or after an event occurred.
766 *
767 * @param data Runtime data struct.
768 * @param threadData Thread data for error handling.
769 * @param solverInfo Storing Runge-Kutta solver data.
770 */
771 ✗ void gbodef_init(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
772 {
773 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
774 SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
775 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
776 ✗ DATA_GBODEF* gbfData = gbData->gbfData;
777 ✗ int nStates = gbfData->nStates;
778 int nStages = gbfData->tableau->nStages;
779
780 int i;
781
782 ✗ gbfData->didEventStep = FALSE;
783 ✗ gbfData->extrapolationBaseTime = INFINITY;
784 ✗ gbfData->extrapolationValid = FALSE;
785 ✗ slowStateCache_invalidate(gbfData->slowStateCache);
786
787 ✗ gbfData->time = gbData->time;
788 ✗ gbfData->stepSize = 0.1*gbData->stepSize*GenericController(&(gbData->err_fast), &(gbData->stepSize), 1, GB_CTRL_I);
789
790 ✗ memcpy(gbfData->yOld, gbData->yOld, sizeof(double) * nStates);
791 ✗ memcpy(gbfData->y, gbData->y, sizeof(double) * nStates);
792
793 ✗ gbfData->timeRight = gbData->timeLeft;
794 ✗ memcpy(gbfData->yRight, gbData->yLeft, sizeof(double) * nStates);
795 ✗ memcpy(gbfData->kRight, gbData->kLeft, sizeof(double) * nStates);
796
797 // set solution ring buffer (extrapolation in case of NLS)
798 ✗ for (i = 0; i < gbfData->ringBufferSize; i++) {
799 ✗ gbfData->tv[i] = gbData->tv[i];
800 ✗ memcpy(gbfData->yv + i * nStates, gbData->yv + i * nStates, nStates * sizeof(double));
801 ✗ memcpy(gbfData->kv + i * nStates, gbData->kv + i * nStates, nStates * sizeof(double));
802 }
803 ✗ }
804
805 /**
806 * @brief Initialize ring buffer and interpolation arrays.
807 *
808 * Called at the beginning of simulation or after an event occurred.
809 *
810 * @param data Runtime data struct.
811 * @param threadData Thread data for error handling.
812 * @param solverInfo Storing Runge-Kutta solver data.
813 */
814 ✗ void gbode_init(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
815 {
816 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
817 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
818 ✗ modelica_real* fODE = &sData->realVars[gbData->nStates];
819 int nStates = gbData->nStates;
820 int i;
821
822 // initialize ring buffer for error and step size control
823 // TODO memset() faster?
824 ✗ for (i = 0; i < gbData->ringBufferSize; i++) {
825 ✗ gbData->errValues[i] = 0;
826 ✗ gbData->stepSizeValues[i] = 0;
827 }
828
829 /* reset statistics, because it is accumulated in solver_main.c */
830 ✗ if (!gbData->isExplicit)
831 ✗ gbData->nlsData->numberOfJEval = 0;
832 ✗ resetSolverStats(&gbData->stats);
833
834 // initialize vector used for interpolation (equidistant time grid)
835 // and for the birate inner integration
836 ✗ gbData->timeRight = gbData->time;
837 ✗ memcpy(gbData->yRight, gbData->yOld, nStates*sizeof(double));
838 ✗ memcpy(gbData->kRight, fODE, nStates*sizeof(double));
839
840 // set solution ring buffer (extrapolation in case of NLS)
841 ✗ for (i = 0; i < gbData->ringBufferSize; i++) {
842 ✗ gbData->tv[i] = gbData->timeRight;
843 ✗ memcpy(gbData->yv + i * nStates, gbData->yRight, nStates * sizeof(double));
844 ✗ memcpy(gbData->kv + i * nStates, gbData->kRight, nStates * sizeof(double));
845 }
846 ✗ gbData->eventTime = DBL_MAX; // reset event time
847 ✗ }
848
849 /*! \fn updateEvalSelection
850 *
851 * updates evalSelectionFast for evaluating gbode_fODE
852 */
853 ✗ static void updateEvalSelection(DATA* data, DATA_GBODE* gbData)
854 {
855 size_t k;
856 ✗ EVAL_SELECTION* selection = gbData->gbfData->evalSelectionFast;
857
858 ✗ clearEvalSelection(selection);
859
860 /* set equations for fast derivatives */
861 ✗ for (k = 0; k < gbData->nFastStates; k++) {
862 ✗ size_t derIdx = gbData->fastStatesIdx[k] + gbData->nStates;
863 ✗ size_t eqnIdx = selection->dag->mapVarToEqNode[derIdx];
864 ✗ if (eqnIdx == (size_t)(-1)) continue;
865 ✗ selection->dag->select[eqnIdx] = TRUE;
866 }
867
868 ✗ activateEvalDependencies(selection);
869
870 /* debug print */
871 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)) {
872 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 1, "updateEvalSelection");
873 enum { bufSize = 40960 };
874 char row_to_print[bufSize];
875 unsigned int ct;
876 ✗ ct = snprintf(row_to_print, bufSize, "%s (time=%g): =", "eqFunctions", data->localData[0]->timeValue);
877 ✗ for (k = 0; k < selection->n; k++) {
878 ✗ ct += snprintf(row_to_print+ct, bufSize-ct, " %zu", selection->idx[k]);
879 }
880 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 0, "%s", row_to_print);
881 ✗ messageClose(OMC_LOG_GBODE_V);
882 }
883 ✗ }
884
885 /*! \fn updateEvalSelectionJacobian
886 *
887 * updates evalSelection for evaluating Jacobian_ODE
888 */
889 ✗ static void updateEvalSelectionJacobian(DATA* data, DATA_GBODE* gbData)
890 {
891 size_t k;
892 ✗ EVAL_SELECTION* selection = gbData->gbfData->jacobian->evalSelection;
893
894 ✗ clearEvalSelection(selection);
895
896 /* set equations for fast derivatives */
897 ✗ for (k = 0; k < gbData->nFastStates; k++) {
898 ✗ size_t derIdx = gbData->fastStatesIdx[k];
899 ✗ size_t eqnIdx = selection->dag->mapVarToEqNode[derIdx];
900 ✗ if (eqnIdx == (size_t)(-1)) continue;
901 ✗ selection->dag->select[eqnIdx] = TRUE;
902 }
903
904 ✗ activateEvalDependencies(selection);
905
906 /* debug print */
907 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)) {
908 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 1, "updateEvalSelectionJacobian");
909 enum { bufSize = 40960 };
910 char row_to_print[bufSize];
911 unsigned int ct;
912 ✗ ct = snprintf(row_to_print, bufSize, "%s (time=%g): =", "Jacobian eqFunctions", data->localData[0]->timeValue);
913 ✗ for (k = 0; k < selection->n; k++) {
914 ✗ ct += snprintf(row_to_print+ct, bufSize-ct, " %zu", selection->idx[k]);
915 }
916 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 0, "%s", row_to_print);
917 ✗ messageClose(OMC_LOG_GBODE_V);
918 }
919 ✗ }
920
921 /*! \fn gbodef_main
922 *
923 * function does one integration step and calculates
924 * next step size by the implicit midpoint rule
925 *
926 * used for solver 'gm'
927 */
928 ✗ int gbodef_main(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, double targetTime)
929 {
930 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA *)data->localData[0];
931 ✗ modelica_real *fODE = sData->realVars + data->modelData->nStates;
932 ✗ DATA_GBODE *gbData = (DATA_GBODE *)solverInfo->solverData;
933 ✗ DATA_GBODEF *gbfData = gbData->gbfData;
934
935 ✗ double stopTime = data->simulationInfo->stopTime;
936
937 double err, eventTime;
938 double tol = data->simulationInfo->tolerance;
939
940 int i, ii, j, jj, l, ll, r, rr;
941 int integrator_step_info;
942
943 ✗ int nStates = gbData->nStates;
944 ✗ int nFastStates = gbData->nFastStates;
945 ✗ int nStages = gbfData->tableau->nStages;
946
947 modelica_boolean fastStatesChange = FALSE;
948 modelica_boolean foundEvent;
949
950 // This is the target time of the main integrator
951 ✗ const double innerTargetTime = fmin(targetTime, gbData->timeRight);
952
953 /* The inner integrator needs to be initialzed, at start time, when an event occured,
954 * and if outer integrations have been done with all states involved
955 * (gbfData->timeRight < gbData->timeLeft)
956 */
957 ✗ if (gbfData->didEventStep || gbfData->timeRight < gbData->timeLeft) {
958 ✗ gbodef_init(data, threadData, solverInfo);
959 }
960
961 ✗ fastStatesChange = checkFastStatesChange(gbData);
962
963 ✗ if (fastStatesChange) {
964 ✗ updateEvalSelection(data, gbData);
965 ✗ gbfData->extrapolationValid = FALSE;
966 ✗ gbfData->fastStateUpdateCount++;
967 }
968
969 ✗ if (fastStatesChange && !gbfData->isExplicit) {
970 ✗ struct dataSolver *solverData = gbfData->nlsData->solverData;
971 // set number of non-linear variables and corresponding nominal values (changes dynamically during simulation)
972 ✗ gbfData->nlsData->size = gbData->nFastStates;
973 ✗ slowStateCache_invalidate(gbfData->slowStateCache);
974
975 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Fast states and corresponding nominal values:");
976 ✗ for (ii = 0; ii < nFastStates; ii++) {
977 ✗ i = gbData->fastStatesIdx[ii];
978 ✗ gbfData->nlsData->nominal[ii] = gbData->nominals[i];
979 ✗ gbfData->nlsData->min[ii] = gbData->mins[i];
980 ✗ gbfData->nlsData->max[ii] = gbData->maxs[i];
981 ✗ infoStreamPrint(OMC_LOG_GBODE, 0, "%s = %g", data->modelData->realVarsData[i].info.name, gbfData->nlsData->nominal[ii]);
982 }
983 ✗ messageClose(OMC_LOG_GBODE);
984
985 ✗ if (gbfData->sparsePattern_NLS) {
986 ✗ updateSparsePattern_GBODEF(data, gbData);
987 ✗ if (gbfData->nlsSolverMethod != GB_NLS_INTERNAL)
988 {
989 ✗ gbfData->jacobian->sizeCols = nFastStates;
990 ✗ gbfData->jacobian->sizeRows = nFastStates;
991 }
992
993 /* TODO don't free and realloc, instead overwrite large enough buffer */
994 ✗ switch (gbfData->nlsSolverMethod)
995 {
996 ✗ case GB_NLS_NEWTON:
997 ✗ ((DATA_NEWTON *)solverData->ordinaryData)->n = gbData->nFastStates;
998 ✗ break;
999 ✗ case GB_NLS_KINSOL:
1000 ✗ nlsKinsolFree(solverData->ordinaryData);
1001 /* Set NLS user data */
1002 ✗ NLS_USERDATA* nlsUserData = initNlsUserData(data, threadData, -1, gbfData->nlsData, gbfData->jacobian);
1003 ✗ nlsUserData->solverData = (void*) gbfData;
1004 ✗ solverData->ordinaryData = (void*) nlsKinsolAllocate(gbfData->nlsData->size, nlsUserData, FALSE, !!gbfData->nlsData->sparsePattern);
1005 ✗ break;
1006 ✗ case GB_NLS_KINSOL_B:
1007 ✗ B_nlsKinsolFree(solverData->ordinaryData);
1008 /* Set NLS user data */
1009 ✗ NLS_USERDATA* B_nlsUserData = initNlsUserData(data, threadData, -1, gbfData->nlsData, gbfData->jacobian);
1010 ✗ B_nlsUserData->solverData = (void*) gbfData;
1011 ✗ solverData->ordinaryData = (void*) B_nlsKinsolAllocate(gbfData->nlsData->size, B_nlsUserData, FALSE, !!gbfData->nlsData->sparsePattern);
1012 ✗ break;
1013 ✗ case GB_NLS_INTERNAL:
1014 // notify internal to update the sparsity + symbolic factorization in the next iteration
1015 ✗ gbInternalScheduleFastStatesUpdate(solverData->ordinaryData);
1016 ✗ break;
1017 ✗ default:
1018 ✗ throwStreamPrint(NULL, "NLS method %s not yet implemented.", GB_NLS_METHOD_NAME[gbfData->nlsSolverMethod]);
1019 }
1020 }
1021 // TODO: -gbnls=internal currently does not use the Jacobian eval selection
1022 ✗ if (gbfData->nlsSolverMethod != GB_NLS_INTERNAL && gbfData->symJacAvailable) {
1023 ✗ updateEvalSelectionJacobian(data, gbData);
1024 }
1025 }
1026
1027 // print informations on the calling details
1028 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "gbodef solver started (fast states/states): %d/%d", gbData->nFastStates,gbData->nStates);
1029 ✗ printIntVector_gb(OMC_LOG_SOLVER, "fast States:", gbData->fastStatesIdx, gbData->nFastStates, gbfData->time);
1030 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "interpolation is done between %10g to %10g (SR-stepsize: %10g)",
1031 gbData->timeLeft, gbData->timeRight, gbData->lastStepSize);
1032
1033 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)) {
1034 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 1, "Interpolation values from outer integration:");
1035 ✗ printVector_gb(OMC_LOG_GBODE_V, "yL", gbData->yLeft, gbData->nStates, gbData->timeLeft);
1036 ✗ printVector_gb(OMC_LOG_GBODE_V, "kL", gbData->kLeft, gbData->nStates, gbData->timeLeft);
1037 ✗ printVector_gb(OMC_LOG_GBODE_V, "yR", gbData->yRight, gbData->nStates, gbData->timeRight);
1038 ✗ printVector_gb(OMC_LOG_GBODE_V, "kR", gbData->kRight, gbData->nStates, gbData->timeRight);
1039 ✗ messageClose(OMC_LOG_GBODE_V);
1040 }
1041
1042 ✗ while (gbfData->time < innerTargetTime) {
1043
1044 // Don't exceed simulation stop time
1045 ✗ if (gbfData->time + gbfData->stepSize > stopTime) {
1046 ✗ gbfData->stepSize = stopTime - gbfData->time;
1047 }
1048
1049 // Synchronize inner integration with outer integration
1050 // Strategy: either set outer step to the inner integration
1051 // or the other way around (depending on, if more or less
1052 // than 2 inner steps required)
1053 ✗ if (gbfData->time + gbfData->stepSize > gbData->timeRight) {
1054 // if (gbfData->time - gbfData->stepSize > gbData->timeLeft) {
1055 // gbData->timeRight = gbfData->timeRight;
1056 // gbData->lastStepSize = gbData->timeRight - gbData->timeLeft;
1057 // messageClose(OMC_LOG_SOLVER); // FIXME what does this belong to?
1058 // return 0;
1059 // } else {
1060 ✗ gbfData->stepSize = gbData->timeRight - gbfData->time;
1061 // }
1062 }
1063
1064 // store left hand data for later interpolation
1065 ✗ gbfData->timeLeft = gbfData->timeRight;
1066 ✗ memcpy(gbfData->yLeft, gbfData->yRight, nStates * sizeof(double));
1067 ✗ memcpy(gbfData->kLeft, gbfData->kRight, nStates * sizeof(double));
1068
1069 // debug the changes of the states and derivatives during integration
1070 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1071 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "states and derivatives at left hand side (inner integration):");
1072 ✗ printVector_gbf(OMC_LOG_GBODE, "yL", gbfData->yLeft, nStates, gbfData->timeLeft, gbData->nFastStates, gbData->fastStatesIdx);
1073 ✗ printVector_gbf(OMC_LOG_GBODE, "kL", gbfData->kLeft, nStates, gbfData->timeLeft, gbData->nFastStates, gbData->fastStatesIdx);
1074 ✗ messageClose(OMC_LOG_GBODE);
1075 }
1076
1077 do {
1078 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER_V)) {
1079 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "States and derivatives of the ring buffer:");
1080 ✗ for (int i = 0; i < gbfData->ringBufferSize; i++) {
1081 ✗ printVector_gbf(OMC_LOG_SOLVER_V, "y", gbfData->yv + i * nStates, nStates, gbfData->tv[i], gbData->nFastStates, gbData->fastStatesIdx);
1082 ✗ printVector_gbf(OMC_LOG_SOLVER_V, "k", gbfData->kv + i * nStates, nStates, gbfData->tv[i], gbData->nFastStates, gbData->fastStatesIdx);
1083 }
1084 ✗ messageClose(OMC_LOG_SOLVER_V);
1085 }
1086
1087 // do one integration step resulting in two different approximations
1088 // results are stored in gbData->y and gbData->yt
1089 ✗ if (gbfData->tableau->richardson) {
1090 ✗ integrator_step_info = gbodef_richardson(data, threadData, solverInfo);
1091 } else {
1092 ✗ integrator_step_info = gbfData->step_fun(data, threadData, solverInfo);
1093 }
1094
1095 // error handling: try half of the step size!
1096 ✗ if (integrator_step_info != 0) {
1097 ✗ (gbfData->stats).nConvergenceTestFailures++;
1098 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) infoStreamPrint(OMC_LOG_SOLVER, 0, "gbodef_main: Failed to calculate step at time = %5g with step size h = %5g.", gbData->time, gbData->stepSize);
1099 ✗ gbfData->stepSize *= 0.5;
1100 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) infoStreamPrint(OMC_LOG_SOLVER, 0, "Try half of the step size = %g", gbfData->stepSize);
1101 ✗ if (gbfData->stepSize < GB_MINIMAL_STEP_SIZE) {
1102 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but error still to large.", GB_MINIMAL_STEP_SIZE);
1103 ✗ messageClose(OMC_LOG_SOLVER); // FIXME what does this belong to?
1104 ✗ return -1;
1105 }
1106 ✗ slowStateCache_invalidate_keep_left(gbfData->slowStateCache);
1107 err = INFINITY;
1108 ✗ continue;
1109 }
1110
1111 ✗ tol = gbScaledErrorTolerance(data->simulationInfo->tolerance, gbfData->tableau->order_b,
1112 ✗ gbfData->currentErrorOrder, gbfData->tableau->richardson);
1113
1114 /* use same error estimate (scaled 2-norm) as for the SR case */
1115 ✗ for (i = 0, err=0; i < nFastStates; i++) {
1116 ✗ ii = gbData->fastStatesIdx[i];
1117 // calculate corresponding values for the error estimator and step size control
1118 ✗ gbfData->errtol[ii] = tol * gbData->nominals[ii] + fmax(fabs(gbfData->yOld[ii]), fabs(gbfData->y[ii])) * tol;
1119 ✗ if (gbfData->tableau->richardson || gbfData->type == MS_TYPE_IMPLICIT) {
1120 ✗ gbfData->errest[ii] = fabs(gbfData->yt[ii]);
1121 }
1122 ✗ gbfData->err[ii] = gbfData->tableau->fac * gbfData->errest[ii] / gbfData->errtol[ii];
1123 ✗ err += gbfData->err[ii] * gbfData->err[ii];
1124 }
1125 ✗ err = sqrt(err / (double) nFastStates);
1126
1127 // debug ring buffer for the states and derviatives of the states
1128 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)) {
1129 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 1, "ring buffer during steps of inner integration");
1130 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 0, "old value:");
1131 ✗ printVector_gbf(OMC_LOG_GBODE_V, "y", gbfData->yOld, nStates, gbfData->time, gbData->nFastStates, gbData->fastStatesIdx);
1132 ✗ debugRingBuffer_gbf(OMC_LOG_GBODE_V, gbfData->x, gbfData->k, nStates, gbfData->tableau, gbfData->time, gbfData->lastStepSize, gbData->nFastStates, gbData->fastStatesIdx);
1133 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 0, "new value:");
1134 ✗ printVector_gbf(OMC_LOG_GBODE_V, "y", gbfData->y, nStates, gbfData->time + gbfData->lastStepSize, gbData->nFastStates, gbData->fastStatesIdx);
1135 ✗ messageClose(OMC_LOG_GBODE_V);
1136 }
1137
1138 // Re-do step, if error is larger than requested
1139 ✗ if (err > 1) {
1140 ✗ gbfData->stats.nErrorTestFailures++;
1141 ✗ gbfData->stepSize *= 0.5;
1142 ✗ slowStateCache_invalidate_keep_left(gbfData->slowStateCache);
1143 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Reject step from %10g to %10g, error %10g, new stepsize %10g",
1144 ✗ gbfData->time, gbfData->time + gbfData->lastStepSize, err, gbfData->stepSize);
1145 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1146 ✗ dumpFastStates_gbf(gbData, gbfData->time + gbfData->lastStepSize, 1);
1147 }
1148 }
1149 ✗ } while (err > 1);
1150
1151 /* remember last time values for dense output extrapolation with yLast, kLast */
1152 ✗ gbfData->extrapolationBaseTime = gbfData->time;
1153 ✗ gbfData->extrapolationStepSize = gbfData->stepSize;
1154 ✗ gbfData->extrapolationValid = TRUE;
1155 ✗ gbData->didFastStep = TRUE;
1156
1157 /* remember kLast and yLast for dense output extrapolation, we have to pack them properly though
1158 TODO: gbfData->k and gbfData->y / yOld should always contain only the fast states packed from
1159 fast index 0, ..., nFastIndex - 1 (so no slow states) - we never need them, except for some interpolations
1160 but that is purely possible with data from gbData itself!! */
1161 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
1162 {
1163 ✗ int slow_idx = gbData->fastStatesIdx[fast_idx];
1164 ✗ gbfData->yLast[fast_idx] = gbfData->yOld[slow_idx];
1165 }
1166
1167 ✗ for (int stage = 0; stage < nStages; stage++)
1168 {
1169 ✗ double *kLast_strided = &gbfData->kLast[stage * nFastStates];
1170 ✗ double *k_strided = &gbfData->k[stage * nStates];
1171 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
1172 {
1173 ✗ int slow_idx = gbData->fastStatesIdx[fast_idx];
1174 ✗ kLast_strided[fast_idx] = k_strided[slow_idx];
1175 }
1176 }
1177
1178 // Count successful integration steps
1179 ✗ gbfData->stats.nStepsTaken += 1;
1180
1181 // pretty sure this is redundant (at least for ESDIRK) -> inspect this further
1182 ✗ slowStateCache_merge_left(gbData, gbfData->slowStateCache, gbfData->yOld);
1183
1184 // interpolation to right boundary
1185 ✗ slowStateCache_merge_right(gbData, gbfData->slowStateCache, gbfData->y);
1186 ✗ slowStateCache_rotate(gbfData->slowStateCache);
1187
1188 // store right hand values for latter interpolation
1189 ✗ gbfData->timeRight = gbfData->time + gbfData->stepSize;
1190 ✗ memcpy(gbfData->yRight, gbfData->y, nStates * sizeof(double));
1191 // update kRight
1192 ✗ if (!gbfData->tableau->isKRightAvailable) {
1193 ✗ sData->timeValue = gbfData->timeRight;
1194 ✗ memcpy(sData->realVars, gbfData->yRight, data->modelData->nStates * sizeof(double));
1195 ✗ gbode_fODE(data, threadData, &(gbData->gbfData->stats.nCallsODE), gbfData->evalSelectionFast);
1196 ✗ memcpy(gbfData->kRight, fODE, nStates * sizeof(double));
1197 }
1198 else
1199 {
1200 // last stage of method already provides the vector
1201 ✗ memcpy(gbfData->kRight, &gbfData->k[nStates * (nStages - 1)], nStates * sizeof(double));
1202 }
1203
1204 ✗ foundEvent = checkForEvents(data, threadData, solverInfo, gbfData->time, gbfData->yOld, gbfData->time + gbfData->stepSize, gbfData->y, TRUE, &eventTime);
1205 ✗ if (foundEvent) {
1206 ✗ solverInfo->currentTime = eventTime;
1207 ✗ sData->timeValue = solverInfo->currentTime;
1208 ✗ gbData->eventHappened = TRUE;
1209
1210 // sData->realVars are the "numerical" values on the right hand side of the event
1211 ✗ gbData->time = eventTime;
1212 ✗ memcpy(gbData->yOld, sData->realVars, gbData->nStates * sizeof(double));
1213
1214 ✗ gbfData->time = eventTime;
1215 ✗ memcpy(gbfData->yOld, sData->realVars, gbData->nStates * sizeof(double));
1216
1217 /* write statistics to the solverInfo data structure */
1218 ✗ memcpy(&solverInfo->solverStatsTmp, &gbfData->stats, sizeof(SOLVERSTATS));
1219
1220 // log the emitted result
1221 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)){
1222 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Emit result (inner integration):");
1223 ✗ printVector_gbf(OMC_LOG_GBODE, " y", sData->realVars, nStates, sData->timeValue, gbData->nFastStates, gbData->fastStatesIdx);
1224 ✗ messageClose(OMC_LOG_GBODE);
1225 }
1226
1227 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1228 ✗ dumpFastStates_gb(gbData, TRUE, eventTime, 0);
1229 }
1230
1231 // Get out of the integration routine for event handling
1232 ✗ messageClose(OMC_LOG_SOLVER); // FIXME what does this belong to?
1233 ✗ return 1;
1234 }
1235
1236 // set previously computed error (is accepted at this point) + predict new step size
1237 ✗ gbData->err_fast = err;
1238
1239 // Rotate and update buffer
1240 ✗ for (i = (gbfData->ringBufferSize - 1); i > 0 ; i--) {
1241 ✗ gbfData->errValues[i] = gbfData->errValues[i - 1];
1242 ✗ gbfData->stepSizeValues[i] = gbfData->stepSizeValues[i - 1];
1243 }
1244
1245 ✗ gbfData->errValues[0] = err;
1246 ✗ gbfData->stepSizeValues[0] = gbfData->stepSize;
1247
1248 /* update time with performed stepSize */
1249 ✗ gbfData->time += gbfData->stepSize;
1250
1251 // Store performed stepSize for adjusting the time in case of latter interpolation
1252 // Call the step size control
1253 ✗ gbfData->lastStepSize = gbfData->stepSize;
1254 ✗ gbfData->stepSize *= GenericController(gbfData->errValues, gbfData->stepSizeValues, gbfData->currentErrorOrder, gbfData->ctrl_method);
1255
1256 // debug the changes of the states and derivatives during integration
1257 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1258 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "States and derivatives at right hand side (inner integration):");
1259 ✗ printVector_gbf(OMC_LOG_GBODE, "yR", gbfData->yRight, nStates, gbfData->timeRight, gbData->nFastStates, gbData->fastStatesIdx);
1260 ✗ printVector_gbf(OMC_LOG_GBODE, "kR", gbfData->kRight, nStates, gbfData->timeRight, gbData->nFastStates, gbData->fastStatesIdx);
1261 ✗ messageClose(OMC_LOG_GBODE);
1262 }
1263
1264 // Rotate ring buffer
1265 ✗ for (i = (gbfData->ringBufferSize - 1); i > 0 ; i--) {
1266 ✗ gbfData->tv[i] = gbfData->tv[i - 1];
1267 ✗ memcpy(gbfData->yv + i * nStates, gbfData->yv + (i - 1) * nStates, nStates * sizeof(double));
1268 ✗ memcpy(gbfData->kv + i * nStates, gbfData->kv + (i - 1) * nStates, nStates * sizeof(double));
1269 }
1270
1271 ✗ gbfData->tv[0] = gbfData->timeRight;
1272 ✗ memcpy(gbfData->yv, gbfData->yRight, nStates * sizeof(double));
1273 ✗ memcpy(gbfData->kv, gbfData->kRight, nStates * sizeof(double));
1274
1275 ✗ debugRingBufferSteps_gbf(OMC_LOG_GBODE, gbfData->yv, gbfData->kv, gbfData->tv, nStates, gbfData->ringBufferSize, gbData->nFastStates, gbData->fastStatesIdx);
1276
1277 /* step is accepted and yOld needs to be updated */
1278 // copyVector_gbf(gbfData->yOld, gbfData->y, nFastStates, gbData->fastStates);
1279 ✗ memcpy(gbfData->yOld, gbfData->y, nStates * sizeof(double));
1280 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Accept step from %10g to %10g, error %10g, new stepsize %10g",
1281 ✗ gbfData->time - gbfData->lastStepSize, gbfData->time, err, gbfData->stepSize);
1282
1283 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1284 ✗ dumpFastStates_gbf(gbData, gbfData->time, 0);
1285 }
1286
1287 /* emit step, if solverNoEquidistantGrid is selected */
1288 ✗ if (solverInfo->solverNoEquidistantGrid) {
1289 ✗ sData->timeValue = gbfData->time;
1290 ✗ solverInfo->currentTime = sData->timeValue;
1291 ✗ memcpy(sData->realVars, gbfData->y, nStates * sizeof(double));
1292 /*
1293 * to emit consistent value we need to update the whole
1294 * continuous system with algebraic variables.
1295 */
1296 ✗ data->callback->updateContinuousSystem(data, threadData);
1297 ✗ sim_result.emit(&sim_result, data, threadData);
1298 // log the emitted result
1299 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)){
1300 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Emit result (inner integration):");
1301 ✗ printVector_gbf(OMC_LOG_GBODE, " y", sData->realVars, nStates, sData->timeValue, gbData->nFastStates, gbData->fastStatesIdx);
1302 ✗ messageClose(OMC_LOG_GBODE);
1303 }
1304 }
1305
1306 ✗ if ((gbData->timeRight - gbfData->time) < GB_MINIMAL_STEP_SIZE || gbData->stepSize < GB_MINIMAL_STEP_SIZE) {
1307 ✗ gbfData->time = gbData->timeRight;
1308 ✗ break;
1309 }
1310 }
1311
1312 // TODO: these 2 full ODE evaluations are a significant performance issue
1313 // I would advocate for removing the kv/yv/tv ring buffer and these 2 ODE evaluations.
1314 // For NLS initial guess extrapolation, yLast + kLast is sufficient as we already
1315 // have the stage derivatives of the last accepted step and a dense output function. If Hermite interpolation to a
1316 // specific point is needed, evaluate kLeft / kRight on demand or with a cache structure similar to slowStateCache.
1317 // For result file output, only state variables need to be set. The fODE is called
1318 // internally by the emit function anyway, so computing it explicitly here just to fill kv is redundant. (only to maintain the ring buffer)
1319
1320 /* update last two entries of ringbuffer with missing values of new fast derivatives */
1321 ✗ for (i = 0; i < 2; i++) {
1322 // TODO actually we only need fast derivatives, but at this point we don't
1323 // yet know which states will become fast so we compute everything.
1324 ✗ sData->timeValue = gbfData->tv[i];
1325 ✗ memcpy(sData->realVars, gbfData->yv + i * nStates, data->modelData->nStates * sizeof(double));
1326 ✗ gbode_fODE(data, threadData, &gbfData->additionalFullODEEvaluations, NULL);
1327
1328 ✗ for (j = 0; j < gbData->nSlowStates; j++)
1329 ✗ (gbfData->kv + i * nStates)[gbData->slowStatesIdx[j]] = fODE[gbData->slowStatesIdx[j]];
1330 }
1331
1332 // copy error and values of the fast states to the outer integrator routine if outer integration time is reached
1333 //gbData->err_fast = gbfData->errValues[0];
1334
1335 ✗ if (!solverInfo->solverNoEquidistantGrid && gbfData->time >= targetTime) {
1336 /* Integrator does large steps and needs to interpolate results with respect to the output grid */
1337 /* Here, only the fast states get updated */
1338 ✗ sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
1339 // solverInfo->currentTime = sData->timeValue;
1340
1341 ✗ gb_interpolation(gbfData->interpolation,
1342 gbfData->timeLeft, gbfData->yLeft, gbfData->kLeft,
1343 gbfData->timeRight, gbfData->yRight, gbfData->kRight,
1344 ✗ sData->timeValue, sData->realVars,
1345 nFastStates, gbData->fastStatesIdx, nStates, gbfData->tableau, gbfData->x, gbfData->k);
1346 }
1347 /* Solver statistics */
1348 ✗ if (!gbfData->isExplicit)
1349 {
1350 ✗ gbfData->stats.nCallsJacobian += gbfData->nlsData->numberOfJEval;
1351 ✗ gbfData->nlsData->numberOfJEval = 0;
1352 }
1353
1354 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "gbodef finished (inner steps).");
1355 ✗ messageClose(OMC_LOG_SOLVER); // FIXME what does this belong to?
1356
1357 ✗ return 0;
1358 }
1359
1360 /**
1361 * @brief Generic Runge-Kutta step.
1362 *
1363 * Do one Runge-Kutta integration step.
1364 * Has step-size control and event handling.
1365 *
1366 * @param data Runtime data struct.
1367 * @param threadData Thread data for error handling.
1368 * @param solverInfo Storing Runge-Kutta solver data.
1369 * @return int Return 0 on success, -1 on failure.
1370 */
1371 ✗ int gbode_main(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
1372 {
1373 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA *)data->localData[0];
1374 ✗ modelica_real *fODE = sData->realVars + data->modelData->nStates;
1375 ✗ DATA_GBODE *gbData = (DATA_GBODE *)solverInfo->solverData;
1376
1377 ✗ double stopTime = data->simulationInfo->stopTime;
1378 double tol = data->simulationInfo->tolerance;
1379
1380 ✗ int nStates = gbData->nStates;
1381 ✗ int nStages = gbData->tableau->nStages;
1382
1383 double targetTime, err;
1384
1385 ✗ const modelica_boolean noConst_intWithErrctrl = gbData->ctrl_method != GB_CTRL_CNST && (gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL || gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL);
1386
1387 int gb_step_info;
1388 int i, retries = 0;
1389 modelica_boolean foundEvent;
1390
1391 double err_states; // error of the (slow, if multirate) states
1392
1393 // root finding will be done in gbode after each accepted step
1394 ✗ solverInfo->solverRootFinding = 1;
1395
1396 /*
1397 * Determine the next target simulation time step.
1398 *
1399 * If the solver is using a non-equidistant grid:
1400 * → The target time is the minimum of the next sample event time
1401 * and the overall stop time.
1402 * Otherwise (equidistant grid):
1403 * → Start from the current time plus the step size,
1404 * but cap it by the stop time and the next scheduled event time.
1405 */
1406 ✗ if (solverInfo->solverNoEquidistantGrid) {
1407 // Non-equidistant grid: next step is driven by the nearest sample event.
1408 ✗ targetTime = fmin(data->simulationInfo->nextSampleEvent, stopTime);
1409 } else {
1410 // Equidistant output grid: targetTime set to the next output time.
1411 ✗ targetTime = solverInfo->currentTime + solverInfo->currentStepSize;
1412
1413 // Ensure we don't run past the stop time.
1414 ✗ targetTime = fmin(targetTime, stopTime);
1415
1416 // Also ensure we don't skip over an event time.
1417 ✗ targetTime = fmin(gbData->eventTime, targetTime);
1418 }
1419
1420 ✗ if (gbData->multi_rate) {
1421 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "Start gbode (birate integration) from %g to %g",
1422 solverInfo->currentTime, targetTime);
1423 } else {
1424 ✗ infoStreamPrint(OMC_LOG_SOLVER, 1, "Start gbode (single-rate integration) from %g to %g",
1425 solverInfo->currentTime, targetTime);
1426 }
1427
1428 ✗ gbData->eventHappened = solverInfo->didEventStep || gbData->isFirstStep;
1429
1430 /*
1431 * Handle step initialization after an event step or at the very first solver step.
1432 *
1433 * This section ensures that the solver’s time, step size, and related buffers
1434 * are correctly initialized before proceeding with integration.
1435 */
1436 ✗ if (solverInfo->didEventStep || gbData->isFirstStep) {
1437 ✗ if (gbData->noRestart && !gbData->isFirstStep) {
1438 /*
1439 * Case: No restart requested after event (-noRestart flag set)
1440 * and we are not at the very first step.
1441 * → Continue from the right boundary of the last interval
1442 * using the optimal step size determined earlier.
1443 */
1444 ✗ gbData->time = gbData->timeRight;
1445 ✗ gbData->stepSize = gbData->optStepSize;
1446 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0,
1447 "Initial step size = %e at time %g",
1448 gbData->stepSize, gbData->time);
1449 } else {
1450 /*
1451 * Case: Either restart is allowed OR this is the very first solver step.
1452 * → Recalculate the initial step size.
1453 * → Reset the ring buffer and solver statistics.
1454 * → Initialize gbData->timeRight, gbData->yRight, and gbData->kRight.
1455 */
1456 ✗ getInitStepSize(data, threadData, gbData, solverInfo);
1457 ✗ gbode_init(data, threadData, solverInfo);
1458 }
1459
1460 // Mark initialization as complete for this step
1461 ✗ gbData->isFirstStep = FALSE;
1462 ✗ solverInfo->didEventStep = FALSE;
1463
1464 // For multi-rate solvers, propagate event-step flag to the fine-level solver
1465 ✗ if (gbData->multi_rate) {
1466 ✗ gbData->gbfData->didEventStep = TRUE;
1467 }
1468 }
1469
1470 ✗ debugRingBufferSteps_gb(OMC_LOG_GBODE, gbData->yv, gbData->kv, gbData->tv, nStates, gbData->ringBufferSize);
1471
1472 /*
1473 * Case: Constant step size control method.
1474 * Use the solver's current step size directly without adjustment.
1475 */
1476 ✗ if (gbData->ctrl_method == GB_CTRL_CNST) {
1477 ✗ gbData->stepSize = solverInfo->currentStepSize;
1478 }
1479
1480
1481 ✗ if (gbData->multi_rate) {
1482 // Check if multirate step is necessary, otherwise the correct values are already stored in sData
1483 ✗ if (gbData->nFastStates > 0 && gbData->gbfData->time < gbData->timeRight && !gbData->gbfData->didEventStep) {
1484 // run multirate step
1485 ✗ gb_step_info = gbodef_main(data, threadData, solverInfo, targetTime);
1486 // synchronize y, yRight , kRight and buffer
1487 ✗ if (fabs(gbData->timeRight - gbData->gbfData->timeRight) < GB_MINIMAL_STEP_SIZE) {
1488 ✗ gbData->time = gbData->timeRight;
1489 ✗ memcpy(gbData->y, gbData->gbfData->y, nStates * sizeof(double));
1490 ✗ memcpy(gbData->yOld, gbData->y, nStates * sizeof(double));
1491 ✗ memcpy(gbData->yRight, gbData->gbfData->yRight, nStates * sizeof(double));
1492 ✗ memcpy(gbData->kRight, gbData->gbfData->kRight, nStates * sizeof(double));
1493 ✗ memcpy(gbData->err, gbData->gbfData->err, nStates * sizeof(double));
1494
1495 // update buffer, rest has already been rotated
1496 ✗ gbData->tv[0] = gbData->timeRight;
1497 ✗ memcpy(gbData->yv, gbData->yRight, nStates * sizeof(double));
1498 ✗ memcpy(gbData->kv, gbData->kRight, nStates * sizeof(double));
1499
1500 /* step is accepted and yOld needs to be updated */
1501 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Accept step from %10g to %10g, error slow states %10g, new stepsize %10g",
1502 ✗ gbData->time - gbData->lastStepSize, gbData->time, gbData->errValues[0], gbData->stepSize);
1503
1504 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1505 // dump fast states in file
1506 ✗ dumpFastStates_gb(gbData, FALSE, gbData->time, 0);
1507 }
1508 }
1509 ✗ if (gb_step_info !=0) {
1510 // get out of here, if an event has happend!
1511 ✗ messageClose(OMC_LOG_SOLVER);
1512 ✗ if (gb_step_info > 0)
1513 return 0;
1514 else
1515 ✗ return gb_step_info;
1516 }
1517 }
1518 }
1519
1520
1521 /* Main integration loop, if gbData->time already greater than targetTime, only the
1522 interpolation is necessary for emitting the output variables (see below) */
1523 ✗ while (gbData->time < targetTime) {
1524 /*
1525 * Limit the step size so we do not overshoot:
1526 * 1. The next sample event time
1527 * 2. The overall simulation stop time
1528 */
1529 ✗ gbData->stepSize = fmin(gbData->stepSize, data->simulationInfo->nextSampleEvent - gbData->time);
1530 ✗ gbData->stepSize = fmin(gbData->stepSize, stopTime - gbData->time);
1531 // TODO maybe easier to use targetTime
1532 //gbData->stepSize = fmin(gbData->stepSize, targetTime - gbData->time);
1533
1534 // Store the “left-hand side” data from the current step
1535 // for later use during interpolation.
1536 // Copies time, states, and derivatives from the “right” (current step)
1537 // to the “left” (previous step).
1538 // FIXME is this comment correct?
1539 ✗ gbData->timeLeft = gbData->timeRight;
1540 ✗ memcpy(gbData->yLeft, gbData->yRight, nStates * sizeof(double));
1541 ✗ memcpy(gbData->kLeft, gbData->kRight, nStates * sizeof(double));
1542
1543 // debug the ring buffer changes of the states and derivatives during integration
1544 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1545 // debug the changes of the states and derivatives during integration
1546 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "States and derivatives at left hand side:");
1547 ✗ printVector_gb(OMC_LOG_GBODE, "yL", gbData->yLeft, nStates, gbData->timeLeft);
1548 ✗ printVector_gb(OMC_LOG_GBODE, "kL", gbData->kLeft, nStates, gbData->timeLeft);
1549 ✗ messageClose(OMC_LOG_GBODE);
1550 }
1551
1552 // Loop will be performed until the error estimate for all states fullfills the
1553 // given tolerance
1554 do {
1555 // set error to INFINITY, in case we break / continue early
1556 err = INFINITY;
1557
1558 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER_V)) {
1559 // debug ring buffer of the states and derivatives during integration
1560 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "States and derivatives of the ring buffer:");
1561 ✗ for (int i=0; i<gbData->ringBufferSize; i++) {
1562 ✗ printVector_gb(OMC_LOG_SOLVER_V, "y", gbData->yv + i * nStates, nStates, gbData->tv[i]);
1563 }
1564 ✗ for (int i=0; i<gbData->ringBufferSize; i++) {
1565 ✗ printVector_gb(OMC_LOG_SOLVER_V, "k", gbData->kv + i * nStates, nStates, gbData->tv[i]);
1566 }
1567 ✗ messageClose(OMC_LOG_SOLVER_V);
1568 }
1569
1570 // Perform one integration step. New error estimators write |error| directly to errest;
1571 // Richardson and MS methods still write a signed error estimate to yt.
1572 // Choose the integration method based on the tableau:
1573 // - If Richardson extrapolation is enabled, use gbode_richardson.
1574 // - Otherwise, use the default step function stored in gbData->step_fun.
1575 ✗ if (gbData->tableau->richardson) {
1576 ✗ gb_step_info = gbode_richardson(data, threadData, solverInfo);
1577 } else {
1578 ✗ gb_step_info = gbData->step_fun(data, threadData, solverInfo);
1579 }
1580
1581 // debug the approximations after performed step
1582 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1583 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Approximations after step calculation:");
1584 ✗ printVector_gb(OMC_LOG_GBODE, " y", gbData->y, nStates, gbData->time + gbData->stepSize);
1585 ✗ if (gbData->tableau->richardson || gbData->type == MS_TYPE_IMPLICIT) {
1586 ✗ printVector_gb(OMC_LOG_GBODE, "yt", gbData->yt, nStates, gbData->time + gbData->stepSize);
1587 } else {
1588 ✗ printVector_gb(OMC_LOG_GBODE, "errest", gbData->errest, nStates, gbData->time + gbData->stepSize);
1589 }
1590 ✗ messageClose(OMC_LOG_GBODE);
1591 }
1592
1593 // Error handling for failed integration step:
1594 // If the step calculation failed (gb_step_info != 0), try reducing the step size by half and retry.
1595 //
1596 // Actions taken on failure:
1597 // - Increment convergence failure statistics counter.
1598 // - Print an informational message about the failure and the current simulation time.
1599 //
1600 // If the solver is using a constant step size control method:
1601 // - Abort the simulation and print an error message since no step size adjustment is possible.
1602 //
1603 // Otherwise (adaptive step size control):
1604 // - Halve the current step size.
1605 // - If multi-rate integration is active and detailed logging is enabled:
1606 // - Reset error metrics for slow, fast, and internal components.
1607 // - Dump the fast states to a file for diagnostics.
1608 // - Print the new reduced step size being tried.
1609 //
1610 // If the step size becomes smaller than the minimal allowed threshold:
1611 // - Abort the simulation with an error indicating minimum step size reached without acceptable error.
1612 //
1613 // If none of the abort conditions occur, the loop continues to retry with the reduced step size.
1614 ✗ if (gb_step_info != 0) {
1615 ✗ gbData->stats.nConvergenceTestFailures++;
1616 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) infoStreamPrint(OMC_LOG_SOLVER, 0, "gbode_main: Failed to calculate step at time = %5g with step size h = %5g.", gbData->time, gbData->stepSize);
1617 ✗ if (gbData->ctrl_method == GB_CTRL_CNST) {
1618 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Simulation aborted since gbode is running with fixed step size and step calculation has failed at time = %5g with step size h = %5g.", gbData->time, gbData->stepSize);
1619 ✗ messageClose(OMC_LOG_SOLVER);
1620 ✗ return -1;
1621 } else {
1622 ✗ if (gbData->eventHappened) {
1623 ✗ gbData->stepSize *= 0.1; // event or initial step rejection: reduce step size to 10% of previous step size
1624 } else {
1625 ✗ gbData->stepSize *= 0.5; // standard rejection: reduce the step size by half to attempt a more accurate integration in the next iteration
1626 }
1627
1628 ✗ if (gbData->multi_rate && OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1629 ✗ gbData->err_slow = 0;
1630 ✗ gbData->err_fast = 0;
1631 ✗ gbData->err_int = 0;
1632 // dump fast states in file
1633 ✗ dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->stepSize, 3);
1634 }
1635 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) infoStreamPrint(OMC_LOG_SOLVER, 0, "Try half of the step size = %g", gbData->stepSize);
1636 ✗ if (gbData->stepSize < GB_MINIMAL_STEP_SIZE) {
1637 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but error still to large.", GB_MINIMAL_STEP_SIZE);
1638 ✗ messageClose(OMC_LOG_SOLVER);
1639 ✗ return -1;
1640 }
1641 ✗ continue;
1642 }
1643 }
1644
1645 // Calculate error estimators and tolerance scaling for each state variable
1646 // Compute error tolerance for the i-th state based on relative and absolute tolerances:
1647 // errtol = Rtol * max(|old state|, |current state|) + Atol * |nominal(state)|
1648 ✗ tol = gbScaledErrorTolerance(data->simulationInfo->tolerance, gbData->tableau->order_b,
1649 ✗ gbData->currentErrorOrder, gbData->tableau->richardson);
1650
1651 ✗ for (i = 0, err=0; i < nStates; i++) {
1652 // calculate corresponding values for the error estimator and step size control
1653 ✗ gbData->errtol[i] = tol * gbData->nominals[i] + fmax(fabs(gbData->yOld[i]), fabs(gbData->y[i])) * tol;
1654 ✗ if (gbData->tableau->richardson || gbData->type == MS_TYPE_IMPLICIT) {
1655 ✗ gbData->errest[i] = fabs(gbData->yt[i]);
1656 }
1657 ✗ gbData->err[i] = gbData->tableau->fac * gbData->errest[i] / gbData->errtol[i];
1658 ✗ err += gbData->err[i] * gbData->err[i];
1659 }
1660
1661 ✗ err = sqrt(err / (double) nStates);
1662
1663 ✗ if (gbData->multi_rate) {
1664 // Multi-rate integration enabled:
1665
1666 // Calculate the error threshold for slow states (used to separate slow and fast states).
1667 ✗ err_states = getErrorThreshold(gbData);
1668 err = err_states;
1669
1670 // Classify states into fast and slow based on the scaled error:
1671 // - States with error >= 1 are considered fast.
1672 // - States with error < 1 are considered slow.
1673 //
1674 // Keep track of the count of fast and slow states,
1675 // and record the maximum error encountered for each group.
1676 ✗ gbData->nFastStates = 0;
1677 ✗ gbData->nSlowStates = 0;
1678 ✗ gbData->err_slow = 0;
1679 ✗ gbData->err_fast = 0;
1680 ✗ gbData->err_int = 0;
1681
1682 ✗ for (i = 0; i < gbData->nStates; i++) {
1683 ✗ if (gbData->err[i] >= 1) {
1684 ✗ gbData->fastStatesIdx[gbData->nFastStates] = i;
1685 ✗ gbData->nFastStates++;
1686 ✗ gbData->err_fast = fmax(gbData->err_fast, gbData->err[i]);
1687 } else {
1688 ✗ gbData->slowStatesIdx[gbData->nSlowStates] = i;
1689 ✗ gbData->nSlowStates++;
1690 ✗ gbData->err_slow = fmax(gbData->err_slow, gbData->err[i]);
1691 }
1692 }
1693 }
1694
1695 // Reject the current integration step if the estimated error exceeds the tolerance,
1696 // and if the solver is not running with a fixed (constant) step size.
1697 ✗ if (err > 1 && gbData->ctrl_method != GB_CTRL_CNST) {
1698
1699 // Logging
1700 ✗ if (gbData->multi_rate) {
1701 // For multi-rate integration, print detailed info about the rejected step,
1702 // including the slow states' error and the reduced step size.
1703 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0,
1704 "Reject step from %.16g to %.16g, error slow states %.16g, new stepsize %.16g",
1705 ✗ gbData->time, gbData->time + gbData->stepSize, err, gbData->stepSize * 0.5);
1706
1707 // If verbose solver logging is enabled, print detailed error information for debugging.
1708 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER_V)) {
1709 ✗ infoStreamPrint(OMC_LOG_SOLVER_V, 1, "Error of the states: threshold = %15.10g", err_states);
1710 ✗ printVector_gb(OMC_LOG_SOLVER_V, "y", gbData->y, nStates, gbData->time + gbData->stepSize);
1711 ✗ printVector_gb(OMC_LOG_SOLVER_V, "er", gbData->err, nStates, gbData->time + gbData->stepSize);
1712 ✗ messageClose(OMC_LOG_SOLVER_V);
1713 }
1714
1715 // If GBODE state logging is active, dump fast state data to file for further analysis.
1716 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1717 ✗ gbData->err_slow = err; // FIXME should this really only happen when logging is active?
1718 ✗ dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->stepSize, 1);
1719 }
1720 } else {
1721 // For single-rate integration, print basic rejection info with the error and new step size.
1722 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0,
1723 "Reject step from %.16g to %.16g, error %.16g, new stepsize %.16g",
1724 ✗ gbData->time, gbData->time + gbData->stepSize, err, gbData->stepSize * 0.5);
1725 }
1726
1727 // Increment the counter for error test failures.
1728 ✗ gbData->stats.nErrorTestFailures++;
1729
1730 ✗ if (gbData->eventHappened)
1731 {
1732 // event or initial step rejection: reduce step size to 10% of previous step size
1733 ✗ gbData->stepSize *= 0.1;
1734 }
1735 else
1736 {
1737 // standard rejection: reduce the step size by half to attempt a more accurate integration in the next iteration.
1738 ✗ gbData->stepSize *= 0.5;
1739 }
1740
1741 // Restart the integration loop with the smaller step size.
1742 ✗ continue;
1743 }
1744
1745 // Store right-hand side values for later interpolation, including event handling:
1746 // - Update gbData->timeRight to the time at the end of the current step.
1747 // - Copy current state values (gbData->y) to gbData->yRight.
1748 //
1749 // Update the derivative estimates gbData->kRight:
1750 // - If the tableau does not provide kRight values directly,
1751 // compute them by evaluating the ODE function at timeRight and current states.
1752 //
1753 // Compute interpolation error estimate (gbData->err_int) if either:
1754 // - Solver logging is enabled, or
1755 // - The control method is not constant step size and
1756 // the interpolation method is one of the error-controlled Hermite or dense output.
1757 //
1758 // For multi-rate integration with fast states, compute the interpolation error only
1759 // for the slow states subset; otherwise, consider all states.
1760 ✗ gbData->timeRight = gbData->time + gbData->stepSize;
1761 ✗ memcpy(gbData->yRight, gbData->y, nStates * sizeof(double));
1762 // update kRight
1763 ✗ if (!gbData->tableau->isKRightAvailable) {
1764 ✗ sData->timeValue = gbData->timeRight;
1765 ✗ memcpy(sData->realVars, gbData->y, data->modelData->nStates * sizeof(double));
1766 ✗ gbode_fODE(data, threadData, &(gbData->stats.nCallsODE), NULL);
1767 ✗ memcpy(gbData->kRight, fODE, nStates * sizeof(double));
1768 }
1769 else
1770 {
1771 // last stage of method already provides the vector (we have some tiny error of the Newton iteration though)
1772 ✗ memcpy(gbData->kRight, &gbData->k[nStates * (nStages - 1)], nStates * sizeof(double));
1773 }
1774
1775 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER) || noConst_intWithErrctrl) {
1776 ✗ if (gbData->multi_rate && gbData->nFastStates>0) {
1777 ✗ gbData->err_int = error_interpolation_gb(gbData, gbData->nSlowStates, gbData->slowStatesIdx, tol);
1778 } else {
1779 ✗ gbData->err_int = error_interpolation_gb(gbData, nStates, NULL, tol);
1780 }
1781 }
1782 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)) {
1783 // debug the changes of the state values during integration
1784 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 1, "Interpolation error of slow states at midpoint:");
1785 ✗ if (gbData->multi_rate) {
1786 ✗ printVector_gbf(OMC_LOG_GBODE_V, "yL", gbData->yLeft, nStates, gbData->timeLeft, gbData->nSlowStates, gbData->slowStatesIdx);
1787 ✗ printVector_gbf(OMC_LOG_GBODE_V, "kL", gbData->kLeft, nStates, gbData->timeLeft, gbData->nSlowStates, gbData->slowStatesIdx);
1788 ✗ printVector_gbf(OMC_LOG_GBODE_V, "yR", gbData->yRight, nStates, gbData->timeRight, gbData->nSlowStates, gbData->slowStatesIdx);
1789 ✗ printVector_gbf(OMC_LOG_GBODE_V, "kR", gbData->kRight, nStates, gbData->timeRight, gbData->nSlowStates, gbData->slowStatesIdx);
1790 ✗ printVector_gbf(OMC_LOG_GBODE_V, "e", gbData->errest, nStates, (gbData->timeLeft + gbData->timeRight)/2, gbData->nSlowStates, gbData->slowStatesIdx);
1791 } else {
1792 ✗ printVector_gb(OMC_LOG_GBODE_V, "yL", gbData->yLeft, nStates, gbData->timeLeft);
1793 ✗ printVector_gb(OMC_LOG_GBODE_V, "yR", gbData->yRight, nStates, gbData->timeRight);
1794 ✗ printVector_gb(OMC_LOG_GBODE_V, "kL", gbData->kLeft, nStates, gbData->timeLeft);
1795 ✗ printVector_gb(OMC_LOG_GBODE_V, "kR", gbData->kRight, nStates, gbData->timeRight);
1796 ✗ printVector_gbf(OMC_LOG_GBODE_V, "e", gbData->errest, nStates, (gbData->timeLeft + gbData->timeRight)/2, gbData->nSlowStates, gbData->slowStatesIdx);
1797 }
1798 ✗ messageClose(OMC_LOG_GBODE_V);
1799 }
1800
1801 /* TODO: do we even need this condition anymore? */
1802
1803 // Adjust the error estimate for step size control by incorporating interpolation error.
1804 // This is done only if:
1805 // - The current error estimate is greater than 0.01,
1806 // - The number of retries is less than 4,
1807 // - The solver is not using a constant step size,
1808 // - And the interpolation method supports error control (Hermite or dense output).
1809 //
1810 // The error used for step size control is set to the maximum of the interpolation error and the current error.
1811 // if ((err > 1e-2) && (retries < 4) && noConst_intWithErrctrl && gbData->multi_rate) {
1812 // err = fmax(gbData->err_int, err);
1813 // }
1814
1815 // Reject the current integration step if the interpolation error exceeds the tolerance,
1816 // provided that the solver is not running with a fixed step size and interpolation error control is enabled.
1817 //
1818 // On rejection:
1819 // - Increment the retry counter and error test failure statistics.
1820 // - Reduce the step size by half to attempt a more accurate integration.
1821 // - Abort the simulation if the step size falls below the minimal allowed threshold.
1822 //
1823 // Logging differs for multi-rate and single-rate integration:
1824 // - For multi-rate, log errors of slow states and interpolation error.
1825 // - For single-rate, log the overall error and interpolation error.
1826 //
1827 // If multi-rate integration and GBODE state logging is active, dump fast states for diagnostics.
1828 //
1829 // If the step is accepted, reset the retry counter.
1830 ✗ if (err > 1 && noConst_intWithErrctrl) {
1831
1832 retries++;
1833 ✗ gbData->stats.nErrorTestFailures++;
1834 ✗ gbData->stepSize *= 0.5;
1835
1836 ✗ if (gbData->stepSize < GB_MINIMAL_STEP_SIZE) {
1837 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
1838 "Simulation aborted! Minimum step size %g reached, but interpolation error still too large.",
1839 GB_MINIMAL_STEP_SIZE);
1840 ✗ messageClose(OMC_LOG_SOLVER);
1841 ✗ return -1;
1842 }
1843
1844 ✗ if (gbData->multi_rate) {
1845 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0,
1846 "Reject step from %.16g to %.16g, error slow states %.16g, error interpolation %.16g, new stepsize %.16g",
1847 ✗ gbData->time, gbData->time + gbData->stepSize, gbData->err_slow, gbData->err_int, gbData->stepSize);
1848 } else {
1849 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0,
1850 "Reject step from %.16g to %.16g, error %.16g, interpolation error %.16g, new stepsize %.16g",
1851 ✗ gbData->time, gbData->time + gbData->stepSize, err_states, gbData->err_int, gbData->stepSize);
1852 }
1853
1854 ✗ if (gbData->multi_rate && OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1855 // Dump fast states to file for further analysis after step rejection.
1856 ✗ dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->stepSize, 2);
1857 }
1858
1859 ✗ continue;
1860 } else {
1861 // Reset retries counter if the step was accepted.
1862 retries = 0;
1863 }
1864
1865 /* Step is accepted from here on, as err <= 1 */
1866
1867 /* remember last time values for dense output extrapolation with yLast, kLast */
1868 ✗ gbData->extrapolationBaseTime = gbData->time;
1869 ✗ gbData->extrapolationStepSize = gbData->stepSize;
1870 ✗ gbData->eventHappened = FALSE;
1871 ✗ gbData->didFastStep = FALSE;
1872
1873 /* remember kLast and yLast for dense output extrapolation */
1874 ✗ memcpy(gbData->kLast, gbData->k, nStates * nStages * sizeof(double));
1875 ✗ memcpy(gbData->yLast, gbData->yOld, nStates * sizeof(double));
1876
1877 // Rotate the error and step size ring buffers to make room for the latest values.
1878 // The oldest entries are shifted one position towards the end,
1879 // and the newest error and step size values are stored at the front (index 0).
1880 // FIXME use actual ring buffer instead of moving data around!
1881 ✗ for (i = (gbData->ringBufferSize - 1); i > 0; i--) {
1882 ✗ gbData->errValues[i] = gbData->errValues[i - 1];
1883 ✗ gbData->stepSizeValues[i] = gbData->stepSizeValues[i - 1];
1884 }
1885 // Store the current error and step size at the beginning of the buffers.
1886 ✗ gbData->errValues[0] = err;
1887 ✗ gbData->stepSizeValues[0] = gbData->stepSize;
1888
1889 // Update the step size using the step size controller
1890 ✗ gbData->lastStepSize = gbData->stepSize; // Save the current step size before updating
1891 // Calculate a new step size based on recent error and step size history,
1892 // the method’s error order, and the control method in use
1893 ✗ gbData->stepSize *= GenericController(gbData->errValues, gbData->stepSizeValues, gbData->currentErrorOrder, gbData->ctrl_method);
1894
1895 // Ensure the new step size does not exceed the user-defined maximum step size (if set)
1896 ✗ if (gbData->maxStepSize > 0 && gbData->maxStepSize < gbData->stepSize)
1897 ✗ gbData->stepSize = gbData->maxStepSize;
1898
1899 // Store the optimized step size for further use
1900 ✗ gbData->optStepSize = gbData->stepSize;
1901
1902 ✗ if (gbData->multi_rate) {
1903 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1904 // debug the changes of the state values during integration
1905 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "States and derivatives at right hand side:");
1906 ✗ printVector_gb(OMC_LOG_GBODE, "yR", gbData->yRight, nStates, gbData->timeRight);
1907 ✗ printVector_gb(OMC_LOG_GBODE, "kR", gbData->kRight, nStates, gbData->timeRight);
1908 ✗ messageClose(OMC_LOG_GBODE);
1909 }
1910
1911 ✗ if (gbData->nFastStates > 0) {
1912 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1913 // debug the error of the states and derivatives after outer integration
1914 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Error of the states before inner integration: threshold = %15.10g", err_states);
1915 ✗ printVector_gb(OMC_LOG_GBODE, "er", gbData->err, nStates, gbData->timeRight);
1916 ✗ messageClose(OMC_LOG_GBODE);
1917 }
1918 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
1919 // dump fast states in file
1920 ✗ dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->lastStepSize, -1);
1921 }
1922 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Refine step from %10g to %10g, error fast states %10g, error interpolation %10g, new stepsize %10g",
1923 ✗ gbData->time, gbData->time + gbData->lastStepSize, gbData->err_fast, error_interpolation_gb(gbData, nStates, NULL, tol), gbData->stepSize);
1924 // run multirate step
1925 ✗ gb_step_info = gbodef_main(data, threadData, solverInfo, targetTime);
1926 // synchronize relevant information
1927 ✗ if (fabs(gbData->timeRight - gbData->gbfData->timeRight) < GB_MINIMAL_STEP_SIZE) {
1928 ✗ memcpy(gbData->y, gbData->gbfData->y, nStates * sizeof(double));
1929 ✗ memcpy(gbData->yRight, gbData->gbfData->yRight, nStates * sizeof(double));
1930 ✗ memcpy(gbData->err, gbData->gbfData->err, nStates * sizeof(double));
1931 ✗ sData->timeValue = gbData->timeRight;
1932 ✗ memcpy(sData->realVars, gbData->yRight, data->modelData->nStates * sizeof(double));
1933 ✗ gbode_fODE(data, threadData, &(gbData->stats.nCallsODE), NULL);
1934 ✗ memcpy(gbData->kRight, fODE, nStates * sizeof(double));
1935 }
1936 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Refined step from %10g to %10g, error fast states %10g, error interpolation %10g, new stepsize %10g",
1937 ✗ gbData->time, gbData->time + gbData->lastStepSize, gbData->err_fast, error_interpolation_gb(gbData, nStates, NULL, tol), gbData->stepSize);
1938 ✗ if (gb_step_info !=0) {
1939 // get out of here, if an event has happend!
1940 ✗ messageClose(OMC_LOG_SOLVER);
1941 ✗ if (gb_step_info>0)
1942 return 0;
1943 else
1944 ✗ return gb_step_info;
1945 }
1946 }
1947
1948 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)) {
1949 // debug the error of the states and derivatives after outer integration
1950 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 1, "Error of the states: threshold = %15.10g", err_states);
1951 ✗ printVector_gb(OMC_LOG_GBODE_V, "er", gbData->err, nStates, gbData->timeRight);
1952 ✗ messageClose(OMC_LOG_GBODE_V);
1953 }
1954 }
1955 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)) {
1956 // debug ring buffer for the states and derviatives of the states
1957 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 1, "Ring buffer during steps of integration");
1958 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 0, "Old value:");
1959 ✗ printVector_gb(OMC_LOG_GBODE_V, "y", gbData->yOld, nStates, gbData->time);
1960 ✗ debugRingBuffer_gb(OMC_LOG_GBODE_V, gbData->x, gbData->k, nStates, gbData->tableau, gbData->time, gbData->lastStepSize);
1961 ✗ infoStreamPrint(OMC_LOG_GBODE_V, 0, "New value:");
1962 ✗ printVector_gb(OMC_LOG_GBODE_V, "y", gbData->y, nStates, gbData->time + gbData->lastStepSize);
1963 ✗ messageClose(OMC_LOG_GBODE_V);
1964 }
1965 ✗ } while (!isfinite(err) || (err > 1 && gbData->ctrl_method != GB_CTRL_CNST));
1966
1967 // count processed steps
1968 ✗ gbData->stats.nStepsTaken++;
1969
1970 // debug the changes of the state values during integration
1971 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1972 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "States and derivatives at right hand side:");
1973 ✗ printVector_gb(OMC_LOG_GBODE, "yR", gbData->yRight, nStates, gbData->timeRight);
1974 ✗ printVector_gb(OMC_LOG_GBODE, "kR", gbData->kRight, nStates, gbData->timeRight);
1975 ✗ messageClose(OMC_LOG_GBODE);
1976 }
1977
1978 // If not using multi-rate integration, or if multi-rate is active but the fast integration time
1979 // is behind the main integrator time, then check for events.
1980 ✗ if (!gbData->multi_rate || (gbData->multi_rate && gbData->gbfData->time < gbData->time)) {
1981
1982 // Check for any events occurring between the previous accepted time (timeLeft) and current time (timeRight).
1983 // The function returns the event time if an event is detected, and sets foundEvent accordingly.
1984 ✗ foundEvent = checkForEvents(data, threadData, solverInfo, gbData->timeLeft, gbData->yLeft, gbData->timeRight, gbData->yRight, FALSE, &(gbData->eventTime));
1985
1986 ✗ if (foundEvent) {
1987 // Clear any pending events in the solver's event list before handling the new event.
1988 ✗ listClear(solverInfo->eventLst);
1989
1990 // Update the current integration time to the event time.
1991 ✗ gbData->time = gbData->eventTime;
1992 ✗ gbData->eventHappened = TRUE;
1993
1994 // Perform interpolation at the event time to estimate states and derivatives accurately.
1995 ✗ gb_interpolation(gbData->interpolation,
1996 gbData->timeLeft, gbData->yLeft, gbData->kLeft,
1997 gbData->timeRight, gbData->yRight, gbData->kRight,
1998 gbData->time, gbData->yOld,
1999 nStates, NULL, nStates, gbData->tableau,
2000 gbData->x, gbData->k);
2001
2002 // Adjust targetTime to not exceed the detected event time,
2003 // ensuring the integrator stops exactly at the event.
2004 ✗ targetTime = fmin(targetTime, gbData->eventTime);
2005
2006 // Exit the integration loop early since an event was detected.
2007 ✗ break;
2008 }
2009 }
2010
2011 ✗ if (gbData->multi_rate) {
2012 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Accept step from %.16g to %.16g, error slow states %.16g, error interpolation %.16g, new stepsize %.16g",
2013 gbData->timeLeft, gbData->timeRight, err_states, gbData->err_int, gbData->stepSize);
2014 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_STATES)) {
2015 // dump fast states in file
2016 ✗ dumpFastStates_gb(gbData, FALSE, gbData->time, 0);
2017 }
2018 } else {
2019 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Accept step from %.16g to %.16g, error %.16g interpolation error %.16g, new stepsize %16g",
2020 gbData->timeLeft, gbData->timeRight, err_states, gbData->err_int, gbData->stepSize);
2021
2022 }
2023
2024 /* update time with performed stepSize */
2025 ✗ gbData->time = gbData->timeRight;
2026
2027 /* step is accepted and yOld needs to be updated */
2028 ✗ memcpy(gbData->yOld, gbData->yRight, nStates * sizeof(double));
2029
2030 // Rotate ring buffer
2031 ✗ for (i = (gbData->ringBufferSize - 1); i > 0 ; i--) {
2032 ✗ gbData->tv[i] = gbData->tv[i - 1];
2033 ✗ memcpy(gbData->yv + i * nStates, gbData->yv + (i - 1) * nStates, nStates * sizeof(double));
2034 ✗ memcpy(gbData->kv + i * nStates, gbData->kv + (i - 1) * nStates, nStates * sizeof(double));
2035 }
2036
2037 // update new values
2038 ✗ gbData->tv[0] = gbData->timeRight;
2039 ✗ memcpy(gbData->yv, gbData->yRight, nStates * sizeof(double));
2040 ✗ memcpy(gbData->kv, gbData->kRight, nStates * sizeof(double));
2041
2042 ✗ debugRingBufferSteps_gb(OMC_LOG_GBODE_V, gbData->yv, gbData->kv, gbData->tv, nStates, gbData->ringBufferSize);
2043
2044 /* emit step, if solverNoEquidistantGrid is selected */
2045 ✗ if (solverInfo->solverNoEquidistantGrid && (!gbData->multi_rate || (gbData->multi_rate && gbData->gbfData->time<gbData->time))) {
2046 ✗ sData->timeValue = gbData->time;
2047 ✗ solverInfo->currentTime = sData->timeValue;
2048 ✗ memcpy(sData->realVars, gbData->y, nStates * sizeof(double));
2049 // log the emitted result
2050 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)){
2051 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Emit result:");
2052 ✗ printVector_gb(OMC_LOG_GBODE, " y", sData->realVars, nStates, sData->timeValue);
2053 ✗ messageClose(OMC_LOG_GBODE);
2054 }
2055 break;
2056 }
2057
2058 // stop, if simulation nearly reached stopTime
2059 ✗ if (stopTime - gbData->time < GB_MINIMAL_STEP_SIZE) {
2060 ✗ gbData->time = stopTime;
2061 ✗ break;
2062 }
2063 } // end of while-loop (gbData->time < targetTime)
2064
2065 ✗ if (gbData->eventTime == targetTime) {
2066
2067 ✗ if (!solverInfo->solverNoEquidistantGrid) {
2068 ✗ foundEvent = checkForEvents(data, threadData, solverInfo, gbData->eventTime, gbData->yOld, gbData->eventTime, gbData->yOld, FALSE, &(gbData->eventTime));
2069 }
2070
2071 ✗ solverInfo->currentTime = gbData->time;
2072 ✗ sData->timeValue = gbData->time;
2073 ✗ memcpy(sData->realVars, gbData->yOld, nStates * sizeof(double));
2074
2075 // if noRestart is set, the right hand side values are stored
2076 ✗ if (gbData->noRestart) {
2077 ✗ gbData->timeRight = gbData->time;
2078 ✗ memcpy(gbData->yRight, gbData->yOld, nStates * sizeof(double));
2079 ✗ gbode_fODE(data, threadData, &(gbData->stats.nCallsODE), NULL);
2080 ✗ memcpy(gbData->kRight, fODE, nStates * sizeof(double));
2081 }
2082
2083 /* write statistics to the solverInfo data structure */
2084 ✗ memcpy(&solverInfo->solverStatsTmp, &gbData->stats, sizeof(SOLVERSTATS));
2085
2086 // log the emitted result
2087 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)){
2088 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Emit result (single-rate integration):");
2089 ✗ printVector_gb(OMC_LOG_GBODE, " y", sData->realVars, nStates, sData->timeValue);
2090 ✗ messageClose(OMC_LOG_GBODE);
2091 }
2092
2093 ✗ listClear(solverInfo->eventLst);
2094 ✗ gbData->eventTime = DBL_MAX; // reset event time, if eventTime is reached
2095
2096 // return to solver main routine for proper event handling (iteration)
2097 ✗ messageClose(OMC_LOG_SOLVER);
2098
2099 ✗ return 0;
2100 }
2101
2102 ✗ if (!solverInfo->solverNoEquidistantGrid) {
2103 /* Integrator does large steps and needs to interpolate results with respect to the output grid */
2104 ✗ sData->timeValue = fmin(solverInfo->currentTime + solverInfo->currentStepSize, targetTime);
2105 ✗ sData->timeValue = fmin(sData->timeValue, stopTime);
2106 ✗ solverInfo->currentTime = sData->timeValue;
2107
2108 ✗ if (gbData->multi_rate) {
2109 // if the inner integration has not been started, the outer values need to be emitted
2110 ✗ if (gbData->gbfData->time >= sData->timeValue) {
2111 ✗ gb_interpolation(gbData->interpolation,
2112 gbData->timeLeft, gbData->yLeft, gbData->kLeft,
2113 gbData->timeRight, gbData->yRight, gbData->kRight,
2114 ✗ sData->timeValue, sData->realVars,
2115 gbData->nSlowStates, gbData->slowStatesIdx, nStates, gbData->tableau, gbData->x, gbData->k);
2116
2117 ✗ DATA_GBODEF *gbfData = gbData->gbfData;
2118 ✗ gb_interpolation(gbData->gbfData->interpolation,
2119 gbfData->timeLeft, gbfData->yLeft, gbfData->kLeft,
2120 gbfData->timeRight, gbfData->yRight, gbfData->kRight,
2121 ✗ sData->timeValue, sData->realVars,
2122 gbfData->nFastStates, gbData->fastStatesIdx, nStates, gbfData->tableau, gbfData->x, gbfData->k);
2123 } else {
2124 ✗ gb_interpolation(gbData->interpolation,
2125 gbData->timeLeft, gbData->yLeft, gbData->kLeft,
2126 gbData->timeRight, gbData->yRight, gbData->kRight,
2127 ✗ sData->timeValue, sData->realVars,
2128 nStates, NULL, nStates, gbData->tableau, gbData->x, gbData->k);
2129 }
2130 } else {
2131 // use chosen interpolation for emitting equidistant output (default hermite)
2132 ✗ if (solverInfo->currentStepSize>0)
2133 ✗ gb_interpolation(gbData->interpolation,
2134 gbData->timeLeft, gbData->yLeft, gbData->kLeft,
2135 gbData->timeRight, gbData->yRight, gbData->kRight,
2136 ✗ sData->timeValue, sData->realVars,
2137 nStates, NULL, nStates, gbData->tableau, gbData->x, gbData->k);
2138 }
2139 // log the emitted result
2140 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)){
2141 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Emit result:");
2142 ✗ printVector_gb(OMC_LOG_GBODE, " y", sData->realVars, nStates, sData->timeValue);
2143 ✗ messageClose(OMC_LOG_GBODE);
2144 }
2145 } else {
2146 // Integrator emits result on the simulation grid (see above)
2147 ✗ sData->timeValue = gbData->time;
2148 ✗ solverInfo->currentTime = gbData->time;
2149 ✗ solverInfo->currentStepSize = gbData->stepSize;
2150 }
2151
2152 /* if a state event occurs than no sample event does need to be activated */
2153 ✗ data->simulationInfo->sampleActivated = data->simulationInfo->sampleActivated
2154 ✗ && solverInfo->currentTime >= data->simulationInfo->nextSampleEvent;
2155
2156 /* Solver statistics */
2157 ✗ if (!gbData->isExplicit && gbData->nlsSolverMethod != GB_NLS_INTERNAL)
2158 ✗ gbData->stats.nCallsJacobian = gbData->nlsData->numberOfJEval;
2159 ✗ if (!solverInfo->solverNoEquidistantGrid && fabs(targetTime - stopTime) < GB_MINIMAL_STEP_SIZE && OMC_ACTIVE_STREAM(OMC_LOG_STATS)) {
2160 ✗ if (gbData->multi_rate) {
2161 ✗ infoStreamPrint(OMC_LOG_STATS, 0, "gbode (birate integration): slow: %s / fast: %s",
2162 ✗ GB_METHOD_NAME[gbData->GM_method], GB_METHOD_NAME[gbData->gbfData->GM_method]);
2163 ✗ logSolverStats(OMC_LOG_STATS, "inner integration", stopTime, stopTime, 0, &gbData->gbfData->stats, &gbData->gbfData->fastStateUpdateCount, &gbData->gbfData->additionalFullODEEvaluations);
2164 ✗ logSolverStats(OMC_LOG_STATS, "outer integration", stopTime, stopTime, 0, &gbData->stats, NULL, NULL);
2165 } else {
2166 ✗ infoStreamPrint(OMC_LOG_STATS, 0, "gbode (single-rate integration): %s", GB_METHOD_NAME[gbData->GM_method]);
2167 }
2168 }
2169 /* Write statistics to the solverInfo data structure */
2170 ✗ logSolverStats(OMC_LOG_SOLVER_V, "gb_singlerate", solverInfo->currentTime, gbData->time, gbData->stepSize, &gbData->stats, NULL, NULL);
2171 ✗ memcpy(&solverInfo->solverStatsTmp, &gbData->stats, sizeof(SOLVERSTATS));
2172
2173 ✗ messageClose(OMC_LOG_SOLVER);
2174
2175 ✗ return 0;
2176 }
2177