1 /*--------------------------------------------------------------- 2 * Programmer(s): Daniel R. Reynolds @ SMU 3 *--------------------------------------------------------------- 4 * SUNDIALS Copyright Start 5 * Copyright (c) 2002-2020, Lawrence Livermore National Security 6 * and Southern Methodist University. 7 * All rights reserved. 8 * 9 * See the top-level LICENSE and NOTICE files for details. 10 * 11 * SPDX-License-Identifier: BSD-3-Clause 12 * SUNDIALS Copyright End 13 *--------------------------------------------------------------- 14 * Implementation header file for ARKode's ARK time stepper 15 * module. 16 *--------------------------------------------------------------*/ 17 18 #ifndef _ARKODE_ARKSTEP_IMPL_H 19 #define _ARKODE_ARKSTEP_IMPL_H 20 21 #include <arkode/arkode_arkstep.h> 22 #include "arkode_impl.h" 23 #include "arkode_ls_impl.h" 24 25 #ifdef __cplusplus /* wrapper to enable C++ usage */ 26 extern "C" { 27 #endif 28 29 /*=============================================================== 30 ARK time step module constants 31 ===============================================================*/ 32 33 #define MAXCOR 3 /* max number of nonlinear iterations */ 34 #define CRDOWN RCONST(0.3) /* constant to estimate the convergence 35 rate for the nonlinear equation */ 36 #define DGMAX RCONST(0.2) /* if |gamma/gammap-1| > DGMAX then call lsetup */ 37 #define RDIV RCONST(2.3) /* declare divergence if ratio del/delp > RDIV */ 38 #define MSBP 20 /* max no. of steps between lsetup calls */ 39 40 /* Default solver tolerance factor */ 41 /* #define NLSCOEF RCONST(0.003) */ /* Hairer & Wanner constant */ 42 /* #define NLSCOEF RCONST(0.2) */ /* CVODE constant */ 43 #define NLSCOEF RCONST(0.1) 44 45 46 47 /*=============================================================== 48 ARK time step module data structure 49 ===============================================================*/ 50 51 /*--------------------------------------------------------------- 52 Types : struct ARKodeARKStepMemRec, ARKodeARKStepMem 53 --------------------------------------------------------------- 54 The type ARKodeARKStepMem is type pointer to struct 55 ARKodeARKStepMemRec. This structure contains fields to 56 perform an additive Runge-Kutta time step. 57 ---------------------------------------------------------------*/ 58 typedef struct ARKodeARKStepMemRec { 59 60 /* ARK problem specification */ 61 ARKRhsFn fe; /* My' = fe(t,y) + fi(t,y) */ 62 ARKRhsFn fi; 63 booleantype linear; /* SUNTRUE if fi is linear */ 64 booleantype linear_timedep; /* SUNTRUE if dfi/dy depends on t */ 65 booleantype explicit; /* SUNTRUE if fe is enabled */ 66 booleantype implicit; /* SUNTRUE if fi is enabled */ 67 68 /* ARK method storage and parameters */ 69 N_Vector *Fe; /* explicit RHS at each stage */ 70 N_Vector *Fi; /* implicit RHS at each stage */ 71 N_Vector sdata; /* old stage data in residual */ 72 N_Vector zpred; /* predicted stage solution */ 73 N_Vector zcor; /* stage correction */ 74 int q; /* method order */ 75 int p; /* embedding order */ 76 int istage; /* current stage */ 77 int stages; /* number of stages */ 78 ARKodeButcherTable Be; /* ERK Butcher table */ 79 ARKodeButcherTable Bi; /* IRK Butcher table */ 80 81 /* User-supplied stage predictor routine */ 82 ARKStepStagePredictFn stage_predict; 83 84 /* (Non)Linear solver parameters & data */ 85 SUNNonlinearSolver NLS; /* generic SUNNonlinearSolver object */ 86 booleantype ownNLS; /* flag indicating ownership of NLS */ 87 realtype gamma; /* gamma = h * A(i,i) */ 88 realtype gammap; /* gamma at the last setup call */ 89 realtype gamrat; /* gamma / gammap */ 90 realtype dgmax; /* call lsetup if |gamma/gammap-1| >= dgmax */ 91 92 int predictor; /* implicit prediction method to use */ 93 realtype crdown; /* nonlinear conv rate estimation constant */ 94 realtype rdiv; /* nonlin divergence if del/delp > rdiv */ 95 realtype crate; /* estimated nonlin convergence rate */ 96 realtype delp; /* norm of previous nonlinear solver update */ 97 realtype eRNrm; /* estimated residual norm, used in nonlin 98 and linear solver convergence tests */ 99 realtype nlscoef; /* coefficient in nonlin. convergence test */ 100 int mnewt; /* internal Newton iteration counter */ 101 102 int msbp; /* positive => max # steps between lsetup 103 negative => call at each Newton iter */ 104 long int nstlp; /* step number of last setup call */ 105 106 int maxcor; /* max num iterations for solving the 107 nonlinear equation */ 108 109 int convfail; /* NLS fail flag (for interface routines) */ 110 booleantype jcur; /* is Jacobian info for lin solver current? */ 111 112 /* Linear Solver Data */ 113 ARKLinsolInitFn linit; 114 ARKLinsolSetupFn lsetup; 115 ARKLinsolSolveFn lsolve; 116 ARKLinsolFreeFn lfree; 117 void *lmem; 118 int lsolve_type; /* interface type: 0=iterative; 1=direct; 2=custom */ 119 120 /* Mass matrix solver data */ 121 ARKMassInitFn minit; 122 ARKMassSetupFn msetup; 123 ARKMassMultFn mmult; 124 ARKMassSolveFn msolve; 125 ARKMassFreeFn mfree; 126 void* mass_mem; 127 realtype msetuptime; /* "t" value at last msetup call */ 128 int msolve_type; /* interface type: 0=iterative; 1=direct; 2=custom */ 129 130 /* Counters */ 131 long int nfe; /* num fe calls */ 132 long int nfi; /* num fi calls */ 133 long int nsetups; /* num setup calls */ 134 135 /* Reusable arrays for fused vector operations */ 136 realtype *cvals; /* scalar array for fused ops */ 137 N_Vector *Xvecs; /* array of vectors for fused ops */ 138 int nfusedopvecs; /* length of cvals and Xvecs arrays */ 139 140 /* Data for using ARKStep with external polynomial forcing */ 141 booleantype expforcing; /* add forcing to explicit RHS */ 142 booleantype impforcing; /* add forcing to implicit RHS */ 143 realtype tshift; /* time normalization shift */ 144 realtype tscale; /* time normalization scaling */ 145 N_Vector* forcing; /* array of forcing vectors */ 146 int nforcing; /* number of forcing vectors */ 147 148 } *ARKodeARKStepMem; 149 150 151 /*=============================================================== 152 ARK time step module private function prototypes 153 ===============================================================*/ 154 155 /* Interface routines supplied to ARKode */ 156 int arkStep_AttachLinsol(void* arkode_mem, ARKLinsolInitFn linit, 157 ARKLinsolSetupFn lsetup, 158 ARKLinsolSolveFn lsolve, 159 ARKLinsolFreeFn lfree, 160 int lsolve_type, void *lmem); 161 int arkStep_AttachMasssol(void* arkode_mem, ARKMassInitFn minit, 162 ARKMassSetupFn msetup, 163 ARKMassMultFn mmult, 164 ARKMassSolveFn msolve, 165 ARKMassFreeFn lfree, 166 int msolve_type, void *mass_mem); 167 void arkStep_DisableLSetup(void* arkode_mem); 168 void arkStep_DisableMSetup(void* arkode_mem); 169 int arkStep_Init(void* arkode_mem, int init_type); 170 void* arkStep_GetLmem(void* arkode_mem); 171 void* arkStep_GetMassMem(void* arkode_mem); 172 ARKRhsFn arkStep_GetImplicitRHS(void* arkode_mem); 173 int arkStep_GetGammas(void* arkode_mem, realtype *gamma, 174 realtype *gamrat, booleantype **jcur, 175 booleantype *dgamma_fail); 176 int arkStep_FullRHS(void* arkode_mem, realtype t, 177 N_Vector y, N_Vector f, int mode); 178 int arkStep_TakeStep(void* arkode_mem, realtype *dsmPtr, int *nflagPtr); 179 180 /* Internal utility routines */ 181 int arkStep_AccessStepMem(void* arkode_mem, const char *fname, 182 ARKodeMem *ark_mem, ARKodeARKStepMem *step_mem); 183 booleantype arkStep_CheckNVector(N_Vector tmpl); 184 int arkStep_SetButcherTables(ARKodeMem ark_mem); 185 int arkStep_CheckButcherTables(ARKodeMem ark_mem); 186 int arkStep_Predict(ARKodeMem ark_mem, int istage, N_Vector yguess); 187 int arkStep_StageSetup(ARKodeMem ark_mem); 188 int arkStep_NlsInit(ARKodeMem ark_mem); 189 int arkStep_Nls(ARKodeMem ark_mem, int nflag); 190 int arkStep_ComputeSolutions(ARKodeMem ark_mem, realtype *dsm); 191 192 /* private functions passed to nonlinear solver */ 193 int arkStep_NlsResidual(N_Vector yy, N_Vector res, void* arkode_mem); 194 int arkStep_NlsFPFunction(N_Vector yy, N_Vector res, void* arkode_mem); 195 int arkStep_NlsLSetup(booleantype jbad, booleantype* jcur, void* arkode_mem); 196 int arkStep_NlsLSolve(N_Vector delta, void* arkode_mem); 197 int arkStep_NlsConvTest(SUNNonlinearSolver NLS, N_Vector y, N_Vector del, 198 realtype tol, N_Vector ewt, void* arkode_mem); 199 200 /* private functions used by MRIStep */ 201 int arkStep_SetInnerForcing(void* arkode_mem, realtype tshift, realtype tscale, 202 N_Vector *f, int nvecs); 203 204 /*=============================================================== 205 Reusable ARKStep Error Messages 206 ===============================================================*/ 207 208 /* Initialization and I/O error messages */ 209 #define MSG_ARKSTEP_NO_MEM "Time step module memory is NULL." 210 #define MSG_NLS_INIT_FAIL "The nonlinear solver's init routine failed." 211 212 /* Other error messages */ 213 #define MSG_ARK_MISSING_FE "Cannot specify that method is explicit without providing a function pointer to fe(t,y)." 214 #define MSG_ARK_MISSING_FI "Cannot specify that method is implicit without providing a function pointer to fi(t,y)." 215 #define MSG_ARK_MISSING_F "Cannot specify that method is ImEx without providing function pointers to fi(t,y) and fe(t,y)." 216 217 #ifdef __cplusplus 218 } 219 #endif 220 221 #endif 222