/[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 913 - (show annotations) (download) (as text)
Sat Oct 28 03:55:19 2006 UTC (13 years, 5 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 /* 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 # 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 /* 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 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
90 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 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
99
100 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 /* 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 enginedata->varlist = NULL;
120 enginedata->safeeval = 1;
121 blsys->enginedata = (void *)enginedata;
122 }
123
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 int size, flag, t_index;
153 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 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
166 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
167 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
168 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 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstol);
197 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 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
217
218 /* 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 /* 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
238 /* set linear solver optional inputs...
239
240 ...nothing here at the moment...
241
242 */
243
244 /* correct initial values, given derivatives */
245 blsys->currentstep=0;
246 t_index=start_index+1;
247 tout1 = samplelist_get(blsys->samples, t_index);
248
249 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
250 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
251 #else
252 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
253 #endif
254
255 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 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
340 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 /* fprintf(stderr,"\n\n"); */
351 /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
352
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 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
363 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 /* if(calc_ok){
378 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
379 }else{
380 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
381 }*/
382 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 /**
400 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
401
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 */
416 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
428 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 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
445 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 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 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 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
485 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 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
491 }
492
493 Jv_i += derivatives[j] * NV_Ith_S(v,variables[j]);
494 }
495
496 NV_Ith_S(Jv,i) = Jv_i;
497 CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i);
498
499 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
511
512
513 return is_error;
514 }
515
516 /*----------------------------------------------
517 ERROR REPORTING
518 */
519 /**
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