/[ascend]/trunk/solvers/ida/idacalc.c
ViewVC logotype

Contents of /trunk/solvers/ida/idacalc.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2555 - (show annotations) (download) (as text)
Fri Feb 3 00:04:07 2012 UTC (10 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 19457 byte(s)
Fix issue with solvers not being loaded on Windows, it was a PATH_MAX limit issue (fixes bug 501).
Fix minor memory leak in idacalc.c.
1 /* ASCEND modelling environment
2 Copyright (C) 2006-2011 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 Calculation of residuals, jacobian, etc, for the ASCEND wrapping of IDA.
21
22 @see http://www.llnl.gov/casc/sundials/
23 *//*
24 by John Pye, May 2006
25 */
26
27 #define _GNU_SOURCE
28
29 #include "idacalc.h"
30 #include "idatypes.h"
31 #include "idalinear.h"
32 #include "idaanalyse.h"
33
34 #include <signal.h>
35 #include <setjmp.h>
36 #include <fenv.h>
37 #include <math.h>
38
39 #ifdef ASC_WITH_MMIO
40 # include <mmio.h>
41 #endif
42
43 #include <ascend/general/platform.h>
44 #include <ascend/utilities/error.h>
45 #include <ascend/utilities/ascSignal.h>
46 #include <ascend/general/panic.h>
47 #include <ascend/compiler/instance_enum.h>
48
49 #include <ascend/system/slv_client.h>
50 #include <ascend/system/relman.h>
51 #include <ascend/system/block.h>
52 #include <ascend/system/slv_stdcalls.h>
53 #include <ascend/system/jacobian.h>
54 #include <ascend/system/bndman.h>
55
56 #include <ascend/utilities/config.h>
57 #include <ascend/integrator/integrator.h>
58
59
60 /* #define FEX_DEBUG */
61 #define JEX_DEBUG
62 /* #define DJEX_DEBUG */
63 /* #define ROOT_DEBUG */
64
65 /*--------------------------------------------------
66 RESIDUALS AND JACOBIAN AND IDAROOTFN
67 */
68
69 #if 0
70 typedef void (SignalHandlerFn)(int);
71 SignalHandlerFn integrator_ida_sig;
72 SignalHandlerFn *integrator_ida_sig_old;
73 jmp_buf integrator_ida_jmp_buf;
74 fenv_t integrator_ida_fenv_old;
75
76
77 void integrator_ida_write_feinfo(){
78 int f;
79 f = fegetexcept();
80 CONSOLE_DEBUG("Locating nature of exception...");
81 if(f & FE_DIVBYZERO)ERROR_REPORTER_HERE(ASC_PROG_ERR,"DIV BY ZERO");
82 if(f & FE_INEXACT)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INEXACT");
83 if(f & FE_INVALID)ERROR_REPORTER_HERE(ASC_PROG_ERR,"INVALID");
84 if(f & FE_OVERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"OVERFLOW");
85 if(f & FE_UNDERFLOW)ERROR_REPORTER_HERE(ASC_PROG_ERR,"UNDERFLOW");
86 if(f==0)ERROR_REPORTER_HERE(ASC_PROG_ERR,"FLAGS ARE CLEAR?!?");
87 }
88
89 void integrator_ida_sig(int sig){
90 /* the wrong signal: rethrow to the default handler */
91 if(sig!=SIGFPE){
92 signal(SIGFPE,SIG_DFL);
93 raise(sig);
94 }
95 integrator_ida_write_feinfo();
96 CONSOLE_DEBUG("Caught SIGFPE=%d (in signal handler). Jumping to...",sig);
97 longjmp(integrator_ida_jmp_buf,sig);
98 }
99 #endif
100
101 /**
102 Function to evaluate system residuals, in the form required for IDA.
103
104 Given tt, yy and yp, we need to evaluate and return rr.
105
106 @param tt current value of indep variable (time)
107 @param yy current values of dependent variable vector
108 @param yp current values of derivatives of dependent variables
109 @param rr the output residual vector (is we're returning data to)
110 @param res_data pointer to our stuff (integ in this case).
111
112 @return 0 on success, positive on recoverable error, and
113 negative on unrecoverable error.
114 */
115 int integrator_ida_fex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *res_data){
116 IntegratorSystem *integ;
117 IntegratorIdaData *enginedata;
118 int i, calc_ok, is_error;
119 struct rel_relation** relptr;
120 double resid;
121 char *relname;
122 #ifdef FEX_DEBUG
123 char *varname;
124 char diffname[30];
125 #endif
126
127 integ = (IntegratorSystem *)res_data;
128 enginedata = integrator_ida_enginedata(integ);
129
130 #ifdef FEX_DEBUG
131 /* fprintf(stderr,"\n\n"); */
132 CONSOLE_DEBUG("EVALUTE RESIDUALS...");
133 #endif
134
135 if(NV_LENGTH_S(rr)!=enginedata->nrels){
136 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid residuals nrels!=length(rr)");
137 return -1; /* unrecoverable */
138 }
139
140 /* pass the values of everything back to the compiler */
141 integrator_set_t(integ, (double)tt);
142 integrator_set_y(integ, NV_DATA_S(yy));
143 integrator_set_ydot(integ, NV_DATA_S(yp));
144
145 /* perform bounds checking on all variables */
146 if(slv_check_bounds(integ->system, 0, -1, NULL)){
147 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
148 return 1;
149 }
150
151 /* evaluate each residual in the rellist */
152 is_error = 0;
153 relptr = enginedata->rellist;
154
155
156 #ifdef ASC_SIGNAL_TRAPS
157 if(enginedata->safeeval){
158 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
159 }else{
160 # ifdef FEX_DEBUG
161 ERROR_REPORTER_HERE(ASC_PROG_ERR,"SETTING TO CATCH SIGFPE...");
162 # endif
163 Asc_SignalHandlerPushDefault(SIGFPE);
164 }
165
166 if (SETJMP(g_fpe_env)==0) {
167 #endif
168
169
170 for(i=0, relptr = enginedata->rellist;
171 i< enginedata->nrels && relptr != NULL;
172 ++i, ++relptr
173 ){
174 resid = relman_eval(*relptr, &calc_ok, enginedata->safeeval);
175
176 NV_Ith_S(rr,i) = resid;
177 if(!calc_ok){
178 relname = rel_make_name(integ->system, *relptr);
179 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
180 ASC_FREE(relname);
181 /* presumable some output already made? */
182 is_error = 1;
183 }/*else{
184 CONSOLE_DEBUG("Calc OK");
185 }*/
186 }
187
188 if(!is_error){
189 for(i=0;i< enginedata->nrels; ++i){
190 if(isnan(NV_Ith_S(rr,i))){
191 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in residual %d",i);
192 is_error=1;
193 }
194 }
195 #ifdef FEX_DEBUG
196 if(!is_error){
197 CONSOLE_DEBUG("No NAN detected");
198 }
199 #endif
200 }
201
202 #ifdef ASC_SIGNAL_TRAPS
203 }else{
204 relname = rel_make_name(integ->system, *relptr);
205 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
206 ASC_FREE(relname);
207 is_error = 1;
208 }
209
210 if(enginedata->safeeval){
211 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
212 }else{
213 Asc_SignalHandlerPopDefault(SIGFPE);
214 }
215 #endif
216
217
218 #ifdef FEX_DEBUG
219 /* output residuals to console */
220 CONSOLE_DEBUG("RESIDUAL OUTPUT");
221 fprintf(stderr,"index\t%25s\t%25s\t%s\n","y","ydot","resid");
222 for(i=0; i<integ->n_y; ++i){
223 varname = var_make_name(integ->system,integ->y[i]);
224 fprintf(stderr,"%d\t%15s=%10f\t",i,varname,NV_Ith_S(yy,i));
225 if(integ->ydot[i]){
226 varname = var_make_name(integ->system,integ->ydot[i]);
227 fprintf(stderr,"%15s=%10f\t",varname,NV_Ith_S(yp,i));
228 }else{
229 snprintf(diffname,99,"diff(%s)",varname);
230 fprintf(stderr,"%15s=%10f\t",diffname,NV_Ith_S(yp,i));
231 }
232 ASC_FREE(varname);
233 relname = rel_make_name(integ->system,enginedata->rellist[i]);
234 fprintf(stderr,"'%s'=%f (%p)\n",relname,NV_Ith_S(rr,i),enginedata->rellist[i]);
235 }
236 #endif
237
238 if(is_error){
239 return 1;
240 }
241
242 #ifdef FEX_DEBUG
243 CONSOLE_DEBUG("RESIDUAL OK");
244 #endif
245 return 0;
246 }
247
248 /**
249 Dense Jacobian evaluation. Only suitable for small problems!
250 Has been seen working for problems up to around 2000 vars, FWIW.
251 */
252 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
253 int integrator_ida_djex(int Neq, realtype tt, realtype c_j
254 , N_Vector yy, N_Vector yp, N_Vector rr
255 , IDA_MTX_T Jac, void *jac_data
256 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
257 ){
258 #else
259 int integrator_ida_djex(long int Neq, realtype tt
260 , N_Vector yy, N_Vector yp, N_Vector rr
261 , realtype c_j, void *jac_data, IDA_MTX_T Jac
262 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
263 ){
264 #endif
265 IntegratorSystem *integ;
266 IntegratorIdaData *enginedata;
267 char *relname;
268 #ifdef DJEX_DEBUG
269 struct var_variable **varlist;
270 char *varname;
271 #endif
272 struct rel_relation **relptr;
273 int i;
274 double *derivatives;
275 struct var_variable **variables;
276 int count, j;
277 int status, is_error = 0;
278
279 integ = (IntegratorSystem *)jac_data;
280 enginedata = integrator_ida_enginedata(integ);
281
282 /* allocate space for returns from relman_diff3 */
283 /** @TODO instead, we should use 'tmp1' and 'tmp2' here... */
284 variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
285 derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
286
287 /* pass the values of everything back to the compiler */
288 integrator_set_t(integ, (double)tt);
289 integrator_set_y(integ, NV_DATA_S(yy));
290 integrator_set_ydot(integ, NV_DATA_S(yp));
291
292 /* perform bounds checking on all variables */
293 if(slv_check_bounds(integ->system, 0, -1, NULL)){
294 /* ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Variable(s) out of bounds"); */
295 return 1;
296 }
297
298 #ifdef DJEX_DEBUG
299 varlist = slv_get_solvers_var_list(integ->system);
300
301 /* print vars */
302 for(i=0; i < integ->n_y; ++i){
303 varname = var_make_name(integ->system, integ->y[i]);
304 CONSOLE_DEBUG("%s = %f",varname,NV_Ith_S(yy,i));
305 asc_assert(NV_Ith_S(yy,i) == var_value(integ->y[i]));
306 ASC_FREE(varname);
307 }
308
309 /* print derivatives */
310 for(i=0; i < integ->n_y; ++i){
311 if(integ->ydot[i]){
312 varname = var_make_name(integ->system, integ->ydot[i]);
313 CONSOLE_DEBUG("%s = %f =%g",varname,NV_Ith_S(yp,i),var_value(integ->ydot[i]));
314 ASC_FREE(varname);
315 }else{
316 varname = var_make_name(integ->system, integ->y[i]);
317 CONSOLE_DEBUG("diff(%s) = %g",varname,NV_Ith_S(yp,i));
318 ASC_FREE(varname);
319 }
320 }
321
322 /* print step size */
323 CONSOLE_DEBUG("<c_j> = %g",c_j);
324 #endif
325
326 /* build up the dense jacobian matrix... */
327 status = 0;
328 for(i=0, relptr = enginedata->rellist;
329 i< enginedata->nrels && relptr != NULL;
330 ++i, ++relptr
331 ){
332 /* get derivatives for this particular relation */
333 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
334
335 if(status){
336 relname = rel_make_name(integ->system, *relptr);
337 CONSOLE_DEBUG("ERROR calculating derivatives for relation '%s'",relname);
338 ASC_FREE(relname);
339 is_error = 1;
340 break;
341 }
342
343 /* output what's going on here ... */
344 #ifdef DJEX_DEBUG
345 relname = rel_make_name(integ->system, *relptr);
346 fprintf(stderr,"%d: '%s': ",i,relname);
347 for(j=0;j<count;++j){
348 varname = var_make_name(integ->system, variables[j]);
349 if(var_deriv(variables[j])){
350 fprintf(stderr," '%s'=",varname);
351 fprintf(stderr,"ydot[%d]",integrator_ida_diffindex(integ,variables[j]));
352 }else{
353 fprintf(stderr," '%s'=y[%d]",varname,var_sindex(variables[j]));
354 }
355 ASC_FREE(varname);
356 }
357 /* relname is freed further down */
358 fprintf(stderr,"\n");
359 #endif
360
361 /* insert values into the Jacobian row in appropriate spots (can assume Jac starts with zeros -- IDA manual) */
362 for(j=0; j < count; ++j){
363 #ifdef DJEX_DEBUG
364 varname = var_make_name(integ->system,variables[j]);
365 fprintf(stderr,"d(%s)/d(%s) = %g",relname,varname,derivatives[j]);
366 ASC_FREE(varname);
367 #endif
368 if(!var_deriv(variables[j])){
369 #ifdef DJEX_DEBUG
370 fprintf(stderr," --> J[%d,%d] += %g\n", i,j,derivatives[j]);
371 asc_assert(var_sindex(variables[j]) >= 0);
372 ASC_ASSERT_LT(var_sindex(variables[j]) , Neq);
373 #endif
374 DENSE_ELEM(Jac,i,var_sindex(variables[j])) += derivatives[j];
375 }else{
376 DENSE_ELEM(Jac,i,integrator_ida_diffindex(integ,variables[j])) += derivatives[j] * c_j;
377 #ifdef DJEX_DEBUG
378 fprintf(stderr," --> * c_j --> J[%d,%d] += %g\n", i,j,derivatives[j] * c_j);
379 #endif
380 }
381 }
382 }
383
384 #ifdef DJEX_DEBUG
385 ASC_FREE(relname);
386 CONSOLE_DEBUG("PRINTING JAC");
387 fprintf(stderr,"\t");
388 for(j=0; j < integ->n_y; ++j){
389 if(j)fprintf(stderr,"\t");
390 varname = var_make_name(integ->system,integ->y[j]);
391 fprintf(stderr,"%11s",varname);
392 ASC_FREE(varname);
393 }
394 fprintf(stderr,"\n");
395 for(i=0; i < enginedata->nrels; ++i){
396 relname = rel_make_name(integ->system, enginedata->rellist[i]);
397 fprintf(stderr,"%s\t",relname);
398 ASC_FREE(relname);
399
400 for(j=0; j < integ->n_y; ++j){
401 if(j)fprintf(stderr,"\t");
402 fprintf(stderr,"%11.2e",DENSE_ELEM(Jac,i,j));
403 }
404 fprintf(stderr,"\n");
405 }
406 #endif
407
408 /* test for NANs */
409 if(!is_error){
410 for(i=0;i< enginedata->nrels; ++i){
411 for(j=0;j<integ->n_y;++j){
412 if(isnan(DENSE_ELEM(Jac,i,j))){
413 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NAN detected in jacobian J[%d,%d]",i,j);
414 is_error=1;
415 }
416 }
417 }
418 #ifdef DJEX_DEBUG
419 if(!is_error){
420 CONSOLE_DEBUG("No NAN detected");
421 }
422 #endif
423 }
424
425 /* if(integrator_ida_check_diffindex(integ)){
426 is_error = 1;
427 }*/
428
429 ASC_FREE(variables);
430 ASC_FREE(derivatives);
431
432 if(is_error){
433 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There were derivative evaluation errors in the dense jacobian");
434 return 1;
435 }
436
437 #ifdef DJEX_DEBUG
438 CONSOLE_DEBUG("DJEX RETURNING 0");
439 /* ASC_PANIC("Quitting"); */
440 #endif
441 return 0;
442 }
443
444 /**
445 Function to evaluate the product J*v, in the form required for IDA (see IDASpilsSetJacTimesVecFn)
446
447 Given tt, yy, yp, rr and v, we need to evaluate and return Jv.
448
449 @param tt current value of the independent variable (time, t)
450 @param yy current value of the dependent variable vector, y(t).
451 @param yp current value of y'(t).
452 @param rr current value of the residual vector F(t, y, y').
453 @param v the vector by which the Jacobian must be multiplied to the right.
454 @param Jv the output vector computed
455 @param c_j the scalar in the system Jacobian, proportional to the inverse of the step size ($ \alpha$ in Eq. (3.5) ).
456 @param jac_data pointer to our stuff (integ in this case, passed into IDA via IDASp*SetJacTimesVecFn.)
457 @param tmp1 @see tmp2
458 @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.
459 @return 0 on success
460 */
461 int integrator_ida_jvex(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr
462 , N_Vector v, N_Vector Jv, realtype c_j
463 , void *jac_data, N_Vector tmp1, N_Vector tmp2
464 ){
465 IntegratorSystem *integ;
466 IntegratorIdaData *enginedata;
467 int i, j, is_error=0;
468 struct rel_relation** relptr = 0;
469 char *relname;
470 int status;
471 double Jv_i;
472
473 struct var_variable **variables;
474 double *derivatives;
475 int count;
476 struct var_variable **varlist;
477 #ifdef JEX_DEBUG
478
479 CONSOLE_DEBUG("EVALUATING JACOBIAN...");
480 #endif
481
482 integ = (IntegratorSystem *)jac_data;
483 enginedata = integrator_ida_enginedata(integ);
484 varlist = slv_get_solvers_var_list(integ->system);
485
486 /* pass the values of everything back to the compiler */
487 integrator_set_t(integ, (double)tt);
488 integrator_set_y(integ, NV_DATA_S(yy));
489 integrator_set_ydot(integ, NV_DATA_S(yp));
490 /* no real use for residuals (rr) here, I don't think? */
491
492 /* allocate space for returns from relman_diff2: we *should* be able to use 'tmp1' and 'tmp2' here... */
493
494 i = NV_LENGTH_S(yy) * 2;
495 #ifdef JEX_DEBUG
496 CONSOLE_DEBUG("Allocating 'variables' with length %d",i);
497 #endif
498 variables = ASC_NEW_ARRAY(struct var_variable*, i);
499 derivatives = ASC_NEW_ARRAY(double, i);
500
501 /* evaluate the derivatives... */
502 /* J = dG_dy = dF_dy + alpha * dF_dyp */
503
504 #ifdef ASC_SIGNAL_TRAPS
505 Asc_SignalHandlerPushDefault(SIGFPE);
506 if (SETJMP(g_fpe_env)==0) {
507 #endif
508 for(i=0, relptr = enginedata->rellist;
509 i< enginedata->nrels && relptr != NULL;
510 ++i, ++relptr
511 ){
512 /* get derivatives for this particular relation */
513 status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
514 #ifdef JEX_DEBUG
515 CONSOLE_DEBUG("Got derivatives against %d matching variables, status = %d", count,status);
516 #endif
517
518 if(status){
519 relname = rel_make_name(integ->system, *relptr);
520 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Calculation error in rel '%s'",relname);
521 ASC_FREE(relname);
522 is_error = 1;
523 break;
524 }
525
526 /*
527 Now we have the derivatives wrt each alg/diff variable in the
528 present equation. variables[] points into the varlist. need
529 a mapping from the varlist to the y and ydot lists.
530 */
531
532 Jv_i = 0;
533 for(j=0; j < count; ++j){
534 /* CONSOLE_DEBUG("j = %d, variables[j] = %d, n_y = %ld", j, variables[j], integ->n_y);
535 varname = var_make_name(integ->system, enginedata->varlist[variables[j]]);
536 if(varname){
537 CONSOLE_DEBUG("Variable %d '%s' derivative = %f", variables[j],varname,derivatives[j]);
538 ASC_FREE(varname);
539 }else{
540 CONSOLE_DEBUG("Variable %d (UNKNOWN!): derivative = %f",variables[j],derivatives[j]);
541 }
542 */
543
544 /* we don't calculate derivatives wrt indep var */
545 asc_assert(variables[j]>=0);
546 if(variables[j] == integ->x) continue;
547 #ifdef JEX_DEBUG
548 CONSOLE_DEBUG("j = %d: variables[j] = %d",j,var_sindex(variables[j]));
549 #endif
550 if(var_deriv(variables[j])){
551 #define DIFFINDEX integrator_ida_diffindex(integ,variables[j])
552 #ifdef JEX_DEBUG
553 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dydot[%d] = %f, v[%d] = %f)\n", i
554 , derivatives[j] * NV_Ith_S(v,DIFFINDEX)
555 , i, DIFFINDEX, derivatives[j]
556 , DIFFINDEX, NV_Ith_S(v,DIFFINDEX)
557 );
558 #endif
559 asc_assert(integ->ydot[DIFFINDEX]==variables[j]);
560 Jv_i += derivatives[j] * NV_Ith_S(v,DIFFINDEX) * c_j;
561 #undef DIFFINDEX
562 }else{
563 #define VARINDEX var_sindex(variables[j])
564 #ifdef JEX_DEBUG
565 asc_assert(integ->y[VARINDEX]==variables[j]);
566 fprintf(stderr,"Jv[%d] += %f (dF[%d]/dy[%d] = %f, v[%d] = %f)\n"
567 , i
568 , derivatives[j] * NV_Ith_S(v,VARINDEX)
569 , i, VARINDEX, derivatives[j]
570 , VARINDEX, NV_Ith_S(v,VARINDEX)
571 );
572 #endif
573 Jv_i += derivatives[j] * NV_Ith_S(v,VARINDEX);
574 #undef VARINDEX
575 }
576 }
577
578 NV_Ith_S(Jv,i) = Jv_i;
579 #ifdef JEX_DEBUG
580 CONSOLE_DEBUG("rel = %p",*relptr);
581 relname = rel_make_name(integ->system, *relptr);
582 CONSOLE_DEBUG("'%s': Jv[%d] = %f", relname, i, NV_Ith_S(Jv,i));
583 ASC_FREE(relname);
584 return 1;
585 #endif
586 }
587 #ifdef ASC_SIGNAL_TRAPS
588 }else{
589 relname = rel_make_name(integ->system, *relptr);
590 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error (SIGFPE) in rel '%s'",relname);
591 ASC_FREE(relname);
592 is_error = 1;
593 }
594 Asc_SignalHandlerPopDefault(SIGFPE);
595 #endif
596
597 if(is_error){
598 CONSOLE_DEBUG("SOME ERRORS FOUND IN EVALUATION");
599 return 1;
600 }
601 return 0;
602 }
603
604 /* sparse jacobian evaluation for IDAASCEND sparse direct linear solver */
605 int integrator_ida_sjex(long int Neq, realtype tt
606 , N_Vector yy, N_Vector yp, N_Vector rr
607 , realtype c_j, void *jac_data, mtx_matrix_t Jac
608 , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
609 ){
610 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
611 return -1;
612 }
613
614 /* root finding function */
615
616 int integrator_ida_rootfn(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *g_data){
617 IntegratorSystem *integ;
618 IntegratorIdaData *enginedata;
619 int i;
620 #ifdef ROOT_DEBUG
621 char *relname;
622 #endif
623
624 asc_assert(g_data!=NULL);
625 integ = (IntegratorSystem *)g_data;
626 enginedata = integrator_ida_enginedata(integ);
627
628 /* pass the values of everything back to the compiler */
629 integrator_set_t(integ, (double)tt);
630 integrator_set_y(integ, NV_DATA_S(yy));
631 integrator_set_ydot(integ, NV_DATA_S(yp));
632
633 asc_assert(gout!=NULL);
634
635 #ifdef ROOT_DEBUG
636 CONSOLE_DEBUG("t = %f",tt);
637 #endif
638
639 /* evaluate the residuals for each of the boundaries */
640 for(i=0; i < enginedata->nbnds; ++i){
641 switch(bnd_kind(enginedata->bndlist[i])){
642 case e_bnd_rel: /* real-valued boundary relation */
643 gout[i] = bndman_real_eval(enginedata->bndlist[i]);
644 #ifdef ROOT_DEBUG
645 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
646 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
647 ASC_FREE(relname);
648 #endif
649 break;
650 case e_bnd_logrel:
651 if(bndman_log_eval(enginedata->bndlist[i])){
652 CONSOLE_DEBUG("bnd[%d] = TRUE",i);
653 #ifdef ROOT_DEBUG
654 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
655 CONSOLE_DEBUG("gout[%d] = %f (boundary '%s')", i, gout[i], relname);
656 ASC_FREE(relname);
657 #endif
658 gout[i] = +1.0;
659 }else{
660 CONSOLE_DEBUG("bnd[%d] = FALSE",i);
661 gout[i] = -1.0;
662 }
663 break;
664 case e_bnd_undefined:
665 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid boundary type e_bnd_undefined");
666 return 1;
667 }
668 }
669
670 return 0; /* no way to detect errors in bndman_*_eval at this stage */
671 }
672

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