/[ascend]/trunk/base/generic/integrator/ida.c
ViewVC logotype

Contents of /trunk/base/generic/integrator/ida.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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