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 / 212
Functions: 0.0% 0 / 0 / 5
Branches: 0.0% 0 / 0 / 68
OMCompiler/SimulationRuntime/c/simulation/solver/sym_solver_ssc.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 sym_solver_ssc.c
29 */
30
31 #include <string.h>
32
33 #include "../../util/omc_error.h"
34 #include "model_help.h"
35
36 #include "sym_solver_ssc.h"
37 #include "external_input.h"
38
39
40 int first_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
41 int generateTwoApproximationsOfDifferentOrder(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
42
43
44 /*! \fn allocateSymEulerImp
45 *
46 * Function allocates memory needed for implicit symbolic euler with step size control.
47 *
48 *
49 */
50 ✗ int allocateSymSolverSsc(SOLVER_INFO* solverInfo, int size)
51 {
52 ✗ DATA_SYM_SOLVER_SSC* userdata = (DATA_SYM_SOLVER_SSC*) malloc(sizeof(DATA_SYM_SOLVER_SSC));
53 ✗ solverInfo->solverData = (void*) userdata;
54
55 ✗ userdata->firstStep = 1;
56 ✗ userdata->y05= malloc(sizeof(double)*size);
57 ✗ userdata->y1 = malloc(sizeof(double)*size);
58 ✗ userdata->y2 = malloc(sizeof(double)*size);
59 ✗ userdata->radauVarsOld = malloc(sizeof(double)*size);
60 ✗ userdata->radauVars = malloc(sizeof(double)*size);
61 ✗ userdata->der_x0 = malloc(sizeof(double)*size);
62
63 /* initialize stats */
64 ✗ userdata->stepsDone = 0;
65 ✗ userdata->evalFunctionODE = 0;
66
67 ✗ userdata->radauStepSizeOld = 0;
68 ✗ return 0;
69 }
70
71 /*! \fn freeSymEulerImp
72 *
73 * Memory needed for solver is set free.
74 */
75 ✗ int freeSymSolverSsc(SOLVER_INFO* solverInfo)
76 {
77 ✗ DATA_SYM_SOLVER_SSC* userdata = (DATA_SYM_SOLVER_SSC*) solverInfo->solverData;
78
79 ✗ free(userdata->y05);
80 ✗ free(userdata->y1);
81 ✗ free(userdata->y2);
82 ✗ free(userdata->radauVarsOld);
83 ✗ free(userdata->radauVars);
84
85 ✗ return 0;
86 }
87
88 /*! \fn sym_euler_im_with_step_size_control_step
89 *
90 * Function does one implicit euler step
91 * and calculates step size for the next step
92 * using the implicit midpoint rule
93 *
94 */
95 ✗ int sym_solver_ssc_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
96 {
97 int retVal = 0;
98 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
99 ✗ SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
100 ✗ DATA_SYM_SOLVER_SSC* userdata = (DATA_SYM_SOLVER_SSC*)solverInfo->solverData;
101 ✗ modelica_real* stateDer = sDataOld->realVars + data->modelData->nStates;
102
103 double sc, err, a, b, diff;
104 ✗ double Atol = data->simulationInfo->tolerance, Rtol = data->simulationInfo->tolerance;
105 int i,j;
106 double fac = 0.9;
107 double facmax = 3.5;
108 double facmin = 0.3;
109 ✗ double saveTime = sDataOld->timeValue;
110 ✗ double targetTime = sDataOld->timeValue + solverInfo->currentStepSize;
111
112
113 ✗ if (userdata->firstStep || solverInfo->didEventStep == 1)
114 {
115 ✗ retVal = first_step(data, threadData, solverInfo);
116 ✗ userdata->radauStepSizeOld = 0;
117
118 ✗ if (retVal != 0)
119 {
120 return -1;
121 }
122 }
123
124 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "new step: time=%e", userdata->radauTime);
125 ✗ while (userdata->radauTime < targetTime)
126 {
127 do
128 {
129 ✗ retVal = generateTwoApproximationsOfDifferentOrder(data, threadData, solverInfo);
130
131 ✗ for (i=0; i<data->modelData->nStates; i++)
132 {
133 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "y1[%d]=%e", i, userdata->y1[i]);
134 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "y2[%d]=%e", i, userdata->y2[i]);
135 }
136
137 /*** calculate error ***/
138 ✗ for (i=0, err=0.0; i<data->modelData->nStates; i++)
139 {
140 ✗ sc = Atol + fmax(fabs(userdata->y2[i]),fabs(userdata->y1[i]))*Rtol;
141 ✗ diff = userdata->y2[i]-userdata->y1[i];
142 ✗ err += (diff*diff)/(sc*sc);
143 }
144
145 ✗ err /= data->modelData->nStates;
146
147 ✗ userdata->stepsDone += 1;
148 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "err = %e", err);
149 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "min(facmax, max(facmin, fac*sqrt(1/err))) = %e", fmin(facmax, fmax(facmin, fac*pow(1.0/err, 4))));
150
151
152 /* update step size */
153 ✗ userdata->radauStepSizeOld = userdata->radauStepSize;
154 ✗ userdata->radauStepSize *= fmin(facmax, fmax(facmin, fac*sqrt(1.0/err)));
155
156 ✗ if (isnan(userdata->radauStepSize) || userdata->radauStepSize < 1e-13)
157 {
158 ✗ userdata->radauStepSize = 1e-13;
159 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Desired step to small try next one");
160 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Interpolate linear");
161
162 /* explicit euler step*/
163 ✗ for(i = 0; i < data->modelData->nStates; i++)
164 {
165 ✗ sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
166 }
167 ✗ sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
168 ✗ solverInfo->currentTime = sData->timeValue;
169
170 ✗ userdata->radauTimeOld = userdata->radauTime;
171 ✗ userdata->radauTime += userdata->radauStepSizeOld;
172
173 ✗ memcpy(userdata->radauVarsOld, userdata->radauVars, data->modelData->nStates*sizeof(double));
174 ✗ memcpy(userdata->radauVars, userdata->y2, data->modelData->nStates*sizeof(double));
175
176 break;
177 }
178
179 ✗ } while (err > 1.0 );
180
181 ✗ userdata->radauTimeOld = userdata->radauTime;
182
183 ✗ userdata->radauTime += userdata->radauStepSizeOld;
184
185 ✗ memcpy(userdata->radauVarsOld, userdata->radauVars, data->modelData->nStates*sizeof(double));
186 ✗ memcpy(userdata->radauVars, userdata->y2, data->modelData->nStates*sizeof(double));
187 }
188
189 ✗ sDataOld->timeValue = saveTime;
190 ✗ solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
191 ✗ sData->timeValue = solverInfo->currentTime;
192
193 ✗ if (userdata->radauTime - userdata->radauTimeOld > 1e-13 && userdata->radauStepSizeOld > 1e-13)
194 {
195 /* linear interpolation */
196 ✗ for (i=0; i<data->modelData->nStates; i++)
197 {
198 ✗ sData->realVars[i] = (userdata->radauVars[i] * (sData->timeValue - userdata->radauTimeOld) + userdata->radauVarsOld[i] * (userdata->radauTime - sData->timeValue))/(userdata->radauTime - userdata->radauTimeOld);
199 }
200
201 /* update first derivative */
202 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "Time %e", sData->timeValue);
203 ✗ for(i=0; i<data->modelData->nStates; ++i)
204 {
205 ✗ a = 4.0 * (userdata->y2[i] - 2.0 * userdata->y05[i] + userdata->radauVarsOld[i]) / (userdata->radauStepSizeOld * userdata->radauStepSizeOld);
206 ✗ b = 2.0 * (userdata->y2[i] - userdata->y05[i])/userdata->radauStepSizeOld - userdata->radauTime * a;
207 ✗ stateDer[i] = a * sData->timeValue + b;
208 }
209 }
210 else
211 {
212 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Desired step to small try next one");
213 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Interpolate linear");
214
215 /* explicit euler step*/
216 ✗ for(i = 0; i < data->modelData->nStates; i++)
217 {
218 ✗ sData->realVars[i] = sDataOld->realVars[i] + stateDer[i] * solverInfo->currentStepSize;
219 }
220 ✗ sData->timeValue = solverInfo->currentTime + solverInfo->currentStepSize;
221 ✗ solverInfo->currentTime = sData->timeValue;
222
223 ✗ userdata->radauTimeOld = userdata->radauTime;
224 ✗ userdata->radauTime += userdata->radauStepSizeOld;
225
226 ✗ memcpy(userdata->radauVarsOld, userdata->radauVars, data->modelData->nStates*sizeof(double));
227 ✗ memcpy(userdata->radauVars, userdata->y2, data->modelData->nStates*sizeof(double));
228 }
229
230 /* update step size */
231 ✗ data->simulationInfo->inlineData->dt = userdata->radauStepSize;
232 ✗ userdata->solverStepSize = userdata->radauStepSizeOld;
233 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "Step done to %f with step size = %e", sData->timeValue, userdata->solverStepSize);
234
235
236 ✗ return retVal;
237 }
238
239 /*! \fn first_step
240 *
241 * function initializes values and sets
242 * initial step size
243 *
244 */
245 ✗ int first_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
246 {
247 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
248 ✗ SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
249 ✗ DATA_SYM_SOLVER_SSC* userdata = (DATA_SYM_SOLVER_SSC*)solverInfo->solverData;
250 ✗ const int n = data->modelData->nStates;
251 ✗ modelica_real* stateDer = sData->realVars + data->modelData->nStates;
252 ✗ modelica_real* stateDerOld = sDataOld->realVars + data->modelData->nStates;
253 double sc, d, d0 = 0.0, d1 = 0.0, d2 = 0.0, h0, h1, delta_ti, infNorm, sum = 0;
254 ✗ double Atol = data->simulationInfo->tolerance, Rtol = data->simulationInfo->tolerance;
255 int i,j,retVal;
256 /* it seems that jacobian is not used yet!
257 #if defined(_MSC_VER) // handle crap compilers
258 double *jacobian = (double*)malloc(n*n*sizeof(double));
259 #else
260 double jacobian[n*n];
261 #endif
262 */
263
264 /* initialize radau values */
265 ✗ for (i=0; i<data->modelData->nStates; i++)
266 {
267 ✗ userdata->radauVars[i] = sData->realVars[i];
268 ✗ userdata->radauVarsOld[i] = sDataOld->realVars[i];
269 }
270
271 ✗ userdata->radauTime = sDataOld->timeValue;
272 ✗ userdata->radauTimeOld = sDataOld->timeValue;
273
274 ✗ userdata->firstStep = 0;
275 ✗ solverInfo->didEventStep = 0;
276
277 ✗ if (compiledWithSymSolver == 2) /* compiled with symSolver - explicit euler*/
278 {
279 /*** calculate starting step size 1st Version ***/
280
281 /* update step size */
282 ✗ data->simulationInfo->inlineData->dt = 1e-8;
283
284 /* evaluate function */
285 ✗ externalInputUpdate(data);
286 ✗ data->callback->input_function(data, threadData);
287 ✗ retVal = data->callback->symbolicInlineSystems(data, threadData);
288
289 ✗ for (i=0; i<data->modelData->nStates; i++)
290 {
291 ✗ stateDer[i] = (sData->realVars[i] - sDataOld->realVars[i])/data->simulationInfo->inlineData->dt;
292 }
293
294 ✗ if(retVal != 0){
295 return -1;
296 }
297
298 ✗ for (i=0; i<data->modelData->nStates; i++)
299 {
300 ✗ sc = Atol + fabs(sDataOld->realVars[i])*Rtol;
301 ✗ d0 += ((sDataOld->realVars[i] * sDataOld->realVars[i])/(sc*sc));
302 ✗ d1 += ((stateDer[i] * stateDer[i]) / (sc*sc));
303 }
304 ✗ d0 /= data->modelData->nStates;
305 ✗ d1 /= data->modelData->nStates;
306
307 ✗ d0 = sqrt(d0);
308 ✗ d1 = sqrt(d1);
309
310
311 ✗ for (i=0; i<data->modelData->nStates; i++)
312 {
313 ✗ userdata->der_x0[i] = stateDer[i];
314 }
315
316 ✗ if (d0 < 1e-5 || d1 < 1e-5)
317 {
318 h0 = 1e-6;
319 }
320 else
321 {
322 ✗ h0 = 0.01 * d0/d1;
323 }
324
325
326 ✗ for (i=0; i<data->modelData->nStates; i++)
327 {
328 ✗ sData->realVars[i] = userdata->radauVars[i] + stateDer[i] * h0;
329 }
330 ✗ sData->timeValue += h0;
331
332 /* update step size */
333 ✗ data->simulationInfo->inlineData->dt = h0;
334
335 /* evaluate function */
336 ✗ externalInputUpdate(data);
337 ✗ data->callback->input_function(data, threadData);
338 ✗ retVal = data->callback->symbolicInlineSystems(data, threadData);
339
340 ✗ for (i=0; i<data->modelData->nStates; i++)
341 {
342 ✗ stateDer[i] = (sData->realVars[i] - sDataOld->realVars[i])/data->simulationInfo->inlineData->dt;
343 }
344
345 ✗ for (i=0; i<data->modelData->nStates; i++)
346 {
347 ✗ sc = Atol + fabs(userdata->radauVars[i])*Rtol;
348 ✗ d2 += ((stateDer[i]-userdata->der_x0[i])*(stateDer[i]-userdata->der_x0[i])/(sc*sc));
349 }
350
351 ✗ d2 = sqrt(d2);
352 ✗ d2 /= h0;
353
354 ✗ d = fmax(d1,d2);
355
356 ✗ if (d > 1e-15)
357 {
358 ✗ h1 = sqrt(0.01/d);
359 }
360 else
361 {
362 ✗ h1 = fmax(1e-6, h0*1e-3);
363 }
364
365 ✗ userdata->radauStepSize = 0.5*fmin(100*h0,h1);
366 ✗ data->simulationInfo->inlineData->dt = userdata->radauStepSize;
367
368 /* end calculation new step size */
369 }
370 else
371 {
372 ✗ userdata->radauStepSize = 0.5*solverInfo->currentStepSize;
373 }
374 /*
375 #if defined(_MSC_VER) // handle crap compilers
376 free(jacobian)
377 #endif
378 */
379 return 0;
380 }
381
382
383 /*! \fn generateTwoApproximationsOfDifferentOrder
384 *
385 * Function generates two approximations of
386 * different convergence order for step
387 * size control (stored in userdata->y1, userdata->y2)
388 *
389 */
390 ✗ int generateTwoApproximationsOfDifferentOrder(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
391 {
392 int retVal = 0;
393 int i;
394 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
395 ✗ SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
396 ✗ DATA_SYM_SOLVER_SSC* userdata = (DATA_SYM_SOLVER_SSC*)solverInfo->solverData;
397 modelica_real* stateDer = sDataOld->realVars + data->modelData->nStates;
398 ✗ if (compiledWithSymSolver == 1) /* compiled with implicit symbolic euler */
399 {
400 /*** do one step with half step size ***/
401 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "radauStepSize = %e", userdata->radauStepSize);
402
403 /* update step size */
404 ✗ userdata->radauStepSize /= 2;
405 ✗ data->simulationInfo->inlineData->dt = userdata->radauStepSize;
406
407 /* update time */
408 ✗ sDataOld->timeValue = userdata->radauTime;
409 ✗ solverInfo->currentTime = userdata->radauTime + userdata->radauStepSize;
410 ✗ sData->timeValue = solverInfo->currentTime;
411
412 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "first system time = %e", sData->timeValue);
413
414 /* update algebraicOld values */
415 ✗ memcpy(data->simulationInfo->inlineData->algOldVars, userdata->radauVars, data->modelData->nStates * sizeof(double));
416
417 /* evaluate function */
418 ✗ externalInputUpdate(data);
419 ✗ data->callback->input_function(data, threadData);
420 ✗ retVal = data->callback->symbolicInlineSystems(data, threadData);
421
422 ✗ if(retVal != 0){
423 return -1;
424 }
425
426 /* save values in y05 */
427 ✗ memcpy(userdata->y05, sData->realVars, data->modelData->nStates*sizeof(double));
428
429 /*** extrapolate values in y1 (= y0 + h * f(y(t+h/2),t+h/2)) ***/
430 ✗ for (i=0; i<data->modelData->nStates; i++)
431 {
432 ✗ userdata->y1[i] = 2.0 * userdata->y05[i] - userdata->radauVars[i];
433 }
434
435 /*** do another step with half step size ***/
436 ✗ memcpy(data->simulationInfo->inlineData->algOldVars, userdata->y05, data->modelData->nStates * sizeof(double));
437
438 /* update time */
439 ✗ sDataOld->timeValue = userdata->radauTime + userdata->radauStepSize;
440 ✗ solverInfo->currentTime = userdata->radauTime + 2.0 * userdata->radauStepSize;
441 ✗ sData->timeValue = solverInfo->currentTime;
442
443 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "second system time = %e", sData->timeValue);
444
445 /* update step size */
446 ✗ data->simulationInfo->inlineData->dt = userdata->radauStepSize;
447
448 /* evaluate function ODE */
449 ✗ externalInputUpdate(data);
450 ✗ data->callback->input_function(data, threadData);
451 ✗ data->callback->symbolicInlineSystems(data, threadData);
452
453
454 ✗ solverInfo->solverStatsTmp.nStepsTaken += 1;
455 ✗ solverInfo->solverStatsTmp.nCallsODE += 2;
456
457 /* save values in y2 */
458 ✗ memcpy(userdata->y2, sData->realVars, data->modelData->nStates*sizeof(double));
459
460 ✗ userdata->radauStepSize *= 2;
461 }
462 ✗ else if (compiledWithSymSolver == 2) /* compiled with explicit symbolic euler */
463 {
464 /*** do one step with half step size***/
465 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "radauStepSize = %e", userdata->radauStepSize);
466
467 /* update step size */
468 ✗ userdata->radauStepSize /= 2;
469 ✗ data->simulationInfo->inlineData->dt = userdata->radauStepSize;
470
471 /* update algOldVars */
472 ✗ memcpy(data->simulationInfo->inlineData->algOldVars, userdata->radauVars, data->modelData->nStates * sizeof(double));
473
474 /* update time */
475 ✗ sDataOld->timeValue = userdata->radauTime;
476 ✗ solverInfo->currentTime = userdata->radauTime + userdata->radauStepSize;
477 ✗ sData->timeValue = solverInfo->currentTime;
478
479 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "first system time = %e", sData->timeValue);
480
481 /* evaluate function */
482 ✗ externalInputUpdate(data);
483 ✗ data->callback->input_function(data, threadData);
484 ✗ retVal = data->callback->symbolicInlineSystems(data, threadData);
485
486 ✗ if(retVal != 0){
487 return -1;
488 }
489
490 /* save values in y05 */
491 ✗ memcpy(userdata->y05, sData->realVars, data->modelData->nStates*sizeof(double));
492
493 /*** extrapolate values in y1 (= y0 + h * f(y,t)) ***/
494 ✗ for (i=0; i<data->modelData->nStates; i++)
495 {
496 ✗ userdata->y1[i] = 2.0 * userdata->y05[i] - userdata->radauVars[i];
497 }
498
499 /*** do another step with half step size ***/
500 ✗ memcpy(data->simulationInfo->inlineData->algOldVars, userdata->y05, data->modelData->nStates * sizeof(double));
501
502 /* update time */
503 ✗ sDataOld->timeValue = userdata->radauTime + userdata->radauStepSize;
504 ✗ solverInfo->currentTime = userdata->radauTime + 2.0 * userdata->radauStepSize;
505 ✗ sData->timeValue = solverInfo->currentTime;
506
507 ✗ infoStreamPrint(OMC_LOG_SOLVER,0, "second system time = %e", sData->timeValue);
508
509 /* update step size */
510 ✗ data->simulationInfo->inlineData->dt = userdata->radauStepSize;
511
512 /* evaluate function ODE */
513 ✗ externalInputUpdate(data);
514 ✗ data->callback->input_function(data, threadData);
515 ✗ data->callback->symbolicInlineSystems(data, threadData);
516
517
518 ✗ solverInfo->solverStatsTmp.nStepsTaken += 1;
519 ✗ solverInfo->solverStatsTmp.nCallsODE += 2;
520
521 /* save values in y2 */
522 ✗ memcpy(userdata->y2, sData->realVars, data->modelData->nStates*sizeof(double));
523
524 /*** generate solution of higher order via richardson extrapolation */
525 ✗ for (i=0; i<data->modelData->nStates; i++)
526 {
527 ✗ userdata->y1[i] = 2.0 * userdata->y2[i] - userdata->y1[i];
528 }
529
530 ✗ userdata->radauStepSize *= 2;
531
532 }
533
534 return 0;
535
536 }
537