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

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