/[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 944 - (show annotations) (download) (as text)
Sat Nov 25 10:46:13 2006 UTC (14 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 20843 byte(s)
Implemented ATOLVECT, ATOL, RTOL parameters for the IDA integrator.
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_PARAM_RTOL
152 ,IDA_PARAM_ATOL
153 ,IDA_PARAM_ATOLVECT
154 ,IDA_PARAMS_SIZE
155 };
156
157 /**
158 Here the full set of parameters is defined, along with upper/lower bounds,
159 etc. The values are stuck into the blsys->params structure.
160
161 @return 0 on success
162 */
163 int integrator_ida_params_default(IntegratorSystem *blsys){
164 asc_assert(blsys!=NULL);
165 asc_assert(blsys->engine==INTEG_IDA);
166 slv_parameters_t *p;
167 p = &(blsys->params);
168
169 slv_destroy_parms(p);
170
171 if(p->parms==NULL){
172 CONSOLE_DEBUG("params NULL");
173 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
174 if(p->parms==NULL)return -1;
175 p->dynamic_parms = 1;
176 }else{
177 CONSOLE_DEBUG("params not NULL");
178 }
179
180 /* reset the number of parameters to zero so that we can check it at the end */
181 p->num_parms = 0;
182
183 slv_param_bool(p,IDA_PARAM_AUTODIFF
184 ,(SlvParameterInitBool){{"autodiff"
185 ,"Use auto-diff?",1
186 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
187 }, TRUE}
188 );
189
190 slv_param_bool(p,IDA_PARAM_ATOLVECT
191 ,(SlvParameterInitBool){{"atolvect"
192 ,"Use 'ode_atol' values as specified?",1
193 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
194 " in the integration. If FALSE, a scalar absolute tolerance value"
195 " is shared by all variables. See IDA manual, section 5.5.1"
196 }, TRUE }
197 );
198
199 slv_param_real(p,IDA_PARAM_ATOL
200 ,(SlvParameterInitReal){{"atol"
201 ,"Scalar absolute error tolerance",1
202 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
203 " See IDA manual, section 5.5.1"
204 }, 1e-5, DBL_MIN, DBL_MAX }
205 );
206
207 slv_param_real(p,IDA_PARAM_RTOL
208 ,(SlvParameterInitReal){{"rtol"
209 ,"Scalar relative error tolerance",1
210 ,"Value of the scalar relative error tolerance."
211 " See IDA manual, section 5.5.1"
212 }, 1e-5, DBL_MIN, DBL_MAX }
213 );
214
215 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
216
217 CONSOLE_DEBUG("Created %d params", p->num_parms);
218
219 return 0;
220 }
221
222 /*-------------------------------------------------------------
223 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
224 */
225
226 /* return 1 on success */
227 int integrator_ida_solve(
228 IntegratorSystem *blsys
229 , unsigned long start_index
230 , unsigned long finish_index
231 ){
232 void *ida_mem;
233 int size, flag, t_index;
234 realtype t0, reltol, abstol, t, tret, tout1;
235 N_Vector y0, yp0, abstolvect, ypret, yret;
236 IntegratorIdaData *enginedata;
237
238 CONSOLE_DEBUG("STARTING IDA...");
239
240 enginedata = integrator_ida_enginedata(blsys);
241 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
242
243 /* store reference to list of relations (in enginedata) */
244 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
245 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
246 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
247 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
248 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
249 size = blsys->n_y;
250
251 if(enginedata->nrels!=size){
252 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
253 return 0; /* failure */
254 }
255
256 /* retrieve initial values from the system */
257
258 /** @TODO fix this, the starting time != first sample */
259 t0 = integrator_get_t(blsys);
260 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
261
262 CONSOLE_DEBUG("RETRIEVING y0");
263
264 y0 = N_VNew_Serial(size);
265 integrator_get_y(blsys,NV_DATA_S(y0));
266
267 CONSOLE_DEBUG("RETRIEVING yp0");
268
269 yp0 = N_VNew_Serial(size);
270 integrator_get_ydot(blsys,NV_DATA_S(yp0));
271
272 N_VPrint_Serial(yp0);
273 CONSOLE_DEBUG("yp0 is at %p",&yp0);
274
275 /* create IDA object */
276 ida_mem = IDACreate();
277
278 /* relative error tolerance */
279 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
280
281 /* allocate internal memory */
282 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
283 /* vector of absolute tolerances */
284 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
285 abstolvect = N_VNew_Serial(size);
286 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
287
288 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
289 }else{
290 /* scalar absolute tolerance (one value for all) */
291 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
292 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
293 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
294
295 }
296
297 if(flag==IDA_MEM_NULL){
298 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
299 return 0;
300 }else if(flag==IDA_MEM_FAIL){
301 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
302 return 0;
303 }else if(flag==IDA_ILL_INPUT){
304 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
305 return 0;
306 }/* else success */
307
308 /* set optional inputs... */
309 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
310 IDASetRdata(ida_mem, (void *)blsys);
311 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
312 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
313 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
314 /* there's no capability for setting *minimum* step size in IDA */
315
316 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
317
318 /* attach linear solver module, using the default value of maxl */
319 flag = IDASpgmr(ida_mem, 0);
320 if(flag==IDASPILS_MEM_NULL){
321 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
322 return 0;
323 }else if(flag==IDASPILS_MEM_FAIL){
324 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
325 return 0;
326 }/* else success */
327
328 /* assign the J*v function */
329 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
330 CONSOLE_DEBUG("USING AUTODIFF");
331 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
332 if(flag==IDASPILS_MEM_NULL){
333 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
334 return 0;
335 }else if(flag==IDASPILS_LMEM_NULL){
336 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
337 return 0;
338 }/* else success */
339 }else{
340 CONSOLE_DEBUG("USING NUMERICAL DIFF");
341 }
342
343 /* set linear solver optional inputs...
344
345 ...nothing here at the moment...
346
347 */
348
349 /* correct initial values, given derivatives */
350 blsys->currentstep=0;
351 t_index=start_index;
352 tout1 = samplelist_get(blsys->samples, t_index);
353
354 /* CONSOLE_DEBUG("Giving t value %f to IDACalcIC", tout1);*/
355
356 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
357 /* note the new API from version 2.3 and onwards */
358 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
359 #else
360 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
361 #endif
362
363 if(flag!=IDA_SUCCESS){
364 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
365 return 0;
366 }/* else success */
367
368 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
369
370 /* optionally, specify ROO-FINDING PROBLEM */
371
372 /* -- set up the IntegratorReporter */
373 integrator_output_init(blsys);
374
375 /* -- store the initial values of all the stuff */
376 integrator_output_write(blsys);
377 integrator_output_write_obs(blsys);
378
379 /* specify where the returned values should be stored */
380 yret = y0;
381 ypret = yp0;
382
383 /* advance solution in time, return values as yret and derivatives as ypret */
384 blsys->currentstep=1;
385 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
386 t = samplelist_get(blsys->samples, t_index);
387
388 /* CONSOLE_DEBUG("SOLVING UP TO t = %f", t); */
389
390 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
391
392 /* pass the values of everything back to the compiler */
393 integrator_set_t(blsys, (double)tret);
394 integrator_set_y(blsys, NV_DATA_S(yret));
395 integrator_set_ydot(blsys, NV_DATA_S(ypret));
396
397 if(flag<0){
398 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
399 break;
400 }
401
402 /* -- do something so that blsys knows the values of tret, yret and ypret */
403
404 /* -- store the current values of all the stuff */
405 integrator_output_write(blsys);
406 integrator_output_write_obs(blsys);
407 }
408
409 /* -- close the IntegratorReporter */
410 integrator_output_close(blsys);
411
412 if(flag < 0){
413 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
414 return 0;
415 }
416
417 /* get optional outputs */
418
419 /* free solution memory */
420 N_VDestroy_Serial(yret);
421 N_VDestroy_Serial(ypret);
422
423 /* free solver memory */
424 IDAFree(ida_mem);
425
426 /* all done */
427 return 1;
428 }
429
430 /*--------------------------------------------------
431 RESIDUALS AND JACOBIAN
432 */
433 /**
434 Function to evaluate system residuals, in the form required for IDA.
435
436 Given tt, yy and yp, we need to evaluate and return rr.
437
438 @param tt current value of indep variable (time)
439 @param yy current values of dependent variable vector
440 @param yp current values of derivatives of dependent variables
441 @param rr the output residual vector (is we're returning data to)
442 @param res_data pointer to our stuff (blsys in this case).
443
444 @return 0 on success, positive on recoverable error, and
445 negative on unrecoverable error.
446 */
447 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
448 IntegratorSystem *blsys;
449 IntegratorIdaData *enginedata;
450 int i, calc_ok, is_error;
451 struct rel_relation** relptr;
452 double resid;
453 char *relname;
454
455 blsys = (IntegratorSystem *)res_data;
456 enginedata = integrator_ida_enginedata(blsys);
457
458 /* fprintf(stderr,"\n\n"); */
459 /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
460
461 /* pass the values of everything back to the compiler */
462 integrator_set_t(blsys, (double)tt);
463 integrator_set_y(blsys, NV_DATA_S(yy));
464 integrator_set_ydot(blsys, NV_DATA_S(yp));
465
466 /* revaluate the system residuals using the new data */
467 is_error = 0;
468 relptr = enginedata->rellist;
469
470 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
471 if(NV_LENGTH_S(rr)!=enginedata->nrels){
472 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
473 return -1; /* unrecoverable */
474 }
475
476 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
477 if (setjmp(g_fpe_env)==0) {
478 for(i=0, relptr = enginedata->rellist;
479 i< enginedata->nrels && relptr != NULL;
480 ++i, ++relptr
481 ){
482 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
483
484 relname = rel_make_name(blsys->system, *relptr);
485 /* if(calc_ok){
486 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
487 }else{
488 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
489 }*/
490 ASC_FREE(relname);
491
492 NV_Ith_S(rr,i) = resid;
493 if(!calc_ok){
494 /* presumable some output already made? */
495 is_error = 1;
496 }
497 }
498 }else{
499 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
500 }
501 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
502
503 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
504 return is_error;
505 }
506
507 /**
508 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
509
510 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
511
512 @param tt current value of the independent variable (time, t)
513 @param yy current value of the dependent variable vector, y(t).
514 @param yp current value of y'(t).
515 @param rr current value of the residual vector F(t, y, y').
516 @param v the vector by which the Jacobian must be multiplied to the right.
517 @param Jv the output vector computed
518 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
519 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
520 @param tmp1 @see tmp2
521 @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.
522 @return 0 on success
523 */
524 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
525 , N_Vector v, N_Vector Jv, realtype c_j
526 , void *jac_data, N_Vector tmp1, N_Vector tmp2
527 ){
528 IntegratorSystem *blsys;
529 IntegratorIdaData *enginedata;
530 int i, j, is_error=0;
531 struct rel_relation** relptr;
532 char *relname, *varname;
533 int status;
534 double Jv_i;
535 int var_yindex;
536
537 int *variables;
538 double *derivatives;
539 var_filter_t filter;
540 int count;
541
542 /* fprintf(stderr,"\n--------------\n"); */
543 /* CONSOLE_DEBUG("EVALUTING JACOBIAN..."); */
544
545 blsys = (IntegratorSystem *)jac_data;
546 enginedata = integrator_ida_enginedata(blsys);
547
548 /* pass the values of everything back to the compiler */
549 integrator_set_t(blsys, (double)tt);
550 integrator_set_y(blsys, NV_DATA_S(yy));
551 integrator_set_ydot(blsys, NV_DATA_S(yp));
552 /* no real use for residuals (rr) here, I don't think? */
553
554 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
555 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
556 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
557
558 /* evaluate the derivatives... */
559 /* J = dG_dy = dF_dy + alpha * dF_dyp */
560
561 filter.matchbits = VAR_SVAR;
562 filter.matchvalue = VAR_SVAR;
563
564 /* CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v)); */
565 /* for(i=0; i<NV_LENGTH_S(v); ++i){
566 CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));
567 }*/
568
569 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
570 if (setjmp(g_fpe_env)==0) {
571 for(i=0, relptr = enginedata->rellist;
572 i< enginedata->nrels && relptr != NULL;
573 ++i, ++relptr
574 ){
575 /* fprintf(stderr,"\n"); */
576 relname = rel_make_name(blsys->system, *relptr);
577 /* CONSOLE_DEBUG("RELATION %d '%s'",i,relname); */
578 ASC_FREE(relname);
579
580 /* get derivatives for this particular relation */
581 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
582 /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
583
584 for(j=0;j<count;++j){
585 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
586 /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */
587 ASC_FREE(varname);
588 }
589
590 if(!status){
591 /* CONSOLE_DEBUG("Derivatives for relation %d OK",i); */
592 }else{
593 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
594 break;
595 }
596
597 /*
598 Now we have the derivatives wrt each alg/diff variable in the
599 present equation. variables[] points into the varlist. need
600 a mapping from the varlist to the y and ydot lists.
601 */
602
603 Jv_i = 0;
604 for(j=0; j < count; ++j){
605 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
606 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
607 if(varname){
608 /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]); */
609 ASC_FREE(varname);
610 }else{
611 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
612 }
613
614 var_yindex = blsys->y_id[variables[j]];
615 /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
616
617 if(var_yindex >= 0){
618 /* CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
619 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
620 }else{
621 var_yindex = -var_yindex-1;
622 /* CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
623 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
624 }
625 }
626
627 NV_Ith_S(Jv,i) = Jv_i;
628 /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */
629
630 if(status){
631 /* presumably some error_reporter will already have been made*/
632 is_error = 1;
633 }
634 }
635 }else{
636 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
637 }
638 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
639
640 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
641
642
643
644 return is_error;
645 }
646
647 /*----------------------------------------------
648 ERROR REPORTING
649 */
650 /**
651 Error message reporter function to be passed to IDA. All error messages
652 will trigger a call to this function, so we should find everything
653 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
654 panel (in the case of PyGTK).
655 */
656 void integrator_ida_error(int error_code
657 , const char *module, const char *function
658 , char *msg, void *eh_data
659 ){
660 IntegratorSystem *blsys;
661 error_severity_t sev;
662
663 /* cast back the IntegratorSystem, just in case we need it */
664 blsys = (IntegratorSystem *)eh_data;
665
666 /* severity depends on the sign of the error_code value */
667 if(error_code <= 0){
668 sev = ASC_PROG_ERR;
669 }else{
670 sev = ASC_PROG_WARNING;
671 }
672
673 /* use our all-purpose error reporting to get stuff back to the GUI */
674 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
675 }

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