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

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