/[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 948 - (show annotations) (download) (as text)
Sun Nov 26 01:36:49 2006 UTC (15 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 27876 byte(s)
Fixed bug with integrator_ida_djex, and switched to using *solvers* varlist in integrator 'visit' routine (needs checking with LSODE!)
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_dense.h>
64 # ifndef IDA_SUCCESS
65 # error "Failed to include SUNDIALS IDA header file"
66 # endif
67 #endif
68
69 /*
70 for the benefit of build tools that didn't sniff the SUNDIALS version, we
71 assume version 2.2.x (and thence possible errors).
72 */
73 #ifndef SUNDIALS_VERSION_MINOR
74 # ifdef __GNUC__
75 # warning "GUESSING SUNDIALS VERSION 2.2"
76 # endif
77 # define SUNDIALS_VERSION_MINOR 2
78 #endif
79 #ifndef SUNDIALS_VERSION_MAJOR
80 # define SUNDIALS_VERSION_MAJOR 2
81 #endif
82
83 /* check that we've got what we expect now */
84 #ifndef ASC_IDA_H
85 # error "Failed to include ASCEND IDA header file"
86 #endif
87
88 /**
89 Struct containing any stuff that IDA needs that doesn't fit into the
90 common IntegratorSystem struct.
91 */
92 typedef struct{
93 struct rel_relation **rellist; /**< NULL terminated list of rels */
94 struct var_variable **varlist; /**< NULL terminated list of vars. ONLY USED FOR DEBUGGING -- get rid of it! */
95 int nrels;
96 int safeeval; /**< whether to pass the 'safe' flag to relman_eval */
97 } IntegratorIdaData;
98
99 /*-------------------------------------------------------------
100 FORWARD DECLS
101 */
102 /* residual function forward declaration */
103 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data);
104
105 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
106 , N_Vector v, N_Vector Jv, realtype c_j
107 , void *jac_data, N_Vector tmp1, N_Vector tmp2
108 );
109
110 /* error handler forward declaration */
111 void integrator_ida_error(int error_code
112 , const char *module, const char *function
113 , char *msg, void *eh_data
114 );
115
116 int integrator_ida_djex(long int Neq, realtype tt
117 , N_Vector yy, N_Vector yp, N_Vector rr
118 , realtype c_j, void *jac_data, DenseMat Jac
119 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
120 );
121
122 /*-------------------------------------------------------------
123 SETUP/TEARDOWN ROUTINES
124 */
125 void integrator_ida_create(IntegratorSystem *blsys){
126 CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
127 IntegratorIdaData *enginedata;
128 enginedata = ASC_NEW(IntegratorIdaData);
129 enginedata->rellist = NULL;
130 enginedata->varlist = NULL;
131 enginedata->safeeval = 1;
132 blsys->enginedata = (void *)enginedata;
133 integrator_ida_params_default(blsys);
134 }
135
136 void integrator_ida_free(void *enginedata){
137 CONSOLE_DEBUG("DELETING IDA ENGINE DATA");
138 IntegratorIdaData *d = (IntegratorIdaData *)enginedata;
139 /* note, we don't own the rellist, so don't need to free it */
140 ASC_FREE(d);
141 }
142
143 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *blsys){
144 IntegratorIdaData *d;
145 assert(blsys!=NULL);
146 assert(blsys->enginedata!=NULL);
147 assert(blsys->engine==INTEG_IDA);
148 d = ((IntegratorIdaData *)(blsys->enginedata));
149 assert(d->safeeval = 1);
150 return d;
151 }
152
153 /*-------------------------------------------------------------
154 PARAMETERS FOR IDA
155 */
156
157 enum ida_parameters{
158 IDA_PARAM_LINSOLVER
159 ,IDA_PARAM_AUTODIFF
160 ,IDA_PARAM_RTOL
161 ,IDA_PARAM_ATOL
162 ,IDA_PARAM_ATOLVECT
163 ,IDA_PARAMS_SIZE
164 };
165
166 /**
167 Here the full set of parameters is defined, along with upper/lower bounds,
168 etc. The values are stuck into the blsys->params structure.
169
170 @return 0 on success
171 */
172 int integrator_ida_params_default(IntegratorSystem *blsys){
173 asc_assert(blsys!=NULL);
174 asc_assert(blsys->engine==INTEG_IDA);
175 slv_parameters_t *p;
176 p = &(blsys->params);
177
178 slv_destroy_parms(p);
179
180 if(p->parms==NULL){
181 CONSOLE_DEBUG("params NULL");
182 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
183 if(p->parms==NULL)return -1;
184 p->dynamic_parms = 1;
185 }else{
186 CONSOLE_DEBUG("params not NULL");
187 }
188
189 /* reset the number of parameters to zero so that we can check it at the end */
190 p->num_parms = 0;
191
192 slv_param_bool(p,IDA_PARAM_AUTODIFF
193 ,(SlvParameterInitBool){{"autodiff"
194 ,"Use auto-diff?",1
195 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
196 }, TRUE}
197 );
198
199 slv_param_bool(p,IDA_PARAM_ATOLVECT
200 ,(SlvParameterInitBool){{"atolvect"
201 ,"Use 'ode_atol' values as specified?",1
202 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
203 " in the integration. If FALSE, a scalar absolute tolerance value"
204 " is shared by all variables. See IDA manual, section 5.5.1"
205 }, TRUE }
206 );
207
208 slv_param_real(p,IDA_PARAM_ATOL
209 ,(SlvParameterInitReal){{"atol"
210 ,"Scalar absolute error tolerance",1
211 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
212 " See IDA manual, section 5.5.1"
213 }, 1e-5, DBL_MIN, DBL_MAX }
214 );
215
216 slv_param_real(p,IDA_PARAM_RTOL
217 ,(SlvParameterInitReal){{"rtol"
218 ,"Scalar relative error tolerance",1
219 ,"Value of the scalar relative error tolerance."
220 " See IDA manual, section 5.5.1"
221 }, 1e-4, DBL_MIN, DBL_MAX }
222 );
223
224 slv_param_char(p,IDA_PARAM_LINSOLVER
225 ,(SlvParameterInitChar){{"linsolver"
226 ,"Linear solver",1
227 ,"See IDA manual, section 5.5.3."
228 }, "SPGMR"}, (char *[]){"DENSE","SPGMR",NULL}
229 );
230
231 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
232
233 CONSOLE_DEBUG("Created %d params", p->num_parms);
234
235 return 0;
236 }
237
238 /*-------------------------------------------------------------
239 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
240 */
241
242 /* return 1 on success */
243 int integrator_ida_solve(
244 IntegratorSystem *blsys
245 , unsigned long start_index
246 , unsigned long finish_index
247 ){
248 void *ida_mem;
249 int size, flag, t_index;
250 realtype t0, reltol, abstol, t, tret, tout1;
251 N_Vector y0, yp0, abstolvect, ypret, yret;
252 IntegratorIdaData *enginedata;
253 char *linsolver;
254
255 CONSOLE_DEBUG("STARTING IDA...");
256
257 enginedata = integrator_ida_enginedata(blsys);
258 CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
259
260 /* store reference to list of relations (in enginedata) */
261 enginedata->nrels = slv_get_num_solvers_rels(blsys->system);
262 enginedata->rellist = slv_get_solvers_rel_list(blsys->system);
263 enginedata->varlist = slv_get_solvers_var_list(blsys->system);
264 CONSOLE_DEBUG("Number of relations: %d",enginedata->nrels);
265 CONSOLE_DEBUG("Number of dependent vars: %ld",blsys->n_y);
266 size = blsys->n_y;
267
268 if(enginedata->nrels!=size){
269 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Integration problem is not square (%d rels, %d vars)", enginedata->nrels, size);
270 return 0; /* failure */
271 }
272
273 /* retrieve initial values from the system */
274
275 /** @TODO fix this, the starting time != first sample */
276 t0 = integrator_get_t(blsys);
277 CONSOLE_DEBUG("RETRIEVED t0 = %f",t0);
278
279 CONSOLE_DEBUG("RETRIEVING y0");
280
281 y0 = N_VNew_Serial(size);
282 integrator_get_y(blsys,NV_DATA_S(y0));
283
284 CONSOLE_DEBUG("RETRIEVING yp0");
285
286 yp0 = N_VNew_Serial(size);
287 integrator_get_ydot(blsys,NV_DATA_S(yp0));
288
289 N_VPrint_Serial(yp0);
290 CONSOLE_DEBUG("yp0 is at %p",&yp0);
291
292 /* create IDA object */
293 ida_mem = IDACreate();
294
295 /* relative error tolerance */
296 reltol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_RTOL);
297
298 /* allocate internal memory */
299 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_ATOLVECT)){
300 /* vector of absolute tolerances */
301 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
302 abstolvect = N_VNew_Serial(size);
303 integrator_get_atol(blsys,NV_DATA_S(abstolvect));
304
305 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV, reltol, abstolvect);
306 }else{
307 /* scalar absolute tolerance (one value for all) */
308 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
309 abstol = SLV_PARAM_REAL(&(blsys->params),IDA_PARAM_ATOL);
310 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS, reltol, &abstol);
311 }
312
313 if(flag==IDA_MEM_NULL){
314 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
315 return 0;
316 }else if(flag==IDA_MEM_FAIL){
317 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
318 return 0;
319 }else if(flag==IDA_ILL_INPUT){
320 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
321 return 0;
322 }/* else success */
323
324 /* set optional inputs... */
325 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *)blsys);
326 IDASetRdata(ida_mem, (void *)blsys);
327 IDASetMaxStep(ida_mem, integrator_get_maxstep(blsys));
328 IDASetInitStep(ida_mem, integrator_get_stepzero(blsys));
329 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(blsys));
330 /* there's no capability for setting *minimum* step size in IDA */
331
332
333 /* attach linear solver module, using the default value of maxl */
334 linsolver = SLV_PARAM_CHAR(&(blsys->params),IDA_PARAM_LINSOLVER);
335 CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
336 if(strcmp(linsolver,"SPGMR")==0){
337 CONSOLE_DEBUG("USING 'SCALED PRECONDITIONER GMRES' LINEAR SOLVER");
338 flag = IDASpgmr(ida_mem, 0);
339 if(flag==IDASPILS_MEM_NULL){
340 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
341 return 0;
342 }else if(flag==IDASPILS_MEM_FAIL){
343 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
344 return 0;
345 }/* else success */
346
347 /* assign the J*v function */
348 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
349 CONSOLE_DEBUG("USING AUTODIFF");
350 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex, (void *)blsys);
351 if(flag==IDASPILS_MEM_NULL){
352 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
353 return 0;
354 }else if(flag==IDASPILS_LMEM_NULL){
355 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
356 return 0;
357 }/* else success */
358 }else{
359 CONSOLE_DEBUG("USING NUMERICAL DIFF");
360 }
361 }else if(strcmp(linsolver,"DENSE")==0){
362 CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",size);
363 flag = IDADense(ida_mem, size);
364 switch(flag){
365 case IDADENSE_SUCCESS: break;
366 case IDADENSE_MEM_NULL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL"); return 0;
367 case IDADENSE_ILL_INPUT: ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module"); return 0;
368 case IDADENSE_MEM_FAIL: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE"); return 0;
369 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return"); return 0;
370 }
371
372 if(SLV_PARAM_BOOL(&(blsys->params),IDA_PARAM_AUTODIFF)){
373 CONSOLE_DEBUG("USING AUTODIFF");
374 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex, (void *)blsys);
375 switch(flag){
376 case IDADENSE_SUCCESS: break;
377 default: ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn"); return 0;
378 }
379 }else{
380 CONSOLE_DEBUG("USING NUMERICAL DIFF");
381 }
382 }else{
383 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
384 return 0;
385 }
386
387 /* set linear solver optional inputs...
388
389 ...nothing here at the moment...
390
391 */
392
393 /* correct initial values, given derivatives */
394 blsys->currentstep=0;
395 t_index=start_index;
396 tout1 = samplelist_get(blsys->samples, t_index);
397
398 /* CONSOLE_DEBUG("Giving t value %f to IDACalcIC", tout1);*/
399
400 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==3
401 /* note the new API from version 2.3 and onwards */
402 flag = IDACalcIC(ida_mem, IDA_Y_INIT, tout1);
403 #else
404 flag = IDACalcIC(ida_mem, t0, y0, yp0, IDA_Y_INIT, tout1);
405 #endif
406
407 if(flag!=IDA_SUCCESS){
408 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to solve initial values (IDACalcIC)");
409 return 0;
410 }/* else success */
411
412 CONSOLE_DEBUG("INITIAL CONDITIONS SOLVED :-)");
413
414 /* optionally, specify ROO-FINDING PROBLEM */
415
416 /* -- set up the IntegratorReporter */
417 integrator_output_init(blsys);
418
419 /* -- store the initial values of all the stuff */
420 integrator_output_write(blsys);
421 integrator_output_write_obs(blsys);
422
423 /* specify where the returned values should be stored */
424 yret = y0;
425 ypret = yp0;
426
427 /* advance solution in time, return values as yret and derivatives as ypret */
428 blsys->currentstep=1;
429 for(t_index=start_index+1;t_index <= finish_index;++t_index, ++blsys->currentstep){
430 t = samplelist_get(blsys->samples, t_index);
431
432 /* CONSOLE_DEBUG("SOLVING UP TO t = %f", t); */
433
434 flag = IDASolve(ida_mem, t, &tret, yret, ypret, IDA_NORMAL);
435
436 /* pass the values of everything back to the compiler */
437 integrator_set_t(blsys, (double)tret);
438 integrator_set_y(blsys, NV_DATA_S(yret));
439 integrator_set_ydot(blsys, NV_DATA_S(ypret));
440
441 if(flag<0){
442 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", t, flag);
443 break;
444 }
445
446 /* -- do something so that blsys knows the values of tret, yret and ypret */
447
448 /* -- store the current values of all the stuff */
449 integrator_output_write(blsys);
450 integrator_output_write_obs(blsys);
451 }
452
453 /* -- close the IntegratorReporter */
454 integrator_output_close(blsys);
455
456 if(flag < 0){
457 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", t);
458 return 0;
459 }
460
461 /* get optional outputs */
462
463 /* free solution memory */
464 N_VDestroy_Serial(yret);
465 N_VDestroy_Serial(ypret);
466
467 /* free solver memory */
468 IDAFree(ida_mem);
469
470 /* all done */
471 return 1;
472 }
473
474 /*--------------------------------------------------
475 RESIDUALS AND JACOBIAN
476 */
477 /**
478 Function to evaluate system residuals, in the form required for IDA.
479
480 Given tt, yy and yp, we need to evaluate and return rr.
481
482 @param tt current value of indep variable (time)
483 @param yy current values of dependent variable vector
484 @param yp current values of derivatives of dependent variables
485 @param rr the output residual vector (is we're returning data to)
486 @param res_data pointer to our stuff (blsys in this case).
487
488 @return 0 on success, positive on recoverable error, and
489 negative on unrecoverable error.
490 */
491 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
492 IntegratorSystem *blsys;
493 IntegratorIdaData *enginedata;
494 int i, calc_ok, is_error;
495 struct rel_relation** relptr;
496 double resid;
497 char *relname;
498
499 blsys = (IntegratorSystem *)res_data;
500 enginedata = integrator_ida_enginedata(blsys);
501
502 /* fprintf(stderr,"\n\n"); */
503 /* CONSOLE_DEBUG("ABOUT TO EVALUTE RESIDUALS..."); */
504
505 /* pass the values of everything back to the compiler */
506 integrator_set_t(blsys, (double)tt);
507 integrator_set_y(blsys, NV_DATA_S(yy));
508 integrator_set_ydot(blsys, NV_DATA_S(yp));
509
510 /* revaluate the system residuals using the new data */
511 is_error = 0;
512 relptr = enginedata->rellist;
513
514 /* CONSOLE_DEBUG("IDA requests residuals of length %lu",NV_LENGTH_S(rr)); */
515 if(NV_LENGTH_S(rr)!=enginedata->nrels){
516 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
517 return -1; /* unrecoverable */
518 }
519
520 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
521 if (setjmp(g_fpe_env)==0) {
522 for(i=0, relptr = enginedata->rellist;
523 i< enginedata->nrels && relptr != NULL;
524 ++i, ++relptr
525 ){
526 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
527
528 relname = rel_make_name(blsys->system, *relptr);
529 /* if(calc_ok){
530 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f",i,relname,resid);
531 }else{
532 CONSOLE_DEBUG("residual[%d:\"%s\"] = %f (ERROR)",i,relname,resid);
533 }*/
534 ASC_FREE(relname);
535
536 NV_Ith_S(rr,i) = resid;
537 if(!calc_ok){
538 /* presumable some output already made? */
539 is_error = 1;
540 }
541 }
542 }else{
543 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
544 }
545 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
546
547 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
548 return is_error;
549 }
550
551 /**
552 Dense Jacobian evaluation. Only suitable for small problems!
553 */
554 int integrator_ida_djex(long int Neq, realtype tt
555 , N_Vector yy, N_Vector yp, N_Vector rr
556 , realtype c_j, void *jac_data, DenseMat Jac
557 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
558 ){
559 IntegratorSystem *blsys;
560 IntegratorIdaData *enginedata;
561 char *relname, *varname;
562 int status;
563 struct rel_relation **relptr;
564 int i;
565 var_filter_t filter = {VAR_SVAR, VAR_SVAR};
566 double *derivatives;
567 int *variables;
568 int count, j, var_yindex;
569
570 blsys = (IntegratorSystem *)jac_data;
571 enginedata = integrator_ida_enginedata(blsys);
572
573 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
574 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
575 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
576
577 /* pass the values of everything back to the compiler */
578 integrator_set_t(blsys, (double)tt);
579 integrator_set_y(blsys, NV_DATA_S(yy));
580 integrator_set_ydot(blsys, NV_DATA_S(yp));
581
582 /* print vars */
583 for(i=0; i < blsys->n_y; ++i){
584 varname = var_make_name(blsys->system, blsys->y[i]);
585 CONSOLE_DEBUG("%s = %f = %f",varname,NV_Ith_S(yy,i),var_value(blsys->y[i]));
586 ASC_FREE(varname);
587 }
588
589 /* print derivatives */
590 for(i=0; i < blsys->n_y; ++i){
591 if(blsys->ydot[i]){
592 varname = var_make_name(blsys->system, blsys->ydot[i]);
593 CONSOLE_DEBUG("%s = %f =%f",varname,NV_Ith_S(yp,i),var_value(blsys->ydot[i]));
594 ASC_FREE(varname);
595 }else{
596 varname = var_make_name(blsys->system, blsys->y[i]);
597 CONSOLE_DEBUG("diff(%s) = %f",varname,NV_Ith_S(yp,i));
598 ASC_FREE(varname);
599 }
600 }
601
602 /* print step size */
603 CONSOLE_DEBUG("<c_j> = %f",c_j);
604
605 for(i=0, relptr = enginedata->rellist;
606 i< enginedata->nrels && relptr != NULL;
607 ++i, ++relptr
608 ){
609 relname = rel_make_name(blsys->system, *relptr);
610 CONSOLE_DEBUG("%d: '%s'",i,relname);
611 }
612
613 /* build up the dense jacobian matrix... */
614 status = 0;
615 for(i=0, relptr = enginedata->rellist;
616 i< enginedata->nrels && relptr != NULL;
617 ++i, ++relptr
618 ){
619 /* get derivatives for this particular relation */
620 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
621
622 if(status){
623 relname = rel_make_name(blsys->system, *relptr);
624 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
625 ASC_FREE(relname);
626 break;
627 }
628
629 /* output what's going on here ... */
630 relname = rel_make_name(blsys->system, *relptr);
631 CONSOLE_DEBUG("RELATION %d '%s'",i,relname);
632 fprintf(stderr,"%d: '%s': ",i,relname);
633 ASC_FREE(relname);
634 for(j=0;j<count;++j){
635 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
636 var_yindex = blsys->y_id[variables[j]];
637 if(var_yindex >=0){
638 fprintf(stderr," var[%d]='%s'=y[%d]",variables[j],varname,var_yindex);
639 }else{
640 fprintf(stderr," var[%d]='%s'=ydot[%d]",variables[j],varname,-var_yindex-1);
641 }
642 ASC_FREE(varname);
643 }
644 fprintf(stderr,"\n");
645
646 /* zero the Jacobian row */
647 for(j=0; j < enginedata->nrels; ++j){
648 DENSE_ELEM(Jac,i,j) = 0;
649 }
650
651 /* insert values into the Jacobian row in appropriate spots */
652 for(j=0; j < count; ++j){
653 var_yindex = blsys->y_id[variables[j]];
654 if(var_yindex >= 0){
655 asc_assert(blsys->y[var_yindex]==enginedata->varlist[variables[j]]);
656 DENSE_ELEM(Jac,i,var_yindex) += derivatives[j];
657 }else{
658 asc_assert(blsys->ydot[-var_yindex-1]==enginedata->varlist[variables[j]]);
659 DENSE_ELEM(Jac,i,-var_yindex-1) += derivatives[j] * c_j;
660 }
661 }
662 }
663
664 CONSOLE_DEBUG("PRINTING JAC");
665 fprintf(stderr,"\t");
666 for(j=0; j < blsys->n_y; ++j){
667 if(j)fprintf(stderr,"\t");
668 varname = var_make_name(blsys->system,blsys->y[j]);
669 fprintf(stderr,"%11s",varname);
670 ASC_FREE(varname);
671 }
672 fprintf(stderr,"\n");
673 for(i=0; i < enginedata->nrels; ++i){
674 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
675 fprintf(stderr,"%s\t",relname);
676 ASC_FREE(relname);
677
678 for(j=0; j < blsys->n_y; ++j){
679 if(j)fprintf(stderr,"\t");
680 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
681 }
682 fprintf(stderr,"\n");
683 }
684
685 #if 0
686 double *yval = NV_DATA_S(yy);
687 /* AWFUL HACK! In an attempt to prove that 'IDADense' works as expected,
688 I've pulled the jacobians straight from idadenx.c (an IDA example) */
689 #define A(REL,I,J,DESIRED) if(fabs((DENSE_ELEM(Jac,I,J) - DESIRED)/DESIRED) > 1e-5){ \
690 CONSOLE_DEBUG("Value for ('%s'=%d,%d) got %e, expected %e",REL,I,J,DENSE_ELEM(Jac,I,J),DESIRED); \
691 status=1; \
692 }else{ \
693 CONSOLE_DEBUG("Value for (%d,%d) is OK",I,J); \
694 }
695 char *eqs[] = {"eq1","eq2","eq3"};
696 for(i=0;i<enginedata->nrels;++i){
697 relname = rel_make_name(blsys->system, enginedata->rellist[i]);
698 if(strcmp(relname,"eq1")==0){
699 A(relname,i,0, -0.04 - c_j );
700 A(relname,i,1, 1.0e4*yval[2] );
701 A(relname,i,2, 1.0e4*yval[1] );
702 }else if(strcmp(relname,"eq2")==0){
703 A(relname,i,0, 0.04);
704 A(relname,i,1, RCONST(-1.0e4)*yval[2] - RCONST(6.0e7)*yval[1] - c_j);
705 A(relname,i,2, RCONST(-1.0e4)*yval[1]);
706 }else if(strcmp(relname,"eq3")==0){
707 A(relname,i,0, -1.0);
708 A(relname,i,1, -1.0);
709 A(relname,i,2, -1.0);
710 }
711 }
712 #undef A
713 #endif
714
715 if(status){
716 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
717 return 1;
718 }
719
720 return 0;
721 }
722
723 /**
724 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
725
726 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
727
728 @param tt current value of the independent variable (time, t)
729 @param yy current value of the dependent variable vector, y(t).
730 @param yp current value of y'(t).
731 @param rr current value of the residual vector F(t, y, y').
732 @param v the vector by which the Jacobian must be multiplied to the right.
733 @param Jv the output vector computed
734 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
735 @param jac_data pointer to our stuff (blsys in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
736 @param tmp1 @see tmp2
737 @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.
738 @return 0 on success
739 */
740 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
741 , N_Vector v, N_Vector Jv, realtype c_j
742 , void *jac_data, N_Vector tmp1, N_Vector tmp2
743 ){
744 IntegratorSystem *blsys;
745 IntegratorIdaData *enginedata;
746 int i, j, is_error=0;
747 struct rel_relation** relptr;
748 char *relname, *varname;
749 int status;
750 double Jv_i;
751 int var_yindex;
752
753 int *variables;
754 double *derivatives;
755 var_filter_t filter;
756 int count;
757
758 /* fprintf(stderr,"\n--------------\n"); */
759 /* CONSOLE_DEBUG("EVALUTING JACOBIAN..."); */
760
761 blsys = (IntegratorSystem *)jac_data;
762 enginedata = integrator_ida_enginedata(blsys);
763
764 /* pass the values of everything back to the compiler */
765 integrator_set_t(blsys, (double)tt);
766 integrator_set_y(blsys, NV_DATA_S(yy));
767 integrator_set_ydot(blsys, NV_DATA_S(yp));
768 /* no real use for residuals (rr) here, I don't think? */
769
770 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
771 variables = ASC_NEW_ARRAY(int, NV_LENGTH_S(yy) * 2);
772 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
773
774 /* evaluate the derivatives... */
775 /* J = dG_dy = dF_dy + alpha * dF_dyp */
776
777 filter.matchbits = VAR_SVAR;
778 filter.matchvalue = VAR_SVAR;
779
780 /* CONSOLE_DEBUG("PRINTING VALUES OF 'v' VECTOR (length %ld)",NV_LENGTH_S(v)); */
781 /* for(i=0; i<NV_LENGTH_S(v); ++i){
782 CONSOLE_DEBUG("v[%d] = %f",i,NV_Ith_S(v,i));
783 }*/
784
785 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
786 if (setjmp(g_fpe_env)==0) {
787 for(i=0, relptr = enginedata->rellist;
788 i< enginedata->nrels && relptr != NULL;
789 ++i, ++relptr
790 ){
791 /* fprintf(stderr,"\n"); */
792 relname = rel_make_name(blsys->system, *relptr);
793 /* CONSOLE_DEBUG("RELATION %d '%s'",i,relname); */
794 ASC_FREE(relname);
795
796 /* get derivatives for this particular relation */
797 status = relman_diff2(*relptr, &filter, derivatives, variables, &count, enginedata->safeeval);
798 /* CONSOLE_DEBUG("Got derivatives against %d matching variables", count); */
799
800 for(j=0;j<count;++j){
801 /* varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]);
802 CONSOLE_DEBUG("derivatives[%d] = %f (variable %d, '%s')",j,derivatives[j],variables[j],varname);
803 ASC_FREE(varname); */
804 }
805
806 if(!status){
807 /* CONSOLE_DEBUG("Derivatives for relation %d OK",i); */
808 }else{
809 CONSOLE_DEBUG("ERROR calculating derivatives for relation %d",i);
810 break;
811 }
812
813 /*
814 Now we have the derivatives wrt each alg/diff variable in the
815 present equation. variables[] points into the varlist. need
816 a mapping from the varlist to the y and ydot lists.
817 */
818
819 Jv_i = 0;
820 for(j=0; j < count; ++j){
821 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], blsys->n_y);
822 varname = var_make_name(blsys->system, enginedata->varlist[variables[j]]); */
823 if(varname){
824 /* CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
825 ASC_FREE(varname); */
826 }else{
827 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
828 }
829
830 var_yindex = blsys->y_id[variables[j] - 1];
831 /* CONSOLE_DEBUG("j = %d: variables[j] = %d, y_id = %d",j,variables[j],var_yindex); */
832
833 if(var_yindex >= 0){
834 /* CONSOLE_DEBUG("j = %d: algebraic, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
835 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex);
836 }else{
837 var_yindex = -var_yindex-1;
838 /* CONSOLE_DEBUG("j = %d: differential, deriv[j] = %f, v[%d] = %f",j,derivatives[j], var_yindex, NV_Ith_S(v,var_yindex)); */
839 Jv_i += derivatives[j] * NV_Ith_S(v,var_yindex) / c_j;
840 }
841 }
842
843 NV_Ith_S(Jv,i) = Jv_i;
844 /* CONSOLE_DEBUG("(J*v)[%d] = %f", i, Jv_i); */
845 }
846 if(status){
847 CONSOLE_DEBUG("JVEX ERROR RESULT");
848 /* presumably some error_reporter will already have been made*/
849 is_error = 1;
850 }
851 }else{
852 CONSOLE_DEBUG("FLOATING POINT ERROR WITH i=%d",i);
853 }
854 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
855
856 if(is_error)CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
857
858
859
860 return is_error;
861 }
862
863 /*----------------------------------------------
864 ERROR REPORTING
865 */
866 /**
867 Error message reporter function to be passed to IDA. All error messages
868 will trigger a call to this function, so we should find everything
869 appearing on the console (in the case of Tcl/Tk) or in the errors/warnings
870 panel (in the case of PyGTK).
871 */
872 void integrator_ida_error(int error_code
873 , const char *module, const char *function
874 , char *msg, void *eh_data
875 ){
876 IntegratorSystem *blsys;
877 error_severity_t sev;
878
879 /* cast back the IntegratorSystem, just in case we need it */
880 blsys = (IntegratorSystem *)eh_data;
881
882 /* severity depends on the sign of the error_code value */
883 if(error_code <= 0){
884 sev = ASC_PROG_ERR;
885 }else{
886 sev = ASC_PROG_WARNING;
887 }
888
889 /* use our all-purpose error reporting to get stuff back to the GUI */
890 error_reporter(sev,module,0,function,"%s (error %d)",msg,error_code);
891 }

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