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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2111 - (show annotations) (download) (as text)
Tue Dec 8 03:05:36 2009 UTC (10 years, 7 months ago) by jpye
File MIME type: text/x-csrc
File size: 72165 byte(s)
CO2 correlation looks to be working now.
Work on support for SUNDIALS on Windows.
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
1065 // FIXME this stuff has not been tested yet, and is very incomplete.
1066
1067 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
1068 /* vector of absolute tolerances */
1069 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
1070 abstolvect = N_VNew_Serial(sys->n_y);
1071 integrator_get_atol(sys,NV_DATA_S(abstolvect));
1072 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SV, reltol, abstolvect);
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 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
1078 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SS, reltol, &abstol);
1079 }
1080 #else
1081 /* allocate internal memory */
1082 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
1083 /* vector of absolute tolerances */
1084 abstolvect = N_VNew_Serial(sys->n_y);
1085 integrator_get_atol(sys,NV_DATA_S(abstolvect));
1086 if(IDA_SUCCESS != IDAReInit(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect)){
1087 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
1088 }
1089 N_VDestroy_Serial(abstolvect);
1090 }else{
1091 /* scalar absolute tolerance (one value for all) */
1092 abstol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_ATOL);
1093 if(IDA_SUCCESS != IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol)){
1094 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
1095 }
1096 }
1097 #endif
1098 }
1099 }else{
1100 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");
1101 }
1102 ASC_FREE(rootsfound);
1103 }
1104
1105
1106 /* pass the values of everything back to the compiler */
1107 integrator_set_t(sys, (double)tret);
1108 integrator_set_y(sys, NV_DATA_S(yret));
1109 integrator_set_ydot(sys, NV_DATA_S(ypret));
1110
1111 if(flag<0){
1112 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
1113 break;
1114 }
1115
1116 /* -- do something so that sys knows the values of tret, yret and ypret */
1117
1118 /* -- store the current values of all the stuff */
1119 integrator_output_write(sys);
1120 integrator_output_write_obs(sys);
1121
1122 }/* loop through next sample timestep */
1123
1124 /* -- close the IntegratorReporter */
1125 integrator_output_close(sys);
1126
1127 /* get optional outputs */
1128 #ifdef STATS_DEBUG
1129 IntegratorIdaStats stats;
1130 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
1131 integrator_ida_write_stats(&stats);
1132 }else{
1133 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
1134 }
1135 #endif
1136
1137 /* free solution memory */
1138 N_VDestroy_Serial(yret);
1139 N_VDestroy_Serial(ypret);
1140
1141 /* free solver memory */
1142 IDAFree(ida_mem);
1143
1144 if(flag < -500){
1145 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", t);
1146 return -flag;
1147 }
1148
1149 if(flag < 0){
1150 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1151 return 14;
1152 }
1153
1154 /* all done, success */
1155 return 0;
1156 }
1157
1158 /*--------------------------------------------------
1159 RESIDUALS AND JACOBIAN AND IDAROOTFN
1160 */
1161
1162 #if 0
1163 typedef void (SignalHandlerFn)(int);
1164 SignalHandlerFn integrator_ida_sig;
1165 SignalHandlerFn *integrator_ida_sig_old;
1166 jmp_buf integrator_ida_jmp_buf;
1167 fenv_t integrator_ida_fenv_old;
1168
1169
1170 void integrator_ida_write_feinfo(){
1171 int f;
1172 f = fegetexcept();
1173 CONSOLE_DEBUG("Locating nature of exception...");
1174 if(f & FE_DIVBYZERO)ERROR_REPORTER_HERE(ASC_PROG_ERR,"DIV BY ZERO");
1175 if(f & FE_INEXACT)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INEXACT");
1176 if(f & FE_INVALID)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INVALID");
1177 if(f & FE_OVERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"OVERFLOW");
1178 if(f & FE_UNDERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"UNDERFLOW");
1179 if(f==0)ERROR_REPORTER_HERE(ASC_PROG_ERR,"FLAGS ARE CLEAR?!?");
1180 }
1181
1182 void integrator_ida_sig(int sig){
1183 /* the wrong signal: rethrow to the default handler */
1184 if(sig!=SIGFPE){
1185 signal(SIGFPE,SIG_DFL);
1186 raise(sig);
1187 }
1188 integrator_ida_write_feinfo();
1189 CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);
1190 longjmp(integrator_ida_jmp_buf,sig);
1191 }
1192 #endif
1193
1194 /**
1195 Function to evaluate system residuals, in the form required for IDA.
1196
1197 Given tt, yy and yp, we need to evaluate and return rr.
1198
1199 @param tt current value of indep variable (time)
1200 @param yy current values of dependent variable vector
1201 @param yp current values of derivatives of dependent variables
1202 @param rr the output residual vector (is we're returning data to)
1203 @param res_data pointer to our stuff (sys in this case).
1204
1205 @return 0 on success, positive on recoverable error, and
1206 negative on unrecoverable error.
1207 */
1208 static int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1209 IntegratorSystem *sys;
1210 IntegratorIdaData *enginedata;
1211 int i, calc_ok, is_error;
1212 struct rel_relation** relptr;
1213 double resid;
1214 char *relname;
1215 #ifdef FEX_DEBUG
1216 char *varname;
1217 char diffname[30];
1218 #endif
1219
1220 sys = (IntegratorSystem *)res_data;
1221 enginedata = integrator_ida_enginedata(sys);
1222
1223 #ifdef FEX_DEBUG
1224 /* fprintf(stderr,"\n\n"); */
1225 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1226 #endif
1227
1228 if(NV_LENGTH_S(rr)!=enginedata->nrels){
1229 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
1230 return -1; /* unrecoverable */
1231 }
1232
1233 /* pass the values of everything back to the compiler */
1234 integrator_set_t(sys, (double)tt);
1235 integrator_set_y(sys, NV_DATA_S(yy));
1236 integrator_set_ydot(sys, NV_DATA_S(yp));
1237
1238 /* perform bounds checking on all variables */
1239 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1240 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1241 return 1;
1242 }
1243
1244 /* evaluate each residual in the rellist */
1245 is_error = 0;
1246 relptr = enginedata->rellist;
1247
1248
1249 #ifdef ASC_SIGNAL_TRAPS
1250 if(enginedata->safeeval){
1251 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1252 }else{
1253 # ifdef FEX_DEBUG
1254 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1255 # endif
1256 Asc_SignalHandlerPushDefault(SIGFPE);
1257 }
1258
1259 if (SETJMP(g_fpe_env)==0) {
1260 #endif
1261
1262
1263 for(i=0, relptr = enginedata->rellist;
1264 i< enginedata->nrels && relptr != NULL;
1265 ++i, ++relptr
1266 ){
1267 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
1268
1269 NV_Ith_S(rr,i) = resid;
1270 if(!calc_ok){
1271 relname = rel_make_name(sys->system, *relptr);
1272 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1273 ASC_FREE(relname);
1274 /* presumable some output already made? */
1275 is_error = 1;
1276 }/*else{
1277 CONSOLE_DEBUG("Calc OK");
1278 }*/
1279 }
1280
1281 if(!is_error){
1282 for(i=0;i< enginedata->nrels; ++i){
1283 if(isnan(NV_Ith_S(rr,i))){
1284 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in residual %d",i);
1285 is_error=1;
1286 }
1287 }
1288 #ifdef FEX_DEBUG
1289 if(!is_error){
1290 CONSOLE_DEBUG("No NAN detected");
1291 }
1292 #endif
1293 }
1294
1295 #ifdef ASC_SIGNAL_TRAPS
1296 }else{
1297 relname = rel_make_name(sys->system, *relptr);
1298 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1299 ASC_FREE(relname);
1300 is_error = 1;
1301 }
1302
1303 if(enginedata->safeeval){
1304 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1305 }else{
1306 Asc_SignalHandlerPopDefault(SIGFPE);
1307 }
1308 #endif
1309
1310
1311 #ifdef FEX_DEBUG
1312 /* output residuals to console */
1313 CONSOLE_DEBUG("RESIDUAL OUTPUT");
1314 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1315 for(i=0; i<sys->n_y; ++i){
1316 varname = var_make_name(sys->system,sys->y[i]);
1317 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1318 if(sys->ydot[i]){
1319 varname = var_make_name(sys->system,sys->ydot[i]);
1320 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1321 }else{
1322 snprintf(diffname,99,"diff(%s)",varname);
1323 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1324 }
1325 ASC_FREE(varname);
1326 relname = rel_make_name(sys->system,enginedata->rellist[i]);
1327 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1328 }
1329 #endif
1330
1331 if(is_error){
1332 return 1;
1333 }
1334
1335 #ifdef FEX_DEBUG
1336 CONSOLE_DEBUG("RESIDUAL OK");
1337 #endif
1338 return 0;
1339 }
1340
1341 /**
1342 Dense Jacobian evaluation. Only suitable for small problems!
1343 Has been seen working for problems up to around 2000 vars, FWIW.
1344 */
1345 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
1346 static int integrator_ida_djex(int Neq, realtype tt, realtype c_j
1347 , N_Vector yy, N_Vector yp, N_Vector rr
1348 , IDA_MTX_T Jac, void *jac_data
1349 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1350 ){
1351 #else
1352 static int integrator_ida_djex(long int Neq, realtype tt
1353 , N_Vector yy, N_Vector yp, N_Vector rr
1354 , realtype c_j, void *jac_data, IDA_MTX_T Jac
1355 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1356 ){
1357 #endif
1358 IntegratorSystem *sys;
1359 IntegratorIdaData *enginedata;
1360 char *relname;
1361 #ifdef DJEX_DEBUG
1362 struct var_variable **varlist;
1363 char *varname;
1364 #endif
1365 struct rel_relation **relptr;
1366 int i;
1367 double *derivatives;
1368 struct var_variable **variables;
1369 int count, j;
1370 int status, is_error = 0;
1371
1372 sys = (IntegratorSystem *)jac_data;
1373 enginedata = integrator_ida_enginedata(sys);
1374
1375 /* allocate space for returns from relman_diff3 */
1376 /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
1377 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1378 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1379
1380 /* pass the values of everything back to the compiler */
1381 integrator_set_t(sys, (double)tt);
1382 integrator_set_y(sys, NV_DATA_S(yy));
1383 integrator_set_ydot(sys, NV_DATA_S(yp));
1384
1385 /* perform bounds checking on all variables */
1386 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1387 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1388 return 1;
1389 }
1390
1391 #ifdef DJEX_DEBUG
1392 varlist = slv_get_solvers_var_list(sys->system);
1393
1394 /* print vars */
1395 for(i=0; i < sys->n_y; ++i){
1396 varname = var_make_name(sys->system, sys->y[i]);
1397 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1398 asc_assert(NV_Ith_S(yy,i) == var_value(sys->y[i]));
1399 ASC_FREE(varname);
1400 }
1401
1402 /* print derivatives */
1403 for(i=0; i < sys->n_y; ++i){
1404 if(sys->ydot[i]){
1405 varname = var_make_name(sys->system, sys->ydot[i]);
1406 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(sys->ydot[i]));
1407 ASC_FREE(varname);
1408 }else{
1409 varname = var_make_name(sys->system, sys->y[i]);
1410 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1411 ASC_FREE(varname);
1412 }
1413 }
1414
1415 /* print step size */
1416 CONSOLE_DEBUG("<c_j> = %g",c_j);
1417 #endif
1418
1419 /* build up the dense jacobian matrix... */
1420 status = 0;
1421 for(i=0, relptr = enginedata->rellist;
1422 i< enginedata->nrels && relptr != NULL;
1423 ++i, ++relptr
1424 ){
1425 /* get derivatives for this particular relation */
1426 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1427
1428 if(status){
1429 relname = rel_make_name(sys->system, *relptr);
1430 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1431 ASC_FREE(relname);
1432 is_error = 1;
1433 break;
1434 }
1435
1436 /* output what's going on here ... */
1437 #ifdef DJEX_DEBUG
1438 relname = rel_make_name(sys->system, *relptr);
1439 fprintf(stderr,"%d: '%s': ",i,relname);
1440 for(j=0;j<count;++j){
1441 varname = var_make_name(sys->system, variables[j]);
1442 if(var_deriv(variables[j])){
1443 fprintf(stderr," '%s'=",varname);
1444 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(sys,variables[j]));
1445 }else{
1446 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1447 }
1448 ASC_FREE(varname);
1449 }
1450 /* relname is freed further down */
1451 fprintf(stderr,"\n");
1452 #endif
1453
1454 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1455 for(j=0; j < count; ++j){
1456 #ifdef DJEX_DEBUG
1457 varname = var_make_name(sys->system,variables[j]);
1458 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1459 ASC_FREE(varname);
1460 #endif
1461 if(!var_deriv(variables[j])){
1462 #ifdef DJEX_DEBUG
1463 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
1464 asc_assert(var_sindex(variables[j]) >= 0);
1465 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1466 #endif
1467 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1468 }else{
1469 DENSE_ELEM(Jac,i,integrator_ida_diffindex(sys,variables[j])) += derivatives[j] * c_j;
1470 #ifdef DJEX_DEBUG
1471 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1472 #endif
1473 }
1474 }
1475 }
1476
1477 #ifdef DJEX_DEBUG
1478 ASC_FREE(relname);
1479 CONSOLE_DEBUG("PRINTING JAC");
1480 fprintf(stderr,"\t");
1481 for(j=0; j < sys->n_y; ++j){
1482 if(j)fprintf(stderr,"\t");
1483 varname = var_make_name(sys->system,sys->y[j]);
1484 fprintf(stderr,"%11s",varname);
1485 ASC_FREE(varname);
1486 }
1487 fprintf(stderr,"\n");
1488 for(i=0; i < enginedata->nrels; ++i){
1489 relname = rel_make_name(sys->system, enginedata->rellist[i]);
1490 fprintf(stderr,"%s\t",relname);
1491 ASC_FREE(relname);
1492
1493 for(j=0; j < sys->n_y; ++j){
1494 if(j)fprintf(stderr,"\t");
1495 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1496 }
1497 fprintf(stderr,"\n");
1498 }
1499 #endif
1500
1501 /* test for NANs */
1502 if(!is_error){
1503 for(i=0;i< enginedata->nrels; ++i){
1504 for(j=0;j<sys->n_y;++j){
1505 if(isnan(DENSE_ELEM(Jac,i,j))){
1506 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
1507 is_error=1;
1508 }
1509 }
1510 }
1511 #ifdef DJEX_DEBUG
1512 if(!is_error){
1513 CONSOLE_DEBUG("No NAN detected");
1514 }
1515 #endif
1516 }
1517
1518 /* if(integrator_ida_check_diffindex(sys)){
1519 is_error = 1;
1520 }*/
1521
1522 if(is_error){
1523 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1524 return 1;
1525 }
1526
1527 #ifdef DJEX_DEBUG
1528 CONSOLE_DEBUG("DJEX RETURNING 0");
1529 /* ASC_PANIC("Quitting"); */
1530 #endif
1531 return 0;
1532 }
1533
1534 /**
1535 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1536
1537 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1538
1539 @param tt current value of the independent variable (time, t)
1540 @param yy current value of the dependent variable vector, y(t).
1541 @param yp current value of y'(t).
1542 @param rr current value of the residual vector F(t, y, y').
1543 @param v the vector by which the Jacobian must be multiplied to the right.
1544 @param Jv the output vector computed
1545 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1546 @param jac_data pointer to our stuff (sys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1547 @param tmp1 @see tmp2
1548 @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.
1549 @return 0 on success
1550 */
1551 static int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1552 , N_Vector v, N_Vector Jv, realtype c_j
1553 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1554 ){
1555 IntegratorSystem *sys;
1556 IntegratorIdaData *enginedata;
1557 int i, j, is_error=0;
1558 struct rel_relation** relptr = 0;
1559 char *relname;
1560 int status;
1561 double Jv_i;
1562
1563 struct var_variable **variables;
1564 double *derivatives;
1565 int count;
1566 struct var_variable **varlist;
1567 #ifdef JEX_DEBUG
1568
1569 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1570 #endif
1571
1572 sys = (IntegratorSystem *)jac_data;
1573 enginedata = integrator_ida_enginedata(sys);
1574 varlist = slv_get_solvers_var_list(sys->system);
1575
1576 /* pass the values of everything back to the compiler */
1577 integrator_set_t(sys, (double)tt);
1578 integrator_set_y(sys, NV_DATA_S(yy));
1579 integrator_set_ydot(sys, NV_DATA_S(yp));
1580 /* no real use for residuals (rr) here, I don't think? */
1581
1582 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1583
1584 i = NV_LENGTH_S(yy) * 2;
1585 #ifdef JEX_DEBUG
1586 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1587 #endif
1588 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1589 derivatives = ASC_NEW_ARRAY(double, i);
1590
1591 /* evaluate the derivatives... */
1592 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1593
1594 #ifdef ASC_SIGNAL_TRAPS
1595 Asc_SignalHandlerPushDefault(SIGFPE);
1596 if (SETJMP(g_fpe_env)==0) {
1597 #endif
1598 for(i=0, relptr = enginedata->rellist;
1599 i< enginedata->nrels && relptr != NULL;
1600 ++i, ++relptr
1601 ){
1602 /* get derivatives for this particular relation */
1603 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1604 #ifdef JEX_DEBUG
1605 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1606 #endif
1607
1608 if(status){
1609 relname = rel_make_name(sys->system, *relptr);
1610 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1611 ASC_FREE(relname);
1612 is_error = 1;
1613 break;
1614 }
1615
1616 /*
1617 Now we have the derivatives wrt each alg/diff variable in the
1618 present equation. variables[] points into the varlist. need
1619 a mapping from the varlist to the y and ydot lists.
1620 */
1621
1622 Jv_i = 0;
1623 for(j=0; j < count; ++j){
1624 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], sys->n_y);
1625 varname = var_make_name(sys->system, enginedata->varlist[variables[j]]);
1626 if(varname){
1627 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1628 ASC_FREE(varname);
1629 }else{
1630 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1631 }
1632 */
1633
1634 /* we don't calculate derivatives wrt indep var */
1635 asc_assert(variables[j]>=0);
1636 if(variables[j] == sys->x) continue;
1637 #ifdef JEX_DEBUG
1638 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1639 #endif
1640 if(var_deriv(variables[j])){
1641 #define DIFFINDEX integrator_ida_diffindex(sys,variables[j])
1642 #ifdef JEX_DEBUG
1643 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
1644 , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
1645 , i, DIFFINDEX, derivatives[j]
1646 , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
1647 );
1648 #endif
1649 asc_assert(sys->ydot[DIFFINDEX]==variables[j]);
1650 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1651 #undef DIFFINDEX
1652 }else{
1653 #define VARINDEX var_sindex(variables[j])
1654 #ifdef JEX_DEBUG
1655 asc_assert(sys->y[VARINDEX]==variables[j]);
1656 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
1657 , i
1658 , derivatives[j] * NV_Ith_S(v,VARINDEX)
1659 , i, VARINDEX, derivatives[j]
1660 , VARINDEX, NV_Ith_S(v,VARINDEX)
1661 );
1662 #endif
1663 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1664 #undef VARINDEX
1665 }
1666 }
1667
1668 NV_Ith_S(Jv,i) = Jv_i;
1669 #ifdef JEX_DEBUG
1670 CONSOLE_DEBUG("rel = %p",*relptr);
1671 relname = rel_make_name(sys->system, *relptr);
1672 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1673 ASC_FREE(relname);
1674 return 1;
1675 #endif
1676 }
1677 #ifdef ASC_SIGNAL_TRAPS
1678 }else{
1679 relname = rel_make_name(sys->system, *relptr);
1680 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1681 ASC_FREE(relname);
1682 is_error = 1;
1683 }
1684 Asc_SignalHandlerPopDefault(SIGFPE);
1685 #endif
1686
1687 if(is_error){
1688 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1689 return 1;
1690 }
1691 return 0;
1692 }
1693
1694 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1695 static int integrator_ida_sjex(long int Neq, realtype tt
1696 , N_Vector yy, N_Vector yp, N_Vector rr
1697 , realtype c_j, void *jac_data, mtx_matrix_t Jac
1698 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1699 ){
1700 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1701 return -1;
1702 }
1703
1704 /* root finding function */
1705
1706 int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data){
1707 IntegratorSystem *sys;
1708 IntegratorIdaData *enginedata;
1709 int i;
1710 #ifdef ROOT_DEBUG
1711 char *relname;
1712 #endif
1713
1714 asc_assert(g_data!=NULL);
1715 sys = (IntegratorSystem *)g_data;
1716 enginedata = integrator_ida_enginedata(sys);
1717
1718 /* pass the values of everything back to the compiler */
1719 integrator_set_t(sys, (double)tt);
1720 integrator_set_y(sys, NV_DATA_S(yy));
1721 integrator_set_ydot(sys, NV_DATA_S(yp));
1722
1723 asc_assert(gout!=NULL);
1724
1725 #ifdef ROOT_DEBUG
1726 CONSOLE_DEBUG("t = %f",tt);
1727 #endif
1728
1729 /* evaluate the residuals for each of the boundaries */
1730 for(i=0; i < enginedata->nbnds; ++i){
1731 switch(bnd_kind(enginedata->bndlist[i])){
1732 case e_bnd_rel: /* real-valued boundary relation */
1733 gout[i] = bndman_real_eval(enginedata->bndlist[i]);
1734 #ifdef ROOT_DEBUG
1735 relname = bnd_make_name(sys->system,enginedata->bndlist[i]);
1736 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
1737 ASC_FREE(relname);
1738 #endif
1739 break;
1740 case e_bnd_logrel:
1741 if(bndman_log_eval(enginedata->bndlist[i])){
1742 CONSOLE_DEBUG("bnd[%d] = TRUE",i);
1743 #ifdef ROOT_DEBUG
1744 relname = bnd_make_name(sys->system,enginedata->bndlist[i]);
1745 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
1746 ASC_FREE(relname);
1747 #endif
1748 gout[i] = +1.0;
1749 }else{
1750 CONSOLE_DEBUG("bnd[%d] = FALSE",i);
1751 gout[i] = -1.0;
1752 }
1753 break;
1754 case e_bnd_undefined:
1755 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid boundary type e_bnd_undefined");
1756 return 1;
1757 }
1758 }
1759
1760 return 0; /* no way to detect errors in bndman_*_eval at this stage */
1761 }
1762
1763
1764 /*----------------------------------------------
1765 FULL JACOBIAN PRECONDITIONER -- EXPERIMENTAL.
1766 */
1767
1768 static void integrator_ida_pcreate_jacobian(IntegratorSystem *sys){
1769 IntegratorIdaData *enginedata =sys->enginedata;
1770 IntegratorIdaPrecDataJacobian *precdata;
1771 precdata = ASC_NEW(IntegratorIdaPrecDataJacobian);
1772 mtx_matrix_t P;
1773 asc_assert(sys->n_y);
1774 precdata->L = linsolqr_create_default();
1775
1776 /* allocate matrix to be used by linsolqr */
1777 P = mtx_create();
1778 mtx_set_order(P, sys->n_y);
1779 linsolqr_set_matrix(precdata->L, P);
1780
1781 enginedata->pfree = &integrator_ida_pfree_jacobian;
1782 enginedata->precdata = precdata;
1783 CONSOLE_DEBUG("Allocated memory for Full Jacobian preconditioner");
1784 }
1785
1786 static void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata){
1787 mtx_matrix_t P;
1788 IntegratorIdaPrecDataJacobian *precdata;
1789
1790 if(enginedata->precdata){
1791 precdata = (IntegratorIdaPrecDataJacobian *)enginedata->precdata;
1792 P = linsolqr_get_matrix(precdata->L);
1793 mtx_destroy(P);
1794 linsolqr_destroy(precdata->L);
1795 ASC_FREE(precdata);
1796 enginedata->precdata = NULL;
1797
1798 CONSOLE_DEBUG("Freed memory for Full Jacobian preconditioner");
1799 }
1800 enginedata->pfree = NULL;
1801 }
1802
1803 /**
1804 EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1805
1806 'setup' function.
1807 */
1808 static int integrator_ida_psetup_jacobian(realtype tt,
1809 N_Vector yy, N_Vector yp, N_Vector rr,
1810 realtype c_j, void *p_data,
1811 N_Vector tmp1, N_Vector tmp2,
1812 N_Vector tmp3
1813 ){
1814 int i, j, res;
1815 IntegratorSystem *sys;
1816 IntegratorIdaData *enginedata;
1817 IntegratorIdaPrecDataJacobian *precdata;
1818 linsolqr_system_t L;
1819 mtx_matrix_t P;
1820 struct rel_relation **relptr;
1821
1822 sys = (IntegratorSystem *)p_data;
1823 enginedata = sys->enginedata;
1824 precdata = (IntegratorIdaPrecDataJacobian *)(enginedata->precdata);
1825 double *derivatives;
1826 struct var_variable **variables;
1827 int count, status;
1828 char *relname;
1829 mtx_coord_t C;
1830
1831 L = precdata->L;
1832 P = linsolqr_get_matrix(L);
1833 mtx_clear(P);
1834
1835 CONSOLE_DEBUG("Setting up Jacobian preconditioner");
1836
1837 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1838 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1839
1840 /**
1841 @TODO FIXME here we are using the very inefficient and contorted approach
1842 of calculating the whole jacobian, then extracting just the diagonal elements.
1843 */
1844
1845 for(i=0, relptr = enginedata->rellist;
1846 i< enginedata->nrels && relptr != NULL;
1847 ++i, ++relptr
1848 ){
1849 /* get derivatives for this particular relation */
1850 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1851 if(status){
1852 relname = rel_make_name(sys->system, *relptr);
1853 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1854 ASC_FREE(relname);
1855 break;
1856 }
1857 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1858 /* find the diagonal elements */
1859 for(j=0; j<count; ++j){
1860 if(var_deriv(variables[j])){
1861 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), c_j * derivatives[j]);
1862 }else{
1863 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), derivatives[j]);
1864 }
1865 }
1866 }
1867
1868 mtx_assemble(P);
1869
1870 if(status){
1871 CONSOLE_DEBUG("Error found when evaluating derivatives");
1872 res = 1; goto finish; /* recoverable */
1873 }
1874
1875 integrator_ida_write_incidence(sys);
1876
1877 res = 0;
1878 finish:
1879 ASC_FREE(variables);
1880 ASC_FREE(derivatives);
1881 return res;
1882 };
1883
1884 /**
1885 EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1886
1887 'solve' function.
1888 */
1889 static int integrator_ida_psolve_jacobian(realtype tt,
1890 N_Vector yy, N_Vector yp, N_Vector rr,
1891 N_Vector rvec, N_Vector zvec,
1892 realtype c_j, realtype delta, void *p_data,
1893 N_Vector tmp
1894 ){
1895 IntegratorSystem *sys;
1896 IntegratorIdaData *data;
1897 IntegratorIdaPrecDataJacobian *precdata;
1898 sys = (IntegratorSystem *)p_data;
1899 data = sys->enginedata;
1900 precdata = (IntegratorIdaPrecDataJacobian *)(data->precdata);
1901 linsolqr_system_t L = precdata->L;
1902
1903 linsolqr_add_rhs(L,NV_DATA_S(rvec),FALSE);
1904
1905 mtx_region_t R;
1906 R.row.low = R.col.low = 0;
1907 R.row.high = R.col.high = mtx_order(linsolqr_get_matrix(L)) - 1;
1908 linsolqr_set_region(L,R);
1909
1910 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
1911 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
1912
1913 /// @TODO more here
1914
1915 linsolqr_remove_rhs(L,NV_DATA_S(rvec));
1916
1917 CONSOLE_DEBUG("Solving Jacobian preconditioner (c_j = %f)",c_j);
1918 return 0;
1919 };
1920
1921
1922 /*----------------------------------------------
1923 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1924 */
1925
1926 static void integrator_ida_pcreate_jacobi(IntegratorSystem *sys){
1927 IntegratorIdaData *enginedata =sys->enginedata;
1928 IntegratorIdaPrecDataJacobi *precdata;
1929 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1930
1931 asc_assert(sys->n_y);
1932 precdata->PIii = N_VNew_Serial(sys->n_y);
1933
1934 enginedata->pfree = &integrator_ida_pfree_jacobi;
1935 enginedata->precdata = precdata;
1936 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1937 }
1938
1939 static void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1940 if(enginedata->precdata){
1941 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1942 N_VDestroy_Serial(precdata->PIii);
1943
1944 ASC_FREE(precdata);
1945 enginedata->precdata = NULL;
1946 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1947 }
1948 enginedata->pfree = NULL;
1949 }
1950
1951 /**
1952 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1953
1954 'setup' function.
1955 */
1956 static int integrator_ida_psetup_jacobi(realtype tt,
1957 N_Vector yy, N_Vector yp, N_Vector rr,
1958 realtype c_j, void *p_data,
1959 N_Vector tmp1, N_Vector tmp2,
1960 N_Vector tmp3
1961 ){
1962 int i, j, res;
1963 IntegratorSystem *sys;
1964 IntegratorIdaData *enginedata;
1965 IntegratorIdaPrecDataJacobi *precdata;
1966 struct rel_relation **relptr;
1967
1968 sys = (IntegratorSystem *)p_data;
1969 enginedata = sys->enginedata;
1970 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1971 double *derivatives;
1972 struct var_variable **variables;
1973 int count, status;
1974 char *relname;
1975
1976 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1977
1978 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1979 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1980
1981 /**
1982 @TODO FIXME here we are using the very inefficient and contorted approach
1983 of calculating the whole jacobian, then extracting just the diagonal elements.
1984 */
1985
1986 for(i=0, relptr = enginedata->rellist;
1987 i< enginedata->nrels && relptr != NULL;
1988 ++i, ++relptr
1989 ){
1990
1991 /* get derivatives for this particular relation */
1992 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1993 if(status){
1994 relname = rel_make_name(sys->system, *relptr);
1995 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1996 ASC_FREE(relname);
1997 break;
1998 }
1999 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
2000 /* find the diagonal elements */
2001 for(j=0; j<count; ++j){
2002 if(var_sindex(variables[j])==i){
2003 if(var_deriv(variables[j])){
2004 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
2005 }else{
2006 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
2007 }
2008
2009 }
2010 }
2011 #ifdef PREC_DEBUG
2012 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
2013 #endif
2014 }
2015
2016 if(status){
2017 CONSOLE_DEBUG("Error found when evaluating derivatives");
2018 res = 1; goto finish; /* recoverable */
2019 }
2020
2021 integrator_ida_write_incidence(sys);
2022
2023 res = 0;
2024 finish:
2025 ASC_FREE(variables);
2026 ASC_FREE(derivatives);
2027 return res;
2028 };
2029
2030 /**
2031 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
2032
2033 'solve' function.
2034 */
2035 static int integrator_ida_psolve_jacobi(realtype tt,
2036 N_Vector yy, N_Vector yp, N_Vector rr,
2037 N_Vector rvec, N_Vector zvec,
2038 realtype c_j, realtype delta, void *p_data,
2039 N_Vector tmp
2040 ){
2041 IntegratorSystem *sys;
2042 IntegratorIdaData *data;
2043 IntegratorIdaPrecDataJacobi *precdata;
2044 sys = (IntegratorSystem *)p_data;
2045 data = sys->enginedata;
2046 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
2047
2048 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
2049 N_VProd(precdata->PIii, rvec, zvec);
2050 return 0;
2051 };
2052
2053 /*----------------------------------------------
2054 STATS
2055 */
2056
2057 /**
2058 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
2059 status in a struct instead of separately.
2060
2061 @return IDA_SUCCESS on success.
2062 */
2063 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
2064
2065 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
2066
2067 int res;
2068
2069 /*
2070 There is an error in the documentation for this function in Sundials 2.2.
2071 According the the header file, the hinused stat is not provided.
2072 */
2073 res = IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups,
2074 &s->netfails, &s->qlast, &s->qcur, &s->hlast, &s->hcur,
2075 &s->tcur
2076 );
2077
2078 /* get the missing statistic */
2079 IDAGetActualInitStep(ida_mem, &s->hinused);
2080
2081 return res;
2082 #else
2083
2084 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
2085 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
2086 ,&s->hlast, &s->hcur, &s->tcur
2087 );
2088
2089 #endif
2090 }
2091
2092 /**
2093 This routine just outputs the stats to the CONSOLE_DEBUG routine.
2094
2095 @TODO provide a GUI way of stats reporting from IDA.
2096 */
2097 static void integrator_ida_write_stats(IntegratorIdaStats *stats){
2098 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
2099 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
2100 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
2101 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
2102 SI(qlast); SI(qcur);
2103 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
2104 # undef SL
2105 # undef SI
2106 # undef SR
2107 }
2108
2109 /*------------------------------------------------------------------------------
2110 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
2111 */
2112
2113 /**
2114 Here we construct the local transfer matrix. It's a bit of a naive
2115 approach; probably far more efficient ways can be worked out. But it will
2116 hopefully be a useful way to examine stability of some spatial
2117 discretisation schemes for PDAE systems.
2118
2119 http://ascendserver.cheme.cmu.edu/wiki/index.php/IDA#Stability
2120 */
2121 static int integrator_ida_transfer_matrix(const IntegratorSystem *sys, struct SystemJacobianStruct *J){
2122 int i=0, res;
2123 enum submat{II_GA=0, II_GD, II_FA, II_FD, II_FDP, II_NUM};
2124
2125 const var_filter_t *matvf[II_NUM] = {
2126 &system_vfilter_algeb
2127 ,&system_vfilter_diff
2128 ,&system_vfilter_algeb
2129 ,&system_vfilter_diff
2130 ,&system_vfilter_deriv
2131 };
2132
2133 const rel_filter_t *matrf[II_NUM] = {
2134 &system_rfilter_algeb
2135 ,&system_rfilter_algeb
2136 ,&system_rfilter_diff
2137 ,&system_rfilter_diff
2138 ,&system_rfilter_diff
2139 };
2140
2141 struct SystemJacobianStruct D[II_NUM];
2142
2143 for(i=0;i<II_NUM;++i){
2144 res = system_jacobian(sys->system, matrf[i], matvf[i], 1/*safe*/ ,&(D[i]));
2145 }
2146
2147 /* compute inverses for matrices that need it */
2148 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
2149 return 1;
2150 }
2151
2152
2153 /**
2154 Our task here is to write the matrices that IDA *should* be seeing. We
2155 are actually making calls to relman_diff in order to do that, so we're
2156 really going back to the variables in the actualy system and computing
2157 row by row what the values are. This should mean just a single call to
2158 each blackbox present in the system (if blackbox caching is working
2159 correctly).
2160 */
2161 static int integrator_ida_write_matrix(const IntegratorSystem *sys, FILE *f, const char *type){
2162 /* IntegratorIdaData *enginedata; */
2163 struct SystemJacobianStruct J = {NULL,NULL,NULL,0,0};
2164 int status=1;
2165 mtx_region_t R;
2166
2167 if(type==NULL)type = "dx'/dx";
2168
2169 if(0==strcmp(type,"dg/dz")){
2170 CONSOLE_DEBUG("Calculating dg/dz...");
2171 status = system_jacobian(sys->system
2172 , &system_rfilter_algeb, &system_vfilter_algeb
2173 , 1 /* safe */
2174 , &J
2175 );
2176 }else if(0==strcmp(type,"dg/dx")){
2177 CONSOLE_DEBUG("Calculating dg/dx...");
2178 status = system_jacobian(sys->system
2179 , &system_rfilter_algeb, &system_vfilter_diff
2180 , 1 /* safe */
2181 , &J
2182 );
2183 }else if(0==strcmp(type,"df/dx'")){
2184 CONSOLE_DEBUG("Calculating df/dx'...");
2185 status = system_jacobian(sys->system
2186 , &system_rfilter_diff, &system_vfilter_deriv
2187 , 1 /* safe */
2188 , &J
2189 );
2190 }else if(0==strcmp(type,"df/dz")){
2191 CONSOLE_DEBUG("Calculating df/dz...");
2192 status = system_jacobian(sys->system
2193 , &system_rfilter_diff, &system_vfilter_algeb
2194 , 1 /* safe */
2195 , &J
2196 );
2197 }else if(0==strcmp(type,"df/dx")){
2198 CONSOLE_DEBUG("Calculating df/dx...");
2199 status = system_jacobian(sys->system
2200 , &system_rfilter_diff, &system_vfilter_diff
2201 , 1 /* safe */
2202 , &J
2203 );
2204 }else if(0==strcmp(type,"dF/dy")){
2205 CONSOLE_DEBUG("Calculating dF/dy...");
2206 status = system_jacobian(sys->system
2207 , &system_rfilter_all, &system_vfilter_nonderiv
2208 , 1 /* safe */
2209 , &J
2210 );
2211 }else if(0==strcmp(type,"dF/dy'")){
2212 CONSOLE_DEBUG("Calculating dF/dy'...");
2213 status = system_jacobian(sys->system
2214 , &system_rfilter_all, &system_vfilter_deriv
2215 , 1 /* safe */
2216 , &J
2217 );
2218 }else if(0==strcmp(type,"dx'/dx")){
2219 /* system state transfer matrix dyd'/dyd */
2220 status = integrator_ida_transfer_matrix(sys, &J);
2221 }else{
2222 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid matrix type '%s'",type);
2223 return 1;
2224 }
2225
2226 if(status){
2227 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating matrix");
2228 }else{
2229 /* send the region explicitly, so that we handle non-square correctly */
2230 R.row.low = 0; R.col.low = 0;
2231 R.row.high = J.n_rels - 1; R.col.high = J.n_vars - 1;
2232 /* note that we're not fussy about empty matrices here... */
2233 mtx_write_region_mmio(f,J.M,&R);
2234 }
2235
2236 if(J.vars)ASC_FREE(J.vars);
2237 if(J.rels)ASC_FREE(J.rels);
2238 if(J.M)mtx_destroy(J.M);
2239
2240 return status;
2241 }
2242
2243 /**
2244 This routine outputs matrix structure in a crude text format, for the sake
2245 of debugging.
2246 */
2247 static void integrator_ida_write_incidence(IntegratorSystem *sys){
2248 int i, j;
2249 struct rel_relation **relptr;
2250 IntegratorIdaData *enginedata = sys->enginedata;
2251 double *derivatives;
2252 struct var_variable **variables;
2253 int count, status;
2254 char *relname;
2255
2256 if(enginedata->nrels > 100){
2257 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
2258 return;
2259 }
2260
2261 variables = ASC_NEW_ARRAY(struct var_variable *, sys->n_y * 2);
2262 derivatives = ASC_NEW_ARRAY(double, sys->n_y * 2);
2263
2264 CONSOLE_DEBUG("Outputting incidence information to console...");
2265
2266 for(i=0, relptr = enginedata->rellist;
2267 i< enginedata->nrels && relptr != NULL;
2268 ++i, ++relptr
2269 ){
2270 relname = rel_make_name(sys->system, *relptr);
2271
2272 /* get derivatives for this particular relation */
2273 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
2274 if(status){
2275 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
2276 ASC_FREE(relname);
2277 break;
2278 }
2279
2280 fprintf(stderr,"%3d:%-15s:",i,relname);
2281 ASC_FREE(relname);
2282
2283 for(j=0; j<count; ++j){
2284 if(var_deriv(variables[j])){
2285 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(sys,variables[j]));
2286 }else{
2287 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
2288 }
2289 }
2290 fprintf(stderr,"\n");
2291 }
2292 ASC_FREE(variables);
2293 ASC_FREE(derivatives);
2294 }
2295
2296 /* @return 0 on success */
2297 int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
2298 char *varname, *relname;
2299 struct var_variable **vlist, *var;
2300 struct rel_relation **rlist, *rel;
2301 long vlen, rlen;
2302 long i;
2303 long di;
2304
2305 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",sys->n_y);
2306
2307 /* if(integrator_sort_obs_vars(sys))return 10; */
2308
2309 if(sys->y && sys->ydot){
2310 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
2311 fprintf(fp,"index\t%-15s\tydot\n","y");
2312 fprintf(fp,"-----\t%-15s\t-----\n","-----");
2313 for(i=0;i<sys->n_y;++i){
2314 varname = var_make_name(sys->system, sys->y[i]);
2315 fprintf(fp,"%ld\t%-15s\t",i,varname);
2316 if(sys->ydot[i]){
2317 ASC_FREE(varname);
2318 varname = var_make_name(sys->system, sys->ydot[i]);
2319 fprintf(fp,"%s\n",varname);
2320 ASC_FREE(varname);
2321 }else{
2322 fprintf(fp,".\n");
2323 ASC_FREE(varname);
2324 }
2325 }
2326 }else{
2327 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
2328 }
2329
2330 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
2331 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
2332 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
2333
2334
2335 /* visit all the slv_system_t master var lists to collect vars */
2336 /* find the vars mostly in this one */
2337 vlist = slv_get_solvers_var_list(sys->system);
2338 vlen = slv_get_num_solvers_vars(sys->system);
2339 for(i=0;i<vlen;i++){
2340 var = vlist[i];
2341
2342 varname = var_make_name(sys->system, var);
2343 fprintf(fp,"%ld\t%-15s\t",i,varname);
2344
2345 if(var_fixed(var)){
2346 // it's fixed, so not really a DAE var
2347 fprintf(fp,"(fixed)\n");
2348 }else if(!var_active(var)){
2349 // inactive
2350 fprintf(fp,"(inactive)\n");
2351 }else if(!var_incident(var)){
2352 // not incident
2353 fprintf(fp,"(not incident)\n");
2354 }else{
2355 if(var_deriv(var)){
2356 if(sys->y_id){
2357 di = integrator_ida_diffindex1(sys,var);
2358 if(di>=0){
2359 ASC_FREE(varname);
2360 varname = var_make_name(sys->system,vlist[di]);
2361 fprintf(fp,".\tdiff(%ld='%s')\n",di,varname);
2362 }else{
2363 fprintf(fp,".\tdiff(???,err=%ld)\n",di);
2364 }
2365 }else{
2366 fprintf(fp,".\tderiv... of??\n");
2367 }
2368 }else{
2369 fprintf(fp,"%d\t.\n",var_sindex(var));
2370 }
2371 }
2372 ASC_FREE(varname);
2373 }
2374
2375 /* let's write out the relations too */
2376 rlist = slv_get_solvers_rel_list(sys->system);
2377 rlen = slv_get_num_solvers_rels(sys->system);
2378
2379 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
2380 fprintf(fp,"index\tname\n");
2381 fprintf(fp,"-----\t----\n");
2382 for(i=0; i<rlen; ++i){
2383 rel = rlist[i];
2384 relname = rel_make_name(sys->system,rel);
2385 fprintf(fp,"%ld\t%s\n",i,relname);
2386 ASC_FREE(relname);
2387 }
2388
2389 /* write out the derivative chains */
2390 fprintf(fp,"\nDERIVATIVE CHAINS\n");
2391 if(integrator_ida_analyse_debug(sys,stderr)){
2392 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
2393 return 340;
2394 }
2395 fprintf(fp,"\n");
2396
2397 /* and lets write block debug output */
2398 system_block_debug(sys->system, fp);
2399
2400 return 0; /* success */
2401 }
2402
2403 /*----------------------------------------------
2404 ERROR REPORTING
2405 */
2406 /**
2407 Error message reporter function to be passed to IDA. All error messages
2408 will trigger a call to this function, so we should find everything
2409 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
2410 panel (in the case of PyGTK).
2411 */
2412 static void integrator_ida_error(int error_code
2413 , const char *module, const char *function
2414 , char *msg, void *eh_data
2415 ){
2416 IntegratorSystem *sys;
2417 error_severity_t sev;
2418
2419 /* cast back the IntegratorSystem, just in case we need it */
2420 sys = (IntegratorSystem *)eh_data;
2421
2422 /* severity depends on the sign of the error_code value */
2423 if(error_code <= 0){
2424 sev = ASC_PROG_ERR;
2425 }else{
2426 sev = ASC_PROG_WARNING;
2427 }
2428
2429 /* use our all-purpose error reporting to get stuff back to the GUI */
2430 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
2431 }

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