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

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