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

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