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

Contents of /trunk/solvers/ida/ida.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: 32295 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 static IntegratorCreateFn integrator_ida_create;
134 static IntegratorParamsDefaultFn integrator_ida_params_default;
135 static IntegratorSolveFn integrator_ida_solve;
136 static IntegratorFreeFn integrator_ida_free;
137
138 /**
139 Everthing that the outside world needs to know about IDA
140 */
141 static const IntegratorInternals integrator_ida_internals = {
142 integrator_ida_create
143 ,integrator_ida_params_default
144 ,integrator_ida_analyse
145 ,integrator_ida_solve
146 ,integrator_ida_write_matrix
147 ,integrator_ida_debug
148 ,integrator_ida_free
149 ,INTEG_IDA
150 ,"IDA"
151 };
152
153 extern ASC_EXPORT int ida_register(void){
154 CONSOLE_DEBUG("Registering IDA...");
155 return integrator_register(&integrator_ida_internals);
156 }
157
158 /*-------------------------------------------------------------
159 FORWARD DECLS
160 */
161
162 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *integ, struct var_variable *var, const int *varindx);
163
164 /*static IntegratorVarVisitorFn integrator_dae_classify_var;
165 static void integrator_visit_system_vars(IntegratorSystem *integ,IntegratorVarVisitorFn *visitor);
166 static void integrator_dae_show_var(IntegratorSystem *integ, struct var_variable *var, const int *varindx); */
167
168 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
169
170 /*-------------------------------------------------------------
171 SETUP/TEARDOWN ROUTINES
172 */
173 static void integrator_ida_create(IntegratorSystem *integ){
174 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
175 IntegratorIdaData *enginedata;
176 enginedata = ASC_NEW(IntegratorIdaData);
177 CONSOLE_DEBUG("enginedata = %p",enginedata);
178 enginedata->rellist = NULL;
179 enginedata->safeeval = 0;
180 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED;
181 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;
182 enginedata->pfree = NULL;
183
184 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
185 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
186
187 integ->enginedata = (void *)enginedata;
188
189 integrator_ida_params_default(integ);
190 }
191
192 static void integrator_ida_free(void *enginedata){
193 #ifdef DESTROY_DEBUG
194 CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);
195 #endif
196 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
197 asc_assert(d);
198 if(d->pfree){
199 CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);
200 /* free the preconditioner data, whatever it happens to be */
201 (d->pfree)(enginedata);
202 }
203
204 ASC_FREE(d->rellist);
205
206 #ifdef DESTROY_DEBUG
207 CONSOLE_DEBUG("Now destroying the enginedata");
208 #endif
209 ASC_FREE(d);
210 #ifdef DESTROY_DEBUG
211 CONSOLE_DEBUG("enginedata freed");
212 #endif
213 }
214
215 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *integ){
216 IntegratorIdaData *d;
217 assert(integ!=NULL);
218 assert(integ->enginedata!=NULL);
219 assert(integ->engine==INTEG_IDA);
220 d = ((IntegratorIdaData *)(integ->enginedata));
221 return d;
222 }
223
224 /*-------------------------------------------------------------
225 PARAMETERS FOR IDA
226 */
227
228 enum ida_parameters{
229 IDA_PARAM_LINSOLVER
230 ,IDA_PARAM_MAXL
231 ,IDA_PARAM_MAXORD
232 ,IDA_PARAM_AUTODIFF
233 ,IDA_PARAM_CALCIC
234 ,IDA_PARAM_SAFEEVAL
235 ,IDA_PARAM_RTOL
236 ,IDA_PARAM_ATOL
237 ,IDA_PARAM_ATOLVECT
238 ,IDA_PARAM_GSMODIFIED
239 ,IDA_PARAM_MAXNCF
240 ,IDA_PARAM_PREC
241 ,IDA_PARAMS_SIZE
242 };
243
244 /**
245 Here the full set of parameters is defined, along with upper/lower bounds,
246 etc. The values are stuck into the integ->params structure.
247
248 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
249 list. Then add a slv_param_*(...) statement below to define the type, description and range
250 for the new parameter.
251
252 @return 0 on success
253 */
254 static int integrator_ida_params_default(IntegratorSystem *integ){
255 asc_assert(integ!=NULL);
256 asc_assert(integ->engine==INTEG_IDA);
257 slv_parameters_t *p;
258 p = &(integ->params);
259
260 slv_destroy_parms(p);
261
262 if(p->parms==NULL){
263 CONSOLE_DEBUG("params NULL");
264 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
265 if(p->parms==NULL)return -1;
266 p->dynamic_parms = 1;
267 }else{
268 CONSOLE_DEBUG("params not NULL");
269 }
270
271 /* reset the number of parameters to zero so that we can check it at the end */
272 p->num_parms = 0;
273
274 slv_param_bool(p,IDA_PARAM_AUTODIFF
275 ,(SlvParameterInitBool){{"autodiff"
276 ,"Use auto-diff?",1
277 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
278 }, TRUE}
279 );
280
281 slv_param_char(p,IDA_PARAM_CALCIC
282 ,(SlvParameterInitChar){{"calcic"
283 ,"Initial conditions calcuation",1
284 ,"Use specified values of ydot to solve for inital y (Y),"
285 " or use the the values of the differential variables (yd) to solve"
286 " for the pure algebraic variables (ya) along with the derivatives"
287 " of the differential variables (yddot) (YA_YDP), or else don't solve"
288 " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
289 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
290 );
291
292 slv_param_bool(p,IDA_PARAM_SAFEEVAL
293 ,(SlvParameterInitBool){{"safeeval"
294 ,"Use safe evaluation?",1
295 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
296 "throw SIGFPE errors which will then halt integration."
297 }, FALSE}
298 );
299
300
301 slv_param_bool(p,IDA_PARAM_ATOLVECT
302 ,(SlvParameterInitBool){{"atolvect"
303 ,"Use 'ode_atol' values as specified?",1
304 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
305 " in the integration. If FALSE, a scalar absolute tolerance value"
306 " is shared by all variables. See IDA manual, section 5.5.1"
307 }, TRUE }
308 );
309
310 slv_param_real(p,IDA_PARAM_ATOL
311 ,(SlvParameterInitReal){{"atol"
312 ,"Scalar absolute error tolerance",1
313 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
314 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
315 }, 1e-5, 0, 1e10 }
316 );
317
318 slv_param_real(p,IDA_PARAM_RTOL
319 ,(SlvParameterInitReal){{"rtol"
320 ,"Scalar relative error tolerance",1
321 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
322 " it's not possible to set per-variable relative tolerances as it is"
323 " with LSODE)."
324 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
325 }, 1e-4, 0, 1 }
326 );
327
328 slv_param_char(p,IDA_PARAM_LINSOLVER
329 ,(SlvParameterInitChar){{"linsolver"
330 ,"Linear solver",1
331 ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"
332 " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"
333 " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"
334 " or SPTFQMR (which still need preconditioners to be implemented"
335 " before they can be very useful."
336 }, "DENSE"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
337 );
338
339 slv_param_int(p,IDA_PARAM_MAXL
340 ,(SlvParameterInitInt){{"maxl"
341 ,"Maximum Krylov dimension",0
342 ,"The maximum dimension of Krylov space used by the linear solver"
343 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
344 " The default of 0 results in IDA using its internal default, which"
345 " is currently a value of 5."
346 }, 0, 0, 20 }
347 );
348
349 slv_param_int(p,IDA_PARAM_MAXORD
350 ,(SlvParameterInitInt){{"maxord"
351 ,"Maximum order of linear multistep method",0
352 ,"The maximum order of the linear multistep method with IDA. See"
353 " IDA manual p 38."
354 }, 5, 1, 5 }
355 );
356
357 slv_param_bool(p,IDA_PARAM_GSMODIFIED
358 ,(SlvParameterInitBool){{"gsmodified"
359 ,"Gram-Schmidt Orthogonalisation Scheme", 2
360 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
361 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
362 }, TRUE}
363 );
364
365 slv_param_int(p,IDA_PARAM_MAXNCF
366 ,(SlvParameterInitInt){{"maxncf"
367 ,"Max nonlinear solver convergence failures per step", 2
368 ,"Maximum number of allowable nonlinear solver convergence failures"
369 " on one step. See IDA manual section 5.5.6.1."
370 }, 10,0,1000 }
371 );
372
373 slv_param_char(p,IDA_PARAM_PREC
374 ,(SlvParameterInitChar){{"prec"
375 ,"Preconditioner",1
376 ,"See IDA manual, section section 5.6.8."
377 },"NONE"}, (char *[]){"NONE","DIAG",NULL}
378 );
379
380 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
381
382 CONSOLE_DEBUG("Created %d params", p->num_parms);
383
384 return 0;
385 }
386
387 /*-------------------------------------------------------------
388 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
389 */
390
391 typedef int IdaFlagFn(void *,int *);
392 typedef char *IdaFlagNameFn(int);
393
394 /* return 0 on success */
395 static int integrator_ida_solve(
396 IntegratorSystem *integ
397 , unsigned long start_index
398 , unsigned long finish_index
399 ){
400 void *ida_mem;
401 int flag, flag1, t_index;
402 realtype t0, reltol, abstol, t, tret, tout1;
403 N_Vector y0, yp0, abstolvect, ypret, yret, id;
404 IntegratorIdaData *enginedata;
405 char *linsolver;
406 int maxl;
407 IdaFlagFn *flagfn;
408 IdaFlagNameFn *flagnamefn;
409 const char *flagfntype;
410 char *pname = NULL;
411 struct rel_relation **rels;
412 int *rootsfound;
413 char havecrossed;
414
415 #ifdef SOLVE_DEBUG
416 char *varname, *relname;
417 #endif
418
419 int i,j,n_activerels,n_solverrels;
420 const IntegratorIdaPrec *prec = NULL;
421 int icopt; /* initial conditions strategy */
422
423 CONSOLE_DEBUG("STARTING IDA...");
424
425 enginedata = integrator_ida_enginedata(integ);
426
427 enginedata->safeeval = SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_SAFEEVAL);
428 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
429
430 /* store reference to list of relations (in enginedata) */
431 n_solverrels = slv_get_num_solvers_rels(integ->system);
432
433 n_activerels = slv_count_solvers_rels(integ->system, &integrator_ida_rel);
434
435 enginedata->bndlist = slv_get_solvers_bnd_list(integ->system);
436 enginedata->nbnds = slv_get_num_solvers_bnds(integ->system);
437
438 enginedata->rellist = ASC_NEW_ARRAY(struct rel_relation *, n_activerels);
439
440 rels = slv_get_solvers_rel_list(integ->system);
441
442 j=0;
443 for(i=0; i < n_solverrels; ++i){
444 if(rel_apply_filter(rels[i], &integrator_ida_rel)){
445 #ifdef SOLVE_DEBUG
446 relname = rel_make_name(integ->system, rels[i]);
447 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));
448 ASC_FREE(relname);
449 #endif
450 enginedata->rellist[j++]=rels[i];
451 }
452 }
453 asc_assert(j==n_activerels);
454
455 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
456 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
457
458 CONSOLE_DEBUG("Number of relations: %d",n_solverrels);
459 CONSOLE_DEBUG("Number of active relations: %d",n_activerels);
460 CONSOLE_DEBUG("Number of dependent vars: %d",integ->n_y);
461 CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);
462
463 enginedata->nrels = n_activerels;
464
465 if(enginedata->nrels != integ->n_y){
466 ERROR_REPORTER_HERE(ASC_USER_ERROR
467 ,"Integration problem is not square (%d active rels, %d vars)"
468 ,n_activerels, integ->n_y
469 );
470 return 1; /* failure */
471 }
472
473 #ifdef SOLVE_DEBUG
474 integrator_ida_debug(integ,stderr);
475 #endif
476
477 /* retrieve initial values from the system */
478
479 /** @TODO fix this, the starting time != first sample */
480 t0 = integrator_get_t(integ);
481 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
482
483 CONSOLE_DEBUG("RETRIEVING y0");
484
485 y0 = N_VNew_Serial(integ->n_y);
486 integrator_get_y(integ,NV_DATA_S(y0));
487
488 #ifdef SOLVE_DEBUG
489 CONSOLE_DEBUG("RETRIEVING yp0");
490 #endif
491
492 yp0 = N_VNew_Serial(integ->n_y);
493 integrator_get_ydot(integ,NV_DATA_S(yp0));
494
495 #ifdef SOLVE_DEBUG
496 N_VPrint_Serial(yp0);
497 CONSOLE_DEBUG("yp0 is at %p",&yp0);
498 #endif
499
500 /* create IDA object */
501 ida_mem = IDACreate();
502
503 /* relative error tolerance */
504 reltol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_RTOL);
505 CONSOLE_DEBUG("rtol = %8.2e",reltol);
506
507
508 /* allocate internal memory */
509 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
510 flag = IDAInit(ida_mem, &integrator_ida_fex, t0, y0 ,yp0);
511 #else
512 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
513 /* vector of absolute tolerances */
514 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
515 abstolvect = N_VNew_Serial(integ->n_y);
516 integrator_get_atol(integ,NV_DATA_S(abstolvect));
517
518 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
519
520 N_VDestroy_Serial(abstolvect);
521 }else{
522 /* scalar absolute tolerance (one value for all) */
523 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
524 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
525 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
526 }
527 #endif
528
529 if(flag==IDA_MEM_NULL){
530 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
531 return 2;
532 }else if(flag==IDA_MEM_FAIL){
533 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
534 return 3;
535 }else if(flag==IDA_ILL_INPUT){
536 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
537 return 4;
538 }/* else success */
539
540
541 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
542 CONSOLE_DEBUG("Assigning tolerances...");
543 /* assign tolerances */
544 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
545 CONSOLE_DEBUG("using vector of atol values");
546 abstolvect = N_VNew_Serial(integ->n_y);
547 integrator_get_atol(integ,NV_DATA_S(abstolvect));
548 IDASVtolerances(ida_mem, reltol, abstolvect);
549 N_VDestroy_Serial(abstolvect);
550 }else{
551 /* scalar tolerances */
552 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
553 CONSOLE_DEBUG("using scalar atol value = %8.2e",abstol);
554 IDASStolerances(ida_mem, reltol, abstol);
555 }
556 #endif
557
558 /* set optional inputs... */
559 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)integ);
560 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
561 IDASetUserData(ida_mem, (void *)integ);
562 #else
563 IDASetRdata(ida_mem, (void *)integ);
564 #endif
565 IDASetMaxStep(ida_mem, integrator_get_maxstep(integ));
566 IDASetInitStep(ida_mem, integrator_get_stepzero(integ));
567 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(integ));
568 if(integrator_get_minstep(integ)>0){
569 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
570 }
571
572 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
573 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
574
575 CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
576 IDASetMaxOrd(ida_mem,SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
577
578 /* there's no capability for setting *minimum* step size in IDA */
579
580
581 /* attach linear solver module, using the default value of maxl */
582 linsolver = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_LINSOLVER);
583 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
584 if(strcmp(linsolver,"ASCEND")==0){
585 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",integ->n_y);
586 IDAASCEND(ida_mem,integ->n_y);
587 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)integ);
588
589 flagfntype = "IDAASCEND";
590 flagfn = &IDAASCENDGetLastFlag;
591 flagnamefn = &IDAASCENDGetReturnFlagName;
592
593 }else if(strcmp(linsolver,"DENSE")==0){
594 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",integ->n_y);
595 flag = IDADense(ida_mem, integ->n_y);
596 switch(flag){
597 case IDADENSE_SUCCESS: break;
598 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 5;
599 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 5;
600 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 5;
601 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
602 }
603
604 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)){
605 CONSOLE_DEBUG("USING AUTODIFF");
606 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
607 flag = IDADlsSetDenseJacFn(ida_mem, &integrator_ida_djex);
608 #else
609 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)integ);
610 #endif
611 switch(flag){
612 case IDADENSE_SUCCESS: break;
613 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
614 }
615 }else{
616 CONSOLE_DEBUG("USING NUMERICAL DIFF");
617 }
618
619 flagfntype = "IDADENSE";
620 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
621 flagfn = &IDADlsGetLastFlag;
622 flagnamefn = &IDADlsGetReturnFlagName;
623 #else
624 flagfn = &IDADenseGetLastFlag;
625 flagnamefn = &IDADenseGetReturnFlagName;
626 #endif
627 }else{
628 /* remaining methods are all SPILS */
629 CONSOLE_DEBUG("IDA SPILS");
630
631 maxl = SLV_PARAM_INT(&(integ->params),IDA_PARAM_MAXL);
632 CONSOLE_DEBUG("maxl = %d",maxl);
633
634 /* what preconditioner? */
635 pname = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_PREC);
636 if(strcmp(pname,"NONE")==0){
637 prec = NULL;
638 }else if(strcmp(pname,"JACOBI")==0){
639 prec = &prec_jacobi;
640 }else{
641 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
642 return 7;
643 }
644
645 /* which SPILS linear solver? */
646 if(strcmp(linsolver,"SPGMR")==0){
647 CONSOLE_DEBUG("IDA SPGMR");
648 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
649 }else if(strcmp(linsolver,"SPBCG")==0){
650 CONSOLE_DEBUG("IDA SPBCG");
651 flag = IDASpbcg(ida_mem, maxl);
652 }else if(strcmp(linsolver,"SPTFQMR")==0){
653 CONSOLE_DEBUG("IDA SPTFQMR");
654 flag = IDASptfqmr(ida_mem,maxl);
655 }else{
656 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
657 return 8;
658 }
659
660 if(prec){
661 /* assign the preconditioner to the linear solver */
662 (prec->pcreate)(integ);
663 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
664 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve);
665 #else
666 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)integ);
667 #endif
668 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
669 }else{
670 CONSOLE_DEBUG("No preconditioner");
671 }
672
673 flagfntype = "IDASPILS";
674 flagfn = &IDASpilsGetLastFlag;
675 flagnamefn = &IDASpilsGetReturnFlagName;
676
677 if(flag==IDASPILS_MEM_NULL){
678 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
679 return 9;
680 }else if(flag==IDASPILS_MEM_FAIL){
681 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
682 return 9;
683 }/* else success */
684
685 /* assign the J*v function */
686 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)){
687 CONSOLE_DEBUG("USING AUTODIFF");
688 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
689 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex);
690 #else
691 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)integ);
692 #endif
693 if(flag==IDASPILS_MEM_NULL){
694 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
695 return 10;
696 }else if(flag==IDASPILS_LMEM_NULL){
697 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
698 return 10;
699 }/* else success */
700 }else{
701 CONSOLE_DEBUG("USING NUMERICAL DIFF");
702 }
703
704 if(strcmp(linsolver,"SPGMR")==0){
705 /* select Gram-Schmidt orthogonalisation */
706 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_GSMODIFIED)){
707 CONSOLE_DEBUG("USING MODIFIED GS");
708 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
709 if(flag!=IDASPILS_SUCCESS){
710 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
711 return 11;
712 }
713 }else{
714 CONSOLE_DEBUG("USING CLASSICAL GS");
715 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
716 if(flag!=IDASPILS_SUCCESS){
717 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
718 return 11;
719 }
720 }
721 }
722 }
723
724 /* set linear solver optional inputs...
725 ...nothing here at the moment...
726 */
727
728 /* calculate initial conditions */
729 icopt = 0;
730 if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"Y")==0){
731 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
732 icopt = IDA_Y_INIT;
733 asc_assert(icopt!=0);
734 }else if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
735 CONSOLE_DEBUG("Solving initial conditions using values of yd");
736 icopt = IDA_YA_YDP_INIT;
737 asc_assert(icopt!=0);
738 id = N_VNew_Serial(integ->n_y);
739 for(i=0; i < integ->n_y; ++i){
740 if(integ->ydot[i] == NULL){
741 NV_Ith_S(id,i) = 0.0;
742 #ifdef SOLVE_DEBUG
743 varname = var_make_name(integ->system,integ->y[i]);
744 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
745 ASC_FREE(varname);
746 #endif
747 }else{
748 #ifdef SOLVE_DEBUG
749 CONSOLE_DEBUG("y[%d] is differential",i);
750 #endif
751 NV_Ith_S(id,i) = 1.0;
752 }
753 }
754 IDASetId(ida_mem, id);
755 N_VDestroy_Serial(id);
756 }else if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"NONE")==0){
757 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
758 }else{
759 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
760 }
761
762 if(icopt){
763 integ->currentstep=0;
764 t_index=start_index + 1;
765 tout1 = samplelist_get(integ->samples, t_index);
766
767 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
768
769 #ifdef ASC_SIGNAL_TRAPS
770 /* catch SIGFPE if desired to */
771 if(enginedata->safeeval){
772 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
773 Asc_SignalHandlerPush(SIGFPE,SIG_DFL);
774 }else{
775 # ifdef FEX_DEBUG
776 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
777 # endif
778 Asc_SignalHandlerPushDefault(SIGFPE);
779 }
780 if (setjmp(g_fpe_env)==0) {
781 #endif
782
783 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=3
784 flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3 */
785 # else
786 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
787 # endif
788 /* check flags and output status */
789 switch(flag){
790 case IDA_SUCCESS:
791 CONSOLE_DEBUG("Initial conditions solved OK");
792 break;
793
794 case IDA_LSETUP_FAIL:
795 case IDA_LINIT_FAIL:
796 case IDA_LSOLVE_FAIL:
797 case IDA_NO_RECOVERY:
798 flag1 = -999;
799 flag = (flagfn)(ida_mem,&flag1);
800 if(flag){
801 ERROR_REPORTER_HERE(ASC_PROG_ERR
802 ,"Unable to retrieve error code from %s (err %d)"
803 ,flagfntype,flag
804 );
805 return 12;
806 }
807 ERROR_REPORTER_HERE(ASC_PROG_ERR
808 ,"%s returned flag '%s' (value = %d)"
809 ,flagfntype,(flagnamefn)(flag1),flag1
810 );
811 return 12;
812
813 default:
814 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
815 return 12;
816 }
817 #ifdef ASC_SIGNAL_TRAPS
818 }else{
819 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
820 return 13;
821 }
822
823 if(enginedata->safeeval){
824 Asc_SignalHandlerPop(SIGFPE,SIG_DFL);
825 }else{
826 CONSOLE_DEBUG("pop...");
827 Asc_SignalHandlerPopDefault(SIGFPE);
828 CONSOLE_DEBUG("...pop");
829 }
830 #endif
831 }/* icopt */
832
833 /* optionally, specify ROO-FINDING problem */
834 if(enginedata->nbnds){
835 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
836 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn);
837 #else
838 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)integ);
839 #endif
840 }
841
842 /* -- set up the IntegratorReporter */
843 integrator_output_init(integ);
844
845 /* -- store the initial values of all the stuff */
846 integrator_output_write(integ);
847 integrator_output_write_obs(integ);
848
849 /* specify where the returned values should be stored */
850 yret = y0;
851 ypret = yp0;
852
853 /* advance solution in time, return values as yret and derivatives as ypret */
854 integ->currentstep=1;
855 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++integ->currentstep){
856 t = samplelist_get(integ->samples, t_index);
857 t0 = integrator_get_t(integ);
858 asc_assert(t > t0);
859
860 #ifdef SOLVE_DEBUG
861 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
862 #endif
863
864 #ifdef ASC_SIGNAL_TRAPS
865 Asc_SignalHandlerPushDefault(SIGINT);
866 if(setjmp(g_int_env)==0) {
867 #endif
868 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
869 #ifdef ASC_SIGNAL_TRAPS
870 }else{
871 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Caught interrupt");
872 flag = -555;
873 }
874 Asc_SignalHandlerPopDefault(SIGINT);
875 #endif
876
877 /* this seems to work, so we can use it to avoid needing 'havecrossed'? */
878 if(flag == IDA_ROOT_RETURN){
879 CONSOLE_DEBUG("IDA reports root found!");
880 }
881
882 /* so we will check for roots found explicitly */
883 if(enginedata->nbnds){
884 rootsfound = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
885 havecrossed = 0;
886 if(IDA_SUCCESS == IDAGetRootInfo(ida_mem, rootsfound)){
887 for(i=0; i < enginedata->nbnds; ++i){
888 if(rootsfound[i]){
889 havecrossed = 1;
890 #ifdef SOLVE_DEBUG
891 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
892 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary '%s' crossed%s",relname,rootsfound[i]>0?" (increasing)":" (decreasing)");
893 ASC_FREE(relname);
894 #else
895 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary %d crossed!%s", i, ,rootsfound[i]>0?" (increasing)":" (decreasing)");
896 #endif
897 }
898 }
899
900 /* so, now we need to restart the integration. we will assume that
901 everything changes: number of variables, etc, etc, etc. */
902
903 if(havecrossed){
904 CONSOLE_DEBUG("Boundaries were crossed; need to reinitialise solver...");
905 /** try resetting the boundary states now? */
906 //IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)integ);
907
908 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
909 IDAReInit(ida_mem, tret, yret, ypret);
910 #elif SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR<3
911 /* TODO find out what version needs the following line... not sure. */
912 //IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret);
913
914 // FIXME this stuff has not been tested yet, and is very incomplete.
915
916 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
917 /* vector of absolute tolerances */
918 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
919 abstolvect = N_VNew_Serial(integ->n_y);
920 integrator_get_atol(integ,NV_DATA_S(abstolvect));
921 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SV, reltol, abstolvect);
922 N_VDestroy_Serial(abstolvect);
923 }else{
924 /* scalar absolute tolerance (one value for all) */
925 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
926 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
927 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SS, reltol, &abstol);
928 }
929 #else
930 /* allocate internal memory */
931 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
932 /* vector of absolute tolerances */
933 abstolvect = N_VNew_Serial(integ->n_y);
934 integrator_get_atol(integ,NV_DATA_S(abstolvect));
935 if(IDA_SUCCESS != IDAReInit(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect)){
936 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
937 }
938 N_VDestroy_Serial(abstolvect);
939 }else{
940 /* scalar absolute tolerance (one value for all) */
941 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
942 if(IDA_SUCCESS != IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol)){
943 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
944 }
945 }
946 #endif
947 }
948 }else{
949 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");
950 }
951 ASC_FREE(rootsfound);
952 }
953
954
955 /* pass the values of everything back to the compiler */
956 integrator_set_t(integ, (double)tret);
957 integrator_set_y(integ, NV_DATA_S(yret));
958 integrator_set_ydot(integ, NV_DATA_S(ypret));
959
960 if(flag<0){
961 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
962 break;
963 }
964
965 /* -- do something so that integ knows the values of tret, yret and ypret */
966
967 /* -- store the current values of all the stuff */
968 integrator_output_write(integ);
969 integrator_output_write_obs(integ);
970
971 }/* loop through next sample timestep */
972
973 /* -- close the IntegratorReporter */
974 integrator_output_close(integ);
975
976 /* get optional outputs */
977 #ifdef STATS_DEBUG
978 IntegratorIdaStats stats;
979 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
980 integrator_ida_write_stats(&stats);
981 }else{
982 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
983 }
984 #endif
985
986 /* free solution memory */
987 N_VDestroy_Serial(yret);
988 N_VDestroy_Serial(ypret);
989
990 /* free solver memory */
991 IDAFree(ida_mem);
992
993 if(flag < -500){
994 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", t);
995 return -flag;
996 }
997
998 if(flag < 0){
999 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1000 return 14;
1001 }
1002
1003 /* all done, success */
1004 return 0;
1005 }
1006
1007 /*----------------------------------------------
1008 STATS
1009 */
1010
1011 /**
1012 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1013 status in a struct instead of separately.
1014
1015 @return IDA_SUCCESS on success.
1016 */
1017 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1018
1019 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
1020
1021 int res;
1022
1023 /*
1024 There is an error in the documentation for this function in Sundials 2.2.
1025 According the the header file, the hinused stat is not provided.
1026 */
1027 res = IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups,
1028 &s->netfails, &s->qlast, &s->qcur, &s->hlast, &s->hcur,
1029 &s->tcur
1030 );
1031
1032 /* get the missing statistic */
1033 IDAGetActualInitStep(ida_mem, &s->hinused);
1034
1035 return res;
1036 #else
1037
1038 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1039 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1040 ,&s->hlast, &s->hcur, &s->tcur
1041 );
1042
1043 #endif
1044 }
1045

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