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

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