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 / 495
Functions: 0.0% 0 / 0 / 8
Branches: 0.0% 0 / 0 / 312
OMCompiler/SimulationRuntime/c/simulation/solver/gbode_step.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_step.c
29 */
30
31 #include "gbode_main.h"
32 #include "gbode_nls.h"
33 #include "gbode_internal_nls.h"
34 #include "gbode_err.h"
35 #include "gbode_util.h"
36
37 #include "kinsolSolver.h"
38
39 #include <math.h>
40
41 /**
42 * @brief Generic multi-step function.
43 *
44 * Internal non-linear equation system will be solved with non-linear solver specified during setup.
45 * Results will be saved in y and the signed error estimate in yt.
46 *
47 * @param data Runtime data struct.
48 * @param threadData Thread data for error handling.
49 * @param solverInfo Storing Runge-Kutta solver data.
50 * @return int Return 0 on success, -1 on failure.
51 */
52 ✗ int full_implicit_MS(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
53 {
54 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
55 ✗ modelica_real* fODE = sData->realVars + data->modelData->nStates;
56 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
57
58 int i;
59 int stage;
60 ✗ int nStates = data->modelData->nStates;
61 ✗ int nStages = gbData->tableau->nStages;
62 NLS_SOLVER_STATUS solved = NLS_FAILED;
63
64 /* Predictor step */
65 ✗ for (i = 0; i < nStates; i++) {
66 ✗ gbData->yt[i] = 0;
67 ✗ for (stage = 0; stage < nStages-1; stage++) {
68 ✗ gbData->yt[i] += -gbData->yv[stage * nStates + i] * gbData->tableau->c[stage] +
69 ✗ gbData->kv[stage * nStates + i] * gbData->tableau->bt[stage] * gbData->stepSize;
70 }
71 ✗ gbData->yt[i] += gbData->kv[stage * nStates + i] * gbData->tableau->bt[stage] * gbData->stepSize;
72 ✗ gbData->yt[i] /= gbData->tableau->c[stage];
73 }
74
75
76 /* Constant part of the multi-step method */
77 ✗ for (i = 0; i < nStates; i++) {
78 ✗ gbData->res_const[i] = 0;
79 ✗ for (stage = 0; stage < nStages-1; stage++) {
80 ✗ gbData->res_const[i] += -gbData->yv[stage * nStates + i] * gbData->tableau->c[stage] +
81 ✗ gbData->kv[stage * nStates + i] * gbData->tableau->b[stage] * gbData->stepSize;
82 }
83 }
84 // printVector_gb("res_const: ", gbData->res_const, nStates, gbData->time);
85
86 /* Compute intermediate step k, explicit if diagonal element is zero, implicit otherwise
87 * k[i] = f(tOld + c[i]*h, yOld + h*sum(A[i,j]*k[j], i=j..i)) */
88 // here, it yields: stage == stage_, and stage * nStages + stage_ is index of the diagonal element
89
90 // set simulation time with respect to the current stage
91 ✗ sData->timeValue = gbData->time + gbData->stepSize;
92
93 // solve for x: 0 = yold-x + h*(sum(A[i,j]*k[j], i=1..j-1) + A[i,i]*f(t + c[i]*h, x))
94 ✗ NONLINEAR_SYSTEM_DATA* nlsData = gbData->nlsData;
95
96 // Set start vectors for the noblinear solver
97 ✗ memcpy(nlsData->nlsx, gbData->yt, nStates*sizeof(modelica_real));
98 ✗ memcpy(nlsData->nlsxOld, nlsData->nlsx, nStates*sizeof(modelica_real));
99 ✗ memcpy(nlsData->nlsxExtrapolation, nlsData->nlsx, nStates*sizeof(modelica_real));
100
101 ✗ solved = solveNLS_gb(data, threadData, nlsData, gbData, FALSE);
102
103 ✗ if (solved != NLS_SOLVED) {
104 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "gbode error: Failed to solve NLS in full_implicit_MS at time t=%g", gbData->time);
105 ✗ return -1;
106 }
107
108 ✗ memcpy(gbData->kv + stage * nStates, fODE, nStates*sizeof(double));
109
110 /* Corrector step */
111 ✗ for (i = 0; i < nStates; i++) {
112 ✗ gbData->y[i] = 0;
113 ✗ for (stage = 0; stage < nStages-1; stage++) {
114 ✗ gbData->y[i] += -gbData->yv[stage * nStates + i] * gbData->tableau->c[stage] +
115 ✗ gbData->kv[stage * nStates + i] * gbData->tableau->b[stage] * gbData->stepSize;
116 }
117 ✗ gbData->y[i] += gbData->kv[stage * nStates + i] * gbData->tableau->b[stage] * gbData->stepSize;
118 ✗ gbData->y[i] /= gbData->tableau->c[stage];
119 ✗ gbData->yt[i] = gbData->y[i] - gbData->yt[i];
120 }
121
122 return 0;
123 }
124
125 /**
126 * @brief Generic multi-step function.
127 *
128 * Internal non-linear equation system will be solved with non-linear solver specified during setup.
129 * Results will be saved in y and the signed error estimate in yt.
130 *
131 * @param data Runtime data struct.
132 * @param threadData Thread data for error handling.
133 * @param solverInfo Storing Runge-Kutta solver data.
134 * @return int Return 0 on success, -1 on failure.
135 */
136 ✗ int full_implicit_MS_MR(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
137 {
138 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
139 ✗ modelica_real* fODE = sData->realVars + data->modelData->nStates;
140 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
141 ✗ DATA_GBODEF* gbfData = gbData->gbfData;
142
143 int i, ii;
144 int stage;
145 ✗ int nStates = data->modelData->nStates;
146 ✗ int nStages = gbfData->tableau->nStages;
147 NLS_SOLVER_STATUS solved = NLS_FAILED;
148
149 /* Predictor step */
150 ✗ for (ii = 0; ii < gbData->nFastStates; ii++)
151 {
152 ✗ i = gbData->fastStatesIdx[ii];
153 ✗ gbfData->yt[i] = 0;
154 ✗ for (stage = 0; stage < nStages-1; stage++)
155 {
156 ✗ gbfData->yt[i] += -gbfData->yv[stage * nStates + i] * gbfData->tableau->c[stage] +
157 ✗ gbfData->kv[stage * nStates + i] * gbfData->tableau->bt[stage] * gbfData->stepSize;
158 }
159 ✗ gbfData->yt[i] += gbfData->kv[stage * nStates + i] * gbfData->tableau->bt[stage] * gbfData->stepSize;
160 ✗ gbfData->yt[i] /= gbfData->tableau->c[stage];
161 }
162
163
164 /* Constant part of the multi-step method */
165 ✗ for (ii = 0; ii < gbData->nFastStates; ii++)
166 {
167 ✗ i = gbData->fastStatesIdx[ii];
168 ✗ gbfData->res_const[i] = 0;
169 ✗ for (stage = 0; stage < nStages-1; stage++)
170 {
171 ✗ gbfData->res_const[i] += -gbfData->yv[stage * nStates + i] * gbfData->tableau->c[stage] +
172 ✗ gbfData->kv[stage * nStates + i] * gbfData->tableau->b[stage] * gbfData->stepSize;
173 }
174 }
175 // printVector_gb("res_const: ", gbData->res_const, nStates, gbData->time);
176
177 /* Compute intermediate step k, explicit if diagonal element is zero, implicit otherwise
178 * k[i] = f(tOld + c[i]*h, yOld + h*sum(A[i,j]*k[j], i=j..i)) */
179 // here, it yields: stage == stage_, and stage * nStages + stage_ is index of the diagonal element
180
181 // set simulation time with respect to the current stage
182 ✗ sData->timeValue = gbfData->time + gbfData->stepSize;
183 // interpolate the slow states on the current time of gbfData->yOld for correct evaluation of gbfData->res_const
184 ✗ gb_interpolation(gbData->interpolation,
185 gbData->timeLeft, gbData->yLeft, gbData->kLeft,
186 gbData->timeRight, gbData->yRight, gbData->kRight,
187 sData->timeValue, sData->realVars,
188 gbData->nSlowStates, gbData->slowStatesIdx, nStates, gbData->tableau, gbData->x, gbData->k);
189
190 // solve for x: 0 = yold-x + h*(sum(A[i,j]*k[j], i=1..j-1) + A[i,i]*f(t + c[i]*h, x))
191 ✗ NONLINEAR_SYSTEM_DATA* nlsData = gbfData->nlsData;
192
193 ✗ projVector_gbf(nlsData->nlsx, gbfData->yt, gbData->nFastStates, gbData->fastStatesIdx);
194 ✗ memcpy(nlsData->nlsxOld, nlsData->nlsx, nStates*sizeof(modelica_real));
195 ✗ memcpy(nlsData->nlsxExtrapolation, nlsData->nlsx, nStates*sizeof(modelica_real));
196
197 ✗ solved = solveNLS_gb(data, threadData, nlsData, gbData, TRUE);
198
199 ✗ if (solved != NLS_SOLVED) {
200 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "gbodef error: Failed to solve NLS in full_implicit_MS_MR at time t=%g", gbfData->time);
201 ✗ return -1;
202 }
203
204 ✗ memcpy(gbfData->kv + stage * nStates, fODE, nStates*sizeof(double));
205
206 /* Corrector step */
207 ✗ for (ii = 0; ii < gbData->nFastStates; ii++)
208 {
209 ✗ i = gbData->fastStatesIdx[ii];
210 ✗ gbfData->y[i] = 0;
211 ✗ for (stage = 0; stage < nStages-1; stage++)
212 {
213 ✗ gbfData->y[i] += -gbfData->yv[stage * nStates + i] * gbfData->tableau->c[stage] +
214 ✗ gbfData->kv[stage * nStates + i] * gbfData->tableau->b[stage] * gbfData->stepSize;
215 }
216 ✗ gbfData->y[i] += gbfData->kv[stage * nStates + i] * gbfData->tableau->b[stage] * gbfData->stepSize;
217 ✗ gbfData->y[i] /= gbfData->tableau->c[stage];
218 ✗ gbfData->yt[i] = gbfData->y[i] - gbfData->yt[i];
219 }
220
221 return 0;
222 }
223
224 /**
225 * @brief Generic diagonal implicit Runge-Kutta step function.
226 *
227 * Internal non-linear equation system will be solved with non-linear solver specified during setup.
228 * Results are saved in y. The selected error estimator writes |error| to errest.
229 *
230 * @param data Runtime data struct.
231 * @param threadData Thread data for error handling.
232 * @param solverInfo Storing Runge-Kutta solver data.
233 * @return int Return 0 on success, -1 on failure.
234 */
235 ✗ int expl_diag_impl_RK(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
236 {
237 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
238 ✗ modelica_real* fODE = sData->realVars + data->modelData->nStates;
239 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
240
241 int i;
242 int stage, stage_;
243 ✗ int nStates = data->modelData->nStates;
244 ✗ int nStages = gbData->tableau->nStages;
245 NLS_SOLVER_STATUS solved = NLS_FAILED;
246 ✗ GB_ERROR_CONTEXT error_context = {data, threadData, gbData, NULL, FALSE};
247
248 ✗ if (!gbData->isExplicit && OMC_ACTIVE_STREAM(OMC_LOG_GBODE_NLS_V)) {
249 // NLS - used values for extrapolation
250 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS_V, 1, "NLS - used values for extrapolation:");
251 ✗ printVector_gb(OMC_LOG_GBODE_NLS_V, "xL", gbData->yv + nStates, nStates, gbData->tv[1]);
252 ✗ printVector_gb(OMC_LOG_GBODE_NLS_V, "kL", gbData->kv + nStates, nStates, gbData->tv[1]);
253 ✗ printVector_gb(OMC_LOG_GBODE_NLS_V, "xR", gbData->yv, nStates, gbData->tv[0]);
254 ✗ printVector_gb(OMC_LOG_GBODE_NLS_V, "kR", gbData->kv, nStates, gbData->tv[0]);
255 ✗ messageClose(OMC_LOG_GBODE_NLS_V);
256 }
257
258 /* Runge-Kutta step */
259 ✗ for (stage = 0; stage < nStages; stage++)
260 {
261 ✗ gbData->act_stage = stage;
262
263 /* Set constant part or residual input
264 * res = f(tOld + c[i]*h, yOld + h*sum(A[i,j]*k[j], i=j..stage-1)) */
265 ✗ for (i = 0; i < nStates; i++)
266 {
267 ✗ gbData->res_const[i] = gbData->yOld[i];
268 ✗ for (stage_ = 0; stage_ < stage; stage_++)
269 {
270 ✗ gbData->res_const[i] += gbData->stepSize * gbData->tableau->A[stage * nStages + stage_] * gbData->k[stage_ * nStates + i];
271 }
272 }
273
274 /* Compute intermediate step k, explicit if diagonal element is zero, implicit otherwise
275 * k[i] = f(tOld + c[i]*h, yOld + h*sum(A[i,j]*k[j], i=j..i)) */
276 // here, it yields: stage == stage_, and stage * nStages + stage_ is index of the diagonal element
277
278 // set simulation time with respect to the current stage
279 ✗ sData->timeValue = gbData->time + gbData->tableau->c[stage_]*gbData->stepSize;
280
281 // if the diagonal element is zero, an explicit step has to be performed
282 ✗ if (gbData->tableau->A[stage * nStages + stage_] == 0) {
283 // Store values in the ring buffer
284 ✗ memcpy(gbData->x + stage_ * nStates, gbData->res_const, nStates*sizeof(double));
285
286 ✗ if (gbData->tableau->isKLeftAvailable && !gbData->didFastStep && (stage == 0)) {
287 ✗ memcpy(fODE, gbData->kLeft, nStates*sizeof(double));
288 } else {
289 ✗ memcpy(sData->realVars, gbData->res_const, nStates*sizeof(double));
290 ✗ gbode_fODE(data, threadData, &(gbData->stats.nCallsODE), NULL);
291 }
292 } else {
293 // solve for x: 0 = yold-x + h*(sum(A[i,j]*k[j], i=1..j-1) + A[i,i]*f(t + c[i]*h, x))
294 ✗ NONLINEAR_SYSTEM_DATA* nlsData = gbData->nlsData;
295 struct dataSolver * solverData = (struct dataSolver *)nlsData->solverData;
296 NLS_KINSOL_DATA* kin_mem = ((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory;
297
298 // Set start vector
299 ✗ memcpy(nlsData->nlsx, gbData->yOld, nStates*sizeof(modelica_real));
300 ✗ memcpy(nlsData->nlsxExtrapolation, gbData->yOld, nStates*sizeof(modelica_real));
301
302 // is the last solution valid and do we use internal nls
303 ✗ modelica_boolean dense_output_valid = (gbData->time != data->simulationInfo->startTime && !gbData->eventHappened
304 ✗ && gbData->nlsSolverMethod == GB_NLS_INTERNAL && gbData->extrapolationBaseTime != INFINITY);
305
306 // for MR integration: start values of fast states are chosen as left boundary y0; avoids poor extrapolation
307 ✗ modelica_boolean do_zero_order_hold_fast_states = (gbData->multi_rate && gbData->nFastStates > 0);
308
309 ✗ if (gbData->tableau->svp != NULL && gbData->tableau->svp->type[stage_] == SVP_LINEAR_COMBINATION)
310 {
311 /* linear combination stage-value-predictors (highest priority) */
312 ✗ gbInternalLinearCombinationSVP(gbData->tableau->svp, stage_, nStates, gbData->stepSize, gbData->k, gbData->yOld, nlsData->nlsxOld);
313
314 // never do 0 order hold if we do sophisticated SVPs
315 do_zero_order_hold_fast_states = FALSE;
316 }
317 ✗ else if (dense_output_valid && gbData->tableau->svp != NULL && gbData->tableau->svp->type[stage_] == SVP_DENSE_OUTPUT)
318 ✗ {
319 /* dense output stage-value-predictor */
320 ✗ double theta = (gbData->time + gbData->tableau->c[stage_] * gbData->stepSize - gbData->extrapolationBaseTime) / gbData->extrapolationStepSize;
321 ✗ gbData->tableau->svp->dense_output_predictor(gbData->tableau, gbData->yLast, NULL, gbData->kLast,
322 ✗ theta, gbData->extrapolationStepSize, nlsData->nlsxOld, 0, NULL, nStates);
323 }
324 ✗ else if (dense_output_valid && gbData->tableau->withDenseOutput)
325 ✗ {
326 /* standard dense output if available / possible */
327 ✗ double theta = (gbData->time + gbData->tableau->c[stage_] * gbData->stepSize - gbData->extrapolationBaseTime) / gbData->extrapolationStepSize;
328 ✗ gbData->tableau->dense_output(gbData->tableau, gbData->yLast, NULL, gbData->kLast,
329 ✗ theta, gbData->extrapolationStepSize, nlsData->nlsxOld, 0, NULL, nStates);
330 }
331 ✗ else if (stage>1)
332 {
333 /* perform hermite to interpolate between two stages */
334 ✗ extrapolation_hermite_gb(nlsData->nlsxOld, gbData->nStates, gbData->time + gbData->tableau->c[stage_-2] * gbData->stepSize, gbData->x + (stage_-2) * nStates, gbData->k + (stage_-2) * nStates,
335 ✗ gbData->time + gbData->tableau->c[stage_-1] * gbData->stepSize, gbData->x + (stage_-1) * nStates, gbData->k + (stage_-1) * nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
336 }
337 else
338 {
339 /* generic extrapolation */
340 ✗ extrapolation_gb(gbData, nlsData->nlsxOld, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
341 }
342
343 // zero order hold for all fast states
344 ✗ if (do_zero_order_hold_fast_states)
345 {
346 ✗ for (int fast = 0; fast < gbData->nFastStates; fast++)
347 {
348 ✗ int full = gbData->fastStatesIdx[fast];
349 ✗ nlsData->nlsxOld[full] = gbData->yOld[full];
350 }
351 }
352
353 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS_V, 0, "Solving NLS of stage %d at time %g", stage_+1, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
354 ✗ solved = solveNLS_gb(data, threadData, nlsData, gbData, FALSE);
355
356 ✗ if (solved != NLS_SOLVED) {
357 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "gbode error: Failed to solve NLS in expl_diag_impl_RK in stage %d at time t=%g", stage_+1, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
358 ✗ return -1;
359 }
360
361 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_NLS_V)) {
362 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS_V, 1, "NLS - start values and solution of the NLS:");
363 ✗ printVector_gb(OMC_LOG_GBODE_NLS_V, "x0", nlsData->nlsxOld, nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
364 ✗ printVector_gb(OMC_LOG_GBODE_NLS_V, "xS", nlsData->nlsxExtrapolation, nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
365 ✗ printVector_gb(OMC_LOG_GBODE_NLS_V, "xL", nlsData->nlsx, nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
366 ✗ messageClose(OMC_LOG_GBODE_NLS_V);
367 }
368
369 ✗ memcpy(gbData->x + stage_ * nStates, nlsData->nlsx, nStates*sizeof(double));
370 ✗ if (/* non explicit stage of (E)SDIRK integrator */ (stage_ != 0 || gbData->tableau->A[0] != 0) && gbData->nlsSolverMethod == GB_NLS_INTERNAL)
371 {
372 // reconstruct k_{stage_} from the solution, avoids repeated call to functionODE()
373 ✗ double ifac = 1.0 / (gbData->stepSize * gbData->tableau->A[stage_ * nStages + stage_]);
374 ✗ for (int i = 0; i < nStates; i++)
375 {
376 ✗ fODE[i] = ifac * (nlsData->nlsx[i] - gbData->res_const[i]);
377 }
378 }
379 }
380 // copy last calculation of fODE, which should coincide with k[i], here, it yields stage == stage_
381 ✗ memcpy(gbData->k + stage_ * nStates, fODE, nStates*sizeof(double));
382 }
383 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS_V, 0, "GBODE: all stages done.");
384
385 // Apply RK-scheme for determining the approximation at (gbData->time + gbData->stepSize)
386 // y = yold + h * sum(b[stage_] * k[stage_], stage_=1..nStages);
387
388 ✗ for (i=0; i<nStates; i++)
389 {
390 ✗ gbData->y[i] = gbData->yOld[i];
391 ✗ for (stage_=0; stage_<nStages; stage_++)
392 {
393 ✗ gbData->y[i] += gbData->stepSize * gbData->tableau->b[stage_] * (gbData->k + stage_ * nStates)[i];
394 }
395 }
396
397 ✗ if (gbEstimateError(&error_context, &gbData->tableau->error.active) < 0)
398 {
399 ✗ return -1;
400 }
401
402 return 0;
403 }
404
405 /**
406 * @brief Generic diagonal implicit Runge-Kutta step function.
407 *
408 * Only for the fast states (inner integration).
409 *
410 * Internal non-linear equation system will be solved with non-linear solver specified during setup.
411 * Results are saved in y. The selected error estimator writes |error| to errest.
412 *
413 * @param data Runtime data struct.
414 * @param threadData Thread data for error handling.
415 * @param solverInfo Storing Runge-Kutta solver data.
416 * @return int Return 0 on success, -1 on failure.
417 */
418 ✗ int expl_diag_impl_RK_MR(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
419 {
420 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
421 ✗ modelica_real* fODE = sData->realVars + data->modelData->nStates;
422 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
423 ✗ DATA_GBODEF* gbfData = gbData->gbfData;
424
425 ✗ int nStates = gbData->nStates;
426 ✗ int nFastStates = gbData->nFastStates;
427 ✗ int nStages = gbfData->tableau->nStages;
428 NLS_SOLVER_STATUS solved = NLS_FAILED;
429
430 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_NLS)) {
431 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS, 1, "NLS - used values for extrapolation:");
432 ✗ printVector_gbf(OMC_LOG_GBODE_NLS, "xL", gbfData->yv + nStates, nStates, gbfData->tv[1], gbData->nFastStates, gbData->fastStatesIdx);
433 ✗ printVector_gbf(OMC_LOG_GBODE_NLS, "kL", gbfData->kv + nStates, nStates, gbfData->tv[1], gbData->nFastStates, gbData->fastStatesIdx);
434 ✗ printVector_gbf(OMC_LOG_GBODE_NLS, "xR", gbfData->yv, nStates, gbfData->tv[0], gbData->nFastStates, gbData->fastStatesIdx);
435 ✗ printVector_gbf(OMC_LOG_GBODE_NLS, "kR", gbfData->kv, nStates, gbfData->tv[0], gbData->nFastStates, gbData->fastStatesIdx);
436 ✗ messageClose(OMC_LOG_GBODE_NLS);
437 }
438
439 ✗ slowStateCache_merge_left(gbData, gbfData->slowStateCache, gbfData->yOld);
440
441 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
442 {
443 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
444 ✗ gbfData->yOldPacked[fast_idx] = gbfData->yOld[full_idx];
445 }
446
447 ✗ for (int stage = 0; stage < nStages; stage++) {
448 ✗ gbfData->act_stage = stage;
449
450 // set simulation time with respect to the current stage
451 // t = t_0 + c[j]*h
452 ✗ sData->timeValue = gbfData->time + gbfData->tableau->c[stage]*gbfData->stepSize;
453
454 // k[i] = f(tOld + c[i]*h, yOld + h*sum(a[i,j]*k[j], i=j..i))
455 // res = f(tOld + c[i]*h, yOld + h*sum(a[i,j]*k[j], i=j..i-1))
456
457 // check for explicit stage
458 ✗ if (gbfData->tableau->A[stage * nStages + stage] == 0)
459 {
460 // check if kLeft is available and potentially reuse said value
461 ✗ if (gbfData->tableau->isKLeftAvailable && (stage == 0) && gbData->didFastStep)
462 {
463 ✗ copyVector_gbf(fODE, gbfData->kLeft, nFastStates, gbData->fastStatesIdx);
464 }
465 else
466 {
467 // for explicit stages, we update the full res_const buffer as we evaluate the ODE at that point (may be optimized further)
468 ✗ memcpy(gbfData->res_const, gbfData->yOld, nStates * sizeof(double));
469
470 ✗ for (int full_idx = 0; full_idx < nStates; full_idx++)
471 {
472 ✗ for (int s = 0; s < stage; s++)
473 {
474 ✗ gbfData->res_const[full_idx] += gbfData->stepSize * gbfData->tableau->A[stage * nStages + s] * gbfData->k[s * nStates + full_idx];
475 }
476 }
477
478 // calculate the fODE values for the explicit stage
479 ✗ memcpy(sData->realVars, gbfData->res_const, nStates * sizeof(double));
480 ✗ gbode_fODE(data, threadData, &(gbfData->stats.nCallsODE), gbfData->evalSelectionFast);
481 }
482 }
483 else
484 {
485 // for implicit stages, only set the fast states for the NLS
486 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
487 {
488 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
489 ✗ gbfData->res_const[full_idx] = gbfData->yOld[full_idx];
490 ✗ for (int s = 0; s < stage; s++)
491 {
492 ✗ gbfData->res_const[full_idx] += gbfData->stepSize * gbfData->tableau->A[stage * nStages + s] * gbfData->k[s * nStates + full_idx];
493 }
494 }
495
496 // interpolate the slow states on the time of the current stage
497 ✗ slowStateCache_overwrite_stage(gbData, gbfData->slowStateCache,stage, sData->realVars);
498
499 // setting the start vector for the newton step
500 // solve for x: 0 = yold-x + h*(sum(A[i,j]*k[j], i=1..j-1) + A[i,i]*f(t + c[i]*h, x))
501 ✗ NONLINEAR_SYSTEM_DATA* nlsData = gbfData->nlsData;
502
503 ✗ projVector_gbf(nlsData->nlsx, gbfData->yOld, nFastStates, gbData->fastStatesIdx);
504 ✗ memcpy(nlsData->nlsxOld, nlsData->nlsx, nFastStates*sizeof(modelica_real));
505
506 // use help vector gbData->y1 for security reasons
507 ✗ extrapolation_gbf(gbData, gbData->y1, gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize);
508 ✗ projVector_gbf(nlsData->nlsxExtrapolation, gbData->y1, nFastStates, gbData->fastStatesIdx);
509
510 // is the last solution valid and do we use internal nls
511 ✗ modelica_boolean dense_output_valid = (gbfData->extrapolationValid && gbfData->nlsSolverMethod == GB_NLS_INTERNAL);
512
513 ✗ if (gbfData->tableau->svp != NULL && gbfData->tableau->svp->type[stage] == SVP_LINEAR_COMBINATION)
514 {
515 /* linear combination stage-value-predictors (highest priority) */
516 ✗ gbInternalLinearCombinationSVP(gbfData->tableau->svp, stage, nFastStates, gbfData->stepSize, gbfData->kCurrPacked, gbfData->yOldPacked, nlsData->nlsxOld);
517 }
518 ✗ else if (dense_output_valid && gbfData->tableau->svp != NULL && gbfData->tableau->svp->type[stage] == SVP_DENSE_OUTPUT)
519 ✗ {
520 /* dense output stage-value-predictor */
521 ✗ double theta = (gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize - gbfData->extrapolationBaseTime) / gbfData->extrapolationStepSize;
522 ✗ gbfData->tableau->svp->dense_output_predictor(gbfData->tableau, gbfData->yLast, NULL, gbfData->kLast,
523 ✗ theta, gbfData->extrapolationStepSize, nlsData->nlsxOld, 0, NULL, nFastStates);
524 }
525 ✗ else if (dense_output_valid && gbfData->tableau->withDenseOutput)
526 {
527 /* standard dense output if available / possible */
528 ✗ double theta = (gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize - gbfData->extrapolationBaseTime) / gbfData->extrapolationStepSize;
529 ✗ gbfData->tableau->dense_output(gbfData->tableau, gbfData->yLast, NULL, gbfData->kLast,
530 ✗ theta, gbfData->extrapolationStepSize, nlsData->nlsxOld, 0, NULL, nFastStates);
531 }
532
533 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS_V, 0, "Solving NLS of gbf stage %d at time %g", stage+1, gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize);
534 ✗ solved = solveNLS_gb(data, threadData, nlsData, gbData, TRUE);
535
536 ✗ if (solved != NLS_SOLVED) {
537 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "gbodef error: Failed to solve NLS in expl_diag_impl_RK_MR in stage %d at time t=%g", stage+1, gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize);
538 ✗ return -1;
539 }
540
541 // debug residuals
542 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_NLS)) {
543 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS, 1, "NLS - start values and solution of the NLS:");
544 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "xS", nlsData->nlsxExtrapolation, nFastStates, gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize);
545 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "xL", nlsData->nlsx, nFastStates, gbfData->time + gbfData->tableau->c[stage] * gbfData->stepSize);
546 ✗ messageClose(OMC_LOG_GBODE_NLS);
547 }
548
549 ✗ if (/* non explicit stage of (E)SDIRK integrator */ (stage != 0 || gbfData->tableau->A[0] != 0) && gbData->nlsSolverMethod == GB_NLS_INTERNAL)
550 {
551 // reconstruct k_{stage} from the solution, avoids repeated call to functionODE()
552 ✗ double ifac = 1.0 / (gbfData->stepSize * gbfData->tableau->A[stage * nStages + stage]);
553 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
554 {
555 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
556 ✗ fODE[full_idx] = ifac * (nlsData->nlsx[fast_idx] - gbfData->res_const[full_idx]);
557 ✗ sData->realVars[full_idx] = nlsData->nlsx[fast_idx];
558 }
559 }
560 }
561
562 // TODO: make k and y only contain fast states. Almost all structures depend on this full vector: this is a todo for a rewrite of GBODE
563 ✗ if (gbfData->nlsSolverMethod == GB_NLS_INTERNAL)
564 {
565 ✗ int stageOffset = nFastStates * stage;
566
567 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
568 {
569 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
570 ✗ gbfData->kCurrPacked[stageOffset + fast_idx] = fODE[full_idx];
571 }
572 }
573
574 // copy last values of sData->realVars and fODE, which should coincide with x[i] and k[i]
575 // TODO: Make the fast state structures only contains the current flat k's
576 // => change the interpolation routines accordingly
577 // in the interpolation routines gbfData->k is also only used with a fastState mapping, so
578 // this is used effectively anyway
579 ✗ memcpy(gbfData->x + stage * nStates, sData->realVars, nStates*sizeof(double));
580 ✗ memcpy(gbfData->k + stage * nStates, fODE, nStates*sizeof(double));
581 }
582
583 // Apply RK-scheme for determining the approximation at (gbData->time + gbData->stepSize)
584 // y = yold + h * sum(b[stage] * k[stage], stage=1..nStages);
585 // for the fast states only!
586 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++) {
587 ✗ int full_idx = gbData->fastStatesIdx[fast_idx];
588 // y is the new approximation
589 ✗ gbfData->y[full_idx] = gbfData->yOld[full_idx];
590 ✗ for (int stage = 0; stage < nStages; stage++) {
591 ✗ gbfData->y[full_idx] += gbfData->stepSize * gbfData->tableau->b[stage] * (gbfData->k + stage * nStates)[full_idx];
592 }
593 }
594
595 ✗ GB_ERROR_CONTEXT error_context = {data, threadData, gbData, gbfData, TRUE};
596 ✗ if (gbEstimateError(&error_context, &gbfData->tableau->error.active) < 0)
597 {
598 ✗ return -1;
599 }
600
601 return 0;
602 }
603
604 /**
605 * @brief Single implicit Runge-Kutta step.
606 *
607 * @param data Runtime data struct.
608 * @param threadData Thread data for error handling.
609 * @param solverInfo Storing Runge-Kutta solver data.
610 * @return int Return 0 on success, -1 on failure.
611 */
612 ✗ int full_implicit_RK(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
613 {
614 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
615 ✗ modelica_real* fODE = sData->realVars + data->modelData->nStates;
616 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
617
618 ✗ NONLINEAR_SYSTEM_DATA* nlsData = gbData->nlsData;
619
620 int i;
621 ✗ int nStates = data->modelData->nStates;
622 ✗ int nStages = gbData->tableau->nStages;
623
624 NLS_SOLVER_STATUS solved = NLS_FAILED;
625
626 // NLS - used values for extrapolation
627 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_NLS)) {
628 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS, 1, "NLS - used values for extrapolation:");
629 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "xL", gbData->yv + nStates, nStates, gbData->tv[1]);
630 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "kL", gbData->kv + nStates, nStates, gbData->tv[1]);
631 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "xR", gbData->yv, nStates, gbData->tv[0]);
632 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "kR", gbData->kv, nStates, gbData->tv[0]);
633 ✗ messageClose(OMC_LOG_GBODE_NLS);
634 }
635
636 /* Set start values for non-linear solver by extrapolation */
637 ✗ for (int stage = 0; stage < nStages; stage++) {
638 ✗ memcpy(nlsData->nlsx + stage*nStates, gbData->yOld, nStates*sizeof(modelica_real));
639 ✗ memcpy(nlsData->nlsxOld + stage*nStates, gbData->yOld, nStates*sizeof(modelica_real));
640
641 ✗ extrapolation_gb(gbData, nlsData->nlsxExtrapolation + stage*nStates, gbData->time + gbData->tableau->c[stage] * gbData->stepSize);
642 }
643
644 // use dense output extrapolation for all slow states (if SR: then for all states)
645 ✗ if (gbData->time != data->simulationInfo->startTime && !gbData->eventHappened
646 ✗ && gbData->tableau->withDenseOutput && gbData->nlsSolverMethod == GB_NLS_INTERNAL
647 ✗ && gbData->extrapolationBaseTime != INFINITY)
648 {
649 ✗ for (int stage = 0; stage < nStages; stage++) {
650 ✗ double theta = (gbData->time + gbData->tableau->c[stage] * gbData->stepSize - gbData->extrapolationBaseTime) / gbData->extrapolationStepSize;
651 ✗ gbData->tableau->dense_output(gbData->tableau, gbData->yLast, NULL, gbData->kLast,
652 ✗ theta, gbData->extrapolationStepSize, nlsData->nlsxOld + stage*nStates, 0, NULL, nStates);
653 }
654 }
655
656 // zero order hold for all fast states
657 ✗ if (gbData->multi_rate && gbData->nFastStates > 0)
658 {
659 ✗ for (int stage = 0; stage < nStages; stage++)
660 {
661 ✗ for (int fast = 0; fast < gbData->nFastStates; fast++)
662 {
663 ✗ int full = gbData->fastStatesIdx[fast];
664 ✗ nlsData->nlsxOld[stage * nStates + full] = gbData->yOld[full];
665 }
666 }
667 }
668
669 ✗ solved = solveNLS_gb(data, threadData, nlsData, gbData, FALSE);
670
671 ✗ if (solved != NLS_SOLVED) {
672 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "gbode error: Failed to solve NLS in full_implicit_RK at time t=%g", gbData->time);
673 ✗ return -1;
674 }
675
676 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_NLS)) {
677 ✗ infoStreamPrint(OMC_LOG_GBODE_NLS, 1, "NLS - start values and solution of the NLS:");
678 ✗ for (int stage = 0; stage < nStages; stage++) {
679 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "xS", nlsData->nlsxExtrapolation + stage*nStates, nStates, gbData->time + gbData->tableau->c[stage] * gbData->stepSize);
680 ✗ printVector_gb(OMC_LOG_GBODE_NLS, "xL", nlsData->nlsx + stage*nStates, nStates, gbData->time + gbData->tableau->c[stage] * gbData->stepSize);
681 }
682 ✗ messageClose(OMC_LOG_GBODE_NLS);
683 }
684
685
686 // Apply RK-scheme for determining the approximation at (gbData->time + gbData->stepSize)
687 // y = yold + h * sum(b[stage_] * k[stage_], stage_=1..nStages);
688
689 // calculate y(t_n+1)
690 ✗ for (i = 0; i < nStates; i++) {
691 ✗ gbData->y[i] = gbData->yOld[i];
692 ✗ for (int stage = 0; stage < nStages; stage++) {
693 ✗ gbData->y[i] += gbData->stepSize * gbData->tableau->b[stage] * (gbData->k + stage * nStates)[i];
694 }
695 }
696
697 ✗ GB_ERROR_CONTEXT error_context = {data, threadData, gbData, NULL, FALSE};
698 ✗ if (gbEstimateError(&error_context, &gbData->tableau->error.active) < 0)
699 {
700 return -1;
701 }
702
703 // copy the whole solution vector to the inner buffer (for latter extrapolation and dense output)
704 ✗ memcpy(gbData->x, nlsData->nlsx, nlsData->size*sizeof(double));
705
706 ✗ return 0;
707 }
708
709 /**
710 * @brief Single implicit Runge-Kutta step.
711 *
712 * @param data Runtime data struct.
713 * @param threadData Thread data for error handling.
714 * @param solverInfo Storing Runge-Kutta solver data.
715 * @return int Return 0 on success, -1 on failure.
716 */
717 ✗ int full_implicit_RK_MR(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
718 {
719 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
720 modelica_real* fODE = sData->realVars + data->modelData->nStates;
721 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
722 ✗ DATA_GBODEF* gbfData = gbData->gbfData;
723 ✗ BUTCHER_TABLEAU *tableau = gbfData->tableau;
724
725 ✗ NONLINEAR_SYSTEM_DATA* nlsData = gbfData->nlsData;
726
727 ✗ int nStates = gbData->nStates;
728 ✗ int nStages = tableau->nStages;
729 ✗ int nFastStates = gbData->nFastStates;
730 ✗ int *fastStatesIdx = gbData->fastStatesIdx;
731
732 NLS_SOLVER_STATUS solved = NLS_FAILED;
733
734 // Attention: as currently all structures in GBODEF_DATA rely on indirect indexing to fast states, e.g.
735 // y(fast_state_i) = gbfData->y[fastStatesIdx[i]], instead of direct (flat) access gbfData->y[i]
736 // usage in gbnls=internal is very inconvenient. Therefore, we use the fields gbfData->yOldPacked and
737 // gbfData->kCurrPacked which represent the yOld and k fields but packed as described above
738 // Thus, internal NLS writes the solutions to kCurrPacked and uses the packed values yOldPacked, so
739 // differences between fast and slow steps is minimal, we can use BLAS routines to full extend and the code is not
740 // that nested. However, we need to extract these solutions to x and k fields at the end though.
741
742 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
743 {
744 ✗ int full_idx = fastStatesIdx[fast_idx];
745 ✗ gbfData->yOldPacked[fast_idx] = gbfData->yOld[full_idx];
746 }
747
748 /* Set start values for non-linear solver by extrapolation */
749 ✗ for (int stage = 0; stage < nStages; stage++) {
750 ✗ int offset = stage * nFastStates;
751 ✗ memcpy(&nlsData->nlsx[offset], gbfData->yOldPacked, nFastStates * sizeof(double));
752 ✗ memcpy(&nlsData->nlsxOld[offset], gbfData->yOldPacked, nFastStates * sizeof(double));
753 }
754
755 ✗ if (gbfData->tableau->withDenseOutput && gbfData->extrapolationValid && gbfData->nlsSolverMethod == GB_NLS_INTERNAL)
756 {
757 ✗ for (int stage = 0; stage < nStages; stage++)
758 {
759 ✗ int offset = stage * nFastStates;
760 ✗ double theta = (gbfData->time + tableau->c[stage] * gbfData->stepSize - gbfData->extrapolationBaseTime) / gbfData->extrapolationStepSize;
761 ✗ tableau->dense_output(tableau, gbfData->yLast, NULL, gbfData->kLast,
762 ✗ theta, gbfData->extrapolationStepSize, &nlsData->nlsxOld[offset], 0, NULL, nFastStates);
763 }
764 }
765
766 ✗ solved = solveNLS_gb(data, threadData, nlsData, gbData, TRUE);
767
768 ✗ if (solved != NLS_SOLVED) {
769 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "gbode error: Failed to solve NLS in full_implicit_RK_MR at time t=%g", gbData->time);
770 ✗ return -1;
771 }
772
773 // calculate x, y, k
774 ✗ for (int fast_idx = 0; fast_idx < nFastStates; fast_idx++)
775 {
776 ✗ int full_idx = fastStatesIdx[fast_idx];
777 ✗ gbfData->y[full_idx] = gbfData->yOld[full_idx];
778 ✗ for (int stage = 0; stage < nStages; stage++)
779 {
780 ✗ int offset_fast = stage * nFastStates;
781 ✗ int offset_full = stage * nStates;
782 ✗ gbfData->x[offset_full + full_idx] = nlsData->nlsx[offset_fast + fast_idx];
783 ✗ gbfData->y[full_idx] += gbfData->stepSize * gbfData->tableau->b[stage] * gbfData->kCurrPacked[offset_fast + fast_idx];
784 ✗ gbfData->k[offset_full + full_idx] = gbfData->kCurrPacked[offset_fast + fast_idx];
785 }
786 }
787
788 ✗ GB_ERROR_CONTEXT error_context = {data, threadData, gbData, gbfData, TRUE};
789 ✗ if (gbEstimateError(&error_context, &gbfData->tableau->error.active) < 0)
790 {
791 ✗ return -1;
792 }
793
794 return 0;
795 }
796
797
798 /**
799 * @brief
800 *
801 * @param data
802 * @param threadData
803 * @param solverInfo
804 * @return int
805 */
806 ✗ int gbodef_richardson(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
807 {
808 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
809 ✗ modelica_real* fODE = sData->realVars + data->modelData->nStates;
810 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
811 ✗ DATA_GBODEF* gbfData = gbData->gbfData;
812
813 double stepSize, lastStepSize, timeValue;
814 int step_info, p;
815 ✗ int nStates = gbfData->nStates;
816 int i;
817
818 // assumption yLeft and yOld coincide!!!
819 ✗ timeValue = gbfData->time;
820 ✗ stepSize = gbfData->stepSize;
821 ✗ lastStepSize = gbfData->lastStepSize;
822 ✗ p = gbfData->tableau->order_b;
823
824 ✗ if (!gbfData->isExplicit) {
825 // Store relevant part of the ring buffer, which is used for extrapolation
826 ✗ for (i = 0; i < 2; i++) {
827 ✗ gbData->tr[i] = gbfData->tv[i];
828 ✗ memcpy(gbData->yr + i * nStates, gbfData->yv + i * nStates, nStates * sizeof(double));
829 ✗ memcpy(gbData->kr + i * nStates, gbfData->kv + i * nStates, nStates * sizeof(double));
830 }
831 }
832
833 ✗ gbfData->stepSize = gbfData->stepSize/2;
834 ✗ step_info = gbfData->step_fun(data, threadData, solverInfo);
835 ✗ if (step_info != 0) {
836 ✗ stepSize = stepSize/2;
837 ✗ lastStepSize = lastStepSize/2;
838 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "Failure: gbode Richardson extrapolation (first half step)");
839 } else {
840 // debug the approximations after performed step
841 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
842 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Richardson extrapolation (first 1/2 step) approximation:");
843 ✗ printVector_gb(OMC_LOG_GBODE, " y", gbfData->y, nStates, gbfData->time + gbfData->stepSize);
844 ✗ printVector_gb(OMC_LOG_GBODE, "yt", gbfData->yt, nStates, gbfData->time + gbfData->stepSize);
845 ✗ messageClose(OMC_LOG_GBODE);
846 }
847 ✗ gbfData->time += gbfData->stepSize;
848 ✗ gbfData->lastStepSize = gbfData->stepSize;
849 ✗ memcpy(gbfData->yOld, gbfData->y, nStates * sizeof(double));
850
851 // prepare for the extrapolation
852 ✗ if (!gbfData->isExplicit) {
853 ✗ sData->timeValue = gbfData->time;
854 ✗ memcpy(sData->realVars, gbfData->y, nStates*sizeof(double));
855 ✗ gbode_fODE(data, threadData, &(gbfData->stats.nCallsODE), gbfData->evalSelectionFast);
856 ✗ gbfData->tv[1] = gbfData->tv[0];
857 ✗ memcpy(gbfData->yv + nStates, gbfData->yv, nStates * sizeof(double));
858 ✗ memcpy(gbfData->kv + nStates, gbfData->kv, nStates * sizeof(double));
859 ✗ gbfData->tv[0] = gbfData->time;
860 ✗ memcpy(gbfData->yv, gbfData->y, nStates * sizeof(double));
861 ✗ memcpy(gbfData->kv, fODE, nStates * sizeof(double));
862 }
863
864 ✗ step_info = gbfData->step_fun(data, threadData, solverInfo);
865 ✗ if (step_info != 0) {
866 ✗ stepSize = stepSize/2;
867 ✗ lastStepSize = lastStepSize/2;
868 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "Failure: gbode Richardson extrapolation (second half step)");
869 } else {
870 // debug the approximations after performed step
871 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
872 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Richardson extrapolation (second 1/2 step) approximation:");
873 ✗ printVector_gb(OMC_LOG_GBODE, " y", gbfData->y, nStates, gbfData->time + gbfData->stepSize);
874 ✗ printVector_gb(OMC_LOG_GBODE, "yt", gbfData->yt, nStates, gbfData->time + gbfData->stepSize);
875 ✗ messageClose(OMC_LOG_GBODE);
876 }
877 ✗ memcpy(gbfData->y1, gbfData->y, nStates * sizeof(double));
878
879 // prepare for the extrapolation
880 ✗ if (!gbfData->isExplicit) {
881 ✗ sData->timeValue = gbfData->time + gbfData->stepSize;
882 ✗ memcpy(sData->realVars, gbfData->y, nStates*sizeof(double));
883 ✗ gbode_fODE(data, threadData, &(gbfData->stats.nCallsODE), gbfData->evalSelectionFast);
884 ✗ gbfData->tv[0] = gbfData->time;
885 ✗ memcpy(gbfData->yv, gbfData->y, nStates * sizeof(double));
886 ✗ memcpy(gbfData->kv, fODE, nStates * sizeof(double));
887 }
888
889 // restore yOld
890 ✗ gbfData->time = timeValue;
891 ✗ gbfData->stepSize = stepSize;
892 ✗ gbfData->lastStepSize = lastStepSize;
893 ✗ memcpy(gbfData->yOld, gbfData->yLeft, nStates * sizeof(double));
894 ✗ step_info = gbfData->step_fun(data, threadData, solverInfo);
895 ✗ if (step_info != 0) {
896 ✗ stepSize = stepSize/2;
897 ✗ lastStepSize = lastStepSize/2;
898 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "Failure: gbode Richardson extrapolation (full step)");
899 } else {
900 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
901 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Richardson extrapolation (full step) approximation");
902 ✗ printVector_gb(OMC_LOG_GBODE, " y", gbfData->y, nStates, gbfData->time + gbfData->stepSize);
903 ✗ printVector_gb(OMC_LOG_GBODE, "yt", gbfData->yt, nStates, gbfData->time + gbfData->stepSize);
904 ✗ messageClose(OMC_LOG_GBODE);
905 }
906 }
907 }
908 }
909
910 // Restore time values and step size
911 ✗ gbfData->time = timeValue;
912 ✗ gbfData->stepSize = stepSize;
913 ✗ gbfData->lastStepSize = lastStepSize;
914 ✗ memcpy(gbfData->yOld, gbfData->yLeft, nStates * sizeof(double));
915 ✗ if (!gbfData->isExplicit) {
916 // Restore ring buffer
917 ✗ for (i = 0; i < 2; i++) {
918 ✗ gbfData->tv[i] = gbData->tr[i];
919 ✗ memcpy(gbfData->yv + i * nStates, gbData->yr + i * nStates, nStates * sizeof(double));
920 ✗ memcpy(gbfData->kv + i * nStates, gbData->kr + i * nStates, nStates * sizeof(double));
921 }
922 }
923 ✗ if (!step_info) {
924 // Extrapolate values based on order of the scheme
925 ✗ double richardsonFactor = pow(2., p);
926 ✗ for (i = 0; i < nStates; i++) {
927 ✗ double y_extrapolated = (richardsonFactor * gbfData->y1[i] - gbfData->y[i]) / (richardsonFactor - 1);
928 ✗ gbfData->yt[i] = gbfData->y[i] - y_extrapolated;
929 }
930 }
931
932 ✗ return step_info;
933 }
934
935 /**
936 * @brief
937 *
938 * @param data
939 * @param threadData
940 * @param solverInfo
941 * @return int
942 */
943 ✗ int gbode_richardson(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
944 {
945 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
946 ✗ modelica_real* fODE = sData->realVars + data->modelData->nStates;
947 ✗ DATA_GBODE* gbData = (DATA_GBODE*)solverInfo->solverData;
948
949 double stepSize, lastStepSize, timeValue;
950 int step_info, p;
951 ✗ int nStates = gbData->nStates;
952 int i;
953
954 // assumption yLeft and yOld coincide!!!
955 ✗ timeValue = gbData->time;
956 ✗ stepSize = gbData->stepSize;
957 ✗ lastStepSize = gbData->lastStepSize;
958 ✗ p = gbData->tableau->order_b;
959
960 ✗ if (!gbData->isExplicit) {
961 // Store relevant part of the ring buffer, which is used for extrapolation
962 ✗ for (i = 0; i < 2; i++) {
963 ✗ gbData->tr[i] = gbData->tv[i];
964 ✗ memcpy(gbData->yr + i * nStates, gbData->yv + i * nStates, nStates * sizeof(double));
965 ✗ memcpy(gbData->kr + i * nStates, gbData->kv + i * nStates, nStates * sizeof(double));
966 }
967 }
968
969 ✗ gbData->stepSize = gbData->stepSize/2;
970 ✗ step_info = gbData->step_fun(data, threadData, solverInfo);
971 ✗ if (step_info != 0) {
972 ✗ stepSize = stepSize/2;
973 ✗ lastStepSize = lastStepSize/2;
974 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "Failure: gbode Richardson extrapolation (first half step)");
975 } else {
976 // debug the approximations after performed step
977 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
978 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Richardson extrapolation (first 1/2 step) approximation:");
979 ✗ printVector_gb(OMC_LOG_GBODE, " y", gbData->y, nStates, gbData->time + gbData->stepSize);
980 ✗ printVector_gb(OMC_LOG_GBODE, "yt", gbData->yt, nStates, gbData->time + gbData->stepSize);
981 ✗ messageClose(OMC_LOG_GBODE);
982 }
983 ✗ gbData->time += gbData->stepSize;
984 ✗ gbData->lastStepSize = gbData->stepSize;
985 ✗ memcpy(gbData->yOld, gbData->y, nStates * sizeof(double));
986
987 // prepare for the extrapolation
988 ✗ if (!gbData->isExplicit) {
989 ✗ sData->timeValue = gbData->time;
990 ✗ memcpy(sData->realVars, gbData->y, nStates*sizeof(double));
991 ✗ gbode_fODE(data, threadData, &(gbData->stats.nCallsODE), NULL);
992 ✗ gbData->tv[1] = gbData->tv[0];
993 ✗ memcpy(gbData->yv + nStates, gbData->yv, nStates * sizeof(double));
994 ✗ memcpy(gbData->kv + nStates, gbData->kv, nStates * sizeof(double));
995 ✗ gbData->tv[0] = gbData->time;
996 ✗ memcpy(gbData->yv, gbData->y, nStates * sizeof(double));
997 ✗ memcpy(gbData->kv, fODE, nStates * sizeof(double));
998 }
999
1000 ✗ step_info = gbData->step_fun(data, threadData, solverInfo);
1001 ✗ if (step_info != 0) {
1002 ✗ stepSize = stepSize/2;
1003 ✗ lastStepSize = lastStepSize/2;
1004 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "Failure: gbode Richardson extrapolation (second half step)");
1005 } else {
1006 // debug the approximations after performed step
1007 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1008 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Richardson extrapolation (second 1/2 step) approximation:");
1009 ✗ printVector_gb(OMC_LOG_GBODE, " y", gbData->y, nStates, gbData->time + gbData->stepSize);
1010 ✗ printVector_gb(OMC_LOG_GBODE, "yt", gbData->yt, nStates, gbData->time + gbData->stepSize);
1011 ✗ messageClose(OMC_LOG_GBODE);
1012 }
1013 ✗ memcpy(gbData->y1, gbData->y, nStates * sizeof(double));
1014
1015 // prepare for the extrapolation
1016 ✗ if (!gbData->isExplicit) {
1017 ✗ sData->timeValue = gbData->time + gbData->stepSize;
1018 ✗ memcpy(sData->realVars, gbData->y, nStates*sizeof(double));
1019 ✗ gbode_fODE(data, threadData, &(gbData->stats.nCallsODE), NULL);
1020 ✗ gbData->tv[0] = gbData->time;
1021 ✗ memcpy(gbData->yv, gbData->y, nStates * sizeof(double));
1022 ✗ memcpy(gbData->kv, fODE, nStates * sizeof(double));
1023 }
1024
1025 // restore yOld
1026 ✗ gbData->time = timeValue;
1027 ✗ gbData->stepSize = stepSize;
1028 ✗ gbData->lastStepSize = lastStepSize;
1029 ✗ memcpy(gbData->yOld, gbData->yLeft, nStates * sizeof(double));
1030 ✗ step_info = gbData->step_fun(data, threadData, solverInfo);
1031 ✗ if (step_info != 0) {
1032 ✗ stepSize = stepSize/2;
1033 ✗ lastStepSize = lastStepSize/2;
1034 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_SOLVER)) warningStreamPrint(OMC_LOG_SOLVER, 0, "Failure: gbode Richardson extrapolation (full step)");
1035 } else {
1036 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE)) {
1037 ✗ infoStreamPrint(OMC_LOG_GBODE, 1, "Richardson extrapolation (full step) approximation");
1038 ✗ printVector_gb(OMC_LOG_GBODE, " y", gbData->y, nStates, gbData->time + gbData->stepSize);
1039 ✗ printVector_gb(OMC_LOG_GBODE, "yt", gbData->yt, nStates, gbData->time + gbData->stepSize);
1040 ✗ messageClose(OMC_LOG_GBODE);
1041 }
1042 }
1043 }
1044 }
1045
1046 // Restore time values and step size
1047 ✗ gbData->time = timeValue;
1048 ✗ gbData->stepSize = stepSize;
1049 ✗ gbData->lastStepSize = lastStepSize;
1050 ✗ memcpy(gbData->yOld, gbData->yLeft, nStates * sizeof(double));
1051
1052 ✗ if (!gbData->isExplicit) {
1053 // Restore ring buffer
1054 ✗ for (i = 0; i < 2; i++) {
1055 ✗ gbData->tv[i] = gbData->tr[i];
1056 ✗ memcpy(gbData->yv + i * nStates, gbData->yr + i * nStates, nStates * sizeof(double));
1057 ✗ memcpy(gbData->kv + i * nStates, gbData->kr + i * nStates, nStates * sizeof(double));
1058 }
1059 }
1060
1061 ✗ if (!step_info) {
1062 // Extrapolate values based on order of the scheme
1063 ✗ double richardsonFactor = pow(2., p);
1064 ✗ for (i = 0; i < nStates; i++) {
1065 ✗ double y_extrapolated = (richardsonFactor * gbData->y1[i] - gbData->y[i]) / (richardsonFactor - 1);
1066 ✗ gbData->yt[i] = gbData->y[i] - y_extrapolated;
1067 }
1068 }
1069
1070 ✗ return step_info;
1071 }
1072