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

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