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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1395 - (hide annotations) (download) (as text)
Sat Apr 21 15:11:45 2007 UTC (15 years, 4 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 johnpye 669 /* ASCEND modelling environment
2 johnpye 1129 Copyright (C) 2006-2007 Carnegie Mellon University
3 johnpye 669
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 johnpye 993
23 johnpye 1068 IDA provides the non-linear parts, as well as a number of pluggable linear
24     solvers: dense, banded and krylov types.
25 johnpye 993
26 johnpye 1068 We also implement here an EXPERIMENTAL direct sparse linear solver for IDA
27     using the ASCEND linsolqr routines.
28    
29 johnpye 669 @see http://www.llnl.gov/casc/sundials/
30     *//*
31     by John Pye, May 2006
32     */
33    
34 johnpye 1232 #define _GNU_SOURCE
35    
36 johnpye 1183 #include "ida.h"
37    
38 johnpye 1232 #include <signal.h>
39     #include <setjmp.h>
40     #include <fenv.h>
41 johnpye 1237 #include <math.h>
42 johnpye 1232
43 johnpye 669 /* SUNDIALS includes */
44     #ifdef ASC_WITH_IDA
45 johnpye 913 # include <sundials/sundials_config.h>
46 johnpye 945 # include <sundials/sundials_dense.h>
47 johnpye 782 # include <ida/ida.h>
48     # include <nvector/nvector_serial.h>
49     # include <ida/ida_spgmr.h>
50 johnpye 969 # include <ida/ida_spbcgs.h>
51     # include <ida/ida_sptfqmr.h>
52 johnpye 945 # include <ida/ida_dense.h>
53 johnpye 669 # ifndef IDA_SUCCESS
54     # error "Failed to include SUNDIALS IDA header file"
55     # endif
56 johnpye 1181 #else
57 johnpye 1213 # error "If you're building this file, you should have ASC_WITH_IDA"
58 johnpye 669 #endif
59    
60 johnpye 994 #ifdef ASC_WITH_MMIO
61     # include <mmio.h>
62     #endif
63    
64 johnpye 1183 #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 johnpye 1181
70 johnpye 1316 #include <system/slv_client.h>
71     #include <system/relman.h>
72     #include <system/block.h>
73     #include <system/slv_stdcalls.h>
74 johnpye 1317 #include <system/jacobian.h>
75 jpye 1337 #include <system/bndman.h>
76 johnpye 1181
77 johnpye 1183 #include "idalinear.h"
78 johnpye 1239 #include "idaanalyse.h"
79 johnpye 1309 #include "ida_impl.h"
80 johnpye 1183
81 johnpye 913 /*
82 johnpye 969 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
83 johnpye 913 */
84     #ifndef SUNDIALS_VERSION_MINOR
85 johnpye 921 # ifdef __GNUC__
86     # warning "GUESSING SUNDIALS VERSION 2.2"
87     # endif
88 johnpye 913 # define SUNDIALS_VERSION_MINOR 2
89     #endif
90     #ifndef SUNDIALS_VERSION_MAJOR
91     # define SUNDIALS_VERSION_MAJOR 2
92     #endif
93    
94 johnpye 669 /* 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 johnpye 1246 /* #define FEX_DEBUG */
100 johnpye 997 #define JEX_DEBUG
101 johnpye 1246 /* #define DJEX_DEBUG */
102 jpye 1339 /* #define SOLVE_DEBUG */
103 johnpye 993 #define STATS_DEBUG
104     #define PREC_DEBUG
105 jpye 1339 /* #define ROOT_DEBUG */
106    
107 johnpye 1237 /* #define DIFFINDEX_DEBUG */
108 johnpye 1279 /* #define ANALYSE_DEBUG */
109 johnpye 1303 /* #define DESTROY_DEBUG */
110 johnpye 1309 /* #define MATRIX_DEBUG */
111 johnpye 1303
112 johnpye 993 /**
113     Everthing that the outside world needs to know about IDA
114     */
115 johnpye 977 const IntegratorInternals integrator_ida_internals = {
116     integrator_ida_create
117     ,integrator_ida_params_default
118 johnpye 1068 ,integrator_ida_analyse
119 johnpye 977 ,integrator_ida_solve
120 johnpye 1306 ,integrator_ida_write_matrix
121 johnpye 1196 ,integrator_ida_debug
122 johnpye 977 ,integrator_ida_free
123     ,INTEG_IDA
124     ,"IDA"
125     };
126    
127 johnpye 993 /*-------------------------------------------------------------
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 johnpye 1306 typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *sys);
136 johnpye 993 typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);
137    
138 jpye 1337
139 johnpye 669 /**
140 johnpye 980 Struct containing any stuff that IDA needs that doesn't fit into the
141 johnpye 669 common IntegratorSystem struct.
142     */
143 johnpye 993 typedef struct IntegratorIdaDataStruct{
144 jpye 1337
145     struct rel_relation **rellist; /**< NULL terminated list of ACTIVE rels */
146     int nrels; /* number of ACTIVE rels */
147    
148 johnpye 972 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
149 jpye 1337 int nbnds; /* number of boundaries */
150    
151 johnpye 669 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
152 johnpye 1196 var_filter_t vfilter;
153 jpye 1337 rel_filter_t rfilter; /**< Used to filter relations from solver's rellist (@TODO needs work) */
154 johnpye 993 void *precdata; /**< For use by the preconditioner */
155     IntegratorIdaPrecFreeFn *pfree; /**< Store instructions here on how to free precdata */
156 jpye 1337
157 johnpye 669 } IntegratorIdaData;
158    
159 jpye 1337
160 johnpye 1060 typedef struct IntegratorIdaPrecDJStruct{
161 johnpye 993 N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */
162     } IntegratorIdaPrecDataJacobi;
163    
164 jpye 1395 typedef struct IntegratorIdaPrecDJFStruct{
165     linsolqr_system_t L;
166     } IntegratorIdaPrecDataJacobian;
167    
168 johnpye 993 /**
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 johnpye 1263
173 johnpye 993 @note IDA uses a different convention for function pointer types, so no '*'.
174 johnpye 669 */
175 johnpye 1060 typedef struct IntegratorIdaPrecStruct{
176 johnpye 993 IntegratorIdaPrecCreateFn *pcreate;
177     IDASpilsPrecSetupFn psetup;
178     IDASpilsPrecSolveFn psolve;
179     } IntegratorIdaPrec;
180    
181 johnpye 669 /* residual function forward declaration */
182 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
183 johnpye 669
184 johnpye 909 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 johnpye 669 /* 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 johnpye 1012 /* dense jacobian evaluation for IDADense dense direct linear solver */
196 johnpye 945 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 johnpye 1012 /* sparse jacobian evaluation for ASCEND's sparse direct solver */
203     IntegratorSparseJacFn integrator_ida_sjex;
204    
205 jpye 1337 /* boundary-detection function */
206     int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
207    
208 johnpye 1060 typedef struct IntegratorIdaStatsStruct{
209 johnpye 993 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 johnpye 1239 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
219    
220 johnpye 1240 /*static IntegratorVarVisitorFn integrator_dae_classify_var;
221 johnpye 1239 static void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
222 johnpye 1240 static void integrator_dae_show_var(IntegratorSystem *sys, struct var_variable *var, const int *varindx); */
223 johnpye 1239
224 johnpye 993 int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
225     void integrator_ida_write_stats(IntegratorIdaStats *stats);
226 johnpye 1306 void integrator_ida_write_incidence(IntegratorSystem *sys);
227 johnpye 1217
228 johnpye 993 /*------
229 jpye 1395 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 johnpye 1263 Jacobi preconditioner -- experimental
258 johnpye 993 */
259    
260     int integrator_ida_psetup_jacobi(realtype tt,
261 johnpye 992 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 johnpye 993 int integrator_ida_psolve_jacobi(realtype tt,
268 johnpye 992 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 johnpye 1306 void integrator_ida_pcreate_jacobi(IntegratorSystem *sys);
275 johnpye 992
276 johnpye 993 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 johnpye 669 /*-------------------------------------------------------------
285     SETUP/TEARDOWN ROUTINES
286     */
287 johnpye 1306 void integrator_ida_create(IntegratorSystem *sys){
288 johnpye 669 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
289     IntegratorIdaData *enginedata;
290     enginedata = ASC_NEW(IntegratorIdaData);
291 johnpye 1221 CONSOLE_DEBUG("enginedata = %p",enginedata);
292 johnpye 669 enginedata->rellist = NULL;
293 johnpye 953 enginedata->safeeval = 0;
294 johnpye 1213 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED;
295 johnpye 1221 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;
296 johnpye 1007 enginedata->pfree = NULL;
297 johnpye 993
298 johnpye 994 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
299     enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
300    
301 johnpye 1306 sys->enginedata = (void *)enginedata;
302 johnpye 993
303 johnpye 1306 integrator_ida_params_default(sys);
304 ben.allan 704 }
305 johnpye 669
306     void integrator_ida_free(void *enginedata){
307 johnpye 1303 #ifdef DESTROY_DEBUG
308 johnpye 1221 CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);
309 johnpye 1303 #endif
310 johnpye 669 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
311 johnpye 1221 asc_assert(d);
312 johnpye 993 if(d->pfree){
313 johnpye 1221 CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);
314 johnpye 993 /* free the preconditioner data, whatever it happens to be */
315     (d->pfree)(enginedata);
316     }
317 jpye 1337
318     ASC_FREE(d->rellist);
319    
320 johnpye 1303 #ifdef DESTROY_DEBUG
321 johnpye 1221 CONSOLE_DEBUG("Now destroying the enginedata");
322 johnpye 1303 #endif
323 johnpye 669 ASC_FREE(d);
324 johnpye 1303 #ifdef DESTROY_DEBUG
325 johnpye 1221 CONSOLE_DEBUG("enginedata freed");
326 johnpye 1303 #endif
327 johnpye 669 }
328    
329 johnpye 1306 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *sys){
330 johnpye 669 IntegratorIdaData *d;
331 johnpye 1306 assert(sys!=NULL);
332     assert(sys->enginedata!=NULL);
333     assert(sys->engine==INTEG_IDA);
334     d = ((IntegratorIdaData *)(sys->enginedata));
335 johnpye 669 return d;
336     }
337    
338     /*-------------------------------------------------------------
339 johnpye 942 PARAMETERS FOR IDA
340     */
341    
342     enum ida_parameters{
343 johnpye 945 IDA_PARAM_LINSOLVER
344 johnpye 970 ,IDA_PARAM_MAXL
345 johnpye 1279 ,IDA_PARAM_MAXORD
346 johnpye 945 ,IDA_PARAM_AUTODIFF
347 johnpye 972 ,IDA_PARAM_CALCIC
348 johnpye 953 ,IDA_PARAM_SAFEEVAL
349 johnpye 944 ,IDA_PARAM_RTOL
350     ,IDA_PARAM_ATOL
351     ,IDA_PARAM_ATOLVECT
352 johnpye 951 ,IDA_PARAM_GSMODIFIED
353 johnpye 991 ,IDA_PARAM_MAXNCF
354 johnpye 992 ,IDA_PARAM_PREC
355 johnpye 951 ,IDA_PARAMS_SIZE
356 johnpye 942 };
357    
358     /**
359     Here the full set of parameters is defined, along with upper/lower bounds,
360 johnpye 1306 etc. The values are stuck into the sys->params structure.
361 johnpye 942
362 johnpye 992 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 johnpye 942 @return 0 on success
367     */
368 johnpye 1306 int integrator_ida_params_default(IntegratorSystem *sys){
369     asc_assert(sys!=NULL);
370     asc_assert(sys->engine==INTEG_IDA);
371 johnpye 942 slv_parameters_t *p;
372 johnpye 1306 p = &(sys->params);
373 johnpye 942
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 jpye 1337 ,(SlvParameterInitBool){{"autodiff"
390 johnpye 942 ,"Use auto-diff?",1
391     ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
392     }, TRUE}
393     );
394    
395 johnpye 993 slv_param_char(p,IDA_PARAM_CALCIC
396 jpye 1337 ,(SlvParameterInitChar){{"calcic"
397 johnpye 993 ,"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 johnpye 1133 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
404 johnpye 972 );
405    
406 johnpye 953 slv_param_bool(p,IDA_PARAM_SAFEEVAL
407 jpye 1337 ,(SlvParameterInitBool){{"safeeval"
408 johnpye 953 ,"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 johnpye 944 slv_param_bool(p,IDA_PARAM_ATOLVECT
416 jpye 1337 ,(SlvParameterInitBool){{"atolvect"
417 johnpye 944 ,"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 jpye 1337 ,(SlvParameterInitReal){{"atol"
426 johnpye 944 ,"Scalar absolute error tolerance",1
427     ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
428 johnpye 997 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
429 johnpye 1251 }, 1e-5, 0, 1e10 }
430 johnpye 944 );
431    
432     slv_param_real(p,IDA_PARAM_RTOL
433 jpye 1337 ,(SlvParameterInitReal){{"rtol"
434 johnpye 944 ,"Scalar relative error tolerance",1
435 johnpye 963 ,"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 johnpye 997 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
439 johnpye 1141 }, 1e-4, 0, 1 }
440 johnpye 944 );
441    
442 johnpye 945 slv_param_char(p,IDA_PARAM_LINSOLVER
443 jpye 1337 ,(SlvParameterInitChar){{"linsolver"
444 johnpye 945 ,"Linear solver",1
445 johnpye 1012 ,"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 johnpye 1300 }, "DENSE"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
451 johnpye 945 );
452    
453 johnpye 970 slv_param_int(p,IDA_PARAM_MAXL
454 jpye 1337 ,(SlvParameterInitInt){{"maxl"
455 johnpye 970 ,"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 johnpye 1279 slv_param_int(p,IDA_PARAM_MAXORD
464 jpye 1337 ,(SlvParameterInitInt){{"maxord"
465 johnpye 1279 ,"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 johnpye 951 slv_param_bool(p,IDA_PARAM_GSMODIFIED
472 jpye 1337 ,(SlvParameterInitBool){{"gsmodified"
473 johnpye 991 ,"Gram-Schmidt Orthogonalisation Scheme", 2
474 johnpye 969 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
475     " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
476 johnpye 951 }, TRUE}
477     );
478    
479 johnpye 991 slv_param_int(p,IDA_PARAM_MAXNCF
480 jpye 1337 ,(SlvParameterInitInt){{"maxncf"
481 johnpye 991 ,"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 johnpye 992 slv_param_char(p,IDA_PARAM_PREC
488 jpye 1337 ,(SlvParameterInitChar){{"prec"
489 johnpye 992 ,"Preconditioner",1
490     ,"See IDA manual, section section 5.6.8."
491     },"NONE"}, (char *[]){"NONE","DIAG",NULL}
492     );
493    
494 johnpye 942 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 johnpye 669 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
503     */
504    
505 johnpye 1060 /*static double div1(double a, double b){
506 johnpye 997 return a/b;
507 johnpye 1060 }*/
508 johnpye 997
509 johnpye 985 typedef int IdaFlagFn(void *,int *);
510 johnpye 991 typedef char *IdaFlagNameFn(int);
511 johnpye 985
512 johnpye 1049 /* return 0 on success */
513 johnpye 669 int integrator_ida_solve(
514 johnpye 1306 IntegratorSystem *sys
515 johnpye 669 , unsigned long start_index
516     , unsigned long finish_index
517     ){
518     void *ida_mem;
519 jpye 1337 int flag, flag1, t_index;
520 johnpye 972 realtype t0, reltol, abstol, t, tret, tout1;
521 johnpye 993 N_Vector y0, yp0, abstolvect, ypret, yret, id;
522 johnpye 669 IntegratorIdaData *enginedata;
523 johnpye 945 char *linsolver;
524 johnpye 971 int maxl;
525 johnpye 985 IdaFlagFn *flagfn;
526     IdaFlagNameFn *flagnamefn;
527 johnpye 986 const char *flagfntype;
528 johnpye 993 char *pname = NULL;
529 jpye 1337 struct rel_relation **rels;
530     int *rootsfound;
531    
532 johnpye 1306 #ifdef SOLVE_DEBUG
533 jpye 1337 char *varname, *relname;
534 johnpye 1306 #endif
535 jpye 1337
536     int i,j,n_activerels,n_solverrels;
537 johnpye 993 const IntegratorIdaPrec *prec = NULL;
538     int icopt; /* initial conditions strategy */
539 johnpye 669
540     CONSOLE_DEBUG("STARTING IDA...");
541    
542 johnpye 1306 enginedata = integrator_ida_enginedata(sys);
543 johnpye 953
544 johnpye 1306 enginedata->safeeval = SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_SAFEEVAL);
545 johnpye 669 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
546    
547     /* store reference to list of relations (in enginedata) */
548 jpye 1337 n_solverrels = slv_get_num_solvers_rels(sys->system);
549    
550     n_activerels = slv_count_solvers_rels(sys->system, &integrator_ida_rel);
551    
552 johnpye 1306 enginedata->bndlist = slv_get_solvers_bnd_list(sys->system);
553 jpye 1337 enginedata->nbnds = slv_get_num_solvers_bnds(sys->system);
554 johnpye 972
555 jpye 1337 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 johnpye 1306 CONSOLE_DEBUG("Number of dependent vars: %d",sys->n_y);
578 jpye 1337 CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);
579 johnpye 669
580 jpye 1337 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 johnpye 1049 return 1; /* failure */
588 johnpye 669 }
589    
590 johnpye 1279 #ifdef SOLVE_DEBUG
591 johnpye 1306 integrator_ida_debug(sys,stderr);
592 johnpye 1279 #endif
593 johnpye 1240
594 johnpye 669 /* retrieve initial values from the system */
595 johnpye 980
596 johnpye 944 /** @TODO fix this, the starting time != first sample */
597 johnpye 1306 t0 = integrator_get_t(sys);
598 johnpye 944 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
599 johnpye 669
600 johnpye 928 CONSOLE_DEBUG("RETRIEVING y0");
601    
602 jpye 1337 y0 = N_VNew_Serial(sys->n_y);
603 johnpye 1306 integrator_get_y(sys,NV_DATA_S(y0));
604 johnpye 669
605 johnpye 991 #ifdef SOLVE_DEBUG
606 johnpye 928 CONSOLE_DEBUG("RETRIEVING yp0");
607 johnpye 991 #endif
608 johnpye 928
609 jpye 1337 yp0 = N_VNew_Serial(sys->n_y);
610 johnpye 1306 integrator_get_ydot(sys,NV_DATA_S(yp0));
611 johnpye 669
612 johnpye 991 #ifdef SOLVE_DEBUG
613 johnpye 669 N_VPrint_Serial(yp0);
614     CONSOLE_DEBUG("yp0 is at %p",&yp0);
615 johnpye 991 #endif
616 johnpye 669
617     /* create IDA object */
618     ida_mem = IDACreate();
619    
620 johnpye 980 /* relative error tolerance */
621 johnpye 1306 reltol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_RTOL);
622 johnpye 949 CONSOLE_DEBUG("rtol = %8.2e",reltol);
623 johnpye 669
624     /* allocate internal memory */
625 johnpye 1306 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_ATOLVECT)){
626 johnpye 944 /* vector of absolute tolerances */
627     CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
628 jpye 1337 abstolvect = N_VNew_Serial(sys->n_y);
629 johnpye 1306 integrator_get_atol(sys,NV_DATA_S(abstolvect));
630 johnpye 944
631     flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
632 johnpye 949
633     N_VDestroy_Serial(abstolvect);
634 johnpye 944 }else{
635     /* scalar absolute tolerance (one value for all) */
636 johnpye 1306 abstol = SLV_PARAM_REAL(&(sys->params),IDA_PARAM_ATOL);
637 johnpye 944 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 johnpye 669 if(flag==IDA_MEM_NULL){
642     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
643 johnpye 1049 return 2;
644 johnpye 669 }else if(flag==IDA_MEM_FAIL){
645     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
646 johnpye 1049 return 3;
647 johnpye 669 }else if(flag==IDA_ILL_INPUT){
648     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
649 johnpye 1049 return 4;
650 johnpye 669 }/* else success */
651    
652     /* set optional inputs... */
653 johnpye 1306 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 johnpye 951 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
660 johnpye 950 }
661 johnpye 991
662 johnpye 1306 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 johnpye 991
665 johnpye 1306 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 johnpye 1279
668 johnpye 669 /* there's no capability for setting *minimum* step size in IDA */
669    
670 johnpye 912
671 johnpye 669 /* attach linear solver module, using the default value of maxl */
672 johnpye 1306 linsolver = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_LINSOLVER);
673 johnpye 946 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
674 johnpye 1012 if(strcmp(linsolver,"ASCEND")==0){
675 jpye 1337 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",sys->n_y);
676     IDAASCEND(ida_mem,sys->n_y);
677 johnpye 1306 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)sys);
678 johnpye 1012
679     flagfntype = "IDAASCEND";
680     flagfn = &IDAASCENDGetLastFlag;
681     flagnamefn = &IDAASCENDGetReturnFlagName;
682    
683     }else if(strcmp(linsolver,"DENSE")==0){
684 jpye 1337 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",sys->n_y);
685     flag = IDADense(ida_mem, sys->n_y);
686 johnpye 969 switch(flag){
687     case IDADENSE_SUCCESS: break;
688 johnpye 1049 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 johnpye 969 }
693 johnpye 980
694 johnpye 1306 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
695 johnpye 969 CONSOLE_DEBUG("USING AUTODIFF");
696 johnpye 1306 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)sys);
697 johnpye 969 switch(flag){
698     case IDADENSE_SUCCESS: break;
699 johnpye 1049 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
700 johnpye 969 }
701     }else{
702     CONSOLE_DEBUG("USING NUMERICAL DIFF");
703     }
704 johnpye 985
705 johnpye 986 flagfntype = "IDADENSE";
706 johnpye 985 flagfn = &IDADenseGetLastFlag;
707     flagnamefn = &IDADenseGetReturnFlagName;
708 johnpye 969 }else{
709     /* remaining methods are all SPILS */
710     CONSOLE_DEBUG("IDA SPILS");
711    
712 johnpye 1306 maxl = SLV_PARAM_INT(&(sys->params),IDA_PARAM_MAXL);
713 johnpye 971 CONSOLE_DEBUG("maxl = %d",maxl);
714    
715 johnpye 993 /* what preconditioner? */
716 johnpye 1306 pname = SLV_PARAM_CHAR(&(sys->params),IDA_PARAM_PREC);
717 johnpye 992 if(strcmp(pname,"NONE")==0){
718 johnpye 993 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 johnpye 1049 return 7;
724 johnpye 992 }
725    
726 johnpye 993 /* which SPILS linear solver? */
727 johnpye 969 if(strcmp(linsolver,"SPGMR")==0){
728     CONSOLE_DEBUG("IDA SPGMR");
729 johnpye 971 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
730 johnpye 969 }else if(strcmp(linsolver,"SPBCG")==0){
731     CONSOLE_DEBUG("IDA SPBCG");
732 johnpye 971 flag = IDASpbcg(ida_mem, maxl);
733 johnpye 969 }else if(strcmp(linsolver,"SPTFQMR")==0){
734     CONSOLE_DEBUG("IDA SPTFQMR");
735 johnpye 971 flag = IDASptfqmr(ida_mem,maxl);
736 johnpye 969 }else{
737     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
738 johnpye 1049 return 8;
739 johnpye 980 }
740    
741 johnpye 993 if(prec){
742     /* assign the preconditioner to the linear solver */
743 johnpye 1306 (prec->pcreate)(sys);
744     IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)sys);
745 johnpye 993 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
746     }else{
747     CONSOLE_DEBUG("No preconditioner");
748     }
749 johnpye 992
750 johnpye 986 flagfntype = "IDASPILS";
751 johnpye 985 flagfn = &IDASpilsGetLastFlag;
752     flagnamefn = &IDASpilsGetReturnFlagName;
753    
754 johnpye 942 if(flag==IDASPILS_MEM_NULL){
755     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
756 johnpye 1049 return 9;
757 johnpye 945 }else if(flag==IDASPILS_MEM_FAIL){
758     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
759 johnpye 1049 return 9;
760 johnpye 942 }/* else success */
761 johnpye 945
762     /* assign the J*v function */
763 johnpye 1306 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_AUTODIFF)){
764 johnpye 945 CONSOLE_DEBUG("USING AUTODIFF");
765 johnpye 1306 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)sys);
766 johnpye 945 if(flag==IDASPILS_MEM_NULL){
767     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
768 johnpye 1049 return 10;
769 johnpye 945 }else if(flag==IDASPILS_LMEM_NULL){
770     ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
771 johnpye 1049 return 10;
772 johnpye 945 }/* else success */
773     }else{
774     CONSOLE_DEBUG("USING NUMERICAL DIFF");
775     }
776 johnpye 951
777 johnpye 969 if(strcmp(linsolver,"SPGMR")==0){
778     /* select Gram-Schmidt orthogonalisation */
779 johnpye 1306 if(SLV_PARAM_BOOL(&(sys->params),IDA_PARAM_GSMODIFIED)){
780 johnpye 969 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 johnpye 1049 return 11;
785 johnpye 969 }
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 johnpye 1049 return 11;
792 johnpye 969 }
793 johnpye 951 }
794     }
795 johnpye 945 }
796 johnpye 669
797 johnpye 909 /* set linear solver optional inputs...
798 johnpye 980 ...nothing here at the moment...
799 johnpye 909 */
800    
801 johnpye 972 /* calculate initial conditions */
802 johnpye 993 icopt = 0;
803 johnpye 1306 if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"Y")==0){
804 johnpye 993 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
805     icopt = IDA_Y_INIT;
806     asc_assert(icopt!=0);
807 johnpye 1306 }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
808 johnpye 993 CONSOLE_DEBUG("Solving initial conditions using values of yd");
809     icopt = IDA_YA_YDP_INIT;
810     asc_assert(icopt!=0);
811 johnpye 1306 id = N_VNew_Serial(sys->n_y);
812     for(i=0; i < sys->n_y; ++i){
813     if(sys->ydot[i] == NULL){
814 johnpye 993 NV_Ith_S(id,i) = 0.0;
815 johnpye 1279 #ifdef SOLVE_DEBUG
816 johnpye 1306 varname = var_make_name(sys->system,sys->y[i]);
817 johnpye 993 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
818     ASC_FREE(varname);
819 johnpye 1279 #endif
820 johnpye 993 }else{
821 johnpye 1279 #ifdef SOLVE_DEBUG
822 johnpye 993 CONSOLE_DEBUG("y[%d] is differential",i);
823 johnpye 1279 #endif
824 johnpye 993 NV_Ith_S(id,i) = 1.0;
825     }
826     }
827     IDASetId(ida_mem, id);
828     N_VDestroy_Serial(id);
829 johnpye 1306 }else if(strcmp(SLV_PARAM_CHAR(&sys->params,IDA_PARAM_CALCIC),"NONE")==0){
830 johnpye 1228 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
831 johnpye 993 }else{
832 johnpye 1228 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
833 johnpye 993 }
834 johnpye 913
835 johnpye 993 if(icopt){
836 johnpye 1306 sys->currentstep=0;
837 johnpye 1016 t_index=start_index + 1;
838 johnpye 1306 tout1 = samplelist_get(sys->samples, t_index);
839 johnpye 944
840 johnpye 972 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
841 johnpye 913
842 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
843 johnpye 996 /* catch SIGFPE if desired to */
844     if(enginedata->safeeval){
845 johnpye 997 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
846 johnpye 1246 Asc_SignalHandlerPush(SIGFPE,SIG_DFL);
847 johnpye 996 }else{
848 johnpye 1228 # ifdef FEX_DEBUG
849 johnpye 997 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
850 johnpye 1228 # endif
851 johnpye 996 Asc_SignalHandlerPushDefault(SIGFPE);
852     }
853 johnpye 1009 if (setjmp(g_fpe_env)==0) {
854 johnpye 1142 #endif
855 johnpye 1263
856 johnpye 993 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
857 johnpye 1254 flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3 */
858 johnpye 993 # else
859 johnpye 1254 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
860 johnpye 993 # endif
861 johnpye 1254 /* check flags and output status */
862     switch(flag){
863     case IDA_SUCCESS:
864     CONSOLE_DEBUG("Initial conditions solved OK");
865     break;
866 johnpye 669
867 johnpye 1254 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 johnpye 1049 return 12;
879 johnpye 1254 }
880     ERROR_REPORTER_HERE(ASC_PROG_ERR
881     ,"%s returned flag '%s' (value = %d)"
882     ,flagfntype,(flagnamefn)(flag1),flag1
883     );
884     return 12;
885 johnpye 669
886 johnpye 1254 default:
887     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
888     return 12;
889     }
890 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
891 johnpye 997 }else{
892     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
893 johnpye 1049 return 13;
894 johnpye 997 }
895 johnpye 996
896     if(enginedata->safeeval){
897 johnpye 1246 Asc_SignalHandlerPop(SIGFPE,SIG_DFL);
898 johnpye 996 }else{
899 johnpye 997 CONSOLE_DEBUG("pop...");
900 johnpye 996 Asc_SignalHandlerPopDefault(SIGFPE);
901 johnpye 997 CONSOLE_DEBUG("...pop");
902 johnpye 991 }
903 johnpye 1142 #endif
904 johnpye 1254 }/* icopt */
905 johnpye 996
906 jpye 1337 /* optionally, specify ROO-FINDING problem */
907     if(enginedata->nbnds){
908     IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn, (void *)sys);
909     }
910 johnpye 669
911     /* -- set up the IntegratorReporter */
912 johnpye 1306 integrator_output_init(sys);
913 johnpye 669
914     /* -- store the initial values of all the stuff */
915 johnpye 1306 integrator_output_write(sys);
916     integrator_output_write_obs(sys);
917 johnpye 669
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 johnpye 1306 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 johnpye 972 asc_assert(t > t0);
928 johnpye 980
929 johnpye 979 #ifdef SOLVE_DEBUG
930 johnpye 1009 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
931 johnpye 979 #endif
932 johnpye 980
933 johnpye 1310 #ifdef ASC_SIGNAL_TRAPS
934     Asc_SignalHandlerPushDefault(SIGINT);
935     if(setjmp(g_int_env)==0) {
936     #endif
937 johnpye 669 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
938 johnpye 1310 #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 johnpye 669
946 jpye 1339 /* 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 johnpye 669 /* pass the values of everything back to the compiler */
971 johnpye 1306 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 johnpye 669
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 johnpye 1306 /* -- do something so that sys knows the values of tret, yret and ypret */
981 johnpye 669
982     /* -- store the current values of all the stuff */
983 johnpye 1306 integrator_output_write(sys);
984     integrator_output_write_obs(sys);
985 johnpye 669
986 jpye 1339 }/* loop through next sample timestep */
987    
988 johnpye 669 /* -- close the IntegratorReporter */
989 johnpye 1306 integrator_output_close(sys);
990 johnpye 669
991 johnpye 993 /* 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 jpye 1337 }else{
997     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
998 johnpye 669 }
999 johnpye 993 #endif
1000 johnpye 669
1001     /* free solution memory */
1002     N_VDestroy_Serial(yret);
1003     N_VDestroy_Serial(ypret);
1004 johnpye 980
1005 johnpye 669 /* free solver memory */
1006     IDAFree(ida_mem);
1007    
1008 johnpye 1311 if(flag < -500){
1009 johnpye 1310 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", t);
1010     return -flag;
1011     }
1012    
1013 johnpye 993 if(flag < 0){
1014     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1015 johnpye 1049 return 14;
1016 johnpye 993 }
1017    
1018     /* all done, success */
1019 johnpye 1049 return 0;
1020 johnpye 669 }
1021    
1022     /*--------------------------------------------------
1023 jpye 1337 RESIDUALS AND JACOBIAN AND IDAROOTFN
1024 johnpye 669 */
1025 johnpye 1232
1026 johnpye 1263 #if 0
1027 johnpye 1232 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 johnpye 1237 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 johnpye 1263 #endif
1057 johnpye 1237
1058 johnpye 669 /**
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 johnpye 1306 @param res_data pointer to our stuff (sys in this case).
1068 johnpye 669
1069     @return 0 on success, positive on recoverable error, and
1070     negative on unrecoverable error.
1071     */
1072 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1073 johnpye 1306 IntegratorSystem *sys;
1074 johnpye 669 IntegratorIdaData *enginedata;
1075     int i, calc_ok, is_error;
1076     struct rel_relation** relptr;
1077     double resid;
1078     char *relname;
1079 johnpye 949 #ifdef FEX_DEBUG
1080     char *varname;
1081 johnpye 991 char diffname[30];
1082 johnpye 949 #endif
1083 johnpye 669
1084 johnpye 1306 sys = (IntegratorSystem *)res_data;
1085     enginedata = integrator_ida_enginedata(sys);
1086 johnpye 669
1087 johnpye 949 #ifdef FEX_DEBUG
1088 johnpye 903 /* fprintf(stderr,"\n\n"); */
1089 johnpye 949 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1090     #endif
1091 johnpye 669
1092 johnpye 1271 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 johnpye 669 /* pass the values of everything back to the compiler */
1098 johnpye 1306 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 johnpye 669
1102 johnpye 1271 /* perform bounds checking on all variables */
1103 johnpye 1306 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1104 johnpye 1271 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1105     return 1;
1106 johnpye 669 }
1107    
1108 johnpye 949 /* evaluate each residual in the rellist */
1109 johnpye 980 is_error = 0;
1110 johnpye 949 relptr = enginedata->rellist;
1111    
1112 johnpye 1236
1113 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
1114 johnpye 968 if(enginedata->safeeval){
1115     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1116     }else{
1117 johnpye 1142 # ifdef FEX_DEBUG
1118 johnpye 968 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1119 johnpye 1142 # endif
1120 johnpye 968 Asc_SignalHandlerPushDefault(SIGFPE);
1121     }
1122    
1123 johnpye 997 if (SETJMP(g_fpe_env)==0) {
1124 johnpye 1142 #endif
1125 johnpye 1236
1126 johnpye 1263
1127 jpye 1337 for(i=0, relptr = enginedata->rellist;
1128 johnpye 669 i< enginedata->nrels && relptr != NULL;
1129     ++i, ++relptr
1130 jpye 1337 ){
1131     resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
1132 johnpye 669
1133 jpye 1337 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 johnpye 1232 }
1144 jpye 1337
1145 johnpye 1232 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 johnpye 1237 #ifdef FEX_DEBUG
1153 johnpye 1232 if(!is_error){
1154     CONSOLE_DEBUG("No NAN detected");
1155     }
1156 johnpye 1237 #endif
1157 johnpye 1232 }
1158 johnpye 1142
1159     #ifdef ASC_SIGNAL_TRAPS
1160 johnpye 669 }else{
1161 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
1162 johnpye 949 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1163     ASC_FREE(relname);
1164     is_error = 1;
1165 johnpye 669 }
1166    
1167 johnpye 968 if(enginedata->safeeval){
1168     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1169     }else{
1170     Asc_SignalHandlerPopDefault(SIGFPE);
1171     }
1172 johnpye 1142 #endif
1173 johnpye 968
1174 johnpye 1236
1175 johnpye 949 #ifdef FEX_DEBUG
1176     /* output residuals to console */
1177     CONSOLE_DEBUG("RESIDUAL OUTPUT");
1178 johnpye 991 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1179 johnpye 1306 for(i=0; i<sys->n_y; ++i){
1180     varname = var_make_name(sys->system,sys->y[i]);
1181 johnpye 991 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1182 johnpye 1306 if(sys->ydot[i]){
1183     varname = var_make_name(sys->system,sys->ydot[i]);
1184 johnpye 991 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1185 johnpye 949 }else{
1186 johnpye 991 snprintf(diffname,99,"diff(%s)",varname);
1187     fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1188 johnpye 949 }
1189     ASC_FREE(varname);
1190 johnpye 1306 relname = rel_make_name(sys->system,enginedata->rellist[i]);
1191 johnpye 979 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1192 johnpye 949 }
1193     #endif
1194    
1195     if(is_error){
1196     return 1;
1197     }
1198    
1199 johnpye 980 #ifdef FEX_DEBUG
1200 johnpye 949 CONSOLE_DEBUG("RESIDUAL OK");
1201     #endif
1202     return 0;
1203 johnpye 669 }
1204    
1205 johnpye 903 /**
1206 johnpye 945 Dense Jacobian evaluation. Only suitable for small problems!
1207 johnpye 1280 Has been seen working for problems up to around 2000 vars, FWIW.
1208 johnpye 945 */
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 johnpye 1306 IntegratorSystem *sys;
1215 johnpye 945 IntegratorIdaData *enginedata;
1216 johnpye 949 char *relname;
1217 johnpye 1009 #ifdef DJEX_DEBUG
1218 johnpye 1196 struct var_variable **varlist;
1219 johnpye 949 char *varname;
1220     #endif
1221 johnpye 946 struct rel_relation **relptr;
1222     int i;
1223     double *derivatives;
1224 johnpye 1221 struct var_variable **variables;
1225 johnpye 976 int count, j;
1226 johnpye 1237 int status, is_error = 0;
1227 johnpye 945
1228 johnpye 1306 sys = (IntegratorSystem *)jac_data;
1229     enginedata = integrator_ida_enginedata(sys);
1230 johnpye 945
1231 johnpye 1221 /* 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 johnpye 946 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1235    
1236 johnpye 945 /* pass the values of everything back to the compiler */
1237 johnpye 1306 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 johnpye 945
1241 johnpye 1271 /* perform bounds checking on all variables */
1242 johnpye 1306 if(slv_check_bounds(sys->system, 0, -1, NULL)){
1243 johnpye 1271 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
1244     return 1;
1245     }
1246    
1247 johnpye 1009 #ifdef DJEX_DEBUG
1248 johnpye 1306 varlist = slv_get_solvers_var_list(sys->system);
1249 johnpye 1196
1250 johnpye 946 /* print vars */
1251 johnpye 1306 for(i=0; i < sys->n_y; ++i){
1252     varname = var_make_name(sys->system, sys->y[i]);
1253 johnpye 1228 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
1254 johnpye 1306 asc_assert(NV_Ith_S(yy,i) == var_value(sys->y[i]));
1255 johnpye 946 ASC_FREE(varname);
1256     }
1257 johnpye 945
1258 johnpye 946 /* print derivatives */
1259 johnpye 1306 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 johnpye 946 ASC_FREE(varname);
1264     }else{
1265 johnpye 1306 varname = var_make_name(sys->system, sys->y[i]);
1266 johnpye 1223 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1267 johnpye 946 ASC_FREE(varname);
1268     }
1269     }
1270    
1271     /* print step size */
1272 johnpye 1223 CONSOLE_DEBUG("<c_j> = %g",c_j);
1273 johnpye 949 #endif
1274 johnpye 980
1275 johnpye 946 /* 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 johnpye 1221 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1283 johnpye 946
1284     if(status){
1285 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
1286 johnpye 946 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1287     ASC_FREE(relname);
1288 johnpye 1237 is_error = 1;
1289 johnpye 946 break;
1290     }
1291    
1292 johnpye 948 /* output what's going on here ... */
1293 johnpye 1009 #ifdef DJEX_DEBUG
1294 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
1295 johnpye 948 fprintf(stderr,"%d: '%s': ",i,relname);
1296     for(j=0;j<count;++j){
1297 johnpye 1306 varname = var_make_name(sys->system, variables[j]);
1298 johnpye 1221 if(var_deriv(variables[j])){
1299     fprintf(stderr," '%s'=",varname);
1300 johnpye 1306 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(sys,variables[j]));
1301 johnpye 946 }else{
1302 johnpye 1221 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1303 johnpye 946 }
1304 johnpye 948 ASC_FREE(varname);
1305     }
1306 johnpye 1228 /* relname is freed further down */
1307 johnpye 948 fprintf(stderr,"\n");
1308 johnpye 980 #endif
1309 johnpye 1221
1310 johnpye 949 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1311 johnpye 980 for(j=0; j < count; ++j){
1312 johnpye 1237 #ifdef DJEX_DEBUG
1313 johnpye 1306 varname = var_make_name(sys->system,variables[j]);
1314 johnpye 1228 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
1315     ASC_FREE(varname);
1316 johnpye 1237 #endif
1317 johnpye 1221 if(!var_deriv(variables[j])){
1318 johnpye 1237 #ifdef DJEX_DEBUG
1319     fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
1320 johnpye 1221 asc_assert(var_sindex(variables[j]) >= 0);
1321     ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1322 johnpye 1237 #endif
1323 johnpye 1221 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1324 johnpye 946 }else{
1325 johnpye 1306 DENSE_ELEM(Jac,i,integrator_ida_diffindex(sys,variables[j])) += derivatives[j] * c_j;
1326 johnpye 1237 #ifdef DJEX_DEBUG
1327     fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
1328     #endif
1329 johnpye 946 }
1330     }
1331     }
1332    
1333 johnpye 1009 #ifdef DJEX_DEBUG
1334 johnpye 1228 ASC_FREE(relname);
1335 johnpye 946 CONSOLE_DEBUG("PRINTING JAC");
1336 johnpye 948 fprintf(stderr,"\t");
1337 johnpye 1306 for(j=0; j < sys->n_y; ++j){
1338 johnpye 948 if(j)fprintf(stderr,"\t");
1339 johnpye 1306 varname = var_make_name(sys->system,sys->y[j]);
1340 johnpye 948 fprintf(stderr,"%11s",varname);
1341     ASC_FREE(varname);
1342     }
1343     fprintf(stderr,"\n");
1344     for(i=0; i < enginedata->nrels; ++i){
1345 johnpye 1306 relname = rel_make_name(sys->system, enginedata->rellist[i]);
1346 johnpye 948 fprintf(stderr,"%s\t",relname);
1347     ASC_FREE(relname);
1348    
1349 johnpye 1306 for(j=0; j < sys->n_y; ++j){
1350 johnpye 948 if(j)fprintf(stderr,"\t");
1351     fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1352 johnpye 946 }
1353     fprintf(stderr,"\n");
1354 johnpye 980 }
1355 johnpye 949 #endif
1356 johnpye 946
1357 johnpye 1237 /* test for NANs */
1358     if(!is_error){
1359     for(i=0;i< enginedata->nrels; ++i){
1360 johnpye 1306 for(j=0;j<sys->n_y;++j){
1361 johnpye 1237 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 johnpye 1306 /* if(integrator_ida_check_diffindex(sys)){
1375 johnpye 1237 is_error = 1;
1376 johnpye 1239 }*/
1377 johnpye 1237
1378     if(is_error){
1379 johnpye 946 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1380     return 1;
1381     }
1382    
1383 johnpye 1009 #ifdef DJEX_DEBUG
1384 johnpye 949 CONSOLE_DEBUG("DJEX RETURNING 0");
1385 johnpye 1228 /* ASC_PANIC("Quitting"); */
1386 johnpye 949 #endif
1387 johnpye 946 return 0;
1388 johnpye 945 }
1389    
1390     /**
1391 johnpye 912 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1392 johnpye 909
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 johnpye 1306 @param jac_data pointer to our stuff (sys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1403 johnpye 909 @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 johnpye 903 */
1407 johnpye 909 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 johnpye 1306 IntegratorSystem *sys;
1412 johnpye 909 IntegratorIdaData *enginedata;
1413     int i, j, is_error=0;
1414     struct rel_relation** relptr;
1415 johnpye 969 char *relname;
1416 johnpye 909 int status;
1417     double Jv_i;
1418 johnpye 903
1419 johnpye 1221 struct var_variable **variables;
1420 johnpye 909 double *derivatives;
1421     int count;
1422 johnpye 1196 struct var_variable **varlist;
1423 johnpye 951 #ifdef JEX_DEBUG
1424 johnpye 1196
1425 johnpye 951 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1426     #endif
1427 johnpye 909
1428 johnpye 1306 sys = (IntegratorSystem *)jac_data;
1429     enginedata = integrator_ida_enginedata(sys);
1430     varlist = slv_get_solvers_var_list(sys->system);
1431 johnpye 909
1432     /* pass the values of everything back to the compiler */
1433 johnpye 1306 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 johnpye 909 /* no real use for residuals (rr) here, I don't think? */
1437    
1438 johnpye 912 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1439 johnpye 980
1440 johnpye 979 i = NV_LENGTH_S(yy) * 2;
1441     #ifdef JEX_DEBUG
1442     CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1443     #endif
1444 johnpye 1221 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1445 johnpye 979 derivatives = ASC_NEW_ARRAY(double, i);
1446 johnpye 909
1447     /* evaluate the derivatives... */
1448     /* J = dG_dy = dF_dy + alpha * dF_dyp */
1449    
1450 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
1451 johnpye 968 Asc_SignalHandlerPushDefault(SIGFPE);
1452 johnpye 997 if (SETJMP(g_fpe_env)==0) {
1453 johnpye 1142 #endif
1454 johnpye 909 for(i=0, relptr = enginedata->rellist;
1455     i< enginedata->nrels && relptr != NULL;
1456     ++i, ++relptr
1457     ){
1458     /* get derivatives for this particular relation */
1459 johnpye 1221 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1460 johnpye 979 #ifdef JEX_DEBUG
1461     CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1462     #endif
1463 johnpye 928
1464 johnpye 951 if(status){
1465 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
1466 johnpye 951 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1467     ASC_FREE(relname);
1468     is_error = 1;
1469 johnpye 909 break;
1470     }
1471    
1472 johnpye 928 /*
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 johnpye 909 Jv_i = 0;
1479     for(j=0; j < count; ++j){
1480 johnpye 1306 /* 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 johnpye 909 if(varname){
1483 johnpye 949 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1484     ASC_FREE(varname);
1485 johnpye 909 }else{
1486 johnpye 912 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1487 johnpye 909 }
1488 johnpye 949 */
1489 johnpye 979
1490     /* we don't calculate derivatives wrt indep var */
1491     asc_assert(variables[j]>=0);
1492 johnpye 1306 if(variables[j] == sys->x) continue;
1493 johnpye 979 #ifdef JEX_DEBUG
1494 johnpye 1221 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
1495 johnpye 1263 #endif
1496 johnpye 1221 if(var_deriv(variables[j])){
1497 johnpye 1306 #define DIFFINDEX integrator_ida_diffindex(sys,variables[j])
1498 johnpye 951 #ifdef JEX_DEBUG
1499 johnpye 1221 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 johnpye 951 );
1504     #endif
1505 johnpye 1306 asc_assert(sys->ydot[DIFFINDEX]==variables[j]);
1506 johnpye 1221 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1507     #undef DIFFINDEX
1508 johnpye 928 }else{
1509 johnpye 1221 #define VARINDEX var_sindex(variables[j])
1510 johnpye 951 #ifdef JEX_DEBUG
1511 johnpye 1306 asc_assert(sys->y[VARINDEX]==variables[j]);
1512 johnpye 1221 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 johnpye 951 );
1518     #endif
1519 johnpye 1221 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1520     #undef VARINDEX
1521 johnpye 928 }
1522 johnpye 909 }
1523    
1524     NV_Ith_S(Jv,i) = Jv_i;
1525 johnpye 951 #ifdef JEX_DEBUG
1526 johnpye 979 CONSOLE_DEBUG("rel = %p",*relptr);
1527 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
1528 johnpye 951 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1529 johnpye 1221 ASC_FREE(relname);
1530 johnpye 980 return 1;
1531 johnpye 951 #endif
1532 johnpye 909 }
1533 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
1534 johnpye 909 }else{
1535 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
1536 johnpye 951 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1537     ASC_FREE(relname);
1538     is_error = 1;
1539 johnpye 909 }
1540 johnpye 968 Asc_SignalHandlerPopDefault(SIGFPE);
1541 johnpye 1142 #endif
1542 johnpye 909
1543 johnpye 951 if(is_error){
1544     CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1545     return 1;
1546     }
1547     return 0;
1548 johnpye 909 }
1549    
1550 johnpye 1012 /* 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 johnpye 1013 ){
1556     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1557     return -1;
1558     }
1559 johnpye 1012
1560 jpye 1337 /* 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 jpye 1339 #ifdef ROOT_DEBUG
1579 jpye 1337 CONSOLE_DEBUG("t = %f",tt);
1580 jpye 1339 #endif
1581 jpye 1337
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 jpye 1339 #ifdef ROOT_DEBUG
1588 jpye 1337 CONSOLE_DEBUG("gout[%d] = %f",i,gout[i]);
1589 jpye 1339 #endif
1590 jpye 1337 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 jpye 1395 case e_bnd_undefined:
1601     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid boundary type e_bnd_undefined");
1602     return 1;
1603 jpye 1337 }
1604     }
1605    
1606     return 0; /* no way to detect errors in bndman_*_eval at this stage */
1607     }
1608    
1609 jpye 1395
1610 johnpye 669 /*----------------------------------------------
1611 jpye 1395 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 johnpye 993 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1771 johnpye 992 */
1772    
1773 johnpye 1306 void integrator_ida_pcreate_jacobi(IntegratorSystem *sys){
1774 jpye 1337 IntegratorIdaData *enginedata =sys->enginedata;
1775 johnpye 993 IntegratorIdaPrecDataJacobi *precdata;
1776     precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1777 johnpye 992
1778 johnpye 1306 asc_assert(sys->n_y);
1779     precdata->PIii = N_VNew_Serial(sys->n_y);
1780 johnpye 993
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 johnpye 992 /**
1799 johnpye 993 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1800 johnpye 992
1801     'setup' function.
1802     */
1803 johnpye 993 int integrator_ida_psetup_jacobi(realtype tt,
1804 johnpye 992 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 johnpye 1221 int i, j, res;
1810 johnpye 1306 IntegratorSystem *sys;
1811 johnpye 993 IntegratorIdaData *enginedata;
1812     IntegratorIdaPrecDataJacobi *precdata;
1813     struct rel_relation **relptr;
1814    
1815 johnpye 1306 sys = (IntegratorSystem *)p_data;
1816     enginedata = sys->enginedata;
1817 johnpye 993 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1818     double *derivatives;
1819 johnpye 1221 struct var_variable **variables;
1820 johnpye 993 int count, status;
1821     char *relname;
1822 johnpye 992
1823 johnpye 993 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1824 johnpye 992
1825 johnpye 1221 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1826 johnpye 993 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 johnpye 1221 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1840 johnpye 993 if(status){
1841 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
1842 johnpye 993 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 johnpye 1221 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 johnpye 993 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1854     }
1855 johnpye 1221
1856 johnpye 993 }
1857     }
1858     #ifdef PREC_DEBUG
1859     CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1860     #endif
1861 johnpye 992 }
1862 johnpye 993
1863     if(status){
1864     CONSOLE_DEBUG("Error found when evaluating derivatives");
1865 johnpye 1221 res = 1; goto finish; /* recoverable */
1866 johnpye 993 }
1867    
1868 johnpye 1306 integrator_ida_write_incidence(sys);
1869 johnpye 993
1870 johnpye 1221 res = 0;
1871     finish:
1872 johnpye 993 ASC_FREE(variables);
1873     ASC_FREE(derivatives);
1874 johnpye 1221 return res;
1875 johnpye 992 };
1876    
1877     /**
1878 johnpye 993 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1879 johnpye 992
1880     'solve' function.
1881     */
1882 johnpye 993 int integrator_ida_psolve_jacobi(realtype tt,
1883 johnpye 992 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 johnpye 1306 IntegratorSystem *sys;
1889 johnpye 992 IntegratorIdaData *data;
1890 johnpye 993 IntegratorIdaPrecDataJacobi *precdata;
1891 johnpye 1306 sys = (IntegratorSystem *)p_data;
1892     data = sys->enginedata;
1893 johnpye 993 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1894 johnpye 992
1895 johnpye 993 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1896     N_VProd(precdata->PIii, rvec, zvec);
1897 johnpye 992 return 0;
1898     };
1899    
1900     /*----------------------------------------------
1901 johnpye 993 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 johnpye 1196 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
1936 johnpye 993 */
1937    
1938 jpye 1320 /**
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 johnpye 994
1944 jpye 1320 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 johnpye 993 /**
1979 johnpye 1309 Our task here is to write the matrices that IDA *should* be seeing. We
1980 johnpye 1306 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 johnpye 994 */
1986 johnpye 1306 int integrator_ida_write_matrix(const IntegratorSystem *sys, FILE *f, const char *type){
1987 jpye 1320 /* IntegratorIdaData *enginedata; */
1988 jpye 1325 struct SystemJacobianStruct J = {NULL,NULL,NULL,0,0};
1989 jpye 1320 int status=1;
1990 jpye 1327 mtx_region_t R;
1991 johnpye 994
1992 jpye 1331 if(type==NULL)type = "dx'/dx";
1993 johnpye 1317
1994 jpye 1331 if(0==strcmp(type,"dg/dz")){
1995     CONSOLE_DEBUG("Calculating dg/dz...");
1996 johnpye 1317 status = system_jacobian(sys->system
1997     , &system_rfilter_algeb, &system_vfilter_algeb
1998     , 1 /* safe */
1999     , &J
2000     );
2001 jpye 1331 }else if(0==strcmp(type,"dg/dx")){
2002     CONSOLE_DEBUG("Calculating dg/dx...");
2003 johnpye 1317 status = system_jacobian(sys->system
2004     , &system_rfilter_algeb, &system_vfilter_diff
2005     , 1 /* safe */
2006     , &J
2007     );
2008 jpye 1331 }else if(0==strcmp(type,"df/dx'")){
2009     CONSOLE_DEBUG("Calculating df/dx'...");
2010 johnpye 1317 status = system_jacobian(sys->system
2011     , &system_rfilter_diff, &system_vfilter_deriv
2012     , 1 /* safe */
2013     , &J
2014     );
2015 jpye 1331 }else if(0==strcmp(type,"df/dz")){
2016     CONSOLE_DEBUG("Calculating df/dz...");
2017 johnpye 1317 status = system_jacobian(sys->system
2018     , &system_rfilter_diff, &system_vfilter_algeb
2019     , 1 /* safe */
2020     , &J
2021     );
2022 jpye 1331 }else if(0==strcmp(type,"df/dx")){
2023     CONSOLE_DEBUG("Calculating df/dx...");
2024 johnpye 1317 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 jpye 1331 }else if(0==strcmp(type,"dF/dy'")){
2037 johnpye 1317 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 jpye 1331 }else if(0==strcmp(type,"dx'/dx")){
2044 jpye 1320 /* system state transfer matrix dyd'/dyd */
2045     status = integrator_ida_transfer_matrix(sys, &J);
2046 johnpye 1309 }else{
2047 johnpye 1317 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid matrix type '%s'",type);
2048     return 1;
2049 johnpye 1306 }
2050    
2051 johnpye 1317 if(status){
2052     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating matrix");
2053     }else{
2054 jpye 1327 /* 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 jpye 1331 /* note that we're not fussy about empty matrices here... */
2058 jpye 1327 mtx_write_region_mmio(f,J.M,&R);
2059 johnpye 1306 }
2060    
2061 johnpye 1317 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 johnpye 994 }
2067    
2068     /**
2069 johnpye 993 This routine outputs matrix structure in a crude text format, for the sake
2070     of debugging.
2071     */
2072 johnpye 1306 void integrator_ida_write_incidence(IntegratorSystem *sys){
2073 johnpye 993 int i, j;
2074     struct rel_relation **relptr;
2075 johnpye 1306 IntegratorIdaData *enginedata = sys->enginedata;
2076 johnpye 993 double *derivatives;
2077 johnpye 1221 struct var_variable **variables;
2078 johnpye 993 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 johnpye 1306 variables = ASC_NEW_ARRAY(struct var_variable *, sys->n_y * 2);
2087     derivatives = ASC_NEW_ARRAY(double, sys->n_y * 2);
2088 johnpye 993
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 johnpye 1306 relname = rel_make_name(sys->system, *relptr);
2096 johnpye 993
2097     /* get derivatives for this particular relation */
2098 johnpye 1221 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
2099 johnpye 993 if(status){
2100     CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
2101     ASC_FREE(relname);
2102     break;
2103     }
2104    
2105 johnpye 1263 fprintf(stderr,"%3d:%-15s:",i,relname);
2106 johnpye 993 ASC_FREE(relname);
2107    
2108     for(j=0; j<count; ++j){
2109 johnpye 1221 if(var_deriv(variables[j])){
2110 johnpye 1306 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(sys,variables[j]));
2111 johnpye 993 }else{
2112 johnpye 1221 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
2113 johnpye 993 }
2114     }
2115     fprintf(stderr,"\n");
2116     }
2117     ASC_FREE(variables);
2118     ASC_FREE(derivatives);
2119     }
2120    
2121 johnpye 1196 /* @return 0 on success */
2122     int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
2123 johnpye 1205 char *varname, *relname;
2124 johnpye 1196 struct var_variable **vlist, *var;
2125 johnpye 1205 struct rel_relation **rlist, *rel;
2126     long vlen, rlen;
2127 johnpye 1221 long i;
2128 jpye 1387 long di;
2129 johnpye 1196
2130 johnpye 1240 fprintf(fp,"THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM\n\n",sys->n_y);
2131 johnpye 1196
2132     /* if(integrator_sort_obs_vars(sys))return 10; */
2133    
2134 johnpye 1246 if(sys->y && sys->ydot){
2135     fprintf(fp,"CONTENTS OF THE 'Y' AND 'YDOT' LISTS\n\n");
2136 jpye 1387 fprintf(fp,"index\t%-15s\tydot\n","y");
2137     fprintf(fp,"-----\t%-15s\t-----\n","-----");
2138 johnpye 1246 for(i=0;i<sys->n_y;++i){
2139     varname = var_make_name(sys->system, sys->y[i]);
2140 jpye 1387 fprintf(fp,"%ld\t%-15s\t",i,varname);
2141 johnpye 1246 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 johnpye 1196 }
2151 johnpye 1246 }else{
2152     fprintf(fp,"'Y' and 'YDOT' LISTS ARE NOT SET!\n");
2153 johnpye 1196 }
2154    
2155 johnpye 1239 fprintf(fp,"\n\nCONTENTS OF THE VAR_FLAGS AND VAR_SINDEX\n\n");
2156 johnpye 1221 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
2157     fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
2158 johnpye 1196
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 johnpye 1205 for(i=0;i<vlen;i++){
2165 johnpye 1196 var = vlist[i];
2166    
2167     varname = var_make_name(sys->system, var);
2168 johnpye 1199 fprintf(fp,"%ld\t%-15s\t",i,varname);
2169 johnpye 1196
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 johnpye 1221 }else if(!var_incident(var)){
2177     // not incident
2178     fprintf(fp,"(not incident)\n");
2179 johnpye 1196 }else{
2180 johnpye 1221 if(var_deriv(var)){
2181 johnpye 1249 if(sys->y_id){
2182 jpye 1387 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 johnpye 1249 }else{
2191     fprintf(fp,".\tderiv... of??\n");
2192     }
2193 johnpye 1196 }else{
2194 johnpye 1239 fprintf(fp,"%d\t.\n",var_sindex(var));
2195 johnpye 1196 }
2196     }
2197     ASC_FREE(varname);
2198     }
2199    
2200 johnpye 1205 /* 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 johnpye 1263
2204 johnpye 1239 fprintf(fp,"\nALL RELATIONS IN THE SOLVER'S LIST (%ld)\n\n",rlen);
2205 johnpye 1205 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 johnpye 1208
2214 johnpye 1238 /* write out the derivative chains */
2215     fprintf(fp,"\nDERIVATIVE CHAINS\n");
2216 johnpye 1239 if(integrator_ida_analyse_debug(sys,stderr)){
2217 johnpye 1238 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error getting diffvars debug info");
2218     return 340;
2219     }
2220     fprintf(fp,"\n");
2221 johnpye 1263
2222 johnpye 1208 /* and lets write block debug output */
2223 johnpye 1240 system_block_debug(sys->system, fp);
2224 johnpye 1263
2225 johnpye 1196 return 0; /* success */
2226     }
2227    
2228 johnpye 993 /*----------------------------------------------
2229 johnpye 669 ERROR REPORTING
2230 johnpye 782 */
2231 johnpye 669 /**
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 johnpye 1306 IntegratorSystem *sys;
2242 johnpye 669 error_severity_t sev;
2243    
2244     /* cast back the IntegratorSystem, just in case we need it */
2245 johnpye 1306 sys = (IntegratorSystem *)eh_data;
2246 johnpye 669
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