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

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