/[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 1251 - (show annotations) (download) (as text)
Sat Jan 27 05:04:49 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 55325 byte(s)
Fixed 'hires' test (needed to rack down the tolerances a bit)
Removed newton test (until we make some changes to the solvers var list)
All of TestIDADENSE passes now and will be added to the buildbot list.
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 //CONSOLE_DEBUG("Raising signal...");
750 //CONSOLE_DEBUG("1/0 = %f", div1(1.0,0.0));
751 //CONSOLE_DEBUG("Still here...");
752
753 /* correct initial values, given derivatives */
754 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
755 /* note the new API from version 2.3 and onwards */
756 flag = IDACalcIC(ida_mem, icopt, tout1);
757 # else
758 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
759 # endif
760
761 switch(flag){
762 case IDA_SUCCESS:
763 CONSOLE_DEBUG("Initial conditions solved OK");
764 break;
765
766 case IDA_LSETUP_FAIL:
767 case IDA_LINIT_FAIL:
768 case IDA_LSOLVE_FAIL:
769 case IDA_NO_RECOVERY:
770 flag1 = -999;
771 flag = (flagfn)(ida_mem,&flag1);
772 if(flag){
773 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);
774 return 12;
775 }
776 ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
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
798 }
799
800 /* optionally, specify ROO-FINDING PROBLEM */
801
802 /* -- set up the IntegratorReporter */
803 integrator_output_init(blsys);
804
805 /* -- store the initial values of all the stuff */
806 integrator_output_write(blsys);
807 integrator_output_write_obs(blsys);
808
809 /* specify where the returned values should be stored */
810 yret = y0;
811 ypret = yp0;
812
813 /* advance solution in time, return values as yret and derivatives as ypret */
814 blsys->currentstep=1;
815 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
816 t = samplelist_get(blsys->samples, t_index);
817 t0 = integrator_get_t(blsys);
818 asc_assert(t > t0);
819
820 #ifdef SOLVE_DEBUG
821 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
822 #endif
823
824 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
825
826 /* pass the values of everything back to the compiler */
827 integrator_set_t(blsys, (double)tret);
828 integrator_set_y(blsys, NV_DATA_S(yret));
829 integrator_set_ydot(blsys, NV_DATA_S(ypret));
830
831 if(flag<0){
832 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
833 break;
834 }
835
836 /* -- do something so that blsys knows the values of tret, yret and ypret */
837
838 /* -- store the current values of all the stuff */
839 integrator_output_write(blsys);
840 integrator_output_write_obs(blsys);
841 }
842
843 /* -- close the IntegratorReporter */
844 integrator_output_close(blsys);
845
846 /* get optional outputs */
847 #ifdef STATS_DEBUG
848 IntegratorIdaStats stats;
849 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
850 integrator_ida_write_stats(&stats);
851 }
852 #endif
853
854 /* free solution memory */
855 N_VDestroy_Serial(yret);
856 N_VDestroy_Serial(ypret);
857
858 /* free solver memory */
859 IDAFree(ida_mem);
860
861 if(flag < 0){
862 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
863 return 14;
864 }
865
866 /* all done, success */
867 return 0;
868 }
869
870 /*--------------------------------------------------
871 RESIDUALS AND JACOBIAN
872 */
873
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
904 /**
905 Function to evaluate system residuals, in the form required for IDA.
906
907 Given tt, yy and yp, we need to evaluate and return rr.
908
909 @param tt current value of indep variable (time)
910 @param yy current values of dependent variable vector
911 @param yp current values of derivatives of dependent variables
912 @param rr the output residual vector (is we're returning data to)
913 @param res_data pointer to our stuff (blsys in this case).
914
915 @return 0 on success, positive on recoverable error, and
916 negative on unrecoverable error.
917 */
918 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
919 IntegratorSystem *blsys;
920 IntegratorIdaData *enginedata;
921 int i, calc_ok, is_error;
922 struct rel_relation** relptr;
923 double resid;
924 char *relname;
925 #ifdef FEX_DEBUG
926 char *varname;
927 char diffname[30];
928 #endif
929
930 blsys = (IntegratorSystem *)res_data;
931 enginedata = integrator_ida_enginedata(blsys);
932
933 #ifdef FEX_DEBUG
934 /* fprintf(stderr,"\n\n"); */
935 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
936 #endif
937
938 /* pass the values of everything back to the compiler */
939 integrator_set_t(blsys, (double)tt);
940 integrator_set_y(blsys, NV_DATA_S(yy));
941 integrator_set_ydot(blsys, NV_DATA_S(yp));
942
943 if(NV_LENGTH_S(rr)!=enginedata->nrels){
944 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
945 return -1; /* unrecoverable */
946 }
947
948 /**
949 @TODO does this function (fex) do bounds checking already?
950 */
951
952 /* evaluate each residual in the rellist */
953 is_error = 0;
954 relptr = enginedata->rellist;
955
956
957 #ifdef ASC_SIGNAL_TRAPS
958 if(enginedata->safeeval){
959 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
960 }else{
961 # ifdef FEX_DEBUG
962 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
963 # endif
964 Asc_SignalHandlerPushDefault(SIGFPE);
965 }
966
967 if (SETJMP(g_fpe_env)==0) {
968 #endif
969
970
971 /*
972 integrator_ida_sig_old = signal(SIGFPE,&integrator_ida_sig);
973 if(fegetenv(&integrator_ida_fenv_old))ASC_PANIC("fenv problem");
974 if(feenableexcept(FE_ALL_EXCEPT))ASC_PANIC("fenv problem");
975 switch(setjmp(integrator_ida_jmp_buf)){
976 case 0:
977 */
978 for(i=0, relptr = enginedata->rellist;
979 i< enginedata->nrels && relptr != NULL;
980 ++i, ++relptr
981 ){
982 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
983
984 NV_Ith_S(rr,i) = resid;
985 if(!calc_ok){
986 relname = rel_make_name(blsys->system, *relptr);
987 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
988 ASC_FREE(relname);
989 /* presumable some output already made? */
990 is_error = 1;
991 }/*else{
992 CONSOLE_DEBUG("Calc OK");
993 }*/
994 }
995 /* break;
996 case SIGFPE:
997 ERROR_REPORTER_HERE(ASC_PROG_ERR,"...here! Caught SIGFPE");
998 integrator_ida_write_feinfo();
999 is_error=1;
1000 break;
1001 default:
1002 ERROR_REPORTER_HERE(ASC_PROG_ERR,"What exception?");
1003 is_error=1;
1004 }
1005 if(signal(SIGFPE,integrator_ida_sig_old)!=&integrator_ida_sig)ASC_PANIC("signal problem");
1006 if(fesetenv(&integrator_ida_fenv_old))ASC_PANIC("fenv problem");
1007 feraiseexcept(FE_DIVBYZERO);
1008 */
1009 if(!is_error){
1010 for(i=0;i< enginedata->nrels; ++i){
1011 if(isnan(NV_Ith_S(rr,i))){
1012 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in residual %d",i);
1013 is_error=1;
1014 }
1015 }
1016 #ifdef FEX_DEBUG
1017 if(!is_error){
1018 CONSOLE_DEBUG("No NAN detected");
1019 }
1020 #endif
1021 }
1022
1023 #ifdef ASC_SIGNAL_TRAPS
1024 }else{
1025 relname = rel_make_name(blsys->system, *relptr);
1026 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1027 ASC_FREE(relname);
1028 is_error = 1;
1029 }
1030
1031 if(enginedata->safeeval){
1032 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1033 }else{
1034 Asc_SignalHandlerPopDefault(SIGFPE);
1035 }
1036 #endif
1037
1038
1039 #ifdef FEX_DEBUG
1040 /* output residuals to console */
1041 CONSOLE_DEBUG("RESIDUAL OUTPUT");
1042 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1043 for(i=0; i<blsys->n_y; ++i){
1044 varname = var_make_name(blsys->system,blsys->y[i]);
1045 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1046 if(blsys->ydot[i]){
1047 varname = var_make_name(blsys->system,blsys->ydot[i]);
1048 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1049 }else{
1050 snprintf(diffname,99,"diff(%s)",varname);
1051 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1052 }
1053 ASC_FREE(varname);
1054 relname = rel_make_name(blsys->system,enginedata->rellist[i]);
1055 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1056 }
1057 #endif
1058
1059 if(is_error){
1060 return 1;
1061 }
1062
1063 #ifdef FEX_DEBUG
1064 CONSOLE_DEBUG("RESIDUAL OK");
1065 #endif
1066 return 0;
1067 }
1068
1069 /**
1070 Dense Jacobian evaluation. Only suitable for small problems!
1071 */
1072 int integrator_ida_djex(long int Neq, realtype tt
1073 , N_Vector yy, N_Vector yp, N_Vector rr
1074 , realtype c_j, void *jac_data, DenseMat Jac
1075 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1076 ){
1077 IntegratorSystem *blsys;
1078 IntegratorIdaData *enginedata;
1079 char *relname;
1080 #ifdef DJEX_DEBUG
1081 struct var_variable **varlist;
1082 char *varname;
1083 #endif
1084 struct rel_relation **relptr;
1085 int i;
1086 double *derivatives;
1087 struct var_variable **variables;
1088 int count, j;
1089 int status, is_error = 0;
1090
1091 blsys = (IntegratorSystem *)jac_data;
1092 enginedata = integrator_ida_enginedata(blsys);
1093
1094 /* allocate space for returns from relman_diff3 */
1095 /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
1096 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1097 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1098
1099 /* pass the values of everything back to the compiler */
1100 integrator_set_t(blsys, (double)tt);
1101 integrator_set_y(blsys, NV_DATA_S(yy));
1102 integrator_set_ydot(blsys, NV_DATA_S(yp));
1103
1104 #ifdef DJEX_DEBUG
1105 varlist = slv_get_solvers_var_list(blsys->system);
1106
1107 /* print vars */
1108 for(i=0; i < blsys->n_y; ++i){
1109 varname = var_make_name(blsys->system, blsys->y[i]);
1110 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1111 asc_assert(NV_Ith_S(yy,i) == var_value(blsys->y[i]));
1112 ASC_FREE(varname);
1113 }
1114
1115 /* print derivatives */
1116 for(i=0; i < blsys->n_y; ++i){
1117 if(blsys->ydot[i]){
1118 varname = var_make_name(blsys->system, blsys->ydot[i]);
1119 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
1120 ASC_FREE(varname);
1121 }else{
1122 varname = var_make_name(blsys->system, blsys->y[i]);
1123 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1124 ASC_FREE(varname);
1125 }
1126 }
1127
1128 /* print step size */
1129 CONSOLE_DEBUG("<c_j> = %g",c_j);
1130 #endif
1131
1132 /* build up the dense jacobian matrix... */
1133 status = 0;
1134 for(i=0, relptr = enginedata->rellist;
1135 i< enginedata->nrels && relptr != NULL;
1136 ++i, ++relptr
1137 ){
1138 /* get derivatives for this particular relation */
1139 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1140
1141 if(status){
1142 relname = rel_make_name(blsys->system, *relptr);
1143 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1144 ASC_FREE(relname);
1145 is_error = 1;
1146 break;
1147 }
1148
1149 /* output what's going on here ... */
1150 #ifdef DJEX_DEBUG
1151 relname = rel_make_name(blsys->system, *relptr);
1152 fprintf(stderr,"%d: '%s': ",i,relname);
1153 for(j=0;j<count;++j){
1154 varname = var_make_name(blsys->system, variables[j]);
1155 if(var_deriv(variables[j])){
1156 fprintf(stderr," '%s'=",varname);
1157 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(blsys,variables[j]));
1158 }else{
1159 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1160 }
1161 ASC_FREE(varname);
1162 }
1163 /* relname is freed further down */
1164 fprintf(stderr,"\n");
1165 #endif
1166
1167 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1168 for(j=0; j < count; ++j){
1169 #ifdef DJEX_DEBUG
1170 varname = var_make_name(blsys->system,variables[j]);
1171 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1172 ASC_FREE(varname);
1173 #endif
1174 if(!var_deriv(variables[j])){
1175 #ifdef DJEX_DEBUG
1176 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
1177 asc_assert(var_sindex(variables[j]) >= 0);
1178 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1179 #endif
1180 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1181 }else{
1182 DENSE_ELEM(Jac,i,integrator_ida_diffindex(blsys,variables[j])) += derivatives[j] * c_j;
1183 #ifdef DJEX_DEBUG
1184 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1185 #endif
1186 }
1187 }
1188 }
1189
1190 #ifdef DJEX_DEBUG
1191 ASC_FREE(relname);
1192 CONSOLE_DEBUG("PRINTING JAC");
1193 fprintf(stderr,"\t");
1194 for(j=0; j < blsys->n_y; ++j){
1195 if(j)fprintf(stderr,"\t");
1196 varname = var_make_name(blsys->system,blsys->y[j]);
1197 fprintf(stderr,"%11s",varname);
1198 ASC_FREE(varname);
1199 }
1200 fprintf(stderr,"\n");
1201 for(i=0; i < enginedata->nrels; ++i){
1202 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
1203 fprintf(stderr,"%s\t",relname);
1204 ASC_FREE(relname);
1205
1206 for(j=0; j < blsys->n_y; ++j){
1207 if(j)fprintf(stderr,"\t");
1208 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1209 }
1210 fprintf(stderr,"\n");
1211 }
1212 #endif
1213
1214 /* test for NANs */
1215 if(!is_error){
1216 for(i=0;i< enginedata->nrels; ++i){
1217 for(j=0;j<blsys->n_y;++j){
1218 if(isnan(DENSE_ELEM(Jac,i,j))){
1219 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
1220 is_error=1;
1221 }
1222 }
1223 }
1224 #ifdef DJEX_DEBUG
1225 if(!is_error){
1226 CONSOLE_DEBUG("No NAN detected");
1227 }
1228 #endif
1229 }
1230
1231 /* if(integrator_ida_check_diffindex(blsys)){
1232 is_error = 1;
1233 }*/
1234
1235 if(is_error){
1236 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1237 return 1;
1238 }
1239
1240 #ifdef DJEX_DEBUG
1241 CONSOLE_DEBUG("DJEX RETURNING 0");
1242 /* ASC_PANIC("Quitting"); */
1243 #endif
1244 return 0;
1245 }
1246
1247 /**
1248 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1249
1250 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1251
1252 @param tt current value of the independent variable (time, t)
1253 @param yy current value of the dependent variable vector, y(t).
1254 @param yp current value of y'(t).
1255 @param rr current value of the residual vector F(t, y, y').
1256 @param v the vector by which the Jacobian must be multiplied to the right.
1257 @param Jv the output vector computed
1258 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1259 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1260 @param tmp1 @see tmp2
1261 @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.
1262 @return 0 on success
1263 */
1264 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1265 , N_Vector v, N_Vector Jv, realtype c_j
1266 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1267 ){
1268 IntegratorSystem *blsys;
1269 IntegratorIdaData *enginedata;
1270 int i, j, is_error=0;
1271 struct rel_relation** relptr;
1272 char *relname;
1273 int status;
1274 double Jv_i;
1275
1276 struct var_variable **variables;
1277 double *derivatives;
1278 int count;
1279 struct var_variable **varlist;
1280 #ifdef JEX_DEBUG
1281
1282 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1283 #endif
1284
1285 blsys = (IntegratorSystem *)jac_data;
1286 enginedata = integrator_ida_enginedata(blsys);
1287 varlist = slv_get_solvers_var_list(blsys->system);
1288
1289 /* pass the values of everything back to the compiler */
1290 integrator_set_t(blsys, (double)tt);
1291 integrator_set_y(blsys, NV_DATA_S(yy));
1292 integrator_set_ydot(blsys, NV_DATA_S(yp));
1293 /* no real use for residuals (rr) here, I don't think? */
1294
1295 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1296
1297 i = NV_LENGTH_S(yy) * 2;
1298 #ifdef JEX_DEBUG
1299 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1300 #endif
1301 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1302 derivatives = ASC_NEW_ARRAY(double, i);
1303
1304 /* evaluate the derivatives... */
1305 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1306
1307 #ifdef ASC_SIGNAL_TRAPS
1308 Asc_SignalHandlerPushDefault(SIGFPE);
1309 if (SETJMP(g_fpe_env)==0) {
1310 #endif
1311 for(i=0, relptr = enginedata->rellist;
1312 i< enginedata->nrels && relptr != NULL;
1313 ++i, ++relptr
1314 ){
1315 /* get derivatives for this particular relation */
1316 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1317 #ifdef JEX_DEBUG
1318 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1319 #endif
1320
1321 if(status){
1322 relname = rel_make_name(blsys->system, *relptr);
1323 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1324 ASC_FREE(relname);
1325 is_error = 1;
1326 break;
1327 }
1328
1329 /*
1330 Now we have the derivatives wrt each alg/diff variable in the
1331 present equation. variables[] points into the varlist. need
1332 a mapping from the varlist to the y and ydot lists.
1333 */
1334
1335 Jv_i = 0;
1336 for(j=0; j < count; ++j){
1337 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
1338 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1339 if(varname){
1340 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1341 ASC_FREE(varname);
1342 }else{
1343 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1344 }
1345 */
1346
1347 /* we don't calculate derivatives wrt indep var */
1348 asc_assert(variables[j]>=0);
1349 if(variables[j] == blsys->x) continue;
1350 #ifdef JEX_DEBUG
1351 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1352 #endif
1353 if(var_deriv(variables[j])){
1354 #define DIFFINDEX integrator_ida_diffindex(blsys,variables[j])
1355 #ifdef JEX_DEBUG
1356 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
1357 , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
1358 , i, DIFFINDEX, derivatives[j]
1359 , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
1360 );
1361 #endif
1362 asc_assert(blsys->ydot[DIFFINDEX]==variables[j]);
1363 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1364 #undef DIFFINDEX
1365 }else{
1366 #define VARINDEX var_sindex(variables[j])
1367 #ifdef JEX_DEBUG
1368 asc_assert(blsys->y[VARINDEX]==variables[j]);
1369 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
1370 , i
1371 , derivatives[j] * NV_Ith_S(v,VARINDEX)
1372 , i, VARINDEX, derivatives[j]
1373 , VARINDEX, NV_Ith_S(v,VARINDEX)
1374 );
1375 #endif
1376 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1377 #undef VARINDEX
1378 }
1379 }
1380
1381 NV_Ith_S(Jv,i) = Jv_i;
1382 #ifdef JEX_DEBUG
1383 CONSOLE_DEBUG("rel = %p",*relptr);
1384 relname = rel_make_name(blsys->system, *relptr);
1385 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1386 ASC_FREE(relname);
1387 return 1;
1388 #endif
1389 }
1390 #ifdef ASC_SIGNAL_TRAPS
1391 }else{
1392 relname = rel_make_name(blsys->system, *relptr);
1393 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1394 ASC_FREE(relname);
1395 is_error = 1;
1396 }
1397 Asc_SignalHandlerPopDefault(SIGFPE);
1398 #endif
1399
1400 if(is_error){
1401 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1402 return 1;
1403 }
1404 return 0;
1405 }
1406
1407 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1408 int integrator_ida_sjex(long int Neq, realtype tt
1409 , N_Vector yy, N_Vector yp, N_Vector rr
1410 , realtype c_j, void *jac_data, mtx_matrix_t Jac
1411 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1412 ){
1413 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1414 return -1;
1415 }
1416
1417 /*----------------------------------------------
1418 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1419 */
1420
1421 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys){
1422 IntegratorIdaData * enginedata =blsys->enginedata;
1423 IntegratorIdaPrecDataJacobi *precdata;
1424 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1425
1426 asc_assert(blsys->n_y);
1427 precdata->PIii = N_VNew_Serial(blsys->n_y);
1428
1429 enginedata->pfree = &integrator_ida_pfree_jacobi;
1430 enginedata->precdata = precdata;
1431 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1432 }
1433
1434 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1435 if(enginedata->precdata){
1436 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1437 N_VDestroy_Serial(precdata->PIii);
1438
1439 ASC_FREE(precdata);
1440 enginedata->precdata = NULL;
1441 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1442 }
1443 enginedata->pfree = NULL;
1444 }
1445
1446 /**
1447 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1448
1449 'setup' function.
1450 */
1451 int integrator_ida_psetup_jacobi(realtype tt,
1452 N_Vector yy, N_Vector yp, N_Vector rr,
1453 realtype c_j, void *p_data,
1454 N_Vector tmp1, N_Vector tmp2,
1455 N_Vector tmp3
1456 ){
1457 int i, j, res;
1458 IntegratorSystem *blsys;
1459 IntegratorIdaData *enginedata;
1460 IntegratorIdaPrecDataJacobi *precdata;
1461 struct rel_relation **relptr;
1462
1463 blsys = (IntegratorSystem *)p_data;
1464 enginedata = blsys->enginedata;
1465 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1466 double *derivatives;
1467 struct var_variable **variables;
1468 int count, status;
1469 char *relname;
1470
1471 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1472
1473 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1474 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1475
1476 /**
1477 @TODO FIXME here we are using the very inefficient and contorted approach
1478 of calculating the whole jacobian, then extracting just the diagonal elements.
1479 */
1480
1481 for(i=0, relptr = enginedata->rellist;
1482 i< enginedata->nrels && relptr != NULL;
1483 ++i, ++relptr
1484 ){
1485
1486 /* get derivatives for this particular relation */
1487 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1488 if(status){
1489 relname = rel_make_name(blsys->system, *relptr);
1490 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1491 ASC_FREE(relname);
1492 break;
1493 }
1494 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1495 /* find the diagonal elements */
1496 for(j=0; j<count; ++j){
1497 if(var_sindex(variables[j])==i){
1498 if(var_deriv(variables[j])){
1499 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
1500 }else{
1501 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1502 }
1503
1504 }
1505 }
1506 #ifdef PREC_DEBUG
1507 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1508 #endif
1509 }
1510
1511 if(status){
1512 CONSOLE_DEBUG("Error found when evaluating derivatives");
1513 res = 1; goto finish; /* recoverable */
1514 }
1515
1516 integrator_ida_write_incidence(blsys);
1517
1518 res = 0;
1519 finish:
1520 ASC_FREE(variables);
1521 ASC_FREE(derivatives);
1522 return res;
1523 };
1524
1525 /**
1526 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1527
1528 'solve' function.
1529 */
1530 int integrator_ida_psolve_jacobi(realtype tt,
1531 N_Vector yy, N_Vector yp, N_Vector rr,
1532 N_Vector rvec, N_Vector zvec,
1533 realtype c_j, realtype delta, void *p_data,
1534 N_Vector tmp
1535 ){
1536 IntegratorSystem *blsys;
1537 IntegratorIdaData *data;
1538 IntegratorIdaPrecDataJacobi *precdata;
1539 blsys = (IntegratorSystem *)p_data;
1540 data = blsys->enginedata;
1541 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1542
1543 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1544 N_VProd(precdata->PIii, rvec, zvec);
1545 return 0;
1546 };
1547
1548 /*----------------------------------------------
1549 STATS
1550 */
1551
1552 /**
1553 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1554 status in a struct instead of separately.
1555
1556 @return IDA_SUCCESS on sucess.
1557 */
1558 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1559 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1560 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1561 ,&s->hlast, &s->hcur, &s->tcur
1562 );
1563 }
1564
1565 /**
1566 This routine just outputs the stats to the CONSOLE_DEBUG routine.
1567
1568 @TODO provide a GUI way of stats reporting from IDA.
1569 */
1570 void integrator_ida_write_stats(IntegratorIdaStats *stats){
1571 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
1572 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
1573 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
1574 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
1575 SI(qlast); SI(qcur);
1576 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
1577 # undef SL
1578 # undef SI
1579 # undef SR
1580 }
1581
1582 /*------------------------------------------------------------------------------
1583 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
1584 */
1585
1586 enum integrator_ida_write_jac_enum{
1587 II_WRITE_Y
1588 , II_WRITE_YDOT
1589 };
1590
1591 /**
1592 @TODO COMPLETE THIS...
1593 */
1594 void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){
1595 IntegratorIdaData *enginedata;
1596 MM_typecode matcode;
1597 int nnz, rhomax;
1598 double *derivatives;
1599 struct var_variable **variables;
1600 struct rel_relation **relptr;
1601 int i, j, status, count;
1602 char *relname;
1603
1604 var_filter_t vfilter = {
1605 VAR_SVAR | VAR_FIXED | VAR_DERIV
1606 ,VAR_SVAR | 0 | 0
1607 };
1608 enginedata = (IntegratorIdaData *)blsys->enginedata;
1609
1610 /* number of non-zeros for all the non-FIXED solver_vars,
1611 in all the active included equality relations.
1612 */
1613 nnz = relman_jacobian_count(enginedata->rellist, enginedata->nrels
1614 , &enginedata->vfilter, &enginedata->rfilter
1615 , &rhomax
1616 );
1617
1618 /* we must have found the same number of relations */
1619 asc_assert(rhomax == enginedata->nrels);
1620
1621 /* output the mmio file header, now that we know our size*/
1622 /* note that we are asserting that our problem is square */
1623 mm_initialize_typecode(&matcode);
1624 mm_set_matrix(&matcode);
1625 mm_set_coordinate(&matcode);
1626 mm_set_real(&matcode);
1627 mm_write_banner(f, matcode);
1628 mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);
1629
1630 variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);
1631 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1632
1633 CONSOLE_DEBUG("Writing sparse Jacobian to file...");
1634
1635 for(i=0, relptr = enginedata->rellist;
1636 i< enginedata->nrels && relptr != NULL;
1637 ++i, ++relptr
1638 ){
1639 relname = rel_make_name(blsys->system, *relptr);
1640
1641 /* get derivatives of y */
1642 status = relman_diff3(*relptr, &vfilter, derivatives, variables, &count, enginedata->safeeval);
1643 if(status){
1644 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1645 ASC_FREE(relname);
1646 break;
1647 }
1648
1649 if(type==II_WRITE_YDOT){
1650 for(j=0; j<count; ++j){
1651 if(var_deriv(variables[j])){
1652 fprintf(f, "%d %d %10.3g\n", i, integrator_ida_diffindex(blsys,variables[j]), derivatives[j]);
1653 }
1654 }
1655 }else if(type==II_WRITE_Y){
1656 for(j=0; j<count; ++j){
1657 if(!var_deriv(variables[j])){
1658 fprintf(f, "%d %d %10.3g\n", i, var_sindex(variables[j]), derivatives[j]);
1659 }
1660 }
1661 }
1662 }
1663 ASC_FREE(variables);
1664 ASC_FREE(derivatives);
1665 }
1666
1667 /**
1668 This routine outputs matrix structure in a crude text format, for the sake
1669 of debugging.
1670 */
1671 void integrator_ida_write_incidence(IntegratorSystem *blsys){
1672 int i, j;
1673 struct rel_relation **relptr;
1674 IntegratorIdaData *enginedata = blsys->enginedata;
1675 double *derivatives;
1676 struct var_variable **variables;
1677 int count, status;
1678 char *relname;
1679
1680 if(enginedata->nrels > 100){
1681 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
1682 return;
1683 }
1684
1685 variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);
1686 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1687
1688 CONSOLE_DEBUG("Outputting incidence information to console...");
1689
1690 for(i=0, relptr = enginedata->rellist;
1691 i< enginedata->nrels && relptr != NULL;
1692 ++i, ++relptr
1693 ){
1694 relname = rel_make_name(blsys->system, *relptr);
1695
1696 /* get derivatives for this particular relation */
1697 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1698 if(status){
1699 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1700 ASC_FREE(relname);
1701 break;
1702 }
1703
1704 fprintf(stderr,"%3d:%-15s:",i,relname);
1705 ASC_FREE(relname);
1706
1707 for(j=0; j<count; ++j){
1708 if(var_deriv(variables[j])){
1709 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(blsys,variables[j]));
1710 }else{
1711 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
1712 }
1713 }
1714 fprintf(stderr,"\n");
1715 }
1716 ASC_FREE(variables);
1717 ASC_FREE(derivatives);
1718 }
1719
1720 /* @return 0 on success */
1721 int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
1722 char *varname, *relname;
1723 struct var_variable **vlist, *var;
1724 struct rel_relation **rlist, *rel;
1725 long vlen, rlen;
1726 long i;
1727
1728 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",sys->n_y);
1729
1730 /* if(integrator_sort_obs_vars(sys))return 10; */
1731
1732 if(sys->y && sys->ydot){
1733 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
1734 fprintf(fp,"index\ty\tydot\n");
1735 fprintf(fp,"-----\t-----\t-----\n");
1736 for(i=0;i<sys->n_y;++i){
1737 varname = var_make_name(sys->system, sys->y[i]);
1738 fprintf(fp,"%ld\t%s\t",i,varname);
1739 if(sys->ydot[i]){
1740 ASC_FREE(varname);
1741 varname = var_make_name(sys->system, sys->ydot[i]);
1742 fprintf(fp,"%s\n",varname);
1743 ASC_FREE(varname);
1744 }else{
1745 fprintf(fp,".\n");
1746 ASC_FREE(varname);
1747 }
1748 }
1749 }else{
1750 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
1751 }
1752
1753 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
1754 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
1755 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
1756
1757
1758 /* visit all the slv_system_t master var lists to collect vars */
1759 /* find the vars mostly in this one */
1760 vlist = slv_get_solvers_var_list(sys->system);
1761 vlen = slv_get_num_solvers_vars(sys->system);
1762 for(i=0;i<vlen;i++){
1763 var = vlist[i];
1764
1765 varname = var_make_name(sys->system, var);
1766 fprintf(fp,"%ld\t%-15s\t",i,varname);
1767
1768 if(var_fixed(var)){
1769 // it's fixed, so not really a DAE var
1770 fprintf(fp,"(fixed)\n");
1771 }else if(!var_active(var)){
1772 // inactive
1773 fprintf(fp,"(inactive)\n");
1774 }else if(!var_incident(var)){
1775 // not incident
1776 fprintf(fp,"(not incident)\n");
1777 }else{
1778 if(var_deriv(var)){
1779 if(sys->y_id){
1780 ASC_FREE(varname);
1781 varname = var_make_name(sys->system,vlist[integrator_ida_diffindex(sys,var)]);
1782 fprintf(fp,".\tdiff(%d='%s')\n",integrator_ida_diffindex(sys,var),varname);
1783 }else{
1784 fprintf(fp,".\tderiv... of??\n");
1785 }
1786 }else{
1787 fprintf(fp,"%d\t.\n",var_sindex(var));
1788 }
1789 }
1790 ASC_FREE(varname);
1791 }
1792
1793 /* let's write out the relations too */
1794 rlist = slv_get_solvers_rel_list(sys->system);
1795 rlen = slv_get_num_solvers_rels(sys->system);
1796
1797 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
1798 fprintf(fp,"index\tname\n");
1799 fprintf(fp,"-----\t----\n");
1800 for(i=0; i<rlen; ++i){
1801 rel = rlist[i];
1802 relname = rel_make_name(sys->system,rel);
1803 fprintf(fp,"%ld\t%s\n",i,relname);
1804 ASC_FREE(relname);
1805 }
1806
1807 /* write out the derivative chains */
1808 fprintf(fp,"\nDERIVATIVE CHAINS\n");
1809 if(integrator_ida_analyse_debug(sys,stderr)){
1810 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
1811 return 340;
1812 }
1813 fprintf(fp,"\n");
1814
1815 /* and lets write block debug output */
1816 system_block_debug(sys->system, fp);
1817
1818 return 0; /* success */
1819 }
1820
1821 /*----------------------------------------------
1822 ERROR REPORTING
1823 */
1824 /**
1825 Error message reporter function to be passed to IDA. All error messages
1826 will trigger a call to this function, so we should find everything
1827 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1828 panel (in the case of PyGTK).
1829 */
1830 void integrator_ida_error(int error_code
1831 , const char *module, const char *function
1832 , char *msg, void *eh_data
1833 ){
1834 IntegratorSystem *blsys;
1835 error_severity_t sev;
1836
1837 /* cast back the IntegratorSystem, just in case we need it */
1838 blsys = (IntegratorSystem *)eh_data;
1839
1840 /* severity depends on the sign of the error_code value */
1841 if(error_code <= 0){
1842 sev = ASC_PROG_ERR;
1843 }else{
1844 sev = ASC_PROG_WARNING;
1845 }
1846
1847 /* use our all-purpose error reporting to get stuff back to the GUI */
1848 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1849 }

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