/[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 985 - (show annotations) (download) (as text)
Thu Dec 21 04:06:02 2006 UTC (15 years, 6 months ago) by johnpye
File MIME type: text/x-csrc
File size: 32093 byte(s)
Trying to get idakryx example working.
Better reporting of failures from IDACalcIC.
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 typedef int IdaFlagFn(void *,int *);
298 typedef const char *IdaFlagNameFn(int);
299
300 /* return 1 on success */
301 int integrator_ida_solve(
302 IntegratorSystem *blsys
303 , unsigned long start_index
304 , unsigned long finish_index
305 ){
306 void *ida_mem;
307 int size, flag, t_index;
308 realtype t0, reltol, abstol, t, tret, tout1;
309 N_Vector y0, yp0, abstolvect, ypret, yret;
310 IntegratorIdaData *enginedata;
311 char *linsolver;
312 int maxl;
313 IdaFlagFn *flagfn;
314 IdaFlagNameFn *flagnamefn;
315
316 CONSOLE_DEBUG("STARTING IDA...");
317
318 enginedata = integrator_ida_enginedata(blsys);
319
320 enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);
321 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
322
323 /* store reference to list of relations (in enginedata) */
324 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
325 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
326 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
327 enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);
328
329 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
330 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
331 size = blsys->n_y;
332
333 if(enginedata->nrels!=size){
334 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
335 return 0; /* failure */
336 }
337
338 /* retrieve initial values from the system */
339
340 /** @TODO fix this, the starting time != first sample */
341 t0 = integrator_get_t(blsys);
342 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
343
344 CONSOLE_DEBUG("RETRIEVING y0");
345
346 y0 = N_VNew_Serial(size);
347 integrator_get_y(blsys,NV_DATA_S(y0));
348
349 CONSOLE_DEBUG("RETRIEVING yp0");
350
351 yp0 = N_VNew_Serial(size);
352 integrator_get_ydot(blsys,NV_DATA_S(yp0));
353
354 N_VPrint_Serial(yp0);
355 CONSOLE_DEBUG("yp0 is at %p",&yp0);
356
357 /* create IDA object */
358 ida_mem = IDACreate();
359
360 /* relative error tolerance */
361 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
362 CONSOLE_DEBUG("rtol = %8.2e",reltol);
363
364 /* allocate internal memory */
365 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
366 /* vector of absolute tolerances */
367 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
368 abstolvect = N_VNew_Serial(size);
369 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
370
371 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
372
373 N_VDestroy_Serial(abstolvect);
374 }else{
375 /* scalar absolute tolerance (one value for all) */
376 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
377 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
378 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
379 }
380
381 if(flag==IDA_MEM_NULL){
382 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
383 return 0;
384 }else if(flag==IDA_MEM_FAIL){
385 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
386 return 0;
387 }else if(flag==IDA_ILL_INPUT){
388 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
389 return 0;
390 }/* else success */
391
392 /* set optional inputs... */
393 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
394 IDASetRdata(ida_mem, (void *)blsys);
395 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
396 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
397 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
398 if(integrator_get_minstep(blsys)>0){
399 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
400 }
401 /* there's no capability for setting *minimum* step size in IDA */
402
403
404 /* attach linear solver module, using the default value of maxl */
405 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
406 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
407 if(strcmp(linsolver,"DENSE")==0){
408 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
409 flag = IDADense(ida_mem, size);
410 switch(flag){
411 case IDADENSE_SUCCESS: break;
412 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 0;
413 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 0;
414 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
415 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
416 }
417
418 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
419 CONSOLE_DEBUG("USING AUTODIFF");
420 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
421 switch(flag){
422 case IDADENSE_SUCCESS: break;
423 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 0;
424 }
425 }else{
426 CONSOLE_DEBUG("USING NUMERICAL DIFF");
427 }
428
429 flagfn = &IDADenseGetLastFlag;
430 flagnamefn = &IDADenseGetReturnFlagName;
431
432 }else{
433 /* remaining methods are all SPILS */
434 CONSOLE_DEBUG("IDA SPILS");
435
436 maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
437 CONSOLE_DEBUG("maxl = %d",maxl);
438
439 if(strcmp(linsolver,"SPGMR")==0){
440 CONSOLE_DEBUG("IDA SPGMR");
441 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
442 }else if(strcmp(linsolver,"SPBCG")==0){
443 CONSOLE_DEBUG("IDA SPBCG");
444 flag = IDASpbcg(ida_mem, maxl);
445 }else if(strcmp(linsolver,"SPTFQMR")==0){
446 CONSOLE_DEBUG("IDA SPTFQMR");
447 flag = IDASptfqmr(ida_mem,maxl);
448 }else{
449 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
450 return 0;
451 }
452
453 flagfn = &IDASpilsGetLastFlag;
454 flagnamefn = &IDASpilsGetReturnFlagName;
455
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_MEM_FAIL){
460 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
461 return 0;
462 }/* else success */
463
464 /* assign the J*v function */
465 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
466 CONSOLE_DEBUG("USING AUTODIFF");
467 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
468 if(flag==IDASPILS_MEM_NULL){
469 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
470 return 0;
471 }else if(flag==IDASPILS_LMEM_NULL){
472 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
473 return 0;
474 }/* else success */
475 }else{
476 CONSOLE_DEBUG("USING NUMERICAL DIFF");
477 }
478
479 if(strcmp(linsolver,"SPGMR")==0){
480 /* select Gram-Schmidt orthogonalisation */
481 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){
482 CONSOLE_DEBUG("USING MODIFIED GS");
483 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
484 if(flag!=IDASPILS_SUCCESS){
485 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
486 return 0;
487 }
488 }else{
489 CONSOLE_DEBUG("USING CLASSICAL GS");
490 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
491 if(flag!=IDASPILS_SUCCESS){
492 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
493 return 0;
494 }
495 }
496 }
497 }
498
499 /* set linear solver optional inputs...
500 ...nothing here at the moment...
501 */
502
503 /* calculate initial conditions */
504 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_CALCIC)){
505 CONSOLE_DEBUG("Solving initial conditions (given derivatives)");
506
507 blsys->currentstep=0;
508 t_index=start_index;
509 tout1 = samplelist_get(blsys->samples, t_index);
510
511 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
512
513 /* correct initial values, given derivatives */
514 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
515 /* note the new API from version 2.3 and onwards */
516 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
517 # else
518 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
519 # endif
520
521 if(flag!=IDA_SUCCESS){
522 (flagfn)(ida_mem,&flag);
523 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC: %s)",(flagnamefn)(flag));
524 return 0;
525 }/* else success */
526
527 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
528 }else{
529 CONSOLE_DEBUG("Not solving initial conditions (see IDA parameter 'calcic')");
530 }
531
532 /* optionally, specify ROO-FINDING PROBLEM */
533
534 /* -- set up the IntegratorReporter */
535 integrator_output_init(blsys);
536
537 /* -- store the initial values of all the stuff */
538 integrator_output_write(blsys);
539 integrator_output_write_obs(blsys);
540
541 /* specify where the returned values should be stored */
542 yret = y0;
543 ypret = yp0;
544
545 /* advance solution in time, return values as yret and derivatives as ypret */
546 blsys->currentstep=1;
547 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
548 t = samplelist_get(blsys->samples, t_index);
549 t0 = integrator_get_t(blsys);
550 asc_assert(t > t0);
551
552 #ifdef SOLVE_DEBUG
553 CONSOLE_DEBUG("Integratoring from t0 = %f to t = %f", t0, t);
554 #endif
555
556 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
557
558 /* pass the values of everything back to the compiler */
559 integrator_set_t(blsys, (double)tret);
560 integrator_set_y(blsys, NV_DATA_S(yret));
561 integrator_set_ydot(blsys, NV_DATA_S(ypret));
562
563 if(flag<0){
564 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
565 break;
566 }
567
568 /* -- do something so that blsys knows the values of tret, yret and ypret */
569
570 /* -- store the current values of all the stuff */
571 integrator_output_write(blsys);
572 integrator_output_write_obs(blsys);
573 }
574
575 /* -- close the IntegratorReporter */
576 integrator_output_close(blsys);
577
578 if(flag < 0){
579 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
580 return 0;
581 }
582
583 /* get optional outputs */
584
585 /* free solution memory */
586 N_VDestroy_Serial(yret);
587 N_VDestroy_Serial(ypret);
588
589 /* free solver memory */
590 IDAFree(ida_mem);
591
592 /* all done */
593 return 1;
594 }
595
596 /*--------------------------------------------------
597 RESIDUALS AND JACOBIAN
598 */
599 /**
600 Function to evaluate system residuals, in the form required for IDA.
601
602 Given tt, yy and yp, we need to evaluate and return rr.
603
604 @param tt current value of indep variable (time)
605 @param yy current values of dependent variable vector
606 @param yp current values of derivatives of dependent variables
607 @param rr the output residual vector (is we're returning data to)
608 @param res_data pointer to our stuff (blsys in this case).
609
610 @return 0 on success, positive on recoverable error, and
611 negative on unrecoverable error.
612 */
613 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
614 IntegratorSystem *blsys;
615 IntegratorIdaData *enginedata;
616 int i, calc_ok, is_error;
617 struct rel_relation** relptr;
618 double resid;
619 char *relname;
620 #ifdef FEX_DEBUG
621 char *varname;
622 #endif
623
624 blsys = (IntegratorSystem *)res_data;
625 enginedata = integrator_ida_enginedata(blsys);
626
627 #ifdef FEX_DEBUG
628 /* fprintf(stderr,"\n\n"); */
629 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
630 #endif
631
632 /* pass the values of everything back to the compiler */
633 integrator_set_t(blsys, (double)tt);
634 integrator_set_y(blsys, NV_DATA_S(yy));
635 integrator_set_ydot(blsys, NV_DATA_S(yp));
636
637 if(NV_LENGTH_S(rr)!=enginedata->nrels){
638 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
639 return -1; /* unrecoverable */
640 }
641
642 /**
643 @TODO does this function (fex) do bounds checking already?
644 */
645
646 /* evaluate each residual in the rellist */
647 is_error = 0;
648 relptr = enginedata->rellist;
649
650 if(enginedata->safeeval){
651 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
652 }else{
653 #ifdef FEX_DEBUG
654 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
655 #endif
656 Asc_SignalHandlerPushDefault(SIGFPE);
657 }
658
659 if (setjmp(g_fpe_env)==0) {
660 for(i=0, relptr = enginedata->rellist;
661 i< enginedata->nrels && relptr != NULL;
662 ++i, ++relptr
663 ){
664 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
665
666 NV_Ith_S(rr,i) = resid;
667 if(!calc_ok){
668 relname = rel_make_name(blsys->system, *relptr);
669 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
670 ASC_FREE(relname);
671 /* presumable some output already made? */
672 is_error = 1;
673 }/*else{
674 CONSOLE_DEBUG("Calc OK");
675 }*/
676 }
677 }else{
678 relname = rel_make_name(blsys->system, *relptr);
679 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
680 ASC_FREE(relname);
681 is_error = 1;
682 }
683
684 if(enginedata->safeeval){
685 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
686 }else{
687 Asc_SignalHandlerPopDefault(SIGFPE);
688 }
689
690 #ifdef FEX_DEBUG
691 /* output residuals to console */
692 CONSOLE_DEBUG("RESIDUAL OUTPUT");
693 fprintf(stderr,"index\t%20s\t%20s\t%s\n","y","ydot","resid");
694 for(i=0; i<blsys->n_y; ++i){
695 varname = var_make_name(blsys->system,blsys->y[i]);
696 fprintf(stderr,"%d\t%10s=%10f\t",i,varname,NV_Ith_S(yy,i));
697 if(blsys->ydot[i]){
698 varname = var_make_name(blsys->system,blsys->ydot[i]);
699 fprintf(stderr,"%10s=%10f\t",varname,NV_Ith_S(yp,i));
700 }else{
701 fprintf(stderr,"diff(%4s)=%10f\t",varname,NV_Ith_S(yp,i));
702 }
703 ASC_FREE(varname);
704 relname = rel_make_name(blsys->system,enginedata->rellist[i]);
705 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
706 }
707 #endif
708
709 if(is_error){
710 return 1;
711 }
712
713 #ifdef FEX_DEBUG
714 CONSOLE_DEBUG("RESIDUAL OK");
715 #endif
716 return 0;
717 }
718
719 /**
720 Dense Jacobian evaluation. Only suitable for small problems!
721 */
722 int integrator_ida_djex(long int Neq, realtype tt
723 , N_Vector yy, N_Vector yp, N_Vector rr
724 , realtype c_j, void *jac_data, DenseMat Jac
725 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
726 ){
727 IntegratorSystem *blsys;
728 IntegratorIdaData *enginedata;
729 char *relname;
730 #ifdef JEX_DEBUG
731 char *varname;
732 #endif
733 int status;
734 struct rel_relation **relptr;
735 int i;
736 var_filter_t filter = {VAR_SVAR | VAR_FIXED, VAR_SVAR};
737 double *derivatives;
738 int *variables;
739 int count, j;
740 long var_yindex;
741
742 blsys = (IntegratorSystem *)jac_data;
743 enginedata = integrator_ida_enginedata(blsys);
744
745 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
746 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
747 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
748
749 /* pass the values of everything back to the compiler */
750 integrator_set_t(blsys, (double)tt);
751 integrator_set_y(blsys, NV_DATA_S(yy));
752 integrator_set_ydot(blsys, NV_DATA_S(yp));
753
754 #ifdef JEX_DEBUG
755 /* print vars */
756 for(i=0; i < blsys->n_y; ++i){
757 varname = var_make_name(blsys->system, blsys->y[i]);
758 CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
759 ASC_FREE(varname);
760 }
761
762 /* print derivatives */
763 for(i=0; i < blsys->n_y; ++i){
764 if(blsys->ydot[i]){
765 varname = var_make_name(blsys->system, blsys->ydot[i]);
766 CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
767 ASC_FREE(varname);
768 }else{
769 varname = var_make_name(blsys->system, blsys->y[i]);
770 CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
771 ASC_FREE(varname);
772 }
773 }
774
775 /* print step size */
776 CONSOLE_DEBUG("<c_j> = %f",c_j);
777 #endif
778
779 /* build up the dense jacobian matrix... */
780 status = 0;
781 for(i=0, relptr = enginedata->rellist;
782 i< enginedata->nrels && relptr != NULL;
783 ++i, ++relptr
784 ){
785 /* get derivatives for this particular relation */
786 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
787
788 if(status){
789 relname = rel_make_name(blsys->system, *relptr);
790 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
791 ASC_FREE(relname);
792 break;
793 }
794
795 /* output what's going on here ... */
796 #ifdef JEX_DEBUG
797 relname = rel_make_name(blsys->system, *relptr);
798 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
799 fprintf(stderr,"%d: '%s': ",i,relname);
800 ASC_FREE(relname);
801 for(j=0;j<count;++j){
802 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
803 var_yindex = blsys->y_id[variables[j]];
804 if(var_yindex >=0){
805 fprintf(stderr," var[%d]='%s'=y[%ld]",variables[j],varname,var_yindex);
806 }else{
807 fprintf(stderr," var[%d]='%s'=ydot[%ld]",variables[j],varname,-var_yindex-1);
808 }
809 ASC_FREE(varname);
810 }
811 fprintf(stderr,"\n");
812 #endif
813 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
814 for(j=0; j < count; ++j){
815 var_yindex = blsys->y_id[variables[j]];
816 #ifndef __WIN32__
817 /* the SUNDIALS headers seem not to store 'N' on Windows */
818 ASC_ASSERT_RANGE(var_yindex, -Jac->N, Jac->N);
819 #endif
820 if(var_yindex >= 0){
821 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
822 DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
823 }else{
824 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
825 DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
826 }
827 }
828 }
829
830 #ifdef JEX_DEBUG
831 CONSOLE_DEBUG("PRINTING JAC");
832 fprintf(stderr,"\t");
833 for(j=0; j < blsys->n_y; ++j){
834 if(j)fprintf(stderr,"\t");
835 varname = var_make_name(blsys->system,blsys->y[j]);
836 fprintf(stderr,"%11s",varname);
837 ASC_FREE(varname);
838 }
839 fprintf(stderr,"\n");
840 for(i=0; i < enginedata->nrels; ++i){
841 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
842 fprintf(stderr,"%s\t",relname);
843 ASC_FREE(relname);
844
845 for(j=0; j < blsys->n_y; ++j){
846 if(j)fprintf(stderr,"\t");
847 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
848 }
849 fprintf(stderr,"\n");
850 }
851 #endif
852
853 if(status){
854 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
855 return 1;
856 }
857
858 #ifdef FEX_DEBUG
859 CONSOLE_DEBUG("DJEX RETURNING 0");
860 #endif
861 return 0;
862 }
863
864 /**
865 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
866
867 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
868
869 @param tt current value of the independent variable (time, t)
870 @param yy current value of the dependent variable vector, y(t).
871 @param yp current value of y'(t).
872 @param rr current value of the residual vector F(t, y, y').
873 @param v the vector by which the Jacobian must be multiplied to the right.
874 @param Jv the output vector computed
875 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
876 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
877 @param tmp1 @see tmp2
878 @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.
879 @return 0 on success
880 */
881 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
882 , N_Vector v, N_Vector Jv, realtype c_j
883 , void *jac_data, N_Vector tmp1, N_Vector tmp2
884 ){
885 IntegratorSystem *blsys;
886 IntegratorIdaData *enginedata;
887 int i, j, is_error=0;
888 struct rel_relation** relptr;
889 char *relname;
890 int status;
891 double Jv_i;
892 long var_yindex;
893
894 int *variables;
895 double *derivatives;
896 var_filter_t filter;
897 int count;
898 #ifdef JEX_DEBUG
899 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
900 #endif
901
902 blsys = (IntegratorSystem *)jac_data;
903 enginedata = integrator_ida_enginedata(blsys);
904
905 /* pass the values of everything back to the compiler */
906 integrator_set_t(blsys, (double)tt);
907 integrator_set_y(blsys, NV_DATA_S(yy));
908 integrator_set_ydot(blsys, NV_DATA_S(yp));
909 /* no real use for residuals (rr) here, I don't think? */
910
911 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
912
913 i = NV_LENGTH_S(yy) * 2;
914 #ifdef JEX_DEBUG
915 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
916 #endif
917 variables = ASC_NEW_ARRAY(int, i);
918 derivatives = ASC_NEW_ARRAY(double, i);
919
920 /* evaluate the derivatives... */
921 /* J = dG_dy = dF_dy + alpha * dF_dyp */
922
923 filter.matchbits = VAR_SVAR | VAR_FIXED;
924 filter.matchvalue = VAR_SVAR;
925
926 Asc_SignalHandlerPushDefault(SIGFPE);
927 if (setjmp(g_fpe_env)==0) {
928 for(i=0, relptr = enginedata->rellist;
929 i< enginedata->nrels && relptr != NULL;
930 ++i, ++relptr
931 ){
932 /* get derivatives for this particular relation */
933 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
934 #ifdef JEX_DEBUG
935 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
936 #endif
937
938 if(status){
939 relname = rel_make_name(blsys->system, *relptr);
940 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
941 ASC_FREE(relname);
942 is_error = 1;
943 break;
944 }
945
946 /*
947 Now we have the derivatives wrt each alg/diff variable in the
948 present equation. variables[] points into the varlist. need
949 a mapping from the varlist to the y and ydot lists.
950 */
951
952 Jv_i = 0;
953 for(j=0; j < count; ++j){
954 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
955 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
956 if(varname){
957 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
958 ASC_FREE(varname);
959 }else{
960 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
961 }
962 */
963
964 /* we don't calculate derivatives wrt indep var */
965 asc_assert(variables[j]>=0);
966 if(enginedata->varlist[variables[j]] == blsys->x) continue;
967
968 var_yindex = blsys->y_id[variables[j]];
969 #ifdef JEX_DEBUG
970 CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %ld",j,variables[j],var_yindex);
971 #endif
972
973 ASC_ASSERT_RANGE(-var_yindex-1, -NV_LENGTH_S(v),NV_LENGTH_S(v));
974
975 if(var_yindex >= 0){
976 #ifdef JEX_DEBUG
977 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
978 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%ld] = %f, v[%ld] = %f)\n", i
979 , derivatives[j] * NV_Ith_S(v,var_yindex)
980 , i, var_yindex, derivatives[j]
981 , var_yindex, NV_Ith_S(v,var_yindex)
982 );
983 #endif
984 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
985 }else{
986 #ifdef JEX_DEBUG
987 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%ld] = %f, v[%ld] = %f)\n", i
988 , derivatives[j] * NV_Ith_S(v,-var_yindex-1)
989 , i, -var_yindex-1, derivatives[j]
990 , -var_yindex-1, NV_Ith_S(v,-var_yindex-1)
991 );
992 #endif
993 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
994 Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
995 }
996 }
997
998 NV_Ith_S(Jv,i) = Jv_i;
999 #ifdef JEX_DEBUG
1000 CONSOLE_DEBUG("rel = %p",*relptr);
1001 relname = rel_make_name(blsys->system, *relptr);
1002 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1003 //ASC_FREE(relname);
1004 return 1;
1005 #endif
1006 }
1007 }else{
1008 relname = rel_make_name(blsys->system, *relptr);
1009 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1010 ASC_FREE(relname);
1011 is_error = 1;
1012 }
1013 Asc_SignalHandlerPopDefault(SIGFPE);
1014
1015 if(is_error){
1016 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1017 return 1;
1018 }
1019 return 0;
1020 }
1021
1022 /*----------------------------------------------
1023 ERROR REPORTING
1024 */
1025 /**
1026 Error message reporter function to be passed to IDA. All error messages
1027 will trigger a call to this function, so we should find everything
1028 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1029 panel (in the case of PyGTK).
1030 */
1031 void integrator_ida_error(int error_code
1032 , const char *module, const char *function
1033 , char *msg, void *eh_data
1034 ){
1035 IntegratorSystem *blsys;
1036 error_severity_t sev;
1037
1038 /* cast back the IntegratorSystem, just in case we need it */
1039 blsys = (IntegratorSystem *)eh_data;
1040
1041 /* severity depends on the sign of the error_code value */
1042 if(error_code <= 0){
1043 sev = ASC_PROG_ERR;
1044 }else{
1045 sev = ASC_PROG_WARNING;
1046 }
1047
1048 /* use our all-purpose error reporting to get stuff back to the GUI */
1049 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1050 }

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