| 1 |
/* ASCEND modelling environment |
| 2 |
Copyright (C) 2007 Carnegie Mellon University |
| 3 |
|
| 4 |
This program is free software; you can redistribute it and/or modify |
| 5 |
it under the terms of the GNU General Public License as published by |
| 6 |
the Free Software Foundation; either version 2, or (at your option) |
| 7 |
any later version. |
| 8 |
|
| 9 |
This program is distributed in the hope that it will be useful, |
| 10 |
but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 11 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 12 |
GNU General Public License for more details. |
| 13 |
|
| 14 |
You should have received a copy of the GNU General Public License |
| 15 |
along with this program; if not, write to the Free Software |
| 16 |
Foundation, Inc., 59 Temple Place - Suite 330, |
| 17 |
Boston, MA 02111-1307, USA. |
| 18 |
*//** @file |
| 19 |
DOPRI5 Runge-Kutta integrator |
| 20 |
|
| 21 |
Based on the implementation of LSODE integrator, but adapted to |
| 22 |
an explicit one-step method. |
| 23 |
*//* |
| 24 |
by John Pye, May 2007. |
| 25 |
*/ |
| 26 |
|
| 27 |
#include <utilities/ascConfig.h> |
| 28 |
#include <utilities/error.h> |
| 29 |
#include <general/ospath.h> |
| 30 |
#include <integrator/integrator.h> |
| 31 |
|
| 32 |
#define INTEG_DOPRI5 5 |
| 33 |
|
| 34 |
IntegratorCreateFn integrator_dopri5_create; |
| 35 |
IntegratorParamsDefaultFn integrator_dopri5_params_default; |
| 36 |
IntegratorSolveFn integrator_dopri5_solve; |
| 37 |
IntegratorFreeFn integrator_dopri5_free; |
| 38 |
IntegratorWriteMatrixFn integrator_dopri5_write_matrix; |
| 39 |
|
| 40 |
const IntegratorInternals integrator_lsode_internals; |
| 41 |
|
| 42 |
const IntegratorInternals integrator_dopri5_internals = |
| 43 |
{ |
| 44 |
integrator_dopri5_create |
| 45 |
,integrator_dopri5_params_default |
| 46 |
,integrator_analyse_ode /* note, this routine is back in integrator.c */ |
| 47 |
,integrator_dopri5_solve |
| 48 |
,integrator_dopri5_write_matrix |
| 49 |
,NULL /* debugfn */ |
| 50 |
,integrator_dopri5_free |
| 51 |
,INTEG_DOPRI5 |
| 52 |
,"DOPRI5" |
| 53 |
}; |
| 54 |
|
| 55 |
extern ASC_EXPORT int dopri5_register(void) |
| 56 |
{ |
| 57 |
CONSOLE_DEBUG("DOPRI5"); |
| 58 |
return 0; |
| 59 |
} |
| 60 |
|
| 61 |
typedef struct IntegratorDopri5DataStruct |
| 62 |
{ |
| 63 |
long n_eqns; /**< dimension of state vector */ |
| 64 |
int *input_indices; /**< vector of state vars indexes */ |
| 65 |
int *output_indices; /**< vector of derivative var indexes */ |
| 66 |
struct var_variable **y_vars; /**< NULL-terminated list of states vars */ |
| 67 |
struct var_variable **ydot_vars; /**< NULL-terminated list of derivative vars*/ |
| 68 |
struct rel_relation **rlist; /**< NULL-terminated list of relevant rels |
| 69 |
to be differentiated */ |
| 70 |
|
| 71 |
char stop; /* stop requested? */ |
| 72 |
int partitioned; /* partioned func evals or not */ |
| 73 |
|
| 74 |
clock_t lastwrite; /* time of last call to the reporter 'write' function */ |
| 75 |
} |
| 76 |
IntegratorDopri5Data; |
| 77 |
|
| 78 |
/*------------------------------------------------------------------------------ |
| 79 |
PARAMETERS |
| 80 |
*/ |
| 81 |
|
| 82 |
enum ida_parameters{ |
| 83 |
DOPRI5_PARAMS_ |
| 84 |
,DOPRI5_PARAMS_RTOL |
| 85 |
,DOPRI5_PARAMS_ATOL |
| 86 |
,DOPRI5_PARAMS_TOLVECT |
| 87 |
|
| 88 |
#if 0 |
| 89 |
// more parameters for adding later... |
| 90 |
SolTrait *solout, /* function providing the numerical solution during integration */ |
| 91 |
int iout, /* switch for calling solout */ |
| 92 |
|
| 93 |
double uround, /* rounding unit */ |
| 94 |
double safe, /* safety factor */ |
| 95 |
double fac1, /* parameters for step size selection */ |
| 96 |
double fac2, |
| 97 |
double beta, /* for stabilized step size control */ |
| 98 |
|
| 99 |
double hmax, /* maximal step size */ |
| 100 |
double h, /* initial step size */ |
| 101 |
long nmax, /* maximal number of allowed steps */ |
| 102 |
|
| 103 |
int meth, /* switch for the choice of the coefficients */ |
| 104 |
long nstiff, /* test for stiffness */ |
| 105 |
#endif |
| 106 |
|
| 107 |
,DOPRI5_PARAMS_SIZE |
| 108 |
}; |
| 109 |
|
| 110 |
/** |
| 111 |
Here the full set of parameters is defined, along with upper/lower bounds, |
| 112 |
etc. The values are stuck into the blsys->params structure. |
| 113 |
|
| 114 |
@return 0 on success |
| 115 |
*/ |
| 116 |
int integrator_dopri5_params_default(IntegratorSystem *blsys){ |
| 117 |
|
| 118 |
asc_assert(blsys!=NULL); |
| 119 |
asc_assert(blsys->engine==INTEG_DOPRI5); |
| 120 |
slv_parameters_t *p; |
| 121 |
p = &(blsys->params); |
| 122 |
|
| 123 |
slv_destroy_parms(p); |
| 124 |
|
| 125 |
if(p->parms==NULL) { |
| 126 |
p->parms = ASC_NEW_ARRAY(struct slv_parameter, LSODE_PARAMS_SIZE); |
| 127 |
if(p->parms==NULL)return -1; |
| 128 |
p->dynamic_parms = 1; |
| 129 |
} else { |
| 130 |
asc_assert(p->num_parms == LSODE_PARAMS_SIZE); |
| 131 |
} |
| 132 |
|
| 133 |
/* reset the number of parameters to zero so that we can check it at the end */ |
| 134 |
p->num_parms = 0; |
| 135 |
|
| 136 |
slv_param_bool(p,DOPRI5_PARAM_TOLVECT |
| 137 |
,(SlvParameterInitBool) { |
| 138 |
{"tolvect" |
| 139 |
,"Use per-variable tolerances?",1 |
| 140 |
,"If TRUE, values of 'ode_rtol' and 'ode_atol' are taken from your" |
| 141 |
" model and used in the integration. If FALSE, scalar values are" |
| 142 |
" used (see 'rtol' and 'atol') and shared by all variables." |
| 143 |
} |
| 144 |
, TRUE |
| 145 |
} |
| 146 |
); |
| 147 |
|
| 148 |
slv_param_real(p,DOPRI5_PARAM_ATOL |
| 149 |
,(SlvParameterInitReal) { |
| 150 |
{"atol" |
| 151 |
,"Scalar absolute error tolerance",1 |
| 152 |
,"Default value of the scalar absolute error tolerance (for cases" |
| 153 |
" where not specified in oda_atol var property. See 'dopri5.h' for" |
| 154 |
" details" |
| 155 |
} |
| 156 |
, 1e-7, 1e-15, 1e10 |
| 157 |
} |
| 158 |
); |
| 159 |
|
| 160 |
slv_param_real(p,DOPRI5_PARAM_RTOL |
| 161 |
,(SlvParameterInitReal) { |
| 162 |
{"rtol" |
| 163 |
,"Scalar relative error tolerance",1 |
| 164 |
,"Default value of the scalar relative error tolerance (for cases" |
| 165 |
" where not specified in oda_rtol var property. See 'dopri5.h' for" |
| 166 |
" details" |
| 167 |
} |
| 168 |
, 1e-7, 1e-15, 1 |
| 169 |
} |
| 170 |
); |
| 171 |
|
| 172 |
/* |
| 173 |
|
| 174 |
slv_param_bool(p,DOPRI5_PARAMS_DENSEREPORTING |
| 175 |
,(SlvParameterInitBool) { |
| 176 |
{"densereporting" |
| 177 |
,"Dense reporting?",1 |
| 178 |
,"If TRUE, output will be made at every sub-timestep" |
| 179 |
" during integration. If false, output is only made |
| 180 |
" according to the samplelist." |
| 181 |
} |
| 182 |
, FALSE |
| 183 |
} |
| 184 |
); |
| 185 |
|
| 186 |
|
| 187 |
slv_param_char(p,LSODE_PARAM_METH |
| 188 |
,(SlvParameterInitChar) { |
| 189 |
{"meth" |
| 190 |
,"Integration method",1 |
| 191 |
,"AM=Adams-Moulton method (for non-stiff problems), BDF=Backwards" |
| 192 |
" Difference Formular (for stiff problems). See 'Description and" |
| 193 |
" Use of LSODE', section 3.1." |
| 194 |
} |
| 195 |
, "BDF" |
| 196 |
} |
| 197 |
, (char *[]) {"AM","BDF",NULL |
| 198 |
} |
| 199 |
); |
| 200 |
|
| 201 |
slv_param_int(p,LSODE_PARAM_MITER |
| 202 |
,(SlvParameterInitInt) { |
| 203 |
{"miter" |
| 204 |
,"Corrector iteration technique", 1 |
| 205 |
,"0=Functional iteration, 1=Modified Newton iteration with user-" |
| 206 |
"supplied analytical Jacobian, 2=Modified Newton iteration with" |
| 207 |
" internally-generated numerical Jacobian, 3=Modified Jacobi-Newton" |
| 208 |
" iteration with internally generated numerical Jacobian. See " |
| 209 |
" 'Description and Use of LSODE', section 3.1. Note that not all" |
| 210 |
" methods described there are available via ASCEND." |
| 211 |
} |
| 212 |
, 1, 0, 3 |
| 213 |
} |
| 214 |
); |
| 215 |
|
| 216 |
slv_param_int(p,LSODE_PARAM_MAXORD |
| 217 |
,(SlvParameterInitInt) { |
| 218 |
{"maxord" |
| 219 |
,"Maximum method order", 1 |
| 220 |
,"See 'Description and Use of LSODE', p 92 and p 8. Limits <=12 for BDF" |
| 221 |
" and <=5 for AM. Higher values are reduced automatically. See notes on" |
| 222 |
" page 92 regarding oscillatory systems." |
| 223 |
} |
| 224 |
, 12, 1, 12 |
| 225 |
} |
| 226 |
); |
| 227 |
|
| 228 |
slv_param_bool(p,LSODE_PARAM_TIMING |
| 229 |
,(SlvParameterInitBool) { |
| 230 |
{"timing" |
| 231 |
,"Output timing statistics?",1 |
| 232 |
,"If TRUE, additional timing statistics will be output to the" |
| 233 |
" console during integration." |
| 234 |
} |
| 235 |
, TRUE |
| 236 |
} |
| 237 |
); |
| 238 |
*/ |
| 239 |
|
| 240 |
asc_assert(p->num_parms == DOPRI5_PARAMS_SIZE); |
| 241 |
return 0; |
| 242 |
} |
| 243 |
|
| 244 |
|
| 245 |
/** |
| 246 |
Allocates, fills, and returns the rtol vector based on LSODE |
| 247 |
|
| 248 |
State variables missing child ode_rtol will be defaulted to RTOL |
| 249 |
|
| 250 |
@param is_r TRUE if we want RTOL, FALSE if we want ATOL. |
| 251 |
*/ |
| 252 |
static double *dopri5_get_artol(IntegratorSystem *blsys, int is_r, int tolvect) { |
| 253 |
|
| 254 |
struct Instance *toli; |
| 255 |
double tolval, *tol; |
| 256 |
int i,len; |
| 257 |
symchar *tolname; |
| 258 |
|
| 259 |
len = blsys->n_y; |
| 260 |
|
| 261 |
if(!tolvect){ |
| 262 |
|
| 263 |
// single tol for all vars |
| 264 |
tolval = SLV_PARAM_REAL(&(blsys->params),DOPRI5_PARAM_RTOL); |
| 265 |
CONSOLE_DEBUG("Using RTOL = %f for all vars", tolval); |
| 266 |
|
| 267 |
tol = ASC_NEW(double); |
| 268 |
if(tol == NULL){ |
| 269 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory"); |
| 270 |
return tol; |
| 271 |
} |
| 272 |
|
| 273 |
*tol = tolval; |
| 274 |
return tolval; |
| 275 |
|
| 276 |
}else{ |
| 277 |
tol = ASC_NEW_ARRAY(double, blsys->n_y+1); |
| 278 |
if (tol == NULL) { |
| 279 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory"); |
| 280 |
return tol; |
| 281 |
} |
| 282 |
|
| 283 |
tolval = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_RTOL) |
| 284 |
|
| 285 |
InitTolNames(); // from where? |
| 286 |
if(is_r)tolname = STATERTOL; |
| 287 |
else tolname = STATEATOL; |
| 288 |
|
| 289 |
for(i=0; i<len; i++){ |
| 290 |
toli = ChildByChar(var_instance(blsys->y[i]),tolname); |
| 291 |
if(toli == NULL || !AtomAssigned(toli)){ |
| 292 |
tol[i] = tolval; |
| 293 |
ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming value %3g" |
| 294 |
"for '%s' child undefined for state variable %ld." |
| 295 |
,tol[i], SCP(tolname), blsys->y_id[i] |
| 296 |
); |
| 297 |
}else{ |
| 298 |
tol[i] = RealAtomValue(toli); |
| 299 |
} |
| 300 |
} |
| 301 |
} |
| 302 |
tol[len] = SLV_PARAM_REAL(&(blsys->params),LSODE_PARAM_RTOL); |
| 303 |
return tol; |
| 304 |
} |
| 305 |
|
| 306 |
/*------------------------------------------------------------------------------ |
| 307 |
FUNCTION EVALUATION |
| 308 |
*/ |
| 309 |
|
| 310 |
FcnEqDiff integrator_dopri5_fex; |
| 311 |
|
| 312 |
/*------------------------------------------------------------------------------ |
| 313 |
REPORTING |
| 314 |
*/ |
| 315 |
|
| 316 |
SolTrait integrator_dopri5_reporter; |
| 317 |
|
| 318 |
/*------------------------------------------------------------------------------ |
| 319 |
SOLVE |
| 320 |
*/ |
| 321 |
int integrator_dopri5_solve(IntegratorSystem *blsys |
| 322 |
, unsigned long start_index, unsigned long finish_index |
| 323 |
){ |
| 324 |
|
| 325 |
slv_status_t status; |
| 326 |
slv_parameters_t params; |
| 327 |
IntegratorLsodeData *d; |
| 328 |
|
| 329 |
double x[2]; |
| 330 |
double xend,xprev; |
| 331 |
unsigned long nsamples, neq; |
| 332 |
long nobs; |
| 333 |
int itol, itask, mf, lrw, liw; |
| 334 |
unsigned long index; |
| 335 |
int istate, iopt; |
| 336 |
double * rwork; |
| 337 |
int * iwork; |
| 338 |
double *y, *atoler, *rtoler, *obs, *dydx; |
| 339 |
int my_neq; |
| 340 |
int reporterstatus; |
| 341 |
const char *method; /* Table 3.1 in D&UoLSODE */ |
| 342 |
int miter; /* Table 3.2 in D&UoLSODE */ |
| 343 |
int maxord; /* page 92 in D&UoLSODE */ |
| 344 |
|
| 345 |
d = (IntegratorLsodeData *)(blsys->enginedata); |
| 346 |
|
| 347 |
/* the numer of equations must be equal to blsys->n_y, the number of states */ |
| 348 |
d->n_eqns = blsys->n_y; |
| 349 |
assert(d->n_eqns>0); |
| 350 |
|
| 351 |
d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns); |
| 352 |
d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns); |
| 353 |
|
| 354 |
d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1); |
| 355 |
d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1); |
| 356 |
|
| 357 |
/* set up the NLA solver here */ |
| 358 |
|
| 359 |
/* here we assume integrators.c is in charge of dynamic loading */ |
| 360 |
|
| 361 |
/* set up parameteers for sending to DOPRI5 */ |
| 362 |
|
| 363 |
#if 0 |
| 364 |
method = SLV_PARAM_CHAR(&(blsys->params),LSODE_PARAM_METH); |
| 365 |
miter = SLV_PARAM_INT(&(blsys->params),LSODE_PARAM_MITER); |
| 366 |
maxord = SLV_PARAM_INT(&(blsys->params),LSODE_PARAM_MAXORD); |
| 367 |
if(miter < 0 || miter > 3) { |
| 368 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unacceptable value '%d' of parameter 'miter'",miter); |
| 369 |
return 5; |
| 370 |
} |
| 371 |
if(strcmp(method,"BDF")==0) { |
| 372 |
CONSOLE_DEBUG("method = BDF"); |
| 373 |
mf = 20 + miter; |
| 374 |
if(maxord > 5) { |
| 375 |
maxord = 5; |
| 376 |
CONSOLE_DEBUG("MAXORD reduced to 5 for BDF"); |
| 377 |
} |
| 378 |
} else if(strcmp(method,"AM")==0) { |
| 379 |
CONSOLE_DEBUG("method = AM"); |
| 380 |
if(maxord > 12) { |
| 381 |
maxord = 12; |
| 382 |
CONSOLE_DEBUG("MAXORD reduced to 12 for AM"); |
| 383 |
} |
| 384 |
mf = 10 + miter; |
| 385 |
} else { |
| 386 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unacceptable value '%d' of parameter 'meth'",method); |
| 387 |
return 5; |
| 388 |
} |
| 389 |
|
| 390 |
CONSOLE_DEBUG("MF = %d",mf); |
| 391 |
#endif |
| 392 |
|
| 393 |
nsamples = integrator_getnsamples(blsys); |
| 394 |
if (nsamples <2) { |
| 395 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system has no end sample time defined."); |
| 396 |
d->status = lsode_nok; |
| 397 |
return 3; |
| 398 |
} |
| 399 |
neq = blsys->n_y; |
| 400 |
nobs = blsys->n_obs; |
| 401 |
|
| 402 |
#if 0 |
| 403 |
unsigned nrdens, /* number of components for which dense outpout is required */ |
| 404 |
unsigned* icont, /* indexes of components for which dense output is required, >= nrdens */ |
| 405 |
unsigned licont /* declared length of icon */ |
| 406 |
#endif |
| 407 |
|
| 408 |
int iout = 1; /* SLV_PARAM_BOOL(&(blsys->params),DOPRI5_PARAM_DENSEREPORTING) */ |
| 409 |
|
| 410 |
int tolvect = SLV_PARAM_BOOL(&(blsys->params),DOPRI5_PARAM_TOLVECT); |
| 411 |
|
| 412 |
/* samplelist_debug(blsys->samples); */ |
| 413 |
|
| 414 |
#if 0 |
| 415 |
/* x[0] = integrator_get_t(blsys); */ |
| 416 |
x[0] = integrator_getsample(blsys, 0); |
| 417 |
x[1] = x[0]-1; /* make sure we don't start with wierd x[1] */ |
| 418 |
|
| 419 |
/* RWORK memory requirements: see D&UoLSODE p 82 */ |
| 420 |
switch(mf) { |
| 421 |
case 10: |
| 422 |
case 20: |
| 423 |
lrw = 20 + neq * (maxord + 1) + 3 * neq; |
| 424 |
break; |
| 425 |
case 11: |
| 426 |
case 12: |
| 427 |
case 21: |
| 428 |
case 22: |
| 429 |
lrw = 22 + neq * (maxord + 1) + 3 * neq + neq*neq; |
| 430 |
break; |
| 431 |
case 13: |
| 432 |
case 23: |
| 433 |
lrw = 22 + neq * (maxord + 1) + 4 * neq; |
| 434 |
break; |
| 435 |
default: |
| 436 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unknown size requirements for this value of 'mf'"); |
| 437 |
return 4; |
| 438 |
} |
| 439 |
|
| 440 |
rwork = ASC_NEW_ARRAY_CLEAR(double, lrw+1); |
| 441 |
liw = 20 + neq; |
| 442 |
|
| 443 |
#endif |
| 444 |
|
| 445 |
iwork = ASC_NEW_ARRAY_CLEAR(int, liw+1); |
| 446 |
y = integrator_get_y(blsys, NULL); |
| 447 |
|
| 448 |
rtoler = lsode_get_artol(blsys,1,tolvect); |
| 449 |
atoler = lsode_get_artol(blsys,0,tolvect); |
| 450 |
|
| 451 |
obs = integrator_get_observations(blsys, NULL); |
| 452 |
|
| 453 |
#if 0 |
| 454 |
dydx = ASC_NEW_ARRAY_CLEAR(double, neq+1); |
| 455 |
if(!y || !obs || !atoler || !rtoler || !rwork || !iwork || !dydx) { |
| 456 |
lsode_free_mem(y,rtoler,atoler,rwork,iwork,obs,dydx); |
| 457 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory for lsode."); |
| 458 |
d->status = lsode_nok; |
| 459 |
return 4; |
| 460 |
} |
| 461 |
#endif |
| 462 |
|
| 463 |
/* |
| 464 |
Prepare args and call lsode. |
| 465 |
*/ |
| 466 |
itol = 4; |
| 467 |
itask = 1; |
| 468 |
istate = 1; |
| 469 |
iopt = 1; |
| 470 |
rwork[4] = integrator_get_stepzero(blsys); |
| 471 |
rwork[5] = integrator_get_maxstep(blsys); |
| 472 |
rwork[6] = integrator_get_minstep(blsys); |
| 473 |
iwork[5] = integrator_get_maxsubsteps(blsys); |
| 474 |
iwork[4] = maxord; |
| 475 |
CONSOLE_DEBUG("MAXORD = %d",maxord); |
| 476 |
|
| 477 |
if(x[0] > integrator_getsample(blsys, 2)) { |
| 478 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid initialisation time: exceeds second timestep value"); |
| 479 |
return 5; |
| 480 |
} |
| 481 |
|
| 482 |
/* put the values from derivative system into the record */ |
| 483 |
integrator_setsample(blsys, start_index, x[0]); |
| 484 |
|
| 485 |
integrator_output_init(blsys); |
| 486 |
|
| 487 |
/* -- store the initial values of all the stuff */ |
| 488 |
integrator_output_write(blsys); |
| 489 |
integrator_output_write_obs(blsys); |
| 490 |
|
| 491 |
my_neq = (int)neq; |
| 492 |
|
| 493 |
/* |
| 494 |
First time entering lsode, x is input. After that, |
| 495 |
lsode uses x as output (y output is y(x)). To drive |
| 496 |
the loop ahead in time, all we need to do is keep upping |
| 497 |
xend. |
| 498 |
*/ |
| 499 |
|
| 500 |
blsys->currentstep = 0; |
| 501 |
for(index = start_index; index < finish_index; index++, blsys->currentstep++) { |
| 502 |
xend = integrator_getsample(blsys, index+1); |
| 503 |
xprev = x[0]; |
| 504 |
asc_assert(xend > xprev); |
| 505 |
/* CONSOLE_DEBUG("LSODE call #%lu: x = [%f,%f]", index,xprev,xend); */ |
| 506 |
|
| 507 |
# ifdef ASC_SIGNAL_TRAPS |
| 508 |
|
| 509 |
Asc_SignalHandlerPushDefault(SIGFPE); |
| 510 |
Asc_SignalHandlerPushDefault(SIGINT); |
| 511 |
|
| 512 |
if(SETJMP(g_fpe_env)==0) { |
| 513 |
# endif /* ASC_SIGNAL_TRAPS */ |
| 514 |
|
| 515 |
/* CONSOLE_DEBUG("Calling LSODE with end-time = %f",xend); */ |
| 516 |
/* |
| 517 |
switch(mf){ |
| 518 |
case 10: |
| 519 |
CONSOLE_DEBUG("Non-stiff (Adams) method; no Jacobian will be used"); break; |
| 520 |
case 21: |
| 521 |
CONSOLE_DEBUG("Stiff (BDF) method, user-supplied full Jacobian"); break; |
| 522 |
case 22: |
| 523 |
CONSOLE_DEBUG("Stiff (BDF) method, internally generated full Jacobian"); break; |
| 524 |
case 24: |
| 525 |
CONSOLE_DEBUG("Stiff (BDF) method, user-supplied banded jacobian"); break; |
| 526 |
case 25: |
| 527 |
CONSOLE_DEBUG("Stiff (BDF) method, internally generated banded jacobian"); break; |
| 528 |
default: |
| 529 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid method id %d for LSODE",mf); |
| 530 |
return 0; * failure * |
| 531 |
} |
| 532 |
*/ |
| 533 |
|
| 534 |
d->lastwrite = clock(); |
| 535 |
|
| 536 |
// tolvect = 0 : scalars |
| 537 |
LSODE(&(LSODE_FEX), &my_neq, y, x, &xend, |
| 538 |
&itol, rtoler, atoler, &itask, &istate, |
| 539 |
&iopt ,rwork, &lrw, iwork, &liw, &(LSODE_JEX), &mf); |
| 540 |
|
| 541 |
res = dopri5 (my_neq, &integrator_dopri5_fex, x, y, xend, rtoler, atoler, itoler, integrator_dopri5_reporter, iout, |
| 542 |
stdout, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, NULL, 0); |
| 543 |
|
| 544 |
# ifdef ASC_SIGNAL_TRAPS |
| 545 |
} else { |
| 546 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to float error in LSODE call."); |
| 547 |
//dopri5_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx); |
| 548 |
return 6; |
| 549 |
} |
| 550 |
Asc_SignalHandlerPopDefault(SIGFPE); |
| 551 |
Asc_SignalHandlerPopDefault(SIGINT); |
| 552 |
|
| 553 |
# endif |
| 554 |
|
| 555 |
/* CONSOLE_DEBUG("AFTER %lu LSODE CALL\n", index); */ |
| 556 |
/* this check is better done in fex,jex, but lsode takes no status */ |
| 557 |
/* if (Solv_C_CheckHalt()) { |
| 558 |
if (istate >= 0) { |
| 559 |
istate=-7; |
| 560 |
} |
| 561 |
} |
| 562 |
*/ |
| 563 |
if(d->stop) { |
| 564 |
istate=-8; |
| 565 |
} |
| 566 |
|
| 567 |
if (istate < 0 ) { |
| 568 |
/* some kind of error occurred... */ |
| 569 |
ERROR_REPORTER_START_HERE(ASC_PROG_ERR); |
| 570 |
//lsode_write_istate(istate); |
| 571 |
FPRINTF(ASCERR, "\nFurthest point reached was t = %g.\n",x[0]); |
| 572 |
error_reporter_end_flush(); |
| 573 |
|
| 574 |
//dopri5_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx); |
| 575 |
integrator_output_close(blsys); |
| 576 |
return 7; |
| 577 |
} |
| 578 |
|
| 579 |
if (d->status==lsode_nok) { |
| 580 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to an error in derivative computations."); |
| 581 |
//lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx); |
| 582 |
//d->status = lsode_ok; /* clean up before we go */ |
| 583 |
//d->lastcall = lsode_none; |
| 584 |
integrator_output_close(blsys); |
| 585 |
return 8; |
| 586 |
} |
| 587 |
|
| 588 |
integrator_setsample(blsys, index+1, x[0]); |
| 589 |
/* record when lsode actually came back */ |
| 590 |
integrator_set_t(blsys, x[0]); |
| 591 |
integrator_set_y(blsys, y); |
| 592 |
/* put x,y in d in case lsode got x,y by interpolation, as it does */ |
| 593 |
|
| 594 |
reporterstatus = integrator_output_write(blsys); |
| 595 |
|
| 596 |
if(reporterstatus==0) { |
| 597 |
ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration cancelled"); |
| 598 |
//lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx); |
| 599 |
//d->status = lsode_ok; |
| 600 |
//d->lastcall = lsode_none; |
| 601 |
integrator_output_close(blsys); |
| 602 |
return 9; |
| 603 |
} |
| 604 |
|
| 605 |
if (nobs > 0) { |
| 606 |
# ifdef ASC_SIGNAL_TRAPS |
| 607 |
if (SETJMP(g_fpe_env)==0) { |
| 608 |
# endif /* ASC_SIGNAL_TRAPS */ |
| 609 |
|
| 610 |
/* solve for obs since d isn't necessarily already |
| 611 |
computed there though lsode's x and y may be. |
| 612 |
Note that since lsode usually steps beyond xend |
| 613 |
x1 usually wouldn't be x0 precisely if the x1/x0 |
| 614 |
scheme worked, which it doesn't anyway. */ |
| 615 |
|
| 616 |
//LSODEDATA_SET(blsys); |
| 617 |
//LSODE_FEX(&my_neq, x, y, dydx); |
| 618 |
//LSODEDATA_RELEASE(); |
| 619 |
|
| 620 |
/* calculate observations, if any, at returned x and y. */ |
| 621 |
obs = integrator_get_observations(blsys, obs); |
| 622 |
|
| 623 |
integrator_output_write_obs(blsys); |
| 624 |
|
| 625 |
# ifdef ASC_SIGNAL_TRAPS |
| 626 |
}else{ |
| 627 |
ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to float error in LSODE FEX call."); |
| 628 |
lsode_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx); |
| 629 |
d->status = lsode_ok; /* clean up before we go */ |
| 630 |
d->lastcall = lsode_none; |
| 631 |
integrator_output_close(blsys); |
| 632 |
return 10; |
| 633 |
} |
| 634 |
# endif /* ASC_SIGNAL_TRAPS */ |
| 635 |
} |
| 636 |
/* CONSOLE_DEBUG("Integration completed from %3g to %3g.",xprev,x[0]); */ |
| 637 |
} |
| 638 |
|
| 639 |
integrator_output_close(blsys); |
| 640 |
|
| 641 |
CONSOLE_DEBUG("--- DOPRI5 done ---"); |
| 642 |
return 0; /* success */ |
| 643 |
} |
| 644 |
|
| 645 |
#endif |
| 646 |
|
Parent Directory
| 
Revision Log