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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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