/[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 928 - (show annotations) (download) (as text)
Wed Nov 22 10:32:18 2006 UTC (15 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 18001 byte(s)
Commented out some stream redirection stuff for simplicity.
The CUnit test suite now works as expected (but without output suppression, for the moment).
Some more effort on IDA (ongoing).
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
21 as part of the GPL-licensed SUNDIALS solver package from LLNL.
22 @see http://www.llnl.gov/casc/sundials/
23 *//*
24 by John Pye, May 2006
25 */
26
27 /*
28 Be careful with the following. This file requires both the 'ida.h' from
29 SUNDIALS as well as the 'ida.h' from ASCEND. Make sure that we're getting
30 both of these; if you get problems check your build tool for the paths being
31 passed to the C preprocessor.
32 */
33
34 /* standard includes */
35 #include <signal.h>
36
37 /* ASCEND includes */
38 #include "ida.h"
39 #include <utilities/error.h>
40 #include <utilities/ascConfig.h>
41 #include <utilities/ascSignal.h>
42 #include <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 CONSOLE_DEBUG("RETRIEVING t0");
178
179 /* retrieve initial values from the system */
180 t0 = samplelist_get(blsys->samples,start_index);
181
182 CONSOLE_DEBUG("RETRIEVING y0");
183
184 y0 = N_VNew_Serial(size);
185 integrator_get_y(blsys,NV_DATA_S(y0));
186
187 CONSOLE_DEBUG("RETRIEVING yp0");
188
189 yp0 = N_VNew_Serial(size);
190 integrator_get_ydot(blsys,NV_DATA_S(yp0));
191
192 N_VPrint_Serial(yp0);
193 CONSOLE_DEBUG("yp0 is at %p",&yp0);
194
195 /* create IDA object */
196 ida_mem = IDACreate();
197
198 /* retrieve the absolute tolerance values for each variable */
199 abstol = N_VNew_Serial(size);
200 N_VConst(0.1,abstol); /** @TODO fill in the abstol values from the model */
201 reltol = 0.001;
202
203 /* allocate internal memory */
204 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstol);
205 if(flag==IDA_MEM_NULL){
206 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
207 return 0;
208 }else if(flag==IDA_MEM_FAIL){
209 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
210 return 0;
211 }else if(flag==IDA_ILL_INPUT){
212 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
213 return 0;
214 }/* else success */
215
216 /* set optional inputs... */
217 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
218 IDASetRdata(ida_mem, (void *)blsys);
219 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
220 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
221 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
222 /* there's no capability for setting *minimum* step size in IDA */
223
224 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
225
226 /* attach linear solver module, using the default value of maxl */
227 flag = IDASpgmr(ida_mem, 0);
228 if(flag==IDASPILS_MEM_NULL){
229 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
230 return 0;
231 }else if(flag==IDASPILS_MEM_FAIL){
232 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
233 return 0;
234 }/* else success */
235
236 /* assign the J*v function */
237 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
238 if(flag==IDASPILS_MEM_NULL){
239 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
240 return 0;
241 }else if(flag==IDASPILS_LMEM_NULL){
242 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
243 return 0;
244 }/* else success */
245
246 /* set linear solver optional inputs...
247
248 ...nothing here at the moment...
249
250 */
251
252 /* correct initial values, given derivatives */
253 blsys->currentstep=0;
254 t_index=start_index+1;
255 tout1 = samplelist_get(blsys->samples, t_index);
256
257 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
258 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
259 #else
260 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
261 #endif
262
263 if(flag!=IDA_SUCCESS){
264 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
265 return 0;
266 }/* else success */
267
268 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
269
270 /* optionally, specify ROO-FINDING PROBLEM */
271
272 /* -- set up the IntegratorReporter */
273 integrator_output_init(blsys);
274
275 /* -- store the initial values of all the stuff */
276 integrator_output_write(blsys);
277 integrator_output_write_obs(blsys);
278
279 /* specify where the returned values should be stored */
280 yret = y0;
281 ypret = yp0;
282
283 /* advance solution in time, return values as yret and derivatives as ypret */
284 blsys->currentstep=1;
285 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
286 t = samplelist_get(blsys->samples, t_index);
287
288 CONSOLE_DEBUG("SOLVING UP TO t = %f", t);
289
290 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
291
292 /* pass the values of everything back to the compiler */
293 integrator_set_t(blsys, (double)tret);
294 integrator_set_y(blsys, NV_DATA_S(yret));
295 integrator_set_ydot(blsys, NV_DATA_S(ypret));
296
297 if(flag<0){
298 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
299 break;
300 }
301
302 /* -- do something so that blsys knows the values of tret, yret and ypret */
303
304 /* -- store the current values of all the stuff */
305 integrator_output_write(blsys);
306 integrator_output_write_obs(blsys);
307 }
308
309 /* -- close the IntegratorReporter */
310 integrator_output_close(blsys);
311
312 if(flag < 0){
313 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
314 return 0;
315 }
316
317 /* get optional outputs */
318
319 /* free solution memory */
320 N_VDestroy_Serial(yret);
321 N_VDestroy_Serial(ypret);
322
323 /* free solver memory */
324 IDAFree(ida_mem);
325
326 /* all done */
327 return 1;
328 }
329
330 /*--------------------------------------------------
331 RESIDUALS AND JACOBIAN
332 */
333 /**
334 Function to evaluate system residuals, in the form required for IDA.
335
336 Given tt, yy and yp, we need to evaluate and return rr.
337
338 @param tt current value of indep variable (time)
339 @param yy current values of dependent variable vector
340 @param yp current values of derivatives of dependent variables
341 @param rr the output residual vector (is we're returning data to)
342 @param res_data pointer to our stuff (blsys in this case).
343
344 @return 0 on success, positive on recoverable error, and
345 negative on unrecoverable error.
346 */
347 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
348 IntegratorSystem *blsys;
349 IntegratorIdaData *enginedata;
350 int i, calc_ok, is_error;
351 struct rel_relation** relptr;
352 double resid;
353 char *relname;
354
355 blsys = (IntegratorSystem *)res_data;
356 enginedata = integrator_ida_enginedata(blsys);
357
358 /* fprintf(stderr,"\n\n"); */
359 /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
360
361 /* pass the values of everything back to the compiler */
362 integrator_set_t(blsys, (double)tt);
363 integrator_set_y(blsys, NV_DATA_S(yy));
364 integrator_set_ydot(blsys, NV_DATA_S(yp));
365
366 /* revaluate the system residuals using the new data */
367 is_error = 0;
368 relptr = enginedata->rellist;
369
370 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
371 if(NV_LENGTH_S(rr)!=enginedata->nrels){
372 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
373 return -1; /* unrecoverable */
374 }
375
376 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
377 if (setjmp(g_fpe_env)==0) {
378 for(i=0, relptr = enginedata->rellist;
379 i< enginedata->nrels && relptr != NULL;
380 ++i, ++relptr
381 ){
382 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
383
384 relname = rel_make_name(blsys->system, *relptr);
385 /* if(calc_ok){
386 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
387 }else{
388 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
389 }*/
390 ASC_FREE(relname);
391
392 NV_Ith_S(rr,i) = resid;
393 if(!calc_ok){
394 /* presumable some output already made? */
395 is_error = 1;
396 }
397 }
398 }else{
399 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
400 }
401 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
402
403 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
404 return is_error;
405 }
406
407 /**
408 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
409
410 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
411
412 @param tt current value of the independent variable (time, t)
413 @param yy current value of the dependent variable vector, y(t).
414 @param yp current value of y'(t).
415 @param rr current value of the residual vector F(t, y, y').
416 @param v the vector by which the Jacobian must be multiplied to the right.
417 @param Jv the output vector computed
418 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
419 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
420 @param tmp1 @see tmp2
421 @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.
422 @return 0 on success
423 */
424 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
425 , N_Vector v, N_Vector Jv, realtype c_j
426 , void *jac_data, N_Vector tmp1, N_Vector tmp2
427 ){
428 IntegratorSystem *blsys;
429 IntegratorIdaData *enginedata;
430 int i, j, is_error=0;
431 struct rel_relation** relptr;
432 char *relname, *varname;
433 int status;
434 double Jv_i;
435 int var_yindex;
436
437 int *variables;
438 double *derivatives;
439 var_filter_t filter;
440 int count;
441
442 fprintf(stderr,"\n--------------\n");
443 CONSOLE_DEBUG("EVALUTING JACOBIAN...");
444
445 blsys = (IntegratorSystem *)jac_data;
446 enginedata = integrator_ida_enginedata(blsys);
447
448 /* pass the values of everything back to the compiler */
449 integrator_set_t(blsys, (double)tt);
450 integrator_set_y(blsys, NV_DATA_S(yy));
451 integrator_set_ydot(blsys, NV_DATA_S(yp));
452 /* no real use for residuals (rr) here, I don't think? */
453
454 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
455 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
456 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
457
458 /* evaluate the derivatives... */
459 /* J = dG_dy = dF_dy + alpha * dF_dyp */
460
461 filter.matchbits = VAR_SVAR;
462 filter.matchvalue = VAR_SVAR;
463
464 CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v));
465 for(i=0; i<NV_LENGTH_S(v); ++i){
466 varname = var_make_name(blsys->system, enginedata->varlist[i]);
467 if(varname==NULL)varname="UNKNOWN";
468 CONSOLE_DEBUG("v[%d] = '%s' = %f",i,varname,NV_Ith_S(v,i));
469 ASC_FREE(varname);
470 }
471
472 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
473 if (setjmp(g_fpe_env)==0) {
474 for(i=0, relptr = enginedata->rellist;
475 i< enginedata->nrels && relptr != NULL;
476 ++i, ++relptr
477 ){
478 fprintf(stderr,"\n");
479 relname = rel_make_name(blsys->system, *relptr);
480 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
481 ASC_FREE(relname);
482
483 /* get derivatives for this particular relation */
484 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
485 CONSOLE_DEBUG("Got derivatives against %d matching variables", count);
486
487 for(j=0;j<count;++j){
488 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
489 CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);
490 ASC_FREE(varname);
491 }
492
493 if(!status){
494 CONSOLE_DEBUG("Derivatives for relation %d OK",i);
495 }else{
496 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
497 break;
498 }
499
500 /*
501 Now we have the derivatives wrt each alg/diff variable in the
502 present equation. variables[] points into the varlist. need
503 a mapping from the varlist to the y and ydot lists.
504 */
505
506 Jv_i = 0;
507 for(j=0; j < count; ++j){
508 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
509 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
510 if(varname){
511 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
512 ASC_FREE(varname);
513 }else{
514 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
515 }
516
517 /* this bit is still not right!!! */
518 var_yindex = blsys->y_id[variables[j]];
519
520 if(var_yindex >= 0){
521 CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));
522 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
523 }else{
524 var_yindex = -var_yindex-1;
525 CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));
526 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
527 }
528 }
529
530 NV_Ith_S(Jv,i) = Jv_i;
531 CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i);
532
533 if(status){
534 /* presumably some error_reporter will already have been made*/
535 is_error = 1;
536 }
537 }
538 }else{
539 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
540 }
541 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
542
543 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
544
545
546
547 return is_error;
548 }
549
550 /*----------------------------------------------
551 ERROR REPORTING
552 */
553 /**
554 Error message reporter function to be passed to IDA. All error messages
555 will trigger a call to this function, so we should find everything
556 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
557 panel (in the case of PyGTK).
558 */
559 void integrator_ida_error(int error_code
560 , const char *module, const char *function
561 , char *msg, void *eh_data
562 ){
563 IntegratorSystem *blsys;
564 error_severity_t sev;
565
566 /* cast back the IntegratorSystem, just in case we need it */
567 blsys = (IntegratorSystem *)eh_data;
568
569 /* severity depends on the sign of the error_code value */
570 if(error_code <= 0){
571 sev = ASC_PROG_ERR;
572 }else{
573 sev = ASC_PROG_WARNING;
574 }
575
576 /* use our all-purpose error reporting to get stuff back to the GUI */
577 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
578 }

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