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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 944 - (hide annotations) (download) (as text)
Sat Nov 25 10:46:13 2006 UTC (14 years, 3 months ago) by johnpye
File MIME type: text/x-csrc
File size: 20843 byte(s)
Implemented ATOLVECT, ATOL, RTOL parameters for the IDA integrator.
1 johnpye 669 /* 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 johnpye 912 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 johnpye 669 */
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 johnpye 942 #include <utilities/ascPanic.h>
43 johnpye 669 #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 johnpye 719 #include "relman.h"
55 johnpye 669
56     /* SUNDIALS includes */
57     #ifdef ASC_WITH_IDA
58 johnpye 913 # include <sundials/sundials_config.h>
59 johnpye 782 # include <ida/ida.h>
60     # include <nvector/nvector_serial.h>
61     # include <ida/ida_spgmr.h>
62 johnpye 669 # ifndef IDA_SUCCESS
63     # error "Failed to include SUNDIALS IDA header file"
64     # endif
65     #endif
66    
67 johnpye 913 /*
68     for the benefit of build tools that didn't sniff the SUNDIALS version, we
69 johnpye 924 assume version 2.2.x (and thence possible errors).
70 johnpye 913 */
71     #ifndef SUNDIALS_VERSION_MINOR
72 johnpye 921 # ifdef __GNUC__
73     # warning "GUESSING SUNDIALS VERSION 2.2"
74     # endif
75 johnpye 913 # define SUNDIALS_VERSION_MINOR 2
76     #endif
77     #ifndef SUNDIALS_VERSION_MAJOR
78     # define SUNDIALS_VERSION_MAJOR 2
79     #endif
80    
81 johnpye 669 /* 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 johnpye 912 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
93 johnpye 669 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 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
102 johnpye 669
103 johnpye 909 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 johnpye 669 /* 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 johnpye 909 enginedata->varlist = NULL;
123 johnpye 669 enginedata->safeeval = 1;
124 ben.allan 704 blsys->enginedata = (void *)enginedata;
125 johnpye 942 integrator_ida_params_default(blsys);
126 ben.allan 704 }
127 johnpye 669
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 johnpye 942 PARAMETERS FOR IDA
147     */
148    
149     enum ida_parameters{
150     IDA_PARAM_AUTODIFF
151 johnpye 944 ,IDA_PARAM_RTOL
152     ,IDA_PARAM_ATOL
153     ,IDA_PARAM_ATOLVECT
154 johnpye 942 ,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 johnpye 944 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 johnpye 942 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 johnpye 669 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 johnpye 719 int size, flag, t_index;
234 johnpye 944 realtype t0, reltol, abstol, t, tret, tout1;
235     N_Vector y0, yp0, abstolvect, ypret, yret;
236 johnpye 669 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 johnpye 909 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
247 johnpye 669 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
248 johnpye 719 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
249 johnpye 669 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 johnpye 944
258     /** @TODO fix this, the starting time != first sample */
259     t0 = integrator_get_t(blsys);
260     CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
261 johnpye 669
262 johnpye 928 CONSOLE_DEBUG("RETRIEVING y0");
263    
264 johnpye 669 y0 = N_VNew_Serial(size);
265     integrator_get_y(blsys,NV_DATA_S(y0));
266    
267 johnpye 928 CONSOLE_DEBUG("RETRIEVING yp0");
268    
269 johnpye 669 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 johnpye 944 /* relative error tolerance */
279     reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
280 johnpye 669
281     /* allocate internal memory */
282 johnpye 944 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 johnpye 669 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 johnpye 912 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER");
317    
318 johnpye 669 /* 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 johnpye 909 /* assign the J*v function */
329 johnpye 942 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 johnpye 669
343 johnpye 909 /* set linear solver optional inputs...
344 johnpye 669
345 johnpye 909 ...nothing here at the moment...
346    
347     */
348    
349 johnpye 669 /* correct initial values, given derivatives */
350     blsys->currentstep=0;
351 johnpye 944 t_index=start_index;
352 johnpye 669 tout1 = samplelist_get(blsys->samples, t_index);
353 johnpye 913
354 johnpye 944 /* CONSOLE_DEBUG("Giving t value %f to IDACalcIC", tout1);*/
355    
356 johnpye 913 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
357 johnpye 944 /* note the new API from version 2.3 and onwards */
358 johnpye 913 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
359     #else
360 johnpye 669 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
361 johnpye 913 #endif
362    
363 johnpye 669 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 johnpye 942 /* CONSOLE_DEBUG("SOLVING UP TO t = %f", t); */
389 johnpye 669
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 johnpye 903 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
448 johnpye 669 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 johnpye 903 /* fprintf(stderr,"\n\n"); */
459     /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
460 johnpye 669
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 johnpye 903 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
471 johnpye 669 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 johnpye 903 /* if(calc_ok){
486 johnpye 669 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
487     }else{
488     CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
489 johnpye 903 }*/
490 johnpye 669 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 johnpye 903 /**
508 johnpye 912 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
509 johnpye 909
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 johnpye 903 */
524 johnpye 909 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 johnpye 928 int var_yindex;
536 johnpye 903
537 johnpye 909 int *variables;
538     double *derivatives;
539     var_filter_t filter;
540     int count;
541    
542 johnpye 944 /* fprintf(stderr,"\n--------------\n"); */
543     /* CONSOLE_DEBUG("EVALUTING JACOBIAN..."); */
544 johnpye 909
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 johnpye 912 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
555 johnpye 909 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 johnpye 944 /* CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v)); */
565     /* for(i=0; i<NV_LENGTH_S(v); ++i){
566 johnpye 930 CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));
567 johnpye 944 }*/
568 johnpye 912
569 johnpye 909 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 johnpye 944 /* fprintf(stderr,"\n"); */
576 johnpye 928 relname = rel_make_name(blsys->system, *relptr);
577 johnpye 944 /* CONSOLE_DEBUG("RELATION %d '%s'",i,relname); */
578 johnpye 928 ASC_FREE(relname);
579    
580 johnpye 909 /* get derivatives for this particular relation */
581     status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
582 johnpye 944 /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
583 johnpye 928
584 johnpye 924 for(j=0;j<count;++j){
585 johnpye 928 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
586 johnpye 944 /* CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname); */
587 johnpye 917 ASC_FREE(varname);
588     }
589 johnpye 909
590     if(!status){
591 johnpye 944 /* CONSOLE_DEBUG("Derivatives for relation %d OK",i); */
592 johnpye 909 }else{
593 johnpye 928 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
594 johnpye 909 break;
595     }
596    
597 johnpye 928 /*
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 johnpye 909 Jv_i = 0;
604     for(j=0; j < count; ++j){
605 johnpye 910 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y); */
606 johnpye 909 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
607     if(varname){
608 johnpye 944 /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]); */
609 johnpye 909 ASC_FREE(varname);
610     }else{
611 johnpye 912 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
612 johnpye 909 }
613 johnpye 912
614 johnpye 928 var_yindex = blsys->y_id[variables[j]];
615 johnpye 944 /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
616 johnpye 928
617     if(var_yindex >= 0){
618 johnpye 944 /* CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
619 johnpye 928 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
620     }else{
621     var_yindex = -var_yindex-1;
622 johnpye 944 /* CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
623 johnpye 928 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
624     }
625 johnpye 909 }
626    
627     NV_Ith_S(Jv,i) = Jv_i;
628 johnpye 944 /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */
629 johnpye 910
630 johnpye 909 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 johnpye 912
642    
643    
644 johnpye 909 return is_error;
645     }
646    
647 johnpye 669 /*----------------------------------------------
648     ERROR REPORTING
649 johnpye 782 */
650 johnpye 669 /**
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