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

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