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

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