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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2367 - (show annotations) (download) (as text)
Mon Jan 31 06:41:54 2011 UTC (13 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 72873 byte(s)
Add direction of boundary detection.
1 /* ASCEND modelling environment
2 Copyright (C) 2006-2011 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/general/platform.h>
74 #include <ascend/utilities/error.h>
75 #include <ascend/utilities/ascSignal.h>
76 #include <ascend/general/panic.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 *integ);
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 *integ, struct var_variable *var, const int *varindx);
257
258 /*static IntegratorVarVisitorFn integrator_dae_classify_var;
259 static void integrator_visit_system_vars(IntegratorSystem *integ,IntegratorVarVisitorFn *visitor);
260 static void integrator_dae_show_var(IntegratorSystem *integ, 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 *integ);
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 *integ);
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 *integ);
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 *integ){
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 integ->enginedata = (void *)enginedata;
340
341 integrator_ida_params_default(integ);
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 *integ){
368 IntegratorIdaData *d;
369 assert(integ!=NULL);
370 assert(integ->enginedata!=NULL);
371 assert(integ->engine==INTEG_IDA);
372 d = ((IntegratorIdaData *)(integ->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 integ->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 *integ){
407 asc_assert(integ!=NULL);
408 asc_assert(integ->engine==INTEG_IDA);
409 slv_parameters_t *p;
410 p = &(integ->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 typedef int IdaFlagFn(void *,int *);
544 typedef char *IdaFlagNameFn(int);
545
546 /* return 0 on success */
547 static int integrator_ida_solve(
548 IntegratorSystem *integ
549 , unsigned long start_index
550 , unsigned long finish_index
551 ){
552 void *ida_mem;
553 int flag, flag1, t_index;
554 realtype t0, reltol, abstol, t, tret, tout1;
555 N_Vector y0, yp0, abstolvect, ypret, yret, id;
556 IntegratorIdaData *enginedata;
557 char *linsolver;
558 int maxl;
559 IdaFlagFn *flagfn;
560 IdaFlagNameFn *flagnamefn;
561 const char *flagfntype;
562 char *pname = NULL;
563 struct rel_relation **rels;
564 int *rootsfound;
565 char havecrossed;
566
567 #ifdef SOLVE_DEBUG
568 char *varname, *relname;
569 #endif
570
571 int i,j,n_activerels,n_solverrels;
572 const IntegratorIdaPrec *prec = NULL;
573 int icopt; /* initial conditions strategy */
574
575 CONSOLE_DEBUG("STARTING IDA...");
576
577 enginedata = integrator_ida_enginedata(integ);
578
579 enginedata->safeeval = SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_SAFEEVAL);
580 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
581
582 /* store reference to list of relations (in enginedata) */
583 n_solverrels = slv_get_num_solvers_rels(integ->system);
584
585 n_activerels = slv_count_solvers_rels(integ->system, &integrator_ida_rel);
586
587 enginedata->bndlist = slv_get_solvers_bnd_list(integ->system);
588 enginedata->nbnds = slv_get_num_solvers_bnds(integ->system);
589
590 enginedata->rellist = ASC_NEW_ARRAY(struct rel_relation *, n_activerels);
591
592 rels = slv_get_solvers_rel_list(integ->system);
593
594 j=0;
595 for(i=0; i < n_solverrels; ++i){
596 if(rel_apply_filter(rels[i], &integrator_ida_rel)){
597 #ifdef SOLVE_DEBUG
598 relname = rel_make_name(integ->system, rels[i]);
599 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));
600 ASC_FREE(relname);
601 #endif
602 enginedata->rellist[j++]=rels[i];
603 }
604 }
605 asc_assert(j==n_activerels);
606
607 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
608 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
609
610 CONSOLE_DEBUG("Number of relations: %d",n_solverrels);
611 CONSOLE_DEBUG("Number of active relations: %d",n_activerels);
612 CONSOLE_DEBUG("Number of dependent vars: %d",integ->n_y);
613 CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);
614
615 enginedata->nrels = n_activerels;
616
617 if(enginedata->nrels != integ->n_y){
618 ERROR_REPORTER_HERE(ASC_USER_ERROR
619 ,"Integration problem is not square (%d active rels, %d vars)"
620 ,n_activerels, integ->n_y
621 );
622 return 1; /* failure */
623 }
624
625 #ifdef SOLVE_DEBUG
626 integrator_ida_debug(integ,stderr);
627 #endif
628
629 /* retrieve initial values from the system */
630
631 /** @TODO fix this, the starting time != first sample */
632 t0 = integrator_get_t(integ);
633 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
634
635 CONSOLE_DEBUG("RETRIEVING y0");
636
637 y0 = N_VNew_Serial(integ->n_y);
638 integrator_get_y(integ,NV_DATA_S(y0));
639
640 #ifdef SOLVE_DEBUG
641 CONSOLE_DEBUG("RETRIEVING yp0");
642 #endif
643
644 yp0 = N_VNew_Serial(integ->n_y);
645 integrator_get_ydot(integ,NV_DATA_S(yp0));
646
647 #ifdef SOLVE_DEBUG
648 N_VPrint_Serial(yp0);
649 CONSOLE_DEBUG("yp0 is at %p",&yp0);
650 #endif
651
652 /* create IDA object */
653 ida_mem = IDACreate();
654
655 /* relative error tolerance */
656 reltol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_RTOL);
657 CONSOLE_DEBUG("rtol = %8.2e",reltol);
658
659
660 /* allocate internal memory */
661 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
662 flag = IDAInit(ida_mem, &integrator_ida_fex, t0, y0 ,yp0);
663 #else
664 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
665 /* vector of absolute tolerances */
666 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
667 abstolvect = N_VNew_Serial(integ->n_y);
668 integrator_get_atol(integ,NV_DATA_S(abstolvect));
669
670 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
671
672 N_VDestroy_Serial(abstolvect);
673 }else{
674 /* scalar absolute tolerance (one value for all) */
675 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
676 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
677 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
678 }
679 #endif
680
681 if(flag==IDA_MEM_NULL){
682 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
683 return 2;
684 }else if(flag==IDA_MEM_FAIL){
685 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
686 return 3;
687 }else if(flag==IDA_ILL_INPUT){
688 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
689 return 4;
690 }/* else success */
691
692
693 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
694 CONSOLE_DEBUG("Assigning tolerances...");
695 /* assign tolerances */
696 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
697 CONSOLE_DEBUG("using vector of atol values");
698 abstolvect = N_VNew_Serial(integ->n_y);
699 integrator_get_atol(integ,NV_DATA_S(abstolvect));
700 IDASVtolerances(ida_mem, reltol, abstolvect);
701 N_VDestroy_Serial(abstolvect);
702 }else{
703 /* scalar tolerances */
704 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
705 CONSOLE_DEBUG("using scalar atol value = %8.2e",abstol);
706 IDASStolerances(ida_mem, reltol, abstol);
707 }
708 #endif
709
710 /* set optional inputs... */
711 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)integ);
712 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
713 IDASetUserData(ida_mem, (void *)integ);
714 #else
715 IDASetRdata(ida_mem, (void *)integ);
716 #endif
717 IDASetMaxStep(ida_mem, integrator_get_maxstep(integ));
718 IDASetInitStep(ida_mem, integrator_get_stepzero(integ));
719 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(integ));
720 if(integrator_get_minstep(integ)>0){
721 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
722 }
723
724 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
725 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
726
727 CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
728 IDASetMaxOrd(ida_mem,SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
729
730 /* there's no capability for setting *minimum* step size in IDA */
731
732
733 /* attach linear solver module, using the default value of maxl */
734 linsolver = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_LINSOLVER);
735 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
736 if(strcmp(linsolver,"ASCEND")==0){
737 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",integ->n_y);
738 IDAASCEND(ida_mem,integ->n_y);
739 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)integ);
740
741 flagfntype = "IDAASCEND";
742 flagfn = &IDAASCENDGetLastFlag;
743 flagnamefn = &IDAASCENDGetReturnFlagName;
744
745 }else if(strcmp(linsolver,"DENSE")==0){
746 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",integ->n_y);
747 flag = IDADense(ida_mem, integ->n_y);
748 switch(flag){
749 case IDADENSE_SUCCESS: break;
750 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 5;
751 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 5;
752 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 5;
753 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
754 }
755
756 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)){
757 CONSOLE_DEBUG("USING AUTODIFF");
758 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
759 flag = IDADlsSetDenseJacFn(ida_mem, &integrator_ida_djex);
760 #else
761 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)integ);
762 #endif
763 switch(flag){
764 case IDADENSE_SUCCESS: break;
765 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
766 }
767 }else{
768 CONSOLE_DEBUG("USING NUMERICAL DIFF");
769 }
770
771 flagfntype = "IDADENSE";
772 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
773 flagfn = &IDADlsGetLastFlag;
774 flagnamefn = &IDADlsGetReturnFlagName;
775 #else
776 flagfn = &IDADenseGetLastFlag;
777 flagnamefn = &IDADenseGetReturnFlagName;
778 #endif
779 }else{
780 /* remaining methods are all SPILS */
781 CONSOLE_DEBUG("IDA SPILS");
782
783 maxl = SLV_PARAM_INT(&(integ->params),IDA_PARAM_MAXL);
784 CONSOLE_DEBUG("maxl = %d",maxl);
785
786 /* what preconditioner? */
787 pname = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_PREC);
788 if(strcmp(pname,"NONE")==0){
789 prec = NULL;
790 }else if(strcmp(pname,"JACOBI")==0){
791 prec = &prec_jacobi;
792 }else{
793 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
794 return 7;
795 }
796
797 /* which SPILS linear solver? */
798 if(strcmp(linsolver,"SPGMR")==0){
799 CONSOLE_DEBUG("IDA SPGMR");
800 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
801 }else if(strcmp(linsolver,"SPBCG")==0){
802 CONSOLE_DEBUG("IDA SPBCG");
803 flag = IDASpbcg(ida_mem, maxl);
804 }else if(strcmp(linsolver,"SPTFQMR")==0){
805 CONSOLE_DEBUG("IDA SPTFQMR");
806 flag = IDASptfqmr(ida_mem,maxl);
807 }else{
808 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
809 return 8;
810 }
811
812 if(prec){
813 /* assign the preconditioner to the linear solver */
814 (prec->pcreate)(integ);
815 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
816 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve);
817 #else
818 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)integ);
819 #endif
820 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
821 }else{
822 CONSOLE_DEBUG("No preconditioner");
823 }
824
825 flagfntype = "IDASPILS";
826 flagfn = &IDASpilsGetLastFlag;
827 flagnamefn = &IDASpilsGetReturnFlagName;
828
829 if(flag==IDASPILS_MEM_NULL){
830 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
831 return 9;
832 }else if(flag==IDASPILS_MEM_FAIL){
833 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
834 return 9;
835 }/* else success */
836
837 /* assign the J*v function */
838 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)){
839 CONSOLE_DEBUG("USING AUTODIFF");
840 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
841 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex);
842 #else
843 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)integ);
844 #endif
845 if(flag==IDASPILS_MEM_NULL){
846 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
847 return 10;
848 }else if(flag==IDASPILS_LMEM_NULL){
849 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
850 return 10;
851 }/* else success */
852 }else{
853 CONSOLE_DEBUG("USING NUMERICAL DIFF");
854 }
855
856 if(strcmp(linsolver,"SPGMR")==0){
857 /* select Gram-Schmidt orthogonalisation */
858 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_GSMODIFIED)){
859 CONSOLE_DEBUG("USING MODIFIED GS");
860 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
861 if(flag!=IDASPILS_SUCCESS){
862 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
863 return 11;
864 }
865 }else{
866 CONSOLE_DEBUG("USING CLASSICAL GS");
867 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
868 if(flag!=IDASPILS_SUCCESS){
869 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
870 return 11;
871 }
872 }
873 }
874 }
875
876 /* set linear solver optional inputs...
877 ...nothing here at the moment...
878 */
879
880 /* calculate initial conditions */
881 icopt = 0;
882 if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"Y")==0){
883 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
884 icopt = IDA_Y_INIT;
885 asc_assert(icopt!=0);
886 }else if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
887 CONSOLE_DEBUG("Solving initial conditions using values of yd");
888 icopt = IDA_YA_YDP_INIT;
889 asc_assert(icopt!=0);
890 id = N_VNew_Serial(integ->n_y);
891 for(i=0; i < integ->n_y; ++i){
892 if(integ->ydot[i] == NULL){
893 NV_Ith_S(id,i) = 0.0;
894 #ifdef SOLVE_DEBUG
895 varname = var_make_name(integ->system,integ->y[i]);
896 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
897 ASC_FREE(varname);
898 #endif
899 }else{
900 #ifdef SOLVE_DEBUG
901 CONSOLE_DEBUG("y[%d] is differential",i);
902 #endif
903 NV_Ith_S(id,i) = 1.0;
904 }
905 }
906 IDASetId(ida_mem, id);
907 N_VDestroy_Serial(id);
908 }else if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC),"NONE")==0){
909 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
910 }else{
911 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
912 }
913
914 if(icopt){
915 integ->currentstep=0;
916 t_index=start_index + 1;
917 tout1 = samplelist_get(integ->samples, t_index);
918
919 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
920
921 #ifdef ASC_SIGNAL_TRAPS
922 /* catch SIGFPE if desired to */
923 if(enginedata->safeeval){
924 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
925 Asc_SignalHandlerPush(SIGFPE,SIG_DFL);
926 }else{
927 # ifdef FEX_DEBUG
928 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
929 # endif
930 Asc_SignalHandlerPushDefault(SIGFPE);
931 }
932 if (setjmp(g_fpe_env)==0) {
933 #endif
934
935 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=3
936 flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3 */
937 # else
938 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
939 # endif
940 /* check flags and output status */
941 switch(flag){
942 case IDA_SUCCESS:
943 CONSOLE_DEBUG("Initial conditions solved OK");
944 break;
945
946 case IDA_LSETUP_FAIL:
947 case IDA_LINIT_FAIL:
948 case IDA_LSOLVE_FAIL:
949 case IDA_NO_RECOVERY:
950 flag1 = -999;
951 flag = (flagfn)(ida_mem,&flag1);
952 if(flag){
953 ERROR_REPORTER_HERE(ASC_PROG_ERR
954 ,"Unable to retrieve error code from %s (err %d)"
955 ,flagfntype,flag
956 );
957 return 12;
958 }
959 ERROR_REPORTER_HERE(ASC_PROG_ERR
960 ,"%s returned flag '%s' (value = %d)"
961 ,flagfntype,(flagnamefn)(flag1),flag1
962 );
963 return 12;
964
965 default:
966 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
967 return 12;
968 }
969 #ifdef ASC_SIGNAL_TRAPS
970 }else{
971 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
972 return 13;
973 }
974
975 if(enginedata->safeeval){
976 Asc_SignalHandlerPop(SIGFPE,SIG_DFL);
977 }else{
978 CONSOLE_DEBUG("pop...");
979 Asc_SignalHandlerPopDefault(SIGFPE);
980 CONSOLE_DEBUG("...pop");
981 }
982 #endif
983 }/* icopt */
984
985 /* optionally, specify ROO-FINDING problem */
986 if(enginedata->nbnds){
987 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
988 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn);
989 #else
990 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)integ);
991 #endif
992 }
993
994 /* -- set up the IntegratorReporter */
995 integrator_output_init(integ);
996
997 /* -- store the initial values of all the stuff */
998 integrator_output_write(integ);
999 integrator_output_write_obs(integ);
1000
1001 /* specify where the returned values should be stored */
1002 yret = y0;
1003 ypret = yp0;
1004
1005 /* advance solution in time, return values as yret and derivatives as ypret */
1006 integ->currentstep=1;
1007 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++integ->currentstep){
1008 t = samplelist_get(integ->samples, t_index);
1009 t0 = integrator_get_t(integ);
1010 asc_assert(t > t0);
1011
1012 #ifdef SOLVE_DEBUG
1013 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
1014 #endif
1015
1016 #ifdef ASC_SIGNAL_TRAPS
1017 Asc_SignalHandlerPushDefault(SIGINT);
1018 if(setjmp(g_int_env)==0) {
1019 #endif
1020 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
1021 #ifdef ASC_SIGNAL_TRAPS
1022 }else{
1023 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Caught interrupt");
1024 flag = -555;
1025 }
1026 Asc_SignalHandlerPopDefault(SIGINT);
1027 #endif
1028
1029 /* this seems to work, so we can use it to avoid needing 'havecrossed'? */
1030 if(flag == IDA_ROOT_RETURN){
1031 CONSOLE_DEBUG("IDA reports root found!");
1032 }
1033
1034 /* so we will check for roots found explicitly */
1035 if(enginedata->nbnds){
1036 rootsfound = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
1037 havecrossed = 0;
1038 if(IDA_SUCCESS == IDAGetRootInfo(ida_mem, rootsfound)){
1039 for(i=0; i < enginedata->nbnds; ++i){
1040 if(rootsfound[i]){
1041 havecrossed = 1;
1042 #ifdef SOLVE_DEBUG
1043 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
1044 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary '%s' crossed%s",relname,rootsfound[i]>0?" (increasing)":" (decreasing)");
1045 ASC_FREE(relname);
1046 #else
1047 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Boundary %d crossed!%s", i, ,rootsfound[i]>0?" (increasing)":" (decreasing)");
1048 #endif
1049 }
1050 }
1051
1052 /* so, now we need to restart the integration. we will assume that
1053 everything changes: number of variables, etc, etc, etc. */
1054
1055 if(havecrossed){
1056 CONSOLE_DEBUG("Boundaries were crossed; need to reinitialise solver...");
1057 /** try resetting the boundary states now? */
1058 //IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)integ);
1059
1060 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
1061 IDAReInit(ida_mem, tret, yret, ypret);
1062 #elif SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR<3
1063 /* TODO find out what version needs the following line... not sure. */
1064 //IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret);
1065
1066 // FIXME this stuff has not been tested yet, and is very incomplete.
1067
1068 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
1069 /* vector of absolute tolerances */
1070 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
1071 abstolvect = N_VNew_Serial(integ->n_y);
1072 integrator_get_atol(integ,NV_DATA_S(abstolvect));
1073 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SV, reltol, abstolvect);
1074 N_VDestroy_Serial(abstolvect);
1075 }else{
1076 /* scalar absolute tolerance (one value for all) */
1077 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
1078 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
1079 flag = IDAReInit(ida_mem, &integrator_ida_fex, tret, yret, ypret, IDA_SS, reltol, &abstol);
1080 }
1081 #else
1082 /* allocate internal memory */
1083 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)){
1084 /* vector of absolute tolerances */
1085 abstolvect = N_VNew_Serial(integ->n_y);
1086 integrator_get_atol(integ,NV_DATA_S(abstolvect));
1087 if(IDA_SUCCESS != IDAReInit(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect)){
1088 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
1089 }
1090 N_VDestroy_Serial(abstolvect);
1091 }else{
1092 /* scalar absolute tolerance (one value for all) */
1093 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
1094 if(IDA_SUCCESS != IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol)){
1095 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to reinitialise IDA");
1096 }
1097 }
1098 #endif
1099 }
1100 }else{
1101 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");
1102 }
1103 ASC_FREE(rootsfound);
1104 }
1105
1106
1107 /* pass the values of everything back to the compiler */
1108 integrator_set_t(integ, (double)tret);
1109 integrator_set_y(integ, NV_DATA_S(yret));
1110 integrator_set_ydot(integ, NV_DATA_S(ypret));
1111
1112 if(flag<0){
1113 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
1114 break;
1115 }
1116
1117 /* -- do something so that integ knows the values of tret, yret and ypret */
1118
1119 /* -- store the current values of all the stuff */
1120 integrator_output_write(integ);
1121 integrator_output_write_obs(integ);
1122
1123 }/* loop through next sample timestep */
1124
1125 /* -- close the IntegratorReporter */
1126 integrator_output_close(integ);
1127
1128 /* get optional outputs */
1129 #ifdef STATS_DEBUG
1130 IntegratorIdaStats stats;
1131 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
1132 integrator_ida_write_stats(&stats);
1133 }else{
1134 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
1135 }
1136 #endif
1137
1138 /* free solution memory */
1139 N_VDestroy_Serial(yret);
1140 N_VDestroy_Serial(ypret);
1141
1142 /* free solver memory */
1143 IDAFree(ida_mem);
1144
1145 if(flag < -500){
1146 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", t);
1147 return -flag;
1148 }
1149
1150 if(flag < 0){
1151 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1152 return 14;
1153 }
1154
1155 /* all done, success */
1156 return 0;
1157 }
1158
1159 /*--------------------------------------------------
1160 RESIDUALS AND JACOBIAN AND IDAROOTFN
1161 */
1162
1163 #if 0
1164 typedef void (SignalHandlerFn)(int);
1165 SignalHandlerFn integrator_ida_sig;
1166 SignalHandlerFn *integrator_ida_sig_old;
1167 jmp_buf integrator_ida_jmp_buf;
1168 fenv_t integrator_ida_fenv_old;
1169
1170
1171 void integrator_ida_write_feinfo(){
1172 int f;
1173 f = fegetexcept();
1174 CONSOLE_DEBUG("Locating nature of exception...");
1175 if(f & FE_DIVBYZERO)ERROR_REPORTER_HERE(ASC_PROG_ERR,"DIV BY ZERO");
1176 if(f & FE_INEXACT)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INEXACT");
1177 if(f & FE_INVALID)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INVALID");
1178 if(f & FE_OVERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"OVERFLOW");
1179 if(f & FE_UNDERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"UNDERFLOW");
1180 if(f==0)ERROR_REPORTER_HERE(ASC_PROG_ERR,"FLAGS ARE CLEAR?!?");
1181 }
1182
1183 void integrator_ida_sig(int sig){
1184 /* the wrong signal: rethrow to the default handler */
1185 if(sig!=SIGFPE){
1186 signal(SIGFPE,SIG_DFL);
1187 raise(sig);
1188 }
1189 integrator_ida_write_feinfo();
1190 CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);
1191 longjmp(integrator_ida_jmp_buf,sig);
1192 }
1193 #endif
1194
1195 /**
1196 Function to evaluate system residuals, in the form required for IDA.
1197
1198 Given tt, yy and yp, we need to evaluate and return rr.
1199
1200 @param tt current value of indep variable (time)
1201 @param yy current values of dependent variable vector
1202 @param yp current values of derivatives of dependent variables
1203 @param rr the output residual vector (is we're returning data to)
1204 @param res_data pointer to our stuff (integ in this case).
1205
1206 @return 0 on success, positive on recoverable error, and
1207 negative on unrecoverable error.
1208 */
1209 static int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1210 IntegratorSystem *integ;
1211 IntegratorIdaData *enginedata;
1212 int i, calc_ok, is_error;
1213 struct rel_relation** relptr;
1214 double resid;
1215 char *relname;
1216 #ifdef FEX_DEBUG
1217 char *varname;
1218 char diffname[30];
1219 #endif
1220
1221 integ = (IntegratorSystem *)res_data;
1222 enginedata = integrator_ida_enginedata(integ);
1223
1224 #ifdef FEX_DEBUG
1225 /* fprintf(stderr,"\n\n"); */
1226 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1227 #endif
1228
1229 if(NV_LENGTH_S(rr)!=enginedata->nrels){
1230 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
1231 return -1; /* unrecoverable */
1232 }
1233
1234 /* pass the values of everything back to the compiler */
1235 integrator_set_t(integ, (double)tt);
1236 integrator_set_y(integ, NV_DATA_S(yy));
1237 integrator_set_ydot(integ, NV_DATA_S(yp));
1238
1239 /* perform bounds checking on all variables */
1240 if(slv_check_bounds(integ->system, 0, -1, NULL)){
1241 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1242 return 1;
1243 }
1244
1245 /* evaluate each residual in the rellist */
1246 is_error = 0;
1247 relptr = enginedata->rellist;
1248
1249
1250 #ifdef ASC_SIGNAL_TRAPS
1251 if(enginedata->safeeval){
1252 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1253 }else{
1254 # ifdef FEX_DEBUG
1255 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1256 # endif
1257 Asc_SignalHandlerPushDefault(SIGFPE);
1258 }
1259
1260 if (SETJMP(g_fpe_env)==0) {
1261 #endif
1262
1263
1264 for(i=0, relptr = enginedata->rellist;
1265 i< enginedata->nrels && relptr != NULL;
1266 ++i, ++relptr
1267 ){
1268 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
1269
1270 NV_Ith_S(rr,i) = resid;
1271 if(!calc_ok){
1272 relname = rel_make_name(integ->system, *relptr);
1273 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1274 ASC_FREE(relname);
1275 /* presumable some output already made? */
1276 is_error = 1;
1277 }/*else{
1278 CONSOLE_DEBUG("Calc OK");
1279 }*/
1280 }
1281
1282 if(!is_error){
1283 for(i=0;i< enginedata->nrels; ++i){
1284 if(isnan(NV_Ith_S(rr,i))){
1285 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in residual %d",i);
1286 is_error=1;
1287 }
1288 }
1289 #ifdef FEX_DEBUG
1290 if(!is_error){
1291 CONSOLE_DEBUG("No NAN detected");
1292 }
1293 #endif
1294 }
1295
1296 #ifdef ASC_SIGNAL_TRAPS
1297 }else{
1298 relname = rel_make_name(integ->system, *relptr);
1299 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1300 ASC_FREE(relname);
1301 is_error = 1;
1302 }
1303
1304 if(enginedata->safeeval){
1305 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1306 }else{
1307 Asc_SignalHandlerPopDefault(SIGFPE);
1308 }
1309 #endif
1310
1311
1312 #ifdef FEX_DEBUG
1313 /* output residuals to console */
1314 CONSOLE_DEBUG("RESIDUAL OUTPUT");
1315 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1316 for(i=0; i<integ->n_y; ++i){
1317 varname = var_make_name(integ->system,integ->y[i]);
1318 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1319 if(integ->ydot[i]){
1320 varname = var_make_name(integ->system,integ->ydot[i]);
1321 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1322 }else{
1323 snprintf(diffname,99,"diff(%s)",varname);
1324 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1325 }
1326 ASC_FREE(varname);
1327 relname = rel_make_name(integ->system,enginedata->rellist[i]);
1328 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1329 }
1330 #endif
1331
1332 if(is_error){
1333 return 1;
1334 }
1335
1336 #ifdef FEX_DEBUG
1337 CONSOLE_DEBUG("RESIDUAL OK");
1338 #endif
1339 return 0;
1340 }
1341
1342 /**
1343 Dense Jacobian evaluation. Only suitable for small problems!
1344 Has been seen working for problems up to around 2000 vars, FWIW.
1345 */
1346 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
1347 static int integrator_ida_djex(int Neq, realtype tt, realtype c_j
1348 , N_Vector yy, N_Vector yp, N_Vector rr
1349 , IDA_MTX_T Jac, void *jac_data
1350 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1351 ){
1352 #else
1353 static int integrator_ida_djex(long int Neq, realtype tt
1354 , N_Vector yy, N_Vector yp, N_Vector rr
1355 , realtype c_j, void *jac_data, IDA_MTX_T Jac
1356 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1357 ){
1358 #endif
1359 IntegratorSystem *integ;
1360 IntegratorIdaData *enginedata;
1361 char *relname;
1362 #ifdef DJEX_DEBUG
1363 struct var_variable **varlist;
1364 char *varname;
1365 #endif
1366 struct rel_relation **relptr;
1367 int i;
1368 double *derivatives;
1369 struct var_variable **variables;
1370 int count, j;
1371 int status, is_error = 0;
1372
1373 integ = (IntegratorSystem *)jac_data;
1374 enginedata = integrator_ida_enginedata(integ);
1375
1376 /* allocate space for returns from relman_diff3 */
1377 /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
1378 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1379 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1380
1381 /* pass the values of everything back to the compiler */
1382 integrator_set_t(integ, (double)tt);
1383 integrator_set_y(integ, NV_DATA_S(yy));
1384 integrator_set_ydot(integ, NV_DATA_S(yp));
1385
1386 /* perform bounds checking on all variables */
1387 if(slv_check_bounds(integ->system, 0, -1, NULL)){
1388 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1389 return 1;
1390 }
1391
1392 #ifdef DJEX_DEBUG
1393 varlist = slv_get_solvers_var_list(integ->system);
1394
1395 /* print vars */
1396 for(i=0; i < integ->n_y; ++i){
1397 varname = var_make_name(integ->system, integ->y[i]);
1398 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1399 asc_assert(NV_Ith_S(yy,i) == var_value(integ->y[i]));
1400 ASC_FREE(varname);
1401 }
1402
1403 /* print derivatives */
1404 for(i=0; i < integ->n_y; ++i){
1405 if(integ->ydot[i]){
1406 varname = var_make_name(integ->system, integ->ydot[i]);
1407 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(integ->ydot[i]));
1408 ASC_FREE(varname);
1409 }else{
1410 varname = var_make_name(integ->system, integ->y[i]);
1411 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1412 ASC_FREE(varname);
1413 }
1414 }
1415
1416 /* print step size */
1417 CONSOLE_DEBUG("<c_j> = %g",c_j);
1418 #endif
1419
1420 /* build up the dense jacobian matrix... */
1421 status = 0;
1422 for(i=0, relptr = enginedata->rellist;
1423 i< enginedata->nrels && relptr != NULL;
1424 ++i, ++relptr
1425 ){
1426 /* get derivatives for this particular relation */
1427 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1428
1429 if(status){
1430 relname = rel_make_name(integ->system, *relptr);
1431 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1432 ASC_FREE(relname);
1433 is_error = 1;
1434 break;
1435 }
1436
1437 /* output what's going on here ... */
1438 #ifdef DJEX_DEBUG
1439 relname = rel_make_name(integ->system, *relptr);
1440 fprintf(stderr,"%d: '%s': ",i,relname);
1441 for(j=0;j<count;++j){
1442 varname = var_make_name(integ->system, variables[j]);
1443 if(var_deriv(variables[j])){
1444 fprintf(stderr," '%s'=",varname);
1445 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(integ,variables[j]));
1446 }else{
1447 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1448 }
1449 ASC_FREE(varname);
1450 }
1451 /* relname is freed further down */
1452 fprintf(stderr,"\n");
1453 #endif
1454
1455 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1456 for(j=0; j < count; ++j){
1457 #ifdef DJEX_DEBUG
1458 varname = var_make_name(integ->system,variables[j]);
1459 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1460 ASC_FREE(varname);
1461 #endif
1462 if(!var_deriv(variables[j])){
1463 #ifdef DJEX_DEBUG
1464 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
1465 asc_assert(var_sindex(variables[j]) >= 0);
1466 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1467 #endif
1468 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1469 }else{
1470 DENSE_ELEM(Jac,i,integrator_ida_diffindex(integ,variables[j])) += derivatives[j] * c_j;
1471 #ifdef DJEX_DEBUG
1472 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1473 #endif
1474 }
1475 }
1476 }
1477
1478 #ifdef DJEX_DEBUG
1479 ASC_FREE(relname);
1480 CONSOLE_DEBUG("PRINTING JAC");
1481 fprintf(stderr,"\t");
1482 for(j=0; j < integ->n_y; ++j){
1483 if(j)fprintf(stderr,"\t");
1484 varname = var_make_name(integ->system,integ->y[j]);
1485 fprintf(stderr,"%11s",varname);
1486 ASC_FREE(varname);
1487 }
1488 fprintf(stderr,"\n");
1489 for(i=0; i < enginedata->nrels; ++i){
1490 relname = rel_make_name(integ->system, enginedata->rellist[i]);
1491 fprintf(stderr,"%s\t",relname);
1492 ASC_FREE(relname);
1493
1494 for(j=0; j < integ->n_y; ++j){
1495 if(j)fprintf(stderr,"\t");
1496 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1497 }
1498 fprintf(stderr,"\n");
1499 }
1500 #endif
1501
1502 /* test for NANs */
1503 if(!is_error){
1504 for(i=0;i< enginedata->nrels; ++i){
1505 for(j=0;j<integ->n_y;++j){
1506 if(isnan(DENSE_ELEM(Jac,i,j))){
1507 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
1508 is_error=1;
1509 }
1510 }
1511 }
1512 #ifdef DJEX_DEBUG
1513 if(!is_error){
1514 CONSOLE_DEBUG("No NAN detected");
1515 }
1516 #endif
1517 }
1518
1519 /* if(integrator_ida_check_diffindex(integ)){
1520 is_error = 1;
1521 }*/
1522
1523 if(is_error){
1524 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1525 return 1;
1526 }
1527
1528 #ifdef DJEX_DEBUG
1529 CONSOLE_DEBUG("DJEX RETURNING 0");
1530 /* ASC_PANIC("Quitting"); */
1531 #endif
1532 return 0;
1533 }
1534
1535 /**
1536 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1537
1538 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1539
1540 @param tt current value of the independent variable (time, t)
1541 @param yy current value of the dependent variable vector, y(t).
1542 @param yp current value of y'(t).
1543 @param rr current value of the residual vector F(t, y, y').
1544 @param v the vector by which the Jacobian must be multiplied to the right.
1545 @param Jv the output vector computed
1546 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1547 @param jac_data pointer to our stuff (integ in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1548 @param tmp1 @see tmp2
1549 @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.
1550 @return 0 on success
1551 */
1552 static int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1553 , N_Vector v, N_Vector Jv, realtype c_j
1554 , void *jac_data, N_Vector tmp1, N_Vector tmp2
1555 ){
1556 IntegratorSystem *integ;
1557 IntegratorIdaData *enginedata;
1558 int i, j, is_error=0;
1559 struct rel_relation** relptr = 0;
1560 char *relname;
1561 int status;
1562 double Jv_i;
1563
1564 struct var_variable **variables;
1565 double *derivatives;
1566 int count;
1567 struct var_variable **varlist;
1568 #ifdef JEX_DEBUG
1569
1570 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1571 #endif
1572
1573 integ = (IntegratorSystem *)jac_data;
1574 enginedata = integrator_ida_enginedata(integ);
1575 varlist = slv_get_solvers_var_list(integ->system);
1576
1577 /* pass the values of everything back to the compiler */
1578 integrator_set_t(integ, (double)tt);
1579 integrator_set_y(integ, NV_DATA_S(yy));
1580 integrator_set_ydot(integ, NV_DATA_S(yp));
1581 /* no real use for residuals (rr) here, I don't think? */
1582
1583 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1584
1585 i = NV_LENGTH_S(yy) * 2;
1586 #ifdef JEX_DEBUG
1587 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1588 #endif
1589 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1590 derivatives = ASC_NEW_ARRAY(double, i);
1591
1592 /* evaluate the derivatives... */
1593 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1594
1595 #ifdef ASC_SIGNAL_TRAPS
1596 Asc_SignalHandlerPushDefault(SIGFPE);
1597 if (SETJMP(g_fpe_env)==0) {
1598 #endif
1599 for(i=0, relptr = enginedata->rellist;
1600 i< enginedata->nrels && relptr != NULL;
1601 ++i, ++relptr
1602 ){
1603 /* get derivatives for this particular relation */
1604 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1605 #ifdef JEX_DEBUG
1606 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1607 #endif
1608
1609 if(status){
1610 relname = rel_make_name(integ->system, *relptr);
1611 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1612 ASC_FREE(relname);
1613 is_error = 1;
1614 break;
1615 }
1616
1617 /*
1618 Now we have the derivatives wrt each alg/diff variable in the
1619 present equation. variables[] points into the varlist. need
1620 a mapping from the varlist to the y and ydot lists.
1621 */
1622
1623 Jv_i = 0;
1624 for(j=0; j < count; ++j){
1625 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], integ->n_y);
1626 varname = var_make_name(integ->system, enginedata->varlist[variables[j]]);
1627 if(varname){
1628 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1629 ASC_FREE(varname);
1630 }else{
1631 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1632 }
1633 */
1634
1635 /* we don't calculate derivatives wrt indep var */
1636 asc_assert(variables[j]>=0);
1637 if(variables[j] == integ->x) continue;
1638 #ifdef JEX_DEBUG
1639 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1640 #endif
1641 if(var_deriv(variables[j])){
1642 #define DIFFINDEX integrator_ida_diffindex(integ,variables[j])
1643 #ifdef JEX_DEBUG
1644 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
1645 , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
1646 , i, DIFFINDEX, derivatives[j]
1647 , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
1648 );
1649 #endif
1650 asc_assert(integ->ydot[DIFFINDEX]==variables[j]);
1651 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1652 #undef DIFFINDEX
1653 }else{
1654 #define VARINDEX var_sindex(variables[j])
1655 #ifdef JEX_DEBUG
1656 asc_assert(integ->y[VARINDEX]==variables[j]);
1657 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
1658 , i
1659 , derivatives[j] * NV_Ith_S(v,VARINDEX)
1660 , i, VARINDEX, derivatives[j]
1661 , VARINDEX, NV_Ith_S(v,VARINDEX)
1662 );
1663 #endif
1664 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1665 #undef VARINDEX
1666 }
1667 }
1668
1669 NV_Ith_S(Jv,i) = Jv_i;
1670 #ifdef JEX_DEBUG
1671 CONSOLE_DEBUG("rel = %p",*relptr);
1672 relname = rel_make_name(integ->system, *relptr);
1673 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1674 ASC_FREE(relname);
1675 return 1;
1676 #endif
1677 }
1678 #ifdef ASC_SIGNAL_TRAPS
1679 }else{
1680 relname = rel_make_name(integ->system, *relptr);
1681 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1682 ASC_FREE(relname);
1683 is_error = 1;
1684 }
1685 Asc_SignalHandlerPopDefault(SIGFPE);
1686 #endif
1687
1688 if(is_error){
1689 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1690 return 1;
1691 }
1692 return 0;
1693 }
1694
1695 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1696 static int integrator_ida_sjex(long int Neq, realtype tt
1697 , N_Vector yy, N_Vector yp, N_Vector rr
1698 , realtype c_j, void *jac_data, mtx_matrix_t Jac
1699 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1700 ){
1701 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1702 return -1;
1703 }
1704
1705 /* root finding function */
1706
1707 int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data){
1708 IntegratorSystem *integ;
1709 IntegratorIdaData *enginedata;
1710 int i;
1711 #ifdef ROOT_DEBUG
1712 char *relname;
1713 #endif
1714
1715 asc_assert(g_data!=NULL);
1716 integ = (IntegratorSystem *)g_data;
1717 enginedata = integrator_ida_enginedata(integ);
1718
1719 /* pass the values of everything back to the compiler */
1720 integrator_set_t(integ, (double)tt);
1721 integrator_set_y(integ, NV_DATA_S(yy));
1722 integrator_set_ydot(integ, NV_DATA_S(yp));
1723
1724 asc_assert(gout!=NULL);
1725
1726 #ifdef ROOT_DEBUG
1727 CONSOLE_DEBUG("t = %f",tt);
1728 #endif
1729
1730 /* evaluate the residuals for each of the boundaries */
1731 for(i=0; i < enginedata->nbnds; ++i){
1732 switch(bnd_kind(enginedata->bndlist[i])){
1733 case e_bnd_rel: /* real-valued boundary relation */
1734 gout[i] = bndman_real_eval(enginedata->bndlist[i]);
1735 #ifdef ROOT_DEBUG
1736 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
1737 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
1738 ASC_FREE(relname);
1739 #endif
1740 break;
1741 case e_bnd_logrel:
1742 if(bndman_log_eval(enginedata->bndlist[i])){
1743 CONSOLE_DEBUG("bnd[%d] = TRUE",i);
1744 #ifdef ROOT_DEBUG
1745 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
1746 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
1747 ASC_FREE(relname);
1748 #endif
1749 gout[i] = +1.0;
1750 }else{
1751 CONSOLE_DEBUG("bnd[%d] = FALSE",i);
1752 gout[i] = -1.0;
1753 }
1754 break;
1755 case e_bnd_undefined:
1756 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid boundary type e_bnd_undefined");
1757 return 1;
1758 }
1759 }
1760
1761 return 0; /* no way to detect errors in bndman_*_eval at this stage */
1762 }
1763
1764
1765 /*----------------------------------------------
1766 FULL JACOBIAN PRECONDITIONER -- EXPERIMENTAL.
1767 */
1768
1769 static void integrator_ida_pcreate_jacobian(IntegratorSystem *integ){
1770 IntegratorIdaData *enginedata =integ->enginedata;
1771 IntegratorIdaPrecDataJacobian *precdata;
1772 precdata = ASC_NEW(IntegratorIdaPrecDataJacobian);
1773 mtx_matrix_t P;
1774 asc_assert(integ->n_y);
1775 precdata->L = linsolqr_create_default();
1776
1777 /* allocate matrix to be used by linsolqr */
1778 P = mtx_create();
1779 mtx_set_order(P, integ->n_y);
1780 linsolqr_set_matrix(precdata->L, P);
1781
1782 enginedata->pfree = &integrator_ida_pfree_jacobian;
1783 enginedata->precdata = precdata;
1784 CONSOLE_DEBUG("Allocated memory for Full Jacobian preconditioner");
1785 }
1786
1787 static void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata){
1788 mtx_matrix_t P;
1789 IntegratorIdaPrecDataJacobian *precdata;
1790
1791 if(enginedata->precdata){
1792 precdata = (IntegratorIdaPrecDataJacobian *)enginedata->precdata;
1793 P = linsolqr_get_matrix(precdata->L);
1794 mtx_destroy(P);
1795 linsolqr_destroy(precdata->L);
1796 ASC_FREE(precdata);
1797 enginedata->precdata = NULL;
1798
1799 CONSOLE_DEBUG("Freed memory for Full Jacobian preconditioner");
1800 }
1801 enginedata->pfree = NULL;
1802 }
1803
1804 /**
1805 EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1806
1807 'setup' function.
1808 */
1809 static int integrator_ida_psetup_jacobian(realtype tt,
1810 N_Vector yy, N_Vector yp, N_Vector rr,
1811 realtype c_j, void *p_data,
1812 N_Vector tmp1, N_Vector tmp2,
1813 N_Vector tmp3
1814 ){
1815 int i, j, res;
1816 IntegratorSystem *integ;
1817 IntegratorIdaData *enginedata;
1818 IntegratorIdaPrecDataJacobian *precdata;
1819 linsolqr_system_t L;
1820 mtx_matrix_t P;
1821 struct rel_relation **relptr;
1822
1823 integ = (IntegratorSystem *)p_data;
1824 enginedata = integ->enginedata;
1825 precdata = (IntegratorIdaPrecDataJacobian *)(enginedata->precdata);
1826 double *derivatives;
1827 struct var_variable **variables;
1828 int count, status;
1829 char *relname;
1830 mtx_coord_t C;
1831
1832 L = precdata->L;
1833 P = linsolqr_get_matrix(L);
1834 mtx_clear(P);
1835
1836 CONSOLE_DEBUG("Setting up Jacobian preconditioner");
1837
1838 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1839 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1840
1841 /**
1842 @TODO FIXME here we are using the very inefficient and contorted approach
1843 of calculating the whole jacobian, then extracting just the diagonal elements.
1844 */
1845
1846 for(i=0, relptr = enginedata->rellist;
1847 i< enginedata->nrels && relptr != NULL;
1848 ++i, ++relptr
1849 ){
1850 /* get derivatives for this particular relation */
1851 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1852 if(status){
1853 relname = rel_make_name(integ->system, *relptr);
1854 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1855 ASC_FREE(relname);
1856 break;
1857 }
1858 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1859 /* find the diagonal elements */
1860 for(j=0; j<count; ++j){
1861 if(var_deriv(variables[j])){
1862 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), c_j * derivatives[j]);
1863 }else{
1864 mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), derivatives[j]);
1865 }
1866 }
1867 }
1868
1869 mtx_assemble(P);
1870
1871 if(status){
1872 CONSOLE_DEBUG("Error found when evaluating derivatives");
1873 res = 1; goto finish; /* recoverable */
1874 }
1875
1876 integrator_ida_write_incidence(integ);
1877
1878 res = 0;
1879 finish:
1880 ASC_FREE(variables);
1881 ASC_FREE(derivatives);
1882 return res;
1883 };
1884
1885 /**
1886 EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1887
1888 'solve' function.
1889 */
1890 static int integrator_ida_psolve_jacobian(realtype tt,
1891 N_Vector yy, N_Vector yp, N_Vector rr,
1892 N_Vector rvec, N_Vector zvec,
1893 realtype c_j, realtype delta, void *p_data,
1894 N_Vector tmp
1895 ){
1896 IntegratorSystem *integ;
1897 IntegratorIdaData *data;
1898 IntegratorIdaPrecDataJacobian *precdata;
1899 integ = (IntegratorSystem *)p_data;
1900 data = integ->enginedata;
1901 precdata = (IntegratorIdaPrecDataJacobian *)(data->precdata);
1902 linsolqr_system_t L = precdata->L;
1903
1904 linsolqr_add_rhs(L,NV_DATA_S(rvec),FALSE);
1905
1906 mtx_region_t R;
1907 R.row.low = R.col.low = 0;
1908 R.row.high = R.col.high = mtx_order(linsolqr_get_matrix(L)) - 1;
1909 linsolqr_set_region(L,R);
1910
1911 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
1912 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
1913
1914 /// @TODO more here
1915
1916 linsolqr_remove_rhs(L,NV_DATA_S(rvec));
1917
1918 CONSOLE_DEBUG("Solving Jacobian preconditioner (c_j = %f)",c_j);
1919 return 0;
1920 };
1921
1922
1923 /*----------------------------------------------
1924 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1925 */
1926
1927 static void integrator_ida_pcreate_jacobi(IntegratorSystem *integ){
1928 IntegratorIdaData *enginedata =integ->enginedata;
1929 IntegratorIdaPrecDataJacobi *precdata;
1930 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1931
1932 asc_assert(integ->n_y);
1933 precdata->PIii = N_VNew_Serial(integ->n_y);
1934
1935 enginedata->pfree = &integrator_ida_pfree_jacobi;
1936 enginedata->precdata = precdata;
1937 CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1938 }
1939
1940 static void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1941 if(enginedata->precdata){
1942 IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1943 N_VDestroy_Serial(precdata->PIii);
1944
1945 ASC_FREE(precdata);
1946 enginedata->precdata = NULL;
1947 CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1948 }
1949 enginedata->pfree = NULL;
1950 }
1951
1952 /**
1953 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1954
1955 'setup' function.
1956 */
1957 static int integrator_ida_psetup_jacobi(realtype tt,
1958 N_Vector yy, N_Vector yp, N_Vector rr,
1959 realtype c_j, void *p_data,
1960 N_Vector tmp1, N_Vector tmp2,
1961 N_Vector tmp3
1962 ){
1963 int i, j, res;
1964 IntegratorSystem *integ;
1965 IntegratorIdaData *enginedata;
1966 IntegratorIdaPrecDataJacobi *precdata;
1967 struct rel_relation **relptr;
1968
1969 integ = (IntegratorSystem *)p_data;
1970 enginedata = integ->enginedata;
1971 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1972 double *derivatives;
1973 struct var_variable **variables;
1974 int count, status;
1975 char *relname;
1976
1977 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1978
1979 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1980 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1981
1982 /**
1983 @TODO FIXME here we are using the very inefficient and contorted approach
1984 of calculating the whole jacobian, then extracting just the diagonal elements.
1985 */
1986
1987 for(i=0, relptr = enginedata->rellist;
1988 i< enginedata->nrels && relptr != NULL;
1989 ++i, ++relptr
1990 ){
1991
1992 /* get derivatives for this particular relation */
1993 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1994 if(status){
1995 relname = rel_make_name(integ->system, *relptr);
1996 CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1997 ASC_FREE(relname);
1998 break;
1999 }
2000 /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
2001 /* find the diagonal elements */
2002 for(j=0; j<count; ++j){
2003 if(var_sindex(variables[j])==i){
2004 if(var_deriv(variables[j])){
2005 NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
2006 }else{
2007 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
2008 }
2009
2010 }
2011 }
2012 #ifdef PREC_DEBUG
2013 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
2014 #endif
2015 }
2016
2017 if(status){
2018 CONSOLE_DEBUG("Error found when evaluating derivatives");
2019 res = 1; goto finish; /* recoverable */
2020 }
2021
2022 integrator_ida_write_incidence(integ);
2023
2024 res = 0;
2025 finish:
2026 ASC_FREE(variables);
2027 ASC_FREE(derivatives);
2028 return res;
2029 };
2030
2031 /**
2032 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
2033
2034 'solve' function.
2035 */
2036 static int integrator_ida_psolve_jacobi(realtype tt,
2037 N_Vector yy, N_Vector yp, N_Vector rr,
2038 N_Vector rvec, N_Vector zvec,
2039 realtype c_j, realtype delta, void *p_data,
2040 N_Vector tmp
2041 ){
2042 IntegratorSystem *integ;
2043 IntegratorIdaData *data;
2044 IntegratorIdaPrecDataJacobi *precdata;
2045 integ = (IntegratorSystem *)p_data;
2046 data = integ->enginedata;
2047 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
2048
2049 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
2050 N_VProd(precdata->PIii, rvec, zvec);
2051 return 0;
2052 };
2053
2054 /*----------------------------------------------
2055 STATS
2056 */
2057
2058 /**
2059 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
2060 status in a struct instead of separately.
2061
2062 @return IDA_SUCCESS on success.
2063 */
2064 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
2065
2066 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
2067
2068 int res;
2069
2070 /*
2071 There is an error in the documentation for this function in Sundials 2.2.
2072 According the the header file, the hinused stat is not provided.
2073 */
2074 res = IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups,
2075 &s->netfails, &s->qlast, &s->qcur, &s->hlast, &s->hcur,
2076 &s->tcur
2077 );
2078
2079 /* get the missing statistic */
2080 IDAGetActualInitStep(ida_mem, &s->hinused);
2081
2082 return res;
2083 #else
2084
2085 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
2086 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
2087 ,&s->hlast, &s->hcur, &s->tcur
2088 );
2089
2090 #endif
2091 }
2092
2093 /**
2094 This routine just outputs the stats to the CONSOLE_DEBUG routine.
2095
2096 @TODO provide a GUI way of stats reporting from IDA.
2097 */
2098 static void integrator_ida_write_stats(IntegratorIdaStats *stats){
2099 # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
2100 # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
2101 # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
2102 SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
2103 SI(qlast); SI(qcur);
2104 SR(hinused); SR(hlast); SR(hcur); SR(tcur);
2105 # undef SL
2106 # undef SI
2107 # undef SR
2108 }
2109
2110 /*------------------------------------------------------------------------------
2111 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
2112 */
2113
2114 /**
2115 Here we construct the local transfer matrix. It's a bit of a naive
2116 approach; probably far more efficient ways can be worked out. But it will
2117 hopefully be a useful way to examine stability of some spatial
2118 discretisation schemes for PDAE systems.
2119
2120 http://ascendserver.cheme.cmu.edu/wiki/index.php/IDA#Stability
2121 */
2122 static int integrator_ida_transfer_matrix(const IntegratorSystem *integ, struct SystemJacobianStruct *J){
2123 int i=0, res;
2124 enum submat{II_GA=0, II_GD, II_FA, II_FD, II_FDP, II_NUM};
2125
2126 const var_filter_t *matvf[II_NUM] = {
2127 &system_vfilter_algeb
2128 ,&system_vfilter_diff
2129 ,&system_vfilter_algeb
2130 ,&system_vfilter_diff
2131 ,&system_vfilter_deriv
2132 };
2133
2134 const rel_filter_t *matrf[II_NUM] = {
2135 &system_rfilter_algeb
2136 ,&system_rfilter_algeb
2137 ,&system_rfilter_diff
2138 ,&system_rfilter_diff
2139 ,&system_rfilter_diff
2140 };
2141
2142 struct SystemJacobianStruct D[II_NUM];
2143
2144 for(i=0;i<II_NUM;++i){
2145 res = system_jacobian(integ->system, matrf[i], matvf[i], 1/*safe*/ ,&(D[i]));
2146 }
2147
2148 /* compute inverses for matrices that need it */
2149 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
2150 return 1;
2151 }
2152
2153
2154 /**
2155 Our task here is to write the matrices that IDA *should* be seeing. We
2156 are actually making calls to relman_diff in order to do that, so we're
2157 really going back to the variables in the actualy system and computing
2158 row by row what the values are. This should mean just a single call to
2159 each blackbox present in the system (if blackbox caching is working
2160 correctly).
2161 */
2162 static int integrator_ida_write_matrix(const IntegratorSystem *integ, FILE *f, const char *type){
2163 /* IntegratorIdaData *enginedata; */
2164 struct SystemJacobianStruct J = {NULL,NULL,NULL,0,0};
2165 int status=1;
2166 mtx_region_t R;
2167
2168 if(type==NULL)type = "dx'/dx";
2169
2170 if(0==strcmp(type,"dg/dz")){
2171 CONSOLE_DEBUG("Calculating dg/dz...");
2172 status = system_jacobian(integ->system
2173 , &system_rfilter_algeb, &system_vfilter_algeb
2174 , 1 /* safe */
2175 , &J
2176 );
2177 }else if(0==strcmp(type,"dg/dx")){
2178 CONSOLE_DEBUG("Calculating dg/dx...");
2179 status = system_jacobian(integ->system
2180 , &system_rfilter_algeb, &system_vfilter_diff
2181 , 1 /* safe */
2182 , &J
2183 );
2184 }else if(0==strcmp(type,"df/dx'")){
2185 CONSOLE_DEBUG("Calculating df/dx'...");
2186 status = system_jacobian(integ->system
2187 , &system_rfilter_diff, &system_vfilter_deriv
2188 , 1 /* safe */
2189 , &J
2190 );
2191 }else if(0==strcmp(type,"df/dz")){
2192 CONSOLE_DEBUG("Calculating df/dz...");
2193 status = system_jacobian(integ->system
2194 , &system_rfilter_diff, &system_vfilter_algeb
2195 , 1 /* safe */
2196 , &J
2197 );
2198 }else if(0==strcmp(type,"df/dx")){
2199 CONSOLE_DEBUG("Calculating df/dx...");
2200 status = system_jacobian(integ->system
2201 , &system_rfilter_diff, &system_vfilter_diff
2202 , 1 /* safe */
2203 , &J
2204 );
2205 }else if(0==strcmp(type,"dF/dy")){
2206 CONSOLE_DEBUG("Calculating dF/dy...");
2207 status = system_jacobian(integ->system
2208 , &system_rfilter_all, &system_vfilter_nonderiv
2209 , 1 /* safe */
2210 , &J
2211 );
2212 }else if(0==strcmp(type,"dF/dy'")){
2213 CONSOLE_DEBUG("Calculating dF/dy'...");
2214 status = system_jacobian(integ->system
2215 , &system_rfilter_all, &system_vfilter_deriv
2216 , 1 /* safe */
2217 , &J
2218 );
2219 }else if(0==strcmp(type,"dx'/dx")){
2220 /* system state transfer matrix dyd'/dyd */
2221 status = integrator_ida_transfer_matrix(integ, &J);
2222 }else{
2223 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid matrix type '%s'",type);
2224 return 1;
2225 }
2226
2227 if(status){
2228 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating matrix");
2229 }else{
2230 /* send the region explicitly, so that we handle non-square correctly */
2231 R.row.low = 0; R.col.low = 0;
2232 R.row.high = J.n_rels - 1; R.col.high = J.n_vars - 1;
2233 /* note that we're not fussy about empty matrices here... */
2234 mtx_write_region_mmio(f,J.M,&R);
2235 }
2236
2237 if(J.vars)ASC_FREE(J.vars);
2238 if(J.rels)ASC_FREE(J.rels);
2239 if(J.M)mtx_destroy(J.M);
2240
2241 return status;
2242 }
2243
2244 /**
2245 This routine outputs matrix structure in a crude text format, for the sake
2246 of debugging.
2247 */
2248 static void integrator_ida_write_incidence(IntegratorSystem *integ){
2249 int i, j;
2250 struct rel_relation **relptr;
2251 IntegratorIdaData *enginedata = integ->enginedata;
2252 double *derivatives;
2253 struct var_variable **variables;
2254 int count, status;
2255 char *relname;
2256
2257 if(enginedata->nrels > 100){
2258 CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
2259 return;
2260 }
2261
2262 variables = ASC_NEW_ARRAY(struct var_variable *, integ->n_y * 2);
2263 derivatives = ASC_NEW_ARRAY(double, integ->n_y * 2);
2264
2265 CONSOLE_DEBUG("Outputting incidence information to console...");
2266
2267 for(i=0, relptr = enginedata->rellist;
2268 i< enginedata->nrels && relptr != NULL;
2269 ++i, ++relptr
2270 ){
2271 relname = rel_make_name(integ->system, *relptr);
2272
2273 /* get derivatives for this particular relation */
2274 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
2275 if(status){
2276 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
2277 ASC_FREE(relname);
2278 break;
2279 }
2280
2281 fprintf(stderr,"%3d:%-15s:",i,relname);
2282 ASC_FREE(relname);
2283
2284 for(j=0; j<count; ++j){
2285 if(var_deriv(variables[j])){
2286 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(integ,variables[j]));
2287 }else{
2288 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
2289 }
2290 }
2291 fprintf(stderr,"\n");
2292 }
2293 ASC_FREE(variables);
2294 ASC_FREE(derivatives);
2295 }
2296
2297 /* @return 0 on success */
2298 int integrator_ida_debug(const IntegratorSystem *integ, FILE *fp){
2299 char *varname, *relname;
2300 struct var_variable **vlist, *var;
2301 struct rel_relation **rlist, *rel;
2302 long vlen, rlen;
2303 long i;
2304 long di;
2305
2306 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",integ->n_y);
2307
2308 /* if(integrator_sort_obs_vars(integ))return 10; */
2309
2310 if(integ->y && integ->ydot){
2311 fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
2312 fprintf(fp,"index\t%-15s\tydot\n","y");
2313 fprintf(fp,"-----\t%-15s\t-----\n","-----");
2314 for(i=0;i<integ->n_y;++i){
2315 varname = var_make_name(integ->system, integ->y[i]);
2316 fprintf(fp,"%ld\t%-15s\t",i,varname);
2317 if(integ->ydot[i]){
2318 ASC_FREE(varname);
2319 varname = var_make_name(integ->system, integ->ydot[i]);
2320 fprintf(fp,"%s\n",varname);
2321 ASC_FREE(varname);
2322 }else{
2323 fprintf(fp,".\n");
2324 ASC_FREE(varname);
2325 }
2326 }
2327 }else{
2328 fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
2329 }
2330
2331 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
2332 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
2333 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
2334
2335
2336 /* visit all the slv_system_t master var lists to collect vars */
2337 /* find the vars mostly in this one */
2338 vlist = slv_get_solvers_var_list(integ->system);
2339 vlen = slv_get_num_solvers_vars(integ->system);
2340 for(i=0;i<vlen;i++){
2341 var = vlist[i];
2342
2343 varname = var_make_name(integ->system, var);
2344 fprintf(fp,"%ld\t%-15s\t",i,varname);
2345
2346 if(var_fixed(var)){
2347 // it's fixed, so not really a DAE var
2348 fprintf(fp,"(fixed)\n");
2349 }else if(!var_active(var)){
2350 // inactive
2351 fprintf(fp,"(inactive)\n");
2352 }else if(!var_incident(var)){
2353 // not incident
2354 fprintf(fp,"(not incident)\n");
2355 }else{
2356 if(var_deriv(var)){
2357 if(integ->y_id){
2358 di = integrator_ida_diffindex1(integ,var);
2359 if(di>=0){
2360 ASC_FREE(varname);
2361 varname = var_make_name(integ->system,vlist[di]);
2362 fprintf(fp,".\tdiff(%ld='%s')\n",di,varname);
2363 }else{
2364 fprintf(fp,".\tdiff(???,err=%ld)\n",di);
2365 }
2366 }else{
2367 fprintf(fp,".\tderiv... of??\n");
2368 }
2369 }else{
2370 fprintf(fp,"%d\t.\n",var_sindex(var));
2371 }
2372 }
2373 ASC_FREE(varname);
2374 }
2375
2376 /* let's write out the relations too */
2377 rlist = slv_get_solvers_rel_list(integ->system);
2378 rlen = slv_get_num_solvers_rels(integ->system);
2379
2380 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
2381 fprintf(fp,"index\tname\n");
2382 fprintf(fp,"-----\t----\n");
2383 for(i=0; i<rlen; ++i){
2384 rel = rlist[i];
2385 relname = rel_make_name(integ->system,rel);
2386 fprintf(fp,"%ld\t%s\n",i,relname);
2387 ASC_FREE(relname);
2388 }
2389
2390 /* write out the derivative chains */
2391 fprintf(fp,"\nDERIVATIVE CHAINS\n");
2392 if(integrator_ida_analyse_debug(integ,stderr)){
2393 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
2394 return 340;
2395 }
2396 fprintf(fp,"\n");
2397
2398 /* and lets write block debug output */
2399 system_block_debug(integ->system, fp);
2400
2401 return 0; /* success */
2402 }
2403
2404 /*----------------------------------------------
2405 ERROR REPORTING
2406 */
2407 /**
2408 Error message reporter function to be passed to IDA. All error messages
2409 will trigger a call to this function, so we should find everything
2410 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
2411 panel (in the case of PyGTK).
2412 */
2413 static void integrator_ida_error(int error_code
2414 , const char *module, const char *function
2415 , char *msg, void *eh_data
2416 ){
2417 IntegratorSystem *integ;
2418 error_severity_t sev;
2419
2420 /* cast back the IntegratorSystem, just in case we need it */
2421 integ = (IntegratorSystem *)eh_data;
2422
2423 /* severity depends on the sign of the error_code value */
2424 if(error_code <= 0){
2425 sev = ASC_PROG_ERR;
2426 }else{
2427 sev = ASC_PROG_WARNING;
2428 }
2429
2430 /* use our all-purpose error reporting to get stuff back to the GUI */
2431 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
2432 }

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