/[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 945 - (show annotations) (download) (as text)
Sat Nov 25 12:41:03 2006 UTC (14 years ago) by johnpye
File MIME type: text/x-csrc
File size: 23688 byte(s)
Fixed a bug in slv_param_char.
More work on 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_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 /**
89 Struct containing any stuff that IDA needs that doesn't fit into the
90 common IntegratorSystem struct.
91 */
92 typedef struct{
93 struct rel_relation **rellist; /**< NULL terminated list of rels */
94 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
95 int nrels;
96 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
97 } IntegratorIdaData;
98
99 /*-------------------------------------------------------------
100 FORWARD DECLS
101 */
102 /* residual function forward declaration */
103 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
104
105 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
106 , N_Vector v, N_Vector Jv, realtype c_j
107 , void *jac_data, N_Vector tmp1, N_Vector tmp2
108 );
109
110 /* error handler forward declaration */
111 void integrator_ida_error(int error_code
112 , const char *module, const char *function
113 , char *msg, void *eh_data
114 );
115
116 int integrator_ida_djex(long int Neq, realtype tt
117 , N_Vector yy, N_Vector yp, N_Vector rr
118 , realtype c_j, void *jac_data, DenseMat Jac
119 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
120 );
121
122 /*-------------------------------------------------------------
123 SETUP/TEARDOWN ROUTINES
124 */
125 void integrator_ida_create(IntegratorSystem *blsys){
126 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
127 IntegratorIdaData *enginedata;
128 enginedata = ASC_NEW(IntegratorIdaData);
129 enginedata->rellist = NULL;
130 enginedata->varlist = NULL;
131 enginedata->safeeval = 1;
132 blsys->enginedata = (void *)enginedata;
133 integrator_ida_params_default(blsys);
134 }
135
136 void integrator_ida_free(void *enginedata){
137 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
138 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
139 /* note, we don't own the rellist, so don't need to free it */
140 ASC_FREE(d);
141 }
142
143 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
144 IntegratorIdaData *d;
145 assert(blsys!=NULL);
146 assert(blsys->enginedata!=NULL);
147 assert(blsys->engine==INTEG_IDA);
148 d = ((IntegratorIdaData *)(blsys->enginedata));
149 assert(d->safeeval = 1);
150 return d;
151 }
152
153 /*-------------------------------------------------------------
154 PARAMETERS FOR IDA
155 */
156
157 enum ida_parameters{
158 IDA_PARAM_LINSOLVER
159 ,IDA_PARAM_AUTODIFF
160 ,IDA_PARAM_RTOL
161 ,IDA_PARAM_ATOL
162 ,IDA_PARAM_ATOLVECT
163 ,IDA_PARAMS_SIZE
164 };
165
166 /**
167 Here the full set of parameters is defined, along with upper/lower bounds,
168 etc. The values are stuck into the blsys->params structure.
169
170 @return 0 on success
171 */
172 int integrator_ida_params_default(IntegratorSystem *blsys){
173 asc_assert(blsys!=NULL);
174 asc_assert(blsys->engine==INTEG_IDA);
175 slv_parameters_t *p;
176 p = &(blsys->params);
177
178 slv_destroy_parms(p);
179
180 if(p->parms==NULL){
181 CONSOLE_DEBUG("params NULL");
182 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
183 if(p->parms==NULL)return -1;
184 p->dynamic_parms = 1;
185 }else{
186 CONSOLE_DEBUG("params not NULL");
187 }
188
189 /* reset the number of parameters to zero so that we can check it at the end */
190 p->num_parms = 0;
191
192 slv_param_bool(p,IDA_PARAM_AUTODIFF
193 ,(SlvParameterInitBool){{"autodiff"
194 ,"Use auto-diff?",1
195 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
196 }, TRUE}
197 );
198
199 slv_param_bool(p,IDA_PARAM_ATOLVECT
200 ,(SlvParameterInitBool){{"atolvect"
201 ,"Use 'ode_atol' values as specified?",1
202 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
203 " in the integration. If FALSE, a scalar absolute tolerance value"
204 " is shared by all variables. See IDA manual, section 5.5.1"
205 }, TRUE }
206 );
207
208 slv_param_real(p,IDA_PARAM_ATOL
209 ,(SlvParameterInitReal){{"atol"
210 ,"Scalar absolute error tolerance",1
211 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
212 " See IDA manual, section 5.5.1"
213 }, 1e-5, DBL_MIN, DBL_MAX }
214 );
215
216 slv_param_real(p,IDA_PARAM_RTOL
217 ,(SlvParameterInitReal){{"rtol"
218 ,"Scalar relative error tolerance",1
219 ,"Value of the scalar relative error tolerance."
220 " See IDA manual, section 5.5.1"
221 }, 1e-4, DBL_MIN, DBL_MAX }
222 );
223
224 slv_param_char(p,IDA_PARAM_LINSOLVER
225 ,(SlvParameterInitChar){{"linsolver"
226 ,"Linear solver",1
227 ,"See IDA manual, section 5.5.3."
228 }, "DENSE"}, (char *[]){"DENSE","SPGMR",NULL}
229 );
230
231 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
232
233 CONSOLE_DEBUG("Created %d params", p->num_parms);
234
235 return 0;
236 }
237
238 /*-------------------------------------------------------------
239 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
240 */
241
242 /* return 1 on success */
243 int integrator_ida_solve(
244 IntegratorSystem *blsys
245 , unsigned long start_index
246 , unsigned long finish_index
247 ){
248 void *ida_mem;
249 int size, flag, t_index;
250 realtype t0, reltol, abstol, t, tret, tout1;
251 N_Vector y0, yp0, abstolvect, ypret, yret;
252 IntegratorIdaData *enginedata;
253 char *linsolver;
254
255 CONSOLE_DEBUG("STARTING IDA...");
256
257 enginedata = integrator_ida_enginedata(blsys);
258 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
259
260 /* store reference to list of relations (in enginedata) */
261 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
262 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
263 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
264 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
265 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
266 size = blsys->n_y;
267
268 if(enginedata->nrels!=size){
269 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
270 return 0; /* failure */
271 }
272
273 /* retrieve initial values from the system */
274
275 /** @TODO fix this, the starting time != first sample */
276 t0 = integrator_get_t(blsys);
277 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
278
279 CONSOLE_DEBUG("RETRIEVING y0");
280
281 y0 = N_VNew_Serial(size);
282 integrator_get_y(blsys,NV_DATA_S(y0));
283
284 CONSOLE_DEBUG("RETRIEVING yp0");
285
286 yp0 = N_VNew_Serial(size);
287 integrator_get_ydot(blsys,NV_DATA_S(yp0));
288
289 N_VPrint_Serial(yp0);
290 CONSOLE_DEBUG("yp0 is at %p",&yp0);
291
292 /* create IDA object */
293 ida_mem = IDACreate();
294
295 /* relative error tolerance */
296 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
297
298 /* allocate internal memory */
299 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
300 /* vector of absolute tolerances */
301 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
302 abstolvect = N_VNew_Serial(size);
303 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
304
305 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
306 }else{
307 /* scalar absolute tolerance (one value for all) */
308 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
309 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
310 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
311 }
312
313 if(flag==IDA_MEM_NULL){
314 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
315 return 0;
316 }else if(flag==IDA_MEM_FAIL){
317 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
318 return 0;
319 }else if(flag==IDA_ILL_INPUT){
320 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
321 return 0;
322 }/* else success */
323
324 /* set optional inputs... */
325 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
326 IDASetRdata(ida_mem, (void *)blsys);
327 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
328 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
329 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
330 /* there's no capability for setting *minimum* step size in IDA */
331
332 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
333
334 /* attach linear solver module, using the default value of maxl */
335 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
336 if(strcmp(linsolver,"SPGMR")==0){
337 flag = IDASpgmr(ida_mem, 0);
338 if(flag==IDASPILS_MEM_NULL){
339 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
340 return 0;
341 }else if(flag==IDASPILS_MEM_FAIL){
342 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
343 return 0;
344 }/* else success */
345
346 /* assign the J*v function */
347 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
348 CONSOLE_DEBUG("USING AUTODIFF");
349 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
350 if(flag==IDASPILS_MEM_NULL){
351 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
352 return 0;
353 }else if(flag==IDASPILS_LMEM_NULL){
354 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
355 return 0;
356 }/* else success */
357 }else{
358 CONSOLE_DEBUG("USING NUMERICAL DIFF");
359 }
360 }else if(strcmp(linsolver,"DENSE")==0){
361 CONSOLE_DEBUG("USING IDADEENSE SOLVER, size = %d",size);
362 flag = IDADense(ida_mem, size);
363 switch(flag){
364 case IDADENSE_SUCCESS: break;
365 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 0;
366 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 0;
367 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
368 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
369 }
370
371 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
372 CONSOLE_DEBUG("USING AUTODIFF");
373 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
374 switch(flag){
375 case IDADENSE_SUCCESS: break;
376 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 0;
377 }
378 }else{
379 CONSOLE_DEBUG("USING NUMERICAL DIFF");
380 }
381 }else{
382 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
383 return 0;
384 }
385
386 /* set linear solver optional inputs...
387
388 ...nothing here at the moment...
389
390 */
391
392 /* correct initial values, given derivatives */
393 blsys->currentstep=0;
394 t_index=start_index;
395 tout1 = samplelist_get(blsys->samples, t_index);
396
397 /* CONSOLE_DEBUG("Giving t value %f to IDACalcIC", tout1);*/
398
399 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
400 /* note the new API from version 2.3 and onwards */
401 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
402 #else
403 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
404 #endif
405
406 if(flag!=IDA_SUCCESS){
407 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
408 return 0;
409 }/* else success */
410
411 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
412
413 /* optionally, specify ROO-FINDING PROBLEM */
414
415 /* -- set up the IntegratorReporter */
416 integrator_output_init(blsys);
417
418 /* -- store the initial values of all the stuff */
419 integrator_output_write(blsys);
420 integrator_output_write_obs(blsys);
421
422 /* specify where the returned values should be stored */
423 yret = y0;
424 ypret = yp0;
425
426 /* advance solution in time, return values as yret and derivatives as ypret */
427 blsys->currentstep=1;
428 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
429 t = samplelist_get(blsys->samples, t_index);
430
431 /* CONSOLE_DEBUG("SOLVING UP TO t = %f", t); */
432
433 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
434
435 /* pass the values of everything back to the compiler */
436 integrator_set_t(blsys, (double)tret);
437 integrator_set_y(blsys, NV_DATA_S(yret));
438 integrator_set_ydot(blsys, NV_DATA_S(ypret));
439
440 if(flag<0){
441 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
442 break;
443 }
444
445 /* -- do something so that blsys knows the values of tret, yret and ypret */
446
447 /* -- store the current values of all the stuff */
448 integrator_output_write(blsys);
449 integrator_output_write_obs(blsys);
450 }
451
452 /* -- close the IntegratorReporter */
453 integrator_output_close(blsys);
454
455 if(flag < 0){
456 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
457 return 0;
458 }
459
460 /* get optional outputs */
461
462 /* free solution memory */
463 N_VDestroy_Serial(yret);
464 N_VDestroy_Serial(ypret);
465
466 /* free solver memory */
467 IDAFree(ida_mem);
468
469 /* all done */
470 return 1;
471 }
472
473 /*--------------------------------------------------
474 RESIDUALS AND JACOBIAN
475 */
476 /**
477 Function to evaluate system residuals, in the form required for IDA.
478
479 Given tt, yy and yp, we need to evaluate and return rr.
480
481 @param tt current value of indep variable (time)
482 @param yy current values of dependent variable vector
483 @param yp current values of derivatives of dependent variables
484 @param rr the output residual vector (is we're returning data to)
485 @param res_data pointer to our stuff (blsys in this case).
486
487 @return 0 on success, positive on recoverable error, and
488 negative on unrecoverable error.
489 */
490 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
491 IntegratorSystem *blsys;
492 IntegratorIdaData *enginedata;
493 int i, calc_ok, is_error;
494 struct rel_relation** relptr;
495 double resid;
496 char *relname;
497
498 blsys = (IntegratorSystem *)res_data;
499 enginedata = integrator_ida_enginedata(blsys);
500
501 /* fprintf(stderr,"\n\n"); */
502 /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
503
504 /* pass the values of everything back to the compiler */
505 integrator_set_t(blsys, (double)tt);
506 integrator_set_y(blsys, NV_DATA_S(yy));
507 integrator_set_ydot(blsys, NV_DATA_S(yp));
508
509 /* revaluate the system residuals using the new data */
510 is_error = 0;
511 relptr = enginedata->rellist;
512
513 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
514 if(NV_LENGTH_S(rr)!=enginedata->nrels){
515 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
516 return -1; /* unrecoverable */
517 }
518
519 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
520 if (setjmp(g_fpe_env)==0) {
521 for(i=0, relptr = enginedata->rellist;
522 i< enginedata->nrels && relptr != NULL;
523 ++i, ++relptr
524 ){
525 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
526
527 relname = rel_make_name(blsys->system, *relptr);
528 /* if(calc_ok){
529 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
530 }else{
531 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
532 }*/
533 ASC_FREE(relname);
534
535 NV_Ith_S(rr,i) = resid;
536 if(!calc_ok){
537 /* presumable some output already made? */
538 is_error = 1;
539 }
540 }
541 }else{
542 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
543 }
544 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
545
546 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
547 return is_error;
548 }
549
550 /**
551 Dense Jacobian evaluation. Only suitable for small problems!
552 */
553 int integrator_ida_djex(long int Neq, realtype tt
554 , N_Vector yy, N_Vector yp, N_Vector rr
555 , realtype c_j, void *jac_data, DenseMat Jac
556 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
557 ){
558 IntegratorSystem *blsys;
559 IntegratorIdaData *enginedata;
560 realtype *yval;
561
562 blsys = (IntegratorSystem *)jac_data;
563 enginedata = integrator_ida_enginedata(blsys);
564
565 /* pass the values of everything back to the compiler */
566 integrator_set_t(blsys, (double)tt);
567 integrator_set_y(blsys, NV_DATA_S(yy));
568 integrator_set_ydot(blsys, NV_DATA_S(yp));
569
570 yval = NV_DATA_S(yy);
571
572 /* AWFUL HACK! In an attempt to prove that 'IDADense' works as expected,
573 I've pulled the jacobians straight from idadenx.c (an IDA example) */
574 # define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1)
575 IJth(Jac,1,1) = RCONST(-0.04) - c_j;
576 IJth(Jac,2,1) = RCONST(0.04);
577 IJth(Jac,3,1) = 1.0;
578 IJth(Jac,1,2) = RCONST(1.0e4)*yval[2];
579 IJth(Jac,2,2) = RCONST(-1.0e4)*yval[2] - RCONST(6.0e7)*yval[1] - c_j;
580 IJth(Jac,3,2) = 1.0;
581 IJth(Jac,1,3) = RCONST(1.0e4)*yval[1];
582 IJth(Jac,2,3) = RCONST(-1.0e4)*yval[1];
583 IJth(Jac,3,3) = 1.0;
584 #undef IJth
585 return(0);
586 }
587
588 /**
589 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
590
591 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
592
593 @param tt current value of the independent variable (time, t)
594 @param yy current value of the dependent variable vector, y(t).
595 @param yp current value of y'(t).
596 @param rr current value of the residual vector F(t, y, y').
597 @param v the vector by which the Jacobian must be multiplied to the right.
598 @param Jv the output vector computed
599 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
600 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
601 @param tmp1 @see tmp2
602 @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.
603 @return 0 on success
604 */
605 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
606 , N_Vector v, N_Vector Jv, realtype c_j
607 , void *jac_data, N_Vector tmp1, N_Vector tmp2
608 ){
609 IntegratorSystem *blsys;
610 IntegratorIdaData *enginedata;
611 int i, j, is_error=0;
612 struct rel_relation** relptr;
613 char *relname, *varname;
614 int status;
615 double Jv_i;
616 int var_yindex;
617
618 int *variables;
619 double *derivatives;
620 var_filter_t filter;
621 int count;
622
623 /* fprintf(stderr,"\n--------------\n"); */
624 /* CONSOLE_DEBUG("EVALUTING JACOBIAN..."); */
625
626 blsys = (IntegratorSystem *)jac_data;
627 enginedata = integrator_ida_enginedata(blsys);
628
629 /* pass the values of everything back to the compiler */
630 integrator_set_t(blsys, (double)tt);
631 integrator_set_y(blsys, NV_DATA_S(yy));
632 integrator_set_ydot(blsys, NV_DATA_S(yp));
633 /* no real use for residuals (rr) here, I don't think? */
634
635 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
636 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
637 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
638
639 /* evaluate the derivatives... */
640 /* J = dG_dy = dF_dy + alpha * dF_dyp */
641
642 filter.matchbits = VAR_SVAR;
643 filter.matchvalue = VAR_SVAR;
644
645 /* CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v)); */
646 /* for(i=0; i<NV_LENGTH_S(v); ++i){
647 CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));
648 }*/
649
650 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
651 if (setjmp(g_fpe_env)==0) {
652 for(i=0, relptr = enginedata->rellist;
653 i< enginedata->nrels && relptr != NULL;
654 ++i, ++relptr
655 ){
656 /* fprintf(stderr,"\n"); */
657 relname = rel_make_name(blsys->system, *relptr);
658 /* CONSOLE_DEBUG("RELATION %d '%s'",i,relname); */
659 ASC_FREE(relname);
660
661 /* get derivatives for this particular relation */
662 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
663 /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
664
665 for(j=0;j<count;++j){
666 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
667 /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */
668 ASC_FREE(varname);
669 }
670
671 if(!status){
672 /* CONSOLE_DEBUG("Derivatives for relation %d OK",i); */
673 }else{
674 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
675 break;
676 }
677
678 /*
679 Now we have the derivatives wrt each alg/diff variable in the
680 present equation. variables[] points into the varlist. need
681 a mapping from the varlist to the y and ydot lists.
682 */
683
684 Jv_i = 0;
685 for(j=0; j < count; ++j){
686 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
687 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
688 if(varname){
689 /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]); */
690 ASC_FREE(varname);
691 }else{
692 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
693 }
694
695 var_yindex = blsys->y_id[variables[j]];
696 /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
697
698 if(var_yindex >= 0){
699 /* CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
700 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
701 }else{
702 var_yindex = -var_yindex-1;
703 /* CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
704 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
705 }
706 }
707
708 NV_Ith_S(Jv,i) = Jv_i;
709 /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */
710
711 if(status){
712 /* presumably some error_reporter will already have been made*/
713 is_error = 1;
714 }
715 }
716 }else{
717 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
718 }
719 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
720
721 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
722
723
724
725 return is_error;
726 }
727
728 /*----------------------------------------------
729 ERROR REPORTING
730 */
731 /**
732 Error message reporter function to be passed to IDA. All error messages
733 will trigger a call to this function, so we should find everything
734 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
735 panel (in the case of PyGTK).
736 */
737 void integrator_ida_error(int error_code
738 , const char *module, const char *function
739 , char *msg, void *eh_data
740 ){
741 IntegratorSystem *blsys;
742 error_severity_t sev;
743
744 /* cast back the IntegratorSystem, just in case we need it */
745 blsys = (IntegratorSystem *)eh_data;
746
747 /* severity depends on the sign of the error_code value */
748 if(error_code <= 0){
749 sev = ASC_PROG_ERR;
750 }else{
751 sev = ASC_PROG_WARNING;
752 }
753
754 /* use our all-purpose error reporting to get stuff back to the GUI */
755 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
756 }

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