/[ascend]/trunk/solvers/dopri5/asc_dopri5.c
ViewVC logotype

Contents of /trunk/solvers/dopri5/asc_dopri5.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1757 - (show annotations) (download) (as text)
Wed Feb 27 05:40:40 2008 UTC (16 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 20970 byte(s)
Add DOPRI5 to installer package (EXPERIMENTAL).
On Windows, use the registry to locate the book.pdf documentation file.
Fix path to test models for DOPRI5 in test.py.
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/ascPanic.h>
29 #include <utilities/ascSignal.h>
30 #include <utilities/error.h>
31 #include <general/ospath.h>
32 #include <integrator/integrator.h>
33 #include <system/slv_stdcalls.h>
34 #include <solver/solver.h>
35 #include "dopri5.h"
36
37 #define INTEG_DOPRI5 5
38
39 /* #define STATS_DEBUG */
40 /* #define DOPRI5_DEBUG */
41
42 IntegratorCreateFn integrator_dopri5_create;
43 IntegratorParamsDefaultFn integrator_dopri5_params_default;
44 IntegratorSolveFn integrator_dopri5_solve;
45 IntegratorFreeFn integrator_dopri5_free;
46 IntegratorWriteMatrixFn integrator_dopri5_write_matrix;
47
48 const IntegratorInternals integrator_dopri5_internals =
49 {
50 integrator_dopri5_create
51 ,integrator_dopri5_params_default
52 ,integrator_analyse_ode /* note, this routine is back in integrator.c */
53 ,integrator_dopri5_solve
54 ,NULL
55 ,NULL /* debugfn */
56 ,integrator_dopri5_free
57 ,INTEG_DOPRI5
58 ,"DOPRI5"
59 };
60
61 extern ASC_EXPORT int dopri5_register(void)
62 {
63 CONSOLE_DEBUG("Registering DOPRI5...");
64 return integrator_register(&integrator_dopri5_internals);
65 }
66
67 enum dopri5_status{
68 DOPRI5_SUCCESS=1
69 ,DOPRI5_INTERRUPT=2
70 ,DOPRI5_BADINPUT=-1
71 ,DOPRI5_ITERLIMIT=-2
72 ,DOPRI5_STEPSMALL=-3
73 ,DOPRI5_STIFF=-4
74 };
75
76 typedef struct IntegratorDopri5DataStruct
77 {
78 long n_eqns; /**< dimension of state vector */
79 int *input_indices; /**< vector of state vars indexes */
80 int *output_indices; /**< vector of derivative var indexes */
81 struct var_variable **y_vars; /**< NULL-terminated list of states vars */
82 struct var_variable **ydot_vars; /**< NULL-terminated list of derivative vars*/
83 struct rel_relation **rlist; /**< NULL-terminated list of relevant rels
84 to be differentiated */
85 long currentsample;
86 char stop; /* stop requested? */
87 int partitioned; /* partioned func evals or not */
88 double *yinter; /* interpolated y values */
89
90 clock_t lastwrite; /* time of last call to the reporter 'write' function */
91 }
92 IntegratorDopri5Data;
93
94 /*------------------------------------------------------------------------------
95 CREATE/FREE
96 */
97
98 void integrator_dopri5_create(struct IntegratorSystemStruct *blsys){
99 IntegratorDopri5Data *d;
100 d = ASC_NEW_CLEAR(IntegratorDopri5Data);
101 d->n_eqns=0;
102 d->input_indices=NULL;
103 d->output_indices=NULL;
104 d->y_vars=NULL;
105 d->ydot_vars=NULL;
106 d->rlist=NULL;
107 blsys->enginedata=(void*)d;
108 integrator_dopri5_params_default(blsys);
109 CONSOLE_DEBUG("CREATED DOPRI5");
110 }
111
112 void integrator_dopri5_free(void *enginedata){
113 IntegratorDopri5Data d;
114 d = *((IntegratorDopri5Data *)enginedata);
115
116 if(d.input_indices)ASC_FREE(d.input_indices);
117 d.input_indices = NULL;
118
119 if(d.output_indices)ASC_FREE(d.output_indices);
120 d.output_indices = NULL;
121
122 if(d.y_vars)ASC_FREE(d.y_vars);
123 d.y_vars = NULL;
124
125 if(d.ydot_vars)ASC_FREE(d.ydot_vars);
126 d.ydot_vars = NULL;
127
128 if(d.rlist)ASC_FREE(d.rlist);
129 d.rlist = NULL;
130
131 if(d.yinter)ASC_FREE(d.yinter);
132 d.yinter = NULL;
133
134 d.n_eqns = 0L;
135 }
136
137 /*------------------------------------------------------------------------------
138 PARAMETERS
139 */
140
141 enum dopri5_parameters{
142 DOPRI5_PARAM_RTOL
143 ,DOPRI5_PARAM_ATOL
144 ,DOPRI5_PARAM_TOLVECT
145 ,DOPRI5_PARAM_NSTIFF
146 ,DOPRI5_PARAMS_SIZE
147 #if 0
148 // more parameters for adding later...
149 SolTrait *solout, /* function providing the numerical solution during integration */
150 int iout, /* switch for calling solout */
151
152 double uround, /* rounding unit */
153 double safe, /* safety factor */
154 double fac1, /* parameters for step size selection */
155 double fac2,
156 double beta, /* for stabilized step size control */
157 #endif
158 };
159
160
161 int integrator_dopri5_params_default(IntegratorSystem *blsys){
162
163 asc_assert(blsys!=NULL);
164 asc_assert(blsys->engine==INTEG_DOPRI5);
165 slv_parameters_t *p;
166 p = &(blsys->params);
167
168 slv_destroy_parms(p);
169
170 if(p->parms==NULL) {
171 p->parms = ASC_NEW_ARRAY(struct slv_parameter, DOPRI5_PARAMS_SIZE);
172 if(p->parms==NULL)return -1;
173 p->dynamic_parms = 1;
174 } else {
175 asc_assert(p->num_parms == DOPRI5_PARAMS_SIZE);
176 }
177
178 /* reset the number of parameters to zero so that we can check it at the end */
179 p->num_parms = 0;
180
181 slv_param_bool(p,DOPRI5_PARAM_TOLVECT
182 ,(SlvParameterInitBool) {
183 {"tolvect"
184 ,"Use per-variable tolerances?",1
185 ,"If TRUE, values of 'ode_rtol' and 'ode_atol' are taken from your"
186 " model and used in the integration. If FALSE, scalar values are"
187 " used (see 'rtol' and 'atol') and shared by all variables. See"
188 " 'itoler' in DOPRI5 code."
189 }
190 , TRUE
191 }
192 );
193
194 slv_param_real(p,DOPRI5_PARAM_ATOL
195 ,(SlvParameterInitReal) {
196 {"atol"
197 ,"Scalar absolute error tolerance",1
198 ,"Default value of the scalar absolute error tolerance (for cases"
199 " where not specified in oda_atol var property. See 'dopri5.h' for"
200 " details"
201 }
202 , 1e-7, 1e-15, 1e10
203 }
204 );
205
206 slv_param_real(p,DOPRI5_PARAM_RTOL
207 ,(SlvParameterInitReal) {
208 {"rtol"
209 ,"Scalar relative error tolerance",1
210 ,"Default value of the scalar relative error tolerance (for cases"
211 " where not specified in oda_rtol var property. See 'dopri5.h' for"
212 " details"
213 }
214 , 1e-7, 1e-15, 1
215 }
216 );
217
218 slv_param_int(p,DOPRI5_PARAM_NSTIFF
219 ,(SlvParameterInitInt) {
220 {"nstiff"
221 ,"Number of steps considered stiff", 1
222 ,"Test for stiffness is activated when the current step number is a"
223 " multiple of nstiff. A negative value means no test."
224 }
225 , 1000, -1, 1e6
226 }
227 );
228
229 /*
230 slv_param_bool(p,DOPRI5_PARAMS_DENSEREPORTING
231 ,(SlvParameterInitBool) {
232 {"densereporting"
233 ,"Dense reporting?",1
234 ,"If TRUE, output will be made at every sub-timestep"
235 " during integration. If false, output is only made
236 " according to the samplelist."
237 }
238 , FALSE
239 }
240 );
241
242 slv_param_char(p,DOPRI5_PARAM_METH
243 ,(SlvParameterInitChar) {
244 {"meth"
245 ,"Integration method",1
246 ,"AM=Adams-Moulton method (for non-stiff problems), BDF=Backwards"
247 " Difference Formular (for stiff problems). See 'Description and"
248 " Use of DOPRI5', section 3.1."
249 }
250 , "BDF"
251 }
252 , (char *[]) {"AM","BDF",NULL
253 }
254 );
255 */
256
257 asc_assert(p->num_parms == DOPRI5_PARAMS_SIZE);
258 return 0;
259 }
260
261 /* solver_var children expected for state variables */
262 static symchar *g_symbols[2];
263 #define STATERTOL g_symbols[0]
264 #define STATEATOL g_symbols[1]
265
266 static void InitTolNames(void){
267 STATERTOL = AddSymbol("ode_rtol");
268 STATEATOL = AddSymbol("ode_atol");
269 }
270
271
272 /**
273 Allocates, fills, and returns the atoler or rtoler data, which may be either
274 a vector of values for each variable, or may be a single value to be
275 shared by all.
276
277 State variables missing child ode_rtol will be defaulted to RTOL
278
279 @param is_r TRUE if we want RTOL, FALSE if we want ATOL.
280 */
281 static double *dopri5_get_artol(IntegratorSystem *blsys, int is_r, int tolvect) {
282
283 struct Instance *toli;
284 double tolval, *tol;
285 int i,len;
286 symchar *tolname;
287
288 len = blsys->n_y;
289
290 if(!tolvect){
291
292 // single tol for all vars
293 tolval = SLV_PARAM_REAL(&(blsys->params),DOPRI5_PARAM_RTOL);
294 CONSOLE_DEBUG("Using RTOL = %g for all vars", tolval);
295
296 tol = ASC_NEW(double);
297 if(tol == NULL){
298 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
299 return tol;
300 }
301
302 *tol = tolval;
303 return tol;
304
305 }else{
306 tol = ASC_NEW_ARRAY(double, blsys->n_y+1);
307 if (tol == NULL) {
308 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
309 return tol;
310 }
311
312 tolval = SLV_PARAM_REAL(&(blsys->params),DOPRI5_PARAM_RTOL);
313
314 InitTolNames(); // from where?
315 if(is_r)tolname = STATERTOL;
316 else tolname = STATEATOL;
317
318 for(i=0; i<len; i++){
319 toli = ChildByChar(var_instance(blsys->y[i]),tolname);
320 if(toli == NULL || !AtomAssigned(toli)){
321 tol[i] = tolval;
322 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Assuming value %3g"
323 "for '%s' child undefined for state variable %ld."
324 ,tol[i], SCP(tolname), blsys->y_id[i]
325 );
326 }else{
327 tol[i] = RealAtomValue(toli);
328 }
329 }
330 }
331 tol[len] = SLV_PARAM_REAL(&(blsys->params),DOPRI5_PARAM_RTOL);
332 return tol;
333 }
334
335 /*------------------------------------------------------------------------------
336 STATS
337 */
338
339 /*
340 Several functions provide access to different values :
341
342 xRead x value for which the solution has been computed (x=xend after
343 successful return).
344
345 hRead Predicted step size of the last accepted step (useful for a
346 subsequent call to dopri5).
347
348 nstepRead Number of used steps.
349 naccptRead Number of accepted steps.
350 nrejctRead Number of rejected steps.
351 nfcnRead Number of function calls.
352 */
353 typedef struct IntegratorDopri5StatsStruct{
354 long nfcn;
355 long nstep;
356 long naccpt;
357 long nrejct;
358 double h;
359 double x;
360 } IntegratorDopri5Stats;
361
362
363 /*----------------------------------------------
364 STATS
365 */
366
367 /**
368 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
369 status in a struct instead of separately.
370
371 @return 0 on success.
372 */
373 int integrator_dopri5_stats(IntegratorSystem *blsys, IntegratorDopri5Stats *s){
374
375 #define S(NAME) s->NAME = NAME##Read()
376 S(nfcn); S(nstep); S(naccpt); S(nrejct);
377 S(h); S(x);
378 #undef S
379
380 return 0;
381 }
382
383 /**
384 This routine just outputs the stats to the CONSOLE_DEBUG routine.
385
386 @TODO provide a GUI way of stats reporting from DOPRI5
387 */
388 void integrator_dopri5_write_stats(IntegratorDopri5Stats *stats){
389 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
390 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
391 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
392 SL(nfcn); SL(nstep); SL(naccpt); SL(nrejct);
393 SR(h); SR(x);
394 # undef SL
395 # undef SI
396 # undef SR
397 }
398
399 /*------------------------------------------------------------------------------
400 FUNCTION EVALUATION
401 */
402
403 static FcnEqDiff integrator_dopri5_fex;
404
405 /**
406 Evaluation function.
407 @param n_eq number of equations, number of y and number of ydot.
408 @param t indep var value
409 @param y input vector of variable values
410 @param ydot return vector of calculated derivatives
411 @param user_data point to whatever we want, in this case the IntegratorSystem
412 */
413 static void integrator_dopri5_fex(
414 unsigned n_eq, double t, double *y, double *ydot
415 , void *user_data
416 ){
417 slv_status_t status;
418 IntegratorSystem *blsys = (IntegratorSystem *)user_data;
419
420 int i;
421 unsigned long res;
422
423 //CONSOLE_DEBUG("Calling for a function evaluation");
424
425 /* pass the time and the unknowns back to the System */
426 integrator_set_t(blsys, t);
427 integrator_set_y(blsys, y);
428 //CONSOLE_DEBUG("t = %f: y[0] = %f",t,y[0]);
429
430 asc_assert(blsys->system);
431
432 slv_resolve(blsys->system);
433
434 if((res = slv_solve(blsys->system))){
435 CONSOLE_DEBUG("solver returns error %ld",res);
436 }
437
438 slv_get_status(blsys->system, &status);
439
440
441 if(slv_check_bounds(blsys->system,0,-1,"")){
442 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Variables went outside boundaries...");
443 // TODO relay that system has gone out of bounds
444 }
445
446 /* pass the NLA solver status to the integrator */
447 res = integrator_checkstatus(status);
448
449 integrator_output_write(blsys);
450
451 if(res){
452 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for derivatives (%d)",res);
453 #if 1
454 ERROR_REPORTER_START_HERE(ASC_PROG_ERR);
455 FPRINTF(ASCERR,"Unable to compute the vector of derivatives with the following values for the state variables:\n");
456 for (i = 0; i< n_eq; i++) {
457 FPRINTF(ASCERR,"y[%4d] = %g\n",i, y[i]);
458 }
459 error_reporter_end_flush();
460 #endif
461 #ifdef ASC_SIGNAL_TRAPS
462 raise(SIGINT);
463 #endif
464 }else{
465 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"lsodedata->status = %d",lsodedata->status); */
466 }
467
468 integrator_get_ydot(blsys, ydot);
469
470 #ifdef DOPRI5_DEBUG
471 CONSOLE_DEBUG("y[0]=%e,y[1]=%e --> ydot[0]=%e,ydot[1]=%e",y[0],y[1],ydot[0],ydot[1]);
472 #endif
473 //CONSOLE_DEBUG("ydot[0] = %f",ydot[0]);
474 // DONE, OK
475 }
476
477 /*------------------------------------------------------------------------------
478 REPORTING
479 */
480
481 static SolTrait integrator_dopri5_reporter;
482
483 static void integrator_dopri5_reporter(
484 long nr, double told, double t, double* y
485 , unsigned n, int* irtrn, void *user_data
486 ){
487 double ts;
488 IntegratorSystem *blsys = (IntegratorSystem *)user_data;
489 IntegratorDopri5Data *d = (IntegratorDopri5Data *)(blsys->enginedata);
490
491 ts = integrator_getsample(blsys,d->currentsample);
492
493 #if 0 /* sub-step output */
494 int i;
495 enum dopri5_status res;
496 slv_status_t status;
497
498 while(t>ts){
499 for(i=0; i<nr; ++i){
500 d->yinter[i] = contd5(i,ts);
501 }
502 integrator_set_y(blsys, d->yinter);
503
504 if((res = slv_solve(blsys->system))){
505 CONSOLE_DEBUG("solver returns error %ld",res);
506 }
507
508 slv_get_status(blsys->system, &status);
509
510 if(slv_check_bounds(blsys->system,0,-1,"")){
511 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Variables went outside boundaries...");
512 // TODO relay that system has gone out of bounds
513 }
514
515 /* pass the NLA solver status to the integrator */
516 res = integrator_checkstatus(status);
517
518 if(res){
519 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve inter-step state (%d)",res);
520
521 #ifdef ASC_SIGNAL_TRAPS
522 raise(SIGINT);
523 #endif
524 }
525
526 integrator_output_write_obs(blsys);
527
528 d->currentsample++;
529 ts = integrator_getsample(blsys,d->currentsample);
530 }
531
532 integrator_output_write(blsys);
533 #else
534 ts = integrator_getsample(blsys,d->currentsample);
535 if(t>ts){
536 //CONSOLE_DEBUG("t=%f > ts=%f (currentsample = %ld",t,ts,d->currentsample);
537 integrator_output_write_obs(blsys);
538 #ifdef DOPRI5_DEBUG
539 CONSOLE_DEBUG("step = %ld", nr-1);
540 #endif
541 while(t>ts){
542 d->currentsample++;
543 blsys->currentstep++;
544 ts = integrator_getsample(blsys,d->currentsample);
545 }
546 }
547
548 //CONSOLE_DEBUG("t = %f, y[0] = %f",t,y[0]);
549 integrator_output_write(blsys);
550 #endif
551 }
552
553 /*------------------------------------------------------------------------------
554 SOLVE
555 */
556
557 #define DOPRI5_FREE CONSOLE_DEBUG("FREE DOPRI5")
558
559 int integrator_dopri5_solve(IntegratorSystem *blsys
560 , unsigned long start_index, unsigned long finish_index
561 ){
562 IntegratorDopri5Data *d;
563 slv_status_t status;
564
565 double x;
566 double xend;
567 unsigned long nsamples, neq;
568 long nobs;
569 //int itol, itask, mf, lrw, liw;
570 //unsigned long index;
571 //int istate, iopt;
572 //double * rwork;
573 //int * iwork;
574 double *y, *atoler, *rtoler, *obs;
575 int my_neq;
576 enum dopri5_status res;
577
578 double hmax, h;
579 long nmax;
580 long nstiff;
581
582 #if 0
583 const char *method; /* Table 3.1 in D&Uo... */
584 int miter; /* Table 3.2 in D&Uo... */
585 int maxord; /* page 92 in D&Uo... */
586 #endif
587
588 d = (IntegratorDopri5Data *)(blsys->enginedata);
589
590 /* the numer of equations must be equal to blsys->n_y, the number of states */
591 d->n_eqns = blsys->n_y;
592 assert(d->n_eqns>0);
593
594 d->input_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
595 d->output_indices = ASC_NEW_ARRAY_CLEAR(int, d->n_eqns);
596
597 d->y_vars = ASC_NEW_ARRAY(struct var_variable *,d->n_eqns+1);
598 d->ydot_vars = ASC_NEW_ARRAY(struct var_variable *, d->n_eqns+1);
599
600 d->yinter = ASC_NEW_ARRAY(double,d->n_eqns);
601
602 /* set up the NLA solver here */
603
604 /*
605 DOPRI5 should be OK to deal with any linsol/linsolqr-based solver.
606 But for the moment we restrict to just QRSlv :-(
607 */
608 if(strcmp(slv_solver_name(slv_get_selected_solver(blsys->system)),"QRSlv") != 0) {
609 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"QRSlv must be selected before integration.");
610 return 1;
611 }
612
613 CONSOLE_DEBUG("Solver selected is '%s'",slv_solver_name(slv_get_selected_solver(blsys->system)));
614
615 slv_get_status(blsys->system, &status);
616
617 if(status.struct_singular){
618 ERROR_REPORTER_HERE(ASC_USER_WARNING
619 ,"The system (according to QRSlv) is structurally singular."
620 " The ODE system may also be singular, but not necessarily."
621 );
622 /* d->status = lsode_nok;
623 return 2;*/
624 }
625
626 /* here we assume integrators.c is in charge of dynamic loading */
627
628 /* set up parameters for sending to DOPRI5 */
629
630 nsamples = integrator_getnsamples(blsys);
631 if (nsamples <2) {
632 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration will not be performed. The system has no end sample time defined.");
633 //d->status = dopri5_nok;
634 return 3;
635 }
636 neq = blsys->n_y;
637 nobs = blsys->n_obs;
638
639 #if 0
640 // TODO implement these:
641 unsigned nrdens, /* number of components for which dense outpout is required */
642 unsigned* icont, /* indexes of components for which dense output is required, >= nrdens */
643 unsigned licont /* declared length of icon */
644 #endif
645
646 int iout = 2; /* SLV_PARAM_BOOL(&(blsys->params),DOPRI5_PARAM_DENSEREPORTING) */
647
648 int tolvect = SLV_PARAM_BOOL(&(blsys->params),DOPRI5_PARAM_TOLVECT);
649
650 /* samplelist_debug(blsys->samples); */
651
652 x = integrator_getsample(blsys, 0);
653 d->currentsample = 1;
654
655 y = integrator_get_y(blsys, NULL);
656
657 rtoler = dopri5_get_artol(blsys,1,tolvect);
658 atoler = dopri5_get_artol(blsys,0,tolvect);
659
660 obs = integrator_get_observations(blsys, NULL);
661
662 // TODO check memory allocations
663
664 h = integrator_get_stepzero(blsys);
665 hmax = integrator_get_maxstep(blsys);
666 CONSOLE_DEBUG("init step = %f, max step = %f", h, hmax);
667
668 /* rwork[6] = integrator_get_minstep(blsys); */ /* ignored */
669 nmax = integrator_get_maxsubsteps(blsys);
670
671 nstiff = SLV_PARAM_INT(&(blsys->params),DOPRI5_PARAM_NSTIFF);
672
673 if(x > integrator_getsample(blsys, 2)) {
674 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid initialisation time: exceeds second timestep value");
675 return 5;
676 }
677
678 /* put the values from derivative system into the record */
679 integrator_setsample(blsys, start_index, x);
680
681 integrator_output_init(blsys);
682
683 /* -- store the initial values of all the stuff */
684 integrator_output_write(blsys);
685 integrator_output_write_obs(blsys);
686
687 my_neq = (int)neq;
688
689 unsigned *icont = NULL;
690 unsigned nrdens = 0;
691 unsigned licont = nrdens;
692 #if 0
693 unsigned licont = 2;
694 icont = ASC_NEW_ARRAY(unsigned, licont);
695 icont[0] = 0;
696 icont[1] = 1;
697 #endif
698
699 blsys->currentstep = 0;
700
701 xend = integrator_getsample(blsys, finish_index);
702
703 # ifdef ASC_SIGNAL_TRAPS
704
705 Asc_SignalHandlerPushDefault(SIGFPE);
706 Asc_SignalHandlerPushDefault(SIGINT);
707
708 if(SETJMP(g_fpe_env)==0) {
709 # endif /* ASC_SIGNAL_TRAPS */
710
711 d->lastwrite = clock();
712
713 res = dopri5 (my_neq, &integrator_dopri5_fex, x, y, xend
714 , rtoler, atoler, tolvect, integrator_dopri5_reporter, iout
715 , stdout, 0.0 /* uround */, 0.0 /*safe*/, 0.0 /*fac1*/, 0.0/*fac2*/
716 , 0.0 /* beta */, hmax, h, nmax, 0
717 , nstiff, nrdens, icont, licont
718 , (void *)blsys
719 );
720
721 # ifdef ASC_SIGNAL_TRAPS
722 }else{
723 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Integration terminated due to float error in DOPRI5 call.");
724 //dopri5_free_mem(y,reltol,abtol,rwork,iwork,obs,dydx);
725 DOPRI5_FREE;
726 return 6;
727 }
728 Asc_SignalHandlerPopDefault(SIGFPE);
729 Asc_SignalHandlerPopDefault(SIGINT);
730
731 # endif
732
733 switch(res){
734 case DOPRI5_SUCCESS:
735 break;
736 case DOPRI5_INTERRUPT:
737 ERROR_REPORTER_HERE(ASC_USER_ERROR,"DOPRI5 interrupted by user");
738 break;
739 case DOPRI5_BADINPUT:
740 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Bad input to DOPRI5");
741 break;
742 case DOPRI5_ITERLIMIT:
743 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Iteration limit exceeded in DOPRI5");
744 break;
745 case DOPRI5_STEPSMALL:
746 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Step size became too small in DOPRI5");
747 break;
748 case DOPRI5_STIFF:
749 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem appears stiff in DOPRI5");
750 break;
751 }
752
753 if(res<0){
754 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Furthest point reached was t = %g.\n",x);
755 DOPRI5_FREE;
756 return 7;
757 }
758
759 /* write final step output */
760 integrator_output_write_obs(blsys);
761
762 #if 0
763 integrator_setsample(blsys, index+1, x);
764 /* record when dopri5 actually came back */
765 integrator_set_t(blsys, x);
766 integrator_set_y(blsys, y);
767 /* put x,y in d in case dopri5 got x,y by interpolation, as it does */
768
769 }
770 #endif
771
772 integrator_output_close(blsys);
773
774 #ifdef STATS_DEBUG
775 IntegratorDopri5Stats stats;
776 if(0 == integrator_dopri5_stats(blsys, &stats)){
777 integrator_dopri5_write_stats(&stats);
778 }else{
779 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
780 }
781 #endif
782
783 CONSOLE_DEBUG("--- DOPRI5 done ---");
784 return 0; /* success */
785 }
786

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22