/[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 913 - (hide annotations) (download) (as text)
Sat Oct 28 03:55:19 2006 UTC (13 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 16869 byte(s)
Added test for SUNDIALS version (2.2.1 and 2.3.0-pre are preferred)
Fixed silly warning message about ignored return values from CONSOLE_DEBUG on GCC.
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     #include <compiler/instance_enum.h>
43     #include "var.h"
44     #include "rel.h"
45     #include "discrete.h"
46     #include "conditional.h"
47     #include "logrel.h"
48     #include "bnd.h"
49     #include "linsol.h"
50     #include "linsolqr.h"
51     #include "slv_common.h"
52     #include "slv_client.h"
53 johnpye 719 #include "relman.h"
54 johnpye 669
55     /* SUNDIALS includes */
56     #ifdef ASC_WITH_IDA
57 johnpye 913 # include <sundials/sundials_config.h>
58 johnpye 782 # include <ida/ida.h>
59     # include <nvector/nvector_serial.h>
60     # include <ida/ida_spgmr.h>
61 johnpye 669 # ifndef IDA_SUCCESS
62     # error "Failed to include SUNDIALS IDA header file"
63     # endif
64     #endif
65    
66 johnpye 913 /*
67     for the benefit of build tools that didn't sniff the SUNDIALS version, we
68     assume version 2.2.x (and thence possible errors)
69     */
70     #ifndef SUNDIALS_VERSION_MINOR
71     # warning "GUESSING SUNDIALS VERSION 2.2"
72     # define SUNDIALS_VERSION_MINOR 2
73     #endif
74     #ifndef SUNDIALS_VERSION_MAJOR
75     # define SUNDIALS_VERSION_MAJOR 2
76     #endif
77    
78 johnpye 669 /* check that we've got what we expect now */
79     #ifndef ASC_IDA_H
80     # error "Failed to include ASCEND IDA header file"
81     #endif
82    
83     /**
84     Struct containing any stuff that IDA needs that doesn't fit into the
85     common IntegratorSystem struct.
86     */
87     typedef struct{
88     struct rel_relation **rellist; /**< NULL terminated list of rels */
89 johnpye 912 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
90 johnpye 669 int nrels;
91     int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
92     } IntegratorIdaData;
93    
94     /*-------------------------------------------------------------
95     FORWARD DECLS
96     */
97     /* residual function forward declaration */
98 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
99 johnpye 669
100 johnpye 909 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
101     , N_Vector v, N_Vector Jv, realtype c_j
102     , void *jac_data, N_Vector tmp1, N_Vector tmp2
103     );
104    
105 johnpye 669 /* error handler forward declaration */
106     void integrator_ida_error(int error_code
107     , const char *module, const char *function
108     , char *msg, void *eh_data
109     );
110    
111     /*-------------------------------------------------------------
112     SETUP/TEARDOWN ROUTINES
113     */
114     void integrator_ida_create(IntegratorSystem *blsys){
115     CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
116     IntegratorIdaData *enginedata;
117     enginedata = ASC_NEW(IntegratorIdaData);
118     enginedata->rellist = NULL;
119 johnpye 909 enginedata->varlist = NULL;
120 johnpye 669 enginedata->safeeval = 1;
121 ben.allan 704 blsys->enginedata = (void *)enginedata;
122     }
123 johnpye 669
124     void integrator_ida_free(void *enginedata){
125     CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
126     IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
127     /* note, we don't own the rellist, so don't need to free it */
128     ASC_FREE(d);
129     }
130    
131     IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
132     IntegratorIdaData *d;
133     assert(blsys!=NULL);
134     assert(blsys->enginedata!=NULL);
135     assert(blsys->engine==INTEG_IDA);
136     d = ((IntegratorIdaData *)(blsys->enginedata));
137     assert(d->safeeval = 1);
138     return d;
139     }
140    
141     /*-------------------------------------------------------------
142     MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
143     */
144    
145     /* return 1 on success */
146     int integrator_ida_solve(
147     IntegratorSystem *blsys
148     , unsigned long start_index
149     , unsigned long finish_index
150     ){
151     void *ida_mem;
152 johnpye 719 int size, flag, t_index;
153 johnpye 669 realtype t0, reltol, t, tret, tout1;
154     N_Vector y0, yp0, abstol, ypret, yret;
155     IntegratorIdaData *enginedata;
156    
157     CONSOLE_DEBUG("STARTING IDA...");
158    
159     enginedata = integrator_ida_enginedata(blsys);
160     CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
161    
162     /* store reference to list of relations (in enginedata) */
163     enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
164     enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
165 johnpye 909 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
166 johnpye 669 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
167 johnpye 719 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
168 johnpye 669 size = blsys->n_y;
169    
170     if(enginedata->nrels!=size){
171     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
172     return 0; /* failure */
173     }
174    
175     /* retrieve initial values from the system */
176     t0 = samplelist_get(blsys->samples,start_index);
177    
178     y0 = N_VNew_Serial(size);
179     integrator_get_y(blsys,NV_DATA_S(y0));
180    
181     yp0 = N_VNew_Serial(size);
182     integrator_get_ydot(blsys,NV_DATA_S(yp0));
183    
184     N_VPrint_Serial(yp0);
185     CONSOLE_DEBUG("yp0 is at %p",&yp0);
186    
187     /* create IDA object */
188     ida_mem = IDACreate();
189    
190     /* retrieve the absolute tolerance values for each variable */
191     abstol = N_VNew_Serial(size);
192     N_VConst(0.1,abstol); /** @TODO fill in the abstol values from the model */
193     reltol = 0.001;
194    
195     /* allocate internal memory */
196 johnpye 903 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstol);
197 johnpye 669 if(flag==IDA_MEM_NULL){
198     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
199     return 0;
200     }else if(flag==IDA_MEM_FAIL){
201     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
202     return 0;
203     }else if(flag==IDA_ILL_INPUT){
204     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
205     return 0;
206     }/* else success */
207    
208     /* set optional inputs... */
209     IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
210     IDASetRdata(ida_mem, (void *)blsys);
211     IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
212     IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
213     IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
214     /* there's no capability for setting *minimum* step size in IDA */
215    
216 johnpye 912 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
217    
218 johnpye 669 /* attach linear solver module, using the default value of maxl */
219     flag = IDASpgmr(ida_mem, 0);
220     if(flag==IDASPILS_MEM_NULL){
221     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
222     return 0;
223     }else if(flag==IDASPILS_MEM_FAIL){
224     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
225     return 0;
226     }/* else success */
227    
228 johnpye 909 /* assign the J*v function */
229     flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
230     if(flag==IDASPILS_MEM_NULL){
231     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
232     return 0;
233     }else if(flag==IDASPILS_LMEM_NULL){
234     ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
235     return 0;
236     }/* else success */
237 johnpye 669
238 johnpye 909 /* set linear solver optional inputs...
239 johnpye 669
240 johnpye 909 ...nothing here at the moment...
241    
242     */
243    
244 johnpye 669 /* correct initial values, given derivatives */
245     blsys->currentstep=0;
246     t_index=start_index+1;
247     tout1 = samplelist_get(blsys->samples, t_index);
248 johnpye 913
249     #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
250     flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
251     #else
252 johnpye 669 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
253 johnpye 913 #endif
254    
255 johnpye 669 if(flag!=IDA_SUCCESS){
256     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
257     return 0;
258     }/* else success */
259    
260     CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
261    
262     /* optionally, specify ROO-FINDING PROBLEM */
263    
264     /* -- set up the IntegratorReporter */
265     integrator_output_init(blsys);
266    
267     /* -- store the initial values of all the stuff */
268     integrator_output_write(blsys);
269     integrator_output_write_obs(blsys);
270    
271     /* specify where the returned values should be stored */
272     yret = y0;
273     ypret = yp0;
274    
275     /* advance solution in time, return values as yret and derivatives as ypret */
276     blsys->currentstep=1;
277     for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
278     t = samplelist_get(blsys->samples, t_index);
279    
280     CONSOLE_DEBUG("SOLVING UP TO t = %f", t);
281    
282     flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
283    
284     /* pass the values of everything back to the compiler */
285     integrator_set_t(blsys, (double)tret);
286     integrator_set_y(blsys, NV_DATA_S(yret));
287     integrator_set_ydot(blsys, NV_DATA_S(ypret));
288    
289     if(flag<0){
290     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
291     break;
292     }
293    
294     /* -- do something so that blsys knows the values of tret, yret and ypret */
295    
296     /* -- store the current values of all the stuff */
297     integrator_output_write(blsys);
298     integrator_output_write_obs(blsys);
299     }
300    
301     /* -- close the IntegratorReporter */
302     integrator_output_close(blsys);
303    
304     if(flag < 0){
305     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
306     return 0;
307     }
308    
309     /* get optional outputs */
310    
311     /* free solution memory */
312     N_VDestroy_Serial(yret);
313     N_VDestroy_Serial(ypret);
314    
315     /* free solver memory */
316     IDAFree(ida_mem);
317    
318     /* all done */
319     return 1;
320     }
321    
322     /*--------------------------------------------------
323     RESIDUALS AND JACOBIAN
324     */
325     /**
326     Function to evaluate system residuals, in the form required for IDA.
327    
328     Given tt, yy and yp, we need to evaluate and return rr.
329    
330     @param tt current value of indep variable (time)
331     @param yy current values of dependent variable vector
332     @param yp current values of derivatives of dependent variables
333     @param rr the output residual vector (is we're returning data to)
334     @param res_data pointer to our stuff (blsys in this case).
335    
336     @return 0 on success, positive on recoverable error, and
337     negative on unrecoverable error.
338     */
339 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
340 johnpye 669 IntegratorSystem *blsys;
341     IntegratorIdaData *enginedata;
342     int i, calc_ok, is_error;
343     struct rel_relation** relptr;
344     double resid;
345     char *relname;
346    
347     blsys = (IntegratorSystem *)res_data;
348     enginedata = integrator_ida_enginedata(blsys);
349    
350 johnpye 903 /* fprintf(stderr,"\n\n"); */
351     /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
352 johnpye 669
353     /* pass the values of everything back to the compiler */
354     integrator_set_t(blsys, (double)tt);
355     integrator_set_y(blsys, NV_DATA_S(yy));
356     integrator_set_ydot(blsys, NV_DATA_S(yp));
357    
358     /* revaluate the system residuals using the new data */
359     is_error = 0;
360     relptr = enginedata->rellist;
361    
362 johnpye 903 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
363 johnpye 669 if(NV_LENGTH_S(rr)!=enginedata->nrels){
364     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
365     return -1; /* unrecoverable */
366     }
367    
368     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
369     if (setjmp(g_fpe_env)==0) {
370     for(i=0, relptr = enginedata->rellist;
371     i< enginedata->nrels && relptr != NULL;
372     ++i, ++relptr
373     ){
374     resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
375    
376     relname = rel_make_name(blsys->system, *relptr);
377 johnpye 903 /* if(calc_ok){
378 johnpye 669 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
379     }else{
380     CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
381 johnpye 903 }*/
382 johnpye 669 ASC_FREE(relname);
383    
384     NV_Ith_S(rr,i) = resid;
385     if(!calc_ok){
386     /* presumable some output already made? */
387     is_error = 1;
388     }
389     }
390     }else{
391     CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
392     }
393     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
394    
395     if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
396     return is_error;
397     }
398    
399 johnpye 903 /**
400 johnpye 912 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
401 johnpye 909
402     Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
403    
404     @param tt current value of the independent variable (time, t)
405     @param yy current value of the dependent variable vector, y(t).
406     @param yp current value of y'(t).
407     @param rr current value of the residual vector F(t, y, y').
408     @param v the vector by which the Jacobian must be multiplied to the right.
409     @param Jv the output vector computed
410     @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
411     @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
412     @param tmp1 @see tmp2
413     @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.
414     @return 0 on success
415 johnpye 903 */
416 johnpye 909 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
417     , N_Vector v, N_Vector Jv, realtype c_j
418     , void *jac_data, N_Vector tmp1, N_Vector tmp2
419     ){
420     IntegratorSystem *blsys;
421     IntegratorIdaData *enginedata;
422     int i, j, is_error=0;
423     struct rel_relation** relptr;
424     char *relname, *varname;
425     int status;
426     double Jv_i;
427 johnpye 903
428 johnpye 909 int *variables;
429     double *derivatives;
430     var_filter_t filter;
431     int count;
432    
433     CONSOLE_DEBUG("EVALUTING JACOBIAN...");
434    
435     blsys = (IntegratorSystem *)jac_data;
436     enginedata = integrator_ida_enginedata(blsys);
437    
438     /* pass the values of everything back to the compiler */
439     integrator_set_t(blsys, (double)tt);
440     integrator_set_y(blsys, NV_DATA_S(yy));
441     integrator_set_ydot(blsys, NV_DATA_S(yp));
442     /* no real use for residuals (rr) here, I don't think? */
443    
444 johnpye 912 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
445 johnpye 909 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
446     derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
447    
448     /* evaluate the derivatives... */
449     /* J = dG_dy = dF_dy + alpha * dF_dyp */
450    
451     filter.matchbits = VAR_SVAR;
452     filter.matchvalue = VAR_SVAR;
453    
454 johnpye 912 CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v));
455     for(i=0; i<NV_LENGTH_S(v); ++i){
456     varname = var_make_name(blsys->system, enginedata->varlist[i]);
457     if(varname==NULL)varname="UNKNOWN";
458     CONSOLE_DEBUG("v[%d] = '%s' = %f",i,varname,NV_Ith_S(v,i));
459     ASC_FREE(varname);
460     }
461    
462 johnpye 909 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
463     if (setjmp(g_fpe_env)==0) {
464     for(i=0, relptr = enginedata->rellist;
465     i< enginedata->nrels && relptr != NULL;
466     ++i, ++relptr
467     ){
468     /* get derivatives for this particular relation */
469     status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
470    
471     CONSOLE_DEBUG("Got derivatives against %d matching variables", count);
472    
473     relname = rel_make_name(blsys->system, *relptr);
474     if(!status){
475     CONSOLE_DEBUG("Derivatives for relation %d '%s' OK",i,relname);
476     }else{
477     CONSOLE_DEBUG("ERROR calculating derivatives for relation %d '%s'",i,relname);
478     break;
479     }
480     ASC_FREE(relname);
481    
482     Jv_i = 0;
483     for(j=0; j < count; ++j){
484 johnpye 910 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
485 johnpye 909 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
486     if(varname){
487     CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
488     ASC_FREE(varname);
489     }else{
490 johnpye 912 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
491 johnpye 909 }
492 johnpye 912
493 johnpye 909 Jv_i += derivatives[j] * NV_Ith_S(v,variables[j]);
494     }
495    
496     NV_Ith_S(Jv,i) = Jv_i;
497 johnpye 910 CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i);
498    
499 johnpye 909 if(status){
500     /* presumably some error_reporter will already have been made*/
501     is_error = 1;
502     }
503     }
504     }else{
505     CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
506     }
507     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
508    
509     if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
510 johnpye 912
511    
512    
513 johnpye 909 return is_error;
514     }
515    
516 johnpye 669 /*----------------------------------------------
517     ERROR REPORTING
518 johnpye 782 */
519 johnpye 669 /**
520     Error message reporter function to be passed to IDA. All error messages
521     will trigger a call to this function, so we should find everything
522     appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
523     panel (in the case of PyGTK).
524     */
525     void integrator_ida_error(int error_code
526     , const char *module, const char *function
527     , char *msg, void *eh_data
528     ){
529     IntegratorSystem *blsys;
530     error_severity_t sev;
531    
532     /* cast back the IntegratorSystem, just in case we need it */
533     blsys = (IntegratorSystem *)eh_data;
534    
535     /* severity depends on the sign of the error_code value */
536     if(error_code <= 0){
537     sev = ASC_PROG_ERR;
538     }else{
539     sev = ASC_PROG_WARNING;
540     }
541    
542     /* use our all-purpose error reporting to get stuff back to the GUI */
543     error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
544     }

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