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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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