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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2369 - (show annotations) (download) (as text)
Mon Jan 31 09:00:44 2011 UTC (13 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 42594 byte(s)
Working to refactor ida.c, which is too big and scary.
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
95 /*
96 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
97 */
98 #ifndef SUNDIALS_VERSION_MINOR
99 # ifdef __GNUC__
100 # warning "GUESSING SUNDIALS VERSION 2.2"
101 # endif
102 # define SUNDIALS_VERSION_MINOR 2
103 #endif
104 #ifndef SUNDIALS_VERSION_MAJOR
105 # define SUNDIALS_VERSION_MAJOR 2
106 #endif
107
108 /* SUNDIALS 2.4.0 introduces new DlsMat in place of DenseMat */
109 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==4
110 # define IDA_MTX_T DlsMat
111 # define IDADENSE_SUCCESS IDADLS_SUCCESS
112 # define IDADENSE_MEM_NULL IDADLS_MEM_NULL
113 # define IDADENSE_ILL_INPUT IDADLS_ILL_INPUT
114 # define IDADENSE_MEM_FAIL IDADLS_MEM_FAIL
115 #else
116 # define IDA_MTX_T DenseMat
117 #endif
118
119 /* #define FEX_DEBUG */
120 #define JEX_DEBUG
121 /* #define DJEX_DEBUG */
122 #define SOLVE_DEBUG
123 #define STATS_DEBUG
124 #define PREC_DEBUG
125 /* #define ROOT_DEBUG */
126
127 /* #define DIFFINDEX_DEBUG */
128 /* #define ANALYSE_DEBUG */
129 /* #define DESTROY_DEBUG */
130 /* #define MATRIX_DEBUG */
131
132 static IntegratorCreateFn integrator_ida_create;
133 static IntegratorParamsDefaultFn integrator_ida_params_default;
134 static IntegratorSolveFn integrator_ida_solve;
135 static IntegratorFreeFn integrator_ida_free;
136 IntegratorDebugFn integrator_ida_debug;
137 static IntegratorWriteMatrixFn integrator_ida_write_matrix;
138
139 /**
140 Everthing that the outside world needs to know about IDA
141 */
142 static const IntegratorInternals integrator_ida_internals = {
143 integrator_ida_create
144 ,integrator_ida_params_default
145 ,integrator_ida_analyse
146 ,integrator_ida_solve
147 ,integrator_ida_write_matrix
148 ,integrator_ida_debug
149 ,integrator_ida_free
150 ,INTEG_IDA
151 ,"IDA"
152 };
153
154 extern ASC_EXPORT int ida_register(void){
155 CONSOLE_DEBUG("Registering IDA...");
156 return integrator_register(&integrator_ida_internals);
157 }
158
159 /*-------------------------------------------------------------
160 FORWARD DECLS
161 */
162
163 /* error handler forward declaration */
164 static void integrator_ida_error(int error_code
165 , const char *module, const char *function
166 , char *msg, void *eh_data
167 );
168
169 typedef struct IntegratorIdaStatsStruct{
170 long nsteps;
171 long nrevals;
172 long nlinsetups;
173 long netfails;
174 int qlast, qcur;
175 realtype hinused, hlast, hcur;
176 realtype tcur;
177 } IntegratorIdaStats;
178
179 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *integ, struct var_variable *var, const int *varindx);
180
181 /*static IntegratorVarVisitorFn integrator_dae_classify_var;
182 static void integrator_visit_system_vars(IntegratorSystem *integ,IntegratorVarVisitorFn *visitor);
183 static void integrator_dae_show_var(IntegratorSystem *integ, struct var_variable *var, const int *varindx); */
184
185 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
186 static void integrator_ida_write_stats(IntegratorIdaStats *stats);
187
188
189 /*-------------------------------------------------------------
190 SETUP/TEARDOWN ROUTINES
191 */
192 static void integrator_ida_create(IntegratorSystem *integ){
193 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
194 IntegratorIdaData *enginedata;
195 enginedata = ASC_NEW(IntegratorIdaData);
196 CONSOLE_DEBUG("enginedata = %p",enginedata);
197 enginedata->rellist = NULL;
198 enginedata->safeeval = 0;
199 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED;
200 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;
201 enginedata->pfree = NULL;
202
203 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
204 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
205
206 integ->enginedata = (void *)enginedata;
207
208 integrator_ida_params_default(integ);
209 }
210
211 static void integrator_ida_free(void *enginedata){
212 #ifdef DESTROY_DEBUG
213 CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);
214 #endif
215 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
216 asc_assert(d);
217 if(d->pfree){
218 CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);
219 /* free the preconditioner data, whatever it happens to be */
220 (d->pfree)(enginedata);
221 }
222
223 ASC_FREE(d->rellist);
224
225 #ifdef DESTROY_DEBUG
226 CONSOLE_DEBUG("Now destroying the enginedata");
227 #endif
228 ASC_FREE(d);
229 #ifdef DESTROY_DEBUG
230 CONSOLE_DEBUG("enginedata freed");
231 #endif
232 }
233
234 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *integ){
235 IntegratorIdaData *d;
236 assert(integ!=NULL);
237 assert(integ->enginedata!=NULL);
238 assert(integ->engine==INTEG_IDA);
239 d = ((IntegratorIdaData *)(integ->enginedata));
240 return d;
241 }
242
243 /*-------------------------------------------------------------
244 PARAMETERS FOR IDA
245 */
246
247 enum ida_parameters{
248 IDA_PARAM_LINSOLVER
249 ,IDA_PARAM_MAXL
250 ,IDA_PARAM_MAXORD
251 ,IDA_PARAM_AUTODIFF
252 ,IDA_PARAM_CALCIC
253 ,IDA_PARAM_SAFEEVAL
254 ,IDA_PARAM_RTOL
255 ,IDA_PARAM_ATOL
256 ,IDA_PARAM_ATOLVECT
257 ,IDA_PARAM_GSMODIFIED
258 ,IDA_PARAM_MAXNCF
259 ,IDA_PARAM_PREC
260 ,IDA_PARAMS_SIZE
261 };
262
263 /**
264 Here the full set of parameters is defined, along with upper/lower bounds,
265 etc. The values are stuck into the integ->params structure.
266
267 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
268 list. Then add a slv_param_*(...) statement below to define the type, description and range
269 for the new parameter.
270
271 @return 0 on success
272 */
273 static int integrator_ida_params_default(IntegratorSystem *integ){
274 asc_assert(integ!=NULL);
275 asc_assert(integ->engine==INTEG_IDA);
276 slv_parameters_t *p;
277 p = &(integ->params);
278
279 slv_destroy_parms(p);
280
281 if(p->parms==NULL){
282 CONSOLE_DEBUG("params NULL");
283 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
284 if(p->parms==NULL)return -1;
285 p->dynamic_parms = 1;
286 }else{
287 CONSOLE_DEBUG("params not NULL");
288 }
289
290 /* reset the number of parameters to zero so that we can check it at the end */
291 p->num_parms = 0;
292
293 slv_param_bool(p,IDA_PARAM_AUTODIFF
294 ,(SlvParameterInitBool){{"autodiff"
295 ,"Use auto-diff?",1
296 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
297 }, TRUE}
298 );
299
300 slv_param_char(p,IDA_PARAM_CALCIC
301 ,(SlvParameterInitChar){{"calcic"
302 ,"Initial conditions calcuation",1
303 ,"Use specified values of ydot to solve for inital y (Y),"
304 " or use the the values of the differential variables (yd) to solve"
305 " for the pure algebraic variables (ya) along with the derivatives"
306 " of the differential variables (yddot) (YA_YDP), or else don't solve"
307 " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
308 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
309 );
310
311 slv_param_bool(p,IDA_PARAM_SAFEEVAL
312 ,(SlvParameterInitBool){{"safeeval"
313 ,"Use safe evaluation?",1
314 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
315 "throw SIGFPE errors which will then halt integration."
316 }, FALSE}
317 );
318
319
320 slv_param_bool(p,IDA_PARAM_ATOLVECT
321 ,(SlvParameterInitBool){{"atolvect"
322 ,"Use 'ode_atol' values as specified?",1
323 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
324 " in the integration. If FALSE, a scalar absolute tolerance value"
325 " is shared by all variables. See IDA manual, section 5.5.1"
326 }, TRUE }
327 );
328
329 slv_param_real(p,IDA_PARAM_ATOL
330 ,(SlvParameterInitReal){{"atol"
331 ,"Scalar absolute error tolerance",1
332 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
333 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
334 }, 1e-5, 0, 1e10 }
335 );
336
337 slv_param_real(p,IDA_PARAM_RTOL
338 ,(SlvParameterInitReal){{"rtol"
339 ,"Scalar relative error tolerance",1
340 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
341 " it's not possible to set per-variable relative tolerances as it is"
342 " with LSODE)."
343 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
344 }, 1e-4, 0, 1 }
345 );
346
347 slv_param_char(p,IDA_PARAM_LINSOLVER
348 ,(SlvParameterInitChar){{"linsolver"
349 ,"Linear solver",1
350 ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"
351 " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"
352 " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"
353 " or SPTFQMR (which still need preconditioners to be implemented"
354 " before they can be very useful."
355 }, "DENSE"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
356 );
357
358 slv_param_int(p,IDA_PARAM_MAXL
359 ,(SlvParameterInitInt){{"maxl"
360 ,"Maximum Krylov dimension",0
361 ,"The maximum dimension of Krylov space used by the linear solver"
362 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
363 " The default of 0 results in IDA using its internal default, which"
364 " is currently a value of 5."
365 }, 0, 0, 20 }
366 );
367
368 slv_param_int(p,IDA_PARAM_MAXORD
369 ,(SlvParameterInitInt){{"maxord"
370 ,"Maximum order of linear multistep method",0
371 ,"The maximum order of the linear multistep method with IDA. See"
372 " IDA manual p 38."
373 }, 5, 1, 5 }
374 );
375
376 slv_param_bool(p,IDA_PARAM_GSMODIFIED
377 ,(SlvParameterInitBool){{"gsmodified"
378 ,"Gram-Schmidt Orthogonalisation Scheme", 2
379 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
380 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
381 }, TRUE}
382 );
383
384 slv_param_int(p,IDA_PARAM_MAXNCF
385 ,(SlvParameterInitInt){{"maxncf"
386 ,"Max nonlinear solver convergence failures per step", 2
387 ,"Maximum number of allowable nonlinear solver convergence failures"
388 " on one step. See IDA manual section 5.5.6.1."
389 }, 10,0,1000 }
390 );
391
392 slv_param_char(p,IDA_PARAM_PREC
393 ,(SlvParameterInitChar){{"prec"
394 ,"Preconditioner",1
395 ,"See IDA manual, section section 5.6.8."
396 },"NONE"}, (char *[]){"NONE","DIAG",NULL}
397 );
398
399 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
400
401 CONSOLE_DEBUG("Created %d params", p->num_parms);
402
403 return 0;
404 }
405
406 /*-------------------------------------------------------------
407 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
408 */
409
410 typedef int IdaFlagFn(void *,int *);
411 typedef char *IdaFlagNameFn(int);
412
413 /* return 0 on success */
414 static int integrator_ida_solve(
415 IntegratorSystem *integ
416 , unsigned long start_index
417 , unsigned long finish_index
418 ){
419 void *ida_mem;
420 int flag, flag1, t_index;
421 realtype t0, reltol, abstol, t, tret, tout1;
422 N_Vector y0, yp0, abstolvect, ypret, yret, id;
423 IntegratorIdaData *enginedata;
424 char *linsolver;
425 int maxl;
426 IdaFlagFn *flagfn;
427 IdaFlagNameFn *flagnamefn;
428 const char *flagfntype;
429 char *pname = NULL;
430 struct rel_relation **rels;
431 int *rootsfound;
432 char havecrossed;
433
434 #ifdef SOLVE_DEBUG
435 char *varname, *relname;
436 #endif
437
438 int i,j,n_activerels,n_solverrels;
439 const IntegratorIdaPrec *prec = NULL;
440 int icopt; /* initial conditions strategy */
441
442 CONSOLE_DEBUG("STARTING IDA...");
443
444 enginedata = integrator_ida_enginedata(integ);
445
446 enginedata->safeeval = SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_SAFEEVAL);
447 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
448
449 /* store reference to list of relations (in enginedata) */
450 n_solverrels = slv_get_num_solvers_rels(integ->system);
451
452 n_activerels = slv_count_solvers_rels(integ->system, &integrator_ida_rel);
453
454 enginedata->bndlist = slv_get_solvers_bnd_list(integ->system);
455 enginedata->nbnds = slv_get_num_solvers_bnds(integ->system);
456
457 enginedata->rellist = ASC_NEW_ARRAY(struct rel_relation *, n_activerels);
458
459 rels = slv_get_solvers_rel_list(integ->system);
460
461 j=0;
462 for(i=0; i < n_solverrels; ++i){
463 if(rel_apply_filter(rels[i], &integrator_ida_rel)){
464 #ifdef SOLVE_DEBUG
465 relname = rel_make_name(integ->system, rels[i]);
466 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));
467 ASC_FREE(relname);
468 #endif
469 enginedata->rellist[j++]=rels[i];
470 }
471 }
472 asc_assert(j==n_activerels);
473
474 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
475 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
476
477 CONSOLE_DEBUG("Number of relations: %d",n_solverrels);
478 CONSOLE_DEBUG("Number of active relations: %d",n_activerels);
479 CONSOLE_DEBUG("Number of dependent vars: %d",integ->n_y);
480 CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);
481
482 enginedata->nrels = n_activerels;
483
484 if(enginedata->nrels != integ->n_y){
485 ERROR_REPORTER_HERE(ASC_USER_ERROR
486 ,"Integration problem is not square (%d active rels, %d vars)"
487 ,n_activerels, integ->n_y
488 );
489 return 1; /* failure */
490 }
491
492 #ifdef SOLVE_DEBUG
493 integrator_ida_debug(integ,stderr);
494 #endif
495
496 /* retrieve initial values from the system */
497
498 /** @TODO fix this, the starting time != first sample */
499 t0 = integrator_get_t(integ);
500 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
501
502 CONSOLE_DEBUG("RETRIEVING y0");
503
504 y0 = N_VNew_Serial(integ->n_y);
505 integrator_get_y(integ,NV_DATA_S(y0));
506
507 #ifdef SOLVE_DEBUG
508 CONSOLE_DEBUG("RETRIEVING yp0");
509 #endif
510
511 yp0 = N_VNew_Serial(integ->n_y);
512 integrator_get_ydot(integ,NV_DATA_S(yp0));
513
514 #ifdef SOLVE_DEBUG
515 N_VPrint_Serial(yp0);
516 CONSOLE_DEBUG("yp0 is at %p",&yp0);
517 #endif
518
519 /* create IDA object */
520 ida_mem = IDACreate();
521
522 /* relative error tolerance */
523 reltol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_RTOL);
524 CONSOLE_DEBUG("rtol = %8.2e",reltol);
525
526
527 /* allocate internal memory */
528 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
529 flag = IDAInit(ida_mem, &integrator_ida_fex, t0, y0 ,yp0);
530 #else
531 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
532 /* vector of absolute tolerances */
533 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
534 abstolvect = N_VNew_Serial(integ->n_y);
535 integrator_get_atol(integ,NV_DATA_S(abstolvect));
536
537 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
538
539 N_VDestroy_Serial(abstolvect);
540 }else{
541 /* scalar absolute tolerance (one value for all) */
542 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
543 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
544 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
545 }
546 #endif
547
548 if(flag==IDA_MEM_NULL){
549 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
550 return 2;
551 }else if(flag==IDA_MEM_FAIL){
552 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
553 return 3;
554 }else if(flag==IDA_ILL_INPUT){
555 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
556 return 4;
557 }/* else success */
558
559
560 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
561 CONSOLE_DEBUG("Assigning tolerances...");
562 /* assign tolerances */
563 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
564 CONSOLE_DEBUG("using vector of atol values");
565 abstolvect = N_VNew_Serial(integ->n_y);
566 integrator_get_atol(integ,NV_DATA_S(abstolvect));
567 IDASVtolerances(ida_mem, reltol, abstolvect);
568 N_VDestroy_Serial(abstolvect);
569 }else{
570 /* scalar tolerances */
571 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
572 CONSOLE_DEBUG("using scalar atol value = %8.2e",abstol);
573 IDASStolerances(ida_mem, reltol, abstol);
574 }
575 #endif
576
577 /* set optional inputs... */
578 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)integ);
579 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
580 IDASetUserData(ida_mem, (void *)integ);
581 #else
582 IDASetRdata(ida_mem, (void *)integ);
583 #endif
584 IDASetMaxStep(ida_mem, integrator_get_maxstep(integ));
585 IDASetInitStep(ida_mem, integrator_get_stepzero(integ));
586 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(integ));
587 if(integrator_get_minstep(integ)>0){
588 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
589 }
590
591 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
592 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
593
594 CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
595 IDASetMaxOrd(ida_mem,SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
596
597 /* there's no capability for setting *minimum* step size in IDA */
598
599
600 /* attach linear solver module, using the default value of maxl */
601 linsolver = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_LINSOLVER);
602 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
603 if(strcmp(linsolver,"ASCEND")==0){
604 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",integ->n_y);
605 IDAASCEND(ida_mem,integ->n_y);
606 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)integ);
607
608 flagfntype = "IDAASCEND";
609 flagfn = &IDAASCENDGetLastFlag;
610 flagnamefn = &IDAASCENDGetReturnFlagName;
611
612 }else if(strcmp(linsolver,"DENSE")==0){
613 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",integ->n_y);
614 flag = IDADense(ida_mem, integ->n_y);
615 switch(flag){
616 case IDADENSE_SUCCESS: break;
617 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 5;
618 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 5;
619 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 5;
620 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
621 }
622
623 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)){
624 CONSOLE_DEBUG("USING AUTODIFF");
625 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
626 flag = IDADlsSetDenseJacFn(ida_mem, &integrator_ida_djex);
627 #else
628 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)integ);
629 #endif
630 switch(flag){
631 case IDADENSE_SUCCESS: break;
632 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
633 }
634 }else{
635 CONSOLE_DEBUG("USING NUMERICAL DIFF");
636 }
637
638 flagfntype = "IDADENSE";
639 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
640 flagfn = &IDADlsGetLastFlag;
641 flagnamefn = &IDADlsGetReturnFlagName;
642 #else
643 flagfn = &IDADenseGetLastFlag;
644 flagnamefn = &IDADenseGetReturnFlagName;
645 #endif
646 }else{
647 /* remaining methods are all SPILS */
648 CONSOLE_DEBUG("IDA SPILS");
649
650 maxl = SLV_PARAM_INT(&(integ->params),IDA_PARAM_MAXL);
651 CONSOLE_DEBUG("maxl = %d",maxl);
652
653 /* what preconditioner? */
654 pname = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_PREC);
655 if(strcmp(pname,"NONE")==0){
656 prec = NULL;
657 }else if(strcmp(pname,"JACOBI")==0){
658 prec = &prec_jacobi;
659 }else{
660 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
661 return 7;
662 }
663
664 /* which SPILS linear solver? */
665 if(strcmp(linsolver,"SPGMR")==0){
666 CONSOLE_DEBUG("IDA SPGMR");
667 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
668 }else if(strcmp(linsolver,"SPBCG")==0){
669 CONSOLE_DEBUG("IDA SPBCG");
670 flag = IDASpbcg(ida_mem, maxl);
671 }else if(strcmp(linsolver,"SPTFQMR")==0){
672 CONSOLE_DEBUG("IDA SPTFQMR");
673 flag = IDASptfqmr(ida_mem,maxl);
674 }else{
675 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
676 return 8;
677 }
678
679 if(prec){
680 /* assign the preconditioner to the linear solver */
681 (prec->pcreate)(integ);
682 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
683 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve);
684 #else
685 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)integ);
686 #endif
687 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
688 }else{
689 CONSOLE_DEBUG("No preconditioner");
690 }
691
692 flagfntype = "IDASPILS";
693 flagfn = &IDASpilsGetLastFlag;
694 flagnamefn = &IDASpilsGetReturnFlagName;
695
696 if(flag==IDASPILS_MEM_NULL){
697 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
698 return 9;
699 }else if(flag==IDASPILS_MEM_FAIL){
700 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
701 return 9;
702 }/* else success */
703
704 /* assign the J*v function */
705 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)){
706 CONSOLE_DEBUG("USING AUTODIFF");
707 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
708 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex);
709 #else
710 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)integ);
711 #endif
712 if(flag==IDASPILS_MEM_NULL){
713 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
714 return 10;
715 }else if(flag==IDASPILS_LMEM_NULL){
716 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
717 return 10;
718 }/* else success */
719 }else{
720 CONSOLE_DEBUG("USING NUMERICAL DIFF");
721 }
722
723 if(strcmp(linsolver,"SPGMR")==0){
724 /* select Gram-Schmidt orthogonalisation */
725 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_GSMODIFIED)){
726 CONSOLE_DEBUG("USING MODIFIED GS");
727 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
728 if(flag!=IDASPILS_SUCCESS){
729 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
730 return 11;
731 }
732 }else{
733 CONSOLE_DEBUG("USING CLASSICAL GS");
734 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
735 if(flag!=IDASPILS_SUCCESS){
736 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
737 return 11;
738 }
739 }
740 }
741 }
742
743 /* set linear solver optional inputs...
744 ...nothing here at the moment...
745 */
746
747 /* calculate initial conditions */
748 icopt = 0;
749 if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"Y")==0){
750 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
751 icopt = IDA_Y_INIT;
752 asc_assert(icopt!=0);
753 }else if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
754 CONSOLE_DEBUG("Solving initial conditions using values of yd");
755 icopt = IDA_YA_YDP_INIT;
756 asc_assert(icopt!=0);
757 id = N_VNew_Serial(integ->n_y);
758 for(i=0; i < integ->n_y; ++i){
759 if(integ->ydot[i] == NULL){
760 NV_Ith_S(id,i) = 0.0;
761 #ifdef SOLVE_DEBUG
762 varname = var_make_name(integ->system,integ->y[i]);
763 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
764 ASC_FREE(varname);
765 #endif
766 }else{
767 #ifdef SOLVE_DEBUG
768 CONSOLE_DEBUG("y[%d] is differential",i);
769 #endif
770 NV_Ith_S(id,i) = 1.0;
771 }
772 }
773 IDASetId(ida_mem, id);
774 N_VDestroy_Serial(id);
775 }else if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"NONE")==0){
776 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
777 }else{
778 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
779 }
780
781 if(icopt){
782 integ->currentstep=0;
783 t_index=start_index + 1;
784 tout1 = samplelist_get(integ->samples, t_index);
785
786 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
787
788 #ifdef ASC_SIGNAL_TRAPS
789 /* catch SIGFPE if desired to */
790 if(enginedata->safeeval){
791 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
792 Asc_SignalHandlerPush(SIGFPE,SIG_DFL);
793 }else{
794 # ifdef FEX_DEBUG
795 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
796 # endif
797 Asc_SignalHandlerPushDefault(SIGFPE);
798 }
799 if (setjmp(g_fpe_env)==0) {
800 #endif
801
802 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=3
803 flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3 */
804 # else
805 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
806 # endif
807 /* check flags and output status */
808 switch(flag){
809 case IDA_SUCCESS:
810 CONSOLE_DEBUG("Initial conditions solved OK");
811 break;
812
813 case IDA_LSETUP_FAIL:
814 case IDA_LINIT_FAIL:
815 case IDA_LSOLVE_FAIL:
816 case IDA_NO_RECOVERY:
817 flag1 = -999;
818 flag = (flagfn)(ida_mem,&flag1);
819 if(flag){
820 ERROR_REPORTER_HERE(ASC_PROG_ERR
821 ,"Unable to retrieve error code from %s (err %d)"
822 ,flagfntype,flag
823 );
824 return 12;
825 }
826 ERROR_REPORTER_HERE(ASC_PROG_ERR
827 ,"%s returned flag '%s' (value = %d)"
828 ,flagfntype,(flagnamefn)(flag1),flag1
829 );
830 return 12;
831
832 default:
833 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
834 return 12;
835 }
836 #ifdef ASC_SIGNAL_TRAPS
837 }else{
838 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
839 return 13;
840 }
841
842 if(enginedata->safeeval){
843 Asc_SignalHandlerPop(SIGFPE,SIG_DFL);
844 }else{
845 CONSOLE_DEBUG("pop...");
846 Asc_SignalHandlerPopDefault(SIGFPE);
847 CONSOLE_DEBUG("...pop");
848 }
849 #endif
850 }/* icopt */
851
852 /* optionally, specify ROO-FINDING problem */
853 if(enginedata->nbnds){
854 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
855 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn);
856 #else
857 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)integ);
858 #endif
859 }
860
861 /* -- set up the IntegratorReporter */
862 integrator_output_init(integ);
863
864 /* -- store the initial values of all the stuff */
865 integrator_output_write(integ);
866 integrator_output_write_obs(integ);
867
868 /* specify where the returned values should be stored */
869 yret = y0;
870 ypret = yp0;
871
872 /* advance solution in time, return values as yret and derivatives as ypret */
873 integ->currentstep=1;
874 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++integ->currentstep){
875 t = samplelist_get(integ->samples, t_index);
876 t0 = integrator_get_t(integ);
877 asc_assert(t > t0);
878
879 #ifdef SOLVE_DEBUG
880 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
881 #endif
882
883 #ifdef ASC_SIGNAL_TRAPS
884 Asc_SignalHandlerPushDefault(SIGINT);
885 if(setjmp(g_int_env)==0) {
886 #endif
887 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
888 #ifdef ASC_SIGNAL_TRAPS
889 }else{
890 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Caught interrupt");
891 flag = -555;
892 }
893 Asc_SignalHandlerPopDefault(SIGINT);
894 #endif
895
896 /* this seems to work, so we can use it to avoid needing 'havecrossed'? */
897 if(flag == IDA_ROOT_RETURN){
898 CONSOLE_DEBUG("IDA reports root found!");
899 }
900
901 /* so we will check for roots found explicitly */
902 if(enginedata->nbnds){
903 rootsfound = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
904 havecrossed = 0;
905 if(IDA_SUCCESS == IDAGetRootInfo(ida_mem, rootsfound)){
906 for(i=0; i < enginedata->nbnds; ++i){
907 if(rootsfound[i]){
908 havecrossed = 1;
909 #ifdef SOLVE_DEBUG
910 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
911 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary '%s' crossed%s",relname,rootsfound[i]>0?" (increasing)":" (decreasing)");
912 ASC_FREE(relname);
913 #else
914 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary %d crossed!%s", i, ,rootsfound[i]>0?" (increasing)":" (decreasing)");
915 #endif
916 }
917 }
918
919 /* so, now we need to restart the integration. we will assume that
920 everything changes: number of variables, etc, etc, etc. */
921
922 if(havecrossed){
923 CONSOLE_DEBUG("Boundaries were crossed; need to reinitialise solver...");
924 /** try resetting the boundary states now? */
925 //IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)integ);
926
927 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
928 IDAReInit(ida_mem, tret, yret, ypret);
929 #elif SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR<3
930 /* TODO find out what version needs the following line... not sure. */
931 //IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret);
932
933 // FIXME this stuff has not been tested yet, and is very incomplete.
934
935 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
936 /* vector of absolute tolerances */
937 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
938 abstolvect = N_VNew_Serial(integ->n_y);
939 integrator_get_atol(integ,NV_DATA_S(abstolvect));
940 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SV, reltol, abstolvect);
941 N_VDestroy_Serial(abstolvect);
942 }else{
943 /* scalar absolute tolerance (one value for all) */
944 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
945 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
946 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SS, reltol, &abstol);
947 }
948 #else
949 /* allocate internal memory */
950 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
951 /* vector of absolute tolerances */
952 abstolvect = N_VNew_Serial(integ->n_y);
953 integrator_get_atol(integ,NV_DATA_S(abstolvect));
954 if(IDA_SUCCESS != IDAReInit(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect)){
955 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
956 }
957 N_VDestroy_Serial(abstolvect);
958 }else{
959 /* scalar absolute tolerance (one value for all) */
960 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
961 if(IDA_SUCCESS != IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol)){
962 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
963 }
964 }
965 #endif
966 }
967 }else{
968 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");
969 }
970 ASC_FREE(rootsfound);
971 }
972
973
974 /* pass the values of everything back to the compiler */
975 integrator_set_t(integ, (double)tret);
976 integrator_set_y(integ, NV_DATA_S(yret));
977 integrator_set_ydot(integ, NV_DATA_S(ypret));
978
979 if(flag<0){
980 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
981 break;
982 }
983
984 /* -- do something so that integ knows the values of tret, yret and ypret */
985
986 /* -- store the current values of all the stuff */
987 integrator_output_write(integ);
988 integrator_output_write_obs(integ);
989
990 }/* loop through next sample timestep */
991
992 /* -- close the IntegratorReporter */
993 integrator_output_close(integ);
994
995 /* get optional outputs */
996 #ifdef STATS_DEBUG
997 IntegratorIdaStats stats;
998 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
999 integrator_ida_write_stats(&stats);
1000 }else{
1001 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
1002 }
1003 #endif
1004
1005 /* free solution memory */
1006 N_VDestroy_Serial(yret);
1007 N_VDestroy_Serial(ypret);
1008
1009 /* free solver memory */
1010 IDAFree(ida_mem);
1011
1012 if(flag < -500){
1013 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", t);
1014 return -flag;
1015 }
1016
1017 if(flag < 0){
1018 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1019 return 14;
1020 }
1021
1022 /* all done, success */
1023 return 0;
1024 }
1025
1026 /*----------------------------------------------
1027 STATS
1028 */
1029
1030 /**
1031 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1032 status in a struct instead of separately.
1033
1034 @return IDA_SUCCESS on success.
1035 */
1036 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1037
1038 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
1039
1040 int res;
1041
1042 /*
1043 There is an error in the documentation for this function in Sundials 2.2.
1044 According the the header file, the hinused stat is not provided.
1045 */
1046 res = IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups,
1047 &s->netfails, &s->qlast, &s->qcur, &s->hlast, &s->hcur,
1048 &s->tcur
1049 );
1050
1051 /* get the missing statistic */
1052 IDAGetActualInitStep(ida_mem, &s->hinused);
1053
1054 return res;
1055 #else
1056
1057 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1058 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1059 ,&s->hlast, &s->hcur, &s->tcur
1060 );
1061
1062 #endif
1063 }
1064
1065 /**
1066 This routine just outputs the stats to the CONSOLE_DEBUG routine.
1067
1068 @TODO provide a GUI way of stats reporting from IDA.
1069 */
1070 static void integrator_ida_write_stats(IntegratorIdaStats *stats){
1071 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
1072 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
1073 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
1074 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
1075 SI(qlast); SI(qcur);
1076 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
1077 # undef SL
1078 # undef SI
1079 # undef SR
1080 }
1081
1082 /*------------------------------------------------------------------------------
1083 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
1084 */
1085
1086 /**
1087 Here we construct the local transfer matrix. It's a bit of a naive
1088 approach; probably far more efficient ways can be worked out. But it will
1089 hopefully be a useful way to examine stability of some spatial
1090 discretisation schemes for PDAE systems.
1091
1092 http://ascendserver.cheme.cmu.edu/wiki/index.php/IDA#Stability
1093 */
1094 static int integrator_ida_transfer_matrix(const IntegratorSystem *integ, struct SystemJacobianStruct *J){
1095 int i=0, res;
1096 enum submat{II_GA=0, II_GD, II_FA, II_FD, II_FDP, II_NUM};
1097
1098 const var_filter_t *matvf[II_NUM] = {
1099 &system_vfilter_algeb
1100 ,&system_vfilter_diff
1101 ,&system_vfilter_algeb
1102 ,&system_vfilter_diff
1103 ,&system_vfilter_deriv
1104 };
1105
1106 const rel_filter_t *matrf[II_NUM] = {
1107 &system_rfilter_algeb
1108 ,&system_rfilter_algeb
1109 ,&system_rfilter_diff
1110 ,&system_rfilter_diff
1111 ,&system_rfilter_diff
1112 };
1113
1114 struct SystemJacobianStruct D[II_NUM];
1115
1116 for(i=0;i<II_NUM;++i){
1117 res = system_jacobian(integ->system, matrf[i], matvf[i], 1/*safe*/ ,&(D[i]));
1118 }
1119
1120 /* compute inverses for matrices that need it */
1121 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1122 return 1;
1123 }
1124
1125
1126 /**
1127 Our task here is to write the matrices that IDA *should* be seeing. We
1128 are actually making calls to relman_diff in order to do that, so we're
1129 really going back to the variables in the actualy system and computing
1130 row by row what the values are. This should mean just a single call to
1131 each blackbox present in the system (if blackbox caching is working
1132 correctly).
1133 */
1134 static int integrator_ida_write_matrix(const IntegratorSystem *integ, FILE *f, const char *type){
1135 /* IntegratorIdaData *enginedata; */
1136 struct SystemJacobianStruct J = {NULL,NULL,NULL,0,0};
1137 int status=1;
1138 mtx_region_t R;
1139
1140 if(type==NULL)type = "dx'/dx";
1141
1142 if(0==strcmp(type,"dg/dz")){
1143 CONSOLE_DEBUG("Calculating dg/dz...");
1144 status = system_jacobian(integ->system
1145 , &system_rfilter_algeb, &system_vfilter_algeb
1146 , 1 /* safe */
1147 , &J
1148 );
1149 }else if(0==strcmp(type,"dg/dx")){
1150 CONSOLE_DEBUG("Calculating dg/dx...");
1151 status = system_jacobian(integ->system
1152 , &system_rfilter_algeb, &system_vfilter_diff
1153 , 1 /* safe */
1154 , &J
1155 );
1156 }else if(0==strcmp(type,"df/dx'")){
1157 CONSOLE_DEBUG("Calculating df/dx'...");
1158 status = system_jacobian(integ->system
1159 , &system_rfilter_diff, &system_vfilter_deriv
1160 , 1 /* safe */
1161 , &J
1162 );
1163 }else if(0==strcmp(type,"df/dz")){
1164 CONSOLE_DEBUG("Calculating df/dz...");
1165 status = system_jacobian(integ->system
1166 , &system_rfilter_diff, &system_vfilter_algeb
1167 , 1 /* safe */
1168 , &J
1169 );
1170 }else if(0==strcmp(type,"df/dx")){
1171 CONSOLE_DEBUG("Calculating df/dx...");
1172 status = system_jacobian(integ->system
1173 , &system_rfilter_diff, &system_vfilter_diff
1174 , 1 /* safe */
1175 , &J
1176 );
1177 }else if(0==strcmp(type,"dF/dy")){
1178 CONSOLE_DEBUG("Calculating dF/dy...");
1179 status = system_jacobian(integ->system
1180 , &system_rfilter_all, &system_vfilter_nonderiv
1181 , 1 /* safe */
1182 , &J
1183 );
1184 }else if(0==strcmp(type,"dF/dy'")){
1185 CONSOLE_DEBUG("Calculating dF/dy'...");
1186 status = system_jacobian(integ->system
1187 , &system_rfilter_all, &system_vfilter_deriv
1188 , 1 /* safe */
1189 , &J
1190 );
1191 }else if(0==strcmp(type,"dx'/dx")){
1192 /* system state transfer matrix dyd'/dyd */
1193 status = integrator_ida_transfer_matrix(integ, &J);
1194 }else{
1195 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid matrix type '%s'",type);
1196 return 1;
1197 }
1198
1199 if(status){
1200 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating matrix");
1201 }else{
1202 /* send the region explicitly, so that we handle non-square correctly */
1203 R.row.low = 0; R.col.low = 0;
1204 R.row.high = J.n_rels - 1; R.col.high = J.n_vars - 1;
1205 /* note that we're not fussy about empty matrices here... */
1206 mtx_write_region_mmio(f,J.M,&R);
1207 }
1208
1209 if(J.vars)ASC_FREE(J.vars);
1210 if(J.rels)ASC_FREE(J.rels);
1211 if(J.M)mtx_destroy(J.M);
1212
1213 return status;
1214 }
1215
1216 /**
1217 This routine outputs matrix structure in a crude text format, for the sake
1218 of debugging.
1219 */
1220 void integrator_ida_write_incidence(IntegratorSystem *integ){
1221 int i, j;
1222 struct rel_relation **relptr;
1223 IntegratorIdaData *enginedata = integ->enginedata;
1224 double *derivatives;
1225 struct var_variable **variables;
1226 int count, status;
1227 char *relname;
1228
1229 if(enginedata->nrels > 100){
1230 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
1231 return;
1232 }
1233
1234 variables = ASC_NEW_ARRAY(struct var_variable *, integ->n_y * 2);
1235 derivatives = ASC_NEW_ARRAY(double, integ->n_y * 2);
1236
1237 CONSOLE_DEBUG("Outputting incidence information to console...");
1238
1239 for(i=0, relptr = enginedata->rellist;
1240 i< enginedata->nrels && relptr != NULL;
1241 ++i, ++relptr
1242 ){
1243 relname = rel_make_name(integ->system, *relptr);
1244
1245 /* get derivatives for this particular relation */
1246 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1247 if(status){
1248 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1249 ASC_FREE(relname);
1250 break;
1251 }
1252
1253 fprintf(stderr,"%3d:%-15s:",i,relname);
1254 ASC_FREE(relname);
1255
1256 for(j=0; j<count; ++j){
1257 if(var_deriv(variables[j])){
1258 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(integ,variables[j]));
1259 }else{
1260 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
1261 }
1262 }
1263 fprintf(stderr,"\n");
1264 }
1265 ASC_FREE(variables);
1266 ASC_FREE(derivatives);
1267 }
1268
1269 /* @return 0 on success */
1270 int integrator_ida_debug(const IntegratorSystem *integ, FILE *fp){
1271 char *varname, *relname;
1272 struct var_variable **vlist, *var;
1273 struct rel_relation **rlist, *rel;
1274 long vlen, rlen;
1275 long i;
1276 long di;
1277
1278 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",integ->n_y);
1279
1280 /* if(integrator_sort_obs_vars(integ))return 10; */
1281
1282 if(integ->y && integ->ydot){
1283 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
1284 fprintf(fp,"index\t%-15s\tydot\n","y");
1285 fprintf(fp,"-----\t%-15s\t-----\n","-----");
1286 for(i=0;i<integ->n_y;++i){
1287 varname = var_make_name(integ->system, integ->y[i]);
1288 fprintf(fp,"%ld\t%-15s\t",i,varname);
1289 if(integ->ydot[i]){
1290 ASC_FREE(varname);
1291 varname = var_make_name(integ->system, integ->ydot[i]);
1292 fprintf(fp,"%s\n",varname);
1293 ASC_FREE(varname);
1294 }else{
1295 fprintf(fp,".\n");
1296 ASC_FREE(varname);
1297 }
1298 }
1299 }else{
1300 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
1301 }
1302
1303 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
1304 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
1305 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
1306
1307
1308 /* visit all the slv_system_t master var lists to collect vars */
1309 /* find the vars mostly in this one */
1310 vlist = slv_get_solvers_var_list(integ->system);
1311 vlen = slv_get_num_solvers_vars(integ->system);
1312 for(i=0;i<vlen;i++){
1313 var = vlist[i];
1314
1315 varname = var_make_name(integ->system, var);
1316 fprintf(fp,"%ld\t%-15s\t",i,varname);
1317
1318 if(var_fixed(var)){
1319 // it's fixed, so not really a DAE var
1320 fprintf(fp,"(fixed)\n");
1321 }else if(!var_active(var)){
1322 // inactive
1323 fprintf(fp,"(inactive)\n");
1324 }else if(!var_incident(var)){
1325 // not incident
1326 fprintf(fp,"(not incident)\n");
1327 }else{
1328 if(var_deriv(var)){
1329 if(integ->y_id){
1330 di = integrator_ida_diffindex1(integ,var);
1331 if(di>=0){
1332 ASC_FREE(varname);
1333 varname = var_make_name(integ->system,vlist[di]);
1334 fprintf(fp,".\tdiff(%ld='%s')\n",di,varname);
1335 }else{
1336 fprintf(fp,".\tdiff(???,err=%ld)\n",di);
1337 }
1338 }else{
1339 fprintf(fp,".\tderiv... of??\n");
1340 }
1341 }else{
1342 fprintf(fp,"%d\t.\n",var_sindex(var));
1343 }
1344 }
1345 ASC_FREE(varname);
1346 }
1347
1348 /* let's write out the relations too */
1349 rlist = slv_get_solvers_rel_list(integ->system);
1350 rlen = slv_get_num_solvers_rels(integ->system);
1351
1352 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
1353 fprintf(fp,"index\tname\n");
1354 fprintf(fp,"-----\t----\n");
1355 for(i=0; i<rlen; ++i){
1356 rel = rlist[i];
1357 relname = rel_make_name(integ->system,rel);
1358 fprintf(fp,"%ld\t%s\n",i,relname);
1359 ASC_FREE(relname);
1360 }
1361
1362 /* write out the derivative chains */
1363 fprintf(fp,"\nDERIVATIVE CHAINS\n");
1364 if(integrator_ida_analyse_debug(integ,stderr)){
1365 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
1366 return 340;
1367 }
1368 fprintf(fp,"\n");
1369
1370 /* and lets write block debug output */
1371 system_block_debug(integ->system, fp);
1372
1373 return 0; /* success */
1374 }
1375
1376 /*----------------------------------------------
1377 ERROR REPORTING
1378 */
1379 /**
1380 Error message reporter function to be passed to IDA. All error messages
1381 will trigger a call to this function, so we should find everything
1382 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1383 panel (in the case of PyGTK).
1384 */
1385 static void integrator_ida_error(int error_code
1386 , const char *module, const char *function
1387 , char *msg, void *eh_data
1388 ){
1389 IntegratorSystem *integ;
1390 error_severity_t sev;
1391
1392 /* cast back the IntegratorSystem, just in case we need it */
1393 integ = (IntegratorSystem *)eh_data;
1394
1395 /* severity depends on the sign of the error_code value */
1396 if(error_code <= 0){
1397 sev = ASC_PROG_ERR;
1398 }else{
1399 sev = ASC_PROG_WARNING;
1400 }
1401
1402 /* use our all-purpose error reporting to get stuff back to the GUI */
1403 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1404 }

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