1 /* 2 * ----------------------------------------------------------------- 3 * $Revision: 4378 $ 4 * $Date: 2015-02-19 10:55:14 -0800 (Thu, 19 Feb 2015) $ 5 * ----------------------------------------------------------------- 6 * Programmer(s): Radu Serban @ LLNL 7 * ----------------------------------------------------------------- 8 * LLNS Copyright Start 9 * Copyright (c) 2014, Lawrence Livermore National Security 10 * This work was performed under the auspices of the U.S. Department 11 * of Energy by Lawrence Livermore National Laboratory in part under 12 * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344. 13 * Produced at the Lawrence Livermore National Laboratory. 14 * All rights reserved. 15 * For details, see the LICENSE file. 16 * LLNS Copyright End 17 * ----------------------------------------------------------------- 18 * This is the header file (private version) for the main IDAS solver. 19 * ----------------------------------------------------------------- 20 */ 21 22 #ifndef _IDAS_IMPL_H 23 #define _IDAS_IMPL_H 24 25 #include <stdarg.h> 26 27 #include <idas/idas.h> 28 #include <sundials/sundials_nvector.h> 29 #include <sundials/sundials_types.h> 30 31 #ifdef __cplusplus /* wrapper to enable C++ usage */ 32 extern "C" { 33 #endif 34 35 /* 36 * ================================================================= 37 * M A I N I N T E G R A T O R M E M O R Y B L O C K 38 * ================================================================= 39 */ 40 41 42 /* Basic IDA constants */ 43 44 #define HMAX_INV_DEFAULT RCONST(0.0) /* hmax_inv default value */ 45 #define MAXORD_DEFAULT 5 /* maxord default value */ 46 #define MXORDP1 6 /* max. number of N_Vectors in phi */ 47 #define MXSTEP_DEFAULT 500 /* mxstep default value */ 48 49 /* itol */ 50 #define IDA_NN 0 51 #define IDA_SS 1 52 #define IDA_SV 2 53 #define IDA_WF 3 54 #define IDA_EE 4 55 56 /* 57 * ----------------------------------------------------------------- 58 * Types: struct IDAMemRec, IDAMem 59 * ----------------------------------------------------------------- 60 * The type IDAMem is type pointer to struct IDAMemRec. 61 * This structure contains fields to keep track of problem state. 62 * ----------------------------------------------------------------- 63 */ 64 65 typedef struct IDAMemRec { 66 67 realtype ida_uround; /* machine unit roundoff */ 68 69 /*-------------------------- 70 Problem Specification Data 71 --------------------------*/ 72 73 IDAResFn ida_res; /* F(t,y(t),y'(t))=0; the function F */ 74 void *ida_user_data; /* user pointer passed to res */ 75 76 int ida_itol; /* itol = IDA_SS, IDA_SV, IDA_WF, IDA_NN */ 77 realtype ida_rtol; /* relative tolerance */ 78 realtype ida_Satol; /* scalar absolute tolerance */ 79 N_Vector ida_Vatol; /* vector absolute tolerance */ 80 booleantype ida_user_efun; /* TRUE if user provides efun */ 81 IDAEwtFn ida_efun; /* function to set ewt */ 82 void *ida_edata; /* user pointer passed to efun */ 83 84 /*----------------------- 85 Quadrature Related Data 86 -----------------------*/ 87 88 booleantype ida_quadr; 89 90 IDAQuadRhsFn ida_rhsQ; 91 void *ida_user_dataQ; 92 93 booleantype ida_errconQ; 94 95 int ida_itolQ; 96 realtype ida_rtolQ; 97 realtype ida_SatolQ; /* scalar absolute tolerance for quadratures */ 98 N_Vector ida_VatolQ; /* vector absolute tolerance for quadratures */ 99 100 /*------------------------ 101 Sensitivity Related Data 102 ------------------------*/ 103 104 booleantype ida_sensi; 105 int ida_Ns; 106 int ida_ism; 107 108 IDASensResFn ida_resS; 109 void *ida_user_dataS; 110 booleantype ida_resSDQ; 111 112 realtype *ida_p; 113 realtype *ida_pbar; 114 int *ida_plist; 115 int ida_DQtype; 116 realtype ida_DQrhomax; 117 118 booleantype ida_errconS; /* TRUE if sensitivities in err. control */ 119 120 int ida_itolS; 121 realtype ida_rtolS; /* relative tolerance for sensitivities */ 122 realtype *ida_SatolS; /* scalar absolute tolerances for sensi. */ 123 N_Vector *ida_VatolS; /* vector absolute tolerances for sensi. */ 124 125 /*----------------------------------- 126 Quadrature Sensitivity Related Data 127 -----------------------------------*/ 128 129 booleantype ida_quadr_sensi; /* TRUE if computing sensitivities of quadrs.*/ 130 131 IDAQuadSensRhsFn ida_rhsQS; /* fQS = (dfQ/dy)*yS + (dfQ/dp) */ 132 void *ida_user_dataQS; /* data pointer passed to fQS */ 133 booleantype ida_rhsQSDQ; /* TRUE if using internal DQ functions */ 134 135 booleantype ida_errconQS; /* TRUE if yQS are considered in err. con. */ 136 137 int ida_itolQS; 138 realtype ida_rtolQS; /* relative tolerance for yQS */ 139 realtype *ida_SatolQS; /* scalar absolute tolerances for yQS */ 140 N_Vector *ida_VatolQS; /* vector absolute tolerances for yQS */ 141 142 /*----------------------------------------------- 143 Divided differences array and associated arrays 144 -----------------------------------------------*/ 145 146 N_Vector ida_phi[MXORDP1]; /* phi = (maxord+1) arrays of divided differences */ 147 148 realtype ida_psi[MXORDP1]; /* differences in t (sums of recent step sizes) */ 149 realtype ida_alpha[MXORDP1]; /* ratios of current stepsize to psi values */ 150 realtype ida_beta[MXORDP1]; /* ratios of current to previous product of psi's */ 151 realtype ida_sigma[MXORDP1]; /* product successive alpha values and factorial */ 152 realtype ida_gamma[MXORDP1]; /* sum of reciprocals of psi values */ 153 154 /*------------------------- 155 N_Vectors for integration 156 -------------------------*/ 157 158 N_Vector ida_ewt; /* error weight vector */ 159 N_Vector ida_yy; /* work space for y vector (= user's yret) */ 160 N_Vector ida_yp; /* work space for y' vector (= user's ypret) */ 161 N_Vector ida_delta; /* residual vector */ 162 N_Vector ida_id; /* bit vector for diff./algebraic components */ 163 N_Vector ida_constraints; /* vector of inequality constraint options */ 164 N_Vector ida_savres; /* saved residual vector (= tempv1) */ 165 N_Vector ida_ee; /* accumulated corrections to y vector, but 166 set equal to estimated local errors upon 167 successful return */ 168 N_Vector ida_mm; /* mask vector in constraints tests (= tempv2) */ 169 N_Vector ida_tempv1; /* work space vector */ 170 N_Vector ida_tempv2; /* work space vector */ 171 N_Vector ida_ynew; /* work vector for y in IDACalcIC (= tempv2) */ 172 N_Vector ida_ypnew; /* work vector for yp in IDACalcIC (= ee) */ 173 N_Vector ida_delnew; /* work vector for delta in IDACalcIC (= phi[2]) */ 174 N_Vector ida_dtemp; /* work vector in IDACalcIC (= phi[3]) */ 175 176 177 /*---------------------------- 178 Quadrature Related N_Vectors 179 ----------------------------*/ 180 181 N_Vector ida_phiQ[MXORDP1]; 182 N_Vector ida_yyQ; 183 N_Vector ida_ypQ; 184 N_Vector ida_ewtQ; 185 N_Vector ida_eeQ; 186 187 /*--------------------------- 188 Sensitivity Related Vectors 189 ---------------------------*/ 190 191 N_Vector *ida_phiS[MXORDP1]; 192 N_Vector *ida_ewtS; 193 194 N_Vector *ida_eeS; /* cumulative sensitivity corrections */ 195 196 N_Vector *ida_yyS; /* allocated and used for: */ 197 N_Vector *ida_ypS; /* ism = SIMULTANEOUS */ 198 N_Vector *ida_deltaS; /* ism = STAGGERED */ 199 200 N_Vector ida_tmpS1; /* work space vectors | tmpS1 = tempv1 */ 201 N_Vector ida_tmpS2; /* for resS | tmpS2 = tempv2 */ 202 N_Vector ida_tmpS3; /* | tmpS3 = allocated */ 203 204 N_Vector *ida_savresS; /* work vector in IDACalcIC for stg (= phiS[2]) */ 205 N_Vector *ida_delnewS; /* work vector in IDACalcIC for stg (= phiS[3]) */ 206 207 N_Vector *ida_yyS0; /* initial yS, ypS vectors allocated and */ 208 N_Vector *ida_ypS0; /* deallocated in IDACalcIC function */ 209 210 N_Vector *ida_yyS0new; /* work vector in IDASensLineSrch (= phiS[4]) */ 211 N_Vector *ida_ypS0new; /* work vector in IDASensLineSrch (= eeS) */ 212 213 /*-------------------------------------- 214 Quadrature Sensitivity Related Vectors 215 --------------------------------------*/ 216 217 N_Vector *ida_phiQS[MXORDP1];/* Mod. div. diffs. for quadr. sensitivities */ 218 N_Vector *ida_ewtQS; /* error weight vectors for sensitivities */ 219 220 N_Vector *ida_eeQS; /* cumulative quadr.sensi.corrections */ 221 222 N_Vector *ida_yyQS; /* Unlike yS, yQS is not allocated by the user */ 223 N_Vector *ida_tempvQS; /* temporary storage vector (~ tempv) */ 224 N_Vector ida_savrhsQ; /* saved quadr. rhs (needed for rhsQS calls) */ 225 226 /*------------------------------ 227 Variables for use by IDACalcIC 228 ------------------------------*/ 229 230 realtype ida_t0; /* initial t */ 231 N_Vector ida_yy0; /* initial y vector (user-supplied). */ 232 N_Vector ida_yp0; /* initial y' vector (user-supplied). */ 233 234 int ida_icopt; /* IC calculation user option */ 235 booleantype ida_lsoff; /* IC calculation linesearch turnoff option */ 236 int ida_maxnh; /* max. number of h tries in IC calculation */ 237 int ida_maxnj; /* max. number of J tries in IC calculation */ 238 int ida_maxnit; /* max. number of Netwon iterations in IC calc. */ 239 int ida_nbacktr; /* number of IC linesearch backtrack operations */ 240 int ida_sysindex; /* computed system index (0 or 1) */ 241 realtype ida_epiccon; /* IC nonlinear convergence test constant */ 242 realtype ida_steptol; /* minimum Newton step size in IC calculation */ 243 realtype ida_tscale; /* time scale factor = abs(tout1 - t0) */ 244 245 /* Tstop information */ 246 247 booleantype ida_tstopset; 248 realtype ida_tstop; 249 250 /* Step Data */ 251 252 int ida_kk; /* current BDF method order */ 253 int ida_knew; /* order for next step from order decrease decision */ 254 int ida_phase; /* flag to trigger step doubling in first few steps */ 255 int ida_ns; /* counts steps at fixed stepsize and order */ 256 257 realtype ida_hin; /* initial step */ 258 realtype ida_hh; /* current step size h */ 259 realtype ida_rr; /* rr = hnext / hused */ 260 realtype ida_tn; /* current internal value of t */ 261 realtype ida_tretlast; /* value of tret previously returned by IDASolve */ 262 realtype ida_cj; /* current value of scalar (-alphas/hh) in Jacobian */ 263 realtype ida_cjlast; /* cj value saved from last successful step */ 264 realtype ida_cjold; /* cj value saved from last call to lsetup */ 265 realtype ida_cjratio; /* ratio of cj values: cj/cjold */ 266 realtype ida_ss; /* scalar used in Newton iteration convergence test */ 267 realtype ida_epsNewt; /* test constant in Newton convergence test */ 268 realtype ida_epcon; /* coeficient of the Newton covergence test */ 269 realtype ida_toldel; /* tolerance in direct test on Newton corrections */ 270 271 realtype ida_ssS; /* scalar ss for staggered sensitivities */ 272 273 /*------ 274 Limits 275 ------*/ 276 277 int ida_maxncf; /* max numer of convergence failures */ 278 int ida_maxcor; /* max number of Newton corrections */ 279 int ida_maxnef; /* max number of error test failures */ 280 281 int ida_maxord; /* max value of method order k: */ 282 int ida_maxord_alloc; /* value of maxord used when allocating memory */ 283 long int ida_mxstep; /* max number of internal steps for one user call */ 284 realtype ida_hmax_inv; /* inverse of max. step size hmax (default = 0.0) */ 285 286 int ida_maxcorS; /* max number of Newton corrections for sensitivity 287 systems (staggered method) */ 288 289 /*-------- 290 Counters 291 --------*/ 292 293 long int ida_nst; /* number of internal steps taken */ 294 295 long int ida_nre; /* number of function (res) calls */ 296 long int ida_nrQe; 297 long int ida_nrSe; 298 long int ida_nrQSe; /* number of fQS calls */ 299 long int ida_nreS; 300 long int ida_nrQeS; /* number of fQ calls from sensi DQ */ 301 302 303 long int ida_ncfn; /* number of corrector convergence failures */ 304 long int ida_ncfnQ; 305 long int ida_ncfnS; 306 307 long int ida_netf; /* number of error test failures */ 308 long int ida_netfQ; 309 long int ida_netfS; 310 long int ida_netfQS; /* number of quadr. sensi. error test failures */ 311 312 long int ida_nni; /* number of Newton iterations performed */ 313 long int ida_nniS; 314 315 long int ida_nsetups; /* number of lsetup calls */ 316 long int ida_nsetupsS; 317 318 /*--------------------------- 319 Space requirements for IDAS 320 ---------------------------*/ 321 322 long int ida_lrw1; /* no. of realtype words in 1 N_Vector */ 323 long int ida_liw1; /* no. of integer words in 1 N_Vector */ 324 long int ida_lrw1Q; 325 long int ida_liw1Q; 326 long int ida_lrw; /* number of realtype words in IDA work vectors */ 327 long int ida_liw; /* no. of integer words in IDA work vectors */ 328 329 330 /*------------------------------------------- 331 Error handler function and error ouput file 332 -------------------------------------------*/ 333 334 IDAErrHandlerFn ida_ehfun; /* Error messages are handled by ehfun */ 335 void *ida_eh_data; /* dats pointer passed to ehfun */ 336 FILE *ida_errfp; /* IDA error messages are sent to errfp */ 337 338 /* Flags to verify correct calling sequence */ 339 340 booleantype ida_SetupDone; /* set to FALSE by IDAInit and IDAReInit 341 set to TRUE by IDACalcIC or IDASolve */ 342 343 booleantype ida_VatolMallocDone; 344 booleantype ida_constraintsMallocDone; 345 booleantype ida_idMallocDone; 346 347 booleantype ida_MallocDone; /* set to FALSE by IDACreate 348 set to TRUE by IDAInit 349 tested by IDAReInit and IDASolve */ 350 351 booleantype ida_VatolQMallocDone; 352 booleantype ida_quadMallocDone; 353 354 booleantype ida_VatolSMallocDone; 355 booleantype ida_SatolSMallocDone; 356 booleantype ida_sensMallocDone; 357 358 booleantype ida_VatolQSMallocDone; 359 booleantype ida_SatolQSMallocDone; 360 booleantype ida_quadSensMallocDone; 361 362 /*------------------ 363 Linear Solver Data 364 ------------------*/ 365 366 /* Linear Solver functions to be called */ 367 368 int (*ida_linit)(struct IDAMemRec *idamem); 369 370 int (*ida_lsetup)(struct IDAMemRec *idamem, N_Vector yyp, 371 N_Vector ypp, N_Vector resp, 372 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3); 373 374 int (*ida_lsolve)(struct IDAMemRec *idamem, N_Vector b, N_Vector weight, 375 N_Vector ycur, N_Vector ypcur, N_Vector rescur); 376 377 int (*ida_lperf)(struct IDAMemRec *idamem, int perftask); 378 379 int (*ida_lfree)(struct IDAMemRec *idamem); 380 381 /* Linear Solver specific memory */ 382 383 void *ida_lmem; 384 385 /* Flag to request a call to the setup routine */ 386 387 booleantype ida_forceSetup; 388 389 /* Flag to indicate successful ida_linit call */ 390 391 booleantype ida_linitOK; 392 393 /*------------ 394 Saved Values 395 ------------*/ 396 397 booleantype ida_setupNonNull; /* Does setup do something? */ 398 booleantype ida_constraintsSet; /* constraints vector present */ 399 booleantype ida_suppressalg; /* TRUE if suppressing algebraic vars. 400 in local error tests */ 401 int ida_kused; /* method order used on last successful step */ 402 realtype ida_h0u; /* actual initial stepsize */ 403 realtype ida_hused; /* step size used on last successful step */ 404 realtype ida_tolsf; /* tolerance scale factor (saved value) */ 405 406 /*---------------- 407 Rootfinding Data 408 ----------------*/ 409 410 IDARootFn ida_gfun; /* Function g for roots sought */ 411 int ida_nrtfn; /* number of components of g */ 412 int *ida_iroots; /* array for root information */ 413 int *ida_rootdir; /* array specifying direction of zero-crossing */ 414 realtype ida_tlo; /* nearest endpoint of interval in root search */ 415 realtype ida_thi; /* farthest endpoint of interval in root search */ 416 realtype ida_trout; /* t return value from rootfinder routine */ 417 realtype *ida_glo; /* saved array of g values at t = tlo */ 418 realtype *ida_ghi; /* saved array of g values at t = thi */ 419 realtype *ida_grout; /* array of g values at t = trout */ 420 realtype ida_toutc; /* copy of tout (if NORMAL mode) */ 421 realtype ida_ttol; /* tolerance on root location */ 422 int ida_taskc; /* copy of parameter itask */ 423 int ida_irfnd; /* flag showing whether last step had a root */ 424 long int ida_nge; /* counter for g evaluations */ 425 booleantype *ida_gactive; /* array with active/inactive event functions */ 426 int ida_mxgnull; /* number of warning messages about possible g==0 */ 427 428 /*------------------------ 429 Adjoint sensitivity data 430 ------------------------*/ 431 432 booleantype ida_adj; /* TRUE if performing ASA */ 433 434 struct IDAadjMemRec *ida_adj_mem; /* Pointer to adjoint memory structure */ 435 436 booleantype ida_adjMallocDone; 437 438 } *IDAMem; 439 440 /* 441 * ================================================================= 442 * A D J O I N T M O D U L E M E M O R Y B L O C K 443 * ================================================================= 444 */ 445 446 /* 447 * ----------------------------------------------------------------- 448 * Forward references for pointers to various structures 449 * ----------------------------------------------------------------- 450 */ 451 452 typedef struct IDAadjMemRec *IDAadjMem; 453 typedef struct CkpntMemRec *CkpntMem; 454 typedef struct DtpntMemRec *DtpntMem; 455 typedef struct IDABMemRec *IDABMem; 456 457 /* 458 * ----------------------------------------------------------------- 459 * Types for functions provided by an interpolation module 460 * ----------------------------------------------------------------- 461 * IDAAMMallocFn: Type for a function that initializes the content 462 * field of the structures in the dt array 463 * IDAAMFreeFn: Type for a function that deallocates the content 464 * field of the structures in the dt array 465 * IDAAGetYFn: Function type for a function that returns the 466 * interpolated forward solution. 467 * IDAAStorePnt: Function type for a function that stores a new 468 * point in the structure d 469 * ----------------------------------------------------------------- 470 */ 471 472 typedef booleantype (*IDAAMMallocFn)(IDAMem IDA_mem); 473 typedef void (*IDAAMFreeFn)(IDAMem IDA_mem); 474 typedef int (*IDAAGetYFn)(IDAMem IDA_mem, realtype t, 475 N_Vector yy, N_Vector yp, 476 N_Vector *yyS, N_Vector *ypS); 477 typedef int (*IDAAStorePntFn)(IDAMem IDA_mem, DtpntMem d); 478 479 /* 480 * ----------------------------------------------------------------- 481 * Types : struct CkpntMemRec, CkpntMem 482 * ----------------------------------------------------------------- 483 * The type CkpntMem is type pointer to struct CkpntMemRec. 484 * This structure contains fields to store all information at a 485 * check point that is needed to 'hot' start IDAS. 486 * ----------------------------------------------------------------- 487 */ 488 489 struct CkpntMemRec { 490 491 /* Integration limits */ 492 realtype ck_t0; 493 realtype ck_t1; 494 495 /* Modified divided difference array */ 496 N_Vector ck_phi[MXORDP1]; 497 498 /* Do we need to carry quadratures? */ 499 booleantype ck_quadr; 500 501 /* Modified divided difference array for quadratures */ 502 N_Vector ck_phiQ[MXORDP1]; 503 504 /* Do we need to carry sensitivities? */ 505 booleantype ck_sensi; 506 507 /* number of sensitivities */ 508 int ck_Ns; 509 510 /* Modified divided difference array for sensitivities */ 511 N_Vector *ck_phiS[MXORDP1]; 512 513 /* Do we need to carry quadrature sensitivities? */ 514 booleantype ck_quadr_sensi; 515 516 /* Modified divided difference array for quadrature sensitivities */ 517 N_Vector *ck_phiQS[MXORDP1]; 518 519 520 /* Step data */ 521 long int ck_nst; 522 realtype ck_tretlast; 523 long int ck_ns; 524 int ck_kk; 525 int ck_kused; 526 int ck_knew; 527 int ck_phase; 528 529 realtype ck_hh; 530 realtype ck_hused; 531 realtype ck_rr; 532 realtype ck_cj; 533 realtype ck_cjlast; 534 realtype ck_cjold; 535 realtype ck_cjratio; 536 realtype ck_ss; 537 realtype ck_ssS; 538 539 realtype ck_psi[MXORDP1]; 540 realtype ck_alpha[MXORDP1]; 541 realtype ck_beta[MXORDP1]; 542 realtype ck_sigma[MXORDP1]; 543 realtype ck_gamma[MXORDP1]; 544 545 /* How many phi, phiS, phiQ and phiQS were allocated? */ 546 int ck_phi_alloc; 547 548 /* Pointer to next structure in list */ 549 struct CkpntMemRec *ck_next; 550 }; 551 552 /* 553 * ----------------------------------------------------------------- 554 * Type : struct DtpntMemRec 555 * ----------------------------------------------------------------- 556 * This structure contains fields to store all information at a 557 * data point that is needed to interpolate solution of forward 558 * simulations. Its content field is interpType-dependent. 559 * ----------------------------------------------------------------- 560 */ 561 562 struct DtpntMemRec { 563 realtype t; /* time */ 564 void *content; /* interpType-dependent content */ 565 }; 566 567 /* Data for cubic Hermite interpolation */ 568 typedef struct HermiteDataMemRec { 569 N_Vector y; 570 N_Vector yd; 571 N_Vector *yS; 572 N_Vector *ySd; 573 } *HermiteDataMem; 574 575 /* Data for polynomial interpolation */ 576 typedef struct PolynomialDataMemRec { 577 N_Vector y; 578 N_Vector *yS; 579 580 /* yd and ySd store the derivative(s) only for the first dt 581 point. NULL otherwise. */ 582 N_Vector yd; 583 N_Vector *ySd; 584 int order; 585 } *PolynomialDataMem; 586 587 /* 588 * ----------------------------------------------------------------- 589 * Type : struct IDABMemRec 590 * ----------------------------------------------------------------- 591 * The type IDABMemRec is a pointer to a structure which stores all 592 * information for ONE backward problem. 593 * The IDAadjMem struct contains a linked list of IDABMem pointers 594 * ----------------------------------------------------------------- 595 */ 596 struct IDABMemRec { 597 598 /* Index of this backward problem */ 599 int ida_index; 600 601 /* Time at which the backward problem is initialized. */ 602 realtype ida_t0; 603 604 /* Memory for this backward problem */ 605 IDAMem IDA_mem; 606 607 /* Flags to indicate that this backward problem's RHS or quad RHS 608 * require forward sensitivities */ 609 booleantype ida_res_withSensi; 610 booleantype ida_rhsQ_withSensi; 611 612 /* Residual function for backward run */ 613 IDAResFnB ida_res; 614 IDAResFnBS ida_resS; 615 616 /* Right hand side quadrature function (fQB) for backward run */ 617 IDAQuadRhsFnB ida_rhsQ; 618 IDAQuadRhsFnBS ida_rhsQS; 619 620 /* User user_data */ 621 void *ida_user_data; 622 623 /* Linear solver's data and functions */ 624 625 /* Memory block for a linear solver's interface to IDAA */ 626 void *ida_lmem; 627 628 /* Function to free any memory allocated by the linear solver */ 629 void (*ida_lfree)(IDABMem IDAB_mem); 630 631 /* Memory block for a preconditioner's module interface to IDAA */ 632 void *ida_pmem; 633 634 /* Function to free any memory allocated by the preconditioner module */ 635 void (*ida_pfree)(IDABMem IDAB_mem); 636 637 /* Time at which to extract solution / quadratures */ 638 realtype ida_tout; 639 640 /* Workspace Nvectors */ 641 N_Vector ida_yy; 642 N_Vector ida_yp; 643 644 /* Link to next structure in list. */ 645 struct IDABMemRec *ida_next; 646 }; 647 648 649 /* 650 * ----------------------------------------------------------------- 651 * Type : struct IDAadjMemRec 652 * ----------------------------------------------------------------- 653 * The type IDAadjMem is type pointer to struct IDAadjMemRec. 654 * This structure contins fields to store all information 655 * necessary for adjoint sensitivity analysis. 656 * ----------------------------------------------------------------- 657 */ 658 659 struct IDAadjMemRec { 660 661 /* -------------------- 662 * Forward problem data 663 * -------------------- */ 664 665 /* Integration interval */ 666 realtype ia_tinitial, ia_tfinal; 667 668 /* Flag for first call to IDASolveF */ 669 booleantype ia_firstIDAFcall; 670 671 /* Flag if IDASolveF was called with TSTOP */ 672 booleantype ia_tstopIDAFcall; 673 realtype ia_tstopIDAF; 674 675 /* ---------------------- 676 * Backward problems data 677 * ---------------------- */ 678 679 /* Storage for backward problems */ 680 struct IDABMemRec *IDAB_mem; 681 682 /* Number of backward problems. */ 683 int ia_nbckpbs; 684 685 /* Address of current backward problem (iterator). */ 686 struct IDABMemRec *ia_bckpbCrt; 687 688 /* Flag for first call to IDASolveB */ 689 booleantype ia_firstIDABcall; 690 691 /* ---------------- 692 * Check point data 693 * ---------------- */ 694 695 /* Storage for check point information */ 696 struct CkpntMemRec *ck_mem; 697 698 /* address of the check point structure for which data is available */ 699 struct CkpntMemRec *ia_ckpntData; 700 701 /* Number of checkpoints. */ 702 int ia_nckpnts; 703 704 /* ------------------ 705 * Interpolation data 706 * ------------------ */ 707 708 /* Number of steps between 2 check points */ 709 long int ia_nsteps; 710 711 /* Storage for data from forward runs */ 712 struct DtpntMemRec **dt_mem; 713 714 /* Actual number of data points saved in current dt_mem */ 715 /* Commonly, np = nsteps+1 */ 716 long int ia_np; 717 718 /* Interpolation type */ 719 int ia_interpType; 720 721 722 /* Functions set by the interpolation module */ 723 IDAAStorePntFn ia_storePnt; /* store a new interpolation point */ 724 IDAAGetYFn ia_getY; /* interpolate forward solution */ 725 IDAAMMallocFn ia_malloc; /* allocate new data point */ 726 IDAAMFreeFn ia_free; /* destroys data point */ 727 728 /* Flags controlling the interpolation module */ 729 booleantype ia_mallocDone; /* IM initialized? */ 730 booleantype ia_newData; /* new data available in dt_mem? */ 731 booleantype ia_storeSensi; /* store sensitivities? */ 732 booleantype ia_interpSensi; /* interpolate sensitivities? */ 733 734 booleantype ia_noInterp; /* interpolations are temporarly */ 735 /* disabled ( IDACalcICB ) */ 736 737 /* Workspace for polynomial interpolation */ 738 N_Vector ia_Y[MXORDP1]; /* pointers phi[i] */ 739 N_Vector *ia_YS[MXORDP1]; /* pointers phiS[i] */ 740 realtype ia_T[MXORDP1]; 741 742 /* Workspace for wrapper functions */ 743 N_Vector ia_yyTmp, ia_ypTmp; 744 N_Vector *ia_yySTmp, *ia_ypSTmp; 745 746 long int ilast; 747 748 }; 749 750 751 /* 752 * ================================================================= 753 * I N T E R F A C E T O L I N E A R S O L V E R S 754 * ================================================================= 755 */ 756 757 /* 758 * ----------------------------------------------------------------- 759 * int (*ida_linit)(IDAMem IDA_mem); 760 * ----------------------------------------------------------------- 761 * The purpose of ida_linit is to allocate memory for the 762 * solver-specific fields in the structure *(idamem->ida_lmem) and 763 * perform any needed initializations of solver-specific memory, 764 * such as counters/statistics. An (*ida_linit) should return 765 * 0 if it has successfully initialized the IDA linear solver and 766 * a non-zero value otherwise. If an error does occur, an 767 * appropriate message should be issued. 768 * ---------------------------------------------------------------- 769 */ 770 771 /* 772 * ----------------------------------------------------------------- 773 * int (*ida_lsetup)(IDAMem IDA_mem, N_Vector yyp, N_Vector ypp, 774 * N_Vector resp, 775 * N_Vector tempv1, N_Vector tempv2, N_Vector tempv3); 776 * ----------------------------------------------------------------- 777 * The job of ida_lsetup is to prepare the linear solver for 778 * subsequent calls to ida_lsolve. Its parameters are as follows: 779 * 780 * idamem - problem memory pointer of type IDAMem. See the big 781 * typedef earlier in this file. 782 * 783 * 784 * yyp - the predicted y vector for the current IDA internal 785 * step. 786 * 787 * ypp - the predicted y' vector for the current IDA internal 788 * step. 789 * 790 * resp - F(tn, yyp, ypp). 791 * 792 * tempv1, tempv2, tempv3 - temporary N_Vectors provided for use 793 * by ida_lsetup. 794 * 795 * The ida_lsetup routine should return 0 if successful, 796 * a positive value for a recoverable error, and a negative value 797 * for an unrecoverable error. 798 * ----------------------------------------------------------------- 799 */ 800 801 /* 802 * ----------------------------------------------------------------- 803 * int (*ida_lsolve)(IDAMem IDA_mem, N_Vector b, N_Vector weight, 804 * N_Vector ycur, N_Vector ypcur, N_Vector rescur); 805 * ----------------------------------------------------------------- 806 * ida_lsolve must solve the linear equation P x = b, where 807 * P is some approximation to the system Jacobian 808 * J = (dF/dy) + cj (dF/dy') 809 * evaluated at (tn,ycur,ypcur) and the RHS vector b is input. 810 * The N-vector ycur contains the solver's current approximation 811 * to y(tn), ypcur contains that for y'(tn), and the vector rescur 812 * contains the N-vector residual F(tn,ycur,ypcur). 813 * The solution is to be returned in the vector b. 814 * 815 * The ida_lsolve routine should return 0 if successful, 816 * a positive value for a recoverable error, and a negative value 817 * for an unrecoverable error. 818 * ----------------------------------------------------------------- 819 */ 820 821 /* 822 * ----------------------------------------------------------------- 823 * int (*ida_lperf)(IDAMem IDA_mem, int perftask); 824 * ----------------------------------------------------------------- 825 * ida_lperf is called two places in IDAS where linear solver 826 * performance data is required by IDAS. For perftask = 0, an 827 * initialization of performance variables is performed, while for 828 * perftask = 1, the performance is evaluated. 829 * ----------------------------------------------------------------- 830 */ 831 832 /* 833 * ================================================================= 834 * I D A S I N T E R N A L F U N C T I O N S 835 * ================================================================= 836 */ 837 838 /* Prototype of internal ewtSet function */ 839 840 int IDAEwtSet(N_Vector ycur, N_Vector weight, void *data); 841 842 /* High level error handler */ 843 844 void IDAProcessError(IDAMem IDA_mem, 845 int error_code, const char *module, const char *fname, 846 const char *msgfmt, ...); 847 848 /* Prototype of internal errHandler function */ 849 850 void IDAErrHandler(int error_code, const char *module, const char *function, 851 char *msg, void *data); 852 853 /* Prototype for internal sensitivity residual DQ function */ 854 855 int IDASensResDQ(int Ns, realtype t, 856 N_Vector yy, N_Vector yp, N_Vector resval, 857 N_Vector *yyS, N_Vector *ypS, N_Vector *resvalS, 858 void *user_dataS, 859 N_Vector ytemp, N_Vector yptemp, N_Vector restemp); 860 861 /* 862 * ================================================================= 863 * I D A S E R R O R M E S S A G E S 864 * ================================================================= 865 */ 866 867 #if defined(SUNDIALS_EXTENDED_PRECISION) 868 869 #define MSG_TIME "t = %Lg, " 870 #define MSG_TIME_H "t = %Lg and h = %Lg, " 871 #define MSG_TIME_INT "t = %Lg is not between tcur - hu = %Lg and tcur = %Lg." 872 #define MSG_TIME_TOUT "tout = %Lg" 873 #define MSG_TIME_TSTOP "tstop = %Lg" 874 875 #elif defined(SUNDIALS_DOUBLE_PRECISION) 876 877 #define MSG_TIME "t = %lg, " 878 #define MSG_TIME_H "t = %lg and h = %lg, " 879 #define MSG_TIME_INT "t = %lg is not between tcur - hu = %lg and tcur = %lg." 880 #define MSG_TIME_TOUT "tout = %lg" 881 #define MSG_TIME_TSTOP "tstop = %lg" 882 883 #else 884 885 #define MSG_TIME "t = %g, " 886 #define MSG_TIME_H "t = %g and h = %g, " 887 #define MSG_TIME_INT "t = %g is not between tcur - hu = %g and tcur = %g." 888 #define MSG_TIME_TOUT "tout = %g" 889 #define MSG_TIME_TSTOP "tstop = %g" 890 891 #endif 892 893 /* General errors */ 894 895 #define MSG_MEM_FAIL "A memory request failed." 896 #define MSG_NO_MEM "ida_mem = NULL illegal." 897 #define MSG_NO_MALLOC "Attempt to call before IDAMalloc." 898 #define MSG_BAD_NVECTOR "A required vector operation is not implemented." 899 900 /* Initialization errors */ 901 902 #define MSG_Y0_NULL "y0 = NULL illegal." 903 #define MSG_YP0_NULL "yp0 = NULL illegal." 904 #define MSG_BAD_ITOL "Illegal value for itol. The legal values are IDA_SS, IDA_SV, and IDA_WF." 905 #define MSG_RES_NULL "res = NULL illegal." 906 #define MSG_BAD_RTOL "rtol < 0 illegal." 907 #define MSG_ATOL_NULL "atol = NULL illegal." 908 #define MSG_BAD_ATOL "Some atol component < 0.0 illegal." 909 #define MSG_ROOT_FUNC_NULL "g = NULL illegal." 910 911 #define MSG_MISSING_ID "id = NULL but suppressalg option on." 912 #define MSG_NO_TOLS "No integration tolerances have been specified." 913 #define MSG_FAIL_EWT "The user-provide EwtSet function failed." 914 #define MSG_BAD_EWT "Some initial ewt component = 0.0 illegal." 915 #define MSG_Y0_FAIL_CONSTR "y0 fails to satisfy constraints." 916 #define MSG_BAD_ISM_CONSTR "Constraints can not be enforced while forward sensitivity is used with simultaneous method." 917 #define MSG_LSOLVE_NULL "The linear solver's solve routine is NULL." 918 #define MSG_LINIT_FAIL "The linear solver's init routine failed." 919 920 #define MSG_NO_QUAD "Illegal attempt to call before calling IDAQuadInit." 921 #define MSG_BAD_EWTQ "Initial ewtQ has component(s) equal to zero (illegal)." 922 #define MSG_BAD_ITOLQ "Illegal value for itolQ. The legal values are IDA_SS and IDA_SV." 923 #define MSG_NO_TOLQ "No integration tolerances for quadrature variables have been specified." 924 #define MSG_NULL_ATOLQ "atolQ = NULL illegal." 925 #define MSG_BAD_RTOLQ "rtolQ < 0 illegal." 926 #define MSG_BAD_ATOLQ "atolQ has negative component(s) (illegal)." 927 928 #define MSG_NO_SENSI "Illegal attempt to call before calling IDASensInit." 929 #define MSG_BAD_EWTS "Initial ewtS has component(s) equal to zero (illegal)." 930 #define MSG_BAD_ITOLS "Illegal value for itolS. The legal values are IDA_SS, IDA_SV, and IDA_EE." 931 #define MSG_NULL_ATOLS "atolS = NULL illegal." 932 #define MSG_BAD_RTOLS "rtolS < 0 illegal." 933 #define MSG_BAD_ATOLS "atolS has negative component(s) (illegal)." 934 #define MSG_BAD_PBAR "pbar has zero component(s) (illegal)." 935 #define MSG_BAD_PLIST "plist has negative component(s) (illegal)." 936 #define MSG_BAD_NS "NS <= 0 illegal." 937 #define MSG_NULL_YYS0 "yyS0 = NULL illegal." 938 #define MSG_NULL_YPS0 "ypS0 = NULL illegal." 939 #define MSG_BAD_ISM "Illegal value for ism. Legal values are: IDA_SIMULTANEOUS and IDA_STAGGERED." 940 #define MSG_BAD_IS "Illegal value for is." 941 #define MSG_NULL_DKYA "dkyA = NULL illegal." 942 #define MSG_BAD_DQTYPE "Illegal value for DQtype. Legal values are: IDA_CENTERED and IDA_FORWARD." 943 #define MSG_BAD_DQRHO "DQrhomax < 0 illegal." 944 945 #define MSG_NULL_ABSTOLQS "abstolQS = NULL illegal parameter." 946 #define MSG_BAD_RELTOLQS "reltolQS < 0 illegal parameter." 947 #define MSG_BAD_ABSTOLQS "abstolQS has negative component(s) (illegal)." 948 #define MSG_NO_QUADSENSI "Forward sensitivity analysis for quadrature variables was not activated." 949 #define MSG_NULL_YQS0 "yQS0 = NULL illegal parameter." 950 951 952 /* IDACalcIC error messages */ 953 954 #define MSG_IC_BAD_ICOPT "icopt has an illegal value." 955 #define MSG_IC_MISSING_ID "id = NULL conflicts with icopt." 956 #define MSG_IC_TOO_CLOSE "tout1 too close to t0 to attempt initial condition calculation." 957 #define MSG_IC_BAD_ID "id has illegal values." 958 #define MSG_IC_BAD_EWT "Some initial ewt component = 0.0 illegal." 959 #define MSG_IC_RES_NONREC "The residual function failed unrecoverably. " 960 #define MSG_IC_RES_FAIL "The residual function failed at the first call. " 961 #define MSG_IC_SETUP_FAIL "The linear solver setup failed unrecoverably." 962 #define MSG_IC_SOLVE_FAIL "The linear solver solve failed unrecoverably." 963 #define MSG_IC_NO_RECOVERY "The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover." 964 #define MSG_IC_FAIL_CONSTR "Unable to satisfy the inequality constraints." 965 #define MSG_IC_FAILED_LINS "The linesearch algorithm failed with too small a step." 966 #define MSG_IC_CONV_FAILED "Newton/Linesearch algorithm failed to converge." 967 968 /* IDASolve error messages */ 969 970 #define MSG_YRET_NULL "yret = NULL illegal." 971 #define MSG_YPRET_NULL "ypret = NULL illegal." 972 #define MSG_TRET_NULL "tret = NULL illegal." 973 #define MSG_BAD_ITASK "itask has an illegal value." 974 #define MSG_TOO_CLOSE "tout too close to t0 to start integration." 975 #define MSG_BAD_HINIT "Initial step is not towards tout." 976 #define MSG_BAD_TSTOP "The value " MSG_TIME_TSTOP " is behind current " MSG_TIME "in the direction of integration." 977 #define MSG_CLOSE_ROOTS "Root found at and very near " MSG_TIME "." 978 #define MSG_MAX_STEPS "At " MSG_TIME ", mxstep steps taken before reaching tout." 979 #define MSG_EWT_NOW_FAIL "At " MSG_TIME "the user-provide EwtSet function failed." 980 #define MSG_EWT_NOW_BAD "At " MSG_TIME "some ewt component has become <= 0.0." 981 #define MSG_TOO_MUCH_ACC "At " MSG_TIME "too much accuracy requested." 982 983 #define MSG_BAD_T "Illegal value for t. " MSG_TIME_INT 984 #define MSG_BAD_TOUT "Trouble interpolating at " MSG_TIME_TOUT ". tout too far back in direction of integration." 985 986 #define MSG_BAD_K "Illegal value for k." 987 #define MSG_NULL_DKY "dky = NULL illegal." 988 #define MSG_NULL_DKYP "dkyp = NULL illegal." 989 990 #define MSG_ERR_FAILS "At " MSG_TIME_H "the error test failed repeatedly or with |h| = hmin." 991 #define MSG_CONV_FAILS "At " MSG_TIME_H "the corrector convergence failed repeatedly or with |h| = hmin." 992 #define MSG_SETUP_FAILED "At " MSG_TIME "the linear solver setup failed unrecoverably." 993 #define MSG_SOLVE_FAILED "At " MSG_TIME "the linear solver solve failed unrecoverably." 994 #define MSG_REP_RES_ERR "At " MSG_TIME "repeated recoverable residual errors." 995 #define MSG_RES_NONRECOV "At " MSG_TIME "the residual function failed unrecoverably." 996 #define MSG_FAILED_CONSTR "At " MSG_TIME "unable to satisfy inequality constraints." 997 #define MSG_RTFUNC_FAILED "At " MSG_TIME ", the rootfinding routine failed in an unrecoverable manner." 998 #define MSG_NO_ROOT "Rootfinding was not initialized." 999 #define MSG_INACTIVE_ROOTS "At the end of the first step, there are still some root functions identically 0. This warning will not be issued again." 1000 1001 #define MSG_EWTQ_NOW_BAD "At " MSG_TIME ", a component of ewtQ has become <= 0." 1002 #define MSG_QRHSFUNC_FAILED "At " MSG_TIME ", the quadrature right-hand side routine failed in an unrecoverable manner." 1003 #define MSG_QRHSFUNC_UNREC "At " MSG_TIME ", the quadrature right-hand side failed in a recoverable manner, but no recovery is possible." 1004 #define MSG_QRHSFUNC_REPTD "At " MSG_TIME "repeated recoverable quadrature right-hand side function errors." 1005 #define MSG_QRHSFUNC_FIRST "The quadrature right-hand side routine failed at the first call." 1006 1007 #define MSG_NULL_P "p = NULL when using internal DQ for sensitivity residual is illegal." 1008 #define MSG_EWTS_NOW_BAD "At " MSG_TIME ", a component of ewtS has become <= 0." 1009 #define MSG_SRHSFUNC_FAILED "At " MSG_TIME ", the sensitivity residual routine failed in an unrecoverable manner." 1010 #define MSG_SRHSFUNC_UNREC "At " MSG_TIME ", the sensitivity residual failed in a recoverable manner, but no recovery is possible." 1011 #define MSG_SRHSFUNC_REPTD "At " MSG_TIME "repeated recoverable sensitivity residual function errors." 1012 1013 #define MSG_NO_TOLQS "No integration tolerances for quadrature sensitivity variables have been specified." 1014 #define MSG_NULL_RHSQ "IDAS is expected to use DQ to evaluate the RHS of quad. sensi., but quadratures were not initialized." 1015 #define MSG_BAD_EWTQS "Initial ewtQS has component(s) equal to zero (illegal)." 1016 #define MSG_EWTQS_NOW_BAD "At " MSG_TIME ", a component of ewtQS has become <= 0." 1017 #define MSG_QSRHSFUNC_FAILED "At " MSG_TIME ", the sensitivity quadrature right-hand side routine failed in an unrecoverable manner." 1018 #define MSG_QSRHSFUNC_FIRST "The quadrature right-hand side routine failed at the first call." 1019 1020 /* IDASet* / IDAGet* error messages */ 1021 #define MSG_NEG_MAXORD "maxord<=0 illegal." 1022 #define MSG_BAD_MAXORD "Illegal attempt to increase maximum order." 1023 #define MSG_NEG_HMAX "hmax < 0 illegal." 1024 #define MSG_NEG_EPCON "epcon <= 0.0 illegal." 1025 #define MSG_BAD_CONSTR "Illegal values in constraints vector." 1026 #define MSG_BAD_EPICCON "epiccon <= 0.0 illegal." 1027 #define MSG_BAD_MAXNH "maxnh <= 0 illegal." 1028 #define MSG_BAD_MAXNJ "maxnj <= 0 illegal." 1029 #define MSG_BAD_MAXNIT "maxnit <= 0 illegal." 1030 #define MSG_BAD_STEPTOL "steptol <= 0.0 illegal." 1031 1032 #define MSG_TOO_LATE "IDAGetConsistentIC can only be called before IDASolve." 1033 1034 /* 1035 * ================================================================= 1036 * I D A A E R R O R M E S S A G E S 1037 * ================================================================= 1038 */ 1039 1040 #define MSGAM_NULL_IDAMEM "ida_mem = NULL illegal." 1041 #define MSGAM_NO_ADJ "Illegal attempt to call before calling IDAadjInit." 1042 #define MSGAM_BAD_INTERP "Illegal value for interp." 1043 #define MSGAM_BAD_STEPS "Steps nonpositive illegal." 1044 #define MSGAM_BAD_WHICH "Illegal value for which." 1045 #define MSGAM_NO_BCK "No backward problems have been defined yet." 1046 #define MSGAM_NO_FWD "Illegal attempt to call before calling IDASolveF." 1047 #define MSGAM_BAD_TB0 "The initial time tB0 is outside the interval over which the forward problem was solved." 1048 #define MSGAM_BAD_SENSI "At least one backward problem requires sensitivities, but they were not stored for interpolation." 1049 #define MSGAM_BAD_ITASKB "Illegal value for itaskB. Legal values are IDA_NORMAL and IDA_ONE_STEP." 1050 #define MSGAM_BAD_TBOUT "The final time tBout is outside the interval over which the forward problem was solved." 1051 #define MSGAM_BACK_ERROR "Error occured while integrating backward problem # %d" 1052 #define MSGAM_BAD_TINTERP "Bad t = %g for interpolation." 1053 #define MSGAM_BAD_T "Bad t for interpolation." 1054 #define MSGAM_WRONG_INTERP "This function cannot be called for the specified interp type." 1055 #define MSGAM_MEM_FAIL "A memory request failed." 1056 #define MSGAM_NO_INITBS "Illegal attempt to call before calling IDAInitBS." 1057 1058 #ifdef __cplusplus 1059 } 1060 #endif 1061 1062 #endif 1063