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

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