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

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