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

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