/[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 782 - (show annotations) (download) (as text)
Tue Jul 25 02:09:12 2006 UTC (14 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 11651 byte(s)
Radu Serban sent me a preview of the new version of IDA which has new header file layout.
This patch updates for the new layout.
Also couple of minor fixes for gcc warnings in numlist and library.cpp.
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 int nrels;
78 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
79 } IntegratorIdaData;
80
81 /*-------------------------------------------------------------
82 FORWARD DECLS
83 */
84 /* residual function forward declaration */
85 int IDA_FEX(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
86
87 /* error handler forward declaration */
88 void integrator_ida_error(int error_code
89 , const char *module, const char *function
90 , char *msg, void *eh_data
91 );
92
93 /*-------------------------------------------------------------
94 SETUP/TEARDOWN ROUTINES
95 */
96 void integrator_ida_create(IntegratorSystem *blsys){
97 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
98 IntegratorIdaData *enginedata;
99 enginedata = ASC_NEW(IntegratorIdaData);
100 enginedata->rellist = NULL;
101 enginedata->safeeval = 1;
102 blsys->enginedata = (void *)enginedata;
103 }
104
105 void integrator_ida_free(void *enginedata){
106 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
107 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
108 /* note, we don't own the rellist, so don't need to free it */
109 ASC_FREE(d);
110 }
111
112 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
113 IntegratorIdaData *d;
114 assert(blsys!=NULL);
115 assert(blsys->enginedata!=NULL);
116 assert(blsys->engine==INTEG_IDA);
117 d = ((IntegratorIdaData *)(blsys->enginedata));
118 assert(d->safeeval = 1);
119 return d;
120 }
121
122 /*-------------------------------------------------------------
123 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
124 */
125
126 /* return 1 on success */
127 int integrator_ida_solve(
128 IntegratorSystem *blsys
129 , unsigned long start_index
130 , unsigned long finish_index
131 ){
132 void *ida_mem;
133 int size, flag, t_index;
134 realtype t0, reltol, t, tret, tout1;
135 N_Vector y0, yp0, abstol, ypret, yret;
136 IntegratorIdaData *enginedata;
137
138 CONSOLE_DEBUG("STARTING IDA...");
139
140 enginedata = integrator_ida_enginedata(blsys);
141 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
142
143 /* store reference to list of relations (in enginedata) */
144 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
145 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
146 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
147 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
148 size = blsys->n_y;
149
150 if(enginedata->nrels!=size){
151 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
152 return 0; /* failure */
153 }
154
155 /* retrieve initial values from the system */
156 t0 = samplelist_get(blsys->samples,start_index);
157
158 y0 = N_VNew_Serial(size);
159 integrator_get_y(blsys,NV_DATA_S(y0));
160
161 yp0 = N_VNew_Serial(size);
162 integrator_get_ydot(blsys,NV_DATA_S(yp0));
163
164 N_VPrint_Serial(yp0);
165 CONSOLE_DEBUG("yp0 is at %p",&yp0);
166
167 /* create IDA object */
168 ida_mem = IDACreate();
169
170 /* retrieve the absolute tolerance values for each variable */
171 abstol = N_VNew_Serial(size);
172 N_VConst(0.1,abstol); /** @TODO fill in the abstol values from the model */
173 reltol = 0.001;
174
175 /* allocate internal memory */
176 flag = IDAMalloc(ida_mem, IDA_FEX, t0, y0, yp0, IDA_SV, reltol, abstol);
177 if(flag==IDA_MEM_NULL){
178 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
179 return 0;
180 }else if(flag==IDA_MEM_FAIL){
181 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
182 return 0;
183 }else if(flag==IDA_ILL_INPUT){
184 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
185 return 0;
186 }/* else success */
187
188 /* set optional inputs... */
189 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
190 IDASetRdata(ida_mem, (void *)blsys);
191 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
192 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
193 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
194 /* there's no capability for setting *minimum* step size in IDA */
195
196 /* attach linear solver module, using the default value of maxl */
197 flag = IDASpgmr(ida_mem, 0);
198 if(flag==IDASPILS_MEM_NULL){
199 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
200 return 0;
201 }else if(flag==IDASPILS_MEM_FAIL){
202 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
203 return 0;
204 }/* else success */
205
206 /* set linear solver optional inputs... */
207
208
209 /* correct initial values, given derivatives */
210 blsys->currentstep=0;
211 t_index=start_index+1;
212 tout1 = samplelist_get(blsys->samples, t_index);
213 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
214 if(flag!=IDA_SUCCESS){
215 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
216 return 0;
217 }/* else success */
218
219 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
220
221 /* optionally, specify ROO-FINDING PROBLEM */
222
223 /* -- set up the IntegratorReporter */
224 integrator_output_init(blsys);
225
226 /* -- store the initial values of all the stuff */
227 integrator_output_write(blsys);
228 integrator_output_write_obs(blsys);
229
230 /* specify where the returned values should be stored */
231 yret = y0;
232 ypret = yp0;
233
234 /* advance solution in time, return values as yret and derivatives as ypret */
235 blsys->currentstep=1;
236 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
237 t = samplelist_get(blsys->samples, t_index);
238
239 CONSOLE_DEBUG("SOLVING UP TO t = %f", t);
240
241 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
242
243 /* pass the values of everything back to the compiler */
244 integrator_set_t(blsys, (double)tret);
245 integrator_set_y(blsys, NV_DATA_S(yret));
246 integrator_set_ydot(blsys, NV_DATA_S(ypret));
247
248 if(flag<0){
249 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
250 break;
251 }
252
253 /* -- do something so that blsys knows the values of tret, yret and ypret */
254
255 /* -- store the current values of all the stuff */
256 integrator_output_write(blsys);
257 integrator_output_write_obs(blsys);
258 }
259
260 /* -- close the IntegratorReporter */
261 integrator_output_close(blsys);
262
263 if(flag < 0){
264 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
265 return 0;
266 }
267
268 /* get optional outputs */
269
270 /* free solution memory */
271 N_VDestroy_Serial(yret);
272 N_VDestroy_Serial(ypret);
273
274 /* free solver memory */
275 IDAFree(ida_mem);
276
277 /* all done */
278 return 1;
279 }
280
281 /*--------------------------------------------------
282 RESIDUALS AND JACOBIAN
283 */
284 /**
285 Function to evaluate system residuals, in the form required for IDA.
286
287 Given tt, yy and yp, we need to evaluate and return rr.
288
289 @param tt current value of indep variable (time)
290 @param yy current values of dependent variable vector
291 @param yp current values of derivatives of dependent variables
292 @param rr the output residual vector (is we're returning data to)
293 @param res_data pointer to our stuff (blsys in this case).
294
295 @return 0 on success, positive on recoverable error, and
296 negative on unrecoverable error.
297 */
298 int IDA_FEX(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
299 IntegratorSystem *blsys;
300 IntegratorIdaData *enginedata;
301 int i, calc_ok, is_error;
302 struct rel_relation** relptr;
303 double resid;
304 char *relname;
305
306 blsys = (IntegratorSystem *)res_data;
307 enginedata = integrator_ida_enginedata(blsys);
308
309 fprintf(stderr,"\n\n");
310 CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS...");
311
312 /* pass the values of everything back to the compiler */
313 integrator_set_t(blsys, (double)tt);
314 integrator_set_y(blsys, NV_DATA_S(yy));
315 integrator_set_ydot(blsys, NV_DATA_S(yp));
316
317 /* revaluate the system residuals using the new data */
318 is_error = 0;
319 relptr = enginedata->rellist;
320
321 CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr));
322 if(NV_LENGTH_S(rr)!=enginedata->nrels){
323 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
324 return -1; /* unrecoverable */
325 }
326
327 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
328 if (setjmp(g_fpe_env)==0) {
329 for(i=0, relptr = enginedata->rellist;
330 i< enginedata->nrels && relptr != NULL;
331 ++i, ++relptr
332 ){
333 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
334
335 relname = rel_make_name(blsys->system, *relptr);
336 if(calc_ok){
337 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
338 }else{
339 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
340 }
341 ASC_FREE(relname);
342
343 NV_Ith_S(rr,i) = resid;
344 if(!calc_ok){
345 /* presumable some output already made? */
346 is_error = 1;
347 }
348 }
349 }else{
350 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
351 }
352 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
353
354 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
355 return is_error;
356 }
357
358 /*----------------------------------------------
359 ERROR REPORTING
360 */
361 /**
362 Error message reporter function to be passed to IDA. All error messages
363 will trigger a call to this function, so we should find everything
364 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
365 panel (in the case of PyGTK).
366 */
367 void integrator_ida_error(int error_code
368 , const char *module, const char *function
369 , char *msg, void *eh_data
370 ){
371 IntegratorSystem *blsys;
372 error_severity_t sev;
373
374 /* cast back the IntegratorSystem, just in case we need it */
375 blsys = (IntegratorSystem *)eh_data;
376
377 /* severity depends on the sign of the error_code value */
378 if(error_code <= 0){
379 sev = ASC_PROG_ERR;
380 }else{
381 sev = ASC_PROG_WARNING;
382 }
383
384 /* use our all-purpose error reporting to get stuff back to the GUI */
385 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
386 }

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