/[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 1279 - (show annotations) (download) (as text)
Sat Feb 10 09:31:55 2007 UTC (17 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 55930 byte(s)
Removed some debug output.
DSGSAT model L=10, n=47
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 */
1092 int integrator_ida_djex(long int Neq, realtype tt
1093 , N_Vector yy, N_Vector yp, N_Vector rr
1094 , realtype c_j, void *jac_data, DenseMat Jac
1095 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1096 ){
1097 IntegratorSystem *blsys;
1098 IntegratorIdaData *enginedata;
1099 char *relname;
1100 #ifdef DJEX_DEBUG
1101 struct var_variable **varlist;
1102 char *varname;
1103 #endif
1104 struct rel_relation **relptr;
1105 int i;
1106 double *derivatives;
1107 struct var_variable **variables;
1108 int count, j;
1109 int status, is_error = 0;
1110
1111 blsys = (IntegratorSystem *)jac_data;
1112 enginedata = integrator_ida_enginedata(blsys);
1113
1114 /* allocate space for returns from relman_diff3 */
1115 /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
1116 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1117 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1118
1119 /* pass the values of everything back to the compiler */
1120 integrator_set_t(blsys, (double)tt);
1121 integrator_set_y(blsys, NV_DATA_S(yy));
1122 integrator_set_ydot(blsys, NV_DATA_S(yp));
1123
1124 /* perform bounds checking on all variables */
1125 if(slv_check_bounds(blsys->system, 0, -1)){
1126 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1127 return 1;
1128 }
1129
1130 #ifdef DJEX_DEBUG
1131 varlist = slv_get_solvers_var_list(blsys->system);
1132
1133 /* print vars */
1134 for(i=0; i < blsys->n_y; ++i){
1135 varname = var_make_name(blsys->system, blsys->y[i]);
1136 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1137 asc_assert(NV_Ith_S(yy,i) == var_value(blsys->y[i]));
1138 ASC_FREE(varname);
1139 }
1140
1141 /* print derivatives */
1142 for(i=0; i < blsys->n_y; ++i){
1143 if(blsys->ydot[i]){
1144 varname = var_make_name(blsys->system, blsys->ydot[i]);
1145 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
1146 ASC_FREE(varname);
1147 }else{
1148 varname = var_make_name(blsys->system, blsys->y[i]);
1149 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1150 ASC_FREE(varname);
1151 }
1152 }
1153
1154 /* print step size */
1155 CONSOLE_DEBUG("<c_j> = %g",c_j);
1156 #endif
1157
1158 /* build up the dense jacobian matrix... */
1159 status = 0;
1160 for(i=0, relptr = enginedata->rellist;
1161 i< enginedata->nrels && relptr != NULL;
1162 ++i, ++relptr
1163 ){
1164 /* get derivatives for this particular relation */
1165 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1166
1167 if(status){
1168 relname = rel_make_name(blsys->system, *relptr);
1169 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1170 ASC_FREE(relname);
1171 is_error = 1;
1172 break;
1173 }
1174
1175 /* output what's going on here ... */
1176 #ifdef DJEX_DEBUG
1177 relname = rel_make_name(blsys->system, *relptr);
1178 fprintf(stderr,"%d: '%s': ",i,relname);
1179 for(j=0;j<count;++j){
1180 varname = var_make_name(blsys->system, variables[j]);
1181 if(var_deriv(variables[j])){
1182 fprintf(stderr," '%s'=",varname);
1183 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(blsys,variables[j]));
1184 }else{
1185 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1186 }
1187 ASC_FREE(varname);
1188 }
1189 /* relname is freed further down */
1190 fprintf(stderr,"\n");
1191 #endif
1192
1193 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1194 for(j=0; j < count; ++j){
1195 #ifdef DJEX_DEBUG
1196 varname = var_make_name(blsys->system,variables[j]);
1197 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1198 ASC_FREE(varname);
1199 #endif
1200 if(!var_deriv(variables[j])){
1201 #ifdef DJEX_DEBUG
1202 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
1203 asc_assert(var_sindex(variables[j]) >= 0);
1204 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1205 #endif
1206 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1207 }else{
1208 DENSE_ELEM(Jac,i,integrator_ida_diffindex(blsys,variables[j])) += derivatives[j] * c_j;
1209 #ifdef DJEX_DEBUG
1210 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1211 #endif
1212 }
1213 }
1214 }
1215
1216 #ifdef DJEX_DEBUG
1217 ASC_FREE(relname);
1218 CONSOLE_DEBUG("PRINTING JAC");
1219 fprintf(stderr,"\t");
1220 for(j=0; j < blsys->n_y; ++j){
1221 if(j)fprintf(stderr,"\t");
1222 varname = var_make_name(blsys->system,blsys->y[j]);
1223 fprintf(stderr,"%11s",varname);
1224 ASC_FREE(varname);
1225 }
1226 fprintf(stderr,"\n");
1227 for(i=0; i < enginedata->nrels; ++i){
1228 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
1229 fprintf(stderr,"%s\t",relname);
1230 ASC_FREE(relname);
1231
1232 for(j=0; j < blsys->n_y; ++j){
1233 if(j)fprintf(stderr,"\t");
1234 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1235 }
1236 fprintf(stderr,"\n");
1237 }
1238 #endif
1239
1240 /* test for NANs */
1241 if(!is_error){
1242 for(i=0;i< enginedata->nrels; ++i){
1243 for(j=0;j<blsys->n_y;++j){
1244 if(isnan(DENSE_ELEM(Jac,i,j))){
1245 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
1246 is_error=1;
1247 }
1248 }
1249 }
1250 #ifdef DJEX_DEBUG
1251 if(!is_error){
1252 CONSOLE_DEBUG("No NAN detected");
1253 }
1254 #endif
1255 }
1256
1257 /* if(integrator_ida_check_diffindex(blsys)){
1258 is_error = 1;
1259 }*/
1260
1261 if(is_error){
1262 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1263 return 1;
1264 }
1265
1266 #ifdef DJEX_DEBUG
1267 CONSOLE_DEBUG("DJEX RETURNING 0");
1268 /* ASC_PANIC("Quitting"); */
1269 #endif
1270 return 0;
1271 }
1272
1273 /**
1274 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1275
1276 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1277
1278 @param tt current value of the independent variable (time, t)
1279 @param yy current value of the dependent variable vector, y(t).
1280 @param yp current value of y'(t).
1281 @param rr current value of the residual vector F(t, y, y').
1282 @param v the vector by which the Jacobian must be multiplied to the right.
1283 @param Jv the output vector computed
1284 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1285 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1286 @param tmp1 @see tmp2
1287 @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.
1288 @return 0 on success
1289 */
1290 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1291 , N_Vector v, N_Vector Jv, realtype c_j
1292 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1293 ){
1294 IntegratorSystem *blsys;
1295 IntegratorIdaData *enginedata;
1296 int i, j, is_error=0;
1297 struct rel_relation** relptr;
1298 char *relname;
1299 int status;
1300 double Jv_i;
1301
1302 struct var_variable **variables;
1303 double *derivatives;
1304 int count;
1305 struct var_variable **varlist;
1306 #ifdef JEX_DEBUG
1307
1308 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1309 #endif
1310
1311 blsys = (IntegratorSystem *)jac_data;
1312 enginedata = integrator_ida_enginedata(blsys);
1313 varlist = slv_get_solvers_var_list(blsys->system);
1314
1315 /* pass the values of everything back to the compiler */
1316 integrator_set_t(blsys, (double)tt);
1317 integrator_set_y(blsys, NV_DATA_S(yy));
1318 integrator_set_ydot(blsys, NV_DATA_S(yp));
1319 /* no real use for residuals (rr) here, I don't think? */
1320
1321 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1322
1323 i = NV_LENGTH_S(yy) * 2;
1324 #ifdef JEX_DEBUG
1325 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1326 #endif
1327 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1328 derivatives = ASC_NEW_ARRAY(double, i);
1329
1330 /* evaluate the derivatives... */
1331 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1332
1333 #ifdef ASC_SIGNAL_TRAPS
1334 Asc_SignalHandlerPushDefault(SIGFPE);
1335 if (SETJMP(g_fpe_env)==0) {
1336 #endif
1337 for(i=0, relptr = enginedata->rellist;
1338 i< enginedata->nrels && relptr != NULL;
1339 ++i, ++relptr
1340 ){
1341 /* get derivatives for this particular relation */
1342 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1343 #ifdef JEX_DEBUG
1344 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1345 #endif
1346
1347 if(status){
1348 relname = rel_make_name(blsys->system, *relptr);
1349 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1350 ASC_FREE(relname);
1351 is_error = 1;
1352 break;
1353 }
1354
1355 /*
1356 Now we have the derivatives wrt each alg/diff variable in the
1357 present equation. variables[] points into the varlist. need
1358 a mapping from the varlist to the y and ydot lists.
1359 */
1360
1361 Jv_i = 0;
1362 for(j=0; j < count; ++j){
1363 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
1364 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1365 if(varname){
1366 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1367 ASC_FREE(varname);
1368 }else{
1369 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1370 }
1371 */
1372
1373 /* we don't calculate derivatives wrt indep var */
1374 asc_assert(variables[j]>=0);
1375 if(variables[j] == blsys->x) continue;
1376 #ifdef JEX_DEBUG
1377 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1378 #endif
1379 if(var_deriv(variables[j])){
1380 #define DIFFINDEX integrator_ida_diffindex(blsys,variables[j])
1381 #ifdef JEX_DEBUG
1382 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
1383 , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
1384 , i, DIFFINDEX, derivatives[j]
1385 , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
1386 );
1387 #endif
1388 asc_assert(blsys->ydot[DIFFINDEX]==variables[j]);
1389 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1390 #undef DIFFINDEX
1391 }else{
1392 #define VARINDEX var_sindex(variables[j])
1393 #ifdef JEX_DEBUG
1394 asc_assert(blsys->y[VARINDEX]==variables[j]);
1395 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
1396 , i
1397 , derivatives[j] * NV_Ith_S(v,VARINDEX)
1398 , i, VARINDEX, derivatives[j]
1399 , VARINDEX, NV_Ith_S(v,VARINDEX)
1400 );
1401 #endif
1402 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1403 #undef VARINDEX
1404 }
1405 }
1406
1407 NV_Ith_S(Jv,i) = Jv_i;
1408 #ifdef JEX_DEBUG
1409 CONSOLE_DEBUG("rel = %p",*relptr);
1410 relname = rel_make_name(blsys->system, *relptr);
1411 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1412 ASC_FREE(relname);
1413 return 1;
1414 #endif
1415 }
1416 #ifdef ASC_SIGNAL_TRAPS
1417 }else{
1418 relname = rel_make_name(blsys->system, *relptr);
1419 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1420 ASC_FREE(relname);
1421 is_error = 1;
1422 }
1423 Asc_SignalHandlerPopDefault(SIGFPE);
1424 #endif
1425
1426 if(is_error){
1427 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1428 return 1;
1429 }
1430 return 0;
1431 }
1432
1433 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1434 int integrator_ida_sjex(long int Neq, realtype tt
1435 , N_Vector yy, N_Vector yp, N_Vector rr
1436 , realtype c_j, void *jac_data, mtx_matrix_t Jac
1437 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1438 ){
1439 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1440 return -1;
1441 }
1442
1443 /*----------------------------------------------
1444 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1445 */
1446
1447 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys){
1448 IntegratorIdaData * enginedata =blsys->enginedata;
1449 IntegratorIdaPrecDataJacobi *precdata;
1450 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1451
1452 asc_assert(blsys->n_y);
1453 precdata->PIii = N_VNew_Serial(blsys->n_y);
1454
1455 enginedata->pfree = &integrator_ida_pfree_jacobi;
1456 enginedata->precdata = precdata;
1457 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1458 }
1459
1460 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1461 if(enginedata->precdata){
1462 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1463 N_VDestroy_Serial(precdata->PIii);
1464
1465 ASC_FREE(precdata);
1466 enginedata->precdata = NULL;
1467 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1468 }
1469 enginedata->pfree = NULL;
1470 }
1471
1472 /**
1473 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1474
1475 'setup' function.
1476 */
1477 int integrator_ida_psetup_jacobi(realtype tt,
1478 N_Vector yy, N_Vector yp, N_Vector rr,
1479 realtype c_j, void *p_data,
1480 N_Vector tmp1, N_Vector tmp2,
1481 N_Vector tmp3
1482 ){
1483 int i, j, res;
1484 IntegratorSystem *blsys;
1485 IntegratorIdaData *enginedata;
1486 IntegratorIdaPrecDataJacobi *precdata;
1487 struct rel_relation **relptr;
1488
1489 blsys = (IntegratorSystem *)p_data;
1490 enginedata = blsys->enginedata;
1491 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1492 double *derivatives;
1493 struct var_variable **variables;
1494 int count, status;
1495 char *relname;
1496
1497 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1498
1499 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1500 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1501
1502 /**
1503 @TODO FIXME here we are using the very inefficient and contorted approach
1504 of calculating the whole jacobian, then extracting just the diagonal elements.
1505 */
1506
1507 for(i=0, relptr = enginedata->rellist;
1508 i< enginedata->nrels && relptr != NULL;
1509 ++i, ++relptr
1510 ){
1511
1512 /* get derivatives for this particular relation */
1513 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1514 if(status){
1515 relname = rel_make_name(blsys->system, *relptr);
1516 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1517 ASC_FREE(relname);
1518 break;
1519 }
1520 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1521 /* find the diagonal elements */
1522 for(j=0; j<count; ++j){
1523 if(var_sindex(variables[j])==i){
1524 if(var_deriv(variables[j])){
1525 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
1526 }else{
1527 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1528 }
1529
1530 }
1531 }
1532 #ifdef PREC_DEBUG
1533 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1534 #endif
1535 }
1536
1537 if(status){
1538 CONSOLE_DEBUG("Error found when evaluating derivatives");
1539 res = 1; goto finish; /* recoverable */
1540 }
1541
1542 integrator_ida_write_incidence(blsys);
1543
1544 res = 0;
1545 finish:
1546 ASC_FREE(variables);
1547 ASC_FREE(derivatives);
1548 return res;
1549 };
1550
1551 /**
1552 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1553
1554 'solve' function.
1555 */
1556 int integrator_ida_psolve_jacobi(realtype tt,
1557 N_Vector yy, N_Vector yp, N_Vector rr,
1558 N_Vector rvec, N_Vector zvec,
1559 realtype c_j, realtype delta, void *p_data,
1560 N_Vector tmp
1561 ){
1562 IntegratorSystem *blsys;
1563 IntegratorIdaData *data;
1564 IntegratorIdaPrecDataJacobi *precdata;
1565 blsys = (IntegratorSystem *)p_data;
1566 data = blsys->enginedata;
1567 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1568
1569 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1570 N_VProd(precdata->PIii, rvec, zvec);
1571 return 0;
1572 };
1573
1574 /*----------------------------------------------
1575 STATS
1576 */
1577
1578 /**
1579 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1580 status in a struct instead of separately.
1581
1582 @return IDA_SUCCESS on sucess.
1583 */
1584 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1585 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1586 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1587 ,&s->hlast, &s->hcur, &s->tcur
1588 );
1589 }
1590
1591 /**
1592 This routine just outputs the stats to the CONSOLE_DEBUG routine.
1593
1594 @TODO provide a GUI way of stats reporting from IDA.
1595 */
1596 void integrator_ida_write_stats(IntegratorIdaStats *stats){
1597 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
1598 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
1599 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
1600 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
1601 SI(qlast); SI(qcur);
1602 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
1603 # undef SL
1604 # undef SI
1605 # undef SR
1606 }
1607
1608 /*------------------------------------------------------------------------------
1609 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
1610 */
1611
1612 enum integrator_ida_write_jac_enum{
1613 II_WRITE_Y
1614 , II_WRITE_YDOT
1615 };
1616
1617 /**
1618 @TODO COMPLETE THIS...
1619 */
1620 void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){
1621 IntegratorIdaData *enginedata;
1622 MM_typecode matcode;
1623 int nnz, rhomax;
1624 double *derivatives;
1625 struct var_variable **variables;
1626 struct rel_relation **relptr;
1627 int i, j, status, count;
1628 char *relname;
1629
1630 var_filter_t vfilter = {
1631 VAR_SVAR | VAR_FIXED | VAR_DERIV
1632 ,VAR_SVAR | 0 | 0
1633 };
1634 enginedata = (IntegratorIdaData *)blsys->enginedata;
1635
1636 /* number of non-zeros for all the non-FIXED solver_vars,
1637 in all the active included equality relations.
1638 */
1639 nnz = relman_jacobian_count(enginedata->rellist, enginedata->nrels
1640 , &enginedata->vfilter, &enginedata->rfilter
1641 , &rhomax
1642 );
1643
1644 /* we must have found the same number of relations */
1645 asc_assert(rhomax == enginedata->nrels);
1646
1647 /* output the mmio file header, now that we know our size*/
1648 /* note that we are asserting that our problem is square */
1649 mm_initialize_typecode(&matcode);
1650 mm_set_matrix(&matcode);
1651 mm_set_coordinate(&matcode);
1652 mm_set_real(&matcode);
1653 mm_write_banner(f, matcode);
1654 mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);
1655
1656 variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);
1657 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1658
1659 CONSOLE_DEBUG("Writing sparse Jacobian to file...");
1660
1661 for(i=0, relptr = enginedata->rellist;
1662 i< enginedata->nrels && relptr != NULL;
1663 ++i, ++relptr
1664 ){
1665 relname = rel_make_name(blsys->system, *relptr);
1666
1667 /* get derivatives of y */
1668 status = relman_diff3(*relptr, &vfilter, derivatives, variables, &count, enginedata->safeeval);
1669 if(status){
1670 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1671 ASC_FREE(relname);
1672 break;
1673 }
1674
1675 if(type==II_WRITE_YDOT){
1676 for(j=0; j<count; ++j){
1677 if(var_deriv(variables[j])){
1678 fprintf(f, "%d %d %10.3g\n", i, integrator_ida_diffindex(blsys,variables[j]), derivatives[j]);
1679 }
1680 }
1681 }else if(type==II_WRITE_Y){
1682 for(j=0; j<count; ++j){
1683 if(!var_deriv(variables[j])){
1684 fprintf(f, "%d %d %10.3g\n", i, var_sindex(variables[j]), derivatives[j]);
1685 }
1686 }
1687 }
1688 }
1689 ASC_FREE(variables);
1690 ASC_FREE(derivatives);
1691 }
1692
1693 /**
1694 This routine outputs matrix structure in a crude text format, for the sake
1695 of debugging.
1696 */
1697 void integrator_ida_write_incidence(IntegratorSystem *blsys){
1698 int i, j;
1699 struct rel_relation **relptr;
1700 IntegratorIdaData *enginedata = blsys->enginedata;
1701 double *derivatives;
1702 struct var_variable **variables;
1703 int count, status;
1704 char *relname;
1705
1706 if(enginedata->nrels > 100){
1707 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
1708 return;
1709 }
1710
1711 variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);
1712 derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1713
1714 CONSOLE_DEBUG("Outputting incidence information to console...");
1715
1716 for(i=0, relptr = enginedata->rellist;
1717 i< enginedata->nrels && relptr != NULL;
1718 ++i, ++relptr
1719 ){
1720 relname = rel_make_name(blsys->system, *relptr);
1721
1722 /* get derivatives for this particular relation */
1723 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1724 if(status){
1725 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1726 ASC_FREE(relname);
1727 break;
1728 }
1729
1730 fprintf(stderr,"%3d:%-15s:",i,relname);
1731 ASC_FREE(relname);
1732
1733 for(j=0; j<count; ++j){
1734 if(var_deriv(variables[j])){
1735 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(blsys,variables[j]));
1736 }else{
1737 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
1738 }
1739 }
1740 fprintf(stderr,"\n");
1741 }
1742 ASC_FREE(variables);
1743 ASC_FREE(derivatives);
1744 }
1745
1746 /* @return 0 on success */
1747 int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
1748 char *varname, *relname;
1749 struct var_variable **vlist, *var;
1750 struct rel_relation **rlist, *rel;
1751 long vlen, rlen;
1752 long i;
1753
1754 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",sys->n_y);
1755
1756 /* if(integrator_sort_obs_vars(sys))return 10; */
1757
1758 if(sys->y && sys->ydot){
1759 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
1760 fprintf(fp,"index\ty\tydot\n");
1761 fprintf(fp,"-----\t-----\t-----\n");
1762 for(i=0;i<sys->n_y;++i){
1763 varname = var_make_name(sys->system, sys->y[i]);
1764 fprintf(fp,"%ld\t%s\t",i,varname);
1765 if(sys->ydot[i]){
1766 ASC_FREE(varname);
1767 varname = var_make_name(sys->system, sys->ydot[i]);
1768 fprintf(fp,"%s\n",varname);
1769 ASC_FREE(varname);
1770 }else{
1771 fprintf(fp,".\n");
1772 ASC_FREE(varname);
1773 }
1774 }
1775 }else{
1776 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
1777 }
1778
1779 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
1780 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
1781 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
1782
1783
1784 /* visit all the slv_system_t master var lists to collect vars */
1785 /* find the vars mostly in this one */
1786 vlist = slv_get_solvers_var_list(sys->system);
1787 vlen = slv_get_num_solvers_vars(sys->system);
1788 for(i=0;i<vlen;i++){
1789 var = vlist[i];
1790
1791 varname = var_make_name(sys->system, var);
1792 fprintf(fp,"%ld\t%-15s\t",i,varname);
1793
1794 if(var_fixed(var)){
1795 // it's fixed, so not really a DAE var
1796 fprintf(fp,"(fixed)\n");
1797 }else if(!var_active(var)){
1798 // inactive
1799 fprintf(fp,"(inactive)\n");
1800 }else if(!var_incident(var)){
1801 // not incident
1802 fprintf(fp,"(not incident)\n");
1803 }else{
1804 if(var_deriv(var)){
1805 if(sys->y_id){
1806 ASC_FREE(varname);
1807 varname = var_make_name(sys->system,vlist[integrator_ida_diffindex(sys,var)]);
1808 fprintf(fp,".\tdiff(%d='%s')\n",integrator_ida_diffindex(sys,var),varname);
1809 }else{
1810 fprintf(fp,".\tderiv... of??\n");
1811 }
1812 }else{
1813 fprintf(fp,"%d\t.\n",var_sindex(var));
1814 }
1815 }
1816 ASC_FREE(varname);
1817 }
1818
1819 /* let's write out the relations too */
1820 rlist = slv_get_solvers_rel_list(sys->system);
1821 rlen = slv_get_num_solvers_rels(sys->system);
1822
1823 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
1824 fprintf(fp,"index\tname\n");
1825 fprintf(fp,"-----\t----\n");
1826 for(i=0; i<rlen; ++i){
1827 rel = rlist[i];
1828 relname = rel_make_name(sys->system,rel);
1829 fprintf(fp,"%ld\t%s\n",i,relname);
1830 ASC_FREE(relname);
1831 }
1832
1833 /* write out the derivative chains */
1834 fprintf(fp,"\nDERIVATIVE CHAINS\n");
1835 if(integrator_ida_analyse_debug(sys,stderr)){
1836 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
1837 return 340;
1838 }
1839 fprintf(fp,"\n");
1840
1841 /* and lets write block debug output */
1842 system_block_debug(sys->system, fp);
1843
1844 return 0; /* success */
1845 }
1846
1847 /*----------------------------------------------
1848 ERROR REPORTING
1849 */
1850 /**
1851 Error message reporter function to be passed to IDA. All error messages
1852 will trigger a call to this function, so we should find everything
1853 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1854 panel (in the case of PyGTK).
1855 */
1856 void integrator_ida_error(int error_code
1857 , const char *module, const char *function
1858 , char *msg, void *eh_data
1859 ){
1860 IntegratorSystem *blsys;
1861 error_severity_t sev;
1862
1863 /* cast back the IntegratorSystem, just in case we need it */
1864 blsys = (IntegratorSystem *)eh_data;
1865
1866 /* severity depends on the sign of the error_code value */
1867 if(error_code <= 0){
1868 sev = ASC_PROG_ERR;
1869 }else{
1870 sev = ASC_PROG_WARNING;
1871 }
1872
1873 /* use our all-purpose error reporting to get stuff back to the GUI */
1874 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1875 }

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