/[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 992 - (show annotations) (download) (as text)
Thu Dec 21 13:37:45 2006 UTC (13 years, 2 months ago) by johnpye
File MIME type: text/x-csrc
File size: 35782 byte(s)
Added support for preconditioner in IDA module. Not quite sure what to put inside it yet though.
1 /* ASCEND modelling environment
2 Copyright (C) 2006 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 @see http://www.llnl.gov/casc/sundials/
23 *//*
24 by John Pye, May 2006
25 */
26
27 /*
28 Be careful with the following. This file requires both the 'ida.h' from
29 SUNDIALS as well as the 'ida.h' from ASCEND. Make sure that we're getting
30 both of these; if you get problems check your build tool for the paths being
31 passed to the C preprocessor.
32 */
33
34 /* standard includes */
35 #include <signal.h>
36
37 /* ASCEND includes */
38 #include "ida.h"
39 #include <utilities/error.h>
40 #include <utilities/ascConfig.h>
41 #include <utilities/ascSignal.h>
42 #include <utilities/ascPanic.h>
43 #include <compiler/instance_enum.h>
44 #include "var.h"
45 #include "rel.h"
46 #include "discrete.h"
47 #include "conditional.h"
48 #include "logrel.h"
49 #include "bnd.h"
50 #include "linsol.h"
51 #include "linsolqr.h"
52 #include "slv_common.h"
53 #include "slv_client.h"
54 #include "relman.h"
55
56 /* SUNDIALS includes */
57 #ifdef ASC_WITH_IDA
58 # include <sundials/sundials_config.h>
59 # include <sundials/sundials_dense.h>
60 # include <ida/ida.h>
61 # include <nvector/nvector_serial.h>
62 # include <ida/ida_spgmr.h>
63 # include <ida/ida_spbcgs.h>
64 # include <ida/ida_sptfqmr.h>
65 # include <ida/ida_dense.h>
66 # ifndef IDA_SUCCESS
67 # error "Failed to include SUNDIALS IDA header file"
68 # endif
69 #endif
70
71 /*
72 for cases where we don't have SUNDIALS_VERSION_MINOR defined, guess version 2.2
73 */
74 #ifndef SUNDIALS_VERSION_MINOR
75 # ifdef __GNUC__
76 # warning "GUESSING SUNDIALS VERSION 2.2"
77 # endif
78 # define SUNDIALS_VERSION_MINOR 2
79 #endif
80 #ifndef SUNDIALS_VERSION_MAJOR
81 # define SUNDIALS_VERSION_MAJOR 2
82 #endif
83
84 /* check that we've got what we expect now */
85 #ifndef ASC_IDA_H
86 # error "Failed to include ASCEND IDA header file"
87 #endif
88
89 /* #define FEX_DEBUG */
90 /* #define JEX_DEBUG */
91 /* #define SOLVE_DEBUG */
92
93 const IntegratorInternals integrator_ida_internals = {
94 integrator_ida_create
95 ,integrator_ida_params_default
96 ,integrator_analyse_dae /* note, this routine is back in integrator.c */
97 ,integrator_ida_solve
98 ,integrator_ida_free
99 ,INTEG_IDA
100 ,"IDA"
101 };
102
103 /**
104 Struct containing any stuff that IDA needs that doesn't fit into the
105 common IntegratorSystem struct.
106 */
107 typedef struct{
108 struct rel_relation **rellist; /**< NULL terminated list of rels */
109 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
110 struct bnd_boundary **bndlist; /**< NULL-terminated list of boundaries, for use in the root-finding code */
111 int nrels;
112 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
113 N_Vector pp; /**< Preconditioner values */
114 } IntegratorIdaData;
115
116 /*-------------------------------------------------------------
117 FORWARD DECLS
118 */
119 /* residual function forward declaration */
120 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
121
122 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
123 , N_Vector v, N_Vector Jv, realtype c_j
124 , void *jac_data, N_Vector tmp1, N_Vector tmp2
125 );
126
127 /* error handler forward declaration */
128 void integrator_ida_error(int error_code
129 , const char *module, const char *function
130 , char *msg, void *eh_data
131 );
132
133 int integrator_ida_djex(long int Neq, realtype tt
134 , N_Vector yy, N_Vector yp, N_Vector rr
135 , realtype c_j, void *jac_data, DenseMat Jac
136 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
137 );
138
139 int integrator_ida_psetup(realtype tt,
140 N_Vector yy, N_Vector yp, N_Vector rr,
141 realtype c_j, void *prec_data,
142 N_Vector tmp1, N_Vector tmp2,
143 N_Vector tmp3
144 );
145
146 int integrator_ida_psolve(realtype tt,
147 N_Vector yy, N_Vector yp, N_Vector rr,
148 N_Vector rvec, N_Vector zvec,
149 realtype c_j, realtype delta, void *prec_data,
150 N_Vector tmp
151 );
152
153 static const IDASpilsPrecSetupFn psetup1 = &integrator_ida_psetup;
154 static const IDASpilsPrecSolveFn psolve1 = &integrator_ida_psolve;
155
156 /*-------------------------------------------------------------
157 SETUP/TEARDOWN ROUTINES
158 */
159 void integrator_ida_create(IntegratorSystem *blsys){
160 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
161 IntegratorIdaData *enginedata;
162 enginedata = ASC_NEW(IntegratorIdaData);
163 enginedata->rellist = NULL;
164 enginedata->varlist = NULL;
165 enginedata->safeeval = 0;
166 enginedata->pp = NULL;
167 blsys->enginedata = (void *)enginedata;
168 integrator_ida_params_default(blsys);
169 }
170
171 void integrator_ida_free(void *enginedata){
172 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
173 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
174 if(d->pp)N_VDestroy_Serial(d->pp);
175 /* note, we don't own the rellist, so don't need to free it */
176 ASC_FREE(d);
177 }
178
179 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
180 IntegratorIdaData *d;
181 assert(blsys!=NULL);
182 assert(blsys->enginedata!=NULL);
183 assert(blsys->engine==INTEG_IDA);
184 d = ((IntegratorIdaData *)(blsys->enginedata));
185 return d;
186 }
187
188 /*-------------------------------------------------------------
189 PARAMETERS FOR IDA
190 */
191
192 enum ida_parameters{
193 IDA_PARAM_LINSOLVER
194 ,IDA_PARAM_MAXL
195 ,IDA_PARAM_AUTODIFF
196 ,IDA_PARAM_CALCIC
197 ,IDA_PARAM_SAFEEVAL
198 ,IDA_PARAM_RTOL
199 ,IDA_PARAM_ATOL
200 ,IDA_PARAM_ATOLVECT
201 ,IDA_PARAM_GSMODIFIED
202 ,IDA_PARAM_MAXNCF
203 ,IDA_PARAM_PREC
204 ,IDA_PARAMS_SIZE
205 };
206
207 /**
208 Here the full set of parameters is defined, along with upper/lower bounds,
209 etc. The values are stuck into the blsys->params structure.
210
211 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
212 list. Then add a slv_param_*(...) statement below to define the type, description and range
213 for the new parameter.
214
215 @return 0 on success
216 */
217 int integrator_ida_params_default(IntegratorSystem *blsys){
218 asc_assert(blsys!=NULL);
219 asc_assert(blsys->engine==INTEG_IDA);
220 slv_parameters_t *p;
221 p = &(blsys->params);
222
223 slv_destroy_parms(p);
224
225 if(p->parms==NULL){
226 CONSOLE_DEBUG("params NULL");
227 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
228 if(p->parms==NULL)return -1;
229 p->dynamic_parms = 1;
230 }else{
231 CONSOLE_DEBUG("params not NULL");
232 }
233
234 /* reset the number of parameters to zero so that we can check it at the end */
235 p->num_parms = 0;
236
237 slv_param_bool(p,IDA_PARAM_AUTODIFF
238 ,(SlvParameterInitBool){{"autodiff"
239 ,"Use auto-diff?",1
240 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
241 }, TRUE}
242 );
243
244 slv_param_bool(p,IDA_PARAM_CALCIC
245 ,(SlvParameterInitBool){{"calcic"
246 ,"Calculate initial conditions?",1
247 ,"Use IDA to calculate initial conditions (1) or else assume that the model will already be solved for this case (0)"
248 }, FALSE}
249 );
250
251 slv_param_bool(p,IDA_PARAM_SAFEEVAL
252 ,(SlvParameterInitBool){{"safeeval"
253 ,"Use safe evaluation?",1
254 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
255 "throw SIGFPE errors which will then halt integration."
256 }, FALSE}
257 );
258
259
260 slv_param_bool(p,IDA_PARAM_ATOLVECT
261 ,(SlvParameterInitBool){{"atolvect"
262 ,"Use 'ode_atol' values as specified?",1
263 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
264 " in the integration. If FALSE, a scalar absolute tolerance value"
265 " is shared by all variables. See IDA manual, section 5.5.1"
266 }, TRUE }
267 );
268
269 slv_param_real(p,IDA_PARAM_ATOL
270 ,(SlvParameterInitReal){{"atol"
271 ,"Scalar absolute error tolerance",1
272 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
273 " See IDA manual, section 5.5.1"
274 }, 1e-5, DBL_MIN, DBL_MAX }
275 );
276
277 slv_param_real(p,IDA_PARAM_RTOL
278 ,(SlvParameterInitReal){{"rtol"
279 ,"Scalar relative error tolerance",1
280 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
281 " it's not possible to set per-variable relative tolerances as it is"
282 " with LSODE)."
283 " See IDA manual, section 5.5.1"
284 }, 1e-4, 0, DBL_MAX }
285 );
286
287 slv_param_char(p,IDA_PARAM_LINSOLVER
288 ,(SlvParameterInitChar){{"linsolver"
289 ,"Linear solver",1
290 ,"See IDA manual, section 5.5.3."
291 }, "SPGMR"}, (char *[]){"DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
292 );
293
294 slv_param_int(p,IDA_PARAM_MAXL
295 ,(SlvParameterInitInt){{"maxl"
296 ,"Maximum Krylov dimension",0
297 ,"The maximum dimension of Krylov space used by the linear solver"
298 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
299 " The default of 0 results in IDA using its internal default, which"
300 " is currently a value of 5."
301 }, 0, 0, 20 }
302 );
303
304 slv_param_bool(p,IDA_PARAM_GSMODIFIED
305 ,(SlvParameterInitBool){{"gsmodified"
306 ,"Gram-Schmidt Orthogonalisation Scheme", 2
307 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
308 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
309 }, TRUE}
310 );
311
312 slv_param_int(p,IDA_PARAM_MAXNCF
313 ,(SlvParameterInitInt){{"maxncf"
314 ,"Max nonlinear solver convergence failures per step", 2
315 ,"Maximum number of allowable nonlinear solver convergence failures"
316 " on one step. See IDA manual section 5.5.6.1."
317 }, 10,0,1000 }
318 );
319
320 slv_param_char(p,IDA_PARAM_PREC
321 ,(SlvParameterInitChar){{"prec"
322 ,"Preconditioner",1
323 ,"See IDA manual, section section 5.6.8."
324 },"NONE"}, (char *[]){"NONE","DIAG",NULL}
325 );
326
327 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
328
329 CONSOLE_DEBUG("Created %d params", p->num_parms);
330
331 return 0;
332 }
333
334 /*-------------------------------------------------------------
335 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
336 */
337
338 typedef int IdaFlagFn(void *,int *);
339 typedef char *IdaFlagNameFn(int);
340
341 /* return 1 on success */
342 int integrator_ida_solve(
343 IntegratorSystem *blsys
344 , unsigned long start_index
345 , unsigned long finish_index
346 ){
347 void *ida_mem;
348 int size, flag, flag1, t_index;
349 realtype t0, reltol, abstol, t, tret, tout1;
350 N_Vector y0, yp0, abstolvect, ypret, yret;
351 IntegratorIdaData *enginedata;
352 char *linsolver;
353 int maxl;
354 IdaFlagFn *flagfn;
355 IdaFlagNameFn *flagnamefn;
356 const char *flagfntype;
357 char *pname;
358 IDASpilsPrecSetupFn psetup;
359 IDASpilsPrecSolveFn psolve;
360
361 CONSOLE_DEBUG("STARTING IDA...");
362
363 enginedata = integrator_ida_enginedata(blsys);
364
365 enginedata->safeeval = SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_SAFEEVAL);
366 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
367
368 /* store reference to list of relations (in enginedata) */
369 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
370 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
371 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
372 enginedata->bndlist = slv_get_solvers_bnd_list(blsys->system);
373
374 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
375 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
376 size = blsys->n_y;
377
378 if(enginedata->nrels!=size){
379 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
380 return 0; /* failure */
381 }
382
383 /* retrieve initial values from the system */
384
385 /** @TODO fix this, the starting time != first sample */
386 t0 = integrator_get_t(blsys);
387 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
388
389 CONSOLE_DEBUG("RETRIEVING y0");
390
391 y0 = N_VNew_Serial(size);
392 integrator_get_y(blsys,NV_DATA_S(y0));
393
394 #ifdef SOLVE_DEBUG
395 CONSOLE_DEBUG("RETRIEVING yp0");
396 #endif
397
398 yp0 = N_VNew_Serial(size);
399 integrator_get_ydot(blsys,NV_DATA_S(yp0));
400
401 #ifdef SOLVE_DEBUG
402 N_VPrint_Serial(yp0);
403 CONSOLE_DEBUG("yp0 is at %p",&yp0);
404 #endif
405
406 /* create IDA object */
407 ida_mem = IDACreate();
408
409 /* relative error tolerance */
410 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
411 CONSOLE_DEBUG("rtol = %8.2e",reltol);
412
413 /* allocate internal memory */
414 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
415 /* vector of absolute tolerances */
416 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
417 abstolvect = N_VNew_Serial(size);
418 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
419
420 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
421
422 N_VDestroy_Serial(abstolvect);
423 }else{
424 /* scalar absolute tolerance (one value for all) */
425 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
426 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
427 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
428 }
429
430 if(flag==IDA_MEM_NULL){
431 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
432 return 0;
433 }else if(flag==IDA_MEM_FAIL){
434 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
435 return 0;
436 }else if(flag==IDA_ILL_INPUT){
437 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
438 return 0;
439 }/* else success */
440
441 /* set optional inputs... */
442 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
443 IDASetRdata(ida_mem, (void *)blsys);
444 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
445 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
446 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
447 if(integrator_get_minstep(blsys)>0){
448 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
449 }
450
451 CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
452 IDASetMaxConvFails(ida_mem,SLV_PARAM_INT(&blsys->params,IDA_PARAM_MAXNCF));
453
454 /* there's no capability for setting *minimum* step size in IDA */
455
456
457 /* attach linear solver module, using the default value of maxl */
458 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
459 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
460 if(strcmp(linsolver,"DENSE")==0){
461 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
462 flag = IDADense(ida_mem, size);
463 switch(flag){
464 case IDADENSE_SUCCESS: break;
465 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 0;
466 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 0;
467 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
468 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
469 }
470
471 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
472 CONSOLE_DEBUG("USING AUTODIFF");
473 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
474 switch(flag){
475 case IDADENSE_SUCCESS: break;
476 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 0;
477 }
478 }else{
479 CONSOLE_DEBUG("USING NUMERICAL DIFF");
480 }
481
482 flagfntype = "IDADENSE";
483 flagfn = &IDADenseGetLastFlag;
484 flagnamefn = &IDADenseGetReturnFlagName;
485 }else{
486 /* remaining methods are all SPILS */
487 CONSOLE_DEBUG("IDA SPILS");
488
489 maxl = SLV_PARAM_INT(&(blsys->params),IDA_PARAM_MAXL);
490 CONSOLE_DEBUG("maxl = %d",maxl);
491
492 pname = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_PREC);
493 if(strcmp(pname,"NONE")==0){
494 psetup=NULL;
495 psolve=NULL;
496 }else if(strcmp(pname,"DIAG")==0){
497 psetup=&integrator_ida_psetup;
498 psolve=&integrator_ida_psolve;
499 }
500
501 if(strcmp(linsolver,"SPGMR")==0){
502 CONSOLE_DEBUG("IDA SPGMR");
503 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
504 }else if(strcmp(linsolver,"SPBCG")==0){
505 CONSOLE_DEBUG("IDA SPBCG");
506 flag = IDASpbcg(ida_mem, maxl);
507 }else if(strcmp(linsolver,"SPTFQMR")==0){
508 CONSOLE_DEBUG("IDA SPTFQMR");
509 flag = IDASptfqmr(ida_mem,maxl);
510 }else{
511 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
512 return 0;
513 }
514
515 enginedata->pp = N_VNew_Serial(blsys->n_y);
516 IDASpilsSetPreconditioner(ida_mem,psetup,psolve,(void *)blsys);
517 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
518
519 flagfntype = "IDASPILS";
520 flagfn = &IDASpilsGetLastFlag;
521 flagnamefn = &IDASpilsGetReturnFlagName;
522
523 if(flag==IDASPILS_MEM_NULL){
524 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
525 return 0;
526 }else if(flag==IDASPILS_MEM_FAIL){
527 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
528 return 0;
529 }/* else success */
530
531 /* assign the J*v function */
532 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
533 CONSOLE_DEBUG("USING AUTODIFF");
534 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
535 if(flag==IDASPILS_MEM_NULL){
536 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
537 return 0;
538 }else if(flag==IDASPILS_LMEM_NULL){
539 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
540 return 0;
541 }/* else success */
542 }else{
543 CONSOLE_DEBUG("USING NUMERICAL DIFF");
544 }
545
546 if(strcmp(linsolver,"SPGMR")==0){
547 /* select Gram-Schmidt orthogonalisation */
548 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_GSMODIFIED)){
549 CONSOLE_DEBUG("USING MODIFIED GS");
550 flag = IDASpilsSetGSType(ida_mem,MODIFIED_GS);
551 if(flag!=IDASPILS_SUCCESS){
552 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
553 return 0;
554 }
555 }else{
556 CONSOLE_DEBUG("USING CLASSICAL GS");
557 flag = IDASpilsSetGSType(ida_mem,CLASSICAL_GS);
558 if(flag!=IDASPILS_SUCCESS){
559 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
560 return 0;
561 }
562 }
563 }
564 }
565
566 /* set linear solver optional inputs...
567 ...nothing here at the moment...
568 */
569
570 /* calculate initial conditions */
571 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_CALCIC)){
572 CONSOLE_DEBUG("Solving initial conditions (given derivatives)");
573
574 blsys->currentstep=0;
575 t_index=start_index;
576 tout1 = samplelist_get(blsys->samples, t_index);
577
578 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
579
580 /* correct initial values, given derivatives */
581 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
582 /* note the new API from version 2.3 and onwards */
583 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
584 # else
585 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
586 # endif
587
588 switch(flag){
589 case IDA_SUCCESS:
590 CONSOLE_DEBUG("Initial conditions solved OK");
591 break;
592
593 case IDA_LSETUP_FAIL:
594 case IDA_LINIT_FAIL:
595 case IDA_LSOLVE_FAIL:
596 case IDA_NO_RECOVERY:
597 flag1 = -999;
598 flag = (flagfn)(ida_mem,&flag1);
599 if(flag){
600 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to retrieve error code from %s (err %d)",flagfntype,flag);
601 return 0;
602 }
603 ERROR_REPORTER_HERE(ASC_PROG_ERR,"%s returned flag '%s' (value = %d)",flagfntype,(flagnamefn)(flag1),flag1);
604 return 0;
605
606 default:
607 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
608 return 0;
609 }
610 }else{
611 CONSOLE_DEBUG("Not solving initial conditions (see IDA parameter 'calcic')");
612 }
613
614 /* optionally, specify ROO-FINDING PROBLEM */
615
616 /* -- set up the IntegratorReporter */
617 integrator_output_init(blsys);
618
619 /* -- store the initial values of all the stuff */
620 integrator_output_write(blsys);
621 integrator_output_write_obs(blsys);
622
623 /* specify where the returned values should be stored */
624 yret = y0;
625 ypret = yp0;
626
627 /* advance solution in time, return values as yret and derivatives as ypret */
628 blsys->currentstep=1;
629 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
630 t = samplelist_get(blsys->samples, t_index);
631 t0 = integrator_get_t(blsys);
632 asc_assert(t > t0);
633
634 #ifdef SOLVE_DEBUG
635 CONSOLE_DEBUG("Integratoring from t0 = %f to t = %f", t0, t);
636 #endif
637
638 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
639
640 /* pass the values of everything back to the compiler */
641 integrator_set_t(blsys, (double)tret);
642 integrator_set_y(blsys, NV_DATA_S(yret));
643 integrator_set_ydot(blsys, NV_DATA_S(ypret));
644
645 if(flag<0){
646 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
647 break;
648 }
649
650 /* -- do something so that blsys knows the values of tret, yret and ypret */
651
652 /* -- store the current values of all the stuff */
653 integrator_output_write(blsys);
654 integrator_output_write_obs(blsys);
655 }
656
657 /* -- close the IntegratorReporter */
658 integrator_output_close(blsys);
659
660 if(flag < 0){
661 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
662 return 0;
663 }
664
665 /* get optional outputs */
666
667 /* free solution memory */
668 N_VDestroy_Serial(yret);
669 N_VDestroy_Serial(ypret);
670
671 /* free solver memory */
672 IDAFree(ida_mem);
673
674 /* all done */
675 return 1;
676 }
677
678 /*--------------------------------------------------
679 RESIDUALS AND JACOBIAN
680 */
681 /**
682 Function to evaluate system residuals, in the form required for IDA.
683
684 Given tt, yy and yp, we need to evaluate and return rr.
685
686 @param tt current value of indep variable (time)
687 @param yy current values of dependent variable vector
688 @param yp current values of derivatives of dependent variables
689 @param rr the output residual vector (is we're returning data to)
690 @param res_data pointer to our stuff (blsys in this case).
691
692 @return 0 on success, positive on recoverable error, and
693 negative on unrecoverable error.
694 */
695 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
696 IntegratorSystem *blsys;
697 IntegratorIdaData *enginedata;
698 int i, calc_ok, is_error;
699 struct rel_relation** relptr;
700 double resid;
701 char *relname;
702 #ifdef FEX_DEBUG
703 char *varname;
704 char diffname[30];
705 #endif
706
707 blsys = (IntegratorSystem *)res_data;
708 enginedata = integrator_ida_enginedata(blsys);
709
710 #ifdef FEX_DEBUG
711 /* fprintf(stderr,"\n\n"); */
712 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
713 #endif
714
715 /* pass the values of everything back to the compiler */
716 integrator_set_t(blsys, (double)tt);
717 integrator_set_y(blsys, NV_DATA_S(yy));
718 integrator_set_ydot(blsys, NV_DATA_S(yp));
719
720 if(NV_LENGTH_S(rr)!=enginedata->nrels){
721 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
722 return -1; /* unrecoverable */
723 }
724
725 /**
726 @TODO does this function (fex) do bounds checking already?
727 */
728
729 /* evaluate each residual in the rellist */
730 is_error = 0;
731 relptr = enginedata->rellist;
732
733 if(enginedata->safeeval){
734 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
735 }else{
736 #ifdef FEX_DEBUG
737 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
738 #endif
739 Asc_SignalHandlerPushDefault(SIGFPE);
740 }
741
742 if (setjmp(g_fpe_env)==0) {
743 for(i=0, relptr = enginedata->rellist;
744 i< enginedata->nrels && relptr != NULL;
745 ++i, ++relptr
746 ){
747 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
748
749 NV_Ith_S(rr,i) = resid;
750 if(!calc_ok){
751 relname = rel_make_name(blsys->system, *relptr);
752 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
753 ASC_FREE(relname);
754 /* presumable some output already made? */
755 is_error = 1;
756 }/*else{
757 CONSOLE_DEBUG("Calc OK");
758 }*/
759 }
760 }else{
761 relname = rel_make_name(blsys->system, *relptr);
762 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
763 ASC_FREE(relname);
764 is_error = 1;
765 }
766
767 if(enginedata->safeeval){
768 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
769 }else{
770 Asc_SignalHandlerPopDefault(SIGFPE);
771 }
772
773 #ifdef FEX_DEBUG
774 /* output residuals to console */
775 CONSOLE_DEBUG("RESIDUAL OUTPUT");
776 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
777 for(i=0; i<blsys->n_y; ++i){
778 varname = var_make_name(blsys->system,blsys->y[i]);
779 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
780 if(blsys->ydot[i]){
781 varname = var_make_name(blsys->system,blsys->ydot[i]);
782 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
783 }else{
784 snprintf(diffname,99,"diff(%s)",varname);
785 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
786 }
787 ASC_FREE(varname);
788 relname = rel_make_name(blsys->system,enginedata->rellist[i]);
789 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
790 }
791 #endif
792
793 if(is_error){
794 return 1;
795 }
796
797 #ifdef FEX_DEBUG
798 CONSOLE_DEBUG("RESIDUAL OK");
799 #endif
800 return 0;
801 }
802
803 /**
804 Dense Jacobian evaluation. Only suitable for small problems!
805 */
806 int integrator_ida_djex(long int Neq, realtype tt
807 , N_Vector yy, N_Vector yp, N_Vector rr
808 , realtype c_j, void *jac_data, DenseMat Jac
809 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
810 ){
811 IntegratorSystem *blsys;
812 IntegratorIdaData *enginedata;
813 char *relname;
814 #ifdef JEX_DEBUG
815 char *varname;
816 #endif
817 int status;
818 struct rel_relation **relptr;
819 int i;
820 var_filter_t filter = {VAR_SVAR | VAR_FIXED, VAR_SVAR};
821 double *derivatives;
822 int *variables;
823 int count, j;
824 long var_yindex;
825
826 blsys = (IntegratorSystem *)jac_data;
827 enginedata = integrator_ida_enginedata(blsys);
828
829 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
830 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
831 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
832
833 /* pass the values of everything back to the compiler */
834 integrator_set_t(blsys, (double)tt);
835 integrator_set_y(blsys, NV_DATA_S(yy));
836 integrator_set_ydot(blsys, NV_DATA_S(yp));
837
838 #ifdef JEX_DEBUG
839 /* print vars */
840 for(i=0; i < blsys->n_y; ++i){
841 varname = var_make_name(blsys->system, blsys->y[i]);
842 CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
843 ASC_FREE(varname);
844 }
845
846 /* print derivatives */
847 for(i=0; i < blsys->n_y; ++i){
848 if(blsys->ydot[i]){
849 varname = var_make_name(blsys->system, blsys->ydot[i]);
850 CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
851 ASC_FREE(varname);
852 }else{
853 varname = var_make_name(blsys->system, blsys->y[i]);
854 CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
855 ASC_FREE(varname);
856 }
857 }
858
859 /* print step size */
860 CONSOLE_DEBUG("<c_j> = %f",c_j);
861 #endif
862
863 /* build up the dense jacobian matrix... */
864 status = 0;
865 for(i=0, relptr = enginedata->rellist;
866 i< enginedata->nrels && relptr != NULL;
867 ++i, ++relptr
868 ){
869 /* get derivatives for this particular relation */
870 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
871
872 if(status){
873 relname = rel_make_name(blsys->system, *relptr);
874 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
875 ASC_FREE(relname);
876 break;
877 }
878
879 /* output what's going on here ... */
880 #ifdef JEX_DEBUG
881 relname = rel_make_name(blsys->system, *relptr);
882 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
883 fprintf(stderr,"%d: '%s': ",i,relname);
884 ASC_FREE(relname);
885 for(j=0;j<count;++j){
886 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
887 var_yindex = blsys->y_id[variables[j]];
888 if(var_yindex >=0){
889 fprintf(stderr," var[%d]='%s'=y[%ld]",variables[j],varname,var_yindex);
890 }else{
891 fprintf(stderr," var[%d]='%s'=ydot[%ld]",variables[j],varname,-var_yindex-1);
892 }
893 ASC_FREE(varname);
894 }
895 fprintf(stderr,"\n");
896 #endif
897 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
898 for(j=0; j < count; ++j){
899 var_yindex = blsys->y_id[variables[j]];
900 #ifndef __WIN32__
901 /* the SUNDIALS headers seem not to store 'N' on Windows */
902 ASC_ASSERT_RANGE(var_yindex, -Jac->N, Jac->N);
903 #endif
904 if(var_yindex >= 0){
905 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
906 DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
907 }else{
908 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
909 DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
910 }
911 }
912 }
913
914 #ifdef JEX_DEBUG
915 CONSOLE_DEBUG("PRINTING JAC");
916 fprintf(stderr,"\t");
917 for(j=0; j < blsys->n_y; ++j){
918 if(j)fprintf(stderr,"\t");
919 varname = var_make_name(blsys->system,blsys->y[j]);
920 fprintf(stderr,"%11s",varname);
921 ASC_FREE(varname);
922 }
923 fprintf(stderr,"\n");
924 for(i=0; i < enginedata->nrels; ++i){
925 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
926 fprintf(stderr,"%s\t",relname);
927 ASC_FREE(relname);
928
929 for(j=0; j < blsys->n_y; ++j){
930 if(j)fprintf(stderr,"\t");
931 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
932 }
933 fprintf(stderr,"\n");
934 }
935 #endif
936
937 if(status){
938 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
939 return 1;
940 }
941
942 #ifdef FEX_DEBUG
943 CONSOLE_DEBUG("DJEX RETURNING 0");
944 #endif
945 return 0;
946 }
947
948 /**
949 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
950
951 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
952
953 @param tt current value of the independent variable (time, t)
954 @param yy current value of the dependent variable vector, y(t).
955 @param yp current value of y'(t).
956 @param rr current value of the residual vector F(t, y, y').
957 @param v the vector by which the Jacobian must be multiplied to the right.
958 @param Jv the output vector computed
959 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
960 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
961 @param tmp1 @see tmp2
962 @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.
963 @return 0 on success
964 */
965 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
966 , N_Vector v, N_Vector Jv, realtype c_j
967 , void *jac_data, N_Vector tmp1, N_Vector tmp2
968 ){
969 IntegratorSystem *blsys;
970 IntegratorIdaData *enginedata;
971 int i, j, is_error=0;
972 struct rel_relation** relptr;
973 char *relname;
974 int status;
975 double Jv_i;
976 long var_yindex;
977
978 int *variables;
979 double *derivatives;
980 var_filter_t filter;
981 int count;
982 #ifdef JEX_DEBUG
983 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
984 #endif
985
986 blsys = (IntegratorSystem *)jac_data;
987 enginedata = integrator_ida_enginedata(blsys);
988
989 /* pass the values of everything back to the compiler */
990 integrator_set_t(blsys, (double)tt);
991 integrator_set_y(blsys, NV_DATA_S(yy));
992 integrator_set_ydot(blsys, NV_DATA_S(yp));
993 /* no real use for residuals (rr) here, I don't think? */
994
995 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
996
997 i = NV_LENGTH_S(yy) * 2;
998 #ifdef JEX_DEBUG
999 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
1000 #endif
1001 variables = ASC_NEW_ARRAY(int, i);
1002 derivatives = ASC_NEW_ARRAY(double, i);
1003
1004 /* evaluate the derivatives... */
1005 /* J = dG_dy = dF_dy + alpha * dF_dyp */
1006
1007 filter.matchbits = VAR_SVAR | VAR_FIXED;
1008 filter.matchvalue = VAR_SVAR;
1009
1010 Asc_SignalHandlerPushDefault(SIGFPE);
1011 if (setjmp(g_fpe_env)==0) {
1012 for(i=0, relptr = enginedata->rellist;
1013 i< enginedata->nrels && relptr != NULL;
1014 ++i, ++relptr
1015 ){
1016 /* get derivatives for this particular relation */
1017 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
1018 #ifdef JEX_DEBUG
1019 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
1020 #endif
1021
1022 if(status){
1023 relname = rel_make_name(blsys->system, *relptr);
1024 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
1025 ASC_FREE(relname);
1026 is_error = 1;
1027 break;
1028 }
1029
1030 /*
1031 Now we have the derivatives wrt each alg/diff variable in the
1032 present equation. variables[] points into the varlist. need
1033 a mapping from the varlist to the y and ydot lists.
1034 */
1035
1036 Jv_i = 0;
1037 for(j=0; j < count; ++j){
1038 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
1039 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
1040 if(varname){
1041 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
1042 ASC_FREE(varname);
1043 }else{
1044 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
1045 }
1046 */
1047
1048 /* we don't calculate derivatives wrt indep var */
1049 asc_assert(variables[j]>=0);
1050 if(enginedata->varlist[variables[j]] == blsys->x) continue;
1051
1052 var_yindex = blsys->y_id[variables[j]];
1053 #ifdef JEX_DEBUG
1054 CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %ld",j,variables[j],var_yindex);
1055 #endif
1056
1057 ASC_ASSERT_RANGE(-var_yindex-1, -NV_LENGTH_S(v),NV_LENGTH_S(v));
1058
1059 if(var_yindex >= 0){
1060 #ifdef JEX_DEBUG
1061 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
1062 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%ld] = %f, v[%ld] = %f)\n", i
1063 , derivatives[j] * NV_Ith_S(v,var_yindex)
1064 , i, var_yindex, derivatives[j]
1065 , var_yindex, NV_Ith_S(v,var_yindex)
1066 );
1067 #endif
1068 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
1069 }else{
1070 #ifdef JEX_DEBUG
1071 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%ld] = %f, v[%ld] = %f)\n", i
1072 , derivatives[j] * NV_Ith_S(v,-var_yindex-1)
1073 , i, -var_yindex-1, derivatives[j]
1074 , -var_yindex-1, NV_Ith_S(v,-var_yindex-1)
1075 );
1076 #endif
1077 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
1078 Jv_i += derivatives[j] * NV_Ith_S(v,-var_yindex-1) * c_j;
1079 }
1080 }
1081
1082 NV_Ith_S(Jv,i) = Jv_i;
1083 #ifdef JEX_DEBUG
1084 CONSOLE_DEBUG("rel = %p",*relptr);
1085 relname = rel_make_name(blsys->system, *relptr);
1086 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
1087 //ASC_FREE(relname);
1088 return 1;
1089 #endif
1090 }
1091 }else{
1092 relname = rel_make_name(blsys->system, *relptr);
1093 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
1094 ASC_FREE(relname);
1095 is_error = 1;
1096 }
1097 Asc_SignalHandlerPopDefault(SIGFPE);
1098
1099 if(is_error){
1100 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
1101 return 1;
1102 }
1103 return 0;
1104 }
1105
1106 /*----------------------------------------------
1107 PRECONDITIONER
1108 */
1109
1110
1111 /**
1112 EXPERIMENTAL. A diagonal preconditioner for use with IDA Krylov solvers
1113
1114 'setup' function.
1115 */
1116 int integrator_ida_psetup(realtype tt,
1117 N_Vector yy, N_Vector yp, N_Vector rr,
1118 realtype c_j, void *p_data,
1119 N_Vector tmp1, N_Vector tmp2,
1120 N_Vector tmp3
1121 ){
1122 int i;
1123 IntegratorSystem *blsys;
1124 IntegratorIdaData *data;
1125 blsys = (IntegratorSystem *)p_data;
1126 data = blsys->enginedata;
1127
1128 double y = 1; /* derivative of y[i] in relation i */
1129 double ydot = 1; /* derivative of ydot[i] in relation i */
1130
1131 N_VConst(1, data->pp);
1132 for(i=0; i<blsys->n_y; ++i){
1133 /* @TODO calculate y, ydot here */
1134 NV_Ith_S(data->pp, i) = 1 / (y + c_j * ydot);
1135 }
1136 return 0;
1137 };
1138
1139 /**
1140 EXPERIMENTAL. A diagonal preconditioner for use with IDA Krylov solvers
1141
1142 'solve' function.
1143 */
1144 int integrator_ida_psolve(realtype tt,
1145 N_Vector yy, N_Vector yp, N_Vector rr,
1146 N_Vector rvec, N_Vector zvec,
1147 realtype c_j, realtype delta, void *p_data,
1148 N_Vector tmp
1149 ){
1150 IntegratorSystem *blsys;
1151 IntegratorIdaData *data;
1152 blsys = (IntegratorSystem *)p_data;
1153 data = blsys->enginedata;
1154
1155 N_VProd(data->pp, rvec, zvec);
1156 return 0;
1157 };
1158
1159 /*----------------------------------------------
1160 ERROR REPORTING
1161 */
1162 /**
1163 Error message reporter function to be passed to IDA. All error messages
1164 will trigger a call to this function, so we should find everything
1165 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
1166 panel (in the case of PyGTK).
1167 */
1168 void integrator_ida_error(int error_code
1169 , const char *module, const char *function
1170 , char *msg, void *eh_data
1171 ){
1172 IntegratorSystem *blsys;
1173 error_severity_t sev;
1174
1175 /* cast back the IntegratorSystem, just in case we need it */
1176 blsys = (IntegratorSystem *)eh_data;
1177
1178 /* severity depends on the sign of the error_code value */
1179 if(error_code <= 0){
1180 sev = ASC_PROG_ERR;
1181 }else{
1182 sev = ASC_PROG_WARNING;
1183 }
1184
1185 /* use our all-purpose error reporting to get stuff back to the GUI */
1186 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
1187 }

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