/[ascend]/trunk/base/generic/solver/ida.c
ViewVC logotype

Contents of /trunk/base/generic/solver/ida.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 980 - (show annotations) (download) (as text)
Wed Dec 20 23:23:50 2006 UTC (17 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 31761 byte(s)
Fixed detection of libpython on Windows.
Fixed TestLSODE.testlotka test.py
Removed a sundials_dense.h assertion under windows.
1 /* ASCEND modelling environment
2 Copyright (C) 2006 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 *//**
19 @file
20 Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
21 as part of the GPL-licensed SUNDIALS solver package from LLNL.
22 @see http://www.llnl.gov/casc/sundials/
23 *//*
24 by John Pye, May 2006
25 */
26
27 /*
28 Be careful with the following. This file requires both the 'ida.h' from
29 SUNDIALS as well as the 'ida.h' from ASCEND. Make sure that we're getting
30 both of these; if you get problems check your build tool for the paths being
31 passed to the C preprocessor.
32 */
33
34 /* standard includes */
35 #include <signal.h>
36
37 /* ASCEND includes */
38 #include "ida.h"
39 #include <utilities/error.h>
40 #include <utilities/ascConfig.h>
41 #include <utilities/ascSignal.h>
42 #include <utilities/ascPanic.h>
43 #include <compiler/instance_enum.h>
44 #include "var.h"
45 #include "rel.h"
46 #include "discrete.h"
47 #include "conditional.h"
48 #include "logrel.h"
49 #include "bnd.h"
50 #include "linsol.h"
51 #include "linsolqr.h"
52 #include "slv_common.h"
53 #include "slv_client.h"
54 #include "relman.h"
55
56 /* SUNDIALS includes */
57 #ifdef ASC_WITH_IDA
58 # include <sundials/sundials_config.h>
59 # include <sundials/sundials_dense.h>
60 # include <ida/ida.h>
61 # include <nvector/nvector_serial.h>
62 # include <ida/ida_spgmr.h>
63 # include <ida/ida_spbcgs.h>
64 # include <ida/ida_sptfqmr.h>
65 # include <ida/ida_dense.h>
66 # ifndef IDA_SUCCESS
67 # error "Failed to include SUNDIALS IDA header file"
68 # endif
69 #endif
70
71 /*
72 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
73 */
74 #ifndef SUNDIALS_VERSION_MINOR
75 # ifdef __GNUC__
76 # warning "GUESSING SUNDIALS VERSION 2.2"
77 # endif
78 # define SUNDIALS_VERSION_MINOR 2
79 #endif
80 #ifndef SUNDIALS_VERSION_MAJOR
81 # define SUNDIALS_VERSION_MAJOR 2
82 #endif
83
84 /* check that we've got what we expect now */
85 #ifndef ASC_IDA_H
86 # error "Failed to include ASCEND IDA header file"
87 #endif
88
89 /* #define FEX_DEBUG */
90 /* #define JEX_DEBUG */
91 /* #define SOLVE_DEBUG */
92
93 const IntegratorInternals integrator_ida_internals = {
94 integrator_ida_create
95 ,integrator_ida_params_default
96 ,integrator_analyse_dae /* note, this routine is back in integrator.c */
97 ,integrator_ida_solve
98 ,integrator_ida_free
99 ,INTEG_IDA
100 ,"IDA"
101 };
102
103 /**
104 Struct containing any stuff that IDA needs that doesn't fit into the
105 common IntegratorSystem struct.
106 */
107 typedef struct{
108 struct rel_relation **rellist; /**< NULL terminated list of rels */
109 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
110 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
111 int nrels;
112 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
113 } IntegratorIdaData;
114
115 /*-------------------------------------------------------------
116 FORWARD DECLS
117 */
118 /* residual function forward declaration */
119 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
120
121 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
122 , N_Vector v, N_Vector Jv, realtype c_j
123 , void *jac_data, N_Vector tmp1, N_Vector tmp2
124 );
125
126 /* error handler forward declaration */
127 void integrator_ida_error(int error_code
128 , const char *module, const char *function
129 , char *msg, void *eh_data
130 );
131
132 int integrator_ida_djex(long int Neq, realtype tt
133 , N_Vector yy, N_Vector yp, N_Vector rr
134 , realtype c_j, void *jac_data, DenseMat Jac
135 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
136 );
137
138 /*-------------------------------------------------------------
139 SETUP/TEARDOWN ROUTINES
140 */
141 void integrator_ida_create(IntegratorSystem *blsys){
142 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
143 IntegratorIdaData *enginedata;
144 enginedata = ASC_NEW(IntegratorIdaData);
145 enginedata->rellist = NULL;
146 enginedata->varlist = NULL;
147 enginedata->safeeval = 0;
148 blsys->enginedata = (void *)enginedata;
149 integrator_ida_params_default(blsys);
150 }
151
152 void integrator_ida_free(void *enginedata){
153 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
154 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
155 /* note, we don't own the rellist, so don't need to free it */
156 ASC_FREE(d);
157 }
158
159 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
160 IntegratorIdaData *d;
161 assert(blsys!=NULL);
162 assert(blsys->enginedata!=NULL);
163 assert(blsys->engine==INTEG_IDA);
164 d = ((IntegratorIdaData *)(blsys->enginedata));
165 return d;
166 }
167
168 /*-------------------------------------------------------------
169 PARAMETERS FOR IDA
170 */
171
172 enum ida_parameters{
173 IDA_PARAM_LINSOLVER
174 ,IDA_PARAM_MAXL
175 ,IDA_PARAM_AUTODIFF
176 ,IDA_PARAM_CALCIC
177 ,IDA_PARAM_SAFEEVAL
178 ,IDA_PARAM_RTOL
179 ,IDA_PARAM_ATOL
180 ,IDA_PARAM_ATOLVECT
181 ,IDA_PARAM_GSMODIFIED
182 ,IDA_PARAMS_SIZE
183 };
184
185 /**
186 Here the full set of parameters is defined, along with upper/lower bounds,
187 etc. The values are stuck into the blsys->params structure.
188
189 @return 0 on success
190 */
191 int integrator_ida_params_default(IntegratorSystem *blsys){
192 asc_assert(blsys!=NULL);
193 asc_assert(blsys->engine==INTEG_IDA);
194 slv_parameters_t *p;
195 p = &(blsys->params);
196
197 slv_destroy_parms(p);
198
199 if(p->parms==NULL){
200 CONSOLE_DEBUG("params NULL");
201 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
202 if(p->parms==NULL)return -1;
203 p->dynamic_parms = 1;
204 }else{
205 CONSOLE_DEBUG("params not NULL");
206 }
207
208 /* reset the number of parameters to zero so that we can check it at the end */
209 p->num_parms = 0;
210
211 slv_param_bool(p,IDA_PARAM_AUTODIFF
212 ,(SlvParameterInitBool){{"autodiff"
213 ,"Use auto-diff?",1
214 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
215 }, TRUE}
216 );
217
218 slv_param_bool(p,IDA_PARAM_CALCIC
219 ,(SlvParameterInitBool){{"calcic"
220 ,"Calculate initial conditions?",1
221 ,"Use IDA to calculate initial conditions (1) or else assume that the model will already be solved for this case (0)"
222 }, FALSE}
223 );
224
225 slv_param_bool(p,IDA_PARAM_SAFEEVAL
226 ,(SlvParameterInitBool){{"safeeval"
227 ,"Use safe evaluation?",1
228 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
229 "throw SIGFPE errors which will then halt integration."
230 }, FALSE}
231 );
232
233
234 slv_param_bool(p,IDA_PARAM_ATOLVECT
235 ,(SlvParameterInitBool){{"atolvect"
236 ,"Use 'ode_atol' values as specified?",1
237 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
238 " in the integration. If FALSE, a scalar absolute tolerance value"
239 " is shared by all variables. See IDA manual, section 5.5.1"
240 }, TRUE }
241 );
242
243 slv_param_real(p,IDA_PARAM_ATOL
244 ,(SlvParameterInitReal){{"atol"
245 ,"Scalar absolute error tolerance",1
246 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
247 " See IDA manual, section 5.5.1"
248 }, 1e-5, DBL_MIN, DBL_MAX }
249 );
250
251 slv_param_real(p,IDA_PARAM_RTOL
252 ,(SlvParameterInitReal){{"rtol"
253 ,"Scalar relative error tolerance",1
254 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
255 " it's not possible to set per-variable relative tolerances as it is"
256 " with LSODE)."
257 " See IDA manual, section 5.5.1"
258 }, 1e-4, 0, DBL_MAX }
259 );
260
261 slv_param_char(p,IDA_PARAM_LINSOLVER
262 ,(SlvParameterInitChar){{"linsolver"
263 ,"Linear solver",1
264 ,"See IDA manual, section 5.5.3."
265 }, "SPGMR"}, (char *[]){"DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
266 );
267
268 slv_param_int(p,IDA_PARAM_MAXL
269 ,(SlvParameterInitInt){{"maxl"
270 ,"Maximum Krylov dimension",0
271 ,"The maximum dimension of Krylov space used by the linear solver"
272 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
273 " The default of 0 results in IDA using its internal default, which"
274 " is currently a value of 5."
275 }, 0, 0, 20 }
276 );
277
278 slv_param_bool(p,IDA_PARAM_GSMODIFIED
279 ,(SlvParameterInitBool){{"gsmodified"
280 ,"Use modified Gram-Schmidt orthogonalisation in SPGMR?",2
281 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
282 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
283 }, TRUE}
284 );
285
286 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
287
288 CONSOLE_DEBUG("Created %d params", p->num_parms);
289
290 return 0;
291 }
292
293 /*-------------------------------------------------------------
294 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
295 */
296
297 /* return 1 on success */
298 int integrator_ida_solve(
299 IntegratorSystem *blsys
300 , unsigned long start_index
301 , unsigned long finish_index
302 ){
303 void *ida_mem;
304 int size, flag, t_index;
305 realtype t0, reltol, abstol, t, tret, tout1;
306 N_Vector y0, yp0, abstolvect, ypret, yret;
307 IntegratorIdaData *enginedata;
308 char *linsolver;
309 int maxl;
310
311 CONSOLE_DEBUG("STARTING IDA...");
312
313 enginedata = integrator_ida_enginedata(blsys);
314
315 enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);
316 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
317
318 /* store reference to list of relations (in enginedata) */
319 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
320 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
321 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
322 enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);
323
324 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
325 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
326 size = blsys->n_y;
327
328 if(enginedata->nrels!=size){
329 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
330 return 0; /* failure */
331 }
332
333 /* retrieve initial values from the system */
334
335 /** @TODO fix this, the starting time != first sample */
336 t0 = integrator_get_t(blsys);
337 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
338
339 CONSOLE_DEBUG("RETRIEVING y0");
340
341 y0 = N_VNew_Serial(size);
342 integrator_get_y(blsys,NV_DATA_S(y0));
343
344 CONSOLE_DEBUG("RETRIEVING yp0");
345
346 yp0 = N_VNew_Serial(size);
347 integrator_get_ydot(blsys,NV_DATA_S(yp0));
348
349 N_VPrint_Serial(yp0);
350 CONSOLE_DEBUG("yp0 is at %p",&yp0);
351
352 /* create IDA object */
353 ida_mem = IDACreate();
354
355 /* relative error tolerance */
356 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
357 CONSOLE_DEBUG("rtol = %8.2e",reltol);
358
359 /* allocate internal memory */
360 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
361 /* vector of absolute tolerances */
362 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
363 abstolvect = N_VNew_Serial(size);
364 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
365
366 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
367
368 N_VDestroy_Serial(abstolvect);
369 }else{
370 /* scalar absolute tolerance (one value for all) */
371 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
372 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
373 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
374 }
375
376 if(flag==IDA_MEM_NULL){
377 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
378 return 0;
379 }else if(flag==IDA_MEM_FAIL){
380 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
381 return 0;
382 }else if(flag==IDA_ILL_INPUT){
383 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
384 return 0;
385 }/* else success */
386
387 /* set optional inputs... */
388 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
389 IDASetRdata(ida_mem, (void *)blsys);
390 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
391 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
392 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
393 if(integrator_get_minstep(blsys)>0){
394 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
395 }
396 /* there's no capability for setting *minimum* step size in IDA */
397
398
399 /* attach linear solver module, using the default value of maxl */
400 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
401 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
402 if(strcmp(linsolver,"DENSE")==0){
403 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
404 flag = IDADense(ida_mem, size);
405 switch(flag){
406 case IDADENSE_SUCCESS: break;
407 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 0;
408 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 0;
409 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
410 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
411 }
412
413 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
414 CONSOLE_DEBUG("USING AUTODIFF");
415 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
416 switch(flag){
417 case IDADENSE_SUCCESS: break;
418 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 0;
419 }
420 }else{
421 CONSOLE_DEBUG("USING NUMERICAL DIFF");
422 }
423 }else{
424 /* remaining methods are all SPILS */
425 CONSOLE_DEBUG("IDA SPILS");
426
427 maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
428 CONSOLE_DEBUG("maxl = %d",maxl);
429
430 if(strcmp(linsolver,"SPGMR")==0){
431 CONSOLE_DEBUG("IDA SPGMR");
432 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
433 }else if(strcmp(linsolver,"SPBCG")==0){
434 CONSOLE_DEBUG("IDA SPBCG");
435 flag = IDASpbcg(ida_mem, maxl);
436 }else if(strcmp(linsolver,"SPTFQMR")==0){
437 CONSOLE_DEBUG("IDA SPTFQMR");
438 flag = IDASptfqmr(ida_mem,maxl);
439 }else{
440 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
441 return 0;
442 }
443
444 if(flag==IDASPILS_MEM_NULL){
445 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
446 return 0;
447 }else if(flag==IDASPILS_MEM_FAIL){
448 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
449 return 0;
450 }/* else success */
451
452 /* assign the J*v function */
453 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
454 CONSOLE_DEBUG("USING AUTODIFF");
455 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
456 if(flag==IDASPILS_MEM_NULL){
457 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
458 return 0;
459 }else if(flag==IDASPILS_LMEM_NULL){
460 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
461 return 0;
462 }/* else success */
463 }else{
464 CONSOLE_DEBUG("USING NUMERICAL DIFF");
465 }
466
467 if(strcmp(linsolver,"SPGMR")==0){
468 /* select Gram-Schmidt orthogonalisation */
469 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){
470 CONSOLE_DEBUG("USING MODIFIED GS");
471 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
472 if(flag!=IDASPILS_SUCCESS){
473 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
474 return 0;
475 }
476 }else{
477 CONSOLE_DEBUG("USING CLASSICAL GS");
478 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
479 if(flag!=IDASPILS_SUCCESS){
480 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
481 return 0;
482 }
483 }
484 }
485 }
486
487 /* set linear solver optional inputs...
488 ...nothing here at the moment...
489 */
490
491 /* calculate initial conditions */
492 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_CALCIC)){
493 CONSOLE_DEBUG("Solving initial conditions (given derivatives)");
494
495 blsys->currentstep=0;
496 t_index=start_index;
497 tout1 = samplelist_get(blsys->samples, t_index);
498
499 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
500
501 /* correct initial values, given derivatives */
502 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
503 /* note the new API from version 2.3 and onwards */
504 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
505 # else
506 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
507 # endif
508
509 if(flag!=IDA_SUCCESS){
510 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
511 return 0;
512 }/* else success */
513
514 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
515 }else{
516 CONSOLE_DEBUG("Not solving initial conditions (see IDA parameter 'calcic')");
517 }
518
519 /* optionally, specify ROO-FINDING PROBLEM */
520
521 /* -- set up the IntegratorReporter */
522 integrator_output_init(blsys);
523
524 /* -- store the initial values of all the stuff */
525 integrator_output_write(blsys);
526 integrator_output_write_obs(blsys);
527
528 /* specify where the returned values should be stored */
529 yret = y0;
530 ypret = yp0;
531
532 /* advance solution in time, return values as yret and derivatives as ypret */
533 blsys->currentstep=1;
534 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
535 t = samplelist_get(blsys->samples, t_index);
536 t0 = integrator_get_t(blsys);
537 asc_assert(t > t0);
538
539 #ifdef SOLVE_DEBUG
540 CONSOLE_DEBUG("Integratoring from t0 = %f to t = %f", t0, t);
541 #endif
542
543 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
544
545 /* pass the values of everything back to the compiler */
546 integrator_set_t(blsys, (double)tret);
547 integrator_set_y(blsys, NV_DATA_S(yret));
548 integrator_set_ydot(blsys, NV_DATA_S(ypret));
549
550 if(flag<0){
551 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
552 break;
553 }
554
555 /* -- do something so that blsys knows the values of tret, yret and ypret */
556
557 /* -- store the current values of all the stuff */
558 integrator_output_write(blsys);
559 integrator_output_write_obs(blsys);
560 }
561
562 /* -- close the IntegratorReporter */
563 integrator_output_close(blsys);
564
565 if(flag < 0){
566 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
567 return 0;
568 }
569
570 /* get optional outputs */
571
572 /* free solution memory */
573 N_VDestroy_Serial(yret);
574 N_VDestroy_Serial(ypret);
575
576 /* free solver memory */
577 IDAFree(ida_mem);
578
579 /* all done */
580 return 1;
581 }
582
583 /*--------------------------------------------------
584 RESIDUALS AND JACOBIAN
585 */
586 /**
587 Function to evaluate system residuals, in the form required for IDA.
588
589 Given tt, yy and yp, we need to evaluate and return rr.
590
591 @param tt current value of indep variable (time)
592 @param yy current values of dependent variable vector
593 @param yp current values of derivatives of dependent variables
594 @param rr the output residual vector (is we're returning data to)
595 @param res_data pointer to our stuff (blsys in this case).
596
597 @return 0 on success, positive on recoverable error, and
598 negative on unrecoverable error.
599 */
600 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
601 IntegratorSystem *blsys;
602 IntegratorIdaData *enginedata;
603 int i, calc_ok, is_error;
604 struct rel_relation** relptr;
605 double resid;
606 char *relname;
607 #ifdef FEX_DEBUG
608 char *varname;
609 #endif
610
611 blsys = (IntegratorSystem *)res_data;
612 enginedata = integrator_ida_enginedata(blsys);
613
614 #ifdef FEX_DEBUG
615 /* fprintf(stderr,"\n\n"); */
616 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
617 #endif
618
619 /* pass the values of everything back to the compiler */
620 integrator_set_t(blsys, (double)tt);
621 integrator_set_y(blsys, NV_DATA_S(yy));
622 integrator_set_ydot(blsys, NV_DATA_S(yp));
623
624 if(NV_LENGTH_S(rr)!=enginedata->nrels){
625 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
626 return -1; /* unrecoverable */
627 }
628
629 /**
630 @TODO does this function (fex) do bounds checking already?
631 */
632
633 /* evaluate each residual in the rellist */
634 is_error = 0;
635 relptr = enginedata->rellist;
636
637 if(enginedata->safeeval){
638 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
639 }else{
640 #ifdef FEX_DEBUG
641 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
642 #endif
643 Asc_SignalHandlerPushDefault(SIGFPE);
644 }
645
646 if (setjmp(g_fpe_env)==0) {
647 for(i=0, relptr = enginedata->rellist;
648 i< enginedata->nrels && relptr != NULL;
649 ++i, ++relptr
650 ){
651 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
652
653 NV_Ith_S(rr,i) = resid;
654 if(!calc_ok){
655 relname = rel_make_name(blsys->system, *relptr);
656 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
657 ASC_FREE(relname);
658 /* presumable some output already made? */
659 is_error = 1;
660 }/*else{
661 CONSOLE_DEBUG("Calc OK");
662 }*/
663 }
664 }else{
665 relname = rel_make_name(blsys->system, *relptr);
666 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
667 ASC_FREE(relname);
668 is_error = 1;
669 }
670
671 if(enginedata->safeeval){
672 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
673 }else{
674 Asc_SignalHandlerPopDefault(SIGFPE);
675 }
676
677 #ifdef FEX_DEBUG
678 /* output residuals to console */
679 CONSOLE_DEBUG("RESIDUAL OUTPUT");
680 fprintf(stderr,"index\t%20s\t%20s\t%s\n","y","ydot","resid");
681 for(i=0; i<blsys->n_y; ++i){
682 varname = var_make_name(blsys->system,blsys->y[i]);
683 fprintf(stderr,"%d\t%10s=%10f\t",i,varname,NV_Ith_S(yy,i));
684 if(blsys->ydot[i]){
685 varname = var_make_name(blsys->system,blsys->ydot[i]);
686 fprintf(stderr,"%10s=%10f\t",varname,NV_Ith_S(yp,i));
687 }else{
688 fprintf(stderr,"diff(%4s)=%10f\t",varname,NV_Ith_S(yp,i));
689 }
690 ASC_FREE(varname);
691 relname = rel_make_name(blsys->system,enginedata->rellist[i]);
692 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
693 }
694 #endif
695
696 if(is_error){
697 return 1;
698 }
699
700 #ifdef FEX_DEBUG
701 CONSOLE_DEBUG("RESIDUAL OK");
702 #endif
703 return 0;
704 }
705
706 /**
707 Dense Jacobian evaluation. Only suitable for small problems!
708 */
709 int integrator_ida_djex(long int Neq, realtype tt
710 , N_Vector yy, N_Vector yp, N_Vector rr
711 , realtype c_j, void *jac_data, DenseMat Jac
712 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
713 ){
714 IntegratorSystem *blsys;
715 IntegratorIdaData *enginedata;
716 char *relname;
717 #ifdef JEX_DEBUG
718 char *varname;
719 #endif
720 int status;
721 struct rel_relation **relptr;
722 int i;
723 var_filter_t filter = {VAR_SVAR | VAR_FIXED, VAR_SVAR};
724 double *derivatives;
725 int *variables;
726 int count, j;
727 long var_yindex;
728
729 blsys = (IntegratorSystem *)jac_data;
730 enginedata = integrator_ida_enginedata(blsys);
731
732 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
733 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
734 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
735
736 /* pass the values of everything back to the compiler */
737 integrator_set_t(blsys, (double)tt);
738 integrator_set_y(blsys, NV_DATA_S(yy));
739 integrator_set_ydot(blsys, NV_DATA_S(yp));
740
741 #ifdef JEX_DEBUG
742 /* print vars */
743 for(i=0; i < blsys->n_y; ++i){
744 varname = var_make_name(blsys->system, blsys->y[i]);
745 CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
746 ASC_FREE(varname);
747 }
748
749 /* print derivatives */
750 for(i=0; i < blsys->n_y; ++i){
751 if(blsys->ydot[i]){
752 varname = var_make_name(blsys->system, blsys->ydot[i]);
753 CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
754 ASC_FREE(varname);
755 }else{
756 varname = var_make_name(blsys->system, blsys->y[i]);
757 CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
758 ASC_FREE(varname);
759 }
760 }
761
762 /* print step size */
763 CONSOLE_DEBUG("<c_j> = %f",c_j);
764 #endif
765
766 /* build up the dense jacobian matrix... */
767 status = 0;
768 for(i=0, relptr = enginedata->rellist;
769 i< enginedata->nrels && relptr != NULL;
770 ++i, ++relptr
771 ){
772 /* get derivatives for this particular relation */
773 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
774
775 if(status){
776 relname = rel_make_name(blsys->system, *relptr);
777 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
778 ASC_FREE(relname);
779 break;
780 }
781
782 /* output what's going on here ... */
783 #ifdef JEX_DEBUG
784 relname = rel_make_name(blsys->system, *relptr);
785 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
786 fprintf(stderr,"%d: '%s': ",i,relname);
787 ASC_FREE(relname);
788 for(j=0;j<count;++j){
789 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
790 var_yindex = blsys->y_id[variables[j]];
791 if(var_yindex >=0){
792 fprintf(stderr," var[%d]='%s'=y[%ld]",variables[j],varname,var_yindex);
793 }else{
794 fprintf(stderr," var[%d]='%s'=ydot[%ld]",variables[j],varname,-var_yindex-1);
795 }
796 ASC_FREE(varname);
797 }
798 fprintf(stderr,"\n");
799 #endif
800 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
801 for(j=0; j < count; ++j){
802 var_yindex = blsys->y_id[variables[j]];
803 #ifndef __WIN32__
804 /* the SUNDIALS headers seem not to store 'N' on Windows */
805 ASC_ASSERT_RANGE(var_yindex, -Jac->N, Jac->N);
806 #endif
807 if(var_yindex >= 0){
808 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
809 DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
810 }else{
811 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
812 DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
813 }
814 }
815 }
816
817 #ifdef JEX_DEBUG
818 CONSOLE_DEBUG("PRINTING JAC");
819 fprintf(stderr,"\t");
820 for(j=0; j < blsys->n_y; ++j){
821 if(j)fprintf(stderr,"\t");
822 varname = var_make_name(blsys->system,blsys->y[j]);
823 fprintf(stderr,"%11s",varname);
824 ASC_FREE(varname);
825 }
826 fprintf(stderr,"\n");
827 for(i=0; i < enginedata->nrels; ++i){
828 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
829 fprintf(stderr,"%s\t",relname);
830 ASC_FREE(relname);
831
832 for(j=0; j < blsys->n_y; ++j){
833 if(j)fprintf(stderr,"\t");
834 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
835 }
836 fprintf(stderr,"\n");
837 }
838 #endif
839
840 if(status){
841 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
842 return 1;
843 }
844
845 #ifdef FEX_DEBUG
846 CONSOLE_DEBUG("DJEX RETURNING 0");
847 #endif
848 return 0;
849 }
850
851 /**
852 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
853
854 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
855
856 @param tt current value of the independent variable (time, t)
857 @param yy current value of the dependent variable vector, y(t).
858 @param yp current value of y'(t).
859 @param rr current value of the residual vector F(t, y, y').
860 @param v the vector by which the Jacobian must be multiplied to the right.
861 @param Jv the output vector computed
862 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
863 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
864 @param tmp1 @see tmp2
865 @param tmp2 (as well as tmp1) pointers to memory allocated for variables of type N_Vector for use here as temporary storage or work space.
866 @return 0 on success
867 */
868 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
869 , N_Vector v, N_Vector Jv, realtype c_j
870 , void *jac_data, N_Vector tmp1, N_Vector tmp2
871 ){
872 IntegratorSystem *blsys;
873 IntegratorIdaData *enginedata;
874 int i, j, is_error=0;
875 struct rel_relation** relptr;
876 char *relname;
877 int status;
878 double Jv_i;
879 long var_yindex;
880
881 int *variables;
882 double *derivatives;
883 var_filter_t filter;
884 int count;
885 #ifdef JEX_DEBUG
886 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
887 #endif
888
889 blsys = (IntegratorSystem *)jac_data;
890 enginedata = integrator_ida_enginedata(blsys);
891
892 /* pass the values of everything back to the compiler */
893 integrator_set_t(blsys, (double)tt);
894 integrator_set_y(blsys, NV_DATA_S(yy));
895 integrator_set_ydot(blsys, NV_DATA_S(yp));
896 /* no real use for residuals (rr) here, I don't think? */
897
898 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
899
900 i = NV_LENGTH_S(yy) * 2;
901 #ifdef JEX_DEBUG
902 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
903 #endif
904 variables = ASC_NEW_ARRAY(int, i);
905 derivatives = ASC_NEW_ARRAY(double, i);
906
907 /* evaluate the derivatives... */
908 /* J = dG_dy = dF_dy + alpha * dF_dyp */
909
910 filter.matchbits = VAR_SVAR | VAR_FIXED;
911 filter.matchvalue = VAR_SVAR;
912
913 Asc_SignalHandlerPushDefault(SIGFPE);
914 if (setjmp(g_fpe_env)==0) {
915 for(i=0, relptr = enginedata->rellist;
916 i< enginedata->nrels && relptr != NULL;
917 ++i, ++relptr
918 ){
919 /* get derivatives for this particular relation */
920 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
921 #ifdef JEX_DEBUG
922 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
923 #endif
924
925 if(status){
926 relname = rel_make_name(blsys->system, *relptr);
927 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
928 ASC_FREE(relname);
929 is_error = 1;
930 break;
931 }
932
933 /*
934 Now we have the derivatives wrt each alg/diff variable in the
935 present equation. variables[] points into the varlist. need
936 a mapping from the varlist to the y and ydot lists.
937 */
938
939 Jv_i = 0;
940 for(j=0; j < count; ++j){
941 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
942 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
943 if(varname){
944 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
945 ASC_FREE(varname);
946 }else{
947 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
948 }
949 */
950
951 /* we don't calculate derivatives wrt indep var */
952 asc_assert(variables[j]>=0);
953 if(enginedata->varlist[variables[j]] == blsys->x) continue;
954
955 var_yindex = blsys->y_id[variables[j]];
956 #ifdef JEX_DEBUG
957 CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %ld",j,variables[j],var_yindex);
958 #endif
959
960 ASC_ASSERT_RANGE(-var_yindex-1, -NV_LENGTH_S(v),NV_LENGTH_S(v));
961
962 if(var_yindex >= 0){
963 #ifdef JEX_DEBUG
964 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
965 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%ld] = %f, v[%ld] = %f)\n", i
966 , derivatives[j] * NV_Ith_S(v,var_yindex)
967 , i, var_yindex, derivatives[j]
968 , var_yindex, NV_Ith_S(v,var_yindex)
969 );
970 #endif
971 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
972 }else{
973 #ifdef JEX_DEBUG
974 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%ld] = %f, v[%ld] = %f)\n", i
975 , derivatives[j] * NV_Ith_S(v,-var_yindex-1)
976 , i, -var_yindex-1, derivatives[j]
977 , -var_yindex-1, NV_Ith_S(v,-var_yindex-1)
978 );
979 #endif
980 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
981 Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
982 }
983 }
984
985 NV_Ith_S(Jv,i) = Jv_i;
986 #ifdef JEX_DEBUG
987 CONSOLE_DEBUG("rel = %p",*relptr);
988 relname = rel_make_name(blsys->system, *relptr);
989 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
990 //ASC_FREE(relname);
991 return 1;
992 #endif
993 }
994 }else{
995 relname = rel_make_name(blsys->system, *relptr);
996 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
997 ASC_FREE(relname);
998 is_error = 1;
999 }
1000 Asc_SignalHandlerPopDefault(SIGFPE);
1001
1002 if(is_error){
1003 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1004 return 1;
1005 }
1006 return 0;
1007 }
1008
1009 /*----------------------------------------------
1010 ERROR REPORTING
1011 */
1012 /**
1013 Error message reporter function to be passed to IDA. All error messages
1014 will trigger a call to this function, so we should find everything
1015 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1016 panel (in the case of PyGTK).
1017 */
1018 void integrator_ida_error(int error_code
1019 , const char *module, const char *function
1020 , char *msg, void *eh_data
1021 ){
1022 IntegratorSystem *blsys;
1023 error_severity_t sev;
1024
1025 /* cast back the IntegratorSystem, just in case we need it */
1026 blsys = (IntegratorSystem *)eh_data;
1027
1028 /* severity depends on the sign of the error_code value */
1029 if(error_code <= 0){
1030 sev = ASC_PROG_ERR;
1031 }else{
1032 sev = ASC_PROG_WARNING;
1033 }
1034
1035 /* use our all-purpose error reporting to get stuff back to the GUI */
1036 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1037 }

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