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

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