/[ascend]/trunk/solvers/ida/idaio.c
ViewVC logotype

Contents of /trunk/solvers/ida/idaio.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2370 - (show annotations) (download) (as text)
Mon Jan 31 11:33:22 2011 UTC (13 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 13362 byte(s)
Split IO routines into extra file.
1 /* ASCEND modelling environment
2 Copyright (C) 2006-2011 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
23 IDA provides the non-linear parts, as well as a number of pluggable linear
24 solvers: dense, banded and krylov types.
25
26 We also implement here an EXPERIMENTAL direct sparse linear solver for IDA
27 using the ASCEND linsolqr routines.
28
29 @see http://www.llnl.gov/casc/sundials/
30 *//*
31 by John Pye, May 2006
32 */
33
34 #define _GNU_SOURCE
35
36 #include <signal.h>
37 #include <setjmp.h>
38 #include <fenv.h>
39 #include <math.h>
40
41 /* SUNDIALS includes */
42 #ifdef ASC_WITH_IDA
43
44 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
45 # include <sundials/sundials_config.h>
46 # include <sundials/sundials_nvector.h>
47 # include <ida/ida_spgmr.h>
48 # include <ida.h>
49 # include <nvector_serial.h>
50 #else
51 # include <sundials/sundials_config.h>
52 # include <nvector/nvector_serial.h>
53 # include <ida/ida.h>
54 #endif
55
56 # include <sundials/sundials_dense.h>
57 # include <ida/ida_spgmr.h>
58 # include <ida/ida_spbcgs.h>
59 # include <ida/ida_sptfqmr.h>
60 # include <ida/ida_dense.h>
61
62 # ifndef IDA_SUCCESS
63 # error "Failed to include SUNDIALS IDA header file"
64 # endif
65 #else
66 # error "If you're building this file, you should have ASC_WITH_IDA"
67 #endif
68
69 #ifdef ASC_WITH_MMIO
70 # include <mmio.h>
71 #endif
72
73 #include <ascend/general/platform.h>
74 #include <ascend/utilities/error.h>
75 #include <ascend/utilities/ascSignal.h>
76 #include <ascend/general/panic.h>
77 #include <ascend/compiler/instance_enum.h>
78
79 #include <ascend/system/slv_client.h>
80 #include <ascend/system/relman.h>
81 #include <ascend/system/block.h>
82 #include <ascend/system/slv_stdcalls.h>
83 #include <ascend/system/jacobian.h>
84 #include <ascend/system/bndman.h>
85
86 #include <ascend/utilities/config.h>
87 #include <ascend/integrator/integrator.h>
88
89 #include "idalinear.h"
90 #include "idaanalyse.h"
91 #include "ida_impl.h"
92 #include "idaprec.h"
93 #include "idacalc.h"
94 #include "idaio.h"
95
96 /*
97 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
98 */
99 #ifndef SUNDIALS_VERSION_MINOR
100 # ifdef __GNUC__
101 # warning "GUESSING SUNDIALS VERSION 2.2"
102 # endif
103 # define SUNDIALS_VERSION_MINOR 2
104 #endif
105 #ifndef SUNDIALS_VERSION_MAJOR
106 # define SUNDIALS_VERSION_MAJOR 2
107 #endif
108
109 /* SUNDIALS 2.4.0 introduces new DlsMat in place of DenseMat */
110 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==4
111 # define IDA_MTX_T DlsMat
112 # define IDADENSE_SUCCESS IDADLS_SUCCESS
113 # define IDADENSE_MEM_NULL IDADLS_MEM_NULL
114 # define IDADENSE_ILL_INPUT IDADLS_ILL_INPUT
115 # define IDADENSE_MEM_FAIL IDADLS_MEM_FAIL
116 #else
117 # define IDA_MTX_T DenseMat
118 #endif
119
120 /* #define FEX_DEBUG */
121 #define JEX_DEBUG
122 /* #define DJEX_DEBUG */
123 #define SOLVE_DEBUG
124 #define STATS_DEBUG
125 #define PREC_DEBUG
126 /* #define ROOT_DEBUG */
127
128 /* #define DIFFINDEX_DEBUG */
129 /* #define ANALYSE_DEBUG */
130 /* #define DESTROY_DEBUG */
131 /* #define MATRIX_DEBUG */
132
133 /**
134 This routine just outputs the stats to the CONSOLE_DEBUG routine.
135
136 @TODO provide a GUI way of stats reporting from IDA.
137 */
138 void integrator_ida_write_stats(IntegratorIdaStats *stats){
139 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
140 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
141 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
142 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
143 SI(qlast); SI(qcur);
144 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
145 # undef SL
146 # undef SI
147 # undef SR
148 }
149
150 /*------------------------------------------------------------------------------
151 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
152 */
153
154 /**
155 Here we construct the local transfer matrix. It's a bit of a naive
156 approach; probably far more efficient ways can be worked out. But it will
157 hopefully be a useful way to examine stability of some spatial
158 discretisation schemes for PDAE systems.
159
160 http://ascendserver.cheme.cmu.edu/wiki/index.php/IDA#Stability
161 */
162 int integrator_ida_transfer_matrix(const IntegratorSystem *integ, struct SystemJacobianStruct *J){
163 int i=0, res;
164 enum submat{II_GA=0, II_GD, II_FA, II_FD, II_FDP, II_NUM};
165
166 const var_filter_t *matvf[II_NUM] = {
167 &system_vfilter_algeb
168 ,&system_vfilter_diff
169 ,&system_vfilter_algeb
170 ,&system_vfilter_diff
171 ,&system_vfilter_deriv
172 };
173
174 const rel_filter_t *matrf[II_NUM] = {
175 &system_rfilter_algeb
176 ,&system_rfilter_algeb
177 ,&system_rfilter_diff
178 ,&system_rfilter_diff
179 ,&system_rfilter_diff
180 };
181
182 struct SystemJacobianStruct D[II_NUM];
183
184 for(i=0;i<II_NUM;++i){
185 res = system_jacobian(integ->system, matrf[i], matvf[i], 1/*safe*/ ,&(D[i]));
186 }
187
188 /* compute inverses for matrices that need it */
189 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
190 return 1;
191 }
192
193 /**
194 Our task here is to write the matrices that IDA *should* be seeing. We
195 are actually making calls to relman_diff in order to do that, so we're
196 really going back to the variables in the actualy system and computing
197 row by row what the values are. This should mean just a single call to
198 each blackbox present in the system (if blackbox caching is working
199 correctly).
200 */
201 int integrator_ida_write_matrix(const IntegratorSystem *integ, FILE *f, const char *type){
202 /* IntegratorIdaData *enginedata; */
203 struct SystemJacobianStruct J = {NULL,NULL,NULL,0,0};
204 int status=1;
205 mtx_region_t R;
206
207 if(type==NULL)type = "dx'/dx";
208
209 if(0==strcmp(type,"dg/dz")){
210 CONSOLE_DEBUG("Calculating dg/dz...");
211 status = system_jacobian(integ->system
212 , &system_rfilter_algeb, &system_vfilter_algeb
213 , 1 /* safe */
214 , &J
215 );
216 }else if(0==strcmp(type,"dg/dx")){
217 CONSOLE_DEBUG("Calculating dg/dx...");
218 status = system_jacobian(integ->system
219 , &system_rfilter_algeb, &system_vfilter_diff
220 , 1 /* safe */
221 , &J
222 );
223 }else if(0==strcmp(type,"df/dx'")){
224 CONSOLE_DEBUG("Calculating df/dx'...");
225 status = system_jacobian(integ->system
226 , &system_rfilter_diff, &system_vfilter_deriv
227 , 1 /* safe */
228 , &J
229 );
230 }else if(0==strcmp(type,"df/dz")){
231 CONSOLE_DEBUG("Calculating df/dz...");
232 status = system_jacobian(integ->system
233 , &system_rfilter_diff, &system_vfilter_algeb
234 , 1 /* safe */
235 , &J
236 );
237 }else if(0==strcmp(type,"df/dx")){
238 CONSOLE_DEBUG("Calculating df/dx...");
239 status = system_jacobian(integ->system
240 , &system_rfilter_diff, &system_vfilter_diff
241 , 1 /* safe */
242 , &J
243 );
244 }else if(0==strcmp(type,"dF/dy")){
245 CONSOLE_DEBUG("Calculating dF/dy...");
246 status = system_jacobian(integ->system
247 , &system_rfilter_all, &system_vfilter_nonderiv
248 , 1 /* safe */
249 , &J
250 );
251 }else if(0==strcmp(type,"dF/dy'")){
252 CONSOLE_DEBUG("Calculating dF/dy'...");
253 status = system_jacobian(integ->system
254 , &system_rfilter_all, &system_vfilter_deriv
255 , 1 /* safe */
256 , &J
257 );
258 }else if(0==strcmp(type,"dx'/dx")){
259 /* system state transfer matrix dyd'/dyd */
260 status = integrator_ida_transfer_matrix(integ, &J);
261 }else{
262 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid matrix type '%s'",type);
263 return 1;
264 }
265
266 if(status){
267 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating matrix");
268 }else{
269 /* send the region explicitly, so that we handle non-square correctly */
270 R.row.low = 0; R.col.low = 0;
271 R.row.high = J.n_rels - 1; R.col.high = J.n_vars - 1;
272 /* note that we're not fussy about empty matrices here... */
273 mtx_write_region_mmio(f,J.M,&R);
274 }
275
276 if(J.vars)ASC_FREE(J.vars);
277 if(J.rels)ASC_FREE(J.rels);
278 if(J.M)mtx_destroy(J.M);
279
280 return status;
281 }
282
283 /**
284 This routine outputs matrix structure in a crude text format, for the sake
285 of debugging.
286 */
287 void integrator_ida_write_incidence(IntegratorSystem *integ){
288 int i, j;
289 struct rel_relation **relptr;
290 IntegratorIdaData *enginedata = integ->enginedata;
291 double *derivatives;
292 struct var_variable **variables;
293 int count, status;
294 char *relname;
295
296 if(enginedata->nrels > 100){
297 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
298 return;
299 }
300
301 variables = ASC_NEW_ARRAY(struct var_variable *, integ->n_y * 2);
302 derivatives = ASC_NEW_ARRAY(double, integ->n_y * 2);
303
304 CONSOLE_DEBUG("Outputting incidence information to console...");
305
306 for(i=0, relptr = enginedata->rellist;
307 i< enginedata->nrels && relptr != NULL;
308 ++i, ++relptr
309 ){
310 relname = rel_make_name(integ->system, *relptr);
311
312 /* get derivatives for this particular relation */
313 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
314 if(status){
315 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
316 ASC_FREE(relname);
317 break;
318 }
319
320 fprintf(stderr,"%3d:%-15s:",i,relname);
321 ASC_FREE(relname);
322
323 for(j=0; j<count; ++j){
324 if(var_deriv(variables[j])){
325 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(integ,variables[j]));
326 }else{
327 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
328 }
329 }
330 fprintf(stderr,"\n");
331 }
332 ASC_FREE(variables);
333 ASC_FREE(derivatives);
334 }
335
336 /* @return 0 on success */
337 int integrator_ida_debug(const IntegratorSystem *integ, FILE *fp){
338 char *varname, *relname;
339 struct var_variable **vlist, *var;
340 struct rel_relation **rlist, *rel;
341 long vlen, rlen;
342 long i;
343 long di;
344
345 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",integ->n_y);
346
347 /* if(integrator_sort_obs_vars(integ))return 10; */
348
349 if(integ->y && integ->ydot){
350 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
351 fprintf(fp,"index\t%-15s\tydot\n","y");
352 fprintf(fp,"-----\t%-15s\t-----\n","-----");
353 for(i=0;i<integ->n_y;++i){
354 varname = var_make_name(integ->system, integ->y[i]);
355 fprintf(fp,"%ld\t%-15s\t",i,varname);
356 if(integ->ydot[i]){
357 ASC_FREE(varname);
358 varname = var_make_name(integ->system, integ->ydot[i]);
359 fprintf(fp,"%s\n",varname);
360 ASC_FREE(varname);
361 }else{
362 fprintf(fp,".\n");
363 ASC_FREE(varname);
364 }
365 }
366 }else{
367 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
368 }
369
370 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
371 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
372 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
373
374
375 /* visit all the slv_system_t master var lists to collect vars */
376 /* find the vars mostly in this one */
377 vlist = slv_get_solvers_var_list(integ->system);
378 vlen = slv_get_num_solvers_vars(integ->system);
379 for(i=0;i<vlen;i++){
380 var = vlist[i];
381
382 varname = var_make_name(integ->system, var);
383 fprintf(fp,"%ld\t%-15s\t",i,varname);
384
385 if(var_fixed(var)){
386 // it's fixed, so not really a DAE var
387 fprintf(fp,"(fixed)\n");
388 }else if(!var_active(var)){
389 // inactive
390 fprintf(fp,"(inactive)\n");
391 }else if(!var_incident(var)){
392 // not incident
393 fprintf(fp,"(not incident)\n");
394 }else{
395 if(var_deriv(var)){
396 if(integ->y_id){
397 di = integrator_ida_diffindex1(integ,var);
398 if(di>=0){
399 ASC_FREE(varname);
400 varname = var_make_name(integ->system,vlist[di]);
401 fprintf(fp,".\tdiff(%ld='%s')\n",di,varname);
402 }else{
403 fprintf(fp,".\tdiff(???,err=%ld)\n",di);
404 }
405 }else{
406 fprintf(fp,".\tderiv... of??\n");
407 }
408 }else{
409 fprintf(fp,"%d\t.\n",var_sindex(var));
410 }
411 }
412 ASC_FREE(varname);
413 }
414
415 /* let's write out the relations too */
416 rlist = slv_get_solvers_rel_list(integ->system);
417 rlen = slv_get_num_solvers_rels(integ->system);
418
419 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
420 fprintf(fp,"index\tname\n");
421 fprintf(fp,"-----\t----\n");
422 for(i=0; i<rlen; ++i){
423 rel = rlist[i];
424 relname = rel_make_name(integ->system,rel);
425 fprintf(fp,"%ld\t%s\n",i,relname);
426 ASC_FREE(relname);
427 }
428
429 /* write out the derivative chains */
430 fprintf(fp,"\nDERIVATIVE CHAINS\n");
431 if(integrator_ida_analyse_debug(integ,stderr)){
432 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
433 return 340;
434 }
435 fprintf(fp,"\n");
436
437 /* and lets write block debug output */
438 system_block_debug(integ->system, fp);
439
440 return 0; /* success */
441 }
442
443 /*----------------------------------------------
444 ERROR REPORTING
445 */
446 /**
447 Error message reporter function to be passed to IDA. All error messages
448 will trigger a call to this function, so we should find everything
449 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
450 panel (in the case of PyGTK).
451 */
452 void integrator_ida_error(int error_code
453 , const char *module, const char *function
454 , char *msg, void *eh_data
455 ){
456 IntegratorSystem *integ;
457 error_severity_t sev;
458
459 /* cast back the IntegratorSystem, just in case we need it */
460 integ = (IntegratorSystem *)eh_data;
461
462 /* severity depends on the sign of the error_code value */
463 if(error_code <= 0){
464 sev = ASC_PROG_ERR;
465 }else{
466 sev = ASC_PROG_WARNING;
467 }
468
469 /* use our all-purpose error reporting to get stuff back to the GUI */
470 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
471 }
472

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