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

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