/[ascend]/trunk/base/generic/solver/ida.c
ViewVC logotype

Contents of /trunk/base/generic/solver/ida.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 996 - (show annotations) (download) (as text)
Sat Dec 23 06:53:04 2006 UTC (15 years, 6 months ago) by johnpye
File MIME type: text/x-csrc
File size: 47175 byte(s)
Added SIGFPE handling to IDACalcIC.
Added ctrl-F9 keystroke in Integrator dialog (PyGTK)
Fixed bun with 'Jac->N' not available under SUNDIALS 2.3.0.
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
21 as part of the GPL-licensed SUNDIALS solver package from LLNL.
22
23 IDA provides dense, banded and sparse solvers. At present this module only
24 implements access to the dense and sparse solvers, using the *serial*
25 vector methods.
26
27 @see http://www.llnl.gov/casc/sundials/
28 *//*
29 by John Pye, May 2006
30 */
31
32 /* standard includes */
33 #include <signal.h>
34
35 /* ASCEND includes */
36 #include "ida.h"
37 #include <utilities/error.h>
38 #include <utilities/ascConfig.h>
39 #include <utilities/ascSignal.h>
40 #include <utilities/ascPanic.h>
41 #include <compiler/instance_enum.h>
42 #include "var.h"
43 #include "rel.h"
44 #include "discrete.h"
45 #include "conditional.h"
46 #include "logrel.h"
47 #include "bnd.h"
48 #include "linsol.h"
49 #include "linsolqr.h"
50 #include "slv_common.h"
51 #include "slv_client.h"
52 #include "relman.h"
53
54 /* SUNDIALS includes */
55 #ifdef ASC_WITH_IDA
56 # include <sundials/sundials_config.h>
57 # include <sundials/sundials_dense.h>
58 # include <ida/ida.h>
59 # include <nvector/nvector_serial.h>
60 # include <ida/ida_spgmr.h>
61 # include <ida/ida_spbcgs.h>
62 # include <ida/ida_sptfqmr.h>
63 # include <ida/ida_dense.h>
64 # ifndef IDA_SUCCESS
65 # error "Failed to include SUNDIALS IDA header file"
66 # endif
67 #endif
68
69 #ifdef ASC_WITH_MMIO
70 # include <mmio.h>
71 #endif
72
73 /*
74 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
75 */
76 #ifndef SUNDIALS_VERSION_MINOR
77 # ifdef __GNUC__
78 # warning "GUESSING SUNDIALS VERSION 2.2"
79 # endif
80 # define SUNDIALS_VERSION_MINOR 2
81 #endif
82 #ifndef SUNDIALS_VERSION_MAJOR
83 # define SUNDIALS_VERSION_MAJOR 2
84 #endif
85
86 /* check that we've got what we expect now */
87 #ifndef ASC_IDA_H
88 # error "Failed to include ASCEND IDA header file"
89 #endif
90
91 /* #define FEX_DEBUG */
92 /* #define JEX_DEBUG */
93 /* #define SOLVE_DEBUG */
94 #define STATS_DEBUG
95 #define PREC_DEBUG
96
97 /**
98 Everthing that the outside world needs to know about IDA
99 */
100 const IntegratorInternals integrator_ida_internals = {
101 integrator_ida_create
102 ,integrator_ida_params_default
103 ,integrator_analyse_dae /* note, this routine is back in integrator.c */
104 ,integrator_ida_solve
105 ,integrator_ida_free
106 ,INTEG_IDA
107 ,"IDA"
108 };
109
110 /*-------------------------------------------------------------
111 FORWARD DECLS
112 */
113
114 /* forward dec needed for IntegratorIdaPrecFreeFn */
115 struct IntegratorIdaDataStruct;
116
117 /* functions for allocating storage for and freeing preconditioner data */
118 typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *blsys);
119 typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);
120
121 /**
122 Struct containing any stuff that IDA needs that doesn't fit into the
123 common IntegratorSystem struct.
124 */
125 typedef struct IntegratorIdaDataStruct{
126 struct rel_relation **rellist; /**< NULL terminated list of rels */
127 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
128 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
129 int nrels;
130 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
131 var_filter_t vfilter; /**< Used to filter variables from varlist in relman_diff2 etc */
132 rel_filter_t rfilter; /**< Used to filter relations from rellist (@TODO needs work) */
133 void *precdata; /**< For use by the preconditioner */
134 IntegratorIdaPrecFreeFn *pfree; /**< Store instructions here on how to free precdata */
135 } IntegratorIdaData;
136
137 typedef struct{
138 N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */
139 } IntegratorIdaPrecDataJacobi;
140
141 /**
142 Hold all the function pointers associated with a particular preconditioner
143 We don't need to store the 'pfree' function here as it is allocated to the enginedata struct
144 by the pcreate function (ensures that corresponding 'free' and 'create' are always used)
145
146 @note IDA uses a different convention for function pointer types, so no '*'.
147 */
148 typedef struct{
149 IntegratorIdaPrecCreateFn *pcreate;
150 IDASpilsPrecSetupFn psetup;
151 IDASpilsPrecSolveFn psolve;
152 } IntegratorIdaPrec;
153
154 /* residual function forward declaration */
155 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
156
157 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
158 , N_Vector v, N_Vector Jv, realtype c_j
159 , void *jac_data, N_Vector tmp1, N_Vector tmp2
160 );
161
162 /* error handler forward declaration */
163 void integrator_ida_error(int error_code
164 , const char *module, const char *function
165 , char *msg, void *eh_data
166 );
167
168 int integrator_ida_djex(long int Neq, realtype tt
169 , N_Vector yy, N_Vector yp, N_Vector rr
170 , realtype c_j, void *jac_data, DenseMat Jac
171 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
172 );
173
174 typedef struct{
175 long nsteps;
176 long nrevals;
177 long nlinsetups;
178 long netfails;
179 int qlast, qcur;
180 realtype hinused, hlast, hcur;
181 realtype tcur;
182 } IntegratorIdaStats;
183
184 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
185 void integrator_ida_write_stats(IntegratorIdaStats *stats);
186 void integrator_ida_write_incidence(IntegratorSystem *blsys);
187 /*------
188 Jacobi preconditioner -- experimental
189 */
190
191 int integrator_ida_psetup_jacobi(realtype tt,
192 N_Vector yy, N_Vector yp, N_Vector rr,
193 realtype c_j, void *prec_data,
194 N_Vector tmp1, N_Vector tmp2,
195 N_Vector tmp3
196 );
197
198 int integrator_ida_psolve_jacobi(realtype tt,
199 N_Vector yy, N_Vector yp, N_Vector rr,
200 N_Vector rvec, N_Vector zvec,
201 realtype c_j, realtype delta, void *prec_data,
202 N_Vector tmp
203 );
204
205 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys);
206
207 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);
208
209 static const IntegratorIdaPrec prec_jacobi = {
210 integrator_ida_pcreate_jacobi
211 , integrator_ida_psetup_jacobi
212 , integrator_ida_psolve_jacobi
213 };
214
215 /*-------------------------------------------------------------
216 SETUP/TEARDOWN ROUTINES
217 */
218 void integrator_ida_create(IntegratorSystem *blsys){
219 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
220 IntegratorIdaData *enginedata;
221 enginedata = ASC_NEW(IntegratorIdaData);
222 enginedata->rellist = NULL;
223 enginedata->varlist = NULL;
224 enginedata->safeeval = 0;
225 enginedata->vfilter.matchbits = VAR_SVAR | VAR_ACTIVE | VAR_FIXED;
226 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_ACTIVE;
227
228 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
229 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
230
231 blsys->enginedata = (void *)enginedata;
232
233 integrator_ida_params_default(blsys);
234 }
235
236 void integrator_ida_free(void *enginedata){
237 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
238 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
239 if(d->pfree){
240 /* free the preconditioner data, whatever it happens to be */
241 (d->pfree)(enginedata);
242 }
243 /* note, we don't own the rellist, so don't need to free it */
244 ASC_FREE(d);
245 }
246
247 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
248 IntegratorIdaData *d;
249 assert(blsys!=NULL);
250 assert(blsys->enginedata!=NULL);
251 assert(blsys->engine==INTEG_IDA);
252 d = ((IntegratorIdaData *)(blsys->enginedata));
253 return d;
254 }
255
256 /*-------------------------------------------------------------
257 PARAMETERS FOR IDA
258 */
259
260 enum ida_parameters{
261 IDA_PARAM_LINSOLVER
262 ,IDA_PARAM_MAXL
263 ,IDA_PARAM_AUTODIFF
264 ,IDA_PARAM_CALCIC
265 ,IDA_PARAM_SAFEEVAL
266 ,IDA_PARAM_RTOL
267 ,IDA_PARAM_ATOL
268 ,IDA_PARAM_ATOLVECT
269 ,IDA_PARAM_GSMODIFIED
270 ,IDA_PARAM_MAXNCF
271 ,IDA_PARAM_PREC
272 ,IDA_PARAMS_SIZE
273 };
274
275 /**
276 Here the full set of parameters is defined, along with upper/lower bounds,
277 etc. The values are stuck into the blsys->params structure.
278
279 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
280 list. Then add a slv_param_*(...) statement below to define the type, description and range
281 for the new parameter.
282
283 @return 0 on success
284 */
285 int integrator_ida_params_default(IntegratorSystem *blsys){
286 asc_assert(blsys!=NULL);
287 asc_assert(blsys->engine==INTEG_IDA);
288 slv_parameters_t *p;
289 p = &(blsys->params);
290
291 slv_destroy_parms(p);
292
293 if(p->parms==NULL){
294 CONSOLE_DEBUG("params NULL");
295 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
296 if(p->parms==NULL)return -1;
297 p->dynamic_parms = 1;
298 }else{
299 CONSOLE_DEBUG("params not NULL");
300 }
301
302 /* reset the number of parameters to zero so that we can check it at the end */
303 p->num_parms = 0;
304
305 slv_param_bool(p,IDA_PARAM_AUTODIFF
306 ,(SlvParameterInitBool){{"autodiff"
307 ,"Use auto-diff?",1
308 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
309 }, TRUE}
310 );
311
312 slv_param_char(p,IDA_PARAM_CALCIC
313 ,(SlvParameterInitChar){{"calcic"
314 ,"Initial conditions calcuation",1
315 ,"Use specified values of ydot to solve for inital y (Y),"
316 " or use the the values of the differential variables (yd) to solve"
317 " for the pure algebraic variables (ya) along with the derivatives"
318 " of the differential variables (yddot) (YA_YDP), or else don't solve"
319 " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
320 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE"}
321 );
322
323 slv_param_bool(p,IDA_PARAM_SAFEEVAL
324 ,(SlvParameterInitBool){{"safeeval"
325 ,"Use safe evaluation?",1
326 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
327 "throw SIGFPE errors which will then halt integration."
328 }, FALSE}
329 );
330
331
332 slv_param_bool(p,IDA_PARAM_ATOLVECT
333 ,(SlvParameterInitBool){{"atolvect"
334 ,"Use 'ode_atol' values as specified?",1
335 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
336 " in the integration. If FALSE, a scalar absolute tolerance value"
337 " is shared by all variables. See IDA manual, section 5.5.1"
338 }, TRUE }
339 );
340
341 slv_param_real(p,IDA_PARAM_ATOL
342 ,(SlvParameterInitReal){{"atol"
343 ,"Scalar absolute error tolerance",1
344 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
345 " See IDA manual, section 5.5.1"
346 }, 1e-5, DBL_MIN, DBL_MAX }
347 );
348
349 slv_param_real(p,IDA_PARAM_RTOL
350 ,(SlvParameterInitReal){{"rtol"
351 ,"Scalar relative error tolerance",1
352 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
353 " it's not possible to set per-variable relative tolerances as it is"
354 " with LSODE)."
355 " See IDA manual, section 5.5.1"
356 }, 1e-4, 0, DBL_MAX }
357 );
358
359 slv_param_char(p,IDA_PARAM_LINSOLVER
360 ,(SlvParameterInitChar){{"linsolver"
361 ,"Linear solver",1
362 ,"See IDA manual, section 5.5.3."
363 }, "SPGMR"}, (char *[]){"DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
364 );
365
366 slv_param_int(p,IDA_PARAM_MAXL
367 ,(SlvParameterInitInt){{"maxl"
368 ,"Maximum Krylov dimension",0
369 ,"The maximum dimension of Krylov space used by the linear solver"
370 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
371 " The default of 0 results in IDA using its internal default, which"
372 " is currently a value of 5."
373 }, 0, 0, 20 }
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 1 on success */
414 int integrator_ida_solve(
415 IntegratorSystem *blsys
416 , unsigned long start_index
417 , unsigned long finish_index
418 ){
419 void *ida_mem;
420 int size, 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 char *varname;
431 int i;
432 const IntegratorIdaPrec *prec = NULL;
433 int icopt; /* initial conditions strategy */
434
435 CONSOLE_DEBUG("STARTING IDA...");
436
437 enginedata = integrator_ida_enginedata(blsys);
438
439 enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);
440 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
441
442 /* store reference to list of relations (in enginedata) */
443 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
444 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
445 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
446 enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);
447
448 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
449 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
450 size = blsys->n_y;
451
452 if(enginedata->nrels!=size){
453 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
454 return 0; /* failure */
455 }
456
457 /* retrieve initial values from the system */
458
459 /** @TODO fix this, the starting time != first sample */
460 t0 = integrator_get_t(blsys);
461 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
462
463 CONSOLE_DEBUG("RETRIEVING y0");
464
465 y0 = N_VNew_Serial(size);
466 integrator_get_y(blsys,NV_DATA_S(y0));
467
468 #ifdef SOLVE_DEBUG
469 CONSOLE_DEBUG("RETRIEVING yp0");
470 #endif
471
472 yp0 = N_VNew_Serial(size);
473 integrator_get_ydot(blsys,NV_DATA_S(yp0));
474
475 #ifdef SOLVE_DEBUG
476 N_VPrint_Serial(yp0);
477 CONSOLE_DEBUG("yp0 is at %p",&yp0);
478 #endif
479
480 /* create IDA object */
481 ida_mem = IDACreate();
482
483 /* relative error tolerance */
484 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
485 CONSOLE_DEBUG("rtol = %8.2e",reltol);
486
487 /* allocate internal memory */
488 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
489 /* vector of absolute tolerances */
490 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
491 abstolvect = N_VNew_Serial(size);
492 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
493
494 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
495
496 N_VDestroy_Serial(abstolvect);
497 }else{
498 /* scalar absolute tolerance (one value for all) */
499 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
500 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
501 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
502 }
503
504 if(flag==IDA_MEM_NULL){
505 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
506 return 0;
507 }else if(flag==IDA_MEM_FAIL){
508 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
509 return 0;
510 }else if(flag==IDA_ILL_INPUT){
511 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
512 return 0;
513 }/* else success */
514
515 /* set optional inputs... */
516 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
517 IDASetRdata(ida_mem, (void *)blsys);
518 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
519 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
520 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
521 if(integrator_get_minstep(blsys)>0){
522 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
523 }
524
525 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
526 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
527
528 /* there's no capability for setting *minimum* step size in IDA */
529
530
531 /* attach linear solver module, using the default value of maxl */
532 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
533 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
534 if(strcmp(linsolver,"DENSE")==0){
535 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
536 flag = IDADense(ida_mem, size);
537 switch(flag){
538 case IDADENSE_SUCCESS: break;
539 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 0;
540 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 0;
541 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
542 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
543 }
544
545 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
546 CONSOLE_DEBUG("USING AUTODIFF");
547 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
548 switch(flag){
549 case IDADENSE_SUCCESS: break;
550 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 0;
551 }
552 }else{
553 CONSOLE_DEBUG("USING NUMERICAL DIFF");
554 }
555
556 flagfntype = "IDADENSE";
557 flagfn = &IDADenseGetLastFlag;
558 flagnamefn = &IDADenseGetReturnFlagName;
559 }else{
560 /* remaining methods are all SPILS */
561 CONSOLE_DEBUG("IDA SPILS");
562
563 maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
564 CONSOLE_DEBUG("maxl = %d",maxl);
565
566 /* what preconditioner? */
567 pname = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_PREC);
568 if(strcmp(pname,"NONE")==0){
569 prec = NULL;
570 }else if(strcmp(pname,"JACOBI")==0){
571 prec = &prec_jacobi;
572 }else{
573 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
574 return 0;
575 }
576
577 /* which SPILS linear solver? */
578 if(strcmp(linsolver,"SPGMR")==0){
579 CONSOLE_DEBUG("IDA SPGMR");
580 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
581 }else if(strcmp(linsolver,"SPBCG")==0){
582 CONSOLE_DEBUG("IDA SPBCG");
583 flag = IDASpbcg(ida_mem, maxl);
584 }else if(strcmp(linsolver,"SPTFQMR")==0){
585 CONSOLE_DEBUG("IDA SPTFQMR");
586 flag = IDASptfqmr(ida_mem,maxl);
587 }else{
588 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
589 return 0;
590 }
591
592 if(prec){
593 /* assign the preconditioner to the linear solver */
594 (prec->pcreate)(blsys);
595 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)blsys);
596 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
597 }else{
598 CONSOLE_DEBUG("No preconditioner");
599 }
600
601 flagfntype = "IDASPILS";
602 flagfn = &IDASpilsGetLastFlag;
603 flagnamefn = &IDASpilsGetReturnFlagName;
604
605 if(flag==IDASPILS_MEM_NULL){
606 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
607 return 0;
608 }else if(flag==IDASPILS_MEM_FAIL){
609 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
610 return 0;
611 }/* else success */
612
613 /* assign the J*v function */
614 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
615 CONSOLE_DEBUG("USING AUTODIFF");
616 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
617 if(flag==IDASPILS_MEM_NULL){
618 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
619 return 0;
620 }else if(flag==IDASPILS_LMEM_NULL){
621 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
622 return 0;
623 }/* else success */
624 }else{
625 CONSOLE_DEBUG("USING NUMERICAL DIFF");
626 }
627
628 if(strcmp(linsolver,"SPGMR")==0){
629 /* select Gram-Schmidt orthogonalisation */
630 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){
631 CONSOLE_DEBUG("USING MODIFIED GS");
632 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
633 if(flag!=IDASPILS_SUCCESS){
634 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
635 return 0;
636 }
637 }else{
638 CONSOLE_DEBUG("USING CLASSICAL GS");
639 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
640 if(flag!=IDASPILS_SUCCESS){
641 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
642 return 0;
643 }
644 }
645 }
646 }
647
648 /* set linear solver optional inputs...
649 ...nothing here at the moment...
650 */
651
652 /* calculate initial conditions */
653 icopt = 0;
654 if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"Y")==0){
655 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
656 icopt = IDA_Y_INIT;
657 asc_assert(icopt!=0);
658 }else if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
659 CONSOLE_DEBUG("Solving initial conditions using values of yd");
660 icopt = IDA_YA_YDP_INIT;
661 asc_assert(icopt!=0);
662 id = N_VNew_Serial(blsys->n_y);
663 for(i=0; i < blsys->n_y; ++i){
664 if(blsys->ydot[i] == NULL){
665 NV_Ith_S(id,i) = 0.0;
666 varname = var_make_name(blsys->system,blsys->y[i]);
667 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
668 ASC_FREE(varname);
669 }else{
670 CONSOLE_DEBUG("y[%d] is differential",i);
671 NV_Ith_S(id,i) = 1.0;
672 }
673 }
674 IDASetId(ida_mem, id);
675 N_VDestroy_Serial(id);
676 }else{
677 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
678 }
679
680 if(icopt){
681 blsys->currentstep=0;
682 t_index=start_index;
683 tout1 = samplelist_get(blsys->samples, t_index);
684
685 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
686
687 /* catch SIGFPE if desired to */
688 if(enginedata->safeeval){
689 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
690 }else{
691 #ifdef FEX_DEBUG
692 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
693 #endif
694 Asc_SignalHandlerPushDefault(SIGFPE);
695 }
696
697 /* correct initial values, given derivatives */
698 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
699 /* note the new API from version 2.3 and onwards */
700 flag = IDACalcIC(ida_mem, icopt, tout1);
701 # else
702 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
703 # endif
704
705 switch(flag){
706 case IDA_SUCCESS:
707 CONSOLE_DEBUG("Initial conditions solved OK");
708 break;
709
710 case IDA_LSETUP_FAIL:
711 case IDA_LINIT_FAIL:
712 case IDA_LSOLVE_FAIL:
713 case IDA_NO_RECOVERY:
714 flag1 = -999;
715 flag = (flagfn)(ida_mem,&flag1);
716 if(flag){
717 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);
718 return 0;
719 }
720 ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
721 return 0;
722
723 default:
724 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
725 return 0;
726 }
727
728 if(enginedata->safeeval){
729 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
730 }else{
731 Asc_SignalHandlerPopDefault(SIGFPE);
732 }
733
734 }
735
736 /* optionally, specify ROO-FINDING PROBLEM */
737
738 /* -- set up the IntegratorReporter */
739 integrator_output_init(blsys);
740
741 /* -- store the initial values of all the stuff */
742 integrator_output_write(blsys);
743 integrator_output_write_obs(blsys);
744
745 /* specify where the returned values should be stored */
746 yret = y0;
747 ypret = yp0;
748
749 /* advance solution in time, return values as yret and derivatives as ypret */
750 blsys->currentstep=1;
751 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
752 t = samplelist_get(blsys->samples, t_index);
753 t0 = integrator_get_t(blsys);
754 asc_assert(t > t0);
755
756 #ifdef SOLVE_DEBUG
757 CONSOLE_DEBUG("Integratoring from t0 = %f to t = %f", t0, t);
758 #endif
759
760 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
761
762 /* pass the values of everything back to the compiler */
763 integrator_set_t(blsys, (double)tret);
764 integrator_set_y(blsys, NV_DATA_S(yret));
765 integrator_set_ydot(blsys, NV_DATA_S(ypret));
766
767 if(flag<0){
768 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
769 break;
770 }
771
772 /* -- do something so that blsys knows the values of tret, yret and ypret */
773
774 /* -- store the current values of all the stuff */
775 integrator_output_write(blsys);
776 integrator_output_write_obs(blsys);
777 }
778
779 /* -- close the IntegratorReporter */
780 integrator_output_close(blsys);
781
782 /* get optional outputs */
783 #ifdef STATS_DEBUG
784 IntegratorIdaStats stats;
785 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
786 integrator_ida_write_stats(&stats);
787 }
788 #endif
789
790 /* free solution memory */
791 N_VDestroy_Serial(yret);
792 N_VDestroy_Serial(ypret);
793
794 /* free solver memory */
795 IDAFree(ida_mem);
796
797 if(flag < 0){
798 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
799 return 0;
800 }
801
802 /* all done, success */
803 return 1;
804 }
805
806 /*--------------------------------------------------
807 RESIDUALS AND JACOBIAN
808 */
809 /**
810 Function to evaluate system residuals, in the form required for IDA.
811
812 Given tt, yy and yp, we need to evaluate and return rr.
813
814 @param tt current value of indep variable (time)
815 @param yy current values of dependent variable vector
816 @param yp current values of derivatives of dependent variables
817 @param rr the output residual vector (is we're returning data to)
818 @param res_data pointer to our stuff (blsys in this case).
819
820 @return 0 on success, positive on recoverable error, and
821 negative on unrecoverable error.
822 */
823 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
824 IntegratorSystem *blsys;
825 IntegratorIdaData *enginedata;
826 int i, calc_ok, is_error;
827 struct rel_relation** relptr;
828 double resid;
829 char *relname;
830 #ifdef FEX_DEBUG
831 char *varname;
832 char diffname[30];
833 #endif
834
835 blsys = (IntegratorSystem *)res_data;
836 enginedata = integrator_ida_enginedata(blsys);
837
838 #ifdef FEX_DEBUG
839 /* fprintf(stderr,"\n\n"); */
840 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
841 #endif
842
843 /* pass the values of everything back to the compiler */
844 integrator_set_t(blsys, (double)tt);
845 integrator_set_y(blsys, NV_DATA_S(yy));
846 integrator_set_ydot(blsys, NV_DATA_S(yp));
847
848 if(NV_LENGTH_S(rr)!=enginedata->nrels){
849 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
850 return -1; /* unrecoverable */
851 }
852
853 /**
854 @TODO does this function (fex) do bounds checking already?
855 */
856
857 /* evaluate each residual in the rellist */
858 is_error = 0;
859 relptr = enginedata->rellist;
860
861 if(enginedata->safeeval){
862 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
863 }else{
864 #ifdef FEX_DEBUG
865 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
866 #endif
867 Asc_SignalHandlerPushDefault(SIGFPE);
868 }
869
870 if (setjmp(g_fpe_env)==0) {
871 for(i=0, relptr = enginedata->rellist;
872 i< enginedata->nrels && relptr != NULL;
873 ++i, ++relptr
874 ){
875 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
876
877 NV_Ith_S(rr,i) = resid;
878 if(!calc_ok){
879 relname = rel_make_name(blsys->system, *relptr);
880 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
881 ASC_FREE(relname);
882 /* presumable some output already made? */
883 is_error = 1;
884 }/*else{
885 CONSOLE_DEBUG("Calc OK");
886 }*/
887 }
888 }else{
889 relname = rel_make_name(blsys->system, *relptr);
890 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
891 ASC_FREE(relname);
892 is_error = 1;
893 }
894
895 if(enginedata->safeeval){
896 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
897 }else{
898 Asc_SignalHandlerPopDefault(SIGFPE);
899 }
900
901 #ifdef FEX_DEBUG
902 /* output residuals to console */
903 CONSOLE_DEBUG("RESIDUAL OUTPUT");
904 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
905 for(i=0; i<blsys->n_y; ++i){
906 varname = var_make_name(blsys->system,blsys->y[i]);
907 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
908 if(blsys->ydot[i]){
909 varname = var_make_name(blsys->system,blsys->ydot[i]);
910 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
911 }else{
912 snprintf(diffname,99,"diff(%s)",varname);
913 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
914 }
915 ASC_FREE(varname);
916 relname = rel_make_name(blsys->system,enginedata->rellist[i]);
917 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
918 }
919 #endif
920
921 if(is_error){
922 return 1;
923 }
924
925 #ifdef FEX_DEBUG
926 CONSOLE_DEBUG("RESIDUAL OK");
927 #endif
928 return 0;
929 }
930
931 /**
932 Dense Jacobian evaluation. Only suitable for small problems!
933 */
934 int integrator_ida_djex(long int Neq, realtype tt
935 , N_Vector yy, N_Vector yp, N_Vector rr
936 , realtype c_j, void *jac_data, DenseMat Jac
937 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
938 ){
939 IntegratorSystem *blsys;
940 IntegratorIdaData *enginedata;
941 char *relname;
942 #ifdef JEX_DEBUG
943 char *varname;
944 #endif
945 int status;
946 struct rel_relation **relptr;
947 int i;
948 double *derivatives;
949 int *variables;
950 int count, j;
951 long var_yindex;
952
953 blsys = (IntegratorSystem *)jac_data;
954 enginedata = integrator_ida_enginedata(blsys);
955
956 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
957 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
958 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
959
960 /* pass the values of everything back to the compiler */
961 integrator_set_t(blsys, (double)tt);
962 integrator_set_y(blsys, NV_DATA_S(yy));
963 integrator_set_ydot(blsys, NV_DATA_S(yp));
964
965 #ifdef JEX_DEBUG
966 /* print vars */
967 for(i=0; i < blsys->n_y; ++i){
968 varname = var_make_name(blsys->system, blsys->y[i]);
969 CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
970 ASC_FREE(varname);
971 }
972
973 /* print derivatives */
974 for(i=0; i < blsys->n_y; ++i){
975 if(blsys->ydot[i]){
976 varname = var_make_name(blsys->system, blsys->ydot[i]);
977 CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
978 ASC_FREE(varname);
979 }else{
980 varname = var_make_name(blsys->system, blsys->y[i]);
981 CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
982 ASC_FREE(varname);
983 }
984 }
985
986 /* print step size */
987 CONSOLE_DEBUG("<c_j> = %f",c_j);
988 #endif
989
990 /* build up the dense jacobian matrix... */
991 status = 0;
992 for(i=0, relptr = enginedata->rellist;
993 i< enginedata->nrels && relptr != NULL;
994 ++i, ++relptr
995 ){
996 /* get derivatives for this particular relation */
997 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
998
999 if(status){
1000 relname = rel_make_name(blsys->system, *relptr);
1001 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1002 ASC_FREE(relname);
1003 break;
1004 }
1005
1006 /* output what's going on here ... */
1007 #ifdef JEX_DEBUG
1008 relname = rel_make_name(blsys->system, *relptr);
1009 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
1010 fprintf(stderr,"%d: '%s': ",i,relname);
1011 ASC_FREE(relname);
1012 for(j=0;j<count;++j){
1013 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1014 var_yindex = blsys->y_id[variables[j]];
1015 if(var_yindex >=0){
1016 fprintf(stderr," var[%d]='%s'=y[%ld]",variables[j],varname,var_yindex);
1017 }else{
1018 fprintf(stderr," var[%d]='%s'=ydot[%ld]",variables[j],varname,-var_yindex-1);
1019 }
1020 ASC_FREE(varname);
1021 }
1022 fprintf(stderr,"\n");
1023 #endif
1024 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1025 for(j=0; j < count; ++j){
1026 var_yindex = blsys->y_id[variables[j]];
1027 /* the SUNDIALS headers seem not to store 'N' on Windows */
1028 ASC_ASSERT_RANGE(var_yindex, -blsys->n_y, blsys->n_y);
1029
1030 if(var_yindex >= 0){
1031 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
1032 DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
1033 }else{
1034 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
1035 DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
1036 }
1037 }
1038 }
1039
1040 #ifdef JEX_DEBUG
1041 CONSOLE_DEBUG("PRINTING JAC");
1042 fprintf(stderr,"\t");
1043 for(j=0; j < blsys->n_y; ++j){
1044 if(j)fprintf(stderr,"\t");
1045 varname = var_make_name(blsys->system,blsys->y[j]);
1046 fprintf(stderr,"%11s",varname);
1047 ASC_FREE(varname);
1048 }
1049 fprintf(stderr,"\n");
1050 for(i=0; i < enginedata->nrels; ++i){
1051 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
1052 fprintf(stderr,"%s\t",relname);
1053 ASC_FREE(relname);
1054
1055 for(j=0; j < blsys->n_y; ++j){
1056 if(j)fprintf(stderr,"\t");
1057 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1058 }
1059 fprintf(stderr,"\n");
1060 }
1061 #endif
1062
1063 if(status){
1064 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1065 return 1;
1066 }
1067
1068 #ifdef FEX_DEBUG
1069 CONSOLE_DEBUG("DJEX RETURNING 0");
1070 #endif
1071 return 0;
1072 }
1073
1074 /**
1075 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1076
1077 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1078
1079 @param tt current value of the independent variable (time, t)
1080 @param yy current value of the dependent variable vector, y(t).
1081 @param yp current value of y'(t).
1082 @param rr current value of the residual vector F(t, y, y').
1083 @param v the vector by which the Jacobian must be multiplied to the right.
1084 @param Jv the output vector computed
1085 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1086 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1087 @param tmp1 @see tmp2
1088 @param tmp2 (as well as tmp1) pointers to memory allocated for variables of type N_Vector for use here as temporary storage or work space.
1089 @return 0 on success
1090 */
1091 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1092 , N_Vector v, N_Vector Jv, realtype c_j
1093 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1094 ){
1095 IntegratorSystem *blsys;
1096 IntegratorIdaData *enginedata;
1097 int i, j, is_error=0;
1098 struct rel_relation** relptr;
1099 char *relname;
1100 int status;
1101 double Jv_i;
1102 long var_yindex;
1103
1104 int *variables;
1105 double *derivatives;
1106 int count;
1107 #ifdef JEX_DEBUG
1108 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1109 #endif
1110
1111 blsys = (IntegratorSystem *)jac_data;
1112 enginedata = integrator_ida_enginedata(blsys);
1113
1114 /* pass the values of everything back to the compiler */
1115 integrator_set_t(blsys, (double)tt);
1116 integrator_set_y(blsys, NV_DATA_S(yy));
1117 integrator_set_ydot(blsys, NV_DATA_S(yp));
1118 /* no real use for residuals (rr) here, I don't think? */
1119
1120 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1121
1122 i = NV_LENGTH_S(yy) * 2;
1123 #ifdef JEX_DEBUG
1124 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1125 #endif
1126 variables = ASC_NEW_ARRAY(int, i);
1127 derivatives = ASC_NEW_ARRAY(double, i);
1128
1129 /* evaluate the derivatives... */
1130 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1131
1132 Asc_SignalHandlerPushDefault(SIGFPE);
1133 if (setjmp(g_fpe_env)==0) {
1134 for(i=0, relptr = enginedata->rellist;
1135 i< enginedata->nrels && relptr != NULL;
1136 ++i, ++relptr
1137 ){
1138 /* get derivatives for this particular relation */
1139 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1140 #ifdef JEX_DEBUG
1141 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1142 #endif
1143
1144 if(status){
1145 relname = rel_make_name(blsys->system, *relptr);
1146 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1147 ASC_FREE(relname);
1148 is_error = 1;
1149 break;
1150 }
1151
1152 /*
1153 Now we have the derivatives wrt each alg/diff variable in the
1154 present equation. variables[] points into the varlist. need
1155 a mapping from the varlist to the y and ydot lists.
1156 */
1157
1158 Jv_i = 0;
1159 for(j=0; j < count; ++j){
1160 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
1161 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1162 if(varname){
1163 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1164 ASC_FREE(varname);
1165 }else{
1166 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1167 }
1168 */
1169
1170 /* we don't calculate derivatives wrt indep var */
1171 asc_assert(variables[j]>=0);
1172 if(enginedata->varlist[variables[j]] == blsys->x) continue;
1173
1174 var_yindex = blsys->y_id[variables[j]];
1175 #ifdef JEX_DEBUG
1176 CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %ld",j,variables[j],var_yindex);
1177 #endif
1178
1179 ASC_ASSERT_RANGE(-var_yindex-1, -NV_LENGTH_S(v),NV_LENGTH_S(v));
1180
1181 if(var_yindex >= 0){
1182 #ifdef JEX_DEBUG
1183 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
1184 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%ld] = %f, v[%ld] = %f)\n", i
1185 , derivatives[j] * NV_Ith_S(v,var_yindex)
1186 , i, var_yindex, derivatives[j]
1187 , var_yindex, NV_Ith_S(v,var_yindex)
1188 );
1189 #endif
1190 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
1191 }else{
1192 #ifdef JEX_DEBUG
1193 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%ld] = %f, v[%ld] = %f)\n", i
1194 , derivatives[j] * NV_Ith_S(v,-var_yindex-1)
1195 , i, -var_yindex-1, derivatives[j]
1196 , -var_yindex-1, NV_Ith_S(v,-var_yindex-1)
1197 );
1198 #endif
1199 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
1200 Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
1201 }
1202 }
1203
1204 NV_Ith_S(Jv,i) = Jv_i;
1205 #ifdef JEX_DEBUG
1206 CONSOLE_DEBUG("rel = %p",*relptr);
1207 relname = rel_make_name(blsys->system, *relptr);
1208 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1209 //ASC_FREE(relname);
1210 return 1;
1211 #endif
1212 }
1213 }else{
1214 relname = rel_make_name(blsys->system, *relptr);
1215 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1216 ASC_FREE(relname);
1217 is_error = 1;
1218 }
1219 Asc_SignalHandlerPopDefault(SIGFPE);
1220
1221 if(is_error){
1222 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1223 return 1;
1224 }
1225 return 0;
1226 }
1227
1228 /*----------------------------------------------
1229 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1230 */
1231
1232 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys){
1233 IntegratorIdaData * enginedata =blsys->enginedata;
1234 IntegratorIdaPrecDataJacobi *precdata;
1235 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1236
1237 asc_assert(blsys->n_y);
1238 precdata->PIii = N_VNew_Serial(blsys->n_y);
1239
1240 enginedata->pfree = &integrator_ida_pfree_jacobi;
1241 enginedata->precdata = precdata;
1242 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1243 }
1244
1245 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1246 if(enginedata->precdata){
1247 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1248 N_VDestroy_Serial(precdata->PIii);
1249
1250 ASC_FREE(precdata);
1251 enginedata->precdata = NULL;
1252 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1253 }
1254 enginedata->pfree = NULL;
1255 }
1256
1257 /**
1258 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1259
1260 'setup' function.
1261 */
1262 int integrator_ida_psetup_jacobi(realtype tt,
1263 N_Vector yy, N_Vector yp, N_Vector rr,
1264 realtype c_j, void *p_data,
1265 N_Vector tmp1, N_Vector tmp2,
1266 N_Vector tmp3
1267 ){
1268 int i, j;
1269 IntegratorSystem *blsys;
1270 IntegratorIdaData *enginedata;
1271 IntegratorIdaPrecDataJacobi *precdata;
1272 struct rel_relation **relptr;
1273
1274 blsys = (IntegratorSystem *)p_data;
1275 enginedata = blsys->enginedata;
1276 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1277 double *derivatives;
1278 int *variables;
1279 int count, status;
1280 char *relname;
1281 int var_yindex;
1282
1283 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1284
1285 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
1286 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1287
1288 /**
1289 @TODO FIXME here we are using the very inefficient and contorted approach
1290 of calculating the whole jacobian, then extracting just the diagonal elements.
1291 */
1292
1293 for(i=0, relptr = enginedata->rellist;
1294 i< enginedata->nrels && relptr != NULL;
1295 ++i, ++relptr
1296 ){
1297
1298 /* get derivatives for this particular relation */
1299 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1300 if(status){
1301 relname = rel_make_name(blsys->system, *relptr);
1302 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1303 ASC_FREE(relname);
1304 break;
1305 }
1306 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1307 /* find the diagonal elements */
1308 for(j=0; j<count; ++j){
1309 if(variables[j]==i){
1310 var_yindex = blsys->y_id[variables[j]];
1311 if(var_yindex >= 0){
1312 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1313 }else{
1314 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
1315 }
1316 }
1317 }
1318 #ifdef PREC_DEBUG
1319 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1320 #endif
1321 }
1322
1323 if(status){
1324 CONSOLE_DEBUG("Error found when evaluating derivatives");
1325 return 1; /* recoverable */
1326 }
1327
1328 integrator_ida_write_incidence(blsys);
1329
1330 ASC_FREE(variables);
1331 ASC_FREE(derivatives);
1332
1333 return 0;
1334 };
1335
1336 /**
1337 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1338
1339 'solve' function.
1340 */
1341 int integrator_ida_psolve_jacobi(realtype tt,
1342 N_Vector yy, N_Vector yp, N_Vector rr,
1343 N_Vector rvec, N_Vector zvec,
1344 realtype c_j, realtype delta, void *p_data,
1345 N_Vector tmp
1346 ){
1347 IntegratorSystem *blsys;
1348 IntegratorIdaData *data;
1349 IntegratorIdaPrecDataJacobi *precdata;
1350 blsys = (IntegratorSystem *)p_data;
1351 data = blsys->enginedata;
1352 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1353
1354 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1355 N_VProd(precdata->PIii, rvec, zvec);
1356 return 0;
1357 };
1358
1359 /*----------------------------------------------
1360 STATS
1361 */
1362
1363 /**
1364 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1365 status in a struct instead of separately.
1366
1367 @return IDA_SUCCESS on sucess.
1368 */
1369 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1370 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1371 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1372 ,&s->hlast, &s->hcur, &s->tcur
1373 );
1374 }
1375
1376 /**
1377 This routine just outputs the stats to the CONSOLE_DEBUG routine.
1378
1379 @TODO provide a GUI way of stats reporting from IDA.
1380 */
1381 void integrator_ida_write_stats(IntegratorIdaStats *stats){
1382 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
1383 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
1384 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
1385 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
1386 SI(qlast); SI(qcur);
1387 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
1388 # undef SL
1389 # undef SI
1390 # undef SR
1391 }
1392
1393 /*------------------------------------------------------------------------------
1394 JACOBIAN / INCIDENCE MATRIX OUTPUT
1395 */
1396
1397 enum integrator_ida_write_jac_enum{
1398 II_WRITE_Y
1399 , II_WRITE_YDOT
1400 };
1401
1402 /**
1403 @TODO COMPLETE THIS...
1404 */
1405 void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){
1406 IntegratorIdaData *enginedata;
1407 MM_typecode matcode;
1408 int nnz, rhomax;
1409 double *derivatives;
1410 int *variables;
1411 struct rel_relation **relptr;
1412 int i, j, status, count, var_yindex;
1413 char *relname;
1414
1415 var_filter_t vfiltery = {
1416 VAR_SVAR | VAR_FIXED | VAR_DERIV
1417 , VAR_SVAR
1418 };
1419 var_filter_t vfilteryd = {
1420 VAR_SVAR | VAR_FIXED | VAR_DERIV
1421 , VAR_SVAR | VAR_DERIV
1422 };
1423
1424 enginedata = (IntegratorIdaData *)blsys->enginedata;
1425
1426 /* number of non-zeros for all the non-FIXED solver_vars,
1427 in all the active included equality relations.
1428 */
1429 nnz = relman_jacobian_count(enginedata->rellist, enginedata->nrels
1430 , &enginedata->vfilter, &enginedata->rfilter
1431 , &rhomax
1432 );
1433
1434 /* we must have found the same number of relations */
1435 asc_assert(rhomax == enginedata->nrels);
1436
1437 /* output the mmio file header, now that we know our size*/
1438 /* note that we are asserting that our problem is square */
1439 mm_initialize_typecode(&matcode);
1440 mm_set_matrix(&matcode);
1441 mm_set_coordinate(&matcode);
1442 mm_set_real(&matcode);
1443 mm_write_banner(f, matcode);
1444 mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);
1445
1446 variables = ASC_NEW_ARRAY(int, blsys->n_y * 2);
1447 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1448
1449 CONSOLE_DEBUG("Writing sparse Jacobian to file...");
1450
1451 for(i=0, relptr = enginedata->rellist;
1452 i< enginedata->nrels && relptr != NULL;
1453 ++i, ++relptr
1454 ){
1455 relname = rel_make_name(blsys->system, *relptr);
1456
1457 /* get derivatives of y */
1458 status = relman_diff2(*relptr, &vfiltery, derivatives, variables, &count, enginedata->safeeval);
1459 if(status){
1460 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1461 ASC_FREE(relname);
1462 break;
1463 }
1464
1465 /* get derivatives of y */
1466 status = relman_diff2(*relptr, &vfilteryd, derivatives, variables, &count, enginedata->safeeval);
1467 if(status){
1468 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1469 ASC_FREE(relname);
1470 break;
1471 }
1472
1473 for(j=0; j<count; ++j){
1474 var_yindex = blsys->y_id[variables[j]];
1475 if(var_yindex >= 0 && type == II_WRITE_Y){
1476 fprintf(f, "%d %d %10.3g\n", i, var_yindex, derivatives[j]);
1477 }else if(var_yindex < 0 && type == II_WRITE_YDOT){
1478 fprintf(f, "%d %d %10.3g\n", i, -var_yindex-1, derivatives[j]);
1479 }
1480 }
1481 }
1482 ASC_FREE(variables);
1483 ASC_FREE(derivatives);
1484 }
1485
1486 /**
1487 This routine outputs matrix structure in a crude text format, for the sake
1488 of debugging.
1489 */
1490 void integrator_ida_write_incidence(IntegratorSystem *blsys){
1491 int i, j;
1492 struct rel_relation **relptr;
1493 IntegratorIdaData *enginedata = blsys->enginedata;
1494 double *derivatives;
1495 int *variables;
1496 int count, status;
1497 char *relname;
1498 int var_yindex;
1499
1500 if(enginedata->nrels > 100){
1501 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
1502 return;
1503 }
1504
1505 variables = ASC_NEW_ARRAY(int, blsys->n_y * 2);
1506 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1507
1508 CONSOLE_DEBUG("Outputting incidence information to console...");
1509
1510 for(i=0, relptr = enginedata->rellist;
1511 i< enginedata->nrels && relptr != NULL;
1512 ++i, ++relptr
1513 ){
1514 relname = rel_make_name(blsys->system, *relptr);
1515
1516 /* get derivatives for this particular relation */
1517 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1518 if(status){
1519 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1520 ASC_FREE(relname);
1521 break;
1522 }
1523
1524 fprintf(stderr,"%3d:%-15s:",i,relname);
1525 ASC_FREE(relname);
1526
1527 for(j=0; j<count; ++j){
1528 var_yindex = blsys->y_id[variables[j]];
1529 if(var_yindex >= 0){
1530 fprintf(stderr," %d:y[%d]",variables[j],var_yindex);
1531 }else{
1532 fprintf(stderr," %d:ydot[%d]",variables[j],-var_yindex-1);
1533 }
1534 }
1535 fprintf(stderr,"\n");
1536 }
1537 ASC_FREE(variables);
1538 ASC_FREE(derivatives);
1539 }
1540
1541 /*----------------------------------------------
1542 ERROR REPORTING
1543 */
1544 /**
1545 Error message reporter function to be passed to IDA. All error messages
1546 will trigger a call to this function, so we should find everything
1547 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1548 panel (in the case of PyGTK).
1549 */
1550 void integrator_ida_error(int error_code
1551 , const char *module, const char *function
1552 , char *msg, void *eh_data
1553 ){
1554 IntegratorSystem *blsys;
1555 error_severity_t sev;
1556
1557 /* cast back the IntegratorSystem, just in case we need it */
1558 blsys = (IntegratorSystem *)eh_data;
1559
1560 /* severity depends on the sign of the error_code value */
1561 if(error_code <= 0){
1562 sev = ASC_PROG_ERR;
1563 }else{
1564 sev = ASC_PROG_WARNING;
1565 }
1566
1567 /* use our all-purpose error reporting to get stuff back to the GUI */
1568 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1569 }

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