/[ascend]/trunk/solvers/ida/ida.c
ViewVC logotype

Contents of /trunk/solvers/ida/ida.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2097 - (show annotations) (download) (as text)
Tue Dec 1 10:18:17 2009 UTC (13 years, 2 months ago) by jpye
File MIME type: text/x-csrc
File size: 71410 byte(s)
Move to 'Science' menu in Ubuntu.
Working to fix support for new IDA.
Fixing errors in building of pygtk/test*.cpp.
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 <signal.h>
37 #include <setjmp.h>
38 #include <fenv.h>
39 #include <math.h>
40
41 /* SUNDIALS includes */
42 #ifdef ASC_WITH_IDA
43
44 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
45 # include <sundials/sundials_config.h>
46 # include <sundials/sundials_nvector.h>
47 # include <ida/ida_spgmr.h>
48 # include <ida.h>
49 # include <nvector_serial.h>
50 #else
51 # include <sundials/sundials_config.h>
52 # include <nvector/nvector_serial.h>
53 # include <ida/ida.h>
54 #endif
55
56 # include <sundials/sundials_dense.h>
57 # include <ida/ida_spgmr.h>
58 # include <ida/ida_spbcgs.h>
59 # include <ida/ida_sptfqmr.h>
60 # include <ida/ida_dense.h>
61
62 # ifndef IDA_SUCCESS
63 # error "Failed to include SUNDIALS IDA header file"
64 # endif
65 #else
66 # error "If you're building this file, you should have ASC_WITH_IDA"
67 #endif
68
69 #ifdef ASC_WITH_MMIO
70 # include <mmio.h>
71 #endif
72
73 #include <ascend/utilities/ascConfig.h>
74 #include <ascend/utilities/error.h>
75 #include <ascend/utilities/ascSignal.h>
76 #include <ascend/utilities/ascPanic.h>
77 #include <ascend/compiler/instance_enum.h>
78
79 #include <ascend/system/slv_client.h>
80 #include <ascend/system/relman.h>
81 #include <ascend/system/block.h>
82 #include <ascend/system/slv_stdcalls.h>
83 #include <ascend/system/jacobian.h>
84 #include <ascend/system/bndman.h>
85
86 #include <ascend/utilities/config.h>
87 #include <ascend/integrator/integrator.h>
88
89 #include "idalinear.h"
90 #include "idaanalyse.h"
91 #include "ida_impl.h"
92
93 /*
94 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
95 */
96 #ifndef SUNDIALS_VERSION_MINOR
97 # ifdef __GNUC__
98 # warning "GUESSING SUNDIALS VERSION 2.2"
99 # endif
100 # define SUNDIALS_VERSION_MINOR 2
101 #endif
102 #ifndef SUNDIALS_VERSION_MAJOR
103 # define SUNDIALS_VERSION_MAJOR 2
104 #endif
105
106 /* SUNDIALS 2.4.0 introduces new DlsMat in place of DenseMat */
107 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==4
108 # define IDA_MTX_T DlsMat
109 # define IDADENSE_SUCCESS IDADLS_SUCCESS
110 # define IDADENSE_MEM_NULL IDADLS_MEM_NULL
111 # define IDADENSE_ILL_INPUT IDADLS_ILL_INPUT
112 # define IDADENSE_MEM_FAIL IDADLS_MEM_FAIL
113 #else
114 # define IDA_MTX_T DenseMat
115 #endif
116
117 /* #define FEX_DEBUG */
118 #define JEX_DEBUG
119 /* #define DJEX_DEBUG */
120 #define SOLVE_DEBUG
121 #define STATS_DEBUG
122 #define PREC_DEBUG
123 /* #define ROOT_DEBUG */
124
125 /* #define DIFFINDEX_DEBUG */
126 /* #define ANALYSE_DEBUG */
127 /* #define DESTROY_DEBUG */
128 /* #define MATRIX_DEBUG */
129
130 static IntegratorCreateFn integrator_ida_create;
131 static IntegratorParamsDefaultFn integrator_ida_params_default;
132 static IntegratorSolveFn integrator_ida_solve;
133 static IntegratorFreeFn integrator_ida_free;
134 IntegratorDebugFn integrator_ida_debug;
135 static IntegratorWriteMatrixFn integrator_ida_write_matrix;
136
137 /**
138 Everthing that the outside world needs to know about IDA
139 */
140 static const IntegratorInternals integrator_ida_internals = {
141 integrator_ida_create
142 ,integrator_ida_params_default
143 ,integrator_ida_analyse
144 ,integrator_ida_solve
145 ,integrator_ida_write_matrix
146 ,integrator_ida_debug
147 ,integrator_ida_free
148 ,INTEG_IDA
149 ,"IDA"
150 };
151
152 extern ASC_EXPORT int ida_register(void){
153 CONSOLE_DEBUG("Registering IDA...");
154 return integrator_register(&integrator_ida_internals);
155 }
156
157 /*-------------------------------------------------------------
158 FORWARD DECLS
159 */
160
161 /* forward dec needed for IntegratorIdaPrecFreeFn */
162 struct IntegratorIdaDataStruct;
163
164 /* functions for allocating storage for and freeing preconditioner data */
165 typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *sys);
166 typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);
167
168
169 /**
170 Struct containing any stuff that IDA needs that doesn't fit into the
171 common IntegratorSystem struct.
172 */
173 typedef struct IntegratorIdaDataStruct{
174
175 struct rel_relation **rellist; /**< NULL terminated list of ACTIVE rels */
176 int nrels; /* number of ACTIVE rels */
177
178 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
179 int nbnds; /* number of boundaries */
180
181 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
182 var_filter_t vfilter;
183 rel_filter_t rfilter; /**< Used to filter relations from solver's rellist (@TODO needs work) */
184 void *precdata; /**< For use by the preconditioner */
185 IntegratorIdaPrecFreeFn *pfree; /**< Store instructions here on how to free precdata */
186
187 } IntegratorIdaData;
188
189
190 typedef struct IntegratorIdaPrecDJStruct{
191 N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */
192 } IntegratorIdaPrecDataJacobi;
193
194 typedef struct IntegratorIdaPrecDJFStruct{
195 linsolqr_system_t L;
196 } IntegratorIdaPrecDataJacobian;
197
198 /**
199 Hold all the function pointers associated with a particular preconditioner
200 We don't need to store the 'pfree' function here as it is allocated to the enginedata struct
201 by the pcreate function (ensures that corresponding 'free' and 'create' are always used)
202
203 @note IDA uses a different convention for function pointer types, so no '*'.
204 */
205 typedef struct IntegratorIdaPrecStruct{
206 IntegratorIdaPrecCreateFn *pcreate;
207 IDASpilsPrecSetupFn psetup;
208 IDASpilsPrecSolveFn psolve;
209 } IntegratorIdaPrec;
210
211 /* residual function forward declaration */
212 static int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
213
214 static int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
215 , N_Vector v, N_Vector Jv, realtype c_j
216 , void *jac_data, N_Vector tmp1, N_Vector tmp2
217 );
218
219 /* error handler forward declaration */
220 static void integrator_ida_error(int error_code
221 , const char *module, const char *function
222 , char *msg, void *eh_data
223 );
224
225 /* dense jacobian evaluation for IDADense dense direct linear solver */
226 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
227 static int integrator_ida_djex(int Neq, realtype tt, realtype c_j
228 , N_Vector yy, N_Vector yp, N_Vector rr
229 , IDA_MTX_T Jac, void *jac_data
230 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
231 );
232 #else
233 static int integrator_ida_djex(long int Neq, realtype tt
234 , N_Vector yy, N_Vector yp, N_Vector rr
235 , realtype c_j, void *jac_data, IDA_MTX_T Jac
236 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
237 );
238 #endif
239
240 /* sparse jacobian evaluation for ASCEND's sparse direct solver */
241 static IntegratorSparseJacFn integrator_ida_sjex;
242
243 /* boundary-detection function */
244 static int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
245
246 typedef struct IntegratorIdaStatsStruct{
247 long nsteps;
248 long nrevals;
249 long nlinsetups;
250 long netfails;
251 int qlast, qcur;
252 realtype hinused, hlast, hcur;
253 realtype tcur;
254 } IntegratorIdaStats;
255
256 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
257
258 /*static IntegratorVarVisitorFn integrator_dae_classify_var;
259 static void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
260 static void integrator_dae_show_var(IntegratorSystem *sys, struct var_variable *var, const int *varindx); */
261
262 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
263 static void integrator_ida_write_stats(IntegratorIdaStats *stats);
264 static void integrator_ida_write_incidence(IntegratorSystem *sys);
265
266 /*------
267 Full jacobian preconditioner -- experimental
268 */
269
270 static int integrator_ida_psetup_jacobian(realtype tt,
271 N_Vector yy, N_Vector yp, N_Vector rr,
272 realtype c_j, void *prec_data,
273 N_Vector tmp1, N_Vector tmp2,
274 N_Vector tmp3
275 );
276
277 static int integrator_ida_psolve_jacobian(realtype tt,
278 N_Vector yy, N_Vector yp, N_Vector rr,
279 N_Vector rvec, N_Vector zvec,
280 realtype c_j, realtype delta, void *prec_data,
281 N_Vector tmp
282 );
283
284 static void integrator_ida_pcreate_jacobian(IntegratorSystem *sys);
285
286 static void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata);
287
288 static const IntegratorIdaPrec prec_jacobian = {
289 integrator_ida_pcreate_jacobian
290 , integrator_ida_psetup_jacobian
291 , integrator_ida_psolve_jacobian
292 };
293
294 /*------
295 Jacobi preconditioner -- experimental
296 */
297
298 static int integrator_ida_psetup_jacobi(realtype tt,
299 N_Vector yy, N_Vector yp, N_Vector rr,
300 realtype c_j, void *prec_data,
301 N_Vector tmp1, N_Vector tmp2,
302 N_Vector tmp3
303 );
304
305 static int integrator_ida_psolve_jacobi(realtype tt,
306 N_Vector yy, N_Vector yp, N_Vector rr,
307 N_Vector rvec, N_Vector zvec,
308 realtype c_j, realtype delta, void *prec_data,
309 N_Vector tmp
310 );
311
312 static void integrator_ida_pcreate_jacobi(IntegratorSystem *sys);
313
314 static void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);
315
316 static const IntegratorIdaPrec prec_jacobi = {
317 integrator_ida_pcreate_jacobi
318 , integrator_ida_psetup_jacobi
319 , integrator_ida_psolve_jacobi
320 };
321
322 /*-------------------------------------------------------------
323 SETUP/TEARDOWN ROUTINES
324 */
325 static void integrator_ida_create(IntegratorSystem *sys){
326 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
327 IntegratorIdaData *enginedata;
328 enginedata = ASC_NEW(IntegratorIdaData);
329 CONSOLE_DEBUG("enginedata = %p",enginedata);
330 enginedata->rellist = NULL;
331 enginedata->safeeval = 0;
332 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED;
333 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;
334 enginedata->pfree = NULL;
335
336 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
337 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
338
339 sys->enginedata = (void *)enginedata;
340
341 integrator_ida_params_default(sys);
342 }
343
344 static void integrator_ida_free(void *enginedata){
345 #ifdef DESTROY_DEBUG
346 CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);
347 #endif
348 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
349 asc_assert(d);
350 if(d->pfree){
351 CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);
352 /* free the preconditioner data, whatever it happens to be */
353 (d->pfree)(enginedata);
354 }
355
356 ASC_FREE(d->rellist);
357
358 #ifdef DESTROY_DEBUG
359 CONSOLE_DEBUG("Now destroying the enginedata");
360 #endif
361 ASC_FREE(d);
362 #ifdef DESTROY_DEBUG
363 CONSOLE_DEBUG("enginedata freed");
364 #endif
365 }
366
367 static IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *sys){
368 IntegratorIdaData *d;
369 assert(sys!=NULL);
370 assert(sys->enginedata!=NULL);
371 assert(sys->engine==INTEG_IDA);
372 d = ((IntegratorIdaData *)(sys->enginedata));
373 return d;
374 }
375
376 /*-------------------------------------------------------------
377 PARAMETERS FOR IDA
378 */
379
380 enum ida_parameters{
381 IDA_PARAM_LINSOLVER
382 ,IDA_PARAM_MAXL
383 ,IDA_PARAM_MAXORD
384 ,IDA_PARAM_AUTODIFF
385 ,IDA_PARAM_CALCIC
386 ,IDA_PARAM_SAFEEVAL
387 ,IDA_PARAM_RTOL
388 ,IDA_PARAM_ATOL
389 ,IDA_PARAM_ATOLVECT
390 ,IDA_PARAM_GSMODIFIED
391 ,IDA_PARAM_MAXNCF
392 ,IDA_PARAM_PREC
393 ,IDA_PARAMS_SIZE
394 };
395
396 /**
397 Here the full set of parameters is defined, along with upper/lower bounds,
398 etc. The values are stuck into the sys->params structure.
399
400 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
401 list. Then add a slv_param_*(...) statement below to define the type, description and range
402 for the new parameter.
403
404 @return 0 on success
405 */
406 static int integrator_ida_params_default(IntegratorSystem *sys){
407 asc_assert(sys!=NULL);
408 asc_assert(sys->engine==INTEG_IDA);
409 slv_parameters_t *p;
410 p = &(sys->params);
411
412 slv_destroy_parms(p);
413
414 if(p->parms==NULL){
415 CONSOLE_DEBUG("params NULL");
416 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
417 if(p->parms==NULL)return -1;
418 p->dynamic_parms = 1;
419 }else{
420 CONSOLE_DEBUG("params not NULL");
421 }
422
423 /* reset the number of parameters to zero so that we can check it at the end */
424 p->num_parms = 0;
425
426 slv_param_bool(p,IDA_PARAM_AUTODIFF
427 ,(SlvParameterInitBool){{"autodiff"
428 ,"Use auto-diff?",1
429 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
430 }, TRUE}
431 );
432
433 slv_param_char(p,IDA_PARAM_CALCIC
434 ,(SlvParameterInitChar){{"calcic"
435 ,"Initial conditions calcuation",1
436 ,"Use specified values of ydot to solve for inital y (Y),"
437 " or use the the values of the differential variables (yd) to solve"
438 " for the pure algebraic variables (ya) along with the derivatives"
439 " of the differential variables (yddot) (YA_YDP), or else don't solve"
440 " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
441 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
442 );
443
444 slv_param_bool(p,IDA_PARAM_SAFEEVAL
445 ,(SlvParameterInitBool){{"safeeval"
446 ,"Use safe evaluation?",1
447 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
448 "throw SIGFPE errors which will then halt integration."
449 }, FALSE}
450 );
451
452
453 slv_param_bool(p,IDA_PARAM_ATOLVECT
454 ,(SlvParameterInitBool){{"atolvect"
455 ,"Use 'ode_atol' values as specified?",1
456 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
457 " in the integration. If FALSE, a scalar absolute tolerance value"
458 " is shared by all variables. See IDA manual, section 5.5.1"
459 }, TRUE }
460 );
461
462 slv_param_real(p,IDA_PARAM_ATOL
463 ,(SlvParameterInitReal){{"atol"
464 ,"Scalar absolute error tolerance",1
465 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
466 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
467 }, 1e-5, 0, 1e10 }
468 );
469
470 slv_param_real(p,IDA_PARAM_RTOL
471 ,(SlvParameterInitReal){{"rtol"
472 ,"Scalar relative error tolerance",1
473 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
474 " it's not possible to set per-variable relative tolerances as it is"
475 " with LSODE)."
476 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
477 }, 1e-4, 0, 1 }
478 );
479
480 slv_param_char(p,IDA_PARAM_LINSOLVER
481 ,(SlvParameterInitChar){{"linsolver"
482 ,"Linear solver",1
483 ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"
484 " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"
485 " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"
486 " or SPTFQMR (which still need preconditioners to be implemented"
487 " before they can be very useful."
488 }, "DENSE"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
489 );
490
491 slv_param_int(p,IDA_PARAM_MAXL
492 ,(SlvParameterInitInt){{"maxl"
493 ,"Maximum Krylov dimension",0
494 ,"The maximum dimension of Krylov space used by the linear solver"
495 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
496 " The default of 0 results in IDA using its internal default, which"
497 " is currently a value of 5."
498 }, 0, 0, 20 }
499 );
500
501 slv_param_int(p,IDA_PARAM_MAXORD
502 ,(SlvParameterInitInt){{"maxord"
503 ,"Maximum order of linear multistep method",0
504 ,"The maximum order of the linear multistep method with IDA. See"
505 " IDA manual p 38."
506 }, 5, 1, 5 }
507 );
508
509 slv_param_bool(p,IDA_PARAM_GSMODIFIED
510 ,(SlvParameterInitBool){{"gsmodified"
511 ,"Gram-Schmidt Orthogonalisation Scheme", 2
512 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
513 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
514 }, TRUE}
515 );
516
517 slv_param_int(p,IDA_PARAM_MAXNCF
518 ,(SlvParameterInitInt){{"maxncf"
519 ,"Max nonlinear solver convergence failures per step", 2
520 ,"Maximum number of allowable nonlinear solver convergence failures"
521 " on one step. See IDA manual section 5.5.6.1."
522 }, 10,0,1000 }
523 );
524
525 slv_param_char(p,IDA_PARAM_PREC
526 ,(SlvParameterInitChar){{"prec"
527 ,"Preconditioner",1
528 ,"See IDA manual, section section 5.6.8."
529 },"NONE"}, (char *[]){"NONE","DIAG",NULL}
530 );
531
532 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
533
534 CONSOLE_DEBUG("Created %d params", p->num_parms);
535
536 return 0;
537 }
538
539 /*-------------------------------------------------------------
540 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
541 */
542
543 /*static double div1(double a, double b){
544 return a/b;
545 }*/
546
547 typedef int IdaFlagFn(void *,int *);
548 typedef char *IdaFlagNameFn(int);
549
550 /* return 0 on success */
551 static int integrator_ida_solve(
552 IntegratorSystem *sys
553 , unsigned long start_index
554 , unsigned long finish_index
555 ){
556 void *ida_mem;
557 int flag, flag1, t_index;
558 realtype t0, reltol, abstol, t, tret, tout1;
559 N_Vector y0, yp0, abstolvect, ypret, yret, id;
560 IntegratorIdaData *enginedata;
561 char *linsolver;
562 int maxl;
563 IdaFlagFn *flagfn;
564 IdaFlagNameFn *flagnamefn;
565 const char *flagfntype;
566 char *pname = NULL;
567 struct rel_relation **rels;
568 int *rootsfound;
569 char havecrossed;
570
571 #ifdef SOLVE_DEBUG
572 char *varname, *relname;
573 #endif
574
575 int i,j,n_activerels,n_solverrels;
576 const IntegratorIdaPrec *prec = NULL;
577 int icopt; /* initial conditions strategy */
578
579 CONSOLE_DEBUG("STARTING IDA...");
580
581 enginedata = integrator_ida_enginedata(sys);
582
583 enginedata->safeeval = SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_SAFEEVAL);
584 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
585
586 /* store reference to list of relations (in enginedata) */
587 n_solverrels = slv_get_num_solvers_rels(sys->system);
588
589 n_activerels = slv_count_solvers_rels(sys->system, &integrator_ida_rel);
590
591 enginedata->bndlist = slv_get_solvers_bnd_list(sys->system);
592 enginedata->nbnds = slv_get_num_solvers_bnds(sys->system);
593
594 enginedata->rellist = ASC_NEW_ARRAY(struct rel_relation *, n_activerels);
595
596 rels = slv_get_solvers_rel_list(sys->system);
597
598 j=0;
599 for(i=0; i < n_solverrels; ++i){
600 if(rel_apply_filter(rels[i], &integrator_ida_rel)){
601 #ifdef SOLVE_DEBUG
602 relname = rel_make_name(sys->system, rels[i]);
603 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));
604 ASC_FREE(relname);
605 #endif
606 enginedata->rellist[j++]=rels[i];
607 }
608 }
609 asc_assert(j==n_activerels);
610
611 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
612 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
613
614 CONSOLE_DEBUG("Number of relations: %d",n_solverrels);
615 CONSOLE_DEBUG("Number of active relations: %d",n_activerels);
616 CONSOLE_DEBUG("Number of dependent vars: %d",sys->n_y);
617 CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);
618
619 enginedata->nrels = n_activerels;
620
621 if(enginedata->nrels != sys->n_y){
622 ERROR_REPORTER_HERE(ASC_USER_ERROR
623 ,"Integration problem is not square (%d active rels, %d vars)"
624 ,n_activerels, sys->n_y
625 );
626 return 1; /* failure */
627 }
628
629 #ifdef SOLVE_DEBUG
630 integrator_ida_debug(sys,stderr);
631 #endif
632
633 /* retrieve initial values from the system */
634
635 /** @TODO fix this, the starting time != first sample */
636 t0 = integrator_get_t(sys);
637 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
638
639 CONSOLE_DEBUG("RETRIEVING y0");
640
641 y0 = N_VNew_Serial(sys->n_y);
642 integrator_get_y(sys,NV_DATA_S(y0));
643
644 #ifdef SOLVE_DEBUG
645 CONSOLE_DEBUG("RETRIEVING yp0");
646 #endif
647
648 yp0 = N_VNew_Serial(sys->n_y);
649 integrator_get_ydot(sys,NV_DATA_S(yp0));
650
651 #ifdef SOLVE_DEBUG
652 N_VPrint_Serial(yp0);
653 CONSOLE_DEBUG("yp0 is at %p",&yp0);
654 #endif
655
656 /* create IDA object */
657 ida_mem = IDACreate();
658
659 /* relative error tolerance */
660 reltol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_RTOL);
661 CONSOLE_DEBUG("rtol = %8.2e",reltol);
662
663
664 /* allocate internal memory */
665 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
666 flag = IDAInit(ida_mem, &integrator_ida_fex, t0, y0 ,yp0);
667 #else
668 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
669 /* vector of absolute tolerances */
670 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
671 abstolvect = N_VNew_Serial(sys->n_y);
672 integrator_get_atol(sys,NV_DATA_S(abstolvect));
673
674 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
675
676 N_VDestroy_Serial(abstolvect);
677 }else{
678 /* scalar absolute tolerance (one value for all) */
679 abstol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_ATOL);
680 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
681 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
682 }
683 #endif
684
685 if(flag==IDA_MEM_NULL){
686 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
687 return 2;
688 }else if(flag==IDA_MEM_FAIL){
689 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
690 return 3;
691 }else if(flag==IDA_ILL_INPUT){
692 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
693 return 4;
694 }/* else success */
695
696
697 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
698 CONSOLE_DEBUG("Assigning tolerances...");
699 /* assign tolerances */
700 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
701 CONSOLE_DEBUG("using vector of atol values");
702 abstolvect = N_VNew_Serial(sys->n_y);
703 integrator_get_atol(sys,NV_DATA_S(abstolvect));
704 IDASVtolerances(ida_mem, reltol, abstolvect);
705 N_VDestroy_Serial(abstolvect);
706 }else{
707 /* scalar tolerances */
708 abstol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_ATOL);
709 CONSOLE_DEBUG("using scalar atol value = %8.2e",abstol);
710 IDASStolerances(ida_mem, reltol, abstol);
711 }
712 #endif
713
714 /* set optional inputs... */
715 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)sys);
716 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
717 IDASetUserData(ida_mem, (void *)sys);
718 #else
719 IDASetRdata(ida_mem, (void *)sys);
720 #endif
721 IDASetMaxStep(ida_mem, integrator_get_maxstep(sys));
722 IDASetInitStep(ida_mem, integrator_get_stepzero(sys));
723 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(sys));
724 if(integrator_get_minstep(sys)>0){
725 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
726 }
727
728 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXNCF));
729 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXNCF));
730
731 CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXORD));
732 IDASetMaxOrd(ida_mem,SLV_PARAM_INT(&sys->params,IDA_PARAM_MAXORD));
733
734 /* there's no capability for setting *minimum* step size in IDA */
735
736
737 /* attach linear solver module, using the default value of maxl */
738 linsolver = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_LINSOLVER);
739 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
740 if(strcmp(linsolver,"ASCEND")==0){
741 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",sys->n_y);
742 IDAASCEND(ida_mem,sys->n_y);
743 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)sys);
744
745 flagfntype = "IDAASCEND";
746 flagfn = &IDAASCENDGetLastFlag;
747 flagnamefn = &IDAASCENDGetReturnFlagName;
748
749 }else if(strcmp(linsolver,"DENSE")==0){
750 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",sys->n_y);
751 flag = IDADense(ida_mem, sys->n_y);
752 switch(flag){
753 case IDADENSE_SUCCESS: break;
754 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 5;
755 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 5;
756 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 5;
757 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
758 }
759
760 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
761 CONSOLE_DEBUG("USING AUTODIFF");
762 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
763 flag = IDADlsSetDenseJacFn(ida_mem, &integrator_ida_djex);
764 #else
765 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)sys);
766 #endif
767 switch(flag){
768 case IDADENSE_SUCCESS: break;
769 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
770 }
771 }else{
772 CONSOLE_DEBUG("USING NUMERICAL DIFF");
773 }
774
775 flagfntype = "IDADENSE";
776 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
777 flagfn = &IDADlsGetLastFlag;
778 flagnamefn = &IDADlsGetReturnFlagName;
779 #else
780 flagfn = &IDADenseGetLastFlag;
781 flagnamefn = &IDADenseGetReturnFlagName;
782 #endif
783 }else{
784 /* remaining methods are all SPILS */
785 CONSOLE_DEBUG("IDA SPILS");
786
787 maxl = SLV_PARAM_INT(&(sys->params),IDA_PARAM_MAXL);
788 CONSOLE_DEBUG("maxl = %d",maxl);
789
790 /* what preconditioner? */
791 pname = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_PREC);
792 if(strcmp(pname,"NONE")==0){
793 prec = NULL;
794 }else if(strcmp(pname,"JACOBI")==0){
795 prec = &prec_jacobi;
796 }else{
797 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
798 return 7;
799 }
800
801 /* which SPILS linear solver? */
802 if(strcmp(linsolver,"SPGMR")==0){
803 CONSOLE_DEBUG("IDA SPGMR");
804 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
805 }else if(strcmp(linsolver,"SPBCG")==0){
806 CONSOLE_DEBUG("IDA SPBCG");
807 flag = IDASpbcg(ida_mem, maxl);
808 }else if(strcmp(linsolver,"SPTFQMR")==0){
809 CONSOLE_DEBUG("IDA SPTFQMR");
810 flag = IDASptfqmr(ida_mem,maxl);
811 }else{
812 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
813 return 8;
814 }
815
816 if(prec){
817 /* assign the preconditioner to the linear solver */
818 (prec->pcreate)(sys);
819 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
820 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve);
821 #else
822 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)sys);
823 #endif
824 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
825 }else{
826 CONSOLE_DEBUG("No preconditioner");
827 }
828
829 flagfntype = "IDASPILS";
830 flagfn = &IDASpilsGetLastFlag;
831 flagnamefn = &IDASpilsGetReturnFlagName;
832
833 if(flag==IDASPILS_MEM_NULL){
834 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
835 return 9;
836 }else if(flag==IDASPILS_MEM_FAIL){
837 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
838 return 9;
839 }/* else success */
840
841 /* assign the J*v function */
842 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
843 CONSOLE_DEBUG("USING AUTODIFF");
844 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
845 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex);
846 #else
847 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)sys);
848 #endif
849 if(flag==IDASPILS_MEM_NULL){
850 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
851 return 10;
852 }else if(flag==IDASPILS_LMEM_NULL){
853 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
854 return 10;
855 }/* else success */
856 }else{
857 CONSOLE_DEBUG("USING NUMERICAL DIFF");
858 }
859
860 if(strcmp(linsolver,"SPGMR")==0){
861 /* select Gram-Schmidt orthogonalisation */
862 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_GSMODIFIED)){
863 CONSOLE_DEBUG("USING MODIFIED GS");
864 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
865 if(flag!=IDASPILS_SUCCESS){
866 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
867 return 11;
868 }
869 }else{
870 CONSOLE_DEBUG("USING CLASSICAL GS");
871 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
872 if(flag!=IDASPILS_SUCCESS){
873 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
874 return 11;
875 }
876 }
877 }
878 }
879
880 /* set linear solver optional inputs...
881 ...nothing here at the moment...
882 */
883
884 /* calculate initial conditions */
885 icopt = 0;
886 if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"Y")==0){
887 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
888 icopt = IDA_Y_INIT;
889 asc_assert(icopt!=0);
890 }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
891 CONSOLE_DEBUG("Solving initial conditions using values of yd");
892 icopt = IDA_YA_YDP_INIT;
893 asc_assert(icopt!=0);
894 id = N_VNew_Serial(sys->n_y);
895 for(i=0; i < sys->n_y; ++i){
896 if(sys->ydot[i] == NULL){
897 NV_Ith_S(id,i) = 0.0;
898 #ifdef SOLVE_DEBUG
899 varname = var_make_name(sys->system,sys->y[i]);
900 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
901 ASC_FREE(varname);
902 #endif
903 }else{
904 #ifdef SOLVE_DEBUG
905 CONSOLE_DEBUG("y[%d] is differential",i);
906 #endif
907 NV_Ith_S(id,i) = 1.0;
908 }
909 }
910 IDASetId(ida_mem, id);
911 N_VDestroy_Serial(id);
912 }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"NONE")==0){
913 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
914 }else{
915 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
916 }
917
918 if(icopt){
919 sys->currentstep=0;
920 t_index=start_index + 1;
921 tout1 = samplelist_get(sys->samples, t_index);
922
923 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
924
925 #ifdef ASC_SIGNAL_TRAPS
926 /* catch SIGFPE if desired to */
927 if(enginedata->safeeval){
928 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
929 Asc_SignalHandlerPush(SIGFPE,SIG_DFL);
930 }else{
931 # ifdef FEX_DEBUG
932 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
933 # endif
934 Asc_SignalHandlerPushDefault(SIGFPE);
935 }
936 if (setjmp(g_fpe_env)==0) {
937 #endif
938
939 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=3
940 flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3 */
941 # else
942 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
943 # endif
944 /* check flags and output status */
945 switch(flag){
946 case IDA_SUCCESS:
947 CONSOLE_DEBUG("Initial conditions solved OK");
948 break;
949
950 case IDA_LSETUP_FAIL:
951 case IDA_LINIT_FAIL:
952 case IDA_LSOLVE_FAIL:
953 case IDA_NO_RECOVERY:
954 flag1 = -999;
955 flag = (flagfn)(ida_mem,&flag1);
956 if(flag){
957 ERROR_REPORTER_HERE(ASC_PROG_ERR
958 ,"Unable to retrieve error code from %s (err %d)"
959 ,flagfntype,flag
960 );
961 return 12;
962 }
963 ERROR_REPORTER_HERE(ASC_PROG_ERR
964 ,"%s returned flag '%s' (value = %d)"
965 ,flagfntype,(flagnamefn)(flag1),flag1
966 );
967 return 12;
968
969 default:
970 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
971 return 12;
972 }
973 #ifdef ASC_SIGNAL_TRAPS
974 }else{
975 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
976 return 13;
977 }
978
979 if(enginedata->safeeval){
980 Asc_SignalHandlerPop(SIGFPE,SIG_DFL);
981 }else{
982 CONSOLE_DEBUG("pop...");
983 Asc_SignalHandlerPopDefault(SIGFPE);
984 CONSOLE_DEBUG("...pop");
985 }
986 #endif
987 }/* icopt */
988
989 /* optionally, specify ROO-FINDING problem */
990 if(enginedata->nbnds){
991 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
992 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn);
993 #else
994 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)sys);
995 #endif
996 }
997
998 /* -- set up the IntegratorReporter */
999 integrator_output_init(sys);
1000
1001 /* -- store the initial values of all the stuff */
1002 integrator_output_write(sys);
1003 integrator_output_write_obs(sys);
1004
1005 /* specify where the returned values should be stored */
1006 yret = y0;
1007 ypret = yp0;
1008
1009 /* advance solution in time, return values as yret and derivatives as ypret */
1010 sys->currentstep=1;
1011 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++sys->currentstep){
1012 t = samplelist_get(sys->samples, t_index);
1013 t0 = integrator_get_t(sys);
1014 asc_assert(t > t0);
1015
1016 #ifdef SOLVE_DEBUG
1017 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
1018 #endif
1019
1020 #ifdef ASC_SIGNAL_TRAPS
1021 Asc_SignalHandlerPushDefault(SIGINT);
1022 if(setjmp(g_int_env)==0) {
1023 #endif
1024 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
1025 #ifdef ASC_SIGNAL_TRAPS
1026 }else{
1027 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Caught interrupt");
1028 flag = -555;
1029 }
1030 Asc_SignalHandlerPopDefault(SIGINT);
1031 #endif
1032
1033 /* check for roots found */
1034 if(enginedata->nbnds){
1035 rootsfound = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
1036 havecrossed = 0;
1037 if(IDA_SUCCESS == IDAGetRootInfo(ida_mem, rootsfound)){
1038 for(i=0; i < enginedata->nbnds; ++i){
1039 if(rootsfound[i]){
1040 havecrossed = 1;
1041 #ifdef SOLVE_DEBUG
1042 relname = bnd_make_name(sys->system,enginedata->bndlist[i]);
1043 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary '%s' crossed",relname);
1044 ASC_FREE(relname);
1045 #else
1046 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary %d crossed!", i);
1047 #endif
1048 }
1049 }
1050
1051 /* so, now we need to restart the integration. we will assume that
1052 everything changes: number of variables, etc, etc, etc. */
1053
1054 if(havecrossed){
1055 CONSOLE_DEBUG("Boundaries were crossed; need to reinitialise solver...");
1056 /** try resetting the boundary states now? */
1057 //IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)sys);
1058
1059 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
1060 IDAReInit(ida_mem, tret, yret, ypret);
1061 #elif SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR<3
1062 /* TODO find out what version needs the following line... not sure. */
1063 IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret);
1064 #else
1065 /* allocate internal memory */
1066 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
1067 /* vector of absolute tolerances */
1068 abstolvect = N_VNew_Serial(sys->n_y);
1069 integrator_get_atol(sys,NV_DATA_S(abstolvect));
1070 if(IDA_SUCCESS != IDAReInit(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect)){
1071 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
1072 }
1073 N_VDestroy_Serial(abstolvect);
1074 }else{
1075 /* scalar absolute tolerance (one value for all) */
1076 abstol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_ATOL);
1077 if(IDA_SUCCESS != IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol)){
1078 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
1079 }
1080 }
1081 #endif
1082 }
1083 }else{
1084 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");
1085 }
1086 ASC_FREE(rootsfound);
1087 }
1088
1089
1090 /* pass the values of everything back to the compiler */
1091 integrator_set_t(sys, (double)tret);
1092 integrator_set_y(sys, NV_DATA_S(yret));
1093 integrator_set_ydot(sys, NV_DATA_S(ypret));
1094
1095 if(flag<0){
1096 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
1097 break;
1098 }
1099
1100 /* -- do something so that sys knows the values of tret, yret and ypret */
1101
1102 /* -- store the current values of all the stuff */
1103 integrator_output_write(sys);
1104 integrator_output_write_obs(sys);
1105
1106 }/* loop through next sample timestep */
1107
1108 /* -- close the IntegratorReporter */
1109 integrator_output_close(sys);
1110
1111 /* get optional outputs */
1112 #ifdef STATS_DEBUG
1113 IntegratorIdaStats stats;
1114 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
1115 integrator_ida_write_stats(&stats);
1116 }else{
1117 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
1118 }
1119 #endif
1120
1121 /* free solution memory */
1122 N_VDestroy_Serial(yret);
1123 N_VDestroy_Serial(ypret);
1124
1125 /* free solver memory */
1126 IDAFree(ida_mem);
1127
1128 if(flag < -500){
1129 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", t);
1130 return -flag;
1131 }
1132
1133 if(flag < 0){
1134 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1135 return 14;
1136 }
1137
1138 /* all done, success */
1139 return 0;
1140 }
1141
1142 /*--------------------------------------------------
1143 RESIDUALS AND JACOBIAN AND IDAROOTFN
1144 */
1145
1146 #if 0
1147 typedef void (SignalHandlerFn)(int);
1148 SignalHandlerFn integrator_ida_sig;
1149 SignalHandlerFn *integrator_ida_sig_old;
1150 jmp_buf integrator_ida_jmp_buf;
1151 fenv_t integrator_ida_fenv_old;
1152
1153
1154 void integrator_ida_write_feinfo(){
1155 int f;
1156 f = fegetexcept();
1157 CONSOLE_DEBUG("Locating nature of exception...");
1158 if(f & FE_DIVBYZERO)ERROR_REPORTER_HERE(ASC_PROG_ERR,"DIV BY ZERO");
1159 if(f & FE_INEXACT)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INEXACT");
1160 if(f & FE_INVALID)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INVALID");
1161 if(f & FE_OVERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"OVERFLOW");
1162 if(f & FE_UNDERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"UNDERFLOW");
1163 if(f==0)ERROR_REPORTER_HERE(ASC_PROG_ERR,"FLAGS ARE CLEAR?!?");
1164 }
1165
1166 void integrator_ida_sig(int sig){
1167 /* the wrong signal: rethrow to the default handler */
1168 if(sig!=SIGFPE){
1169 signal(SIGFPE,SIG_DFL);
1170 raise(sig);
1171 }
1172 integrator_ida_write_feinfo();
1173 CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);
1174 longjmp(integrator_ida_jmp_buf,sig);
1175 }
1176 #endif
1177
1178 /**
1179 Function to evaluate system residuals, in the form required for IDA.
1180
1181 Given tt, yy and yp, we need to evaluate and return rr.
1182
1183 @param tt current value of indep variable (time)
1184 @param yy current values of dependent variable vector
1185 @param yp current values of derivatives of dependent variables
1186 @param rr the output residual vector (is we're returning data to)
1187 @param res_data pointer to our stuff (sys in this case).
1188
1189 @return 0 on success, positive on recoverable error, and
1190 negative on unrecoverable error.
1191 */
1192 static int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1193 IntegratorSystem *sys;
1194 IntegratorIdaData *enginedata;
1195 int i, calc_ok, is_error;
1196 struct rel_relation** relptr;
1197 double resid;
1198 char *relname;
1199 #ifdef FEX_DEBUG
1200 char *varname;
1201 char diffname[30];
1202 #endif
1203
1204 sys = (IntegratorSystem *)res_data;
1205 enginedata = integrator_ida_enginedata(sys);
1206
1207 #ifdef FEX_DEBUG
1208 /* fprintf(stderr,"\n\n"); */
1209 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1210 #endif
1211
1212 if(NV_LENGTH_S(rr)!=enginedata->nrels){
1213 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
1214 return -1; /* unrecoverable */
1215 }
1216
1217 /* pass the values of everything back to the compiler */
1218 integrator_set_t(sys, (double)tt);
1219 integrator_set_y(sys, NV_DATA_S(yy));
1220 integrator_set_ydot(sys, NV_DATA_S(yp));
1221
1222 /* perform bounds checking on all variables */
1223 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1224 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1225 return 1;
1226 }
1227
1228 /* evaluate each residual in the rellist */
1229 is_error = 0;
1230 relptr = enginedata->rellist;
1231
1232
1233 #ifdef ASC_SIGNAL_TRAPS
1234 if(enginedata->safeeval){
1235 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1236 }else{
1237 # ifdef FEX_DEBUG
1238 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1239 # endif
1240 Asc_SignalHandlerPushDefault(SIGFPE);
1241 }
1242
1243 if (SETJMP(g_fpe_env)==0) {
1244 #endif
1245
1246
1247 for(i=0, relptr = enginedata->rellist;
1248 i< enginedata->nrels && relptr != NULL;
1249 ++i, ++relptr
1250 ){
1251 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
1252
1253 NV_Ith_S(rr,i) = resid;
1254 if(!calc_ok){
1255 relname = rel_make_name(sys->system, *relptr);
1256 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1257 ASC_FREE(relname);
1258 /* presumable some output already made? */
1259 is_error = 1;
1260 }/*else{
1261 CONSOLE_DEBUG("Calc OK");
1262 }*/
1263 }
1264
1265 if(!is_error){
1266 for(i=0;i< enginedata->nrels; ++i){
1267 if(isnan(NV_Ith_S(rr,i))){
1268 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in residual %d",i);
1269 is_error=1;
1270 }
1271 }
1272 #ifdef FEX_DEBUG
1273 if(!is_error){
1274 CONSOLE_DEBUG("No NAN detected");
1275 }
1276 #endif
1277 }
1278
1279 #ifdef ASC_SIGNAL_TRAPS
1280 }else{
1281 relname = rel_make_name(sys->system, *relptr);
1282 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1283 ASC_FREE(relname);
1284 is_error = 1;
1285 }
1286
1287 if(enginedata->safeeval){
1288 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1289 }else{
1290 Asc_SignalHandlerPopDefault(SIGFPE);
1291 }
1292 #endif
1293
1294
1295 #ifdef FEX_DEBUG
1296 /* output residuals to console */
1297 CONSOLE_DEBUG("RESIDUAL OUTPUT");
1298 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1299 for(i=0; i<sys->n_y; ++i){
1300 varname = var_make_name(sys->system,sys->y[i]);
1301 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1302 if(sys->ydot[i]){
1303 varname = var_make_name(sys->system,sys->ydot[i]);
1304 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1305 }else{
1306 snprintf(diffname,99,"diff(%s)",varname);
1307 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1308 }
1309 ASC_FREE(varname);
1310 relname = rel_make_name(sys->system,enginedata->rellist[i]);
1311 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1312 }
1313 #endif
1314
1315 if(is_error){
1316 return 1;
1317 }
1318
1319 #ifdef FEX_DEBUG
1320 CONSOLE_DEBUG("RESIDUAL OK");
1321 #endif
1322 return 0;
1323 }
1324
1325 /**
1326 Dense Jacobian evaluation. Only suitable for small problems!
1327 Has been seen working for problems up to around 2000 vars, FWIW.
1328 */
1329 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
1330 static int integrator_ida_djex(int Neq, realtype tt, realtype c_j
1331 , N_Vector yy, N_Vector yp, N_Vector rr
1332 , IDA_MTX_T Jac, void *jac_data
1333 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1334 ){
1335 #else
1336 static int integrator_ida_djex(long int Neq, realtype tt
1337 , N_Vector yy, N_Vector yp, N_Vector rr
1338 , realtype c_j, void *jac_data, IDA_MTX_T Jac
1339 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1340 ){
1341 #endif
1342 IntegratorSystem *sys;
1343 IntegratorIdaData *enginedata;
1344 char *relname;
1345 #ifdef DJEX_DEBUG
1346 struct var_variable **varlist;
1347 char *varname;
1348 #endif
1349 struct rel_relation **relptr;
1350 int i;
1351 double *derivatives;
1352 struct var_variable **variables;
1353 int count, j;
1354 int status, is_error = 0;
1355
1356 sys = (IntegratorSystem *)jac_data;
1357 enginedata = integrator_ida_enginedata(sys);
1358
1359 /* allocate space for returns from relman_diff3 */
1360 /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
1361 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1362 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1363
1364 /* pass the values of everything back to the compiler */
1365 integrator_set_t(sys, (double)tt);
1366 integrator_set_y(sys, NV_DATA_S(yy));
1367 integrator_set_ydot(sys, NV_DATA_S(yp));
1368
1369 /* perform bounds checking on all variables */
1370 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1371 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1372 return 1;
1373 }
1374
1375 #ifdef DJEX_DEBUG
1376 varlist = slv_get_solvers_var_list(sys->system);
1377
1378 /* print vars */
1379 for(i=0; i < sys->n_y; ++i){
1380 varname = var_make_name(sys->system, sys->y[i]);
1381 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1382 asc_assert(NV_Ith_S(yy,i) == var_value(sys->y[i]));
1383 ASC_FREE(varname);
1384 }
1385
1386 /* print derivatives */
1387 for(i=0; i < sys->n_y; ++i){
1388 if(sys->ydot[i]){
1389 varname = var_make_name(sys->system, sys->ydot[i]);
1390 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(sys->ydot[i]));
1391 ASC_FREE(varname);
1392 }else{
1393 varname = var_make_name(sys->system, sys->y[i]);
1394 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1395 ASC_FREE(varname);
1396 }
1397 }
1398
1399 /* print step size */
1400 CONSOLE_DEBUG("<c_j> = %g",c_j);
1401 #endif
1402
1403 /* build up the dense jacobian matrix... */
1404 status = 0;
1405 for(i=0, relptr = enginedata->rellist;
1406 i< enginedata->nrels && relptr != NULL;
1407 ++i, ++relptr
1408 ){
1409 /* get derivatives for this particular relation */
1410 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1411
1412 if(status){
1413 relname = rel_make_name(sys->system, *relptr);
1414 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1415 ASC_FREE(relname);
1416 is_error = 1;
1417 break;
1418 }
1419
1420 /* output what's going on here ... */
1421 #ifdef DJEX_DEBUG
1422 relname = rel_make_name(sys->system, *relptr);
1423 fprintf(stderr,"%d: '%s': ",i,relname);
1424 for(j=0;j<count;++j){
1425 varname = var_make_name(sys->system, variables[j]);
1426 if(var_deriv(variables[j])){
1427 fprintf(stderr," '%s'=",varname);
1428 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(sys,variables[j]));
1429 }else{
1430 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1431 }
1432 ASC_FREE(varname);
1433 }
1434 /* relname is freed further down */
1435 fprintf(stderr,"\n");
1436 #endif
1437
1438 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1439 for(j=0; j < count; ++j){
1440 #ifdef DJEX_DEBUG
1441 varname = var_make_name(sys->system,variables[j]);
1442 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1443 ASC_FREE(varname);
1444 #endif
1445 if(!var_deriv(variables[j])){
1446 #ifdef DJEX_DEBUG
1447 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
1448 asc_assert(var_sindex(variables[j]) >= 0);
1449 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1450 #endif
1451 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1452 }else{
1453 DENSE_ELEM(Jac,i,integrator_ida_diffindex(sys,variables[j])) += derivatives[j] * c_j;
1454 #ifdef DJEX_DEBUG
1455 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1456 #endif
1457 }
1458 }
1459 }
1460
1461 #ifdef DJEX_DEBUG
1462 ASC_FREE(relname);
1463 CONSOLE_DEBUG("PRINTING JAC");
1464 fprintf(stderr,"\t");
1465 for(j=0; j < sys->n_y; ++j){
1466 if(j)fprintf(stderr,"\t");
1467 varname = var_make_name(sys->system,sys->y[j]);
1468 fprintf(stderr,"%11s",varname);
1469 ASC_FREE(varname);
1470 }
1471 fprintf(stderr,"\n");
1472 for(i=0; i < enginedata->nrels; ++i){
1473 relname = rel_make_name(sys->system, enginedata->rellist[i]);
1474 fprintf(stderr,"%s\t",relname);
1475 ASC_FREE(relname);
1476
1477 for(j=0; j < sys->n_y; ++j){
1478 if(j)fprintf(stderr,"\t");
1479 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1480 }
1481 fprintf(stderr,"\n");
1482 }
1483 #endif
1484
1485 /* test for NANs */
1486 if(!is_error){
1487 for(i=0;i< enginedata->nrels; ++i){
1488 for(j=0;j<sys->n_y;++j){
1489 if(isnan(DENSE_ELEM(Jac,i,j))){
1490 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
1491 is_error=1;
1492 }
1493 }
1494 }
1495 #ifdef DJEX_DEBUG
1496 if(!is_error){
1497 CONSOLE_DEBUG("No NAN detected");
1498 }
1499 #endif
1500 }
1501
1502 /* if(integrator_ida_check_diffindex(sys)){
1503 is_error = 1;
1504 }*/
1505
1506 if(is_error){
1507 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1508 return 1;
1509 }
1510
1511 #ifdef DJEX_DEBUG
1512 CONSOLE_DEBUG("DJEX RETURNING 0");
1513 /* ASC_PANIC("Quitting"); */
1514 #endif
1515 return 0;
1516 }
1517
1518 /**
1519 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1520
1521 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1522
1523 @param tt current value of the independent variable (time, t)
1524 @param yy current value of the dependent variable vector, y(t).
1525 @param yp current value of y'(t).
1526 @param rr current value of the residual vector F(t, y, y').
1527 @param v the vector by which the Jacobian must be multiplied to the right.
1528 @param Jv the output vector computed
1529 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1530 @param jac_data pointer to our stuff (sys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1531 @param tmp1 @see tmp2
1532 @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.
1533 @return 0 on success
1534 */
1535 static int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1536 , N_Vector v, N_Vector Jv, realtype c_j
1537 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1538 ){
1539 IntegratorSystem *sys;
1540 IntegratorIdaData *enginedata;
1541 int i, j, is_error=0;
1542 struct rel_relation** relptr = 0;
1543 char *relname;
1544 int status;
1545 double Jv_i;
1546
1547 struct var_variable **variables;
1548 double *derivatives;
1549 int count;
1550 struct var_variable **varlist;
1551 #ifdef JEX_DEBUG
1552
1553 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1554 #endif
1555
1556 sys = (IntegratorSystem *)jac_data;
1557 enginedata = integrator_ida_enginedata(sys);
1558 varlist = slv_get_solvers_var_list(sys->system);
1559
1560 /* pass the values of everything back to the compiler */
1561 integrator_set_t(sys, (double)tt);
1562 integrator_set_y(sys, NV_DATA_S(yy));
1563 integrator_set_ydot(sys, NV_DATA_S(yp));
1564 /* no real use for residuals (rr) here, I don't think? */
1565
1566 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1567
1568 i = NV_LENGTH_S(yy) * 2;
1569 #ifdef JEX_DEBUG
1570 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1571 #endif
1572 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1573 derivatives = ASC_NEW_ARRAY(double, i);
1574
1575 /* evaluate the derivatives... */
1576 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1577
1578 #ifdef ASC_SIGNAL_TRAPS
1579 Asc_SignalHandlerPushDefault(SIGFPE);
1580 if (SETJMP(g_fpe_env)==0) {
1581 #endif
1582 for(i=0, relptr = enginedata->rellist;
1583 i< enginedata->nrels && relptr != NULL;
1584 ++i, ++relptr
1585 ){
1586 /* get derivatives for this particular relation */
1587 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1588 #ifdef JEX_DEBUG
1589 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1590 #endif
1591
1592 if(status){
1593 relname = rel_make_name(sys->system, *relptr);
1594 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1595 ASC_FREE(relname);
1596 is_error = 1;
1597 break;
1598 }
1599
1600 /*
1601 Now we have the derivatives wrt each alg/diff variable in the
1602 present equation. variables[] points into the varlist. need
1603 a mapping from the varlist to the y and ydot lists.
1604 */
1605
1606 Jv_i = 0;
1607 for(j=0; j < count; ++j){
1608 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], sys->n_y);
1609 varname = var_make_name(sys->system, enginedata->varlist[variables[j]]);
1610 if(varname){
1611 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1612 ASC_FREE(varname);
1613 }else{
1614 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1615 }
1616 */
1617
1618 /* we don't calculate derivatives wrt indep var */
1619 asc_assert(variables[j]>=0);
1620 if(variables[j] == sys->x) continue;
1621 #ifdef JEX_DEBUG
1622 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1623 #endif
1624 if(var_deriv(variables[j])){
1625 #define DIFFINDEX integrator_ida_diffindex(sys,variables[j])
1626 #ifdef JEX_DEBUG
1627 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
1628 , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
1629 , i, DIFFINDEX, derivatives[j]
1630 , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
1631 );
1632 #endif
1633 asc_assert(sys->ydot[DIFFINDEX]==variables[j]);
1634 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1635 #undef DIFFINDEX
1636 }else{
1637 #define VARINDEX var_sindex(variables[j])
1638 #ifdef JEX_DEBUG
1639 asc_assert(sys->y[VARINDEX]==variables[j]);
1640 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
1641 , i
1642 , derivatives[j] * NV_Ith_S(v,VARINDEX)
1643 , i, VARINDEX, derivatives[j]
1644 , VARINDEX, NV_Ith_S(v,VARINDEX)
1645 );
1646 #endif
1647 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1648 #undef VARINDEX
1649 }
1650 }
1651
1652 NV_Ith_S(Jv,i) = Jv_i;
1653 #ifdef JEX_DEBUG
1654 CONSOLE_DEBUG("rel = %p",*relptr);
1655 relname = rel_make_name(sys->system, *relptr);
1656 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1657 ASC_FREE(relname);
1658 return 1;
1659 #endif
1660 }
1661 #ifdef ASC_SIGNAL_TRAPS
1662 }else{
1663 relname = rel_make_name(sys->system, *relptr);
1664 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1665 ASC_FREE(relname);
1666 is_error = 1;
1667 }
1668 Asc_SignalHandlerPopDefault(SIGFPE);
1669 #endif
1670
1671 if(is_error){
1672 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1673 return 1;
1674 }
1675 return 0;
1676 }
1677
1678 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1679 static int integrator_ida_sjex(long int Neq, realtype tt
1680 , N_Vector yy, N_Vector yp, N_Vector rr
1681 , realtype c_j, void *jac_data, mtx_matrix_t Jac
1682 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1683 ){
1684 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1685 return -1;
1686 }
1687
1688 /* root finding function */
1689
1690 int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data){
1691 IntegratorSystem *sys;
1692 IntegratorIdaData *enginedata;
1693 int i;
1694 #ifdef ROOT_DEBUG
1695 char *relname;
1696 #endif
1697
1698 asc_assert(g_data!=NULL);
1699 sys = (IntegratorSystem *)g_data;
1700 enginedata = integrator_ida_enginedata(sys);
1701
1702 /* pass the values of everything back to the compiler */
1703 integrator_set_t(sys, (double)tt);
1704 integrator_set_y(sys, NV_DATA_S(yy));
1705 integrator_set_ydot(sys, NV_DATA_S(yp));
1706
1707 asc_assert(gout!=NULL);
1708
1709 #ifdef ROOT_DEBUG
1710 CONSOLE_DEBUG("t = %f",tt);
1711 #endif
1712
1713 /* evaluate the residuals for each of the boundaries */
1714 for(i=0; i < enginedata->nbnds; ++i){
1715 switch(bnd_kind(enginedata->bndlist[i])){
1716 case e_bnd_rel: /* real-valued boundary relation */
1717 gout[i] = bndman_real_eval(enginedata->bndlist[i]);
1718 #ifdef ROOT_DEBUG
1719 relname = bnd_make_name(sys->system,enginedata->bndlist[i]);
1720 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
1721 ASC_FREE(relname);
1722 #endif
1723 break;
1724 case e_bnd_logrel:
1725 if(bndman_log_eval(enginedata->bndlist[i])){
1726 CONSOLE_DEBUG("bnd[%d] = TRUE",i);
1727 #ifdef ROOT_DEBUG
1728 relname = bnd_make_name(sys->system,enginedata->bndlist[i]);
1729 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
1730 ASC_FREE(relname);
1731 #endif
1732 gout[i] = +1.0;
1733 }else{
1734 CONSOLE_DEBUG("bnd[%d] = FALSE",i);
1735 gout[i] = -1.0;
1736 }
1737 break;
1738 case e_bnd_undefined:
1739 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid boundary type e_bnd_undefined");
1740 return 1;
1741 }
1742 }
1743
1744 return 0; /* no way to detect errors in bndman_*_eval at this stage */
1745 }
1746
1747
1748 /*----------------------------------------------
1749 FULL JACOBIAN PRECONDITIONER -- EXPERIMENTAL.
1750 */
1751
1752 static void integrator_ida_pcreate_jacobian(IntegratorSystem *sys){
1753 IntegratorIdaData *enginedata =sys->enginedata;
1754 IntegratorIdaPrecDataJacobian *precdata;
1755 precdata = ASC_NEW(IntegratorIdaPrecDataJacobian);
1756 mtx_matrix_t P;
1757 asc_assert(sys->n_y);
1758 precdata->L = linsolqr_create_default();
1759
1760 /* allocate matrix to be used by linsolqr */
1761 P = mtx_create();
1762 mtx_set_order(P, sys->n_y);
1763 linsolqr_set_matrix(precdata->L, P);
1764
1765 enginedata->pfree = &integrator_ida_pfree_jacobian;
1766 enginedata->precdata = precdata;
1767 CONSOLE_DEBUG("Allocated memory for Full Jacobian preconditioner");
1768 }
1769
1770 static void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata){
1771 mtx_matrix_t P;
1772 IntegratorIdaPrecDataJacobian *precdata;
1773
1774 if(enginedata->precdata){
1775 precdata = (IntegratorIdaPrecDataJacobian *)enginedata->precdata;
1776 P = linsolqr_get_matrix(precdata->L);
1777 mtx_destroy(P);
1778 linsolqr_destroy(precdata->L);
1779 ASC_FREE(precdata);
1780 enginedata->precdata = NULL;
1781
1782 CONSOLE_DEBUG("Freed memory for Full Jacobian preconditioner");
1783 }
1784 enginedata->pfree = NULL;
1785 }
1786
1787 /**
1788 EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1789
1790 'setup' function.
1791 */
1792 static int integrator_ida_psetup_jacobian(realtype tt,
1793 N_Vector yy, N_Vector yp, N_Vector rr,
1794 realtype c_j, void *p_data,
1795 N_Vector tmp1, N_Vector tmp2,
1796 N_Vector tmp3
1797 ){
1798 int i, j, res;
1799 IntegratorSystem *sys;
1800 IntegratorIdaData *enginedata;
1801 IntegratorIdaPrecDataJacobian *precdata;
1802 linsolqr_system_t L;
1803 mtx_matrix_t P;
1804 struct rel_relation **relptr;
1805
1806 sys = (IntegratorSystem *)p_data;
1807 enginedata = sys->enginedata;
1808 precdata = (IntegratorIdaPrecDataJacobian *)(enginedata->precdata);
1809 double *derivatives;
1810 struct var_variable **variables;
1811 int count, status;
1812 char *relname;
1813 mtx_coord_t C;
1814
1815 L = precdata->L;
1816 P = linsolqr_get_matrix(L);
1817 mtx_clear(P);
1818
1819 CONSOLE_DEBUG("Setting up Jacobian preconditioner");
1820
1821 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1822 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1823
1824 /**
1825 @TODO FIXME here we are using the very inefficient and contorted approach
1826 of calculating the whole jacobian, then extracting just the diagonal elements.
1827 */
1828
1829 for(i=0, relptr = enginedata->rellist;
1830 i< enginedata->nrels && relptr != NULL;
1831 ++i, ++relptr
1832 ){
1833 /* get derivatives for this particular relation */
1834 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1835 if(status){
1836 relname = rel_make_name(sys->system, *relptr);
1837 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1838 ASC_FREE(relname);
1839 break;
1840 }
1841 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1842 /* find the diagonal elements */
1843 for(j=0; j<count; ++j){
1844 if(var_deriv(variables[j])){
1845 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), c_j * derivatives[j]);
1846 }else{
1847 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), derivatives[j]);
1848 }
1849 }
1850 }
1851
1852 mtx_assemble(P);
1853
1854 if(status){
1855 CONSOLE_DEBUG("Error found when evaluating derivatives");
1856 res = 1; goto finish; /* recoverable */
1857 }
1858
1859 integrator_ida_write_incidence(sys);
1860
1861 res = 0;
1862 finish:
1863 ASC_FREE(variables);
1864 ASC_FREE(derivatives);
1865 return res;
1866 };
1867
1868 /**
1869 EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1870
1871 'solve' function.
1872 */
1873 static int integrator_ida_psolve_jacobian(realtype tt,
1874 N_Vector yy, N_Vector yp, N_Vector rr,
1875 N_Vector rvec, N_Vector zvec,
1876 realtype c_j, realtype delta, void *p_data,
1877 N_Vector tmp
1878 ){
1879 IntegratorSystem *sys;
1880 IntegratorIdaData *data;
1881 IntegratorIdaPrecDataJacobian *precdata;
1882 sys = (IntegratorSystem *)p_data;
1883 data = sys->enginedata;
1884 precdata = (IntegratorIdaPrecDataJacobian *)(data->precdata);
1885 linsolqr_system_t L = precdata->L;
1886
1887 linsolqr_add_rhs(L,NV_DATA_S(rvec),FALSE);
1888
1889 mtx_region_t R;
1890 R.row.low = R.col.low = 0;
1891 R.row.high = R.col.high = mtx_order(linsolqr_get_matrix(L)) - 1;
1892 linsolqr_set_region(L,R);
1893
1894 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
1895 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
1896
1897 /// @TODO more here
1898
1899 linsolqr_remove_rhs(L,NV_DATA_S(rvec));
1900
1901 CONSOLE_DEBUG("Solving Jacobian preconditioner (c_j = %f)",c_j);
1902 return 0;
1903 };
1904
1905
1906 /*----------------------------------------------
1907 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1908 */
1909
1910 static void integrator_ida_pcreate_jacobi(IntegratorSystem *sys){
1911 IntegratorIdaData *enginedata =sys->enginedata;
1912 IntegratorIdaPrecDataJacobi *precdata;
1913 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1914
1915 asc_assert(sys->n_y);
1916 precdata->PIii = N_VNew_Serial(sys->n_y);
1917
1918 enginedata->pfree = &integrator_ida_pfree_jacobi;
1919 enginedata->precdata = precdata;
1920 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1921 }
1922
1923 static void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1924 if(enginedata->precdata){
1925 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1926 N_VDestroy_Serial(precdata->PIii);
1927
1928 ASC_FREE(precdata);
1929 enginedata->precdata = NULL;
1930 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1931 }
1932 enginedata->pfree = NULL;
1933 }
1934
1935 /**
1936 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1937
1938 'setup' function.
1939 */
1940 static int integrator_ida_psetup_jacobi(realtype tt,
1941 N_Vector yy, N_Vector yp, N_Vector rr,
1942 realtype c_j, void *p_data,
1943 N_Vector tmp1, N_Vector tmp2,
1944 N_Vector tmp3
1945 ){
1946 int i, j, res;
1947 IntegratorSystem *sys;
1948 IntegratorIdaData *enginedata;
1949 IntegratorIdaPrecDataJacobi *precdata;
1950 struct rel_relation **relptr;
1951
1952 sys = (IntegratorSystem *)p_data;
1953 enginedata = sys->enginedata;
1954 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1955 double *derivatives;
1956 struct var_variable **variables;
1957 int count, status;
1958 char *relname;
1959
1960 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1961
1962 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1963 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1964
1965 /**
1966 @TODO FIXME here we are using the very inefficient and contorted approach
1967 of calculating the whole jacobian, then extracting just the diagonal elements.
1968 */
1969
1970 for(i=0, relptr = enginedata->rellist;
1971 i< enginedata->nrels && relptr != NULL;
1972 ++i, ++relptr
1973 ){
1974
1975 /* get derivatives for this particular relation */
1976 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1977 if(status){
1978 relname = rel_make_name(sys->system, *relptr);
1979 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1980 ASC_FREE(relname);
1981 break;
1982 }
1983 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1984 /* find the diagonal elements */
1985 for(j=0; j<count; ++j){
1986 if(var_sindex(variables[j])==i){
1987 if(var_deriv(variables[j])){
1988 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
1989 }else{
1990 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1991 }
1992
1993 }
1994 }
1995 #ifdef PREC_DEBUG
1996 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1997 #endif
1998 }
1999
2000 if(status){
2001 CONSOLE_DEBUG("Error found when evaluating derivatives");
2002 res = 1; goto finish; /* recoverable */
2003 }
2004
2005 integrator_ida_write_incidence(sys);
2006
2007 res = 0;
2008 finish:
2009 ASC_FREE(variables);
2010 ASC_FREE(derivatives);
2011 return res;
2012 };
2013
2014 /**
2015 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
2016
2017 'solve' function.
2018 */
2019 static int integrator_ida_psolve_jacobi(realtype tt,
2020 N_Vector yy, N_Vector yp, N_Vector rr,
2021 N_Vector rvec, N_Vector zvec,
2022 realtype c_j, realtype delta, void *p_data,
2023 N_Vector tmp
2024 ){
2025 IntegratorSystem *sys;
2026 IntegratorIdaData *data;
2027 IntegratorIdaPrecDataJacobi *precdata;
2028 sys = (IntegratorSystem *)p_data;
2029 data = sys->enginedata;
2030 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
2031
2032 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
2033 N_VProd(precdata->PIii, rvec, zvec);
2034 return 0;
2035 };
2036
2037 /*----------------------------------------------
2038 STATS
2039 */
2040
2041 /**
2042 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
2043 status in a struct instead of separately.
2044
2045 @return IDA_SUCCESS on success.
2046 */
2047 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
2048
2049 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
2050
2051 int res;
2052
2053 /*
2054 There is an error in the documentation for this function in Sundials 2.2.
2055 According the the header file, the hinused stat is not provided.
2056 */
2057 res = IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups,
2058 &s->netfails, &s->qlast, &s->qcur, &s->hlast, &s->hcur,
2059 &s->tcur
2060 );
2061
2062 /* get the missing statistic */
2063 IDAGetActualInitStep(ida_mem, &s->hinused);
2064
2065 return res;
2066 #else
2067
2068 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
2069 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
2070 ,&s->hlast, &s->hcur, &s->tcur
2071 );
2072
2073 #endif
2074 }
2075
2076 /**
2077 This routine just outputs the stats to the CONSOLE_DEBUG routine.
2078
2079 @TODO provide a GUI way of stats reporting from IDA.
2080 */
2081 static void integrator_ida_write_stats(IntegratorIdaStats *stats){
2082 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
2083 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
2084 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
2085 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
2086 SI(qlast); SI(qcur);
2087 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
2088 # undef SL
2089 # undef SI
2090 # undef SR
2091 }
2092
2093 /*------------------------------------------------------------------------------
2094 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
2095 */
2096
2097 /**
2098 Here we construct the local transfer matrix. It's a bit of a naive
2099 approach; probably far more efficient ways can be worked out. But it will
2100 hopefully be a useful way to examine stability of some spatial
2101 discretisation schemes for PDAE systems.
2102
2103 http://ascendserver.cheme.cmu.edu/wiki/index.php/IDA#Stability
2104 */
2105 static int integrator_ida_transfer_matrix(const IntegratorSystem *sys, struct SystemJacobianStruct *J){
2106 int i=0, res;
2107 enum submat{II_GA=0, II_GD, II_FA, II_FD, II_FDP, II_NUM};
2108
2109 const var_filter_t *matvf[II_NUM] = {
2110 &system_vfilter_algeb
2111 ,&system_vfilter_diff
2112 ,&system_vfilter_algeb
2113 ,&system_vfilter_diff
2114 ,&system_vfilter_deriv
2115 };
2116
2117 const rel_filter_t *matrf[II_NUM] = {
2118 &system_rfilter_algeb
2119 ,&system_rfilter_algeb
2120 ,&system_rfilter_diff
2121 ,&system_rfilter_diff
2122 ,&system_rfilter_diff
2123 };
2124
2125 struct SystemJacobianStruct D[II_NUM];
2126
2127 for(i=0;i<II_NUM;++i){
2128 res = system_jacobian(sys->system, matrf[i], matvf[i], 1/*safe*/ ,&(D[i]));
2129 }
2130
2131 /* compute inverses for matrices that need it */
2132 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
2133 return 1;
2134 }
2135
2136
2137 /**
2138 Our task here is to write the matrices that IDA *should* be seeing. We
2139 are actually making calls to relman_diff in order to do that, so we're
2140 really going back to the variables in the actualy system and computing
2141 row by row what the values are. This should mean just a single call to
2142 each blackbox present in the system (if blackbox caching is working
2143 correctly).
2144 */
2145 static int integrator_ida_write_matrix(const IntegratorSystem *sys, FILE *f, const char *type){
2146 /* IntegratorIdaData *enginedata; */
2147 struct SystemJacobianStruct J = {NULL,NULL,NULL,0,0};
2148 int status=1;
2149 mtx_region_t R;
2150
2151 if(type==NULL)type = "dx'/dx";
2152
2153 if(0==strcmp(type,"dg/dz")){
2154 CONSOLE_DEBUG("Calculating dg/dz...");
2155 status = system_jacobian(sys->system
2156 , &system_rfilter_algeb, &system_vfilter_algeb
2157 , 1 /* safe */
2158 , &J
2159 );
2160 }else if(0==strcmp(type,"dg/dx")){
2161 CONSOLE_DEBUG("Calculating dg/dx...");
2162 status = system_jacobian(sys->system
2163 , &system_rfilter_algeb, &system_vfilter_diff
2164 , 1 /* safe */
2165 , &J
2166 );
2167 }else if(0==strcmp(type,"df/dx'")){
2168 CONSOLE_DEBUG("Calculating df/dx'...");
2169 status = system_jacobian(sys->system
2170 , &system_rfilter_diff, &system_vfilter_deriv
2171 , 1 /* safe */
2172 , &J
2173 );
2174 }else if(0==strcmp(type,"df/dz")){
2175 CONSOLE_DEBUG("Calculating df/dz...");
2176 status = system_jacobian(sys->system
2177 , &system_rfilter_diff, &system_vfilter_algeb
2178 , 1 /* safe */
2179 , &J
2180 );
2181 }else if(0==strcmp(type,"df/dx")){
2182 CONSOLE_DEBUG("Calculating df/dx...");
2183 status = system_jacobian(sys->system
2184 , &system_rfilter_diff, &system_vfilter_diff
2185 , 1 /* safe */
2186 , &J
2187 );
2188 }else if(0==strcmp(type,"dF/dy")){
2189 CONSOLE_DEBUG("Calculating dF/dy...");
2190 status = system_jacobian(sys->system
2191 , &system_rfilter_all, &system_vfilter_nonderiv
2192 , 1 /* safe */
2193 , &J
2194 );
2195 }else if(0==strcmp(type,"dF/dy'")){
2196 CONSOLE_DEBUG("Calculating dF/dy'...");
2197 status = system_jacobian(sys->system
2198 , &system_rfilter_all, &system_vfilter_deriv
2199 , 1 /* safe */
2200 , &J
2201 );
2202 }else if(0==strcmp(type,"dx'/dx")){
2203 /* system state transfer matrix dyd'/dyd */
2204 status = integrator_ida_transfer_matrix(sys, &J);
2205 }else{
2206 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid matrix type '%s'",type);
2207 return 1;
2208 }
2209
2210 if(status){
2211 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating matrix");
2212 }else{
2213 /* send the region explicitly, so that we handle non-square correctly */
2214 R.row.low = 0; R.col.low = 0;
2215 R.row.high = J.n_rels - 1; R.col.high = J.n_vars - 1;
2216 /* note that we're not fussy about empty matrices here... */
2217 mtx_write_region_mmio(f,J.M,&R);
2218 }
2219
2220 if(J.vars)ASC_FREE(J.vars);
2221 if(J.rels)ASC_FREE(J.rels);
2222 if(J.M)mtx_destroy(J.M);
2223
2224 return status;
2225 }
2226
2227 /**
2228 This routine outputs matrix structure in a crude text format, for the sake
2229 of debugging.
2230 */
2231 static void integrator_ida_write_incidence(IntegratorSystem *sys){
2232 int i, j;
2233 struct rel_relation **relptr;
2234 IntegratorIdaData *enginedata = sys->enginedata;
2235 double *derivatives;
2236 struct var_variable **variables;
2237 int count, status;
2238 char *relname;
2239
2240 if(enginedata->nrels > 100){
2241 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
2242 return;
2243 }
2244
2245 variables = ASC_NEW_ARRAY(struct var_variable *, sys->n_y * 2);
2246 derivatives = ASC_NEW_ARRAY(double, sys->n_y * 2);
2247
2248 CONSOLE_DEBUG("Outputting incidence information to console...");
2249
2250 for(i=0, relptr = enginedata->rellist;
2251 i< enginedata->nrels && relptr != NULL;
2252 ++i, ++relptr
2253 ){
2254 relname = rel_make_name(sys->system, *relptr);
2255
2256 /* get derivatives for this particular relation */
2257 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
2258 if(status){
2259 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
2260 ASC_FREE(relname);
2261 break;
2262 }
2263
2264 fprintf(stderr,"%3d:%-15s:",i,relname);
2265 ASC_FREE(relname);
2266
2267 for(j=0; j<count; ++j){
2268 if(var_deriv(variables[j])){
2269 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(sys,variables[j]));
2270 }else{
2271 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
2272 }
2273 }
2274 fprintf(stderr,"\n");
2275 }
2276 ASC_FREE(variables);
2277 ASC_FREE(derivatives);
2278 }
2279
2280 /* @return 0 on success */
2281 int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
2282 char *varname, *relname;
2283 struct var_variable **vlist, *var;
2284 struct rel_relation **rlist, *rel;
2285 long vlen, rlen;
2286 long i;
2287 long di;
2288
2289 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",sys->n_y);
2290
2291 /* if(integrator_sort_obs_vars(sys))return 10; */
2292
2293 if(sys->y && sys->ydot){
2294 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
2295 fprintf(fp,"index\t%-15s\tydot\n","y");
2296 fprintf(fp,"-----\t%-15s\t-----\n","-----");
2297 for(i=0;i<sys->n_y;++i){
2298 varname = var_make_name(sys->system, sys->y[i]);
2299 fprintf(fp,"%ld\t%-15s\t",i,varname);
2300 if(sys->ydot[i]){
2301 ASC_FREE(varname);
2302 varname = var_make_name(sys->system, sys->ydot[i]);
2303 fprintf(fp,"%s\n",varname);
2304 ASC_FREE(varname);
2305 }else{
2306 fprintf(fp,".\n");
2307 ASC_FREE(varname);
2308 }
2309 }
2310 }else{
2311 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
2312 }
2313
2314 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
2315 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
2316 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
2317
2318
2319 /* visit all the slv_system_t master var lists to collect vars */
2320 /* find the vars mostly in this one */
2321 vlist = slv_get_solvers_var_list(sys->system);
2322 vlen = slv_get_num_solvers_vars(sys->system);
2323 for(i=0;i<vlen;i++){
2324 var = vlist[i];
2325
2326 varname = var_make_name(sys->system, var);
2327 fprintf(fp,"%ld\t%-15s\t",i,varname);
2328
2329 if(var_fixed(var)){
2330 // it's fixed, so not really a DAE var
2331 fprintf(fp,"(fixed)\n");
2332 }else if(!var_active(var)){
2333 // inactive
2334 fprintf(fp,"(inactive)\n");
2335 }else if(!var_incident(var)){
2336 // not incident
2337 fprintf(fp,"(not incident)\n");
2338 }else{
2339 if(var_deriv(var)){
2340 if(sys->y_id){
2341 di = integrator_ida_diffindex1(sys,var);
2342 if(di>=0){
2343 ASC_FREE(varname);
2344 varname = var_make_name(sys->system,vlist[di]);
2345 fprintf(fp,".\tdiff(%ld='%s')\n",di,varname);
2346 }else{
2347 fprintf(fp,".\tdiff(???,err=%ld)\n",di);
2348 }
2349 }else{
2350 fprintf(fp,".\tderiv... of??\n");
2351 }
2352 }else{
2353 fprintf(fp,"%d\t.\n",var_sindex(var));
2354 }
2355 }
2356 ASC_FREE(varname);
2357 }
2358
2359 /* let's write out the relations too */
2360 rlist = slv_get_solvers_rel_list(sys->system);
2361 rlen = slv_get_num_solvers_rels(sys->system);
2362
2363 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
2364 fprintf(fp,"index\tname\n");
2365 fprintf(fp,"-----\t----\n");
2366 for(i=0; i<rlen; ++i){
2367 rel = rlist[i];
2368 relname = rel_make_name(sys->system,rel);
2369 fprintf(fp,"%ld\t%s\n",i,relname);
2370 ASC_FREE(relname);
2371 }
2372
2373 /* write out the derivative chains */
2374 fprintf(fp,"\nDERIVATIVE CHAINS\n");
2375 if(integrator_ida_analyse_debug(sys,stderr)){
2376 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
2377 return 340;
2378 }
2379 fprintf(fp,"\n");
2380
2381 /* and lets write block debug output */
2382 system_block_debug(sys->system, fp);
2383
2384 return 0; /* success */
2385 }
2386
2387 /*----------------------------------------------
2388 ERROR REPORTING
2389 */
2390 /**
2391 Error message reporter function to be passed to IDA. All error messages
2392 will trigger a call to this function, so we should find everything
2393 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
2394 panel (in the case of PyGTK).
2395 */
2396 static void integrator_ida_error(int error_code
2397 , const char *module, const char *function
2398 , char *msg, void *eh_data
2399 ){
2400 IntegratorSystem *sys;
2401 error_severity_t sev;
2402
2403 /* cast back the IntegratorSystem, just in case we need it */
2404 sys = (IntegratorSystem *)eh_data;
2405
2406 /* severity depends on the sign of the error_code value */
2407 if(error_code <= 0){
2408 sev = ASC_PROG_ERR;
2409 }else{
2410 sev = ASC_PROG_WARNING;
2411 }
2412
2413 /* use our all-purpose error reporting to get stuff back to the GUI */
2414 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
2415 }

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