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

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