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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1223 - (show annotations) (download) (as text)
Wed Jan 24 14:25:15 2007 UTC (15 years, 5 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 /* ASCEND modelling environment
2 Copyright (C) 2006-2007 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
21 as part of the GPL-licensed SUNDIALS solver package from LLNL.
22
23 IDA provides the non-linear parts, as well as a number of pluggable linear
24 solvers: dense, banded and krylov types.
25
26 We also implement here an EXPERIMENTAL direct sparse linear solver for IDA
27 using the ASCEND linsolqr routines.
28
29 @see http://www.llnl.gov/casc/sundials/
30 *//*
31 by John Pye, May 2006
32 */
33
34 #include "ida.h"
35
36 #include <signal.h>
37
38 /* SUNDIALS includes */
39 #ifdef ASC_WITH_IDA
40 # include <sundials/sundials_config.h>
41 # include <sundials/sundials_dense.h>
42 # include <ida/ida.h>
43 # include <nvector/nvector_serial.h>
44 # include <ida/ida_spgmr.h>
45 # include <ida/ida_spbcgs.h>
46 # include <ida/ida_sptfqmr.h>
47 # include <ida/ida_dense.h>
48 # ifndef IDA_SUCCESS
49 # error "Failed to include SUNDIALS IDA header file"
50 # endif
51 #else
52 # error "If you're building this file, you should have ASC_WITH_IDA"
53 #endif
54
55 #ifdef ASC_WITH_MMIO
56 # include <mmio.h>
57 #endif
58
59 #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
65 #include <solver/slv_client.h>
66 #include <solver/block.h>
67 #include <solver/slvDOF.h>
68 #include <solver/relman.h>
69 #ifdef ASC_IDA_NEW_ANALYSE
70 # include <solver/diffvars.h>
71 # include <solver/slv_stdcalls.h>
72 #endif
73
74 #include "idalinear.h"
75
76 /*
77 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
78 */
79 #ifndef SUNDIALS_VERSION_MINOR
80 # ifdef __GNUC__
81 # warning "GUESSING SUNDIALS VERSION 2.2"
82 # endif
83 # define SUNDIALS_VERSION_MINOR 2
84 #endif
85 #ifndef SUNDIALS_VERSION_MAJOR
86 # define SUNDIALS_VERSION_MAJOR 2
87 #endif
88
89 /* 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 /* #define FEX_DEBUG */
95 #define JEX_DEBUG
96 #define DJEX_DEBUG
97 #define SOLVE_DEBUG
98 #define STATS_DEBUG
99 #define PREC_DEBUG
100 /* #define ANALYSE_DEBUG */
101 /**
102 Everthing that the outside world needs to know about IDA
103 */
104 const IntegratorInternals integrator_ida_internals = {
105 integrator_ida_create
106 ,integrator_ida_params_default
107 #ifdef ASC_IDA_NEW_ANALYSE
108 ,integrator_ida_analyse
109 #else
110 ,integrator_analyse_dae /* note, this routine is back in integrator.c */
111 #endif
112 ,integrator_ida_solve
113 ,NULL /* writematrixfn */
114 ,integrator_ida_debug
115 ,integrator_ida_free
116 ,INTEG_IDA
117 ,"IDA"
118 };
119
120 /*-------------------------------------------------------------
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 /**
132 Struct containing any stuff that IDA needs that doesn't fit into the
133 common IntegratorSystem struct.
134 */
135 typedef struct IntegratorIdaDataStruct{
136 struct rel_relation **rellist; /**< NULL terminated list of rels */
137 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
138 int nrels;
139 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
140 var_filter_t vfilter;
141 rel_filter_t rfilter; /**< Used to filter relations from rellist (@TODO needs work) */
142 void *precdata; /**< For use by the preconditioner */
143 IntegratorIdaPrecFreeFn *pfree; /**< Store instructions here on how to free precdata */
144 } IntegratorIdaData;
145
146 typedef struct IntegratorIdaPrecDJStruct{
147 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 */
157 typedef struct IntegratorIdaPrecStruct{
158 IntegratorIdaPrecCreateFn *pcreate;
159 IDASpilsPrecSetupFn psetup;
160 IDASpilsPrecSolveFn psolve;
161 } IntegratorIdaPrec;
162
163 /* residual function forward declaration */
164 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
165
166 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 /* 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 /* dense jacobian evaluation for IDADense dense direct linear solver */
178 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 /* sparse jacobian evaluation for ASCEND's sparse direct solver */
185 IntegratorSparseJacFn integrator_ida_sjex;
186
187 typedef struct IntegratorIdaStatsStruct{
188 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
201 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 /*------
205 Jacobi preconditioner -- experimental
206 */
207
208 int integrator_ida_psetup_jacobi(realtype tt,
209 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 int integrator_ida_psolve_jacobi(realtype tt,
216 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 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys);
223
224 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 /*-------------------------------------------------------------
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 CONSOLE_DEBUG("enginedata = %p",enginedata);
240 enginedata->rellist = NULL;
241 enginedata->safeeval = 0;
242 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED;
243 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;
244 enginedata->pfree = NULL;
245
246 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
247 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
248
249 blsys->enginedata = (void *)enginedata;
250
251 integrator_ida_params_default(blsys);
252 }
253
254 void integrator_ida_free(void *enginedata){
255 CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);
256 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
257 asc_assert(d);
258 if(d->pfree){
259 CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);
260 /* free the preconditioner data, whatever it happens to be */
261 (d->pfree)(enginedata);
262 }
263 /* note, we don't own the rellist, so don't need to free it */
264 CONSOLE_DEBUG("Now destroying the enginedata");
265 ASC_FREE(d);
266 CONSOLE_DEBUG("enginedata freed");
267 }
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 PARAMETERS FOR IDA
280 */
281
282 enum ida_parameters{
283 IDA_PARAM_LINSOLVER
284 ,IDA_PARAM_MAXL
285 ,IDA_PARAM_AUTODIFF
286 ,IDA_PARAM_CALCIC
287 ,IDA_PARAM_SAFEEVAL
288 ,IDA_PARAM_RTOL
289 ,IDA_PARAM_ATOL
290 ,IDA_PARAM_ATOLVECT
291 ,IDA_PARAM_GSMODIFIED
292 ,IDA_PARAM_MAXNCF
293 ,IDA_PARAM_PREC
294 ,IDA_PARAMS_SIZE
295 };
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 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 @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 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 }, "YA_YDP"}, (char *[]){"Y", "YA_YDP", "NONE",NULL}
343 );
344
345 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 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 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
368 }, 1e-5, 1e-15, 1.0e15 }
369 );
370
371 slv_param_real(p,IDA_PARAM_RTOL
372 ,(SlvParameterInitReal){{"rtol"
373 ,"Scalar relative error tolerance",1
374 ,"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 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
378 }, 1e-4, 0, 1 }
379 );
380
381 slv_param_char(p,IDA_PARAM_LINSOLVER
382 ,(SlvParameterInitChar){{"linsolver"
383 ,"Linear solver",1
384 ,"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 }, "SPGMR"}, (char *[]){"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
390 );
391
392 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 slv_param_bool(p,IDA_PARAM_GSMODIFIED
403 ,(SlvParameterInitBool){{"gsmodified"
404 ,"Gram-Schmidt Orthogonalisation Scheme", 2
405 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
406 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
407 }, TRUE}
408 );
409
410 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 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 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 /*------------------------------------------------------------------------------
433 ANALYSIS ROUTINE (new implementation)
434 */
435
436 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 var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
444 ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
445 for(i=0;i<nv;++i){
446 v=vlist[i];
447 if(!var_apply_filter(v,&vf))continue;
448 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 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 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 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 static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys){
529 int i;
530 const SolverDiffVarCollection *diffvars;
531 SolverDiffVarSequence seq;
532 char *varname;
533
534 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 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 void integrator_dae_show_var(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
578
579 #ifdef ASC_IDA_NEW_ANALYSE
580 /**
581 Perform additional problem analysis to prepare problem for integration with
582 IDA.
583
584 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
588 We can also assume that the independent variable has been found.
589
590 See mailing list, ~Jan 2007.
591
592 Note, the stuff for identifying the static and output sub-problems should
593 be part of integrator.c, not this file. We will assume this is handled
594
595 @return 0 on success
596 @see integrator_analyse
597 */
598 int integrator_ida_analyse(IntegratorSystem *sys){
599 /* struct var_variable **solversvars;
600 unsigned long nsolversvars;
601 struct rel_relation **solversrels;
602 unsigned long nsolversrels;*/
603 const SolverDiffVarCollection *diffvars;
604 SolverDiffVarSequence seq;
605 long i, j, n_y, n_ydot, n_skipped_diff, n_skipped_alg, n_skipped_deriv;
606 /* int res; */
607
608 struct var_variable *v;
609 char *varname;
610
611 CONSOLE_DEBUG("NEW integrator_ida_analyse------------------>");
612
613 asc_assert(sys->engine==INTEG_IDA);
614
615 #if 0
616 /* partition into static, dynamic and output problems */
617 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 #endif
623
624 /* 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 diffvars = slv_get_diffvars(sys->system);
631 if(diffvars==NULL){
632 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
633 return 6;
634 }
635
636 CONSOLE_DEBUG("Got %ld chains, maxorder = %ld",diffvars->nseqs,diffvars->maxorder);
637
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
647 /* set up the dynamic problem */
648 CONSOLE_DEBUG("Setting up the dynamic problem");
649
650 asc_assert(sys->y == NULL);
651
652 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 n_y = 0; n_ydot = 0;
655 n_skipped_alg = 0; n_skipped_diff = 0; n_skipped_deriv = 0;
656
657 /* add the variables from the derivative chains */
658 CONSOLE_DEBUG("Counting vars (%ld)",n_y);
659 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 n_skipped_alg++;
669 continue;
670 }else if(!var_incident(v)){
671 CONSOLE_DEBUG("'%s' is not incident",varname);
672 ASC_FREE(varname);
673 n_skipped_alg++;
674 continue;
675 }
676 CONSOLE_DEBUG("'%s' is algebraic (%ld)",varname,n_y);
677 }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 n_skipped_diff++;
684 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 CONSOLE_DEBUG("Derivative '%s' SET INACTIVE",varname);
691 ASC_FREE(varname);
692 n_skipped_deriv++;
693 }
694 continue;
695 }else if(var_fixed(seq.vars[1])){
696 /* diff var with fixed derivative */
697 CONSOLE_DEBUG("Derivative of var '%s' is fixed; CONVERTING TO ALGEBRAIC",varname);
698 n_skipped_deriv++;
699 }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 CONSOLE_DEBUG("'%s' is differential (%ld)",varname,n_y);
705 ASC_FREE(varname);
706 v = seq.vars[1];
707 sys->ydot[n_y] = v;
708 n_ydot++;
709 varname = var_make_name(sys->system,v);
710 CONSOLE_DEBUG("'%s' is derivative (%ld)",varname,n_y);
711 ASC_FREE(varname);
712 n_y++;
713 continue;
714 }
715 }
716 /* fall through: v is algebraic */
717 ASC_FREE(varname);
718 asc_assert(var_active(v));
719 sys->y[n_y] = v;
720 sys->ydot[n_y] = NULL;
721 n_y++;
722 continue;
723 }
724 sys->n_y = n_y;
725
726 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 sys->y_id = ASC_NEW_ARRAY_CLEAR(int, slv_get_num_solvers_vars(sys->system));
731 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 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 slv_block_partition(sys->system);
750 res = integrator_ida_block_check(sys);
751 if(res)return 100 + res;
752
753 /* ...(2) FREE the derivatives, FIX the diffvars */
754 CONSOLE_DEBUG("Checking system with differential variables fixed...");
755 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 /* ...(3) restore as it was, FREE the diffvars */
766 for(i=0;i<n_y;++i){
767 if(sys->ydot[i]){
768 var_set_fixed(sys->y[i],0);
769 }
770 }
771 #endif
772
773 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
779 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 /* 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 return 3;
797 }else if(diffvars->nindep<1){
798 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
799 return 4;
800 }else if(diffvars->indep[0]!=sys->x){
801 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
802 return 5;
803 }
804
805 /* get the observations */
806 /** @TODO this should use a 'problem provider' API instead of static data
807 from the system object */
808 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 }
816
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 /* ERROR_REPORTER_HERE(ASC_PROG_ERR,"Implementation incomplete");
824 return -1; */
825
826 return 0;
827 }
828 #endif
829
830 /*-------------------------------------------------------------
831 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
832 */
833
834 /*static double div1(double a, double b){
835 return a/b;
836 }*/
837
838 typedef int IdaFlagFn(void *,int *);
839 typedef char *IdaFlagNameFn(int);
840
841 /* return 0 on success */
842 int integrator_ida_solve(
843 IntegratorSystem *blsys
844 , unsigned long start_index
845 , unsigned long finish_index
846 ){
847 void *ida_mem;
848 int size, flag, flag1, t_index;
849 realtype t0, reltol, abstol, t, tret, tout1;
850 N_Vector y0, yp0, abstolvect, ypret, yret, id;
851 IntegratorIdaData *enginedata;
852 char *linsolver;
853 int maxl;
854 IdaFlagFn *flagfn;
855 IdaFlagNameFn *flagnamefn;
856 const char *flagfntype;
857 char *pname = NULL;
858 char *varname;
859 int i;
860 const IntegratorIdaPrec *prec = NULL;
861 int icopt; /* initial conditions strategy */
862
863 CONSOLE_DEBUG("STARTING IDA...");
864
865 enginedata = integrator_ida_enginedata(blsys);
866
867 enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);
868 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 enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);
874
875 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
876 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
877 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 return 1; /* failure */
882 }
883
884 /* retrieve initial values from the system */
885
886 /** @TODO fix this, the starting time != first sample */
887 t0 = integrator_get_t(blsys);
888 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
889
890 CONSOLE_DEBUG("RETRIEVING y0");
891
892 y0 = N_VNew_Serial(size);
893 integrator_get_y(blsys,NV_DATA_S(y0));
894
895 #ifdef SOLVE_DEBUG
896 CONSOLE_DEBUG("RETRIEVING yp0");
897 #endif
898
899 yp0 = N_VNew_Serial(size);
900 integrator_get_ydot(blsys,NV_DATA_S(yp0));
901
902 #ifdef SOLVE_DEBUG
903 N_VPrint_Serial(yp0);
904 CONSOLE_DEBUG("yp0 is at %p",&yp0);
905 #endif
906
907 /* create IDA object */
908 ida_mem = IDACreate();
909
910 /* relative error tolerance */
911 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
912 CONSOLE_DEBUG("rtol = %8.2e",reltol);
913
914 /* allocate internal memory */
915 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
923 N_VDestroy_Serial(abstolvect);
924 }else{
925 /* scalar absolute tolerance (one value for all) */
926 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
927 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 if(flag==IDA_MEM_NULL){
932 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
933 return 2;
934 }else if(flag==IDA_MEM_FAIL){
935 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
936 return 3;
937 }else if(flag==IDA_ILL_INPUT){
938 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
939 return 4;
940 }/* else success */
941
942 /* set optional inputs... */
943 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
944 IDASetRdata(ida_mem, (void *)blsys);
945 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
946 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
947 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
948 if(integrator_get_minstep(blsys)>0){
949 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
950 }
951
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 /* there's no capability for setting *minimum* step size in IDA */
956
957
958 /* attach linear solver module, using the default value of maxl */
959 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
960 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
961 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 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
972 flag = IDADense(ida_mem, size);
973 switch(flag){
974 case IDADENSE_SUCCESS: break;
975 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 }
980
981 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 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 6;
987 }
988 }else{
989 CONSOLE_DEBUG("USING NUMERICAL DIFF");
990 }
991
992 flagfntype = "IDADENSE";
993 flagfn = &IDADenseGetLastFlag;
994 flagnamefn = &IDADenseGetReturnFlagName;
995 }else{
996 /* remaining methods are all SPILS */
997 CONSOLE_DEBUG("IDA SPILS");
998
999 maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
1000 CONSOLE_DEBUG("maxl = %d",maxl);
1001
1002 /* what preconditioner? */
1003 pname = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_PREC);
1004 if(strcmp(pname,"NONE")==0){
1005 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 return 7;
1011 }
1012
1013 /* which SPILS linear solver? */
1014 if(strcmp(linsolver,"SPGMR")==0){
1015 CONSOLE_DEBUG("IDA SPGMR");
1016 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
1017 }else if(strcmp(linsolver,"SPBCG")==0){
1018 CONSOLE_DEBUG("IDA SPBCG");
1019 flag = IDASpbcg(ida_mem, maxl);
1020 }else if(strcmp(linsolver,"SPTFQMR")==0){
1021 CONSOLE_DEBUG("IDA SPTFQMR");
1022 flag = IDASptfqmr(ida_mem,maxl);
1023 }else{
1024 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
1025 return 8;
1026 }
1027
1028 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
1037 flagfntype = "IDASPILS";
1038 flagfn = &IDASpilsGetLastFlag;
1039 flagnamefn = &IDASpilsGetReturnFlagName;
1040
1041 if(flag==IDASPILS_MEM_NULL){
1042 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
1043 return 9;
1044 }else if(flag==IDASPILS_MEM_FAIL){
1045 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
1046 return 9;
1047 }/* else success */
1048
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 return 10;
1056 }else if(flag==IDASPILS_LMEM_NULL){
1057 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
1058 return 10;
1059 }/* else success */
1060 }else{
1061 CONSOLE_DEBUG("USING NUMERICAL DIFF");
1062 }
1063
1064 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 return 11;
1072 }
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 return 11;
1079 }
1080 }
1081 }
1082 }
1083
1084 /* set linear solver optional inputs...
1085 ...nothing here at the moment...
1086 */
1087
1088 /* calculate initial conditions */
1089 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
1116 if(icopt){
1117 blsys->currentstep=0;
1118 t_index=start_index + 1;
1119 tout1 = samplelist_get(blsys->samples, t_index);
1120
1121 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
1122
1123 #ifdef ASC_SIGNAL_TRAPS
1124 /* catch SIGFPE if desired to */
1125 if(enginedata->safeeval){
1126 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
1127 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1128 }else{
1129 #ifdef FEX_DEBUG
1130 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
1131 #endif
1132 Asc_SignalHandlerPushDefault(SIGFPE);
1133 }
1134 if (setjmp(g_fpe_env)==0) {
1135 #endif
1136
1137 //CONSOLE_DEBUG("Raising signal...");
1138 //CONSOLE_DEBUG("1/0 = %f", div1(1.0,0.0));
1139 //CONSOLE_DEBUG("Still here...");
1140
1141 /* correct initial values, given derivatives */
1142 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
1143 /* note the new API from version 2.3 and onwards */
1144 flag = IDACalcIC(ida_mem, icopt, tout1);
1145 # else
1146 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
1147 # endif
1148
1149 switch(flag){
1150 case IDA_SUCCESS:
1151 CONSOLE_DEBUG("Initial conditions solved OK");
1152 break;
1153
1154 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 return 12;
1163 }
1164 ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
1165 return 12;
1166
1167 default:
1168 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
1169 return 12;
1170 }
1171 #ifdef ASC_SIGNAL_TRAPS
1172 }else{
1173 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
1174 return 13;
1175 }
1176
1177 if(enginedata->safeeval){
1178 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1179 }else{
1180 CONSOLE_DEBUG("pop...");
1181 Asc_SignalHandlerPopDefault(SIGFPE);
1182 CONSOLE_DEBUG("...pop");
1183 }
1184 #endif
1185
1186 }
1187
1188 /* 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 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
1204 t = samplelist_get(blsys->samples, t_index);
1205 t0 = integrator_get_t(blsys);
1206 asc_assert(t > t0);
1207
1208 #ifdef SOLVE_DEBUG
1209 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, t);
1210 #endif
1211
1212 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 /* 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 }
1240 #endif
1241
1242 /* free solution memory */
1243 N_VDestroy_Serial(yret);
1244 N_VDestroy_Serial(ypret);
1245
1246 /* free solver memory */
1247 IDAFree(ida_mem);
1248
1249 if(flag < 0){
1250 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
1251 return 14;
1252 }
1253
1254 /* all done, success */
1255 return 0;
1256 }
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 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
1276 IntegratorSystem *blsys;
1277 IntegratorIdaData *enginedata;
1278 int i, calc_ok, is_error;
1279 struct rel_relation** relptr;
1280 double resid;
1281 char *relname;
1282 #ifdef FEX_DEBUG
1283 char *varname;
1284 char diffname[30];
1285 #endif
1286
1287 blsys = (IntegratorSystem *)res_data;
1288 enginedata = integrator_ida_enginedata(blsys);
1289
1290 #ifdef FEX_DEBUG
1291 /* fprintf(stderr,"\n\n"); */
1292 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
1293 #endif
1294
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 /**
1306 @TODO does this function (fex) do bounds checking already?
1307 */
1308
1309 /* evaluate each residual in the rellist */
1310 is_error = 0;
1311 relptr = enginedata->rellist;
1312
1313 #ifdef ASC_SIGNAL_TRAPS
1314 if(enginedata->safeeval){
1315 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
1316 }else{
1317 # ifdef FEX_DEBUG
1318 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
1319 # endif
1320 Asc_SignalHandlerPushDefault(SIGFPE);
1321 }
1322
1323 if (SETJMP(g_fpe_env)==0) {
1324 #endif
1325 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 relname = rel_make_name(blsys->system, *relptr);
1334 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1335 ASC_FREE(relname);
1336 /* presumable some output already made? */
1337 is_error = 1;
1338 }/*else{
1339 CONSOLE_DEBUG("Calc OK");
1340 }*/
1341 }
1342
1343 #ifdef ASC_SIGNAL_TRAPS
1344 }else{
1345 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 }
1350
1351 if(enginedata->safeeval){
1352 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
1353 }else{
1354 Asc_SignalHandlerPopDefault(SIGFPE);
1355 }
1356 #endif
1357
1358 #ifdef FEX_DEBUG
1359 /* output residuals to console */
1360 CONSOLE_DEBUG("RESIDUAL OUTPUT");
1361 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
1362 for(i=0; i<blsys->n_y; ++i){
1363 varname = var_make_name(blsys->system,blsys->y[i]);
1364 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
1365 if(blsys->ydot[i]){
1366 varname = var_make_name(blsys->system,blsys->ydot[i]);
1367 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
1368 }else{
1369 snprintf(diffname,99,"diff(%s)",varname);
1370 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
1371 }
1372 ASC_FREE(varname);
1373 relname = rel_make_name(blsys->system,enginedata->rellist[i]);
1374 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
1375 }
1376 #endif
1377
1378 if(is_error){
1379 return 1;
1380 }
1381
1382 #ifdef FEX_DEBUG
1383 CONSOLE_DEBUG("RESIDUAL OK");
1384 #endif
1385 return 0;
1386 }
1387
1388 /**
1389 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 char *relname;
1399 #ifdef DJEX_DEBUG
1400 struct var_variable **varlist;
1401 char *varname;
1402 #endif
1403 int status;
1404 struct rel_relation **relptr;
1405 int i;
1406 double *derivatives;
1407 struct var_variable **variables;
1408 int count, j;
1409
1410 blsys = (IntegratorSystem *)jac_data;
1411 enginedata = integrator_ida_enginedata(blsys);
1412
1413 /* 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 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1417
1418 /* 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 #ifdef DJEX_DEBUG
1424 varlist = slv_get_solvers_var_list(blsys->system);
1425
1426 /* print vars */
1427 for(i=0; i < blsys->n_y; ++i){
1428 varname = var_make_name(blsys->system, blsys->y[i]);
1429 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 ASC_FREE(varname);
1432 }
1433
1434 /* 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 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
1439 ASC_FREE(varname);
1440 }else{
1441 varname = var_make_name(blsys->system, blsys->y[i]);
1442 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
1443 ASC_FREE(varname);
1444 }
1445 }
1446
1447 /* print step size */
1448 CONSOLE_DEBUG("<c_j> = %g",c_j);
1449 #endif
1450
1451 /* 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 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1459
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 /* output what's going on here ... */
1468 #ifdef DJEX_DEBUG
1469 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 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 }else{
1478 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
1479 }
1480 ASC_FREE(varname);
1481 }
1482 fprintf(stderr,"\n");
1483 #endif
1484
1485 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
1486 for(j=0; j < count; ++j){
1487 if(!var_deriv(variables[j])){
1488 #ifdef DJEX_DEBUG1
1489 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 varname = var_make_name(blsys->system,variables[j]);
1494 CONSOLE_DEBUG("var name '%s'",varname);
1495 ASC_FREE(varname);
1496 asc_assert(var_sindex(variables[j]) >= 0);
1497 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
1498 #endif
1499 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
1500 }else{
1501 DENSE_ELEM(Jac,i,integrator_ida_diffindex(blsys,variables[j])) += derivatives[j] * c_j;
1502 }
1503 }
1504 }
1505
1506 #ifdef DJEX_DEBUG
1507 CONSOLE_DEBUG("PRINTING JAC");
1508 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 for(j=0; j < blsys->n_y; ++j){
1522 if(j)fprintf(stderr,"\t");
1523 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
1524 }
1525 fprintf(stderr,"\n");
1526 }
1527 #endif
1528
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 #ifdef DJEX_DEBUG
1535 CONSOLE_DEBUG("DJEX RETURNING 0");
1536 #endif
1537 return 0;
1538 }
1539
1540 /**
1541 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
1542
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 */
1557 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 char *relname;
1566 int status;
1567 double Jv_i;
1568
1569 struct var_variable **variables;
1570 double *derivatives;
1571 int count;
1572 struct var_variable **varlist;
1573 #ifdef JEX_DEBUG
1574
1575 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
1576 #endif
1577
1578 blsys = (IntegratorSystem *)jac_data;
1579 enginedata = integrator_ida_enginedata(blsys);
1580 varlist = slv_get_solvers_var_list(blsys->system);
1581
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 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
1589
1590 i = NV_LENGTH_S(yy) * 2;
1591 #ifdef JEX_DEBUG
1592 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1593 #endif
1594 variables = ASC_NEW_ARRAY(struct var_variable*, i);
1595 derivatives = ASC_NEW_ARRAY(double, i);
1596
1597 /* evaluate the derivatives... */
1598 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1599
1600 #ifdef ASC_SIGNAL_TRAPS
1601 Asc_SignalHandlerPushDefault(SIGFPE);
1602 if (SETJMP(g_fpe_env)==0) {
1603 #endif
1604 for(i=0, relptr = enginedata->rellist;
1605 i< enginedata->nrels && relptr != NULL;
1606 ++i, ++relptr
1607 ){
1608 /* get derivatives for this particular relation */
1609 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1610 #ifdef JEX_DEBUG
1611 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1612 #endif
1613
1614 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 break;
1620 }
1621
1622 /*
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 Jv_i = 0;
1629 for(j=0; j < count; ++j){
1630 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
1631 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1632 if(varname){
1633 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1634 ASC_FREE(varname);
1635 }else{
1636 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1637 }
1638 */
1639
1640 /* we don't calculate derivatives wrt indep var */
1641 asc_assert(variables[j]>=0);
1642 if(variables[j] == blsys->x) continue;
1643 #ifdef JEX_DEBUG
1644 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 #ifdef JEX_DEBUG
1649 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 );
1654 #endif
1655 asc_assert(blsys->ydot[DIFFINDEX]==variables[j]);
1656 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
1657 #undef DIFFINDEX
1658 }else{
1659 #define VARINDEX var_sindex(variables[j])
1660 #ifdef JEX_DEBUG
1661 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 );
1668 #endif
1669 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
1670 #undef VARINDEX
1671 }
1672 }
1673
1674 NV_Ith_S(Jv,i) = Jv_i;
1675 #ifdef JEX_DEBUG
1676 CONSOLE_DEBUG("rel = %p",*relptr);
1677 relname = rel_make_name(blsys->system, *relptr);
1678 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1679 ASC_FREE(relname);
1680 return 1;
1681 #endif
1682 }
1683 #ifdef ASC_SIGNAL_TRAPS
1684 }else{
1685 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 }
1690 Asc_SignalHandlerPopDefault(SIGFPE);
1691 #endif
1692
1693 if(is_error){
1694 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1695 return 1;
1696 }
1697 return 0;
1698 }
1699
1700 /* 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 ){
1706 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
1707 return -1;
1708 }
1709
1710 /*----------------------------------------------
1711 JACOBI PRECONDITIONER -- EXPERIMENTAL.
1712 */
1713
1714 void integrator_ida_pcreate_jacobi(IntegratorSystem *blsys){
1715 IntegratorIdaData * enginedata =blsys->enginedata;
1716 IntegratorIdaPrecDataJacobi *precdata;
1717 precdata = ASC_NEW(IntegratorIdaPrecDataJacobi);
1718
1719 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 /**
1740 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1741
1742 'setup' function.
1743 */
1744 int integrator_ida_psetup_jacobi(realtype tt,
1745 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 int i, j, res;
1751 IntegratorSystem *blsys;
1752 IntegratorIdaData *enginedata;
1753 IntegratorIdaPrecDataJacobi *precdata;
1754 struct rel_relation **relptr;
1755
1756 blsys = (IntegratorSystem *)p_data;
1757 enginedata = blsys->enginedata;
1758 precdata = (IntegratorIdaPrecDataJacobi *)(enginedata->precdata);
1759 double *derivatives;
1760 struct var_variable **variables;
1761 int count, status;
1762 char *relname;
1763
1764 CONSOLE_DEBUG("Setting up Jacobi preconditioner");
1765
1766 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1767 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 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1781 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 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 NV_Ith_S(precdata->PIii, i) = 1./derivatives[j];
1795 }
1796
1797 }
1798 }
1799 #ifdef PREC_DEBUG
1800 CONSOLE_DEBUG("PI[%d] = %f",i,NV_Ith_S(precdata->PIii,i));
1801 #endif
1802 }
1803
1804 if(status){
1805 CONSOLE_DEBUG("Error found when evaluating derivatives");
1806 res = 1; goto finish; /* recoverable */
1807 }
1808
1809 integrator_ida_write_incidence(blsys);
1810
1811 res = 0;
1812 finish:
1813 ASC_FREE(variables);
1814 ASC_FREE(derivatives);
1815 return res;
1816 };
1817
1818 /**
1819 EXPERIMENTAL. Jacobi preconditioner for use with IDA Krylov solvers
1820
1821 'solve' function.
1822 */
1823 int integrator_ida_psolve_jacobi(realtype tt,
1824 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 IntegratorIdaPrecDataJacobi *precdata;
1832 blsys = (IntegratorSystem *)p_data;
1833 data = blsys->enginedata;
1834 precdata = (IntegratorIdaPrecDataJacobi *)(data->precdata);
1835
1836 CONSOLE_DEBUG("Solving Jacobi preconditioner (c_j = %f)",c_j);
1837 N_VProd(precdata->PIii, rvec, zvec);
1838 return 0;
1839 };
1840
1841 /*----------------------------------------------
1842 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 OUTPUT OF INTERNALS: JACOBIAN / INCIDENCE MATRIX / DEBUG INFO
1877 */
1878
1879 enum integrator_ida_write_jac_enum{
1880 II_WRITE_Y
1881 , II_WRITE_YDOT
1882 };
1883
1884 /**
1885 @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 struct var_variable **variables;
1893 struct rel_relation **relptr;
1894 int i, j, status, count;
1895 char *relname;
1896
1897 var_filter_t vfilter = {
1898 VAR_SVAR | VAR_FIXED | VAR_DERIV
1899 ,VAR_SVAR | 0 | 0
1900 };
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 variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);
1924 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 status = relman_diff3(*relptr, &vfilter, derivatives, variables, &count, enginedata->safeeval);
1936 if(status){
1937 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
1938 ASC_FREE(relname);
1939 break;
1940 }
1941
1942 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 }
1948 }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 }
1955 }
1956 ASC_FREE(variables);
1957 ASC_FREE(derivatives);
1958 }
1959
1960 /**
1961 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 struct var_variable **variables;
1970 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 variables = ASC_NEW_ARRAY(struct var_variable *, blsys->n_y * 2);
1979 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 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1991 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 if(var_deriv(variables[j])){
2002 fprintf(stderr," %p:ydot[%d]",variables[j],integrator_ida_diffindex(blsys,variables[j]));
2003 }else{
2004 fprintf(stderr," %p:y[%d]",variables[j],var_sindex(variables[j]));
2005 }
2006 }
2007 fprintf(stderr,"\n");
2008 }
2009 ASC_FREE(variables);
2010 ASC_FREE(derivatives);
2011 }
2012
2013 /* @return 0 on success */
2014 int integrator_ida_debug(const IntegratorSystem *sys, FILE *fp){
2015 char *varname, *relname;
2016 struct var_variable **vlist, *var;
2017 struct rel_relation **rlist, *rel;
2018 long vlen, rlen;
2019 long i;
2020
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 fprintf(fp,".\n");
2038 ASC_FREE(varname);
2039 }
2040 }
2041
2042 fprintf(fp,"\n\nCORRESPONDENCE OF SOLVER VARS TO INTEGRATOR VARS\n\n");
2043 fprintf(fp,"sindex\t%-15s\ty \tydot \n","name");
2044 fprintf(fp,"------\t%-15s\t-----\t-----\n","----");
2045
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 for(i=0;i<vlen;i++){
2052 var = vlist[i];
2053
2054 varname = var_make_name(sys->system, var);
2055 fprintf(fp,"%ld\t%-15s\t",i,varname);
2056
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 }else if(!var_incident(var)){
2064 // not incident
2065 fprintf(fp,"(not incident)\n");
2066 }else{
2067 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 }else{
2072 fprintf(fp,"\t%d\t.\n",var_sindex(var));
2073 }
2074 }
2075 ASC_FREE(varname);
2076 }
2077
2078 /* 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
2092 /* and lets write block debug output */
2093 block_debug(sys->system, fp);
2094
2095 return 0; /* success */
2096 }
2097
2098 /*----------------------------------------------
2099 ERROR REPORTING
2100 */
2101 /**
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