/[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 1187 - (hide annotations) (download) (as text)
Sun Jan 21 06:03:58 2007 UTC (15 years, 6 months ago) by johnpye
File MIME type: text/x-csrc
File size: 54330 byte(s)
analyse.c can now pick up 'chains' of derivatives and pass them to the solver.
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 1183 #include "ida.h"
35    
36 johnpye 669 #include <signal.h>
37    
38     /* SUNDIALS includes */
39     #ifdef ASC_WITH_IDA
40 johnpye 913 # include <sundials/sundials_config.h>
41 johnpye 945 # include <sundials/sundials_dense.h>
42 johnpye 782 # include <ida/ida.h>
43     # include <nvector/nvector_serial.h>
44     # include <ida/ida_spgmr.h>
45 johnpye 969 # include <ida/ida_spbcgs.h>
46     # include <ida/ida_sptfqmr.h>
47 johnpye 945 # include <ida/ida_dense.h>
48 johnpye 669 # ifndef IDA_SUCCESS
49     # error "Failed to include SUNDIALS IDA header file"
50     # endif
51 johnpye 1181 #else
52     # error "Where is IDA?"
53 johnpye 669 #endif
54    
55 johnpye 994 #ifdef ASC_WITH_MMIO
56     # include <mmio.h>
57     #endif
58    
59 johnpye 1183 #include <utilities/config.h>
60     #include <utilities/ascConfig.h>
61     #include <utilities/error.h>
62     #include <utilities/ascSignal.h>
63     #include <utilities/ascPanic.h>
64     #include <compiler/instance_enum.h>
65 johnpye 1181
66 johnpye 1187 #include <solver/slv_client.h>
67 johnpye 1183 #include <solver/relman.h>
68 johnpye 1187 #ifdef ASC_IDA_NEW_ANALYSE
69     # include <solver/analyze.h>
70     # include <solver/slv_stdcalls.h>
71     #endif
72 johnpye 1181
73 johnpye 1183 #include "idalinear.h"
74    
75 johnpye 913 /*
76 johnpye 969 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
77 johnpye 913 */
78     #ifndef SUNDIALS_VERSION_MINOR
79 johnpye 921 # ifdef __GNUC__
80     # warning "GUESSING SUNDIALS VERSION 2.2"
81     # endif
82 johnpye 913 # define SUNDIALS_VERSION_MINOR 2
83     #endif
84     #ifndef SUNDIALS_VERSION_MAJOR
85     # define SUNDIALS_VERSION_MAJOR 2
86     #endif
87    
88 johnpye 669 /* check that we've got what we expect now */
89     #ifndef ASC_IDA_H
90     # error "Failed to include ASCEND IDA header file"
91     #endif
92    
93 johnpye 1009 /* #define FEX_DEBUG */
94 johnpye 997 #define JEX_DEBUG
95     #define SOLVE_DEBUG
96 johnpye 993 #define STATS_DEBUG
97     #define PREC_DEBUG
98 johnpye 951
99 johnpye 993 /**
100     Everthing that the outside world needs to know about IDA
101     */
102 johnpye 977 const IntegratorInternals integrator_ida_internals = {
103     integrator_ida_create
104     ,integrator_ida_params_default
105 johnpye 1183 #ifdef ASC_IDA_NEW_ANALYSE
106 johnpye 1068 ,integrator_ida_analyse
107     #else
108 johnpye 977 ,integrator_analyse_dae /* note, this routine is back in integrator.c */
109 johnpye 1068 #endif
110 johnpye 977 ,integrator_ida_solve
111 johnpye 1129 ,NULL /* writematrixfn */
112 johnpye 977 ,integrator_ida_free
113     ,INTEG_IDA
114     ,"IDA"
115     };
116    
117 johnpye 993 /*-------------------------------------------------------------
118     FORWARD DECLS
119     */
120    
121     /* forward dec needed for IntegratorIdaPrecFreeFn */
122     struct IntegratorIdaDataStruct;
123    
124     /* functions for allocating storage for and freeing preconditioner data */
125     typedef void IntegratorIdaPrecCreateFn(IntegratorSystem *blsys);
126     typedef void IntegratorIdaPrecFreeFn(struct IntegratorIdaDataStruct *enginedata);
127    
128 johnpye 669 /**
129 johnpye 980 Struct containing any stuff that IDA needs that doesn't fit into the
130 johnpye 669 common IntegratorSystem struct.
131     */
132 johnpye 993 typedef struct IntegratorIdaDataStruct{
133 johnpye 669 struct rel_relation **rellist; /**< NULL terminated list of rels */
134 johnpye 912 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
135 johnpye 972 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
136 johnpye 669 int nrels;
137     int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
138 johnpye 994 var_filter_t vfilter; /**< Used to filter variables from varlist in relman_diff2 etc */
139     rel_filter_t rfilter; /**< Used to filter relations from rellist (@TODO needs work) */
140 johnpye 993 void *precdata; /**< For use by the preconditioner */
141     IntegratorIdaPrecFreeFn *pfree; /**< Store instructions here on how to free precdata */
142 johnpye 669 } IntegratorIdaData;
143    
144 johnpye 1060 typedef struct IntegratorIdaPrecDJStruct{
145 johnpye 993 N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */
146     } IntegratorIdaPrecDataJacobi;
147    
148     /**
149     Hold all the function pointers associated with a particular preconditioner
150     We don't need to store the 'pfree' function here as it is allocated to the enginedata struct
151     by the pcreate function (ensures that corresponding 'free' and 'create' are always used)
152    
153     @note IDA uses a different convention for function pointer types, so no '*'.
154 johnpye 669 */
155 johnpye 1060 typedef struct IntegratorIdaPrecStruct{
156 johnpye 993 IntegratorIdaPrecCreateFn *pcreate;
157     IDASpilsPrecSetupFn psetup;
158     IDASpilsPrecSolveFn psolve;
159     } IntegratorIdaPrec;
160    
161 johnpye 669 /* residual function forward declaration */
162 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
163 johnpye 669
164 johnpye 909 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
165     , N_Vector v, N_Vector Jv, realtype c_j
166     , void *jac_data, N_Vector tmp1, N_Vector tmp2
167     );
168    
169 johnpye 669 /* error handler forward declaration */
170     void integrator_ida_error(int error_code
171     , const char *module, const char *function
172     , char *msg, void *eh_data
173     );
174    
175 johnpye 1012 /* dense jacobian evaluation for IDADense dense direct linear solver */
176 johnpye 945 int integrator_ida_djex(long int Neq, realtype tt
177     , N_Vector yy, N_Vector yp, N_Vector rr
178     , realtype c_j, void *jac_data, DenseMat Jac
179     , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
180     );
181    
182 johnpye 1012 /* sparse jacobian evaluation for ASCEND's sparse direct solver */
183     IntegratorSparseJacFn integrator_ida_sjex;
184    
185 johnpye 1060 typedef struct IntegratorIdaStatsStruct{
186 johnpye 993 long nsteps;
187     long nrevals;
188     long nlinsetups;
189     long netfails;
190     int qlast, qcur;
191     realtype hinused, hlast, hcur;
192     realtype tcur;
193     } IntegratorIdaStats;
194    
195     int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
196     void integrator_ida_write_stats(IntegratorIdaStats *stats);
197     void integrator_ida_write_incidence(IntegratorSystem *blsys);
198     /*------
199     Jacobi preconditioner -- experimental
200     */
201    
202     int integrator_ida_psetup_jacobi(realtype tt,
203 johnpye 992 N_Vector yy, N_Vector yp, N_Vector rr,
204     realtype c_j, void *prec_data,
205     N_Vector tmp1, N_Vector tmp2,
206     N_Vector tmp3
207     );
208    
209 johnpye 993 int integrator_ida_psolve_jacobi(realtype tt,
210 johnpye 992 N_Vector yy, N_Vector yp, N_Vector rr,
211     N_Vector rvec, N_Vector zvec,
212     realtype c_j, realtype delta, void *prec_data,
213     N_Vector tmp
214     );
215    
216 johnpye 993 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys);
217 johnpye 992
218 johnpye 993 void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata);
219    
220     static const IntegratorIdaPrec prec_jacobi = {
221     integrator_ida_pcreate_jacobi
222     , integrator_ida_psetup_jacobi
223     , integrator_ida_psolve_jacobi
224     };
225    
226 johnpye 669 /*-------------------------------------------------------------
227     SETUP/TEARDOWN ROUTINES
228     */
229     void integrator_ida_create(IntegratorSystem *blsys){
230     CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
231     IntegratorIdaData *enginedata;
232     enginedata = ASC_NEW(IntegratorIdaData);
233     enginedata->rellist = NULL;
234 johnpye 909 enginedata->varlist = NULL;
235 johnpye 953 enginedata->safeeval = 0;
236 johnpye 994 enginedata->vfilter.matchbits = VAR_SVAR | VAR_ACTIVE | VAR_FIXED;
237     enginedata->vfilter.matchvalue = VAR_SVAR | VAR_ACTIVE;
238 johnpye 1007 enginedata->pfree = NULL;
239 johnpye 993
240 johnpye 994 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
241     enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
242    
243 ben.allan 704 blsys->enginedata = (void *)enginedata;
244 johnpye 993
245 johnpye 942 integrator_ida_params_default(blsys);
246 ben.allan 704 }
247 johnpye 669
248     void integrator_ida_free(void *enginedata){
249     CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
250     IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
251 johnpye 993 if(d->pfree){
252     /* free the preconditioner data, whatever it happens to be */
253     (d->pfree)(enginedata);
254     }
255 johnpye 669 /* note, we don't own the rellist, so don't need to free it */
256     ASC_FREE(d);
257     }
258    
259     IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
260     IntegratorIdaData *d;
261     assert(blsys!=NULL);
262     assert(blsys->enginedata!=NULL);
263     assert(blsys->engine==INTEG_IDA);
264     d = ((IntegratorIdaData *)(blsys->enginedata));
265     return d;
266     }
267    
268     /*-------------------------------------------------------------
269 johnpye 942 PARAMETERS FOR IDA
270     */
271    
272     enum ida_parameters{
273 johnpye 945 IDA_PARAM_LINSOLVER
274 johnpye 970 ,IDA_PARAM_MAXL
275 johnpye 945 ,IDA_PARAM_AUTODIFF
276 johnpye 972 ,IDA_PARAM_CALCIC
277 johnpye 953 ,IDA_PARAM_SAFEEVAL
278 johnpye 944 ,IDA_PARAM_RTOL
279     ,IDA_PARAM_ATOL
280     ,IDA_PARAM_ATOLVECT
281 johnpye 951 ,IDA_PARAM_GSMODIFIED
282 johnpye 991 ,IDA_PARAM_MAXNCF
283 johnpye 992 ,IDA_PARAM_PREC
284 johnpye 951 ,IDA_PARAMS_SIZE
285 johnpye 942 };
286    
287     /**
288     Here the full set of parameters is defined, along with upper/lower bounds,
289     etc. The values are stuck into the blsys->params structure.
290    
291 johnpye 992 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
292     list. Then add a slv_param_*(...) statement below to define the type, description and range
293     for the new parameter.
294    
295 johnpye 942 @return 0 on success
296     */
297     int integrator_ida_params_default(IntegratorSystem *blsys){
298     asc_assert(blsys!=NULL);
299     asc_assert(blsys->engine==INTEG_IDA);
300     slv_parameters_t *p;
301     p = &(blsys->params);
302    
303     slv_destroy_parms(p);
304    
305     if(p->parms==NULL){
306     CONSOLE_DEBUG("params NULL");
307     p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
308     if(p->parms==NULL)return -1;
309     p->dynamic_parms = 1;
310     }else{
311     CONSOLE_DEBUG("params not NULL");
312     }
313    
314     /* reset the number of parameters to zero so that we can check it at the end */
315     p->num_parms = 0;
316    
317     slv_param_bool(p,IDA_PARAM_AUTODIFF
318     ,(SlvParameterInitBool){{"autodiff"
319     ,"Use auto-diff?",1
320     ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
321     }, TRUE}
322     );
323    
324 johnpye 993 slv_param_char(p,IDA_PARAM_CALCIC
325     ,(SlvParameterInitChar){{"calcic"
326     ,"Initial conditions calcuation",1
327     ,"Use specified values of ydot to solve for inital y (Y),"
328     " or use the the values of the differential variables (yd) to solve"
329     " for the pure algebraic variables (ya) along with the derivatives"
330     " of the differential variables (yddot) (YA_YDP), or else don't solve"
331     " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
332 johnpye 1133 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
333 johnpye 972 );
334    
335 johnpye 953 slv_param_bool(p,IDA_PARAM_SAFEEVAL
336     ,(SlvParameterInitBool){{"safeeval"
337     ,"Use safe evaluation?",1
338     ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
339     "throw SIGFPE errors which will then halt integration."
340     }, FALSE}
341     );
342    
343    
344 johnpye 944 slv_param_bool(p,IDA_PARAM_ATOLVECT
345     ,(SlvParameterInitBool){{"atolvect"
346     ,"Use 'ode_atol' values as specified?",1
347     ,"If TRUE, values of 'ode_atol' are taken from your model and used "
348     " in the integration. If FALSE, a scalar absolute tolerance value"
349     " is shared by all variables. See IDA manual, section 5.5.1"
350     }, TRUE }
351     );
352    
353     slv_param_real(p,IDA_PARAM_ATOL
354     ,(SlvParameterInitReal){{"atol"
355     ,"Scalar absolute error tolerance",1
356     ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
357 johnpye 997 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
358 johnpye 1141 }, 1e-5, 1e-15, 1.0e15 }
359 johnpye 944 );
360    
361     slv_param_real(p,IDA_PARAM_RTOL
362     ,(SlvParameterInitReal){{"rtol"
363     ,"Scalar relative error tolerance",1
364 johnpye 963 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
365     " it's not possible to set per-variable relative tolerances as it is"
366     " with LSODE)."
367 johnpye 997 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
368 johnpye 1141 }, 1e-4, 0, 1 }
369 johnpye 944 );
370    
371 johnpye 945 slv_param_char(p,IDA_PARAM_LINSOLVER
372     ,(SlvParameterInitChar){{"linsolver"
373     ,"Linear solver",1
374 johnpye 1012 ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"
375     " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"
376     " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"
377     " or SPTFQMR (which still need preconditioners to be implemented"
378     " before they can be very useful."
379 johnpye 1016 }, "SPGMR"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
380 johnpye 945 );
381    
382 johnpye 970 slv_param_int(p,IDA_PARAM_MAXL
383     ,(SlvParameterInitInt){{"maxl"
384     ,"Maximum Krylov dimension",0
385     ,"The maximum dimension of Krylov space used by the linear solver"
386     " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
387     " The default of 0 results in IDA using its internal default, which"
388     " is currently a value of 5."
389     }, 0, 0, 20 }
390     );
391    
392 johnpye 951 slv_param_bool(p,IDA_PARAM_GSMODIFIED
393     ,(SlvParameterInitBool){{"gsmodified"
394 johnpye 991 ,"Gram-Schmidt Orthogonalisation Scheme", 2
395 johnpye 969 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
396     " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
397 johnpye 951 }, TRUE}
398     );
399    
400 johnpye 991 slv_param_int(p,IDA_PARAM_MAXNCF
401     ,(SlvParameterInitInt){{"maxncf"
402     ,"Max nonlinear solver convergence failures per step", 2
403     ,"Maximum number of allowable nonlinear solver convergence failures"
404     " on one step. See IDA manual section 5.5.6.1."
405     }, 10,0,1000 }
406     );
407    
408 johnpye 992 slv_param_char(p,IDA_PARAM_PREC
409     ,(SlvParameterInitChar){{"prec"
410     ,"Preconditioner",1
411     ,"See IDA manual, section section 5.6.8."
412     },"NONE"}, (char *[]){"NONE","DIAG",NULL}
413     );
414    
415 johnpye 942 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
416    
417     CONSOLE_DEBUG("Created %d params", p->num_parms);
418    
419     return 0;
420     }
421    
422 johnpye 1068 /*------------------------------------------------------------------------------
423     ANALYSIS ROUTINE (new implementation)
424     */
425    
426 johnpye 1069 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
427     void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
428     IntegratorVarVisitorFn integrator_dae_classify_var;
429    
430 johnpye 1187 #ifdef ASC_IDA_NEW_ANALYSE
431 johnpye 1069 /**
432     Perform additional problem analysis to prepare problem for integration with
433     IDA.
434    
435     Note, we can assume that analyse_make_problem has already been called.
436    
437     We can also assume that the independent variable has been found.
438    
439 johnpye 1070 See emails ~Jan 8 2006.
440    
441     Note, the stuff for identifying the static and output sub-problems should
442     be part of integrator.c, not this file.
443    
444 johnpye 1069 @return 0 on success
445     @see integrator_analyse
446     */
447 johnpye 1068 int integrator_ida_analyse(struct IntegratorSystemStruct *sys){
448 johnpye 1069 struct var_variable **solversvars;
449 johnpye 1070 const struct var_variable **vlist;
450     unsigned long nsolversvars, i, j, nderivs, nvlist, nfixed;
451     struct rel_relation **solversrels;
452     unsigned long nsolversrels, dynamicstart, dynamicend;
453     const SolverDiffVarCollection *diffvars;
454    
455 johnpye 1068 struct var_variable *v;
456     struct var_variable **derivs, **derivs2;
457 johnpye 1069 char *varname;
458     struct Instance *inst;
459    
460     CONSOLE_DEBUG("NEW integrator_ida_analyse------------------>");
461 johnpye 1068
462     asc_assert(sys->engine==INTEG_IDA);
463    
464     /* get our list of derivative variables */
465 johnpye 1069 solversvars = slv_get_solvers_var_list(sys->system);
466     nsolversvars = slv_get_num_solvers_vars(sys->system);
467     CONSOLE_DEBUG("Solver has %lu vars",nsolversvars);
468    
469 johnpye 1068 derivs = ASC_NEW_ARRAY(struct var_variable *,nsolversvars);
470     if(derivs==NULL){
471 johnpye 1069 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
472     return 3;
473 johnpye 1068 }
474     j = 0;
475 johnpye 1070 nfixed = 0;
476 johnpye 1069 for(i=0; i<nsolversvars; ++i){
477     v = solversvars[i];
478     asc_assert((int)v!=0x1);
479    
480 johnpye 1070 if(var_fixed(v))++nfixed;
481    
482 johnpye 1069 /** @TODO remove this */
483     varname = var_make_name(sys->system,v);
484    
485 johnpye 1068 /* the solver has marked derivs already, so just test */
486     if(var_deriv(v)){
487 johnpye 1069 CONSOLE_DEBUG("Found derivative var '%s'",varname);
488 johnpye 1068 derivs[j++] = v;
489     }
490     }
491     nderivs = j;
492     if(nderivs==0){
493     ERROR_REPORTER_HERE(ASC_USER_ERROR,"System is not a dynamic problem: contains no derivatives");
494     return 1;
495     }
496    
497     /* shrink the list down to size */
498     derivs2 = derivs;
499     derivs = ASC_NEW_ARRAY(struct var_variable *,nderivs);
500     memcpy(derivs,derivs2,nderivs*sizeof(struct var_variable *));
501     ASC_FREE(derivs2);
502    
503     CONSOLE_DEBUG("FOUND %lu DERIV VARS",nderivs);
504    
505     /* partition into static, dynamic and output problems */
506 johnpye 1069 CONSOLE_DEBUG("Block partitioning system...");
507     if(slv_block_partition(sys->system)){
508     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to block-partition system");
509     return 1;
510     }
511 johnpye 1068
512 johnpye 1070 solversrels = slv_get_solvers_rel_list(sys->system);
513     nsolversrels = slv_get_num_solvers_rels(sys->system);
514     CONSOLE_DEBUG("System has %lu rels",nsolversrels);
515    
516     dynamicstart = nsolversrels;
517     dynamicend = 0;
518     for(i=0; i<nsolversrels; ++i){
519     vlist = rel_incidence_list(solversrels[i]);
520     nvlist = rel_n_incidences(solversrels[i]);
521     for(j=0; j<nvlist; ++j){
522     if(dynamicstart > i && var_deriv(vlist[j]) ){
523     CONSOLE_DEBUG("Start dynamic problem in rel %lu",i);
524     dynamicstart = i;
525     }
526     if(dynamicend < i && var_deriv(vlist[j])){
527     CONSOLE_DEBUG("End dynamic problem with rel %lu",i);
528     dynamicend = i + 1;
529     }
530     }
531     }
532    
533     if(dynamicstart > 0){
534     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Static problem is not implemented yet");
535     return 2;
536     }
537    
538     if(dynamicend < nsolversrels){
539     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Output problem is not implemented yet");
540     return 3;
541     }
542    
543 johnpye 1068 /* raise error if any of the above are non-square */
544 johnpye 1070 if(dynamicend - dynamicstart != nsolversvars - nfixed - nderivs){
545     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Dynamic problem is not square:"
546     " %lu rels != %lu vars - %lu fixed - %lu derivs"
547     ,(dynamicend-dynamicstart),nsolversvars,nfixed,nderivs
548     );
549     return 4;
550     }
551     CONSOLE_DEBUG("Dynamic problem is square with %lu rels",dynamicend - dynamicstart);
552 johnpye 1068
553     /* get our list of differential variables */
554 johnpye 1070 /* analyse_generate_diffvars(sys->system); */
555     diffvars = analyse_get_diffvars(sys->system);
556     if(diffvars==NULL){
557     ERROR_REPORTER_HERE(ASC_PROG_ERR,"diffvars structure is NULL");
558     return 5;
559     }
560 johnpye 1068
561 johnpye 1070 /* for(i=0; i<diffvars->n; ++i){
562     diffvars->seqs[i] */
563    
564 johnpye 1068 /* get our list of algebraic varibles */
565 johnpye 1070
566    
567 johnpye 1068
568     /* set up the static problem */
569     /* - rel, var (etc) lists */
570     /* - block decomposition */
571 johnpye 1070 CONSOLE_DEBUG("Skipping the static problem setup");
572 johnpye 1068
573     /* set up the output problem */
574     /* - rel, var (etc) lists */
575     /* - block decomposition */
576 johnpye 1070 CONSOLE_DEBUG("Skipping the output problem setup");
577 johnpye 1068
578     /* set up the dynamic problem */
579 johnpye 1070 CONSOLE_DEBUG("Setting up the dynamic problem");
580 johnpye 1068 /* - 'y' list as [ya|yd] */
581     /* - sparsity pattern for dF/dy and dF/dy' */
582     /* - sparsity pattern for union of above */
583     /* - block decomposition based on above */
584     /* - block decomposition results in reordering of y and y' */
585     /* - boundaries (optional) */
586 johnpye 1070 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Implementation incomplete");
587 johnpye 1068 return 1;
588     }
589 johnpye 1187 #endif
590 johnpye 1068
591 johnpye 942 /*-------------------------------------------------------------
592 johnpye 669 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
593     */
594    
595 johnpye 1060 /*static double div1(double a, double b){
596 johnpye 997 return a/b;
597 johnpye 1060 }*/
598 johnpye 997
599 johnpye 985 typedef int IdaFlagFn(void *,int *);
600 johnpye 991 typedef char *IdaFlagNameFn(int);
601 johnpye 985
602 johnpye 1049 /* return 0 on success */
603 johnpye 669 int integrator_ida_solve(
604     IntegratorSystem *blsys
605     , unsigned long start_index
606     , unsigned long finish_index
607     ){
608     void *ida_mem;
609 johnpye 986 int size, flag, flag1, t_index;
610 johnpye 972 realtype t0, reltol, abstol, t, tret, tout1;
611 johnpye 993 N_Vector y0, yp0, abstolvect, ypret, yret, id;
612 johnpye 669 IntegratorIdaData *enginedata;
613 johnpye 945 char *linsolver;
614 johnpye 971 int maxl;
615 johnpye 985 IdaFlagFn *flagfn;
616     IdaFlagNameFn *flagnamefn;
617 johnpye 986 const char *flagfntype;
618 johnpye 993 char *pname = NULL;
619     char *varname;
620     int i;
621     const IntegratorIdaPrec *prec = NULL;
622     int icopt; /* initial conditions strategy */
623 johnpye 669
624     CONSOLE_DEBUG("STARTING IDA...");
625    
626     enginedata = integrator_ida_enginedata(blsys);
627 johnpye 953
628     enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);
629 johnpye 669 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
630    
631     /* store reference to list of relations (in enginedata) */
632     enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
633     enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
634 johnpye 909 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
635 johnpye 972 enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);
636    
637 johnpye 669 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
638 johnpye 719 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
639 johnpye 669 size = blsys->n_y;
640    
641     if(enginedata->nrels!=size){
642     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
643 johnpye 1049 return 1; /* failure */
644 johnpye 669 }
645    
646     /* retrieve initial values from the system */
647 johnpye 980
648 johnpye 944 /** @TODO fix this, the starting time != first sample */
649     t0 = integrator_get_t(blsys);
650     CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
651 johnpye 669
652 johnpye 928 CONSOLE_DEBUG("RETRIEVING y0");
653    
654 johnpye 669 y0 = N_VNew_Serial(size);
655     integrator_get_y(blsys,NV_DATA_S(y0));
656    
657 johnpye 991 #ifdef SOLVE_DEBUG
658 johnpye 928 CONSOLE_DEBUG("RETRIEVING yp0");
659 johnpye 991 #endif
660 johnpye 928
661 johnpye 669 yp0 = N_VNew_Serial(size);
662     integrator_get_ydot(blsys,NV_DATA_S(yp0));
663    
664 johnpye 991 #ifdef SOLVE_DEBUG
665 johnpye 669 N_VPrint_Serial(yp0);
666     CONSOLE_DEBUG("yp0 is at %p",&yp0);
667 johnpye 991 #endif
668 johnpye 669
669     /* create IDA object */
670     ida_mem = IDACreate();
671    
672 johnpye 980 /* relative error tolerance */
673 johnpye 944 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
674 johnpye 949 CONSOLE_DEBUG("rtol = %8.2e",reltol);
675 johnpye 669
676     /* allocate internal memory */
677 johnpye 944 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
678     /* vector of absolute tolerances */
679     CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
680     abstolvect = N_VNew_Serial(size);
681     integrator_get_atol(blsys,NV_DATA_S(abstolvect));
682    
683     flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
684 johnpye 949
685     N_VDestroy_Serial(abstolvect);
686 johnpye 944 }else{
687     /* scalar absolute tolerance (one value for all) */
688 johnpye 1023 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
689 johnpye 944 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
690     flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
691     }
692    
693 johnpye 669 if(flag==IDA_MEM_NULL){
694     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
695 johnpye 1049 return 2;
696 johnpye 669 }else if(flag==IDA_MEM_FAIL){
697     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
698 johnpye 1049 return 3;
699 johnpye 669 }else if(flag==IDA_ILL_INPUT){
700     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
701 johnpye 1049 return 4;
702 johnpye 669 }/* else success */
703    
704     /* set optional inputs... */
705     IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
706     IDASetRdata(ida_mem, (void *)blsys);
707 johnpye 980 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
708 johnpye 950 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
709     IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
710     if(integrator_get_minstep(blsys)>0){
711 johnpye 951 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
712 johnpye 950 }
713 johnpye 991
714     CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
715     IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
716    
717 johnpye 669 /* there's no capability for setting *minimum* step size in IDA */
718    
719 johnpye 912
720 johnpye 669 /* attach linear solver module, using the default value of maxl */
721 johnpye 945 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
722 johnpye 946 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
723 johnpye 1012 if(strcmp(linsolver,"ASCEND")==0){
724     CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",size);
725     IDAASCEND(ida_mem,size);
726     IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *)blsys);
727    
728     flagfntype = "IDAASCEND";
729     flagfn = &IDAASCENDGetLastFlag;
730     flagnamefn = &IDAASCENDGetReturnFlagName;
731    
732     }else if(strcmp(linsolver,"DENSE")==0){
733 johnpye 969 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
734     flag = IDADense(ida_mem, size);
735     switch(flag){
736     case IDADENSE_SUCCESS: break;
737 johnpye 1049 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 5;
738     case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 5;
739     case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 5;
740     default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 5;
741 johnpye 969 }
742 johnpye 980
743 johnpye 969 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
744     CONSOLE_DEBUG("USING AUTODIFF");
745     flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
746     switch(flag){
747     case IDADENSE_SUCCESS: break;
748 johnpye 1049 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
749 johnpye 969 }
750     }else{
751     CONSOLE_DEBUG("USING NUMERICAL DIFF");
752     }
753 johnpye 985
754 johnpye 986 flagfntype = "IDADENSE";
755 johnpye 985 flagfn = &IDADenseGetLastFlag;
756     flagnamefn = &IDADenseGetReturnFlagName;
757 johnpye 969 }else{
758     /* remaining methods are all SPILS */
759     CONSOLE_DEBUG("IDA SPILS");
760    
761 johnpye 971 maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
762     CONSOLE_DEBUG("maxl = %d",maxl);
763    
764 johnpye 993 /* what preconditioner? */
765 johnpye 992 pname = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_PREC);
766     if(strcmp(pname,"NONE")==0){
767 johnpye 993 prec = NULL;
768     }else if(strcmp(pname,"JACOBI")==0){
769     prec = &prec_jacobi;
770     }else{
771     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
772 johnpye 1049 return 7;
773 johnpye 992 }
774    
775 johnpye 993 /* which SPILS linear solver? */
776 johnpye 969 if(strcmp(linsolver,"SPGMR")==0){
777     CONSOLE_DEBUG("IDA SPGMR");
778 johnpye 971 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
779 johnpye 969 }else if(strcmp(linsolver,"SPBCG")==0){
780     CONSOLE_DEBUG("IDA SPBCG");
781 johnpye 971 flag = IDASpbcg(ida_mem, maxl);
782 johnpye 969 }else if(strcmp(linsolver,"SPTFQMR")==0){
783     CONSOLE_DEBUG("IDA SPTFQMR");
784 johnpye 971 flag = IDASptfqmr(ida_mem,maxl);
785 johnpye 969 }else{
786     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
787 johnpye 1049 return 8;
788 johnpye 980 }
789    
790 johnpye 993 if(prec){
791     /* assign the preconditioner to the linear solver */
792     (prec->pcreate)(blsys);
793     IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve,(void *)blsys);
794     CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
795     }else{
796     CONSOLE_DEBUG("No preconditioner");
797     }
798 johnpye 992
799 johnpye 986 flagfntype = "IDASPILS";
800 johnpye 985 flagfn = &IDASpilsGetLastFlag;
801     flagnamefn = &IDASpilsGetReturnFlagName;
802    
803 johnpye 942 if(flag==IDASPILS_MEM_NULL){
804     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
805 johnpye 1049 return 9;
806 johnpye 945 }else if(flag==IDASPILS_MEM_FAIL){
807     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
808 johnpye 1049 return 9;
809 johnpye 942 }/* else success */
810 johnpye 945
811     /* assign the J*v function */
812     if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
813     CONSOLE_DEBUG("USING AUTODIFF");
814     flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
815     if(flag==IDASPILS_MEM_NULL){
816     ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
817 johnpye 1049 return 10;
818 johnpye 945 }else if(flag==IDASPILS_LMEM_NULL){
819     ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
820 johnpye 1049 return 10;
821 johnpye 945 }/* else success */
822     }else{
823     CONSOLE_DEBUG("USING NUMERICAL DIFF");
824     }
825 johnpye 951
826 johnpye 969 if(strcmp(linsolver,"SPGMR")==0){
827     /* select Gram-Schmidt orthogonalisation */
828     if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){
829     CONSOLE_DEBUG("USING MODIFIED GS");
830     flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
831     if(flag!=IDASPILS_SUCCESS){
832     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
833 johnpye 1049 return 11;
834 johnpye 969 }
835     }else{
836     CONSOLE_DEBUG("USING CLASSICAL GS");
837     flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
838     if(flag!=IDASPILS_SUCCESS){
839     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
840 johnpye 1049 return 11;
841 johnpye 969 }
842 johnpye 951 }
843     }
844 johnpye 945 }
845 johnpye 669
846 johnpye 909 /* set linear solver optional inputs...
847 johnpye 980 ...nothing here at the moment...
848 johnpye 909 */
849    
850 johnpye 972 /* calculate initial conditions */
851 johnpye 993 icopt = 0;
852     if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"Y")==0){
853     CONSOLE_DEBUG("Solving initial conditions using values of yddot");
854     icopt = IDA_Y_INIT;
855     asc_assert(icopt!=0);
856     }else if(strcmp(SLV_PARAM_CHAR(&blsys->params,IDA_PARAM_CALCIC),"YA_YDP")==0){
857     CONSOLE_DEBUG("Solving initial conditions using values of yd");
858     icopt = IDA_YA_YDP_INIT;
859     asc_assert(icopt!=0);
860     id = N_VNew_Serial(blsys->n_y);
861     for(i=0; i < blsys->n_y; ++i){
862     if(blsys->ydot[i] == NULL){
863     NV_Ith_S(id,i) = 0.0;
864     varname = var_make_name(blsys->system,blsys->y[i]);
865     CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
866     ASC_FREE(varname);
867     }else{
868     CONSOLE_DEBUG("y[%d] is differential",i);
869     NV_Ith_S(id,i) = 1.0;
870     }
871     }
872     IDASetId(ida_mem, id);
873     N_VDestroy_Serial(id);
874     }else{
875     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
876     }
877 johnpye 913
878 johnpye 993 if(icopt){
879 johnpye 972 blsys->currentstep=0;
880 johnpye 1016 t_index=start_index + 1;
881 johnpye 972 tout1 = samplelist_get(blsys->samples, t_index);
882 johnpye 944
883 johnpye 972 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
884 johnpye 913
885 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
886 johnpye 996 /* catch SIGFPE if desired to */
887     if(enginedata->safeeval){
888 johnpye 997 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
889 johnpye 996 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
890     }else{
891     #ifdef FEX_DEBUG
892 johnpye 997 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
893 johnpye 996 #endif
894     Asc_SignalHandlerPushDefault(SIGFPE);
895     }
896 johnpye 1009 if (setjmp(g_fpe_env)==0) {
897 johnpye 1142 #endif
898 johnpye 997
899 johnpye 1009 //CONSOLE_DEBUG("Raising signal...");
900     //CONSOLE_DEBUG("1/0 = %f", div1(1.0,0.0));
901     //CONSOLE_DEBUG("Still here...");
902 johnpye 996
903     /* correct initial values, given derivatives */
904 johnpye 993 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
905 johnpye 996 /* note the new API from version 2.3 and onwards */
906     flag = IDACalcIC(ida_mem, icopt, tout1);
907 johnpye 993 # else
908 johnpye 996 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
909 johnpye 993 # endif
910 johnpye 669
911 johnpye 996 switch(flag){
912     case IDA_SUCCESS:
913     CONSOLE_DEBUG("Initial conditions solved OK");
914     break;
915 johnpye 986
916 johnpye 996 case IDA_LSETUP_FAIL:
917     case IDA_LINIT_FAIL:
918     case IDA_LSOLVE_FAIL:
919     case IDA_NO_RECOVERY:
920     flag1 = -999;
921     flag = (flagfn)(ida_mem,&flag1);
922     if(flag){
923     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);
924 johnpye 1049 return 12;
925 johnpye 996 }
926     ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
927 johnpye 1049 return 12;
928 johnpye 669
929 johnpye 996 default:
930     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
931 johnpye 1049 return 12;
932 johnpye 996 }
933 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
934 johnpye 997 }else{
935     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
936 johnpye 1049 return 13;
937 johnpye 997 }
938 johnpye 996
939     if(enginedata->safeeval){
940     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
941     }else{
942 johnpye 997 CONSOLE_DEBUG("pop...");
943 johnpye 996 Asc_SignalHandlerPopDefault(SIGFPE);
944 johnpye 997 CONSOLE_DEBUG("...pop");
945 johnpye 991 }
946 johnpye 1142 #endif
947 johnpye 996
948 johnpye 972 }
949    
950 johnpye 669 /* optionally, specify ROO-FINDING PROBLEM */
951    
952     /* -- set up the IntegratorReporter */
953     integrator_output_init(blsys);
954    
955     /* -- store the initial values of all the stuff */
956     integrator_output_write(blsys);
957     integrator_output_write_obs(blsys);
958    
959     /* specify where the returned values should be stored */
960     yret = y0;
961     ypret = yp0;
962    
963     /* advance solution in time, return values as yret and derivatives as ypret */
964     blsys->currentstep=1;
965 johnpye 972 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
966 johnpye 669 t = samplelist_get(blsys->samples, t_index);
967 johnpye 972 t0 = integrator_get_t(blsys);
968     asc_assert(t > t0);
969 johnpye 980
970 johnpye 979 #ifdef SOLVE_DEBUG
971 johnpye 1009 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
972 johnpye 979 #endif
973 johnpye 980
974 johnpye 669 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
975    
976     /* pass the values of everything back to the compiler */
977     integrator_set_t(blsys, (double)tret);
978     integrator_set_y(blsys, NV_DATA_S(yret));
979     integrator_set_ydot(blsys, NV_DATA_S(ypret));
980    
981     if(flag<0){
982     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
983     break;
984     }
985    
986     /* -- do something so that blsys knows the values of tret, yret and ypret */
987    
988     /* -- store the current values of all the stuff */
989     integrator_output_write(blsys);
990     integrator_output_write_obs(blsys);
991     }
992    
993     /* -- close the IntegratorReporter */
994     integrator_output_close(blsys);
995    
996 johnpye 993 /* get optional outputs */
997     #ifdef STATS_DEBUG
998     IntegratorIdaStats stats;
999     if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)){
1000     integrator_ida_write_stats(&stats);
1001 johnpye 669 }
1002 johnpye 993 #endif
1003 johnpye 669
1004     /* free solution memory */
1005     N_VDestroy_Serial(yret);
1006     N_VDestroy_Serial(ypret);
1007 johnpye 980
1008 johnpye 669 /* free solver memory */
1009     IDAFree(ida_mem);
1010    
1011 johnpye 993 if(flag < 0){
1012     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1013 johnpye 1049 return 14;
1014 johnpye 993 }
1015    
1016     /* all done, success */
1017 johnpye 1049 return 0;
1018 johnpye 669 }
1019    
1020     /*--------------------------------------------------
1021     RESIDUALS AND JACOBIAN
1022     */
1023     /**
1024     Function to evaluate system residuals, in the form required for IDA.
1025    
1026     Given tt, yy and yp, we need to evaluate and return rr.
1027    
1028     @param tt current value of indep variable (time)
1029     @param yy current values of dependent variable vector
1030     @param yp current values of derivatives of dependent variables
1031     @param rr the output residual vector (is we're returning data to)
1032     @param res_data pointer to our stuff (blsys in this case).
1033    
1034     @return 0 on success, positive on recoverable error, and
1035     negative on unrecoverable error.
1036     */
1037 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1038 johnpye 669 IntegratorSystem *blsys;
1039     IntegratorIdaData *enginedata;
1040     int i, calc_ok, is_error;
1041     struct rel_relation** relptr;
1042     double resid;
1043     char *relname;
1044 johnpye 949 #ifdef FEX_DEBUG
1045     char *varname;
1046 johnpye 991 char diffname[30];
1047 johnpye 949 #endif
1048 johnpye 669
1049     blsys = (IntegratorSystem *)res_data;
1050     enginedata = integrator_ida_enginedata(blsys);
1051    
1052 johnpye 949 #ifdef FEX_DEBUG
1053 johnpye 903 /* fprintf(stderr,"\n\n"); */
1054 johnpye 949 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1055     #endif
1056 johnpye 669
1057     /* pass the values of everything back to the compiler */
1058     integrator_set_t(blsys, (double)tt);
1059     integrator_set_y(blsys, NV_DATA_S(yy));
1060     integrator_set_ydot(blsys, NV_DATA_S(yp));
1061    
1062     if(NV_LENGTH_S(rr)!=enginedata->nrels){
1063     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
1064     return -1; /* unrecoverable */
1065     }
1066    
1067 johnpye 949 /**
1068     @TODO does this function (fex) do bounds checking already?
1069     */
1070    
1071     /* evaluate each residual in the rellist */
1072 johnpye 980 is_error = 0;
1073 johnpye 949 relptr = enginedata->rellist;
1074    
1075 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
1076 johnpye 968 if(enginedata->safeeval){
1077     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1078     }else{
1079 johnpye 1142 # ifdef FEX_DEBUG
1080 johnpye 968 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1081 johnpye 1142 # endif
1082 johnpye 968 Asc_SignalHandlerPushDefault(SIGFPE);
1083     }
1084    
1085 johnpye 997 if (SETJMP(g_fpe_env)==0) {
1086 johnpye 1142 #endif
1087 johnpye 669 for(i=0, relptr = enginedata->rellist;
1088     i< enginedata->nrels && relptr != NULL;
1089     ++i, ++relptr
1090     ){
1091     resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
1092    
1093     NV_Ith_S(rr,i) = resid;
1094     if(!calc_ok){
1095 johnpye 949 relname = rel_make_name(blsys->system, *relptr);
1096 johnpye 951 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1097 johnpye 949 ASC_FREE(relname);
1098 johnpye 669 /* presumable some output already made? */
1099     is_error = 1;
1100 johnpye 969 }/*else{
1101 johnpye 968 CONSOLE_DEBUG("Calc OK");
1102 johnpye 969 }*/
1103 johnpye 669 }
1104 johnpye 1142
1105     #ifdef ASC_SIGNAL_TRAPS
1106 johnpye 669 }else{
1107 johnpye 949 relname = rel_make_name(blsys->system, *relptr);
1108     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1109     ASC_FREE(relname);
1110     is_error = 1;
1111 johnpye 669 }
1112    
1113 johnpye 968 if(enginedata->safeeval){
1114     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1115     }else{
1116     Asc_SignalHandlerPopDefault(SIGFPE);
1117     }
1118 johnpye 1142 #endif
1119 johnpye 968
1120 johnpye 949 #ifdef FEX_DEBUG
1121     /* output residuals to console */
1122     CONSOLE_DEBUG("RESIDUAL OUTPUT");
1123 johnpye 991 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1124 johnpye 949 for(i=0; i<blsys->n_y; ++i){
1125     varname = var_make_name(blsys->system,blsys->y[i]);
1126 johnpye 991 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1127 johnpye 949 if(blsys->ydot[i]){
1128     varname = var_make_name(blsys->system,blsys->ydot[i]);
1129 johnpye 991 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1130 johnpye 949 }else{
1131 johnpye 991 snprintf(diffname,99,"diff(%s)",varname);
1132     fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1133 johnpye 949 }
1134     ASC_FREE(varname);
1135     relname = rel_make_name(blsys->system,enginedata->rellist[i]);
1136 johnpye 979 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1137 johnpye 949 }
1138     #endif
1139    
1140     if(is_error){
1141     return 1;
1142     }
1143    
1144 johnpye 980 #ifdef FEX_DEBUG
1145 johnpye 949 CONSOLE_DEBUG("RESIDUAL OK");
1146     #endif
1147     return 0;
1148 johnpye 669 }
1149    
1150 johnpye 903 /**
1151 johnpye 945 Dense Jacobian evaluation. Only suitable for small problems!
1152     */
1153     int integrator_ida_djex(long int Neq, realtype tt
1154     , N_Vector yy, N_Vector yp, N_Vector rr
1155     , realtype c_j, void *jac_data, DenseMat Jac
1156     , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1157     ){
1158     IntegratorSystem *blsys;
1159     IntegratorIdaData *enginedata;
1160 johnpye 949 char *relname;
1161 johnpye 1009 #ifdef DJEX_DEBUG
1162 johnpye 949 char *varname;
1163     #endif
1164 johnpye 946 int status;
1165     struct rel_relation **relptr;
1166     int i;
1167     double *derivatives;
1168     int *variables;
1169 johnpye 976 int count, j;
1170     long var_yindex;
1171 johnpye 945
1172     blsys = (IntegratorSystem *)jac_data;
1173     enginedata = integrator_ida_enginedata(blsys);
1174    
1175 johnpye 946 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1176     variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
1177     derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1178    
1179 johnpye 945 /* pass the values of everything back to the compiler */
1180     integrator_set_t(blsys, (double)tt);
1181     integrator_set_y(blsys, NV_DATA_S(yy));
1182     integrator_set_ydot(blsys, NV_DATA_S(yp));
1183    
1184 johnpye 1009 #ifdef DJEX_DEBUG
1185 johnpye 946 /* print vars */
1186     for(i=0; i < blsys->n_y; ++i){
1187     varname = var_make_name(blsys->system, blsys->y[i]);
1188 johnpye 948 CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
1189 johnpye 946 ASC_FREE(varname);
1190     }
1191 johnpye 945
1192 johnpye 946 /* print derivatives */
1193     for(i=0; i < blsys->n_y; ++i){
1194     if(blsys->ydot[i]){
1195     varname = var_make_name(blsys->system, blsys->ydot[i]);
1196 johnpye 948 CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
1197 johnpye 946 ASC_FREE(varname);
1198     }else{
1199     varname = var_make_name(blsys->system, blsys->y[i]);
1200     CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
1201     ASC_FREE(varname);
1202     }
1203     }
1204    
1205     /* print step size */
1206     CONSOLE_DEBUG("<c_j> = %f",c_j);
1207 johnpye 949 #endif
1208 johnpye 980
1209 johnpye 946 /* build up the dense jacobian matrix... */
1210     status = 0;
1211     for(i=0, relptr = enginedata->rellist;
1212     i< enginedata->nrels && relptr != NULL;
1213     ++i, ++relptr
1214     ){
1215     /* get derivatives for this particular relation */
1216 johnpye 994 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1217 johnpye 946
1218     if(status){
1219     relname = rel_make_name(blsys->system, *relptr);
1220     CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1221     ASC_FREE(relname);
1222     break;
1223     }
1224    
1225 johnpye 948 /* output what's going on here ... */
1226 johnpye 1009 #ifdef DJEX_DEBUG
1227 johnpye 948 relname = rel_make_name(blsys->system, *relptr);
1228     CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
1229     fprintf(stderr,"%d: '%s': ",i,relname);
1230     ASC_FREE(relname);
1231     for(j=0;j<count;++j){
1232 johnpye 946 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1233 johnpye 948 var_yindex = blsys->y_id[variables[j]];
1234     if(var_yindex >=0){
1235 johnpye 976 fprintf(stderr," var[%d]='%s'=y[%ld]",variables[j],varname,var_yindex);
1236 johnpye 946 }else{
1237 johnpye 976 fprintf(stderr," var[%d]='%s'=ydot[%ld]",variables[j],varname,-var_yindex-1);
1238 johnpye 946 }
1239 johnpye 948 ASC_FREE(varname);
1240     }
1241     fprintf(stderr,"\n");
1242 johnpye 980 #endif
1243 johnpye 949 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1244 johnpye 980 for(j=0; j < count; ++j){
1245 johnpye 948 var_yindex = blsys->y_id[variables[j]];
1246 johnpye 980 /* the SUNDIALS headers seem not to store 'N' on Windows */
1247 johnpye 996 ASC_ASSERT_RANGE(var_yindex, -blsys->n_y, blsys->n_y);
1248    
1249 johnpye 946 if(var_yindex >= 0){
1250 johnpye 948 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
1251 johnpye 946 DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
1252     }else{
1253 johnpye 948 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
1254 johnpye 946 DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
1255     }
1256     }
1257     }
1258    
1259 johnpye 1009 #ifdef DJEX_DEBUG
1260 johnpye 946 CONSOLE_DEBUG("PRINTING JAC");
1261 johnpye 948 fprintf(stderr,"\t");
1262     for(j=0; j < blsys->n_y; ++j){
1263     if(j)fprintf(stderr,"\t");
1264     varname = var_make_name(blsys->system,blsys->y[j]);
1265     fprintf(stderr,"%11s",varname);
1266     ASC_FREE(varname);
1267     }
1268     fprintf(stderr,"\n");
1269     for(i=0; i < enginedata->nrels; ++i){
1270     relname = rel_make_name(blsys->system, enginedata->rellist[i]);
1271     fprintf(stderr,"%s\t",relname);
1272     ASC_FREE(relname);
1273    
1274 johnpye 946 for(j=0; j < blsys->n_y; ++j){
1275 johnpye 948 if(j)fprintf(stderr,"\t");
1276     fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1277 johnpye 946 }
1278     fprintf(stderr,"\n");
1279 johnpye 980 }
1280 johnpye 949 #endif
1281 johnpye 946
1282     if(status){
1283     ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
1284     return 1;
1285     }
1286    
1287 johnpye 1009 #ifdef DJEX_DEBUG
1288 johnpye 949 CONSOLE_DEBUG("DJEX RETURNING 0");
1289     #endif
1290 johnpye 946 return 0;
1291 johnpye 945 }
1292    
1293     /**
1294 johnpye 912 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1295 johnpye 909
1296     Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
1297    
1298     @param tt current value of the independent variable (time, t)
1299     @param yy current value of the dependent variable vector, y(t).
1300     @param yp current value of y'(t).
1301     @param rr current value of the residual vector F(t, y, y').
1302     @param v the vector by which the Jacobian must be multiplied to the right.
1303     @param Jv the output vector computed
1304     @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
1305     @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
1306     @param tmp1 @see tmp2
1307     @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.
1308     @return 0 on success
1309 johnpye 903 */
1310 johnpye 909 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
1311     , N_Vector v, N_Vector Jv, realtype c_j
1312     , void *jac_data, N_Vector tmp1, N_Vector tmp2
1313     ){
1314     IntegratorSystem *blsys;
1315     IntegratorIdaData *enginedata;
1316     int i, j, is_error=0;
1317     struct rel_relation** relptr;
1318 johnpye 969 char *relname;
1319 johnpye 909 int status;
1320     double Jv_i;
1321 johnpye 976 long var_yindex;
1322 johnpye 903
1323 johnpye 909 int *variables;
1324     double *derivatives;
1325     int count;
1326 johnpye 951 #ifdef JEX_DEBUG
1327     CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1328     #endif
1329 johnpye 909
1330     blsys = (IntegratorSystem *)jac_data;
1331     enginedata = integrator_ida_enginedata(blsys);
1332    
1333     /* pass the values of everything back to the compiler */
1334     integrator_set_t(blsys, (double)tt);
1335     integrator_set_y(blsys, NV_DATA_S(yy));
1336     integrator_set_ydot(blsys, NV_DATA_S(yp));
1337     /* no real use for residuals (rr) here, I don't think? */
1338    
1339 johnpye 912 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1340 johnpye 980
1341 johnpye 979 i = NV_LENGTH_S(yy) * 2;
1342     #ifdef JEX_DEBUG
1343     CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1344     #endif
1345     variables = ASC_NEW_ARRAY(int, i);
1346     derivatives = ASC_NEW_ARRAY(double, i);
1347 johnpye 909
1348     /* evaluate the derivatives... */
1349     /* J = dG_dy = dF_dy + alpha * dF_dyp */
1350    
1351 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
1352 johnpye 968 Asc_SignalHandlerPushDefault(SIGFPE);
1353 johnpye 997 if (SETJMP(g_fpe_env)==0) {
1354 johnpye 1142 #endif
1355 johnpye 909 for(i=0, relptr = enginedata->rellist;
1356     i< enginedata->nrels && relptr != NULL;
1357     ++i, ++relptr
1358     ){
1359     /* get derivatives for this particular relation */
1360 johnpye 994 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1361 johnpye 979 #ifdef JEX_DEBUG
1362     CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1363     #endif
1364 johnpye 928
1365 johnpye 951 if(status){
1366     relname = rel_make_name(blsys->system, *relptr);
1367     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1368     ASC_FREE(relname);
1369     is_error = 1;
1370 johnpye 909 break;
1371     }
1372    
1373 johnpye 928 /*
1374     Now we have the derivatives wrt each alg/diff variable in the
1375     present equation. variables[] points into the varlist. need
1376     a mapping from the varlist to the y and ydot lists.
1377     */
1378    
1379 johnpye 909 Jv_i = 0;
1380     for(j=0; j < count; ++j){
1381 johnpye 946 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
1382 johnpye 949 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1383 johnpye 909 if(varname){
1384 johnpye 949 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1385     ASC_FREE(varname);
1386 johnpye 909 }else{
1387 johnpye 912 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1388 johnpye 909 }
1389 johnpye 949 */
1390 johnpye 979
1391     /* we don't calculate derivatives wrt indep var */
1392     asc_assert(variables[j]>=0);
1393     if(enginedata->varlist[variables[j]] == blsys->x) continue;
1394    
1395 johnpye 951 var_yindex = blsys->y_id[variables[j]];
1396 johnpye 979 #ifdef JEX_DEBUG
1397 johnpye 977 CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %ld",j,variables[j],var_yindex);
1398 johnpye 979 #endif
1399    
1400 johnpye 977 ASC_ASSERT_RANGE(-var_yindex-1, -NV_LENGTH_S(v),NV_LENGTH_S(v));
1401 johnpye 928
1402     if(var_yindex >= 0){
1403 johnpye 951 #ifdef JEX_DEBUG
1404     asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
1405 johnpye 976 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%ld] = %f, v[%ld] = %f)\n", i
1406 johnpye 951 , derivatives[j] * NV_Ith_S(v,var_yindex)
1407     , i, var_yindex, derivatives[j]
1408     , var_yindex, NV_Ith_S(v,var_yindex)
1409     );
1410     #endif
1411 johnpye 928 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
1412     }else{
1413 johnpye 951 #ifdef JEX_DEBUG
1414 johnpye 976 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%ld] = %f, v[%ld] = %f)\n", i
1415 johnpye 951 , derivatives[j] * NV_Ith_S(v,-var_yindex-1)
1416     , i, -var_yindex-1, derivatives[j]
1417     , -var_yindex-1, NV_Ith_S(v,-var_yindex-1)
1418     );
1419     #endif
1420     asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
1421 johnpye 980 Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
1422 johnpye 928 }
1423 johnpye 909 }
1424    
1425     NV_Ith_S(Jv,i) = Jv_i;
1426 johnpye 951 #ifdef JEX_DEBUG
1427 johnpye 979 CONSOLE_DEBUG("rel = %p",*relptr);
1428 johnpye 951 relname = rel_make_name(blsys->system, *relptr);
1429     CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1430 johnpye 979 //ASC_FREE(relname);
1431 johnpye 980 return 1;
1432 johnpye 951 #endif
1433 johnpye 909 }
1434 johnpye 1142 #ifdef ASC_SIGNAL_TRAPS
1435 johnpye 909 }else{
1436 johnpye 951 relname = rel_make_name(blsys->system, *relptr);
1437     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1438     ASC_FREE(relname);
1439     is_error = 1;
1440 johnpye 909 }
1441 johnpye 968 Asc_SignalHandlerPopDefault(SIGFPE);
1442 johnpye 1142 #endif
1443 johnpye 909
1444 johnpye 951 if(is_error){
1445     CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1446     return 1;
1447     }
1448     return 0;
1449 johnpye 909 }
1450    
1451 johnpye 1012 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
1452     int integrator_ida_sjex(long int Neq, realtype tt
1453     , N_Vector yy, N_Vector yp, N_Vector rr
1454     , realtype c_j, void *jac_data, mtx_matrix_t Jac
1455     , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
1456 johnpye 1013 ){
1457     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1458     return -1;
1459     }
1460 johnpye 1012
1461 johnpye 669 /*----------------------------------------------
1462 johnpye 993 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1463 johnpye 992 */
1464    
1465 johnpye 993 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys){
1466     IntegratorIdaData * enginedata =blsys->enginedata;
1467     IntegratorIdaPrecDataJacobi *precdata;
1468     precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1469 johnpye 992
1470 johnpye 993 asc_assert(blsys->n_y);
1471     precdata->PIii = N_VNew_Serial(blsys->n_y);
1472    
1473     enginedata->pfree = &integrator_ida_pfree_jacobi;
1474     enginedata->precdata = precdata;
1475     CONSOLE_DEBUG("Allocated memory for Jacobi preconditioner");
1476     }
1477    
1478     void integrator_ida_pfree_jacobi(IntegratorIdaData *enginedata){
1479     if(enginedata->precdata){
1480     IntegratorIdaPrecDataJacobi *precdata = (IntegratorIdaPrecDataJacobi *)enginedata->precdata;
1481     N_VDestroy_Serial(precdata->PIii);
1482    
1483     ASC_FREE(precdata);
1484     enginedata->precdata = NULL;
1485     CONSOLE_DEBUG("Freed memory for Jacobi preconditioner");
1486     }
1487     enginedata->pfree = NULL;
1488     }
1489    
1490 johnpye 992 /**
1491 johnpye 993 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1492 johnpye 992
1493     'setup' function.
1494     */
1495 johnpye 993 int integrator_ida_psetup_jacobi(realtype tt,
1496 johnpye 992 N_Vector yy, N_Vector yp, N_Vector rr,
1497     realtype c_j, void *p_data,
1498     N_Vector tmp1, N_Vector tmp2,
1499     N_Vector tmp3
1500     ){
1501 johnpye 993 int i, j;
1502 johnpye 992 IntegratorSystem *blsys;
1503 johnpye 993 IntegratorIdaData *enginedata;
1504     IntegratorIdaPrecDataJacobi *precdata;
1505     struct rel_relation **relptr;
1506    
1507 johnpye 992 blsys = (IntegratorSystem *)p_data;
1508 johnpye 993 enginedata = blsys->enginedata;
1509     precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1510     double *derivatives;
1511     int *variables;
1512     int count, status;
1513     char *relname;
1514     int var_yindex;
1515 johnpye 992
1516 johnpye 993 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1517 johnpye 992
1518 johnpye 993 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
1519     derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1520    
1521     /**
1522     @TODO FIXME here we are using the very inefficient and contorted approach
1523     of calculating the whole jacobian, then extracting just the diagonal elements.
1524     */
1525    
1526     for(i=0, relptr = enginedata->rellist;
1527     i< enginedata->nrels && relptr != NULL;
1528     ++i, ++relptr
1529     ){
1530    
1531     /* get derivatives for this particular relation */
1532 johnpye 994 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1533 johnpye 993 if(status){
1534     relname = rel_make_name(blsys->system, *relptr);
1535     CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1536     ASC_FREE(relname);
1537     break;
1538     }
1539     /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1540     /* find the diagonal elements */
1541     for(j=0; j<count; ++j){
1542     if(variables[j]==i){
1543     var_yindex = blsys->y_id[variables[j]];
1544     if(var_yindex >= 0){
1545     NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1546     }else{
1547     NV_Ith_S(precdata->PIii, i) = 1./(c_j * derivatives[j]);
1548     }
1549     }
1550     }
1551     #ifdef PREC_DEBUG
1552     CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1553     #endif
1554 johnpye 992 }
1555 johnpye 993
1556     if(status){
1557     CONSOLE_DEBUG("Error found when evaluating derivatives");
1558     return 1; /* recoverable */
1559     }
1560    
1561     integrator_ida_write_incidence(blsys);
1562    
1563     ASC_FREE(variables);
1564     ASC_FREE(derivatives);
1565    
1566 johnpye 992 return 0;
1567     };
1568    
1569     /**
1570 johnpye 993 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1571 johnpye 992
1572     'solve' function.
1573     */
1574 johnpye 993 int integrator_ida_psolve_jacobi(realtype tt,
1575 johnpye 992 N_Vector yy, N_Vector yp, N_Vector rr,
1576     N_Vector rvec, N_Vector zvec,
1577     realtype c_j, realtype delta, void *p_data,
1578     N_Vector tmp
1579     ){
1580     IntegratorSystem *blsys;
1581     IntegratorIdaData *data;
1582 johnpye 993 IntegratorIdaPrecDataJacobi *precdata;
1583 johnpye 992 blsys = (IntegratorSystem *)p_data;
1584     data = blsys->enginedata;
1585 johnpye 993 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1586 johnpye 992
1587 johnpye 993 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1588     N_VProd(precdata->PIii, rvec, zvec);
1589 johnpye 992 return 0;
1590     };
1591    
1592     /*----------------------------------------------
1593 johnpye 993 STATS
1594     */
1595    
1596     /**
1597     A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1598     status in a struct instead of separately.
1599    
1600     @return IDA_SUCCESS on sucess.
1601     */
1602     int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s){
1603     return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1604     ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1605     ,&s->hlast, &s->hcur, &s->tcur
1606     );
1607     }
1608    
1609     /**
1610     This routine just outputs the stats to the CONSOLE_DEBUG routine.
1611    
1612     @TODO provide a GUI way of stats reporting from IDA.
1613     */
1614     void integrator_ida_write_stats(IntegratorIdaStats *stats){
1615     # define SL(N) CONSOLE_DEBUG("%s = %ld",#N,stats->N)
1616     # define SI(N) CONSOLE_DEBUG("%s = %d",#N,stats->N)
1617     # define SR(N) CONSOLE_DEBUG("%s = %f",#N,stats->N)
1618     SL(nsteps); SL(nrevals); SL(nlinsetups); SL(netfails);
1619     SI(qlast); SI(qcur);
1620     SR(hinused); SR(hlast); SR(hcur); SR(tcur);
1621     # undef SL
1622     # undef SI
1623     # undef SR
1624     }
1625    
1626     /*------------------------------------------------------------------------------
1627 johnpye 994 JACOBIAN / INCIDENCE MATRIX OUTPUT
1628 johnpye 993 */
1629    
1630 johnpye 994 enum integrator_ida_write_jac_enum{
1631     II_WRITE_Y
1632     , II_WRITE_YDOT
1633     };
1634    
1635 johnpye 993 /**
1636 johnpye 994 @TODO COMPLETE THIS...
1637     */
1638     void integrator_ida_write_jacobian(IntegratorSystem *blsys, realtype c_j, FILE *f, enum integrator_ida_write_jac_enum type){
1639     IntegratorIdaData *enginedata;
1640     MM_typecode matcode;
1641     int nnz, rhomax;
1642     double *derivatives;
1643     int *variables;
1644     struct rel_relation **relptr;
1645     int i, j, status, count, var_yindex;
1646     char *relname;
1647    
1648     var_filter_t vfiltery = {
1649     VAR_SVAR | VAR_FIXED | VAR_DERIV
1650     , VAR_SVAR
1651     };
1652     var_filter_t vfilteryd = {
1653     VAR_SVAR | VAR_FIXED | VAR_DERIV
1654     , VAR_SVAR | VAR_DERIV
1655     };
1656    
1657     enginedata = (IntegratorIdaData *)blsys->enginedata;
1658    
1659     /* number of non-zeros for all the non-FIXED solver_vars,
1660     in all the active included equality relations.
1661     */
1662     nnz = relman_jacobian_count(enginedata->rellist, enginedata->nrels
1663     , &enginedata->vfilter, &enginedata->rfilter
1664     , &rhomax
1665     );
1666    
1667     /* we must have found the same number of relations */
1668     asc_assert(rhomax == enginedata->nrels);
1669    
1670     /* output the mmio file header, now that we know our size*/
1671     /* note that we are asserting that our problem is square */
1672     mm_initialize_typecode(&matcode);
1673     mm_set_matrix(&matcode);
1674     mm_set_coordinate(&matcode);
1675     mm_set_real(&matcode);
1676     mm_write_banner(f, matcode);
1677     mm_write_mtx_crd_size(f, enginedata->nrels, enginedata->nrels, nnz);
1678    
1679     variables = ASC_NEW_ARRAY(int, blsys->n_y * 2);
1680     derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1681    
1682     CONSOLE_DEBUG("Writing sparse Jacobian to file...");
1683    
1684     for(i=0, relptr = enginedata->rellist;
1685     i< enginedata->nrels && relptr != NULL;
1686     ++i, ++relptr
1687     ){
1688     relname = rel_make_name(blsys->system, *relptr);
1689    
1690     /* get derivatives of y */
1691     status = relman_diff2(*relptr, &vfiltery, derivatives, variables, &count, enginedata->safeeval);
1692     if(status){
1693     CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1694     ASC_FREE(relname);
1695     break;
1696     }
1697    
1698     /* get derivatives of y */
1699     status = relman_diff2(*relptr, &vfilteryd, derivatives, variables, &count, enginedata->safeeval);
1700     if(status){
1701     CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1702     ASC_FREE(relname);
1703     break;
1704     }
1705    
1706     for(j=0; j<count; ++j){
1707     var_yindex = blsys->y_id[variables[j]];
1708     if(var_yindex >= 0 && type == II_WRITE_Y){
1709     fprintf(f, "%d %d %10.3g\n", i, var_yindex, derivatives[j]);
1710     }else if(var_yindex < 0 && type == II_WRITE_YDOT){
1711     fprintf(f, "%d %d %10.3g\n", i, -var_yindex-1, derivatives[j]);
1712     }
1713     }
1714     }
1715     ASC_FREE(variables);
1716     ASC_FREE(derivatives);
1717     }
1718    
1719     /**
1720 johnpye 993 This routine outputs matrix structure in a crude text format, for the sake
1721     of debugging.
1722     */
1723     void integrator_ida_write_incidence(IntegratorSystem *blsys){
1724     int i, j;
1725     struct rel_relation **relptr;
1726     IntegratorIdaData *enginedata = blsys->enginedata;
1727     double *derivatives;
1728     int *variables;
1729     int count, status;
1730     char *relname;
1731     int var_yindex;
1732    
1733     if(enginedata->nrels > 100){
1734     CONSOLE_DEBUG("Ignoring call (matrix size too big = %d)",enginedata->nrels);
1735     return;
1736     }
1737    
1738     variables = ASC_NEW_ARRAY(int, blsys->n_y * 2);
1739     derivatives = ASC_NEW_ARRAY(double, blsys->n_y * 2);
1740    
1741     CONSOLE_DEBUG("Outputting incidence information to console...");
1742    
1743     for(i=0, relptr = enginedata->rellist;
1744     i< enginedata->nrels && relptr != NULL;
1745     ++i, ++relptr
1746     ){
1747     relname = rel_make_name(blsys->system, *relptr);
1748    
1749     /* get derivatives for this particular relation */
1750 johnpye 994 status = relman_diff2(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1751 johnpye 993 if(status){
1752     CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1753     ASC_FREE(relname);
1754     break;
1755     }
1756    
1757     fprintf(stderr,"%3d:%-15s:",i,relname);
1758     ASC_FREE(relname);
1759    
1760     for(j=0; j<count; ++j){
1761     var_yindex = blsys->y_id[variables[j]];
1762     if(var_yindex >= 0){
1763     fprintf(stderr," %d:y[%d]",variables[j],var_yindex);
1764     }else{
1765     fprintf(stderr," %d:ydot[%d]",variables[j],-var_yindex-1);
1766     }
1767     }
1768     fprintf(stderr,"\n");
1769     }
1770     ASC_FREE(variables);
1771     ASC_FREE(derivatives);
1772     }
1773    
1774     /*----------------------------------------------
1775 johnpye 669 ERROR REPORTING
1776 johnpye 782 */
1777 johnpye 669 /**
1778     Error message reporter function to be passed to IDA. All error messages
1779     will trigger a call to this function, so we should find everything
1780     appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1781     panel (in the case of PyGTK).
1782     */
1783     void integrator_ida_error(int error_code
1784     , const char *module, const char *function
1785     , char *msg, void *eh_data
1786     ){
1787     IntegratorSystem *blsys;
1788     error_severity_t sev;
1789    
1790     /* cast back the IntegratorSystem, just in case we need it */
1791     blsys = (IntegratorSystem *)eh_data;
1792    
1793     /* severity depends on the sign of the error_code value */
1794     if(error_code <= 0){
1795     sev = ASC_PROG_ERR;
1796     }else{
1797     sev = ASC_PROG_WARNING;
1798     }
1799    
1800     /* use our all-purpose error reporting to get stuff back to the GUI */
1801     error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1802     }

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