/[ascend]/trunk/base/generic/solver/ida.c
ViewVC logotype

Annotation of /trunk/base/generic/solver/ida.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 945 - (hide annotations) (download) (as text)
Sat Nov 25 12:41:03 2006 UTC (14 years, 2 months 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 johnpye 669 /* 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 johnpye 912 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 johnpye 669 */
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 johnpye 942 #include <utilities/ascPanic.h>
43 johnpye 669 #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 johnpye 719 #include "relman.h"
55 johnpye 669
56     /* SUNDIALS includes */
57     #ifdef ASC_WITH_IDA
58 johnpye 913 # include <sundials/sundials_config.h>
59 johnpye 945 # include <sundials/sundials_dense.h>
60 johnpye 782 # include <ida/ida.h>
61     # include <nvector/nvector_serial.h>
62     # include <ida/ida_spgmr.h>
63 johnpye 945 # include <ida/ida_dense.h>
64 johnpye 669 # ifndef IDA_SUCCESS
65     # error "Failed to include SUNDIALS IDA header file"
66     # endif
67     #endif
68    
69 johnpye 913 /*
70     for the benefit of build tools that didn't sniff the SUNDIALS version, we
71 johnpye 924 assume version 2.2.x (and thence possible errors).
72 johnpye 913 */
73     #ifndef SUNDIALS_VERSION_MINOR
74 johnpye 921 # ifdef __GNUC__
75     # warning "GUESSING SUNDIALS VERSION 2.2"
76     # endif
77 johnpye 913 # define SUNDIALS_VERSION_MINOR 2
78     #endif
79     #ifndef SUNDIALS_VERSION_MAJOR
80     # define SUNDIALS_VERSION_MAJOR 2
81     #endif
82    
83 johnpye 669 /* 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 johnpye 912 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
95 johnpye 669 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 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
104 johnpye 669
105 johnpye 909 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 johnpye 669 /* 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 johnpye 945 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 johnpye 669 /*-------------------------------------------------------------
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 johnpye 909 enginedata->varlist = NULL;
131 johnpye 669 enginedata->safeeval = 1;
132 ben.allan 704 blsys->enginedata = (void *)enginedata;
133 johnpye 942 integrator_ida_params_default(blsys);
134 ben.allan 704 }
135 johnpye 669
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 johnpye 942 PARAMETERS FOR IDA
155     */
156    
157     enum ida_parameters{
158 johnpye 945 IDA_PARAM_LINSOLVER
159     ,IDA_PARAM_AUTODIFF
160 johnpye 944 ,IDA_PARAM_RTOL
161     ,IDA_PARAM_ATOL
162     ,IDA_PARAM_ATOLVECT
163 johnpye 942 ,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 johnpye 944 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 johnpye 945 }, 1e-4, DBL_MIN, DBL_MAX }
222 johnpye 944 );
223    
224 johnpye 945 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 johnpye 942 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 johnpye 669 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 johnpye 719 int size, flag, t_index;
250 johnpye 944 realtype t0, reltol, abstol, t, tret, tout1;
251     N_Vector y0, yp0, abstolvect, ypret, yret;
252 johnpye 669 IntegratorIdaData *enginedata;
253 johnpye 945 char *linsolver;
254 johnpye 669
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 johnpye 909 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
264 johnpye 669 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
265 johnpye 719 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
266 johnpye 669 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 johnpye 944
275     /** @TODO fix this, the starting time != first sample */
276     t0 = integrator_get_t(blsys);
277     CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
278 johnpye 669
279 johnpye 928 CONSOLE_DEBUG("RETRIEVING y0");
280    
281 johnpye 669 y0 = N_VNew_Serial(size);
282     integrator_get_y(blsys,NV_DATA_S(y0));
283    
284 johnpye 928 CONSOLE_DEBUG("RETRIEVING yp0");
285    
286 johnpye 669 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 johnpye 944 /* relative error tolerance */
296     reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
297 johnpye 669
298     /* allocate internal memory */
299 johnpye 944 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 johnpye 669 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 johnpye 912 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
333    
334 johnpye 669 /* attach linear solver module, using the default value of maxl */
335 johnpye 945 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
336     if(strcmp(linsolver,"SPGMR")==0){
337     flag = IDASpgmr(ida_mem, 0);
338 johnpye 942 if(flag==IDASPILS_MEM_NULL){
339     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
340     return 0;
341 johnpye 945 }else if(flag==IDASPILS_MEM_FAIL){
342     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
343 johnpye 942 return 0;
344     }/* else success */
345 johnpye 945
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 johnpye 942 }else{
382 johnpye 945 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
383     return 0;
384     }
385 johnpye 669
386 johnpye 909 /* set linear solver optional inputs...
387 johnpye 669
388 johnpye 909 ...nothing here at the moment...
389    
390     */
391    
392 johnpye 669 /* correct initial values, given derivatives */
393     blsys->currentstep=0;
394 johnpye 944 t_index=start_index;
395 johnpye 669 tout1 = samplelist_get(blsys->samples, t_index);
396 johnpye 913
397 johnpye 944 /* CONSOLE_DEBUG("Giving t value %f to IDACalcIC", tout1);*/
398    
399 johnpye 913 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
400 johnpye 944 /* note the new API from version 2.3 and onwards */
401 johnpye 913 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
402     #else
403 johnpye 669 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
404 johnpye 913 #endif
405    
406 johnpye 669 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 johnpye 942 /* CONSOLE_DEBUG("SOLVING UP TO t = %f", t); */
432 johnpye 669
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 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
491 johnpye 669 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 johnpye 903 /* fprintf(stderr,"\n\n"); */
502     /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
503 johnpye 669
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 johnpye 903 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
514 johnpye 669 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 johnpye 903 /* if(calc_ok){
529 johnpye 669 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
530     }else{
531     CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
532 johnpye 903 }*/
533 johnpye 669 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 johnpye 903 /**
551 johnpye 945 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 johnpye 912 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
590 johnpye 909
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 johnpye 903 */
605 johnpye 909 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 johnpye 928 int var_yindex;
617 johnpye 903
618 johnpye 909 int *variables;
619     double *derivatives;
620     var_filter_t filter;
621     int count;
622    
623 johnpye 944 /* fprintf(stderr,"\n--------------\n"); */
624     /* CONSOLE_DEBUG("EVALUTING JACOBIAN..."); */
625 johnpye 909
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 johnpye 912 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
636 johnpye 909 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 johnpye 944 /* CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v)); */
646     /* for(i=0; i<NV_LENGTH_S(v); ++i){
647 johnpye 930 CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));
648 johnpye 944 }*/
649 johnpye 912
650 johnpye 909 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 johnpye 944 /* fprintf(stderr,"\n"); */
657 johnpye 928 relname = rel_make_name(blsys->system, *relptr);
658 johnpye 944 /* CONSOLE_DEBUG("RELATION %d '%s'",i,relname); */
659 johnpye 928 ASC_FREE(relname);
660    
661 johnpye 909 /* get derivatives for this particular relation */
662     status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
663 johnpye 944 /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
664 johnpye 928
665 johnpye 924 for(j=0;j<count;++j){
666 johnpye 928 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
667 johnpye 944 /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */
668 johnpye 917 ASC_FREE(varname);
669     }
670 johnpye 909
671     if(!status){
672 johnpye 944 /* CONSOLE_DEBUG("Derivatives for relation %d OK",i); */
673 johnpye 909 }else{
674 johnpye 928 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
675 johnpye 909 break;
676     }
677    
678 johnpye 928 /*
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 johnpye 909 Jv_i = 0;
685     for(j=0; j < count; ++j){
686 johnpye 910 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
687 johnpye 909 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
688     if(varname){
689 johnpye 944 /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]); */
690 johnpye 909 ASC_FREE(varname);
691     }else{
692 johnpye 912 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
693 johnpye 909 }
694 johnpye 912
695 johnpye 928 var_yindex = blsys->y_id[variables[j]];
696 johnpye 944 /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
697 johnpye 928
698     if(var_yindex >= 0){
699 johnpye 944 /* CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
700 johnpye 928 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
701     }else{
702     var_yindex = -var_yindex-1;
703 johnpye 944 /* CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
704 johnpye 928 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
705     }
706 johnpye 909 }
707    
708     NV_Ith_S(Jv,i) = Jv_i;
709 johnpye 944 /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */
710 johnpye 910
711 johnpye 909 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 johnpye 912
723    
724    
725 johnpye 909 return is_error;
726     }
727    
728 johnpye 669 /*----------------------------------------------
729     ERROR REPORTING
730 johnpye 782 */
731 johnpye 669 /**
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