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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1387 - (show annotations) (download) (as text)
Sat Apr 7 14:43:31 2007 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 60471 byte(s)
Added plot support in Integrator output tabs.
Some other minor debugging for pylab integration and idaanalyse output.
1 /* ASCEND modelling environment
2 Copyright (C) 2006-2007 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
21 as part of the GPL-licensed SUNDIALS solver package from LLNL.
22
23 IDA provides the non-linear parts, as well as a number of pluggable linear
24 solvers: dense, banded and krylov types.
25
26 We also implement here an EXPERIMENTAL direct sparse linear solver for IDA
27 using the ASCEND linsolqr routines.
28
29 @see http://www.llnl.gov/casc/sundials/
30 *//*
31 by John Pye, May 2006
32 */
33
34 #define _GNU_SOURCE
35
36 #include "ida.h"
37
38 #include <signal.h>
39 #include <setjmp.h>
40 #include <fenv.h>
41 #include <math.h>
42
43 /* SUNDIALS includes */
44 #ifdef ASC_WITH_IDA
45 # include <sundials/sundials_config.h>
46 # include <sundials/sundials_dense.h>
47 # include <ida/ida.h>
48 # include <nvector/nvector_serial.h>
49 # include <ida/ida_spgmr.h>
50 # include <ida/ida_spbcgs.h>
51 # include <ida/ida_sptfqmr.h>
52 # include <ida/ida_dense.h>
53 # ifndef IDA_SUCCESS
54 # error "Failed to include SUNDIALS IDA header file"
55 # endif
56 #else
57 # error "If you're building this file, you should have ASC_WITH_IDA"
58 #endif
59
60 #ifdef ASC_WITH_MMIO
61 # include <mmio.h>
62 #endif
63
64 #include <utilities/ascConfig.h>
65 #include <utilities/error.h>
66 #include <utilities/ascSignal.h>
67 #include <utilities/ascPanic.h>
68 #include <compiler/instance_enum.h>
69
70 #include <system/slv_client.h>
71 #include <system/relman.h>
72 #include <system/block.h>
73 #include <system/slv_stdcalls.h>
74 #include <system/jacobian.h>
75 #include <system/bndman.h>
76
77 #include "idalinear.h"
78 #include "idaanalyse.h"
79 #include "ida_impl.h"
80
81 /*
82 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
83 */
84 #ifndef SUNDIALS_VERSION_MINOR
85 # ifdef __GNUC__
86 # warning "GUESSING SUNDIALS VERSION 2.2"
87 # endif
88 # define SUNDIALS_VERSION_MINOR 2
89 #endif
90 #ifndef SUNDIALS_VERSION_MAJOR
91 # define SUNDIALS_VERSION_MAJOR 2
92 #endif
93
94 /* check that we've got what we expect now */
95 #ifndef ASC_IDA_H
96 # error "Failed to include ASCEND IDA header file"
97 #endif
98
99 /* #define FEX_DEBUG */
100 #define JEX_DEBUG
101 /* #define DJEX_DEBUG */
102 /* #define SOLVE_DEBUG */
103 #define STATS_DEBUG
104 #define PREC_DEBUG
105 /* #define ROOT_DEBUG */
106
107 /* #define DIFFINDEX_DEBUG */
108 /* #define ANALYSE_DEBUG */
109 /* #define DESTROY_DEBUG */
110 /* #define MATRIX_DEBUG */
111
112 /**
113 Everthing that the outside world needs to know about IDA
114 */
115 const IntegratorInternals integrator_ida_internals = {
116 integrator_ida_create
117 ,integrator_ida_params_default
118 ,integrator_ida_analyse
119 ,integrator_ida_solve
120 ,integrator_ida_write_matrix
121 ,integrator_ida_debug
122 ,integrator_ida_free
123 ,INTEG_IDA
124 ,"IDA"
125 };
126
127 /*-------------------------------------------------------------
128 FORWARD DECLS
129 */
130
131 /* forward dec needed for IntegratorIdaPrecFreeFn */
132 struct IntegratorIdaDataStruct;
133
134 /* functions for allocating storage for and freeing preconditioner data */
135 typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *sys);
136 typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);
137
138
139 /**
140 Struct containing any stuff that IDA needs that doesn't fit into the
141 common IntegratorSystem struct.
142 */
143 typedef struct IntegratorIdaDataStruct{
144
145 struct rel_relation **rellist; /**< NULL terminated list of ACTIVE rels */
146 int nrels; /* number of ACTIVE rels */
147
148 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
149 int nbnds; /* number of boundaries */
150
151 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
152 var_filter_t vfilter;
153 rel_filter_t rfilter; /**< Used to filter relations from solver's rellist (@TODO needs work) */
154 void *precdata; /**< For use by the preconditioner */
155 IntegratorIdaPrecFreeFn *pfree; /**< Store instructions here on how to free precdata */
156
157 } IntegratorIdaData;
158
159
160 typedef struct IntegratorIdaPrecDJStruct{
161 N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */
162 } IntegratorIdaPrecDataJacobi;
163
164 /**
165 Hold all the function pointers associated with a particular preconditioner
166 We don't need to store the 'pfree' function here as it is allocated to the enginedata struct
167 by the pcreate function (ensures that corresponding 'free' and 'create' are always used)
168
169 @note IDA uses a different convention for function pointer types, so no '*'.
170 */
171 typedef struct IntegratorIdaPrecStruct{
172 IntegratorIdaPrecCreateFn *pcreate;
173 IDASpilsPrecSetupFn psetup;
174 IDASpilsPrecSolveFn psolve;
175 } IntegratorIdaPrec;
176
177 /* residual function forward declaration */
178 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
179
180 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
181 , N_Vector v, N_Vector Jv, realtype c_j
182 , void *jac_data, N_Vector tmp1, N_Vector tmp2
183 );
184
185 /* error handler forward declaration */
186 void integrator_ida_error(int error_code
187 , const char *module, const char *function
188 , char *msg, void *eh_data
189 );
190
191 /* dense jacobian evaluation for IDADense dense direct linear solver */
192 int integrator_ida_djex(long int Neq, realtype tt
193 , N_Vector yy, N_Vector yp, N_Vector rr
194 , realtype c_j, void *jac_data, DenseMat Jac
195 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
196 );
197
198 /* sparse jacobian evaluation for ASCEND's sparse direct solver */
199 IntegratorSparseJacFn integrator_ida_sjex;
200
201 /* boundary-detection function */
202 int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
203
204 typedef struct IntegratorIdaStatsStruct{
205 long nsteps;
206 long nrevals;
207 long nlinsetups;
208 long netfails;
209 int qlast, qcur;
210 realtype hinused, hlast, hcur;
211 realtype tcur;
212 } IntegratorIdaStats;
213
214 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
215
216 /*static IntegratorVarVisitorFn integrator_dae_classify_var;
217 static void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
218 static void integrator_dae_show_var(IntegratorSystem *sys, struct var_variable *var, const int *varindx); */
219
220 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
221 void integrator_ida_write_stats(IntegratorIdaStats *stats);
222 void integrator_ida_write_incidence(IntegratorSystem *sys);
223
224 /*------
225 Jacobi preconditioner -- experimental
226 */
227
228 int integrator_ida_psetup_jacobi(realtype tt,
229 N_Vector yy, N_Vector yp, N_Vector rr,
230 realtype c_j, void *prec_data,
231 N_Vector tmp1, N_Vector tmp2,
232 N_Vector tmp3
233 );
234
235 int integrator_ida_psolve_jacobi(realtype tt,
236 N_Vector yy, N_Vector yp, N_Vector rr,
237 N_Vector rvec, N_Vector zvec,
238 realtype c_j, realtype delta, void *prec_data,
239 N_Vector tmp
240 );
241
242 void integrator_ida_pcreate_jacobi(IntegratorSystem *sys);
243
244 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);
245
246 static const IntegratorIdaPrec prec_jacobi = {
247 integrator_ida_pcreate_jacobi
248 , integrator_ida_psetup_jacobi
249 , integrator_ida_psolve_jacobi
250 };
251
252 /*-------------------------------------------------------------
253 SETUP/TEARDOWN ROUTINES
254 */
255 void integrator_ida_create(IntegratorSystem *sys){
256 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
257 IntegratorIdaData *enginedata;
258 enginedata = ASC_NEW(IntegratorIdaData);
259 CONSOLE_DEBUG("enginedata = %p",enginedata);
260 enginedata->rellist = NULL;
261 enginedata->safeeval = 0;
262 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED;
263 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;
264 enginedata->pfree = NULL;
265
266 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
267 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
268
269 sys->enginedata = (void *)enginedata;
270
271 integrator_ida_params_default(sys);
272 }
273
274 void integrator_ida_free(void *enginedata){
275 #ifdef DESTROY_DEBUG
276 CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);
277 #endif
278 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
279 asc_assert(d);
280 if(d->pfree){
281 CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);
282 /* free the preconditioner data, whatever it happens to be */
283 (d->pfree)(enginedata);
284 }
285
286 ASC_FREE(d->rellist);
287
288 #ifdef DESTROY_DEBUG
289 CONSOLE_DEBUG("Now destroying the enginedata");
290 #endif
291 ASC_FREE(d);
292 #ifdef DESTROY_DEBUG
293 CONSOLE_DEBUG("enginedata freed");
294 #endif
295 }
296
297 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *sys){
298 IntegratorIdaData *d;
299 assert(sys!=NULL);
300 assert(sys->enginedata!=NULL);
301 assert(sys->engine==INTEG_IDA);
302 d = ((IntegratorIdaData *)(sys->enginedata));
303 return d;
304 }
305
306 /*-------------------------------------------------------------
307 PARAMETERS FOR IDA
308 */
309
310 enum ida_parameters{
311 IDA_PARAM_LINSOLVER
312 ,IDA_PARAM_MAXL
313 ,IDA_PARAM_MAXORD
314 ,IDA_PARAM_AUTODIFF
315 ,IDA_PARAM_CALCIC
316 ,IDA_PARAM_SAFEEVAL
317 ,IDA_PARAM_RTOL
318 ,IDA_PARAM_ATOL
319 ,IDA_PARAM_ATOLVECT
320 ,IDA_PARAM_GSMODIFIED
321 ,IDA_PARAM_MAXNCF
322 ,IDA_PARAM_PREC
323 ,IDA_PARAMS_SIZE
324 };
325
326 /**
327 Here the full set of parameters is defined, along with upper/lower bounds,
328 etc. The values are stuck into the sys->params structure.
329
330 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
331 list. Then add a slv_param_*(...) statement below to define the type, description and range
332 for the new parameter.
333
334 @return 0 on success
335 */
336 int integrator_ida_params_default(IntegratorSystem *sys){
337 asc_assert(sys!=NULL);
338 asc_assert(sys->engine==INTEG_IDA);
339 slv_parameters_t *p;
340 p = &(sys->params);
341
342 slv_destroy_parms(p);
343
344 if(p->parms==NULL){
345 CONSOLE_DEBUG("params NULL");
346 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
347 if(p->parms==NULL)return -1;
348 p->dynamic_parms = 1;
349 }else{
350 CONSOLE_DEBUG("params not NULL");
351 }
352
353 /* reset the number of parameters to zero so that we can check it at the end */
354 p->num_parms = 0;
355
356 slv_param_bool(p,IDA_PARAM_AUTODIFF
357 ,(SlvParameterInitBool){{"autodiff"
358 ,"Use auto-diff?",1
359 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
360 }, TRUE}
361 );
362
363 slv_param_char(p,IDA_PARAM_CALCIC
364 ,(SlvParameterInitChar){{"calcic"
365 ,"Initial conditions calcuation",1
366 ,"Use specified values of ydot to solve for inital y (Y),"
367 " or use the the values of the differential variables (yd) to solve"
368 " for the pure algebraic variables (ya) along with the derivatives"
369 " of the differential variables (yddot) (YA_YDP), or else don't solve"
370 " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
371 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
372 );
373
374 slv_param_bool(p,IDA_PARAM_SAFEEVAL
375 ,(SlvParameterInitBool){{"safeeval"
376 ,"Use safe evaluation?",1
377 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
378 "throw SIGFPE errors which will then halt integration."
379 }, FALSE}
380 );
381
382
383 slv_param_bool(p,IDA_PARAM_ATOLVECT
384 ,(SlvParameterInitBool){{"atolvect"
385 ,"Use 'ode_atol' values as specified?",1
386 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
387 " in the integration. If FALSE, a scalar absolute tolerance value"
388 " is shared by all variables. See IDA manual, section 5.5.1"
389 }, TRUE }
390 );
391
392 slv_param_real(p,IDA_PARAM_ATOL
393 ,(SlvParameterInitReal){{"atol"
394 ,"Scalar absolute error tolerance",1
395 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
396 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
397 }, 1e-5, 0, 1e10 }
398 );
399
400 slv_param_real(p,IDA_PARAM_RTOL
401 ,(SlvParameterInitReal){{"rtol"
402 ,"Scalar relative error tolerance",1
403 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
404 " it's not possible to set per-variable relative tolerances as it is"
405 " with LSODE)."
406 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
407 }, 1e-4, 0, 1 }
408 );
409
410 slv_param_char(p,IDA_PARAM_LINSOLVER
411 ,(SlvParameterInitChar){{"linsolver"
412 ,"Linear solver",1
413 ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"
414 " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"
415 " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"
416 " or SPTFQMR (which still need preconditioners to be implemented"
417 " before they can be very useful."
418 }, "DENSE"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
419 );
420
421 slv_param_int(p,IDA_PARAM_MAXL
422 ,(SlvParameterInitInt){{"maxl"
423 ,"Maximum Krylov dimension",0
424 ,"The maximum dimension of Krylov space used by the linear solver"
425 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
426 " The default of 0 results in IDA using its internal default, which"
427 " is currently a value of 5."
428 }, 0, 0, 20 }
429 );
430
431 slv_param_int(p,IDA_PARAM_MAXORD
432 ,(SlvParameterInitInt){{"maxord"
433 ,"Maximum order of linear multistep method",0
434 ,"The maximum order of the linear multistep method with IDA. See"
435 " IDA manual p 38."
436 }, 5, 1, 5 }
437 );
438
439 slv_param_bool(p,IDA_PARAM_GSMODIFIED
440 ,(SlvParameterInitBool){{"gsmodified"
441 ,"Gram-Schmidt Orthogonalisation Scheme", 2
442 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
443 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
444 }, TRUE}
445 );
446
447 slv_param_int(p,IDA_PARAM_MAXNCF
448 ,(SlvParameterInitInt){{"maxncf"
449 ,"Max nonlinear solver convergence failures per step", 2
450 ,"Maximum number of allowable nonlinear solver convergence failures"
451 " on one step. See IDA manual section 5.5.6.1."
452 }, 10,0,1000 }
453 );
454
455 slv_param_char(p,IDA_PARAM_PREC
456 ,(SlvParameterInitChar){{"prec"
457 ,"Preconditioner",1
458 ,"See IDA manual, section section 5.6.8."
459 },"NONE"}, (char *[]){"NONE","DIAG",NULL}
460 );
461
462 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
463
464 CONSOLE_DEBUG("Created %d params", p->num_parms);
465
466 return 0;
467 }
468
469 /*-------------------------------------------------------------
470 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
471 */
472
473 /*static double div1(double a, double b){
474 return a/b;
475 }*/
476
477 typedef int IdaFlagFn(void *,int *);
478 typedef char *IdaFlagNameFn(int);
479
480 /* return 0 on success */
481 int integrator_ida_solve(
482 IntegratorSystem *sys
483 , unsigned long start_index
484 , unsigned long finish_index
485 ){
486 void *ida_mem;
487 int flag, flag1, t_index;
488 realtype t0, reltol, abstol, t, tret, tout1;
489 N_Vector y0, yp0, abstolvect, ypret, yret, id;
490 IntegratorIdaData *enginedata;
491 char *linsolver;
492 int maxl;
493 IdaFlagFn *flagfn;
494 IdaFlagNameFn *flagnamefn;
495 const char *flagfntype;
496 char *pname = NULL;
497 struct rel_relation **rels;
498 int *rootsfound;
499
500 #ifdef SOLVE_DEBUG
501 char *varname, *relname;
502 #endif
503
504 int i,j,n_activerels,n_solverrels;
505 const IntegratorIdaPrec *prec = NULL;
506 int icopt; /* initial conditions strategy */
507
508 CONSOLE_DEBUG("STARTING IDA...");
509
510 enginedata = integrator_ida_enginedata(sys);
511
512 enginedata->safeeval = SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_SAFEEVAL);
513 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
514
515 /* store reference to list of relations (in enginedata) */
516 n_solverrels = slv_get_num_solvers_rels(sys->system);
517
518 n_activerels = slv_count_solvers_rels(sys->system, &integrator_ida_rel);
519
520 enginedata->bndlist = slv_get_solvers_bnd_list(sys->system);
521 enginedata->nbnds = slv_get_num_solvers_bnds(sys->system);
522
523 enginedata->rellist = ASC_NEW_ARRAY(struct rel_relation *, n_activerels);
524
525 rels = slv_get_solvers_rel_list(sys->system);
526
527 j=0;
528 for(i=0; i < n_solverrels; ++i){
529 if(rel_apply_filter(rels[i], &integrator_ida_rel)){
530 #ifdef SOLVE_DEBUG
531 relname = rel_make_name(sys->system, rels[i]);
532 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));
533 ASC_FREE(relname);
534 #endif
535 enginedata->rellist[j++]=rels[i];
536 }
537 }
538 asc_assert(j==n_activerels);
539
540 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
541 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
542
543 CONSOLE_DEBUG("Number of relations: %d",n_solverrels);
544 CONSOLE_DEBUG("Number of active relations: %d",n_activerels);
545 CONSOLE_DEBUG("Number of dependent vars: %d",sys->n_y);
546 CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);
547
548 enginedata->nrels = n_activerels;
549
550 if(enginedata->nrels != sys->n_y){
551 ERROR_REPORTER_HERE(ASC_USER_ERROR
552 ,"Integration problem is not square (%d active rels, %d vars)"
553 ,n_activerels, sys->n_y
554 );
555 return 1; /* failure */
556 }
557
558 #ifdef SOLVE_DEBUG
559 integrator_ida_debug(sys,stderr);
560 #endif
561
562 /* retrieve initial values from the system */
563
564 /** @TODO fix this, the starting time != first sample */
565 t0 = integrator_get_t(sys);
566 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
567
568 CONSOLE_DEBUG("RETRIEVING y0");
569
570 y0 = N_VNew_Serial(sys->n_y);
571 integrator_get_y(sys,NV_DATA_S(y0));
572
573 #ifdef SOLVE_DEBUG
574 CONSOLE_DEBUG("RETRIEVING yp0");
575 #endif
576
577 yp0 = N_VNew_Serial(sys->n_y);
578 integrator_get_ydot(sys,NV_DATA_S(yp0));
579
580 #ifdef SOLVE_DEBUG
581 N_VPrint_Serial(yp0);
582 CONSOLE_DEBUG("yp0 is at %p",&yp0);
583 #endif
584
585 /* create IDA object */
586 ida_mem = IDACreate();
587
588 /* relative error tolerance */
589 reltol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_RTOL);
590 CONSOLE_DEBUG("rtol = %8.2e",reltol);
591
592 /* allocate internal memory */
593 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
594 /* vector of absolute tolerances */
595 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
596 abstolvect = N_VNew_Serial(sys->n_y);
597 integrator_get_atol(sys,NV_DATA_S(abstolvect));
598
599 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
600
601 N_VDestroy_Serial(abstolvect);
602 }else{
603 /* scalar absolute tolerance (one value for all) */
604 abstol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_ATOL);
605 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
606 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
607 }
608
609 if(flag==IDA_MEM_NULL){
610 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
611 return 2;
612 }else if(flag==IDA_MEM_FAIL){
613 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
614 return 3;
615 }else if(flag==IDA_ILL_INPUT){
616 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
617 return 4;
618 }/* else success */
619
620 /* set optional inputs... */
621 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)sys);
622 IDASetRdata(ida_mem, (void *)sys);
623 IDASetMaxStep(ida_mem, integrator_get_maxstep(sys));
624 IDASetInitStep(ida_mem, integrator_get_stepzero(sys));
625 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(sys));
626 if(integrator_get_minstep(sys)>0){
627 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
628 }
629
630 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXNCF));
631 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXNCF));
632
633 CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXORD));
634 IDASetMaxOrd(ida_mem,SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXORD));
635
636 /* there's no capability for setting *minimum* step size in IDA */
637
638
639 /* attach linear solver module, using the default value of maxl */
640 linsolver = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_LINSOLVER);
641 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
642 if(strcmp(linsolver,"ASCEND")==0){
643 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",sys->n_y);
644 IDAASCEND(ida_mem,sys->n_y);
645 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)sys);
646
647 flagfntype = "IDAASCEND";
648 flagfn = &IDAASCENDGetLastFlag;
649 flagnamefn = &IDAASCENDGetReturnFlagName;
650
651 }else if(strcmp(linsolver,"DENSE")==0){
652 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",sys->n_y);
653 flag = IDADense(ida_mem, sys->n_y);
654 switch(flag){
655 case IDADENSE_SUCCESS: break;
656 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 5;
657 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 5;
658 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 5;
659 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
660 }
661
662 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
663 CONSOLE_DEBUG("USING AUTODIFF");
664 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)sys);
665 switch(flag){
666 case IDADENSE_SUCCESS: break;
667 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
668 }
669 }else{
670 CONSOLE_DEBUG("USING NUMERICAL DIFF");
671 }
672
673 flagfntype = "IDADENSE";
674 flagfn = &IDADenseGetLastFlag;
675 flagnamefn = &IDADenseGetReturnFlagName;
676 }else{
677 /* remaining methods are all SPILS */
678 CONSOLE_DEBUG("IDA SPILS");
679
680 maxl = SLV_PARAM_INT(&(sys->params),IDA_PARAM_MAXL);
681 CONSOLE_DEBUG("maxl = %d",maxl);
682
683 /* what preconditioner? */
684 pname = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_PREC);
685 if(strcmp(pname,"NONE")==0){
686 prec = NULL;
687 }else if(strcmp(pname,"JACOBI")==0){
688 prec = &prec_jacobi;
689 }else{
690 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
691 return 7;
692 }
693
694 /* which SPILS linear solver? */
695 if(strcmp(linsolver,"SPGMR")==0){
696 CONSOLE_DEBUG("IDA SPGMR");
697 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
698 }else if(strcmp(linsolver,"SPBCG")==0){
699 CONSOLE_DEBUG("IDA SPBCG");
700 flag = IDASpbcg(ida_mem, maxl);
701 }else if(strcmp(linsolver,"SPTFQMR")==0){
702 CONSOLE_DEBUG("IDA SPTFQMR");
703 flag = IDASptfqmr(ida_mem,maxl);
704 }else{
705 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
706 return 8;
707 }
708
709 if(prec){
710 /* assign the preconditioner to the linear solver */
711 (prec->pcreate)(sys);
712 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)sys);
713 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
714 }else{
715 CONSOLE_DEBUG("No preconditioner");
716 }
717
718 flagfntype = "IDASPILS";
719 flagfn = &IDASpilsGetLastFlag;
720 flagnamefn = &IDASpilsGetReturnFlagName;
721
722 if(flag==IDASPILS_MEM_NULL){
723 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
724 return 9;
725 }else if(flag==IDASPILS_MEM_FAIL){
726 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
727 return 9;
728 }/* else success */
729
730 /* assign the J*v function */
731 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
732 CONSOLE_DEBUG("USING AUTODIFF");
733 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)sys);
734 if(flag==IDASPILS_MEM_NULL){
735 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
736 return 10;
737 }else if(flag==IDASPILS_LMEM_NULL){
738 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
739 return 10;
740 }/* else success */
741 }else{
742 CONSOLE_DEBUG("USING NUMERICAL DIFF");
743 }
744
745 if(strcmp(linsolver,"SPGMR")==0){
746 /* select Gram-Schmidt orthogonalisation */
747 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_GSMODIFIED)){
748 CONSOLE_DEBUG("USING MODIFIED GS");
749 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
750 if(flag!=IDASPILS_SUCCESS){
751 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
752 return 11;
753 }
754 }else{
755 CONSOLE_DEBUG("USING CLASSICAL GS");
756 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
757 if(flag!=IDASPILS_SUCCESS){
758 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
759 return 11;
760 }
761 }
762 }
763 }
764
765 /* set linear solver optional inputs...
766 ...nothing here at the moment...
767 */
768
769 /* calculate initial conditions */
770 icopt = 0;
771 if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"Y")==0){
772 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
773 icopt = IDA_Y_INIT;
774 asc_assert(icopt!=0);
775 }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
776 CONSOLE_DEBUG("Solving initial conditions using values of yd");
777 icopt = IDA_YA_YDP_INIT;
778 asc_assert(icopt!=0);
779 id = N_VNew_Serial(sys->n_y);
780 for(i=0; i < sys->n_y; ++i){
781 if(sys->ydot[i] == NULL){
782 NV_Ith_S(id,i) = 0.0;
783 #ifdef SOLVE_DEBUG
784 varname = var_make_name(sys->system,sys->y[i]);
785 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
786 ASC_FREE(varname);
787 #endif
788 }else{
789 #ifdef SOLVE_DEBUG
790 CONSOLE_DEBUG("y[%d] is differential",i);
791 #endif
792 NV_Ith_S(id,i) = 1.0;
793 }
794 }
795 IDASetId(ida_mem, id);
796 N_VDestroy_Serial(id);
797 }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"NONE")==0){
798 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
799 }else{
800 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
801 }
802
803 if(icopt){
804 sys->currentstep=0;
805 t_index=start_index + 1;
806 tout1 = samplelist_get(sys->samples, t_index);
807
808 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
809
810 #ifdef ASC_SIGNAL_TRAPS
811 /* catch SIGFPE if desired to */
812 if(enginedata->safeeval){
813 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
814 Asc_SignalHandlerPush(SIGFPE,SIG_DFL);
815 }else{
816 # ifdef FEX_DEBUG
817 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
818 # endif
819 Asc_SignalHandlerPushDefault(SIGFPE);
820 }
821 if (setjmp(g_fpe_env)==0) {
822 #endif
823
824 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
825 flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3 */
826 # else
827 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
828 # endif
829 /* check flags and output status */
830 switch(flag){
831 case IDA_SUCCESS:
832 CONSOLE_DEBUG("Initial conditions solved OK");
833 break;
834
835 case IDA_LSETUP_FAIL:
836 case IDA_LINIT_FAIL:
837 case IDA_LSOLVE_FAIL:
838 case IDA_NO_RECOVERY:
839 flag1 = -999;
840 flag = (flagfn)(ida_mem,&flag1);
841 if(flag){
842 ERROR_REPORTER_HERE(ASC_PROG_ERR
843 ,"Unable to retrieve error code from %s (err %d)"
844 ,flagfntype,flag
845 );
846 return 12;
847 }
848 ERROR_REPORTER_HERE(ASC_PROG_ERR
849 ,"%s returned flag '%s' (value = %d)"
850 ,flagfntype,(flagnamefn)(flag1),flag1
851 );
852 return 12;
853
854 default:
855 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
856 return 12;
857 }
858 #ifdef ASC_SIGNAL_TRAPS
859 }else{
860 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
861 return 13;
862 }
863
864 if(enginedata->safeeval){
865 Asc_SignalHandlerPop(SIGFPE,SIG_DFL);
866 }else{
867 CONSOLE_DEBUG("pop...");
868 Asc_SignalHandlerPopDefault(SIGFPE);
869 CONSOLE_DEBUG("...pop");
870 }
871 #endif
872 }/* icopt */
873
874 /* optionally, specify ROO-FINDING problem */
875 if(enginedata->nbnds){
876 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)sys);
877 }
878
879 /* -- set up the IntegratorReporter */
880 integrator_output_init(sys);
881
882 /* -- store the initial values of all the stuff */
883 integrator_output_write(sys);
884 integrator_output_write_obs(sys);
885
886 /* specify where the returned values should be stored */
887 yret = y0;
888 ypret = yp0;
889
890 /* advance solution in time, return values as yret and derivatives as ypret */
891 sys->currentstep=1;
892 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++sys->currentstep){
893 t = samplelist_get(sys->samples, t_index);
894 t0 = integrator_get_t(sys);
895 asc_assert(t > t0);
896
897 #ifdef SOLVE_DEBUG
898 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
899 #endif
900
901 #ifdef ASC_SIGNAL_TRAPS
902 Asc_SignalHandlerPushDefault(SIGINT);
903 if(setjmp(g_int_env)==0) {
904 #endif
905 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
906 #ifdef ASC_SIGNAL_TRAPS
907 }else{
908 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Caught interrupt");
909 flag = -555;
910 }
911 Asc_SignalHandlerPopDefault(SIGINT);
912 #endif
913
914 /* check for roots found */
915 if(enginedata->nbnds){
916 rootsfound = ASC_NEW_ARRAY(int,enginedata->nbnds);
917 if(IDA_SUCCESS == IDAGetRootInfo(ida_mem, rootsfound)){
918 for(i=0; i < enginedata->nbnds; ++i){
919 if(rootsfound[i]){
920 #ifdef SOLVE_DEBUG
921 relname = bnd_make_name(sys->system,enginedata->bndlist[i]);
922 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary '%s' crossed",relname);
923 ASC_FREE(relname);
924 #else
925 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary crossed!");
926 #endif
927 }
928 }
929 }else{
930 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");
931 }
932 ASC_FREE(rootsfound);
933 /* so, now we need to restart the integration. we will assume that
934 everything changes: number of variables, etc, etc, etc. */
935 }
936
937
938 /* pass the values of everything back to the compiler */
939 integrator_set_t(sys, (double)tret);
940 integrator_set_y(sys, NV_DATA_S(yret));
941 integrator_set_ydot(sys, NV_DATA_S(ypret));
942
943 if(flag<0){
944 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
945 break;
946 }
947
948 /* -- do something so that sys knows the values of tret, yret and ypret */
949
950 /* -- store the current values of all the stuff */
951 integrator_output_write(sys);
952 integrator_output_write_obs(sys);
953
954 }/* loop through next sample timestep */
955
956 /* -- close the IntegratorReporter */
957 integrator_output_close(sys);
958
959 /* get optional outputs */
960 #ifdef STATS_DEBUG
961 IntegratorIdaStats stats;
962 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
963 integrator_ida_write_stats(&stats);
964 }else{
965 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
966 }
967 #endif
968
969 /* free solution memory */
970 N_VDestroy_Serial(yret);
971 N_VDestroy_Serial(ypret);
972
973 /* free solver memory */
974 IDAFree(ida_mem);
975
976 if(flag < -500){
977 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", t);
978 return -flag;
979 }
980
981 if(flag < 0){
982 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
983 return 14;
984 }
985
986 /* all done, success */
987 return 0;
988 }
989
990 /*--------------------------------------------------
991 RESIDUALS AND JACOBIAN AND IDAROOTFN
992 */
993
994 #if 0
995 typedef void (SignalHandlerFn)(int);
996 SignalHandlerFn integrator_ida_sig;
997 SignalHandlerFn *integrator_ida_sig_old;
998 jmp_buf integrator_ida_jmp_buf;
999 fenv_t integrator_ida_fenv_old;
1000
1001
1002 void integrator_ida_write_feinfo(){
1003 int f;
1004 f = fegetexcept();
1005 CONSOLE_DEBUG("Locating nature of exception...");
1006 if(f & FE_DIVBYZERO)ERROR_REPORTER_HERE(ASC_PROG_ERR,"DIV BY ZERO");
1007 if(f & FE_INEXACT)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INEXACT");
1008 if(f & FE_INVALID)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INVALID");
1009 if(f & FE_OVERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"OVERFLOW");
1010 if(f & FE_UNDERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"UNDERFLOW");
1011 if(f==0)ERROR_REPORTER_HERE(ASC_PROG_ERR,"FLAGS ARE CLEAR?!?");
1012 }
1013
1014 void integrator_ida_sig(int sig){
1015 /* the wrong signal: rethrow to the default handler */
1016 if(sig!=SIGFPE){
1017 signal(SIGFPE,SIG_DFL);
1018 raise(sig);
1019 }
1020 integrator_ida_write_feinfo();
1021 CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);
1022 longjmp(integrator_ida_jmp_buf,sig);
1023 }
1024 #endif
1025
1026 /**
1027 Function to evaluate system residuals, in the form required for IDA.
1028
1029 Given tt, yy and yp, we need to evaluate and return rr.
1030
1031 @param tt current value of indep variable (time)
1032 @param yy current values of dependent variable vector
1033 @param yp current values of derivatives of dependent variables
1034 @param rr the output residual vector (is we're returning data to)
1035 @param res_data pointer to our stuff (sys in this case).
1036
1037 @return 0 on success, positive on recoverable error, and
1038 negative on unrecoverable error.
1039 */
1040 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1041 IntegratorSystem *sys;
1042 IntegratorIdaData *enginedata;
1043 int i, calc_ok, is_error;
1044 struct rel_relation** relptr;
1045 double resid;
1046 char *relname;
1047 #ifdef FEX_DEBUG
1048 char *varname;
1049 char diffname[30];
1050 #endif
1051
1052 sys = (IntegratorSystem *)res_data;
1053 enginedata = integrator_ida_enginedata(sys);
1054
1055 #ifdef FEX_DEBUG
1056 /* fprintf(stderr,"\n\n"); */
1057 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1058 #endif
1059
1060 if(NV_LENGTH_S(rr)!=enginedata->nrels){
1061 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
1062 return -1; /* unrecoverable */
1063 }
1064
1065 /* pass the values of everything back to the compiler */
1066 integrator_set_t(sys, (double)tt);
1067 integrator_set_y(sys, NV_DATA_S(yy));
1068 integrator_set_ydot(sys, NV_DATA_S(yp));
1069
1070 /* perform bounds checking on all variables */
1071 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1072 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1073 return 1;
1074 }
1075
1076 /* evaluate each residual in the rellist */
1077 is_error = 0;
1078 relptr = enginedata->rellist;
1079
1080
1081 #ifdef ASC_SIGNAL_TRAPS
1082 if(enginedata->safeeval){
1083 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1084 }else{
1085 # ifdef FEX_DEBUG
1086 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1087 # endif
1088 Asc_SignalHandlerPushDefault(SIGFPE);
1089 }
1090
1091 if (SETJMP(g_fpe_env)==0) {
1092 #endif
1093
1094
1095 for(i=0, relptr = enginedata->rellist;
1096 i< enginedata->nrels && relptr != NULL;
1097 ++i, ++relptr
1098 ){
1099 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
1100
1101 NV_Ith_S(rr,i) = resid;
1102 if(!calc_ok){
1103 relname = rel_make_name(sys->system, *relptr);
1104 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1105 ASC_FREE(relname);
1106 /* presumable some output already made? */
1107 is_error = 1;
1108 }/*else{
1109 CONSOLE_DEBUG("Calc OK");
1110 }*/
1111 }
1112
1113 if(!is_error){
1114 for(i=0;i< enginedata->nrels; ++i){
1115 if(isnan(NV_Ith_S(rr,i))){
1116 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in residual %d",i);
1117 is_error=1;
1118 }
1119 }
1120 #ifdef FEX_DEBUG
1121 if(!is_error){
1122 CONSOLE_DEBUG("No NAN detected");
1123 }
1124 #endif
1125 }
1126
1127 #ifdef ASC_SIGNAL_TRAPS
1128 }else{
1129 relname = rel_make_name(sys->system, *relptr);
1130 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1131 ASC_FREE(relname);
1132 is_error = 1;
1133 }
1134
1135 if(enginedata->safeeval){
1136 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1137 }else{
1138 Asc_SignalHandlerPopDefault(SIGFPE);
1139 }
1140 #endif
1141
1142
1143 #ifdef FEX_DEBUG
1144 /* output residuals to console */
1145 CONSOLE_DEBUG("RESIDUAL OUTPUT");
1146 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1147 for(i=0; i<sys->n_y; ++i){
1148 varname = var_make_name(sys->system,sys->y[i]);
1149 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1150 if(sys->ydot[i]){
1151 varname = var_make_name(sys->system,sys->ydot[i]);
1152 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1153 }else{
1154 snprintf(diffname,99,"diff(%s)",varname);
1155 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1156 }
1157 ASC_FREE(varname);
1158 relname = rel_make_name(sys->system,enginedata->rellist[i]);
1159 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1160 }
1161 #endif
1162
1163 if(is_error){
1164 return 1;
1165 }
1166
1167 #ifdef FEX_DEBUG
1168 CONSOLE_DEBUG("RESIDUAL OK");
1169 #endif
1170 return 0;
1171 }
1172
1173 /**
1174 Dense Jacobian evaluation. Only suitable for small problems!
1175 Has been seen working for problems up to around 2000 vars, FWIW.
1176 */
1177 int integrator_ida_djex(long int Neq, realtype tt
1178 , N_Vector yy, N_Vector yp, N_Vector rr
1179 , realtype c_j, void *jac_data, DenseMat Jac
1180 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1181 ){
1182 IntegratorSystem *sys;
1183 IntegratorIdaData *enginedata;
1184 char *relname;
1185 #ifdef DJEX_DEBUG
1186 struct var_variable **varlist;
1187 char *varname;
1188 #endif
1189 struct rel_relation **relptr;
1190 int i;
1191 double *derivatives;
1192 struct var_variable **variables;
1193 int count, j;
1194 int status, is_error = 0;
1195
1196 sys = (IntegratorSystem *)jac_data;
1197 enginedata = integrator_ida_enginedata(sys);
1198
1199 /* allocate space for returns from relman_diff3 */
1200 /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
1201 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1202 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1203
1204 /* pass the values of everything back to the compiler */
1205 integrator_set_t(sys, (double)tt);
1206 integrator_set_y(sys, NV_DATA_S(yy));
1207 integrator_set_ydot(sys, NV_DATA_S(yp));
1208
1209 /* perform bounds checking on all variables */
1210 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1211 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1212 return 1;
1213 }
1214
1215 #ifdef DJEX_DEBUG
1216 varlist = slv_get_solvers_var_list(sys->system);
1217
1218 /* print vars */
1219 for(i=0; i < sys->n_y; ++i){
1220 varname = var_make_name(sys->system, sys->y[i]);
1221 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1222 asc_assert(NV_Ith_S(yy,i) == var_value(sys->y[i]));
1223 ASC_FREE(varname);
1224 }
1225
1226 /* print derivatives */
1227 for(i=0; i < sys->n_y; ++i){
1228 if(sys->ydot[i]){
1229 varname = var_make_name(sys->system, sys->ydot[i]);
1230 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(sys->ydot[i]));
1231 ASC_FREE(varname);
1232 }else{
1233 varname = var_make_name(sys->system, sys->y[i]);
1234 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1235 ASC_FREE(varname);
1236 }
1237 }
1238
1239 /* print step size */
1240 CONSOLE_DEBUG("<c_j> = %g",c_j);
1241 #endif
1242
1243 /* build up the dense jacobian matrix... */
1244 status = 0;
1245 for(i=0, relptr = enginedata->rellist;
1246 i< enginedata->nrels && relptr != NULL;
1247 ++i, ++relptr
1248 ){
1249 /* get derivatives for this particular relation */
1250 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1251
1252 if(status){
1253 relname = rel_make_name(sys->system, *relptr);
1254 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1255 ASC_FREE(relname);
1256 is_error = 1;
1257 break;
1258 }
1259
1260 /* output what's going on here ... */
1261 #ifdef DJEX_DEBUG
1262 relname = rel_make_name(sys->system, *relptr);
1263 fprintf(stderr,"%d: '%s': ",i,relname);
1264 for(j=0;j<count;++j){
1265 varname = var_make_name(sys->system, variables[j]);
1266 if(var_deriv(variables[j])){
1267 fprintf(stderr," '%s'=",varname);
1268 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(sys,variables[j]));
1269 }else{
1270 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1271 }
1272 ASC_FREE(varname);
1273 }
1274 /* relname is freed further down */
1275 fprintf(stderr,"\n");
1276 #endif
1277
1278 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1279 for(j=0; j < count; ++j){
1280 #ifdef DJEX_DEBUG
1281 varname = var_make_name(sys->system,variables[j]);
1282 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1283 ASC_FREE(varname);
1284 #endif
1285 if(!var_deriv(variables[j])){
1286 #ifdef DJEX_DEBUG
1287 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
1288 asc_assert(var_sindex(variables[j]) >= 0);
1289 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1290 #endif
1291 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1292 }else{
1293 DENSE_ELEM(Jac,i,integrator_ida_diffindex(sys,variables[j])) += derivatives[j] * c_j;
1294 #ifdef DJEX_DEBUG
1295 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1296 #endif
1297 }
1298 }
1299 }
1300
1301 #ifdef DJEX_DEBUG
1302 ASC_FREE(relname);
1303 CONSOLE_DEBUG("PRINTING JAC");
1304 fprintf(stderr,"\t");
1305 for(j=0; j < sys->n_y; ++j){
1306 if(j)fprintf(stderr,"\t");
1307 varname = var_make_name(sys->system,sys->y[j]);
1308 fprintf(stderr,"%11s",varname);
1309 ASC_FREE(varname);
1310 }
1311 fprintf(stderr,"\n");
1312 for(i=0; i < enginedata->nrels; ++i){
1313 relname = rel_make_name(sys->system, enginedata->rellist[i]);
1314 fprintf(stderr,"%s\t",relname);
1315 ASC_FREE(relname);
1316
1317 for(j=0; j < sys->n_y; ++j){
1318 if(j)fprintf(stderr,"\t");
1319 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1320 }
1321 fprintf(stderr,"\n");
1322 }
1323 #endif
1324
1325 /* test for NANs */
1326 if(!is_error){
1327 for(i=0;i< enginedata->nrels; ++i){
1328 for(j=0;j<sys->n_y;++j){
1329 if(isnan(DENSE_ELEM(Jac,i,j))){
1330 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
1331 is_error=1;
1332 }
1333 }
1334 }
1335 #ifdef DJEX_DEBUG
1336 if(!is_error){
1337 CONSOLE_DEBUG("No NAN detected");
1338 }
1339 #endif
1340 }
1341
1342 /* if(integrator_ida_check_diffindex(sys)){
1343 is_error = 1;
1344 }*/
1345
1346 if(is_error){
1347 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1348 return 1;
1349 }
1350
1351 #ifdef DJEX_DEBUG
1352 CONSOLE_DEBUG("DJEX RETURNING 0");
1353 /* ASC_PANIC("Quitting"); */
1354 #endif
1355 return 0;
1356 }
1357
1358 /**
1359 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1360
1361 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1362
1363 @param tt current value of the independent variable (time, t)
1364 @param yy current value of the dependent variable vector, y(t).
1365 @param yp current value of y'(t).
1366 @param rr current value of the residual vector F(t, y, y').
1367 @param v the vector by which the Jacobian must be multiplied to the right.
1368 @param Jv the output vector computed
1369 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1370 @param jac_data pointer to our stuff (sys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1371 @param tmp1 @see tmp2
1372 @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.
1373 @return 0 on success
1374 */
1375 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1376 , N_Vector v, N_Vector Jv, realtype c_j
1377 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1378 ){
1379 IntegratorSystem *sys;
1380 IntegratorIdaData *enginedata;
1381 int i, j, is_error=0;
1382 struct rel_relation** relptr;
1383 char *relname;
1384 int status;
1385 double Jv_i;
1386
1387 struct var_variable **variables;
1388 double *derivatives;
1389 int count;
1390 struct var_variable **varlist;
1391 #ifdef JEX_DEBUG
1392
1393 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1394 #endif
1395
1396 sys = (IntegratorSystem *)jac_data;
1397 enginedata = integrator_ida_enginedata(sys);
1398 varlist = slv_get_solvers_var_list(sys->system);
1399
1400 /* pass the values of everything back to the compiler */
1401 integrator_set_t(sys, (double)tt);
1402 integrator_set_y(sys, NV_DATA_S(yy));
1403 integrator_set_ydot(sys, NV_DATA_S(yp));
1404 /* no real use for residuals (rr) here, I don't think? */
1405
1406 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1407
1408 i = NV_LENGTH_S(yy) * 2;
1409 #ifdef JEX_DEBUG
1410 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1411 #endif
1412 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1413 derivatives = ASC_NEW_ARRAY(double, i);
1414
1415 /* evaluate the derivatives... */
1416 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1417
1418 #ifdef ASC_SIGNAL_TRAPS
1419 Asc_SignalHandlerPushDefault(SIGFPE);
1420 if (SETJMP(g_fpe_env)==0) {
1421 #endif
1422 for(i=0, relptr = enginedata->rellist;
1423 i< enginedata->nrels && relptr != NULL;
1424 ++i, ++relptr
1425 ){
1426 /* get derivatives for this particular relation */
1427 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1428 #ifdef JEX_DEBUG
1429 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1430 #endif
1431
1432 if(status){
1433 relname = rel_make_name(sys->system, *relptr);
1434 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1435 ASC_FREE(relname);
1436 is_error = 1;
1437 break;
1438 }
1439
1440 /*
1441 Now we have the derivatives wrt each alg/diff variable in the
1442 present equation. variables[] points into the varlist. need
1443 a mapping from the varlist to the y and ydot lists.
1444 */
1445
1446 Jv_i = 0;
1447 for(j=0; j < count; ++j){
1448 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], sys->n_y);
1449 varname = var_make_name(sys->system, enginedata->varlist[variables[j]]);
1450 if(varname){
1451 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1452 ASC_FREE(varname);
1453 }else{
1454 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1455 }
1456 */
1457
1458 /* we don't calculate derivatives wrt indep var */
1459 asc_assert(variables[j]>=0);
1460 if(variables[j] == sys->x) continue;
1461 #ifdef JEX_DEBUG
1462 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1463 #endif
1464 if(var_deriv(variables[j])){
1465 #define DIFFINDEX integrator_ida_diffindex(sys,variables[j])
1466 #ifdef JEX_DEBUG
1467 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
1468 , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
1469 , i, DIFFINDEX, derivatives[j]
1470 , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
1471 );
1472 #endif
1473 asc_assert(sys->ydot[DIFFINDEX]==variables[j]);
1474 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1475 #undef DIFFINDEX
1476 }else{
1477 #define VARINDEX var_sindex(variables[j])
1478 #ifdef JEX_DEBUG
1479 asc_assert(sys->y[VARINDEX]==variables[j]);
1480 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
1481 , i
1482 , derivatives[j] * NV_Ith_S(v,VARINDEX)
1483 , i, VARINDEX, derivatives[j]
1484 , VARINDEX, NV_Ith_S(v,VARINDEX)
1485 );
1486 #endif
1487 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1488 #undef VARINDEX
1489 }
1490 }
1491
1492 NV_Ith_S(Jv,i) = Jv_i;
1493 #ifdef JEX_DEBUG
1494 CONSOLE_DEBUG("rel = %p",*relptr);
1495 relname = rel_make_name(sys->system, *relptr);
1496 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1497 ASC_FREE(relname);
1498 return 1;
1499 #endif
1500 }
1501 #ifdef ASC_SIGNAL_TRAPS
1502 }else{
1503 relname = rel_make_name(sys->system, *relptr);
1504 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1505 ASC_FREE(relname);
1506 is_error = 1;
1507 }
1508 Asc_SignalHandlerPopDefault(SIGFPE);
1509 #endif
1510
1511 if(is_error){
1512 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1513 return 1;
1514 }
1515 return 0;
1516 }
1517
1518 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1519 int integrator_ida_sjex(long int Neq, realtype tt
1520 , N_Vector yy, N_Vector yp, N_Vector rr
1521 , realtype c_j, void *jac_data, mtx_matrix_t Jac
1522 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1523 ){
1524 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1525 return -1;
1526 }
1527
1528 /* root finding function */
1529
1530 int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data){
1531 IntegratorSystem *sys;
1532 IntegratorIdaData *enginedata;
1533 int i;
1534
1535 asc_assert(g_data!=NULL);
1536 sys = (IntegratorSystem *)g_data;
1537 enginedata = integrator_ida_enginedata(sys);
1538
1539 /* pass the values of everything back to the compiler */
1540 integrator_set_t(sys, (double)tt);
1541 integrator_set_y(sys, NV_DATA_S(yy));
1542 integrator_set_ydot(sys, NV_DATA_S(yp));
1543
1544 asc_assert(gout!=NULL);
1545
1546 #ifdef ROOT_DEBUG
1547 CONSOLE_DEBUG("t = %f",tt);
1548 #endif
1549
1550 /* evaluate the residuals for each of the boundaries */
1551 for(i=0; i < enginedata->nbnds; ++i){
1552 switch(bnd_kind(enginedata->bndlist[i])){
1553 case e_bnd_rel: /* real-valued boundary relation */
1554 gout[i] = bndman_real_eval(enginedata->bndlist[i]);
1555 #ifdef ROOT_DEBUG
1556 CONSOLE_DEBUG("gout[%d] = %f",i,gout[i]);
1557 #endif
1558 break;
1559 case e_bnd_logrel:
1560 if(bndman_log_eval(enginedata->bndlist[i])){
1561 CONSOLE_DEBUG("bnd[%d] = TRUE",i);
1562 gout[i] = +2.0;
1563 }else{
1564 CONSOLE_DEBUG("bnd[%d] = FALSE",i);
1565 gout[i] = -2.0;
1566 }
1567 break;
1568 }
1569 }
1570
1571 return 0; /* no way to detect errors in bndman_*_eval at this stage */
1572 }
1573
1574 /*----------------------------------------------
1575 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1576 */
1577
1578 void integrator_ida_pcreate_jacobi(IntegratorSystem *sys){
1579 IntegratorIdaData *enginedata =sys->enginedata;
1580 IntegratorIdaPrecDataJacobi *precdata;
1581 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1582
1583 asc_assert(sys->n_y);
1584 precdata->PIii = N_VNew_Serial(sys->n_y);
1585
1586 enginedata->pfree = &integrator_ida_pfree_jacobi;
1587 enginedata->precdata = precdata;
1588 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1589 }
1590
1591 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1592 if(enginedata->precdata){
1593 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1594 N_VDestroy_Serial(precdata->PIii);
1595
1596 ASC_FREE(precdata);
1597 enginedata->precdata = NULL;
1598 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1599 }
1600 enginedata->pfree = NULL;
1601 }
1602
1603 /**
1604 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1605
1606 'setup' function.
1607 */
1608 int integrator_ida_psetup_jacobi(realtype tt,
1609 N_Vector yy, N_Vector yp, N_Vector rr,
1610 realtype c_j, void *p_data,
1611 N_Vector tmp1, N_Vector tmp2,
1612 N_Vector tmp3
1613 ){
1614 int i, j, res;
1615 IntegratorSystem *sys;
1616 IntegratorIdaData *enginedata;
1617 IntegratorIdaPrecDataJacobi *precdata;
1618 struct rel_relation **relptr;
1619
1620 sys = (IntegratorSystem *)p_data;
1621 enginedata = sys->enginedata;
1622 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1623 double *derivatives;
1624 struct var_variable **variables;
1625 int count, status;
1626 char *relname;
1627
1628 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1629
1630 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1631 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1632
1633 /**
1634 @TODO FIXME here we are using the very inefficient and contorted approach
1635 of calculating the whole jacobian, then extracting just the diagonal elements.
1636 */
1637
1638 for(i=0, relptr = enginedata->rellist;
1639 i< enginedata->nrels && relptr != NULL;
1640 ++i, ++relptr
1641 ){
1642
1643 /* get derivatives for this particular relation */
1644 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1645 if(status){
1646 relname = rel_make_name(sys->system, *relptr);
1647 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1648 ASC_FREE(relname);
1649 break;
1650 }
1651 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1652 /* find the diagonal elements */
1653 for(j=0; j<count; ++j){
1654 if(var_sindex(variables[j])==i){
1655 if(var_deriv(variables[j])){
1656 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
1657 }else{
1658 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1659 }
1660
1661 }
1662 }
1663 #ifdef PREC_DEBUG
1664 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1665 #endif
1666 }
1667
1668 if(status){
1669 CONSOLE_DEBUG("Error found when evaluating derivatives");
1670 res = 1; goto finish; /* recoverable */
1671 }
1672
1673 integrator_ida_write_incidence(sys);
1674
1675 res = 0;
1676 finish:
1677 ASC_FREE(variables);
1678 ASC_FREE(derivatives);
1679 return res;
1680 };
1681
1682 /**
1683 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1684
1685 'solve' function.
1686 */
1687 int integrator_ida_psolve_jacobi(realtype tt,
1688 N_Vector yy, N_Vector yp, N_Vector rr,
1689 N_Vector rvec, N_Vector zvec,
1690 realtype c_j, realtype delta, void *p_data,
1691 N_Vector tmp
1692 ){
1693 IntegratorSystem *sys;
1694 IntegratorIdaData *data;
1695 IntegratorIdaPrecDataJacobi *precdata;
1696 sys = (IntegratorSystem *)p_data;
1697 data = sys->enginedata;
1698 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1699
1700 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1701 N_VProd(precdata->PIii, rvec, zvec);
1702 return 0;
1703 };
1704
1705 /*----------------------------------------------
1706 STATS
1707 */
1708
1709 /**
1710 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1711 status in a struct instead of separately.
1712
1713 @return IDA_SUCCESS on sucess.
1714 */
1715 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1716 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1717 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1718 ,&s->hlast, &s->hcur, &s->tcur
1719 );
1720 }
1721
1722 /**
1723 This routine just outputs the stats to the CONSOLE_DEBUG routine.
1724
1725 @TODO provide a GUI way of stats reporting from IDA.
1726 */
1727 void integrator_ida_write_stats(IntegratorIdaStats *stats){
1728 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
1729 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
1730 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
1731 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
1732 SI(qlast); SI(qcur);
1733 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
1734 # undef SL
1735 # undef SI
1736 # undef SR
1737 }
1738
1739 /*------------------------------------------------------------------------------
1740 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
1741 */
1742
1743 /**
1744 Here we construct the local transfer matrix. It's a bit of a naive
1745 approach; probably far more efficient ways can be worked out. But it will
1746 hopefully be a useful way to examine stability of some spatial
1747 discretisation schemes for PDAE systems.
1748
1749 http://ascendserver.cheme.cmu.edu/wiki/index.php/IDA#Stability
1750 */
1751 static int integrator_ida_transfer_matrix(const IntegratorSystem *sys, struct SystemJacobianStruct *J){
1752 int i=0, res;
1753 enum submat{II_GA=0, II_GD, II_FA, II_FD, II_FDP, II_NUM};
1754
1755 const var_filter_t *matvf[II_NUM] = {
1756 &system_vfilter_algeb
1757 ,&system_vfilter_diff
1758 ,&system_vfilter_algeb
1759 ,&system_vfilter_diff
1760 ,&system_vfilter_deriv
1761 };
1762
1763 const rel_filter_t *matrf[II_NUM] = {
1764 &system_rfilter_algeb
1765 ,&system_rfilter_algeb
1766 ,&system_rfilter_diff
1767 ,&system_rfilter_diff
1768 ,&system_rfilter_diff
1769 };
1770
1771 struct SystemJacobianStruct D[II_NUM];
1772
1773 for(i=0;i<II_NUM;++i){
1774 res = system_jacobian(sys->system, matrf[i], matvf[i], 1/*safe*/ ,&(D[i]));
1775 }
1776
1777 /* compute inverses for matrices that need it */
1778 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1779 return 1;
1780 }
1781
1782
1783 /**
1784 Our task here is to write the matrices that IDA *should* be seeing. We
1785 are actually making calls to relman_diff in order to do that, so we're
1786 really going back to the variables in the actualy system and computing
1787 row by row what the values are. This should mean just a single call to
1788 each blackbox present in the system (if blackbox caching is working
1789 correctly).
1790 */
1791 int integrator_ida_write_matrix(const IntegratorSystem *sys, FILE *f, const char *type){
1792 /* IntegratorIdaData *enginedata; */
1793 struct SystemJacobianStruct J = {NULL,NULL,NULL,0,0};
1794 int status=1;
1795 mtx_region_t R;
1796
1797 if(type==NULL)type = "dx'/dx";
1798
1799 if(0==strcmp(type,"dg/dz")){
1800 CONSOLE_DEBUG("Calculating dg/dz...");
1801 status = system_jacobian(sys->system
1802 , &system_rfilter_algeb, &system_vfilter_algeb
1803 , 1 /* safe */
1804 , &J
1805 );
1806 }else if(0==strcmp(type,"dg/dx")){
1807 CONSOLE_DEBUG("Calculating dg/dx...");
1808 status = system_jacobian(sys->system
1809 , &system_rfilter_algeb, &system_vfilter_diff
1810 , 1 /* safe */
1811 , &J
1812 );
1813 }else if(0==strcmp(type,"df/dx'")){
1814 CONSOLE_DEBUG("Calculating df/dx'...");
1815 status = system_jacobian(sys->system
1816 , &system_rfilter_diff, &system_vfilter_deriv
1817 , 1 /* safe */
1818 , &J
1819 );
1820 }else if(0==strcmp(type,"df/dz")){
1821 CONSOLE_DEBUG("Calculating df/dz...");
1822 status = system_jacobian(sys->system
1823 , &system_rfilter_diff, &system_vfilter_algeb
1824 , 1 /* safe */
1825 , &J
1826 );
1827 }else if(0==strcmp(type,"df/dx")){
1828 CONSOLE_DEBUG("Calculating df/dx...");
1829 status = system_jacobian(sys->system
1830 , &system_rfilter_diff, &system_vfilter_diff
1831 , 1 /* safe */
1832 , &J
1833 );
1834 }else if(0==strcmp(type,"dF/dy")){
1835 CONSOLE_DEBUG("Calculating dF/dy...");
1836 status = system_jacobian(sys->system
1837 , &system_rfilter_all, &system_vfilter_nonderiv
1838 , 1 /* safe */
1839 , &J
1840 );
1841 }else if(0==strcmp(type,"dF/dy'")){
1842 CONSOLE_DEBUG("Calculating dF/dy'...");
1843 status = system_jacobian(sys->system
1844 , &system_rfilter_all, &system_vfilter_deriv
1845 , 1 /* safe */
1846 , &J
1847 );
1848 }else if(0==strcmp(type,"dx'/dx")){
1849 /* system state transfer matrix dyd'/dyd */
1850 status = integrator_ida_transfer_matrix(sys, &J);
1851 }else{
1852 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid matrix type '%s'",type);
1853 return 1;
1854 }
1855
1856 if(status){
1857 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating matrix");
1858 }else{
1859 /* send the region explicitly, so that we handle non-square correctly */
1860 R.row.low = 0; R.col.low = 0;
1861 R.row.high = J.n_rels - 1; R.col.high = J.n_vars - 1;
1862 /* note that we're not fussy about empty matrices here... */
1863 mtx_write_region_mmio(f,J.M,&R);
1864 }
1865
1866 if(J.vars)ASC_FREE(J.vars);
1867 if(J.rels)ASC_FREE(J.rels);
1868 if(J.M)mtx_destroy(J.M);
1869
1870 return status;
1871 }
1872
1873 /**
1874 This routine outputs matrix structure in a crude text format, for the sake
1875 of debugging.
1876 */
1877 void integrator_ida_write_incidence(IntegratorSystem *sys){
1878 int i, j;
1879 struct rel_relation **relptr;
1880 IntegratorIdaData *enginedata = sys->enginedata;
1881 double *derivatives;
1882 struct var_variable **variables;
1883 int count, status;
1884 char *relname;
1885
1886 if(enginedata->nrels > 100){
1887 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
1888 return;
1889 }
1890
1891 variables = ASC_NEW_ARRAY(struct var_variable *, sys->n_y * 2);
1892 derivatives = ASC_NEW_ARRAY(double, sys->n_y * 2);
1893
1894 CONSOLE_DEBUG("Outputting incidence information to console...");
1895
1896 for(i=0, relptr = enginedata->rellist;
1897 i< enginedata->nrels && relptr != NULL;
1898 ++i, ++relptr
1899 ){
1900 relname = rel_make_name(sys->system, *relptr);
1901
1902 /* get derivatives for this particular relation */
1903 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1904 if(status){
1905 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1906 ASC_FREE(relname);
1907 break;
1908 }
1909
1910 fprintf(stderr,"%3d:%-15s:",i,relname);
1911 ASC_FREE(relname);
1912
1913 for(j=0; j<count; ++j){
1914 if(var_deriv(variables[j])){
1915 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(sys,variables[j]));
1916 }else{
1917 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
1918 }
1919 }
1920 fprintf(stderr,"\n");
1921 }
1922 ASC_FREE(variables);
1923 ASC_FREE(derivatives);
1924 }
1925
1926 /* @return 0 on success */
1927 int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
1928 char *varname, *relname;
1929 struct var_variable **vlist, *var;
1930 struct rel_relation **rlist, *rel;
1931 long vlen, rlen;
1932 long i;
1933 long di;
1934
1935 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",sys->n_y);
1936
1937 /* if(integrator_sort_obs_vars(sys))return 10; */
1938
1939 if(sys->y && sys->ydot){
1940 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
1941 fprintf(fp,"index\t%-15s\tydot\n","y");
1942 fprintf(fp,"-----\t%-15s\t-----\n","-----");
1943 for(i=0;i<sys->n_y;++i){
1944 varname = var_make_name(sys->system, sys->y[i]);
1945 fprintf(fp,"%ld\t%-15s\t",i,varname);
1946 if(sys->ydot[i]){
1947 ASC_FREE(varname);
1948 varname = var_make_name(sys->system, sys->ydot[i]);
1949 fprintf(fp,"%s\n",varname);
1950 ASC_FREE(varname);
1951 }else{
1952 fprintf(fp,".\n");
1953 ASC_FREE(varname);
1954 }
1955 }
1956 }else{
1957 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
1958 }
1959
1960 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
1961 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
1962 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
1963
1964
1965 /* visit all the slv_system_t master var lists to collect vars */
1966 /* find the vars mostly in this one */
1967 vlist = slv_get_solvers_var_list(sys->system);
1968 vlen = slv_get_num_solvers_vars(sys->system);
1969 for(i=0;i<vlen;i++){
1970 var = vlist[i];
1971
1972 varname = var_make_name(sys->system, var);
1973 fprintf(fp,"%ld\t%-15s\t",i,varname);
1974
1975 if(var_fixed(var)){
1976 // it's fixed, so not really a DAE var
1977 fprintf(fp,"(fixed)\n");
1978 }else if(!var_active(var)){
1979 // inactive
1980 fprintf(fp,"(inactive)\n");
1981 }else if(!var_incident(var)){
1982 // not incident
1983 fprintf(fp,"(not incident)\n");
1984 }else{
1985 if(var_deriv(var)){
1986 if(sys->y_id){
1987 di = integrator_ida_diffindex1(sys,var);
1988 if(di>=0){
1989 ASC_FREE(varname);
1990 varname = var_make_name(sys->system,vlist[di]);
1991 fprintf(fp,".\tdiff(%ld='%s')\n",di,varname);
1992 }else{
1993 fprintf(fp,".\tdiff(???,err=%ld)\n",di);
1994 }
1995 }else{
1996 fprintf(fp,".\tderiv... of??\n");
1997 }
1998 }else{
1999 fprintf(fp,"%d\t.\n",var_sindex(var));
2000 }
2001 }
2002 ASC_FREE(varname);
2003 }
2004
2005 /* let's write out the relations too */
2006 rlist = slv_get_solvers_rel_list(sys->system);
2007 rlen = slv_get_num_solvers_rels(sys->system);
2008
2009 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
2010 fprintf(fp,"index\tname\n");
2011 fprintf(fp,"-----\t----\n");
2012 for(i=0; i<rlen; ++i){
2013 rel = rlist[i];
2014 relname = rel_make_name(sys->system,rel);
2015 fprintf(fp,"%ld\t%s\n",i,relname);
2016 ASC_FREE(relname);
2017 }
2018
2019 /* write out the derivative chains */
2020 fprintf(fp,"\nDERIVATIVE CHAINS\n");
2021 if(integrator_ida_analyse_debug(sys,stderr)){
2022 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
2023 return 340;
2024 }
2025 fprintf(fp,"\n");
2026
2027 /* and lets write block debug output */
2028 system_block_debug(sys->system, fp);
2029
2030 return 0; /* success */
2031 }
2032
2033 /*----------------------------------------------
2034 ERROR REPORTING
2035 */
2036 /**
2037 Error message reporter function to be passed to IDA. All error messages
2038 will trigger a call to this function, so we should find everything
2039 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
2040 panel (in the case of PyGTK).
2041 */
2042 void integrator_ida_error(int error_code
2043 , const char *module, const char *function
2044 , char *msg, void *eh_data
2045 ){
2046 IntegratorSystem *sys;
2047 error_severity_t sev;
2048
2049 /* cast back the IntegratorSystem, just in case we need it */
2050 sys = (IntegratorSystem *)eh_data;
2051
2052 /* severity depends on the sign of the error_code value */
2053 if(error_code <= 0){
2054 sev = ASC_PROG_ERR;
2055 }else{
2056 sev = ASC_PROG_WARNING;
2057 }
2058
2059 /* use our all-purpose error reporting to get stuff back to the GUI */
2060 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
2061 }

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