Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/irksco.c
Warning:line 284, column 9
Value stored to 'data' during its initialization is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2014, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
11 * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
14 * ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the OSMC (Open Source Modelica Consortium)
17 * Public License (OSMC-PL) are obtained from OSMC, either from the above
18 * address, from the URLs: http://www.openmodelica.org or
19 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica
20 * distribution. GNU version 3 is obtained from:
21 * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from:
22 * http://www.opensource.org/licenses/BSD-3-Clause.
23 *
24 * This program is distributed WITHOUT ANY WARRANTY; without even the implied
25 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS
26 * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE
27 * CONDITIONS OF OSMC-PL.
28 *
29 */
30
31/*! \file irksco.c
32 * implicit solver for the numerical solution of ordinary differential equations (with step size control)
33 * \author kbalzereit, wbraun
34 */
35
36
37#include <string.h>
38#include <float.h>
39
40#include "simulation/results/simulation_result.h"
41#include "util/omc_error.h"
42#include "util/varinfo.h"
43#include "model_help.h"
44#include "external_input.h"
45#include "newtonIteration.h"
46#include "irksco.h"
47
48int wrapper_fvec_irksco(int* n, double* x, double* f, void* userdata, int fj);
49static int refreshModel(DATA* data, threadData_t *threadData, double* x, double time);
50void irksco_first_step(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo);
51int rk_imp_step(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo, double* y_new);
52
53
54/*! \fn allocateIrksco
55 *
56 * Function allocates memory needed for implicit rk solving methods.
57 *
58 * Integration methods of higher order can be used by increasing userdata->order:
59 * (1) implicit euler method
60 *
61 */
62int allocateIrksco(SOLVER_INFO* solverInfo, int size, int zcSize)
63{
64 DATA_IRKSCO* userdata = (DATA_IRKSCO*) malloc(sizeof(DATA_IRKSCO));
65 solverInfo->solverData = (void*) userdata;
66 userdata->order = 1;
67 userdata->ordersize = 1;
68
69 allocateNewtonData(userdata->ordersize*size, &(userdata->solverData));
70 userdata->firstStep = 1;
71 userdata->y0 = malloc(sizeof(double)*size);
72 userdata->y05= malloc(sizeof(double)*size);
73 userdata->y1 = malloc(sizeof(double)*size);
74 userdata->y2 = malloc(sizeof(double)*size);
75 userdata->der_x0 = malloc(sizeof(double)*size);
76 userdata->radauVarsOld = malloc(sizeof(double)*size);
77 userdata->radauVars = malloc(sizeof(double)*size);
78 userdata->zeroCrossingValues = malloc(sizeof(double)*zcSize);
79 userdata->zeroCrossingValuesOld = malloc(sizeof(double)*zcSize);
80
81 userdata->m = malloc(sizeof(double)*size);
82 userdata->n = malloc(sizeof(double)*size);
83
84 userdata->A = malloc(sizeof(double)*userdata->ordersize*userdata->ordersize);
85 userdata->Ainv = malloc(sizeof(double)*userdata->ordersize*userdata->ordersize);
86 userdata->c = malloc(sizeof(double)*userdata->ordersize);
87 userdata->d = malloc(sizeof(double)*userdata->ordersize);
88
89 /* initialize stats */
90 userdata->stepsDone = 0;
91 userdata->evalFunctionODE = 0;
92 userdata->evalJacobians = 0;
93
94 userdata->radauStepSizeOld = 0;
95 userdata->A[0] = 1;
96 userdata->c[0] = 1;
97 userdata->d[0] = 1;
98
99 return 0;
100}
101
102/*! \fn freeIrksco
103 *
104 * Memory needed for solver is set free.
105 */
106int freeIrksco(SOLVER_INFO* solverInfo)
107{
108 DATA_IRKSCO* userdata = (DATA_IRKSCO*) solverInfo->solverData;
109 freeNewtonData(&(userdata->solverData));
110
111 free(userdata->y0);
112 free(userdata->y05);
113 free(userdata->y1);
114 free(userdata->y2);
115 free(userdata->der_x0);
116 free(userdata->radauVarsOld);
117 free(userdata->radauVars);
118 free(userdata->zeroCrossingValues);
119 free(userdata->zeroCrossingValuesOld);
120
121 return 0;
122}
123
124/*! \fn checkForZeroCrossingsIrksco
125 *
126 * This function checks for ZeroCrossings.
127 */
128int checkForZeroCrossingsIrksco(DATA* data, threadData_t *threadData, DATA_IRKSCO* irkscoData, double *gout)
129{
130 TRACE_PUSH
131
132 /* read input vars */
133 externalInputUpdate(data);
134 data->callback->input_function(data, threadData);
135 /* eval needed equations*/
136 data->callback->function_ZeroCrossingsEquations(data, threadData);
137
138 data->callback->function_ZeroCrossings(data, threadData, gout);
139
140 TRACE_POP
141 return 0;
142}
143
144/*! \fn compareZeroCrossings
145 *
146 * This function compares gout vs. gout_old and return 1,
147 * if they are not equal, otherwise it returns 0,
148 *
149 * \param [ref] [data]
150 * \param [in] [gout]
151 * \param [in] [gout_old]
152 *
153 */
154int compareZeroCrossings(DATA* data, double* gout, double* gout_old)
155{
156 TRACE_PUSH
157 int i;
158
159 for(i=0; i<data->modelData->nZeroCrossings; ++i)
160 if(gout[i] != gout_old[i])
161 return 1;
162
163 TRACE_POP
164 return 0;
165}
166
167/*! \fn rk_imp_step
168 *
169 * function does one implicit euler step with the stepSize given in radauStepSize
170 * function omc_newton is used for solving nonlinear system
171 * results will be saved in y_new
172 *
173 */
174int rk_imp_step(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo, double* y_new)
175{
176 int i, j, n=data->modelData->nStates;
177
178 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
179 modelica_real* stateDer = sData->realVars + data->modelData->nStates;
180 NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo->nonlinearSystemData;
181 DATA_IRKSCO* userdata = (DATA_IRKSCO*)solverInfo->solverData;
182 DATA_NEWTON* solverData = (DATA_NEWTON*) userdata->solverData;
183
184 userdata->data = (void*) data;
185 userdata->threadData = threadData;
186
187 double a,b;
188
189 sData->timeValue = userdata->radauTime + userdata->radauStepSize;
190 solverInfo->currentTime = sData->timeValue;
191
192 solverData->initialized = 1;
193 solverData->numberOfIterations = 0;
194 solverData->numberOfFunctionEvaluations = 0;
195 solverData->n = n*userdata->ordersize;
196
197
198 /* linear extrapolation for start value of newton iteration */
199 for (i=0; i<n; i++)
200 {
201 if (userdata->radauStepSizeOld > 1e-16)
202 {
203 userdata->m[i] = (userdata->radauVars[i] - userdata->radauVarsOld[i]) / userdata->radauStepSizeOld;
204 userdata->n[i] = userdata->radauVars[i] - userdata->radauTime * userdata->m[i];
205 }
206 else
207 {
208 userdata->m[i] = 0;
209 userdata->n[i] = 0;
210 }
211 }
212
213 /* initial guess calculated via linear extrapolation */
214 for (i=0; i<userdata->ordersize; i++)
215 {
216 if (userdata->radauStepSizeOld > 1e-16)
217 {
218 for (j=0; j<n; j++)
219 {
220 solverData->x[i*n+j] = userdata->m[j] * (userdata->radauTimeOld + userdata->c[i] * userdata->radauStepSize )+ userdata->n[j] - userdata->y0[j];
221 }
222 }
223 else
224 {
225 for (j=0; j<n; j++)
226 {
227 solverData->x[i*n+j] = userdata->radauVars[i];
228 }
229 }
230 }
231
232 solverData->newtonStrategy = NEWTON_DAMPED2;
233 _omc_newton(wrapper_fvec_irksco, solverData, (void*)userdata);
234
235 /* if newton solver did not converge, do iteration again but calculate jacobian in every step */
236 if (solverData->info == -1)
237 {
238 for (i=0; i<userdata->ordersize; i++)
239 {
240 for (j=0; j<n; j++)
241 {
242 solverData->x[i*n+j] = userdata->m[j] * (userdata->radauTimeOld + userdata->c[i] * userdata->radauStepSize )+ userdata->n[j] - userdata->y0[j];
243 }
244 }
245 solverData->numberOfIterations = 0;
246 solverData->numberOfFunctionEvaluations = 0;
247 solverData->calculate_jacobian = 1;
248
249 warningStreamPrint(LOG_SOLVER, 0, "nonlinear solver did not converge at time %e, do iteration again with calculating jacobian in every step", solverInfo->currentTime);
250 _omc_newton(wrapper_fvec_irksco, solverData, (void*)userdata);
251
252 solverData->calculate_jacobian = -1;
253 }
254
255 for (j=0; j<n; j++)
256 {
257 y_new[j] = userdata->y0[j];
258 }
259
260 for (i=0; i<userdata->ordersize; i++)
261 {
262 if (userdata->d[i] != 0)
263 {
264 for (j=0; j<n; j++)
265 {
266 y_new[j] += userdata->d[i] * solverData->x[i*n+j];
267 }
268 }
269 }
270
271
272 return 0;
273}
274
275
276/*! \fn wrapper_fvec_irksco
277 *
278 * calculate function values or jacobian matrix
279 * fj = 1 ==> calculate function values
280 * fj = 0 ==> calculate jacobian matrix
281 */
282int wrapper_fvec_irksco(int* n, double* x, double* fvec, void* userdata, int fj)
283{
284 DATA* data = ((DATA_IRKSCO*)userdata)->data;
Value stored to 'data' during its initialization is never read
285 threadData_t* threadData = ((DATA_IRKSCO*)userdata)->threadData;
286 DATA_NEWTON* solverData = (DATA_NEWTON*)((DATA_IRKSCO*)userdata)->solverData;
287 if (fj)
288 {
289 int i, j, k;
290 DATA_IRKSCO* irkscoData = (DATA_IRKSCO*) userdata;
291 DATA* data = irkscoData->data;
292 int n0 = (*n)/irkscoData->ordersize;
293 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
294 modelica_real* stateDer = sData->realVars + data->modelData->nStates;
295
296 ((DATA_IRKSCO*)userdata)->evalFunctionODE++;
297
298 for (k=0; k < irkscoData->ordersize; k++)
299 {
300 for (j=0; j<n0; j++)
301 {
302 fvec[k*n0+j] = x[k*n0+j];
303 }
304 }
305
306 for (i=0; i < irkscoData->ordersize; i++)
307 {
308 sData->timeValue = irkscoData->radauTimeOld + irkscoData->c[i] * irkscoData->radauStepSize;
309
310 for (j=0; j < n0; j++)
311 {
312 sData->realVars[j] = irkscoData->y0[j] + x[n0*i+j];
313 }
314
315
316 externalInputUpdate(data);
317 data->callback->input_function(data, threadData);
318 data->callback->functionODE(data, threadData);
319
320 for (k=0; k < irkscoData->ordersize; k++)
321 {
322 for (j=0; j<n0; j++)
323 {
324 fvec[k*n0+j] -= irkscoData->A[i*irkscoData->ordersize+k] * irkscoData->radauStepSize * stateDer[j];
325 }
326 }
327
328 }
329 }
330 else
331 {
332 double delta_h = sqrt(solverData->epsfcn);
333 double delta_hh;
334 double xsave;
335
336 int i,j,l;
337
338 /* profiling */
339 rt_tick(SIM_TIMER_JACOBIAN5);
340
341 ((DATA_IRKSCO*)userdata)->evalJacobians++;
342
343 for(i = 0; i < *n; i++)
344 {
345 delta_hh = fmax(delta_h * fmax(fabs(x[i]), fabs(fvec[i])), delta_h);
346 delta_hh = ((fvec[i] >= 0) ? delta_hh : -delta_hh);
347 delta_hh = x[i] + delta_hh - x[i];
348 xsave = x[i];
349 x[i] += delta_hh;
350 delta_hh = 1. / delta_hh;
351
352 wrapper_fvec_irksco(n, x, solverData->rwork, userdata, 1);
353 solverData->nfev++;
354
355 for(j = 0; j < *n; j++)
356 {
357 l = i * *n + j;
358 solverData->fjac[l] = (solverData->rwork[j] - fvec[j]) * delta_hh;
359 }
360 x[i] = xsave;
361 }
362
363 /* profiling */
364 rt_accumulate(SIM_TIMER_JACOBIAN5);
365 }
366 return 0;
367}
368
369/*! \fn refreshModel
370 *
371 * function updates values in sData->realVars
372 *
373 * used for solver 'irksco'
374 */
375static
376int refreshModel(DATA* data, threadData_t *threadData, double* x, double time)
377{
378 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
379
380 memcpy(sData->realVars, x, sizeof(double)*data->modelData->nStates);
381 sData->timeValue = time;
382 /* read input vars */
383 externalInputUpdate(data);
384 data->callback->input_function(data, threadData);
385 data->callback->functionODE(data, threadData);
386
387 return 0;
388}
389
390/*! \fn irksco_midpoint_rule
391 *
392 * function does one integration step and calculates
393 * next step size by the implicit midpoint rule
394 *
395 * used for solver 'irksco'
396 */
397int irksco_midpoint_rule(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
398{
399 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
400 SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
401 modelica_real* stateDer = sData->realVars + data->modelData->nStates;
402 DATA_IRKSCO* userdata = (DATA_IRKSCO*)solverInfo->solverData;
403 DATA_NEWTON* solverData = (DATA_NEWTON*)userdata->solverData;
404
405 double sc, err, a, b, diff;
406 double Atol = data->simulationInfo->tolerance, Rtol = data->simulationInfo->tolerance;
407 int i;
408 double fac = 0.9;
409 double facmax = 3.5;
410 double facmin = 0.3;
411 double saveTime = sDataOld->timeValue;
412 double targetTime;
413
414
415
416 /* Calculate steps until targetTime is reached */
417 if (solverInfo->integratorSteps)
418 {
419 if (data->simulationInfo->nextSampleEvent < data->simulationInfo->stopTime)
420 {
421 targetTime = data->simulationInfo->nextSampleEvent;
422 }
423 else
424 {
425 targetTime = data->simulationInfo->stopTime;
426 }
427 }
428 else
429 {
430 targetTime = sDataOld->timeValue + solverInfo->currentStepSize;
431 }
432
433 if (userdata->firstStep || solverInfo->didEventStep == 1)
434 {
435 irksco_first_step(data, threadData, solverInfo);
436 userdata->radauStepSizeOld = 0;
437 }
438
439 memcpy(userdata->y0, sDataOld->realVars, data->modelData->nStates*sizeof(double));
440
441 while (userdata->radauTime < targetTime)
442 {
443 infoStreamPrint(LOG_SOLVER, 1, "new step to %f -> targetTime: %f", userdata->radauTime, targetTime);
444
445 do
446 {
447 /*** do one step with original step size ***/
448 /* set y0 */
449 memcpy(userdata->y0, userdata->radauVars, data->modelData->nStates*sizeof(double));
450
451 /* calculate jacobian once for the first iteration */
452 if (userdata->stepsDone == 0)
453 solverData->calculate_jacobian = 0;
454
455 /* solve nonlinear system */
456 rk_imp_step(data, threadData, solverInfo, userdata->y05);
457
458 /* extrapolate values in y1 */
459 for (i=0; i<data->modelData->nStates; i++)
460 {
461 userdata->y1[i] = 2.0*userdata->y05[i] - userdata->radauVars[i];
462 }
463
464 /*** do another step with original step size ***/
465
466 /* update y0 */
467 memcpy(userdata->y0, userdata->y05, data->modelData->nStates*sizeof(double));
468
469 /* update time */
470 userdata->radauTime += userdata->radauStepSize;
471
472 /* do not calculate jacobian again */
473 solverData->calculate_jacobian = -1;
474
475 /* solve nonlinear system */
476 rk_imp_step(data, threadData, solverInfo, userdata->y2);
477
478 /* reset time */
479 userdata->radauTime -= userdata->radauStepSize;
480
481
482 /*** calculate error ***/
483 for (i=0, err=0.0; i<data->modelData->nStates; i++)
484 {
485 sc = Atol + fmax(fabs(userdata->y2[i]),fabs(userdata->y1[i]))*Rtol;
486 diff = userdata->y2[i]-userdata->y1[i];
487 err += (diff*diff)/(sc*sc);
488 }
489
490 err /= data->modelData->nStates;
491 err = sqrt(err);
492
493 userdata->stepsDone += 1;
494 /* debug
495 infoStreamPrint(LOG_SOLVER, 0, "err = %e", err);
496 infoStreamPrint(LOG_SOLVER, 0, "min(facmax, max(facmin, fac*sqrt(1/err))) = %e", fmin(facmax, fmax(facmin, fac*sqrt(1.0/err))));
497 */
498 /* update step size */
499 userdata->radauStepSizeOld = 2.0*userdata->radauStepSize;
500 userdata->radauStepSize *= fmin(facmax, fmax(facmin, fac*sqrt(1.0/err)));
501
502 if (isnan(userdata->radauStepSize)(sizeof ((userdata->radauStepSize)) == sizeof (float) ? __isnanf
(userdata->radauStepSize) : sizeof ((userdata->radauStepSize
)) == sizeof (double) ? __isnan (userdata->radauStepSize) :
__isnanl (userdata->radauStepSize))
)
503 {
504 userdata->radauStepSize = 1e-6;
505 }
506
507 } while (err > 1.0 );
508
509 userdata->radauTimeOld = userdata->radauTime;
510
511 userdata->radauTime += userdata->radauStepSizeOld;
512
513 memcpy(userdata->radauVarsOld, userdata->radauVars, data->modelData->nStates*sizeof(double));
514 memcpy(userdata->radauVars, userdata->y2, data->modelData->nStates*sizeof(double));
515
516 /* emit step, if integratorSteps is selected */
517 if (solverInfo->integratorSteps)
518 {
519 sData->timeValue = userdata->radauTime;
520 memcpy(sData->realVars, userdata->radauVars, data->modelData->nStates*sizeof(double));
521 /*
522 * to emit consistent value we need to update the whole
523 * continuous system with algebraic variables.
524 */
525 data->callback->updateContinuousSystem(data, threadData);
526 sim_result.emit(&sim_result, data, threadData);
527 }
528 messageClose(LOG_SOLVER);
529 }
530
531 if (!solverInfo->integratorSteps)
532 {
533 solverInfo->currentTime = sDataOld->timeValue + solverInfo->currentStepSize;
534 sData->timeValue = solverInfo->currentTime;
535 /* linear interpolation */
536 for (i=0; i<data->modelData->nStates; i++)
537 {
538 a = (userdata->radauVars[i] - userdata->radauVarsOld[i]) / userdata->radauStepSizeOld;
539 b = userdata->radauVars[i] - userdata->radauTime * a;
540 sData->realVars[i] = a * sData->timeValue + b;
541 }
542 }else{
543 solverInfo->currentTime = userdata->radauTime;
544 }
545
546 /* if a state event occurs than no sample event does need to be activated */
547 if (data->simulationInfo->sampleActivated && solverInfo->currentTime < data->simulationInfo->nextSampleEvent)
548 {
549 data->simulationInfo->sampleActivated = 0;
550 }
551
552
553 if(ACTIVE_STREAM(LOG_SOLVER)(useStream[LOG_SOLVER]))
554 {
555 infoStreamPrint(LOG_SOLVER, 1, "irksco call statistics: ");
556 infoStreamPrint(LOG_SOLVER, 0, "current time value: %0.4g", solverInfo->currentTime);
557 infoStreamPrint(LOG_SOLVER, 0, "current integration time value: %0.4g", userdata->radauTime);
558 infoStreamPrint(LOG_SOLVER, 0, "step size H to be attempted on next step: %0.4g", userdata->radauStepSize);
559 infoStreamPrint(LOG_SOLVER, 0, "number of steps taken so far: %d", userdata->stepsDone);
560 infoStreamPrint(LOG_SOLVER, 0, "number of calls of functionODE() : %d", userdata->evalFunctionODE);
561 infoStreamPrint(LOG_SOLVER, 0, "number of calculation of jacobian : %d", userdata->evalJacobians);
562 messageClose(LOG_SOLVER);
563 }
564
565 /* write stats */
566 solverInfo->solverStatsTmp[0] = userdata->stepsDone;
567 solverInfo->solverStatsTmp[1] = userdata->evalFunctionODE;
568 solverInfo->solverStatsTmp[2] = userdata->evalJacobians;
569
570 infoStreamPrint(LOG_SOLVER, 0, "Finished irksco step.");
571
572 return 0;
573}
574
575/*! \fn irksco_first_step
576 *
577 * function initializes values and calculates
578 * initial step size at the beginning or after an event
579 *
580 */
581void irksco_first_step(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
582{
583 SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
584 SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
585 DATA_IRKSCO* userdata = (DATA_IRKSCO*)solverInfo->solverData;
586 const int n = data->modelData->nStates;
587 modelica_real* stateDer = sData->realVars + data->modelData->nStates;
588 double sc, d, d0 = 0.0, d1 = 0.0, d2 = 0.0, h0, h1, delta_ti, infNorm, sum = 0;
589 double Atol = 1e-6, Rtol = 1e-3;
590
591 int i,j;
592
593 /* initialize radau values */
594 for (i=0; i<data->modelData->nStates; i++)
595 {
596 userdata->radauVars[i] = sData->realVars[i];
597 userdata->radauVarsOld[i] = sDataOld->realVars[i];
598 }
599
600 userdata->radauTime = sDataOld->timeValue;
601 userdata->radauTimeOld = sDataOld->timeValue;
602
603 userdata->firstStep = 0;
604 solverInfo->didEventStep = 0;
605
606
607 /* calculate starting step size 1st Version */
608
609 refreshModel(data, threadData, sDataOld->realVars, sDataOld->timeValue);
610
611
612 for (i=0; i<data->modelData->nStates; i++)
613 {
614 sc = Atol + fabs(sDataOld->realVars[i])*Rtol;
615 d0 += ((sDataOld->realVars[i] * sDataOld->realVars[i])/(sc*sc));
616 d1 += ((stateDer[i] * stateDer[i]) / (sc*sc));
617 }
618 d0 /= data->modelData->nStates;
619 d1 /= data->modelData->nStates;
620
621 d0 = sqrt(d0);
622 d1 = sqrt(d1);
623
624
625 for (i=0; i<data->modelData->nStates; i++)
626 {
627 userdata->der_x0[i] = stateDer[i];
628 }
629
630 if (d0 < 1e-5 || d1 < 1e-5)
631 {
632 h0 = 1e-6;
633 }
634 else
635 {
636 h0 = 0.01 * d0/d1;
637 }
638
639
640 for (i=0; i<data->modelData->nStates; i++)
641 {
642 sData->realVars[i] = userdata->radauVars[i] + stateDer[i] * h0;
643 }
644 sData->timeValue += h0;
645
646 externalInputUpdate(data);
647 data->callback->input_function(data, threadData);
648 data->callback->functionODE(data, threadData);
649
650
651 for (i=0; i<data->modelData->nStates; i++)
652 {
653 sc = Atol + fabs(userdata->radauVars[i])*Rtol;
654 d2 += ((stateDer[i]-userdata->der_x0[i])*(stateDer[i]-userdata->der_x0[i])/(sc*sc));
655 }
656
657 d2 /= h0;
658 d2 = sqrt(d2);
659
660
661 d = fmax(d1,d2);
662
663 if (d > 1e-15)
664 {
665 h1 = sqrt(0.01/d);
666 }
667 else
668 {
669 h1 = fmax(1e-6, h0*1e-3);
670 }
671
672 userdata->radauStepSize = 0.5*fmin(100*h0,h1);
673
674 /* end calculation new step size */
675
676 infoStreamPrint(LOG_SOLVER, 0, "initial step size = %e", userdata->radauStepSize);
677}