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

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