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

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