/[ascend]/branches/ksenija2/solvers/ida/ida.c
ViewVC logotype

Contents of /branches/ksenija2/solvers/ida/ida.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2882 - (show annotations) (download) (as text)
Sun Mar 29 12:20:00 2015 UTC (7 years, 8 months ago) by jpye
File MIME type: text/x-csrc
File size: 42806 byte(s)
more debugging info

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, see <http://www.gnu.org/licenses/>.
16 *//**
17 @file
18 Access to the IDA integrator for ASCEND. IDA is a DAE solver that comes
19 as part of the GPL-licensed SUNDIALS solver package from LLNL.
20
21 IDA provides the non-linear parts, as well as a number of pluggable linear
22 solvers: dense, banded and krylov types.
23
24 We also implement here an EXPERIMENTAL direct sparse linear solver for IDA
25 using the ASCEND linsolqr routines.
26
27 @see http://www.llnl.gov/casc/sundials/
28 *//*
29 by John Pye, May 2006
30 */
31
32 #define _GNU_SOURCE
33
34 #include "ida.h"
35 #include "idalinear.h"
36 #include "idaanalyse.h"
37 #include "idatypes.h"
38 #include "idaprec.h"
39 #include "idacalc.h"
40 #include "idaio.h"
41 #include "idaboundary.h"
42
43 #include <signal.h>
44 #include <setjmp.h>
45 #include <fenv.h>
46 #include <math.h>
47
48 #ifdef ASC_WITH_MMIO
49 # include <mmio.h>
50 #endif
51
52 #include <ascend/general/platform.h>
53 #include <ascend/utilities/error.h>
54 #include <ascend/utilities/ascSignal.h>
55 #include <ascend/general/panic.h>
56 #include <ascend/compiler/instance_enum.h>
57
58
59 #include <ascend/system/slv_client.h>
60 #include <ascend/system/relman.h>
61 #include <ascend/system/block.h>
62 #include <ascend/system/slv_stdcalls.h>
63 #include <ascend/system/jacobian.h>
64 #include <ascend/system/bndman.h>
65
66 #include <ascend/utilities/config.h>
67 #include <ascend/integrator/integrator.h>
68
69 /* #define FEX_DEBUG */
70 /* #define SOLVE_DEBUG */
71 #define IDA_BND_DEBUG
72 /* #define STATS_DEBUG */
73 /* #define DESTROY_DEBUG */
74
75 /*-------------------------------------------------------------
76 SOLVER REGISTRATION
77 */
78
79 /*
80 The following functions are declared static because we don't want them
81 ever to be called directly from outside code. Instead, they are call via
82 function pointers passed through the ida_register function.
83 */
84 static IntegratorCreateFn integrator_ida_create;
85 static IntegratorParamsDefaultFn integrator_ida_params_default;
86 static IntegratorSolveFn integrator_ida_solve;
87 static IntegratorFreeFn integrator_ida_free;
88
89 /**
90 This data structure contains pointers to the various functions that are
91 exposed to libascend in order for ASCEND to drive this solver.
92 */
93 static const IntegratorInternals integrator_ida_internals = {
94 integrator_ida_create, integrator_ida_params_default,
95 integrator_ida_analyse, integrator_ida_solve,
96 integrator_ida_write_matrix, integrator_ida_debug, integrator_ida_free,
97 INTEG_IDA, "IDA" };
98
99 /**
100 This function is accessed by libascend when loading this solver. The
101 function will register the integrator such that it can then be applied
102 to solving problems.
103 */
104 extern ASC_EXPORT int ida_register(void) {
105 //CONSOLE_DEBUG("Registering IDA...");
106 return integrator_register(&integrator_ida_internals);
107 }
108
109 /*-------------------------------------------------------------
110 FORWARD DECLS
111 */
112
113 typedef void ( IntegratorVarVisitorFn)(IntegratorSystem *integ,
114 struct var_variable *var, const int *varindx);
115
116 /*static IntegratorVarVisitorFn integrator_dae_classify_var;
117 static void integrator_visit_system_vars(IntegratorSystem *integ,IntegratorVarVisitorFn *visitor);
118 static void integrator_dae_show_var(IntegratorSystem *integ, struct var_variable *var, const int *varindx); */
119
120 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s);
121
122 /*-------------------------------------------------------------
123 SETUP/TEARDOWN ROUTINES
124 */
125
126 /**
127 Allocate memory inside the IntegratorIdaData structure at
128 the time when IDA is assigned to a particular system as its integrator.
129 */
130 static void integrator_ida_create(IntegratorSystem *integ) {
131 //CONSOLE_DEBUG("ALLOCATING IDA ENGINE DATA");
132 IntegratorIdaData *enginedata;
133 enginedata = ASC_NEW(IntegratorIdaData);
134 //CONSOLE_DEBUG("enginedata = %p",enginedata);
135 enginedata->rellist = NULL;
136 enginedata->safeeval = 0;
137 enginedata->vfilter.matchbits = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE
138 | VAR_FIXED;
139 enginedata->vfilter.matchvalue = VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0;
140 enginedata->pfree = NULL;
141
142 enginedata->rfilter.matchbits = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
143 enginedata->rfilter.matchvalue = REL_EQUALITY | REL_INCLUDED | REL_ACTIVE;
144
145 enginedata->flagfn = NULL;
146 enginedata->flagfntype = NULL;
147 enginedata->flagnamefn = NULL;
148
149 integ->enginedata = (void *) enginedata;
150
151 integrator_ida_params_default(integ);
152 }
153
154 /**
155 Free any memory specifically allocated for use by our ASCEND wrapper for
156 the IDA integrator, other than that which the IDA code itself may be
157 managing.
158 */
159 static void integrator_ida_free(void *enginedata) {
160 #ifdef DESTROY_DEBUG
161 CONSOLE_DEBUG("DESTROYING IDA engine data at %p",enginedata);
162 #endif
163 IntegratorIdaData *d = (IntegratorIdaData *) enginedata;
164 asc_assert(d);
165 if (d->pfree) {
166 CONSOLE_DEBUG("DESTROYING preconditioner data using fn at %p",d->pfree);
167 /* free the preconditioner data, whatever it happens to be */
168 (d->pfree)(enginedata);
169 }
170
171 ASC_FREE(d->rellist);
172
173 #ifdef DESTROY_DEBUG
174 CONSOLE_DEBUG("Now destroying the enginedata");
175 #endif
176 ASC_FREE(d);
177 #ifdef DESTROY_DEBUG
178 CONSOLE_DEBUG("enginedata freed");
179 #endif
180 }
181
182 IntegratorIdaData *integrator_ida_enginedata(IntegratorSystem *integ) {
183 IntegratorIdaData *d;
184 assert(integ!=NULL);
185 assert(integ->enginedata!=NULL);
186 assert(integ->engine==INTEG_IDA);
187 d = ((IntegratorIdaData *) (integ->enginedata));
188 return d;
189 }
190
191 /*-------------------------------------------------------------
192 PARAMETERS FOR IDA
193 */
194
195 enum ida_parameters {
196 IDA_PARAM_LINSOLVER,
197 IDA_PARAM_MAXL,
198 IDA_PARAM_MAXORD,
199 IDA_PARAM_AUTODIFF,
200 IDA_PARAM_CALCIC,
201 IDA_PARAM_SAFEEVAL,
202 IDA_PARAM_RTOL,
203 IDA_PARAM_ATOL,
204 IDA_PARAM_ATOLVECT,
205 IDA_PARAM_GSMODIFIED,
206 IDA_PARAM_MAXNCF,
207 IDA_PARAM_PREC,
208 IDA_PARAMS_SIZE
209 };
210
211 /**
212 Here the full set of parameters is defined, along with upper/lower bounds,
213 etc. The values are stuck into the integ->params structure.
214
215 To add a new parameter, first give it a name IDA_PARAM_* in thge above enum ida_parameters
216 list. Then add a slv_param_*(...) statement below to define the type, description and range
217 for the new parameter.
218
219 @return 0 on success
220 */
221 static int integrator_ida_params_default(IntegratorSystem *integ) {
222 asc_assert(integ!=NULL);
223 asc_assert(integ->engine==INTEG_IDA);
224 slv_parameters_t *p;
225 p = &(integ->params);
226
227 slv_destroy_parms(p);
228
229 if (p->parms == NULL) {
230 //CONSOLE_DEBUG("params NULL");
231 p->parms = ASC_NEW_ARRAY(struct slv_parameter, IDA_PARAMS_SIZE);
232 if (p->parms == NULL)
233 return -1;
234 p->dynamic_parms = 1;
235 } else {
236 //CONSOLE_DEBUG("params not NULL");
237 }
238
239 /* reset the number of parameters to zero so that we can check it at the end */
240 p->num_parms = 0;
241
242 slv_param_bool(
243 p,
244 IDA_PARAM_AUTODIFF,
245 (SlvParameterInitBool) { {"autodiff"
246 ,"Use auto-diff?",1
247 ,"Use automatic differentiation of expressions (1) or use numerical derivatives (0)"
248 }, TRUE}
249 );
250
251 slv_param_char(p,IDA_PARAM_CALCIC
252 ,(SlvParameterInitChar) { {"calcic"
253 ,"Initial conditions calcuation",1
254 ,"Use specified values of ydot to solve for inital y (Y),"
255 " or use the the values of the differential variables (yd) to solve"
256 " for the pure algebraic variables (ya) along with the derivatives"
257 " of the differential variables (yddot) (YA_YDP), or else don't solve"
258 " the intial conditions at all (NONE). See IDA manual p 41 (IDASetId)"
259 }, "YA_YDP"}, (char *[]) {"Y", "YA_YDP", "NONE",NULL}
260 );
261
262 slv_param_bool(p,IDA_PARAM_SAFEEVAL
263 ,(SlvParameterInitBool) { {"safeeval"
264 ,"Use safe evaluation?",1
265 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
266 "throw SIGFPE errors which will then halt integration."
267 }, FALSE}
268 );
269
270 slv_param_bool(p,IDA_PARAM_ATOLVECT
271 ,(SlvParameterInitBool) { {"atolvect"
272 ,"Use 'ode_atol' values as specified?",1
273 ,"If TRUE, values of 'ode_atol' are taken from your model and used "
274 " in the integration. If FALSE, a scalar absolute tolerance value"
275 " is shared by all variables. See IDA manual, section 5.5.1"
276 }, TRUE}
277 );
278
279 slv_param_real(p,IDA_PARAM_ATOL
280 ,(SlvParameterInitReal) { {"atol"
281 ,"Scalar absolute error tolerance",1
282 ,"Value of the scalar absolute error tolerance. See also 'atolvect'."
283 " See IDA manual, sections 5.5.1 and 5.5.2 'Advice on choice and use of tolerances'"
284 }, 1e-5, 0, 1e10}
285 );
286
287 slv_param_real(p,IDA_PARAM_RTOL
288 ,(SlvParameterInitReal) { {"rtol"
289 ,"Scalar relative error tolerance",1
290 ,"Value of the scalar relative error tolerance. (Note that for IDA,"
291 " it's not possible to set per-variable relative tolerances as it is"
292 " with LSODE)."
293 " See IDA manual, section 5.5.2 'Advice on choice and use of tolerances'"
294 }, 1e-4, 0, 1}
295 );
296
297 slv_param_char(p,IDA_PARAM_LINSOLVER
298 ,(SlvParameterInitChar) { {"linsolver"
299 ,"Linear solver",1
300 ,"See IDA manual, section 5.5.3. Choose 'ASCEND' to use the linsolqr"
301 " direct linear solver bundled with ASCEND, 'DENSE' to use the dense"
302 " solver bundled with IDA, or one of the Krylov solvers SPGMR, SPBCG"
303 " or SPTFQMR (which still need preconditioners to be implemented"
304 " before they can be very useful."
305 }, "DENSE"}, (char *[]) {"ASCEND","DENSE","BAND","SPGMR","SPBCG","SPTFQMR",NULL}
306 );
307
308 slv_param_int(p,IDA_PARAM_MAXL
309 ,(SlvParameterInitInt) { {"maxl"
310 ,"Maximum Krylov dimension",0
311 ,"The maximum dimension of Krylov space used by the linear solver"
312 " (for SPGMR, SPBCG, SPTFQMR) with IDA. See IDA manual section 5.5."
313 " The default of 0 results in IDA using its internal default, which"
314 " is currently a value of 5."
315 }, 0, 0, 20}
316 );
317
318 slv_param_int(p,IDA_PARAM_MAXORD
319 ,(SlvParameterInitInt) { {"maxord"
320 ,"Maximum order of linear multistep method",0
321 ,"The maximum order of the linear multistep method with IDA. See"
322 " IDA manual p 38."
323 }, 5, 1, 5}
324 );
325
326 slv_param_bool(p,IDA_PARAM_GSMODIFIED
327 ,(SlvParameterInitBool) { {"gsmodified"
328 ,"Gram-Schmidt Orthogonalisation Scheme", 2
329 ,"TRUE = GS_MODIFIED, FALSE = GS_CLASSICAL. See IDA manual section"
330 " 5.5.6.6. Only applies when linsolve=SPGMR is selected."
331 }, TRUE}
332 );
333
334 slv_param_int(p,IDA_PARAM_MAXNCF
335 ,(SlvParameterInitInt) { {"maxncf"
336 ,"Max nonlinear solver convergence failures per step", 2
337 ,"Maximum number of allowable nonlinear solver convergence failures"
338 " on one step. See IDA manual section 5.5.6.1."
339 }, 10,0,1000}
340 );
341
342 slv_param_char(p,IDA_PARAM_PREC
343 ,(SlvParameterInitChar) { {"prec"
344 ,"Preconditioner",1
345 ,"See IDA manual, section section 5.6.8."
346 },"NONE"}, (char *[]) {"NONE","DIAG",NULL}
347 );
348
349 asc_assert(p->num_parms == IDA_PARAMS_SIZE);
350
351 //CONSOLE_DEBUG("Created %d params", p->num_parms);
352
353 return 0;
354 }
355
356
357
358 /*******************************************
359 * SOLVE SETUP FUNCTIONS
360 *******************************************/
361
362 int ida_load_rellist(IntegratorSystem *integ) {
363 IntegratorIdaData *enginedata;
364 struct rel_relation **rels;
365 int i, j, n_solverrels, n_active_rels;
366
367 enginedata = integrator_ida_enginedata(integ);
368
369 n_solverrels = slv_get_num_solvers_rels(integ->system);
370 n_active_rels = slv_count_solvers_rels(integ->system, &integrator_ida_rel);
371 rels = slv_get_solvers_rel_list(integ->system);
372
373 if (enginedata->rellist != NULL) {
374 ASC_FREE(enginedata->rellist);
375 enginedata->rellist = NULL;
376 }
377
378 enginedata->rellist
379 = ASC_NEW_ARRAY(struct rel_relation *, n_active_rels);
380
381 #ifdef SOLVE_DEBUG
382 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
383 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
384
385 CONSOLE_DEBUG("Number of relations: %d",n_solverrels);
386 CONSOLE_DEBUG("Number of active relations: %d",n_active_rels);
387 CONSOLE_DEBUG("Number of dependent vars: %d",integ->n_y);
388 CONSOLE_DEBUG("Number of boundaries: %d",enginedata->nbnds);
389 #endif
390
391
392 j = 0;
393 for (i = 0; i < n_solverrels; ++i) {
394 if (rel_apply_filter(rels[i], &integrator_ida_rel)) {
395 #ifdef SOLVE_DEBUG
396 {
397 char *relname = rel_make_name(integ->system, rels[i]);
398 CONSOLE_DEBUG("rel '%s': 0x%x", relname, rel_flags(rels[i]));
399 ASC_FREE(relname);
400 }
401 #endif
402 enginedata->rellist[j++] = rels[i];
403 }
404 }
405
406 asc_assert(j == n_active_rels);
407 enginedata->nrels = n_active_rels;
408
409 if (enginedata->nrels != integ->n_y) {
410 ERROR_REPORTER_HERE(ASC_USER_ERROR
411 ,"Integration problem is not square (%d active rels, %d vars)"
412 ,n_active_rels, integ->n_y
413 );
414 return 1; /* failure */
415 }
416
417 return 0;
418 }
419
420 /*
421 * Retrieve initial values from the system
422 */
423 int ida_retrieve_IVs(IntegratorSystem *integ, realtype t0, N_Vector y0,
424 N_Vector yp0) {
425
426 #ifdef SOLVE_DEBUG
427 char *varname;
428 char diffname[30];
429 int i;
430 CONSOLE_DEBUG("RETRIEVING INITIAL VALUES:");
431 CONSOLE_DEBUG("t0 = %f",t0);
432 #endif
433
434 integrator_get_y(integ, NV_DATA_S(y0));
435 integrator_get_ydot(integ, NV_DATA_S(yp0));
436
437 #ifdef SOLVE_DEBUG
438 fprintf(stderr, "index\t%25s\t%25s\n", "y", "ydot");
439 for (i = 0; i < integ->n_y; ++i) {
440 varname = var_make_name(integ->system, integ->y[i]);
441 fprintf(stderr, "%d\t%15s=%10f\t", i, varname, NV_Ith_S(y0,i));
442
443 if (integ->ydot[i]) {
444 ASC_FREE(varname);
445 varname = var_make_name(integ->system, integ->ydot[i]);
446 fprintf(stderr, "%15s=%10f\t\n", varname, NV_Ith_S(yp0,i));
447 } else {
448 snprintf(diffname,99,"diff(%s)",varname);
449 fprintf(stderr,"%15s=%10f\t\n",diffname,NV_Ith_S(yp0,i));
450 }
451 ASC_FREE(varname);
452 }
453 #endif
454
455 return 0;
456 }
457
458 /*
459 * Assign internal memory for the appropriate sundials version and error
460 * tolerance parameter settings.
461 */
462 int ida_malloc(IntegratorSystem *integ, void *ida_mem, realtype t0,
463 N_Vector y0, N_Vector yp0) {
464 int flag;
465 N_Vector abstolvect;
466 realtype reltol, abstol;
467
468 /* relative error tolerance */
469 reltol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_RTOL);
470 //CONSOLE_DEBUG("rtol = %8.2e",reltol);
471
472
473 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
474 flag = IDAInit(ida_mem, &integrator_ida_fex, t0, y0 ,yp0);
475 #else
476 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)) {
477 /* vector of absolute tolerances */
478 CONSOLE_DEBUG("USING VECTOR OF ATOL VALUES");
479 abstolvect = N_VNew_Serial(integ->n_y);
480 integrator_get_atol(integ, NV_DATA_S(abstolvect));
481
482 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SV,
483 reltol, abstolvect);
484
485 N_VDestroy_Serial(abstolvect);
486 }else{
487 /* scalar absolute tolerance (one value for all) */
488 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
489 CONSOLE_DEBUG("USING SCALAR ATOL VALUE = %8.2e",abstol);
490 flag = IDAMalloc(ida_mem, &integrator_ida_fex, t0, y0, yp0, IDA_SS,
491 reltol, &abstol);
492 }
493 #endif
494
495 if(flag == IDA_MEM_NULL) {
496 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
497 return 2;
498 }else if(flag == IDA_MEM_FAIL) {
499 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDAMalloc)");
500 return 3;
501 }else if(flag == IDA_ILL_INPUT) {
502 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid input to IDAMalloc");
503 return 4;
504 }
505
506 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
507 //CONSOLE_DEBUG("Assigning tolerances...");
508 /* assign tolerances */
509 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_ATOLVECT)) {
510 //CONSOLE_DEBUG("using vector of atol values");
511 abstolvect = N_VNew_Serial(integ->n_y);
512 integrator_get_atol(integ,NV_DATA_S(abstolvect));
513 IDASVtolerances(ida_mem, reltol, abstolvect);
514 N_VDestroy_Serial(abstolvect);
515 }else{
516 /* scalar tolerances */
517 abstol = SLV_PARAM_REAL(&(integ->params),IDA_PARAM_ATOL);
518 CONSOLE_DEBUG("using scalar atol value = %8.2e",abstol);
519 IDASStolerances(ida_mem, reltol, abstol);
520 }
521 #endif
522
523 /* success */
524 return 0;
525 }
526
527 /*
528 * Set parameter inputs for step size, the linear solver module and preconditioner
529 */
530 int ida_set_optional_inputs(IntegratorSystem *integ, void *ida_mem) {
531 int flag;
532 char *linsolver;
533 char *pname = NULL;
534 int maxl;
535 const IntegratorIdaPrec *prec = NULL;
536
537 IntegratorIdaData *enginedata = integ->enginedata;
538
539 IDASetErrHandlerFn(ida_mem, &integrator_ida_error, (void *) integ);
540 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
541 IDASetUserData(ida_mem, (void *)integ);
542 #else
543 IDASetRdata(ida_mem, (void *) integ);
544 #endif
545 IDASetMaxStep(ida_mem, integrator_get_maxstep(integ));
546 IDASetInitStep(ida_mem, integrator_get_stepzero(integ));
547 IDASetMaxNumSteps(ida_mem, integrator_get_maxsubsteps(integ));
548 if (integrator_get_minstep(integ) > 0) {
549 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"IDA does not support minstep (ignored)\n");
550 }
551
552 //CONSOLE_DEBUG("MAXNCF = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
553 IDASetMaxConvFails(ida_mem, SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXNCF));
554
555 //CONSOLE_DEBUG("MAXORD = %d",SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
556 IDASetMaxOrd(ida_mem, SLV_PARAM_INT(&integ->params,IDA_PARAM_MAXORD));
557
558 /* there's no capability for setting *minimum* step size in IDA */
559
560 /* attach linear solver module, using the default value of maxl */
561 linsolver = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_LINSOLVER);
562 //CONSOLE_DEBUG("ASSIGNING LINEAR SOLVER '%s'",linsolver);
563 if(strcmp(linsolver, "ASCEND") == 0) {
564 CONSOLE_DEBUG("ASCEND DIRECT SOLVER, size = %d",integ->n_y);
565 IDAASCEND(ida_mem, integ->n_y);
566 IDAASCENDSetJacFn(ida_mem, &integrator_ida_sjex, (void *) integ);
567
568 enginedata->flagfntype = "IDAASCEND";
569 enginedata->flagfn = &IDAASCENDGetLastFlag;
570 enginedata->flagnamefn = &IDAASCENDGetReturnFlagName;
571
572 }else if (strcmp(linsolver, "DENSE") == 0) {
573 //CONSOLE_DEBUG("DENSE DIRECT SOLVER, size = %d",integ->n_y);
574 flag = IDADense(ida_mem, integ->n_y);
575 switch(flag){
576 case IDADENSE_SUCCESS:
577 break;
578 case IDADENSE_MEM_NULL:
579 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
580 return 5;
581 case IDADENSE_ILL_INPUT:
582 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDADENSE is not compatible with current nvector module");
583 return 5;
584 case IDADENSE_MEM_FAIL:
585 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Memory allocation failed for IDADENSE");
586 return 5;
587 default:
588 ERROR_REPORTER_HERE(ASC_PROG_ERR,"bad return");
589 return 5;
590 }
591
592 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)) {
593 //CONSOLE_DEBUG("USING AUTODIFF");
594 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
595 flag = IDADlsSetDenseJacFn(ida_mem, &integrator_ida_djex);
596 #else
597 flag = IDADenseSetJacFn(ida_mem, &integrator_ida_djex,
598 (void *) integ);
599 #endif
600 switch (flag) {
601 case IDADENSE_SUCCESS:
602 break;
603 default:
604 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed IDADenseSetJacFn");
605 return 6;
606 }
607 }else{
608 CONSOLE_DEBUG("USING NUMERICAL DIFF");
609 }
610
611 enginedata->flagfntype = "IDADENSE";
612 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
613 enginedata->flagfn = &IDADlsGetLastFlag;
614 enginedata->flagnamefn = &IDADlsGetReturnFlagName;
615 #else
616 enginedata->flagfn = &IDADenseGetLastFlag;
617 enginedata->flagnamefn = &IDADenseGetReturnFlagName;
618 #endif
619 }else{
620 /* remaining methods are all SPILS */
621 CONSOLE_DEBUG("IDA SPILS");
622
623 maxl = SLV_PARAM_INT(&(integ->params),IDA_PARAM_MAXL);
624 CONSOLE_DEBUG("maxl = %d",maxl);
625
626 /* what preconditioner? */
627 pname = SLV_PARAM_CHAR(&(integ->params),IDA_PARAM_PREC);
628 if(strcmp(pname, "NONE") == 0){
629 prec = NULL;
630 }else if(strcmp(pname, "JACOBI") == 0){
631 prec = &prec_jacobi;
632 }else{
633 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid preconditioner choice '%s'",pname);
634 return 7;
635 }
636
637 /* which SPILS linear solver? */
638 if(strcmp(linsolver, "SPGMR") == 0){
639 CONSOLE_DEBUG("IDA SPGMR");
640 flag = IDASpgmr(ida_mem, maxl); /* 0 means use the default max Krylov dimension of 5 */
641 }else if(strcmp(linsolver, "SPBCG") == 0){
642 CONSOLE_DEBUG("IDA SPBCG");
643 flag = IDASpbcg(ida_mem, maxl);
644 }else if(strcmp(linsolver, "SPTFQMR") == 0){
645 CONSOLE_DEBUG("IDA SPTFQMR");
646 flag = IDASptfqmr(ida_mem, maxl);
647 }else{
648 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown IDA linear solver choice '%s'",linsolver);
649 return 8;
650 }
651
652 if(prec){
653 /* assign the preconditioner to the linear solver */
654 (prec->pcreate)(integ);
655 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
656 IDASpilsSetPreconditioner(ida_mem,prec->psetup,prec->psolve);
657 #else
658 IDASpilsSetPreconditioner(ida_mem, prec->psetup, prec->psolve,
659 (void *) integ);
660 #endif
661 CONSOLE_DEBUG("PRECONDITIONER = %s",pname);
662 }else{
663 CONSOLE_DEBUG("No preconditioner");
664 }
665
666 enginedata->flagfntype = "IDASPILS";
667 enginedata->flagfn = &IDASpilsGetLastFlag;
668 enginedata->flagnamefn = &IDASpilsGetReturnFlagName;
669
670 if(flag == IDASPILS_MEM_NULL){
671 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
672 return 9;
673 }else if (flag == IDASPILS_MEM_FAIL) {
674 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to allocate memory (IDASpgmr)");
675 return 9;
676 }/* else success */
677
678
679 /* assign the J*v function */
680 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_AUTODIFF)) {
681 CONSOLE_DEBUG("USING AUTODIFF");
682 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
683 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex);
684 #else
685 flag = IDASpilsSetJacTimesVecFn(ida_mem, &integrator_ida_jvex,
686 (void *) integ);
687 #endif
688 if(flag == IDASPILS_MEM_NULL) {
689 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ida_mem is NULL");
690 return 10;
691 }else if (flag == IDASPILS_LMEM_NULL) {
692 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IDASPILS linear solver has not been initialized");
693 return 10;
694 }/* else success */
695 }else{
696 CONSOLE_DEBUG("USING NUMERICAL DIFF");
697 }
698
699 if(strcmp(linsolver, "SPGMR") == 0) {
700 /* select Gram-Schmidt orthogonalisation */
701 if(SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_GSMODIFIED)) {
702 CONSOLE_DEBUG("USING MODIFIED GS");
703 flag = IDASpilsSetGSType(ida_mem, MODIFIED_GS);
704 if(flag != IDASPILS_SUCCESS) {
705 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
706 return 11;
707 }
708 }else{
709 CONSOLE_DEBUG("USING CLASSICAL GS");
710 flag = IDASpilsSetGSType(ida_mem, CLASSICAL_GS);
711 if(flag != IDASPILS_SUCCESS) {
712 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to set GS_MODIFIED");
713 return 11;
714 }
715 }
716 }
717 }
718
719 /* set linear solver optional inputs...
720 ...nothing here at the moment...
721 */
722
723 return 0;
724 } /* ida_set_optional_inputs */
725
726 /**
727 * Calculate initial conditions using IDACalcIC
728 */
729
730 int ida_setup_IC(IntegratorSystem *integ, void *ida_mem,
731 realtype tout1, realtype t0, N_Vector y0, N_Vector yp0) {
732 int i, flag, flag1;
733 int icopt; /* initial conditions strategy */
734 N_Vector id;
735
736 IntegratorIdaData *enginedata = integ->enginedata;
737
738 #ifdef SOLVE_DEBUG
739 char *varname;
740 #endif
741
742 #ifdef IDA_BND_DEBUG
743 CONSOLE_DEBUG("Solving initial conditions...");
744 #endif
745
746 icopt = 0;
747 if(strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "Y") == 0) {
748 #ifdef SOLVE_DEBUG
749 CONSOLE_DEBUG("Solving initial conditions using values of yddot");
750 #endif
751 icopt = IDA_Y_INIT;
752 asc_assert(icopt!=0);
753 }else if(0==strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "YA_YDP")){
754 #ifdef SOLVE_DEBUG
755 CONSOLE_DEBUG("Solving initial conditions using values of yd");
756 #endif
757 icopt = IDA_YA_YDP_INIT;
758 asc_assert(icopt!=0);
759 id = N_VNew_Serial(integ->n_y);
760 for(i = 0; i < integ->n_y; ++i) {
761 if(integ->ydot[i] == NULL) {
762 NV_Ith_S(id, i) = 0.0;
763 #ifdef SOLVE_DEBUG
764 varname = var_make_name(integ->system, integ->y[i]);
765 CONSOLE_DEBUG("y[%d] = '%s' is pure algebraic",i,varname);
766 ASC_FREE(varname);
767 #endif
768 }else{
769 #ifdef SOLVE_DEBUG
770 CONSOLE_DEBUG("y[%d] is differential",i);
771 #endif
772 NV_Ith_S(id, i) = 1.0;
773 }
774 }
775 IDASetId(ida_mem, id);
776 N_VDestroy_Serial(id);
777 }else if (strcmp(SLV_PARAM_CHAR(&integ->params,IDA_PARAM_CALCIC), "NONE")
778 == 0) {
779 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Not solving initial conditions: check current residuals");
780 } else {
781 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid 'iccalc' value: check solver parameters.");
782 }
783
784 if (icopt) {
785 #ifdef SOLVE_DEBUG
786 CONSOLE_DEBUG("SOLVING INITIAL CONDITIONS IDACalcIC (tout1 = %f)", tout1);
787 #endif
788
789 #ifdef ASC_SIGNAL_TRAPS
790 /* catch SIGFPE if desired to */
791 if(enginedata->safeeval) {
792 CONSOLE_DEBUG("SETTING TO IGNORE SIGFPE...");
793 Asc_SignalHandlerPush(SIGFPE, SIG_DFL);
794 }else {
795 # ifdef FEX_DEBUG
796 CONSOLE_DEBUG("SETTING TO CATCH SIGFPE...");
797 # endif
798 Asc_SignalHandlerPushDefault(SIGFPE);
799 }
800 if(setjmp(g_fpe_env) == 0) {
801 #endif
802
803 # if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=3
804 flag = IDACalcIC(ida_mem, icopt, tout1);/* new API from v2.3 */
805 # else
806 flag = IDACalcIC(ida_mem, t0, y0, yp0, icopt, tout1);
807 # endif
808 /* check flags and output status */
809 switch(flag){
810 case IDA_SUCCESS:
811 #ifdef SOLVE_DEBUG
812 CONSOLE_DEBUG("Initial conditions solved OK");
813 #endif
814 break;
815
816 case IDA_LSETUP_FAIL:
817 case IDA_LINIT_FAIL:
818 case IDA_LSOLVE_FAIL:
819 case IDA_NO_RECOVERY:
820 flag1 = -999;
821 flag = (enginedata->flagfn)(ida_mem, &flag1);
822 if (flag) {
823 ERROR_REPORTER_HERE(ASC_PROG_ERR
824 ,"Unable to retrieve error code from %s (err %d)"
825 ,enginedata->flagfntype,flag
826 );
827 return 12;
828 }
829 ERROR_REPORTER_HERE(ASC_PROG_ERR
830 ,"%s returned flag '%s' (value = %d)"
831 ,enginedata->flagfntype,(enginedata->flagnamefn)(flag1),flag1
832 );
833 return 12;
834
835 default:
836 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve initial condition (IDACalcIC)");
837 return 12;
838 }
839 #ifdef ASC_SIGNAL_TRAPS
840 }else{
841 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Floating point error while solving initial conditions");
842 return 13;
843 }
844
845 if(enginedata->safeeval){
846 Asc_SignalHandlerPop(SIGFPE, SIG_DFL);
847 }else{
848 /* CONSOLE_DEBUG("pop..."); */
849 Asc_SignalHandlerPopDefault(SIGFPE);
850 /* CONSOLE_DEBUG("...pop"); */
851 }
852 #endif
853 }/* icopt */
854
855 return 0;
856 } /* ida_setup_IC */
857
858
859 int ida_root_init(IntegratorSystem *integ, void *ida_mem) {
860 IntegratorIdaData *enginedata = integ->enginedata;
861
862 if(enginedata->nbnds) {
863 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
864 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn);
865 #else
866 IDARootInit(ida_mem, enginedata->nbnds, &integrator_ida_rootfn,
867 (void *) integ);
868 #endif
869 }
870
871 return 0;
872
873 }
874
875 /**
876 * Allocate memory and prepare for integration.
877 *
878 * @param t the first point at which a solution is desired. Required by IdaCalcIC
879 *
880 * @return 0 on success
881 */
882 int ida_prepare_integrator(IntegratorSystem *integ, void *ida_mem,
883 realtype tout1) {
884 realtype t0;
885 N_Vector y0, yp0;
886
887 y0 = ida_bnd_new_zero_NV(integ->n_y);
888 yp0 = ida_bnd_new_zero_NV(integ->n_y);
889
890 #if 0
891 int i;
892 double val;
893 CONSOLE_DEBUG("Values of the derivatives present in the model");
894 for(i=0; i < integ->n_y; i++) {
895 if(integ->ydot[i]){
896 val = var_value(integ->ydot[i]);
897 CONSOLE_DEBUG("ydot[%d]= %g", i, val);
898 }
899 }
900 #endif
901
902 t0 = integrator_get_t(integ);
903 ida_retrieve_IVs(integ, t0, y0, yp0);
904
905 #ifdef IDA_BND_DEBUG
906 CONSOLE_DEBUG("Initial values BEFORE IDACalcIC: y0 =");
907 N_VPrint_Serial(y0);
908 CONSOLE_DEBUG("yp0 = ");
909 N_VPrint_Serial(yp0);
910 #endif
911
912 /* allocate internal memory */
913 ida_malloc(integ, ida_mem, t0, y0, yp0);
914
915 /* set optional inputs... */
916 ida_set_optional_inputs(integ, ida_mem);
917
918 /* calculate initial conditions */
919 ida_setup_IC(integ, ida_mem, tout1, t0, y0, yp0);
920
921 /* specify ROOT-FINDING problem (if necessary) */
922 ida_root_init(integ, ida_mem);
923
924 /* Clean up */
925 N_VDestroy_Serial(y0);
926 N_VDestroy_Serial(yp0);
927 return 0;
928 }
929
930 /**
931 * Reinitialise IDA after boundary-crossing. We assume that IDAInit has been called.
932 * Works with Sundials 2.4.0 and later versions.
933 *
934 * @param t the first point at which a solution is desired. Required by IdaCalcIC
935 *
936 * @return 0 on success
937 */
938 int ida_reinit_integrator(IntegratorSystem *integ, void *ida_mem,
939 realtype tout1) {
940 realtype t0;
941 N_Vector y0, yp0;
942
943 y0 = ida_bnd_new_zero_NV(integ->n_y);
944 yp0 = ida_bnd_new_zero_NV(integ->n_y);
945
946 t0 = integrator_get_t(integ);
947 ida_retrieve_IVs(integ, t0, y0, yp0);
948
949 int flag;
950
951 flag = IDAReInit(ida_mem, t0, y0, yp0);
952 if(flag!=IDA_SUCCESS){
953 ERROR_REPORTER_HERE(ASC_PROG_ERR, "Reinitialisation failed.");
954 }
955
956 /* calculate initial conditions */
957 ida_setup_IC(integ, ida_mem, tout1, t0, y0, yp0);
958
959 /* Clean up */
960 N_VDestroy_Serial(y0);
961 N_VDestroy_Serial(yp0);
962 return 0;
963 }
964
965 /*-------------------------------------------------------------
966 MAIN IDA SOLVER ROUTINE, see IDA manual, sec 5.4, p. 27 ff.
967 */
968
969 /**
970 Main IDA solver routine. We can assume that the integrator_ida_analyse
971 function will already have been run.
972
973 The presence of 'start_index' and 'finish_index' is a hangover from the
974 Tcl/Tk GUI. We would like to get rid of them eventually.
975
976 @return 0 on success */
977 static int integrator_ida_solve(IntegratorSystem *integ,
978 unsigned long start_index, unsigned long finish_index) {
979 void *ida_mem;
980 int t_index;
981 realtype t0, tout, tret, tol;
982 realtype tmin; /** < The length of a small additional step made after an event is triggered */
983 N_Vector ypret, yret;
984 IntegratorIdaData *enginedata;
985 int i, flag;
986
987 int *rootsfound; /** < IDA rootfinder reports root index in here */
988 int *rootdir = NULL; /** < Used to tell IDA to ignore doulve crossings */
989 int *bnd_cond_states; /** < Record of boundary states so that IDA can tell LRSlv
990 how to evaluate a boundary crossing */
991 int *bnd_not_set;
992 int all_bnds_set = 1;
993
994 int need_to_reconfigure; /** < Flag to indicate system rebuild after crossing */
995 int need_to_reinteg = 0; /** < Flag for when crossings happen on or very close to timesteps */
996 int skipping_output; /** < Flag to skip output to reporter */
997 int qrslv_ind, lrslv_ind;
998
999 int after_root = 0;
1000
1001 #if defined(SOLVE_DEBUG) || defined(IDA_BND_DEBUG)
1002 char *relname;
1003 #endif
1004
1005 //CONSOLE_DEBUG("STARTING IDA...");
1006 /* Setup boundary list */
1007 enginedata = integrator_ida_enginedata(integ);
1008 enginedata->bndlist = slv_get_solvers_bnd_list(integ->system);
1009 enginedata->nbnds = slv_get_num_solvers_bnds(integ->system);
1010 enginedata->safeeval = SLV_PARAM_BOOL(&(integ->params),IDA_PARAM_SAFEEVAL);
1011 //CONSOLE_DEBUG("safeeval = %d",enginedata->safeeval);
1012
1013 qrslv_ind = slv_lookup_client("QRSlv");
1014 lrslv_ind = slv_lookup_client("LRSlv");
1015
1016 bnd_not_set = ASC_NEW_ARRAY(int,enginedata->nbnds);
1017
1018 #ifdef SOLVE_DEBUG
1019 integrator_ida_debug(integ, stderr);
1020 #endif
1021
1022 /* Initialise boundary condition states if appropriate. Reconfigure if necessary */
1023 if(enginedata->nbnds){
1024 CONSOLE_DEBUG("Initialising boundary states");
1025 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR<4
1026 ERROR_REPORTER_HERE(ASC_PROG_WARNING, "Warning: boundary detection is"
1027 "unreliable with SUNDIALS pre version 2.4.0. Please update if you"
1028 "wish to use IDA for conditional integration");
1029 #endif
1030 bnd_cond_states = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
1031
1032 /* identify if we're exactly *on* any boundaries currently */
1033 for(i = 0; i < enginedata->nbnds; i++) {
1034 #ifdef IDA_BND_DEBUG
1035 relname = bnd_make_name(integ->system,enginedata->bndlist[i]);
1036 #endif
1037 bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);
1038 bnd_set_ida_first_cross(enginedata->bndlist[i],1);
1039 if(bndman_real_eval(enginedata->bndlist[i]) == 0) {
1040 /* if the residual for the boundary is zero (ie looks like we are *on* the boundary?) JP */
1041 #ifdef IDA_BND_DEBUG
1042 CONSOLE_DEBUG("Boundary '%s': not set",relname);
1043 #endif
1044 bnd_not_set[i] = 1;
1045 all_bnds_set = 0;
1046 }else{
1047 bnd_not_set[i] = 0;
1048 }
1049 #ifdef IDA_BND_DEBUG
1050 CONSOLE_DEBUG("Boundary '%s' is %d",relname,bnd_cond_states[i]);
1051 ASC_FREE(relname);
1052 #endif
1053 }
1054 CONSOLE_DEBUG("Setting up LRSlv...");
1055 if(ida_setup_lrslv(integ,qrslv_ind,lrslv_ind)){
1056 ERROR_REPORTER_HERE(ASC_USER_ERROR, "Idaanalyse failed.");
1057 return 1;
1058 }
1059
1060 }
1061
1062 /* store reference to list of relations (in enginedata) */
1063 ida_load_rellist(integ);
1064
1065 /* create IDA object */
1066 ida_mem = IDACreate();
1067
1068 /* Setup parameter inputs and initial conditions for IDA. */
1069 tout = samplelist_get(integ->samples, start_index + 1);
1070 /* solve the initial conditions, allocate memory, other stuff... */
1071 ida_prepare_integrator(integ, ida_mem, tout);
1072
1073 tol = 0.0001*(samplelist_get(integ->samples, finish_index)
1074 - samplelist_get(integ->samples, start_index))
1075 /samplelist_length(integ->samples);
1076 tmin = tol;
1077
1078 /* -- set up the IntegratorReporter */
1079 integrator_output_init(integ);
1080
1081 /* -- store the initial values of all the stuff */
1082 integrator_output_write(integ);
1083 integrator_output_write_obs(integ);
1084
1085 /* specify where the returned values should be stored */
1086 yret = ida_bnd_new_zero_NV(integ->n_y);
1087 ypret = ida_bnd_new_zero_NV(integ->n_y);
1088
1089 /* advance solution in time, return values as yret and derivatives as ypret */
1090 integ->currentstep = 1;
1091
1092 rootdir = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
1093
1094 for(t_index = start_index + 1; t_index <= finish_index; ++t_index, ++integ->currentstep) {
1095 tout = samplelist_get(integ->samples, t_index);
1096 t0 = integrator_get_t(integ);
1097
1098 asc_assert(tout > t0);
1099
1100 #ifdef SOLVE_DEBUG
1101 CONSOLE_DEBUG("Integrating from t0 = %f to t = %f", t0, tout);
1102 #endif
1103
1104 if(!after_root){
1105 /* reset the root-crossing-direction flags */
1106 for(i = 0; i < enginedata->nbnds; i++){
1107 rootdir[i] = 0;
1108 #ifdef IDA_BND_DEBUG
1109 char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);
1110 CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d]=%d, bndman_calc_satisfied=%d; (trigger dirn=%s)"
1111 ,n, i, bnd_cond_states[i], bndman_calc_satisfied(enginedata->bndlist[i])
1112 ,rootdir[i]==1?"UP":(rootdir[i]==-1?"DOWN":"both")
1113 );
1114 ASC_FREE(n);
1115 #endif
1116 }
1117
1118 if(enginedata->nbnds) IDASetRootDirection(ida_mem, rootdir);
1119 }
1120
1121 /**
1122 * do { solve routine } while(need_to_reintegrate)
1123 * If IDA detects a boundary sufficiently far from the desired output
1124 * time tout, then the system is reconfigured and IDASolve is recalled
1125 * up to the same tindex.
1126 *
1127 * @TODO "Sufficiently far" is currently tol = 0.0001, i.e complete arbitrary
1128 * Should reflect IDACalcIC's requirements?
1129 */
1130 do{
1131 if(need_to_reinteg && !after_root){
1132 #ifdef IDA_BND_DEBUG
1133 CONSOLE_DEBUG("Resuming integration from %f to %f", integrator_get_t(integ), tout);
1134 #endif
1135 integrator_output_write(integ);
1136 }
1137
1138 #ifdef IDA_BND_DEBUG
1139 if(after_root && (tout - tret) > tol){
1140 CONSOLE_DEBUG("Resuming integration from %f to %f", integrator_get_t(integ), tret + tmin);
1141 }
1142 #endif
1143
1144 /* Control flags for boundary crossings */
1145 need_to_reinteg = 0;
1146 skipping_output = 0;
1147
1148 #ifdef ASC_SIGNAL_TRAPS
1149 Asc_SignalHandlerPushDefault(SIGINT);
1150 if(setjmp(g_int_env) == 0) {
1151 #endif
1152 if(after_root && (tout - tret) > tmin){
1153 CONSOLE_DEBUG("Continuing (slightly after root)...");
1154 /* after a root, we advance time by a tiny amount to ensure
1155 we traverse past the switching poing -- FIXME a bit hacky? */
1156 flag = IDASolve(ida_mem, tret + tmin, &tret, yret, ypret, IDA_NORMAL);
1157 }else{
1158 CONSOLE_DEBUG("Continuing immediately...");
1159 /* if we didn't pass a root, continue immediately from prev end time */
1160 flag = IDASolve(ida_mem, tout, &tret, yret, ypret, IDA_NORMAL);
1161 }
1162 if(after_root){
1163 CONSOLE_DEBUG("Running logic solver after root...");
1164 ida_log_solve(integ, lrslv_ind);
1165 after_root = 0;
1166 }
1167 #ifdef ASC_SIGNAL_TRAPS
1168 }else{
1169 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Caught interrupt");
1170 flag = -555;
1171 }
1172 Asc_SignalHandlerPopDefault(SIGINT);
1173 #endif
1174
1175 if(enginedata->nbnds) {
1176 if(flag == IDA_ROOT_RETURN){
1177 #ifdef IDA_BND_DEBUG
1178 CONSOLE_DEBUG("IDA reports root found!");
1179 #endif
1180 /* Store the root index */
1181 rootsfound = ASC_NEW_ARRAY_CLEAR(int,enginedata->nbnds);
1182
1183 if(IDA_SUCCESS != IDAGetRootInfo(ida_mem, rootsfound)) {
1184 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch boundary-crossing info");
1185 return 14;
1186 }
1187
1188 #ifdef IDA_BND_DEBUG
1189 /* write out the boundaries that were crossed */
1190 for(i = 0; i < enginedata->nbnds; i++) {
1191 if(rootsfound[i]) {
1192 relname = bnd_make_name(integ->system, enginedata->bndlist[i]);
1193 CONSOLE_DEBUG("Boundary '%s' crossed at time x = %f, direction %s"
1194 ,relname,tret,(rootsfound>0?"UP":"DOWN")
1195 );
1196 ASC_FREE(relname);
1197 }
1198 }
1199 #endif
1200
1201
1202 if(all_bnds_set == 0){
1203 #ifdef IDA_BND_DEBUG
1204 CONSOLE_DEBUG("Unset bounds exist; evaluate them explicitly...");
1205 #endif
1206 all_bnds_set = 1;
1207 for(i = 0; i < enginedata->nbnds; i++){
1208 if(bnd_not_set[i]){
1209 if(!rootsfound[i]){
1210 bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);
1211 #ifdef IDA_BND_DEBUG
1212 relname = bnd_make_name(integ->system, enginedata->bndlist[i]);
1213 CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d] = %d"
1214 ,relname,i,bnd_cond_states[i]
1215 );
1216 ASC_FREE(relname);
1217 #endif
1218 }else all_bnds_set = 0;
1219 }
1220 }
1221 }
1222
1223 if(!after_root){
1224 #ifdef IDA_BND_DEBUG
1225 CONSOLE_DEBUG("Just 'after_root'...");
1226 for(i=0;i<enginedata->nbnds;++i){
1227 relname = bnd_make_name(integ->system, enginedata->bndlist[i]);
1228 CONSOLE_DEBUG("Boundary '%s': bnd_cond_states[%d] = %d"
1229 ,relname,i,bnd_cond_states[i]
1230 );
1231 ASC_FREE(relname);
1232 }
1233 #endif
1234 need_to_reconfigure = ida_cross_boundary(integ, rootsfound,
1235 bnd_cond_states, qrslv_ind, lrslv_ind);
1236 }
1237 if(need_to_reconfigure == 2) {
1238 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Analysis after the boundary failed.");
1239 return 1;
1240 }
1241 if(need_to_reconfigure){
1242 after_root = 1;
1243 if (ida_bnd_update_relist(integ) != 0) {
1244 /* system not square, failure */
1245 return 1;
1246 }
1247 #ifdef IDA_BND_DEBUG
1248 CONSOLE_DEBUG("Boundary(ies) crossed at x=%f: need to reinitialize solver",tret);
1249 #endif
1250 /* so, now we need to restart the integration. we will assume that
1251 everything changes: number of variables, etc, etc, etc. */
1252
1253 /* Need to destroy and rebuild system */
1254 /* IDAFree(ida_mem); */
1255 /* ida_mem = IDACreate(); */
1256 /* ida_reinit_integrator(integ, ida_mem, tout); */
1257 /* n_y may have changed */
1258 N_VDestroy_Serial(yret);
1259 N_VDestroy_Serial(ypret);
1260
1261 yret = N_VNew_Serial(integ->n_y);
1262 ypret = N_VNew_Serial(integ->n_y);
1263 } /* need to reconfigure */
1264
1265 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR>=4
1266 /* set rootdir to -1*rootsfound to set IDA
1267 * to ignore double crossings */
1268 for(i = 0; i < enginedata->nbnds; i++) {
1269 rootdir[i] = -1*rootsfound[i];
1270 #ifdef IDA_BND_DEBUG
1271 char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);
1272 CONSOLE_DEBUG("Set direction=%d for boundary '%s'",rootdir[i],n);
1273 ASC_FREE(n);
1274 #endif
1275 }
1276 IDASetRootDirection(ida_mem, rootdir);\
1277 /*^^^ FIXME what about setting root direction for 'normal' roots? */
1278 #endif
1279 ASC_FREE(rootsfound);
1280 ida_reinit_integrator(integ, ida_mem, tout);
1281 } /* IDA_ROOT_RETURN */
1282
1283 } /* nbnds */
1284
1285 /* Are we sufficiently far from tout to continue
1286 * integrating on this timestep? */
1287 if(fabs(tret - tout) > tol) need_to_reinteg = 1;
1288 else if(after_root){
1289 /* Advance timestep, skip writing output once*/
1290 tout = samplelist_get(integ->samples, t_index + 1);
1291 skipping_output = 1;
1292 }
1293
1294 }while (need_to_reinteg); /* end of solve time step */
1295
1296 if(flag != IDA_ROOT_RETURN && all_bnds_set == 0) {
1297 all_bnds_set = 1;
1298 for(i = 0; i < enginedata->nbnds; i++) {
1299 if(bnd_not_set[i]) {
1300 bnd_cond_states[i] = bndman_calc_satisfied(enginedata->bndlist[i]);
1301 #ifdef IDA_BND_DEBUG
1302 char *n = bnd_make_name(integ->system,enginedata->bndlist[i]);
1303 CONSOLE_DEBUG("Boundary '%s' not set; satisfied=%d",n,bnd_cond_states[i]);
1304 ASC_FREE(n);
1305 #endif
1306 }
1307 }
1308 }
1309
1310 if(!skipping_output){
1311 /* pass the values of everything back to the compiler */
1312 integrator_set_t(integ, (double) tret);
1313 integrator_set_y(integ, NV_DATA_S(yret));
1314 integrator_set_ydot(integ, NV_DATA_S(ypret));
1315
1316 /* -- store the current values of all the stuff */
1317 if(integrator_output_write(integ) == 0) {
1318 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Interrupted at t = %f", tout);
1319 break;
1320 }
1321 integrator_output_write_obs(integ);
1322 }
1323
1324 if(flag < 0) {
1325 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve t = %f (IDASolve), error %d", tout, flag);
1326 break;
1327 }
1328 }/* loop through next sample timestep */
1329
1330 ASC_FREE(rootdir);
1331
1332 /* -- close the IntegratorReporter */
1333 integrator_output_close(integ);
1334
1335 /* get optional outputs */
1336 #ifdef STATS_DEBUG
1337 IntegratorIdaStats stats;
1338 if(IDA_SUCCESS == integrator_ida_stats(ida_mem, &stats)) {
1339 integrator_ida_write_stats(&stats);
1340 }else{
1341 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to fetch stats!?!?");
1342 }
1343 #endif
1344
1345 /* free solution memory */
1346 N_VDestroy_Serial(yret);
1347 N_VDestroy_Serial(ypret);
1348
1349 /* free bnd states if appropriate */
1350 if(enginedata->nbnds) {
1351 ASC_FREE(bnd_cond_states);
1352 }
1353
1354 /* free solver memory */
1355 IDAFree(&ida_mem);
1356
1357 if(flag < -500) {
1358 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Interrupted while attempting t = %f", tout);
1359 return -flag;
1360 }
1361
1362 if(flag < 0) {
1363 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Solving aborted while attempting t = %f", tout);
1364 return 14;
1365 }
1366
1367 /* all done, success */
1368 return 0;
1369 }
1370
1371 /*----------------------------------------------
1372 STATS
1373 */
1374
1375 /**
1376 A simple wrapper to the IDAGetIntegratorStats function. Returns all the
1377 status in a struct instead of separately.
1378
1379 @return IDA_SUCCESS on success.
1380 */
1381 static int integrator_ida_stats(void *ida_mem, IntegratorIdaStats *s) {
1382
1383 #if SUNDIALS_VERSION_MAJOR==2 && SUNDIALS_VERSION_MINOR==2
1384
1385 int res;
1386
1387 /*
1388 There is an error in the documentation for this function in Sundials 2.2.
1389 According the the header file, the hinused stat is not provided.
1390 */
1391 res = IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals,
1392 &s->nlinsetups, &s->netfails, &s->qlast, &s->qcur, &s->hlast,
1393 &s->hcur, &s->tcur);
1394
1395 /* get the missing statistic */
1396 IDAGetActualInitStep(ida_mem, &s->hinused);
1397
1398 return res;
1399 #else
1400
1401 return IDAGetIntegratorStats(ida_mem, &s->nsteps, &s->nrevals, &s->nlinsetups
1402 ,&s->netfails, &s->qlast, &s->qcur, &s->hinused
1403 ,&s->hlast, &s->hcur, &s->tcur
1404 );
1405
1406 #endif
1407 }
1408
1409 /* vim: set ts=4: */

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