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

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