/[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 951 - (show annotations) (download) (as text)
Sun Nov 26 05:01:49 2006 UTC (13 years, 11 months ago) by johnpye
File MIME type: text/x-csrc
File size: 28436 byte(s)
Added 'GSMODIFIED' parameter to IDA implementation.
Working on getting the idadenx model to converge using the SPGMR solver.
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, DBL_MIN, 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","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