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

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