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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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