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

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