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

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