/[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 942 - (show annotations) (download) (as text)
Sat Nov 25 05:26:47 2006 UTC (14 years, 2 months ago) by johnpye
File MIME type: text/x-csrc
File size: 19291 byte(s)
Incorporated 'SolverParameters' functionality into 'Integrator', both at C level and C++/Python.
This shouldn't break Tcl/Tk as default parameters will be set and used invisibly.
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 <ida/ida.h>
60 # include <nvector/nvector_serial.h>
61 # include <ida/ida_spgmr.h>
62 # ifndef IDA_SUCCESS
63 # error "Failed to include SUNDIALS IDA header file"
64 # endif
65 #endif
66
67 /*
68 for the benefit of build tools that didn't sniff the SUNDIALS version, we
69 assume version 2.2.x (and thence possible errors).
70 */
71 #ifndef SUNDIALS_VERSION_MINOR
72 # ifdef __GNUC__
73 # warning "GUESSING SUNDIALS VERSION 2.2"
74 # endif
75 # define SUNDIALS_VERSION_MINOR 2
76 #endif
77 #ifndef SUNDIALS_VERSION_MAJOR
78 # define SUNDIALS_VERSION_MAJOR 2
79 #endif
80
81 /* check that we've got what we expect now */
82 #ifndef ASC_IDA_H
83 # error "Failed to include ASCEND IDA header file"
84 #endif
85
86 /**
87 Struct containing any stuff that IDA needs that doesn't fit into the
88 common IntegratorSystem struct.
89 */
90 typedef struct{
91 struct rel_relation **rellist; /**< NULL terminated list of rels */
92 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
93 int nrels;
94 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
95 } IntegratorIdaData;
96
97 /*-------------------------------------------------------------
98 FORWARD DECLS
99 */
100 /* residual function forward declaration */
101 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
102
103 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
104 , N_Vector v, N_Vector Jv, realtype c_j
105 , void *jac_data, N_Vector tmp1, N_Vector tmp2
106 );
107
108 /* error handler forward declaration */
109 void integrator_ida_error(int error_code
110 , const char *module, const char *function
111 , char *msg, void *eh_data
112 );
113
114 /*-------------------------------------------------------------
115 SETUP/TEARDOWN ROUTINES
116 */
117 void integrator_ida_create(IntegratorSystem *blsys){
118 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
119 IntegratorIdaData *enginedata;
120 enginedata = ASC_NEW(IntegratorIdaData);
121 enginedata->rellist = NULL;
122 enginedata->varlist = NULL;
123 enginedata->safeeval = 1;
124 blsys->enginedata = (void *)enginedata;
125 integrator_ida_params_default(blsys);
126 }
127
128 void integrator_ida_free(void *enginedata){
129 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
130 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
131 /* note, we don't own the rellist, so don't need to free it */
132 ASC_FREE(d);
133 }
134
135 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
136 IntegratorIdaData *d;
137 assert(blsys!=NULL);
138 assert(blsys->enginedata!=NULL);
139 assert(blsys->engine==INTEG_IDA);
140 d = ((IntegratorIdaData *)(blsys->enginedata));
141 assert(d->safeeval = 1);
142 return d;
143 }
144
145 /*-------------------------------------------------------------
146 PARAMETERS FOR IDA
147 */
148
149 enum ida_parameters{
150 IDA_PARAM_AUTODIFF
151 ,IDA_PARAMS_SIZE
152 };
153
154 /**
155 Here the full set of parameters is defined, along with upper/lower bounds,
156 etc. The values are stuck into the blsys->params structure.
157
158 @return 0 on success
159 */
160 int integrator_ida_params_default(IntegratorSystem *blsys){
161 asc_assert(blsys!=NULL);
162 asc_assert(blsys->engine==INTEG_IDA);
163 slv_parameters_t *p;
164 p = &(blsys->params);
165
166 slv_destroy_parms(p);
167
168 if(p->parms==NULL){
169 CONSOLE_DEBUG("params NULL");
170 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
171 if(p->parms==NULL)return -1;
172 p->dynamic_parms = 1;
173 }else{
174 CONSOLE_DEBUG("params not NULL");
175 }
176
177 /* reset the number of parameters to zero so that we can check it at the end */
178 p->num_parms = 0;
179
180 slv_param_bool(p,IDA_PARAM_AUTODIFF
181 ,(SlvParameterInitBool){{"autodiff"
182 ,"Use auto-diff?",1
183 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
184 }, TRUE}
185 );
186
187 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
188
189 CONSOLE_DEBUG("Created %d params", p->num_parms);
190
191 return 0;
192 }
193
194 /*-------------------------------------------------------------
195 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
196 */
197
198 /* return 1 on success */
199 int integrator_ida_solve(
200 IntegratorSystem *blsys
201 , unsigned long start_index
202 , unsigned long finish_index
203 ){
204 void *ida_mem;
205 int size, flag, t_index;
206 realtype t0, reltol, t, tret, tout1;
207 N_Vector y0, yp0, abstol, ypret, yret;
208 IntegratorIdaData *enginedata;
209
210 CONSOLE_DEBUG("STARTING IDA...");
211
212 enginedata = integrator_ida_enginedata(blsys);
213 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
214
215 /* store reference to list of relations (in enginedata) */
216 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
217 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
218 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
219 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
220 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
221 size = blsys->n_y;
222
223 if(enginedata->nrels!=size){
224 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
225 return 0; /* failure */
226 }
227
228 CONSOLE_DEBUG("RETRIEVING t0");
229
230 /* retrieve initial values from the system */
231 t0 = samplelist_get(blsys->samples,start_index);
232
233 CONSOLE_DEBUG("RETRIEVING y0");
234
235 y0 = N_VNew_Serial(size);
236 integrator_get_y(blsys,NV_DATA_S(y0));
237
238 CONSOLE_DEBUG("RETRIEVING yp0");
239
240 yp0 = N_VNew_Serial(size);
241 integrator_get_ydot(blsys,NV_DATA_S(yp0));
242
243 N_VPrint_Serial(yp0);
244 CONSOLE_DEBUG("yp0 is at %p",&yp0);
245
246 /* create IDA object */
247 ida_mem = IDACreate();
248
249 /* retrieve the absolute tolerance values for each variable */
250 abstol = N_VNew_Serial(size);
251 N_VConst(0.1,abstol); /** @TODO fill in the abstol values from the model */
252 reltol = 0.001;
253
254 /* allocate internal memory */
255 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstol);
256 if(flag==IDA_MEM_NULL){
257 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
258 return 0;
259 }else if(flag==IDA_MEM_FAIL){
260 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
261 return 0;
262 }else if(flag==IDA_ILL_INPUT){
263 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
264 return 0;
265 }/* else success */
266
267 /* set optional inputs... */
268 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
269 IDASetRdata(ida_mem, (void *)blsys);
270 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
271 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
272 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
273 /* there's no capability for setting *minimum* step size in IDA */
274
275 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
276
277 /* attach linear solver module, using the default value of maxl */
278 flag = IDASpgmr(ida_mem, 0);
279 if(flag==IDASPILS_MEM_NULL){
280 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
281 return 0;
282 }else if(flag==IDASPILS_MEM_FAIL){
283 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
284 return 0;
285 }/* else success */
286
287 /* assign the J*v function */
288 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
289 CONSOLE_DEBUG("USING AUTODIFF");
290 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
291 if(flag==IDASPILS_MEM_NULL){
292 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
293 return 0;
294 }else if(flag==IDASPILS_LMEM_NULL){
295 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
296 return 0;
297 }/* else success */
298 }else{
299 CONSOLE_DEBUG("USING NUMERICAL DIFF");
300 }
301
302 /* set linear solver optional inputs...
303
304 ...nothing here at the moment...
305
306 */
307
308 /* correct initial values, given derivatives */
309 blsys->currentstep=0;
310 t_index=start_index+1;
311 tout1 = samplelist_get(blsys->samples, t_index);
312
313 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
314 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
315 #else
316 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
317 #endif
318
319 if(flag!=IDA_SUCCESS){
320 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
321 return 0;
322 }/* else success */
323
324 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
325
326 /* optionally, specify ROO-FINDING PROBLEM */
327
328 /* -- set up the IntegratorReporter */
329 integrator_output_init(blsys);
330
331 /* -- store the initial values of all the stuff */
332 integrator_output_write(blsys);
333 integrator_output_write_obs(blsys);
334
335 /* specify where the returned values should be stored */
336 yret = y0;
337 ypret = yp0;
338
339 /* advance solution in time, return values as yret and derivatives as ypret */
340 blsys->currentstep=1;
341 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
342 t = samplelist_get(blsys->samples, t_index);
343
344 /* CONSOLE_DEBUG("SOLVING UP TO t = %f", t); */
345
346 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
347
348 /* pass the values of everything back to the compiler */
349 integrator_set_t(blsys, (double)tret);
350 integrator_set_y(blsys, NV_DATA_S(yret));
351 integrator_set_ydot(blsys, NV_DATA_S(ypret));
352
353 if(flag<0){
354 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
355 break;
356 }
357
358 /* -- do something so that blsys knows the values of tret, yret and ypret */
359
360 /* -- store the current values of all the stuff */
361 integrator_output_write(blsys);
362 integrator_output_write_obs(blsys);
363 }
364
365 /* -- close the IntegratorReporter */
366 integrator_output_close(blsys);
367
368 if(flag < 0){
369 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
370 return 0;
371 }
372
373 /* get optional outputs */
374
375 /* free solution memory */
376 N_VDestroy_Serial(yret);
377 N_VDestroy_Serial(ypret);
378
379 /* free solver memory */
380 IDAFree(ida_mem);
381
382 /* all done */
383 return 1;
384 }
385
386 /*--------------------------------------------------
387 RESIDUALS AND JACOBIAN
388 */
389 /**
390 Function to evaluate system residuals, in the form required for IDA.
391
392 Given tt, yy and yp, we need to evaluate and return rr.
393
394 @param tt current value of indep variable (time)
395 @param yy current values of dependent variable vector
396 @param yp current values of derivatives of dependent variables
397 @param rr the output residual vector (is we're returning data to)
398 @param res_data pointer to our stuff (blsys in this case).
399
400 @return 0 on success, positive on recoverable error, and
401 negative on unrecoverable error.
402 */
403 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
404 IntegratorSystem *blsys;
405 IntegratorIdaData *enginedata;
406 int i, calc_ok, is_error;
407 struct rel_relation** relptr;
408 double resid;
409 char *relname;
410
411 blsys = (IntegratorSystem *)res_data;
412 enginedata = integrator_ida_enginedata(blsys);
413
414 /* fprintf(stderr,"\n\n"); */
415 /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
416
417 /* pass the values of everything back to the compiler */
418 integrator_set_t(blsys, (double)tt);
419 integrator_set_y(blsys, NV_DATA_S(yy));
420 integrator_set_ydot(blsys, NV_DATA_S(yp));
421
422 /* revaluate the system residuals using the new data */
423 is_error = 0;
424 relptr = enginedata->rellist;
425
426 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
427 if(NV_LENGTH_S(rr)!=enginedata->nrels){
428 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
429 return -1; /* unrecoverable */
430 }
431
432 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
433 if (setjmp(g_fpe_env)==0) {
434 for(i=0, relptr = enginedata->rellist;
435 i< enginedata->nrels && relptr != NULL;
436 ++i, ++relptr
437 ){
438 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
439
440 relname = rel_make_name(blsys->system, *relptr);
441 /* if(calc_ok){
442 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
443 }else{
444 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
445 }*/
446 ASC_FREE(relname);
447
448 NV_Ith_S(rr,i) = resid;
449 if(!calc_ok){
450 /* presumable some output already made? */
451 is_error = 1;
452 }
453 }
454 }else{
455 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
456 }
457 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
458
459 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
460 return is_error;
461 }
462
463 /**
464 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
465
466 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
467
468 @param tt current value of the independent variable (time, t)
469 @param yy current value of the dependent variable vector, y(t).
470 @param yp current value of y'(t).
471 @param rr current value of the residual vector F(t, y, y').
472 @param v the vector by which the Jacobian must be multiplied to the right.
473 @param Jv the output vector computed
474 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
475 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
476 @param tmp1 @see tmp2
477 @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.
478 @return 0 on success
479 */
480 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
481 , N_Vector v, N_Vector Jv, realtype c_j
482 , void *jac_data, N_Vector tmp1, N_Vector tmp2
483 ){
484 IntegratorSystem *blsys;
485 IntegratorIdaData *enginedata;
486 int i, j, is_error=0;
487 struct rel_relation** relptr;
488 char *relname, *varname;
489 int status;
490 double Jv_i;
491 int var_yindex;
492
493 int *variables;
494 double *derivatives;
495 var_filter_t filter;
496 int count;
497
498 fprintf(stderr,"\n--------------\n");
499 CONSOLE_DEBUG("EVALUTING JACOBIAN...");
500
501 blsys = (IntegratorSystem *)jac_data;
502 enginedata = integrator_ida_enginedata(blsys);
503
504 /* pass the values of everything back to the compiler */
505 integrator_set_t(blsys, (double)tt);
506 integrator_set_y(blsys, NV_DATA_S(yy));
507 integrator_set_ydot(blsys, NV_DATA_S(yp));
508 /* no real use for residuals (rr) here, I don't think? */
509
510 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
511 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
512 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
513
514 /* evaluate the derivatives... */
515 /* J = dG_dy = dF_dy + alpha * dF_dyp */
516
517 filter.matchbits = VAR_SVAR;
518 filter.matchvalue = VAR_SVAR;
519
520 CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v));
521 for(i=0; i<NV_LENGTH_S(v); ++i){
522 CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));
523 }
524
525 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
526 if (setjmp(g_fpe_env)==0) {
527 for(i=0, relptr = enginedata->rellist;
528 i< enginedata->nrels && relptr != NULL;
529 ++i, ++relptr
530 ){
531 fprintf(stderr,"\n");
532 relname = rel_make_name(blsys->system, *relptr);
533 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
534 ASC_FREE(relname);
535
536 /* get derivatives for this particular relation */
537 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
538 CONSOLE_DEBUG("Got derivatives against %d matching variables", count);
539
540 for(j=0;j<count;++j){
541 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
542 CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);
543 ASC_FREE(varname);
544 }
545
546 if(!status){
547 CONSOLE_DEBUG("Derivatives for relation %d OK",i);
548 }else{
549 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
550 break;
551 }
552
553 /*
554 Now we have the derivatives wrt each alg/diff variable in the
555 present equation. variables[] points into the varlist. need
556 a mapping from the varlist to the y and ydot lists.
557 */
558
559 Jv_i = 0;
560 for(j=0; j < count; ++j){
561 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
562 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
563 if(varname){
564 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
565 ASC_FREE(varname);
566 }else{
567 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
568 }
569
570 var_yindex = blsys->y_id[variables[j]];
571 CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex);
572
573 if(var_yindex >= 0){
574 CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));
575 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
576 }else{
577 var_yindex = -var_yindex-1;
578 CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex));
579 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
580 }
581 }
582
583 NV_Ith_S(Jv,i) = Jv_i;
584 CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i);
585
586 if(status){
587 /* presumably some error_reporter will already have been made*/
588 is_error = 1;
589 }
590 }
591 }else{
592 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
593 }
594 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
595
596 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
597
598
599
600 return is_error;
601 }
602
603 /*----------------------------------------------
604 ERROR REPORTING
605 */
606 /**
607 Error message reporter function to be passed to IDA. All error messages
608 will trigger a call to this function, so we should find everything
609 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
610 panel (in the case of PyGTK).
611 */
612 void integrator_ida_error(int error_code
613 , const char *module, const char *function
614 , char *msg, void *eh_data
615 ){
616 IntegratorSystem *blsys;
617 error_severity_t sev;
618
619 /* cast back the IntegratorSystem, just in case we need it */
620 blsys = (IntegratorSystem *)eh_data;
621
622 /* severity depends on the sign of the error_code value */
623 if(error_code <= 0){
624 sev = ASC_PROG_ERR;
625 }else{
626 sev = ASC_PROG_WARNING;
627 }
628
629 /* use our all-purpose error reporting to get stuff back to the GUI */
630 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
631 }

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