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 * This is the implementation file for the optional input and
15 * output functions for the ARKode ARKStep time stepper module.
16 *
17 * NOTE: many functions currently in arkode_io.c will move here,
18 * with slightly different names. The code transition will be
19 * minimal, but the documentation changes will be significant.
20 *--------------------------------------------------------------*/
21
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include "arkode_arkstep_impl.h"
26 #include <sundials/sundials_math.h>
27 #include <sundials/sundials_types.h>
28
29 #if defined(SUNDIALS_EXTENDED_PRECISION)
30 #define RSYM "Lg"
31 #else
32 #define RSYM "g"
33 #endif
34
35
36 /*===============================================================
37 ARKStep Optional input functions (wrappers for generic ARKode
38 utility routines). All are documented in arkode_io.c.
39 ===============================================================*/
ARKStepSetDenseOrder(void * arkode_mem,int dord)40 int ARKStepSetDenseOrder(void *arkode_mem, int dord) {
41 return(ARKStepSetInterpolantDegree(arkode_mem, dord)); }
ARKStepSetInterpolantDegree(void * arkode_mem,int degree)42 int ARKStepSetInterpolantDegree(void *arkode_mem, int degree) {
43 if (degree < 0) degree = ARK_INTERP_MAX_DEGREE;
44 return(arkSetInterpolantDegree(arkode_mem, degree)); }
ARKStepSetInterpolantType(void * arkode_mem,int itype)45 int ARKStepSetInterpolantType(void *arkode_mem, int itype) {
46 return(arkSetInterpolantType(arkode_mem, itype)); }
ARKStepSetErrHandlerFn(void * arkode_mem,ARKErrHandlerFn ehfun,void * eh_data)47 int ARKStepSetErrHandlerFn(void *arkode_mem, ARKErrHandlerFn ehfun,
48 void *eh_data) {
49 return(arkSetErrHandlerFn(arkode_mem, ehfun, eh_data)); }
ARKStepSetErrFile(void * arkode_mem,FILE * errfp)50 int ARKStepSetErrFile(void *arkode_mem, FILE *errfp) {
51 return(arkSetErrFile(arkode_mem, errfp)); }
ARKStepSetDiagnostics(void * arkode_mem,FILE * diagfp)52 int ARKStepSetDiagnostics(void *arkode_mem, FILE *diagfp) {
53 return(arkSetDiagnostics(arkode_mem, diagfp)); }
ARKStepSetMaxNumSteps(void * arkode_mem,long int mxsteps)54 int ARKStepSetMaxNumSteps(void *arkode_mem, long int mxsteps) {
55 return(arkSetMaxNumSteps(arkode_mem, mxsteps)); }
ARKStepSetMaxHnilWarns(void * arkode_mem,int mxhnil)56 int ARKStepSetMaxHnilWarns(void *arkode_mem, int mxhnil) {
57 return(arkSetMaxHnilWarns(arkode_mem, mxhnil)); }
ARKStepSetInitStep(void * arkode_mem,realtype hin)58 int ARKStepSetInitStep(void *arkode_mem, realtype hin) {
59 return(arkSetInitStep(arkode_mem, hin)); }
ARKStepSetMinStep(void * arkode_mem,realtype hmin)60 int ARKStepSetMinStep(void *arkode_mem, realtype hmin) {
61 return(arkSetMinStep(arkode_mem, hmin)); }
ARKStepSetMaxStep(void * arkode_mem,realtype hmax)62 int ARKStepSetMaxStep(void *arkode_mem, realtype hmax) {
63 return(arkSetMaxStep(arkode_mem, hmax)); }
ARKStepSetStopTime(void * arkode_mem,realtype tstop)64 int ARKStepSetStopTime(void *arkode_mem, realtype tstop) {
65 return(arkSetStopTime(arkode_mem, tstop)); }
ARKStepSetRootDirection(void * arkode_mem,int * rootdir)66 int ARKStepSetRootDirection(void *arkode_mem, int *rootdir) {
67 return(arkSetRootDirection(arkode_mem, rootdir)); }
ARKStepSetNoInactiveRootWarn(void * arkode_mem)68 int ARKStepSetNoInactiveRootWarn(void *arkode_mem) {
69 return(arkSetNoInactiveRootWarn(arkode_mem)); }
ARKStepSetConstraints(void * arkode_mem,N_Vector constraints)70 int ARKStepSetConstraints(void *arkode_mem, N_Vector constraints) {
71 return(arkSetConstraints(arkode_mem, constraints)); }
ARKStepSetMaxNumConstrFails(void * arkode_mem,int maxfails)72 int ARKStepSetMaxNumConstrFails(void *arkode_mem, int maxfails) {
73 return(arkSetMaxNumConstrFails(arkode_mem, maxfails)); }
ARKStepSetPostprocessStepFn(void * arkode_mem,ARKPostProcessFn ProcessStep)74 int ARKStepSetPostprocessStepFn(void *arkode_mem,
75 ARKPostProcessFn ProcessStep) {
76 return(arkSetPostprocessStepFn(arkode_mem, ProcessStep)); }
ARKStepSetPostprocessStageFn(void * arkode_mem,ARKPostProcessFn ProcessStage)77 int ARKStepSetPostprocessStageFn(void *arkode_mem,
78 ARKPostProcessFn ProcessStage) {
79 return(arkSetPostprocessStageFn(arkode_mem, ProcessStage)); }
ARKStepSetCFLFraction(void * arkode_mem,realtype cfl_frac)80 int ARKStepSetCFLFraction(void *arkode_mem, realtype cfl_frac) {
81 return(arkSetCFLFraction(arkode_mem, cfl_frac)); }
ARKStepSetSafetyFactor(void * arkode_mem,realtype safety)82 int ARKStepSetSafetyFactor(void *arkode_mem, realtype safety) {
83 return(arkSetSafetyFactor(arkode_mem, safety)); }
ARKStepSetErrorBias(void * arkode_mem,realtype bias)84 int ARKStepSetErrorBias(void *arkode_mem, realtype bias) {
85 return(arkSetErrorBias(arkode_mem, bias)); }
ARKStepSetMaxGrowth(void * arkode_mem,realtype mx_growth)86 int ARKStepSetMaxGrowth(void *arkode_mem, realtype mx_growth) {
87 return(arkSetMaxGrowth(arkode_mem, mx_growth)); }
ARKStepSetMinReduction(void * arkode_mem,realtype eta_min)88 int ARKStepSetMinReduction(void *arkode_mem, realtype eta_min) {
89 return(arkSetMinReduction(arkode_mem, eta_min)); }
ARKStepSetFixedStepBounds(void * arkode_mem,realtype lb,realtype ub)90 int ARKStepSetFixedStepBounds(void *arkode_mem, realtype lb, realtype ub) {
91 return(arkSetFixedStepBounds(arkode_mem, lb, ub)); }
ARKStepSetAdaptivityMethod(void * arkode_mem,int imethod,int idefault,int pq,realtype adapt_params[3])92 int ARKStepSetAdaptivityMethod(void *arkode_mem, int imethod, int idefault,
93 int pq, realtype adapt_params[3]) {
94 return(arkSetAdaptivityMethod(arkode_mem, imethod, idefault, pq, adapt_params)); }
ARKStepSetAdaptivityFn(void * arkode_mem,ARKAdaptFn hfun,void * h_data)95 int ARKStepSetAdaptivityFn(void *arkode_mem, ARKAdaptFn hfun, void *h_data) {
96 return(arkSetAdaptivityFn(arkode_mem, hfun, h_data)); }
ARKStepSetMaxFirstGrowth(void * arkode_mem,realtype etamx1)97 int ARKStepSetMaxFirstGrowth(void *arkode_mem, realtype etamx1) {
98 return(arkSetMaxFirstGrowth(arkode_mem, etamx1)); }
ARKStepSetMaxEFailGrowth(void * arkode_mem,realtype etamxf)99 int ARKStepSetMaxEFailGrowth(void *arkode_mem, realtype etamxf) {
100 return(arkSetMaxEFailGrowth(arkode_mem, etamxf)); }
ARKStepSetSmallNumEFails(void * arkode_mem,int small_nef)101 int ARKStepSetSmallNumEFails(void *arkode_mem, int small_nef) {
102 return(arkSetSmallNumEFails(arkode_mem, small_nef)); }
ARKStepSetMaxCFailGrowth(void * arkode_mem,realtype etacf)103 int ARKStepSetMaxCFailGrowth(void *arkode_mem, realtype etacf) {
104 return(arkSetMaxCFailGrowth(arkode_mem, etacf)); }
ARKStepSetStabilityFn(void * arkode_mem,ARKExpStabFn EStab,void * estab_data)105 int ARKStepSetStabilityFn(void *arkode_mem, ARKExpStabFn EStab, void *estab_data) {
106 return(arkSetStabilityFn(arkode_mem, EStab, estab_data)); }
ARKStepSetMaxErrTestFails(void * arkode_mem,int maxnef)107 int ARKStepSetMaxErrTestFails(void *arkode_mem, int maxnef) {
108 return(arkSetMaxErrTestFails(arkode_mem, maxnef)); }
ARKStepSetMaxConvFails(void * arkode_mem,int maxncf)109 int ARKStepSetMaxConvFails(void *arkode_mem, int maxncf) {
110 return(arkSetMaxConvFails(arkode_mem, maxncf)); }
ARKStepSetFixedStep(void * arkode_mem,realtype hfixed)111 int ARKStepSetFixedStep(void *arkode_mem, realtype hfixed) {
112 return(arkSetFixedStep(arkode_mem, hfixed)); }
113
114
115 /*===============================================================
116 ARKStep Optional input functions (customized wrappers for
117 generic ARKode utility routines). All are documented in
118 arkode_io.c and arkode_ls.c.
119 ===============================================================*/
120
ARKStepSetUserData(void * arkode_mem,void * user_data)121 int ARKStepSetUserData(void *arkode_mem, void *user_data)
122 {
123 ARKodeMem ark_mem;
124 ARKodeARKStepMem step_mem;
125 int retval;
126
127 /* access ARKodeARKStepMem structure */
128 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetUserData",
129 &ark_mem, &step_mem);
130 if (retval != ARK_SUCCESS) return(retval);
131
132 /* set user_data in ARKode mem */
133 retval = arkSetUserData(arkode_mem, user_data);
134 if (retval != ARK_SUCCESS) return(retval);
135
136 /* set user data in ARKodeLS mem */
137 if (step_mem->lmem != NULL) {
138 retval = arkLSSetUserData(arkode_mem, user_data);
139 if (retval != ARKLS_SUCCESS) return(retval);
140 }
141
142 /* set user data in ARKodeLSMass mem */
143 if (step_mem->mass_mem != NULL) {
144 retval = arkLSSetMassUserData(arkode_mem, user_data);
145 if (retval != ARKLS_SUCCESS) return(retval);
146 }
147
148 return(ARK_SUCCESS);
149 }
150
151
152 /*---------------------------------------------------------------
153 ARKStepSetStagePredictFn: Specifies a user-provided step
154 predictor function having type ARKStepStagePredictFn. A
155 NULL input function disables calls to this routine.
156 ---------------------------------------------------------------*/
ARKStepSetStagePredictFn(void * arkode_mem,ARKStepStagePredictFn PredictStage)157 int ARKStepSetStagePredictFn(void *arkode_mem,
158 ARKStepStagePredictFn PredictStage)
159 {
160 ARKodeMem ark_mem;
161 ARKodeARKStepMem step_mem;
162 int retval;
163
164 /* access ARKodeARKStepMem structure and set function pointer */
165 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetStagePredictFn",
166 &ark_mem, &step_mem);
167 if (retval != ARK_SUCCESS) return(retval);
168
169 /* override predictor method 5 if non-NULL PredictStage is supplied */
170 if ((step_mem->predictor == 5) && (PredictStage != NULL)) {
171 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
172 "ARKStepSetStagePredictFn",
173 "User-supplied predictor is incompatible with predictor method 5");
174 return(ARK_ILL_INPUT);
175 }
176
177 step_mem->stage_predict = PredictStage;
178 return(ARK_SUCCESS);
179 }
180
181
182 /*---------------------------------------------------------------
183 These wrappers for ARKLs module 'set' routines all are
184 documented in arkode_arkstep.h.
185 ---------------------------------------------------------------*/
ARKStepSetLinearSolver(void * arkode_mem,SUNLinearSolver LS,SUNMatrix A)186 int ARKStepSetLinearSolver(void *arkode_mem, SUNLinearSolver LS,
187 SUNMatrix A) {
188 return(arkLSSetLinearSolver(arkode_mem, LS, A)); }
ARKStepSetMassLinearSolver(void * arkode_mem,SUNLinearSolver LS,SUNMatrix M,booleantype time_dep)189 int ARKStepSetMassLinearSolver(void *arkode_mem, SUNLinearSolver LS,
190 SUNMatrix M, booleantype time_dep) {
191 return(arkLSSetMassLinearSolver(arkode_mem, LS, M, time_dep)); }
ARKStepSetJacFn(void * arkode_mem,ARKLsJacFn jac)192 int ARKStepSetJacFn(void *arkode_mem, ARKLsJacFn jac) {
193 return(arkLSSetJacFn(arkode_mem, jac)); }
ARKStepSetMassFn(void * arkode_mem,ARKLsMassFn mass)194 int ARKStepSetMassFn(void *arkode_mem, ARKLsMassFn mass) {
195 return(arkLSSetMassFn(arkode_mem, mass)); }
ARKStepSetMaxStepsBetweenJac(void * arkode_mem,long int msbj)196 int ARKStepSetMaxStepsBetweenJac(void *arkode_mem, long int msbj) {
197 return(arkLSSetMaxStepsBetweenJac(arkode_mem, msbj)); }
ARKStepSetLinearSolutionScaling(void * arkode_mem,booleantype onoff)198 int ARKStepSetLinearSolutionScaling(void *arkode_mem, booleantype onoff) {
199 return(arkLSSetLinearSolutionScaling(arkode_mem, onoff)); }
ARKStepSetEpsLin(void * arkode_mem,realtype eplifac)200 int ARKStepSetEpsLin(void *arkode_mem, realtype eplifac) {
201 return(arkLSSetEpsLin(arkode_mem, eplifac)); }
ARKStepSetMassEpsLin(void * arkode_mem,realtype eplifac)202 int ARKStepSetMassEpsLin(void *arkode_mem, realtype eplifac) {
203 return(arkLSSetMassEpsLin(arkode_mem, eplifac)); }
ARKStepSetPreconditioner(void * arkode_mem,ARKLsPrecSetupFn psetup,ARKLsPrecSolveFn psolve)204 int ARKStepSetPreconditioner(void *arkode_mem, ARKLsPrecSetupFn psetup,
205 ARKLsPrecSolveFn psolve) {
206 return(arkLSSetPreconditioner(arkode_mem, psetup, psolve)); }
ARKStepSetMassPreconditioner(void * arkode_mem,ARKLsMassPrecSetupFn psetup,ARKLsMassPrecSolveFn psolve)207 int ARKStepSetMassPreconditioner(void *arkode_mem, ARKLsMassPrecSetupFn psetup,
208 ARKLsMassPrecSolveFn psolve) {
209 return(arkLSSetMassPreconditioner(arkode_mem, psetup, psolve)); }
ARKStepSetJacTimes(void * arkode_mem,ARKLsJacTimesSetupFn jtsetup,ARKLsJacTimesVecFn jtimes)210 int ARKStepSetJacTimes(void *arkode_mem, ARKLsJacTimesSetupFn jtsetup,
211 ARKLsJacTimesVecFn jtimes) {
212 return(arkLSSetJacTimes(arkode_mem, jtsetup, jtimes)); }
ARKStepSetJacTimesRhsFn(void * arkode_mem,ARKRhsFn jtimesRhsFn)213 int ARKStepSetJacTimesRhsFn(void *arkode_mem, ARKRhsFn jtimesRhsFn) {
214 return(arkLSSetJacTimesRhsFn(arkode_mem, jtimesRhsFn)); }
ARKStepSetMassTimes(void * arkode_mem,ARKLsMassTimesSetupFn msetup,ARKLsMassTimesVecFn mtimes,void * mtimes_data)215 int ARKStepSetMassTimes(void *arkode_mem, ARKLsMassTimesSetupFn msetup,
216 ARKLsMassTimesVecFn mtimes, void *mtimes_data) {
217 return(arkLSSetMassTimes(arkode_mem, msetup, mtimes, mtimes_data)); }
ARKStepSetLinSysFn(void * arkode_mem,ARKLsLinSysFn linsys)218 int ARKStepSetLinSysFn(void *arkode_mem, ARKLsLinSysFn linsys) {
219 return(arkLSSetLinSysFn(arkode_mem, linsys)); }
220
221
222 /*===============================================================
223 ARKStep Optional output functions (wrappers for generic ARKode
224 utility routines). All are documented in arkode_io.c.
225 ===============================================================*/
ARKStepGetNumStepAttempts(void * arkode_mem,long int * nstep_attempts)226 int ARKStepGetNumStepAttempts(void *arkode_mem, long int *nstep_attempts) {
227 return(arkGetNumStepAttempts(arkode_mem, nstep_attempts)); }
ARKStepGetNumSteps(void * arkode_mem,long int * nsteps)228 int ARKStepGetNumSteps(void *arkode_mem, long int *nsteps) {
229 return(arkGetNumSteps(arkode_mem, nsteps)); }
ARKStepGetActualInitStep(void * arkode_mem,realtype * hinused)230 int ARKStepGetActualInitStep(void *arkode_mem, realtype *hinused) {
231 return(arkGetActualInitStep(arkode_mem, hinused)); }
ARKStepGetLastStep(void * arkode_mem,realtype * hlast)232 int ARKStepGetLastStep(void *arkode_mem, realtype *hlast) {
233 return(arkGetLastStep(arkode_mem, hlast)); }
ARKStepGetCurrentStep(void * arkode_mem,realtype * hcur)234 int ARKStepGetCurrentStep(void *arkode_mem, realtype *hcur) {
235 return(arkGetCurrentStep(arkode_mem, hcur)); }
ARKStepGetCurrentTime(void * arkode_mem,realtype * tcur)236 int ARKStepGetCurrentTime(void *arkode_mem, realtype *tcur) {
237 return(arkGetCurrentTime(arkode_mem, tcur)); }
ARKStepGetCurrentState(void * arkode_mem,N_Vector * ycur)238 int ARKStepGetCurrentState(void *arkode_mem, N_Vector *ycur) {
239 return(arkGetCurrentState(arkode_mem, ycur)); }
ARKStepGetTolScaleFactor(void * arkode_mem,realtype * tolsfact)240 int ARKStepGetTolScaleFactor(void *arkode_mem, realtype *tolsfact) {
241 return(arkGetTolScaleFactor(arkode_mem, tolsfact)); }
ARKStepGetErrWeights(void * arkode_mem,N_Vector eweight)242 int ARKStepGetErrWeights(void *arkode_mem, N_Vector eweight) {
243 return(arkGetErrWeights(arkode_mem, eweight)); }
ARKStepGetResWeights(void * arkode_mem,N_Vector rweight)244 int ARKStepGetResWeights(void *arkode_mem, N_Vector rweight) {
245 return(arkGetResWeights(arkode_mem, rweight)); }
ARKStepGetWorkSpace(void * arkode_mem,long int * lenrw,long int * leniw)246 int ARKStepGetWorkSpace(void *arkode_mem, long int *lenrw, long int *leniw) {
247 return(arkGetWorkSpace(arkode_mem, lenrw, leniw)); }
ARKStepGetNumGEvals(void * arkode_mem,long int * ngevals)248 int ARKStepGetNumGEvals(void *arkode_mem, long int *ngevals) {
249 return(arkGetNumGEvals(arkode_mem, ngevals)); }
ARKStepGetRootInfo(void * arkode_mem,int * rootsfound)250 int ARKStepGetRootInfo(void *arkode_mem, int *rootsfound) {
251 return(arkGetRootInfo(arkode_mem, rootsfound)); }
ARKStepGetStepStats(void * arkode_mem,long int * nsteps,realtype * hinused,realtype * hlast,realtype * hcur,realtype * tcur)252 int ARKStepGetStepStats(void *arkode_mem, long int *nsteps,
253 realtype *hinused, realtype *hlast,
254 realtype *hcur, realtype *tcur) {
255 return(arkGetStepStats(arkode_mem, nsteps, hinused, hlast, hcur, tcur)); }
ARKStepGetNumConstrFails(void * arkode_mem,long int * nconstrfails)256 int ARKStepGetNumConstrFails(void *arkode_mem, long int *nconstrfails) {
257 return(arkGetNumConstrFails(arkode_mem, nconstrfails)); }
ARKStepGetNumExpSteps(void * arkode_mem,long int * nsteps)258 int ARKStepGetNumExpSteps(void *arkode_mem, long int *nsteps) {
259 return(arkGetNumExpSteps(arkode_mem, nsteps)); }
ARKStepGetNumAccSteps(void * arkode_mem,long int * nsteps)260 int ARKStepGetNumAccSteps(void *arkode_mem, long int *nsteps) {
261 return(arkGetNumAccSteps(arkode_mem, nsteps)); }
ARKStepGetNumErrTestFails(void * arkode_mem,long int * netfails)262 int ARKStepGetNumErrTestFails(void *arkode_mem, long int *netfails) {
263 return(arkGetNumErrTestFails(arkode_mem, netfails)); }
ARKStepGetReturnFlagName(long int flag)264 char *ARKStepGetReturnFlagName(long int flag) {
265 return(arkGetReturnFlagName(flag)); }
266
267 /*---------------------------------------------------------------
268 These wrappers for ARKLs module 'get' routines all are
269 documented in arkode_arkstep.h.
270 ---------------------------------------------------------------*/
ARKStepGetLinWorkSpace(void * arkode_mem,long int * lenrwLS,long int * leniwLS)271 int ARKStepGetLinWorkSpace(void *arkode_mem, long int *lenrwLS, long int *leniwLS) {
272 return(arkLSGetWorkSpace(arkode_mem, lenrwLS, leniwLS)); }
ARKStepGetNumJacEvals(void * arkode_mem,long int * njevals)273 int ARKStepGetNumJacEvals(void *arkode_mem, long int *njevals) {
274 return(arkLSGetNumJacEvals(arkode_mem, njevals)); }
ARKStepGetNumPrecEvals(void * arkode_mem,long int * npevals)275 int ARKStepGetNumPrecEvals(void *arkode_mem, long int *npevals) {
276 return(arkLSGetNumPrecEvals(arkode_mem, npevals)); }
ARKStepGetNumPrecSolves(void * arkode_mem,long int * npsolves)277 int ARKStepGetNumPrecSolves(void *arkode_mem, long int *npsolves) {
278 return(arkLSGetNumPrecSolves(arkode_mem, npsolves)); }
ARKStepGetNumLinIters(void * arkode_mem,long int * nliters)279 int ARKStepGetNumLinIters(void *arkode_mem, long int *nliters) {
280 return(arkLSGetNumLinIters(arkode_mem, nliters)); }
ARKStepGetNumLinConvFails(void * arkode_mem,long int * nlcfails)281 int ARKStepGetNumLinConvFails(void *arkode_mem, long int *nlcfails) {
282 return(arkLSGetNumConvFails(arkode_mem, nlcfails)); }
ARKStepGetNumJTSetupEvals(void * arkode_mem,long int * njtsetups)283 int ARKStepGetNumJTSetupEvals(void *arkode_mem, long int *njtsetups) {
284 return(arkLSGetNumJTSetupEvals(arkode_mem, njtsetups)); }
ARKStepGetNumJtimesEvals(void * arkode_mem,long int * njvevals)285 int ARKStepGetNumJtimesEvals(void *arkode_mem, long int *njvevals) {
286 return(arkLSGetNumJtimesEvals(arkode_mem, njvevals)); }
ARKStepGetNumLinRhsEvals(void * arkode_mem,long int * nfevalsLS)287 int ARKStepGetNumLinRhsEvals(void *arkode_mem, long int *nfevalsLS) {
288 return(arkLSGetNumRhsEvals(arkode_mem, nfevalsLS)); }
ARKStepGetLastLinFlag(void * arkode_mem,long int * flag)289 int ARKStepGetLastLinFlag(void *arkode_mem, long int *flag) {
290 return(arkLSGetLastFlag(arkode_mem, flag)); }
291
ARKStepGetMassWorkSpace(void * arkode_mem,long int * lenrwMLS,long int * leniwMLS)292 int ARKStepGetMassWorkSpace(void *arkode_mem, long int *lenrwMLS, long int *leniwMLS) {
293 return(arkLSGetMassWorkSpace(arkode_mem, lenrwMLS, leniwMLS)); }
ARKStepGetNumMassSetups(void * arkode_mem,long int * nmsetups)294 int ARKStepGetNumMassSetups(void *arkode_mem, long int *nmsetups) {
295 return(arkLSGetNumMassSetups(arkode_mem, nmsetups)); }
ARKStepGetNumMassMultSetups(void * arkode_mem,long int * nmvsetups)296 int ARKStepGetNumMassMultSetups(void *arkode_mem, long int *nmvsetups) {
297 return(arkLSGetNumMassMatvecSetups(arkode_mem, nmvsetups)); }
ARKStepGetNumMassMult(void * arkode_mem,long int * nmvevals)298 int ARKStepGetNumMassMult(void *arkode_mem, long int *nmvevals) {
299 return(arkLSGetNumMassMult(arkode_mem, nmvevals)); }
ARKStepGetNumMassSolves(void * arkode_mem,long int * nmsolves)300 int ARKStepGetNumMassSolves(void *arkode_mem, long int *nmsolves) {
301 return(arkLSGetNumMassSolves(arkode_mem, nmsolves)); }
ARKStepGetNumMassPrecEvals(void * arkode_mem,long int * nmpevals)302 int ARKStepGetNumMassPrecEvals(void *arkode_mem, long int *nmpevals) {
303 return(arkLSGetNumMassPrecEvals(arkode_mem, nmpevals)); }
ARKStepGetNumMassPrecSolves(void * arkode_mem,long int * nmpsolves)304 int ARKStepGetNumMassPrecSolves(void *arkode_mem, long int *nmpsolves) {
305 return(arkLSGetNumMassPrecSolves(arkode_mem, nmpsolves)); }
ARKStepGetNumMassIters(void * arkode_mem,long int * nmiters)306 int ARKStepGetNumMassIters(void *arkode_mem, long int *nmiters) {
307 return(arkLSGetNumMassIters(arkode_mem, nmiters)); }
ARKStepGetNumMassConvFails(void * arkode_mem,long int * nmcfails)308 int ARKStepGetNumMassConvFails(void *arkode_mem, long int *nmcfails) {
309 return(arkLSGetNumMassConvFails(arkode_mem, nmcfails)); }
ARKStepGetNumMTSetups(void * arkode_mem,long int * nmtsetups)310 int ARKStepGetNumMTSetups(void *arkode_mem, long int *nmtsetups) {
311 return(arkLSGetNumMTSetups(arkode_mem, nmtsetups)); }
ARKStepGetLastMassFlag(void * arkode_mem,long int * flag)312 int ARKStepGetLastMassFlag(void *arkode_mem, long int *flag) {
313 return(arkLSGetLastMassFlag(arkode_mem, flag)); }
314
ARKStepGetLinReturnFlagName(long int flag)315 char *ARKStepGetLinReturnFlagName(long int flag) {
316 return(arkLSGetReturnFlagName(flag)); }
317
318
319
320 /*===============================================================
321 ARKStep optional input functions -- stepper-specific
322 ===============================================================*/
323
324 /*---------------------------------------------------------------
325 ARKStepSetDefaults:
326
327 Resets all ARKStep optional inputs to their default values.
328 Does not change problem-defining function pointers or
329 user_data pointer. Also leaves alone any data
330 structures/options related to the ARKode infrastructure itself
331 (e.g., root-finding and post-process step).
332 ---------------------------------------------------------------*/
ARKStepSetDefaults(void * arkode_mem)333 int ARKStepSetDefaults(void* arkode_mem)
334 {
335 ARKodeMem ark_mem;
336 ARKodeARKStepMem step_mem;
337 int retval;
338
339 /* access ARKodeARKStepMem structure */
340 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetDefaults",
341 &ark_mem, &step_mem);
342 if (retval != ARK_SUCCESS) return(retval);
343
344 /* Set default ARKode infrastructure parameters */
345 retval = arkSetDefaults(ark_mem);
346 if (retval != ARK_SUCCESS) {
347 arkProcessError(NULL, ARK_MEM_NULL, "ARKode::ARKStep",
348 "ARKStepSetDefaults",
349 "Error setting ARKode infrastructure defaults");
350 return(retval);
351 }
352
353 /* Set default values for integrator optional inputs */
354 step_mem->q = Q_DEFAULT; /* method order */
355 step_mem->p = 0; /* embedding order */
356 step_mem->predictor = 0; /* trivial predictor */
357 step_mem->linear = SUNFALSE; /* nonlinear problem */
358 step_mem->linear_timedep = SUNTRUE; /* dfi/dy depends on t */
359 step_mem->explicit = SUNTRUE; /* fe(t,y) will be used */
360 step_mem->implicit = SUNTRUE; /* fi(t,y) will be used */
361 step_mem->maxcor = MAXCOR; /* max nonlinear iters/stage */
362 step_mem->nlscoef = NLSCOEF; /* nonlinear tolerance coefficient */
363 step_mem->crdown = CRDOWN; /* nonlinear convergence estimate coeff. */
364 step_mem->rdiv = RDIV; /* nonlinear divergence tolerance */
365 step_mem->dgmax = DGMAX; /* max step change before recomputing J or P */
366 step_mem->msbp = MSBP; /* max steps between updates to J or P */
367 step_mem->stages = 0; /* no stages */
368 step_mem->istage = 0; /* current stage */
369 step_mem->Be = NULL; /* no Butcher tables */
370 step_mem->Bi = NULL;
371 step_mem->NLS = NULL; /* no nonlinear solver object */
372 step_mem->jcur = SUNFALSE;
373 step_mem->convfail = ARK_NO_FAILURES;
374 step_mem->stage_predict = NULL; /* no user-supplied stage predictor */
375 return(ARK_SUCCESS);
376 }
377
378
379 /*---------------------------------------------------------------
380 ARKStepSetOptimalParams:
381
382 Sets all adaptivity and solver parameters to our 'best guess'
383 values, for a given ARKStep integration method (ERK, DIRK, ARK),
384 a given method order, and a given nonlinear solver type. Should
385 only be called after the method order, solver, and integration
386 method have been set, and only if time step adaptivity is
387 enabled.
388 ---------------------------------------------------------------*/
ARKStepSetOptimalParams(void * arkode_mem)389 int ARKStepSetOptimalParams(void *arkode_mem)
390 {
391 ARKodeMem ark_mem;
392 ARKodeARKStepMem step_mem;
393 ARKodeHAdaptMem hadapt_mem;
394 int retval;
395
396 /* access ARKodeARKStepMem structure */
397 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetOptimalParams",
398 &ark_mem, &step_mem);
399 if (retval != ARK_SUCCESS) return(retval);
400
401 /* access ARKodeHAdaptMem structure */
402 if (ark_mem->hadapt_mem == NULL) {
403 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
404 "ARKStepSetOptimalParams",
405 MSG_ARKADAPT_NO_MEM);
406 return(ARK_MEM_NULL);
407 }
408 hadapt_mem = ark_mem->hadapt_mem;
409
410 /* Choose values based on method, order */
411
412 /* explicit */
413 if (step_mem->explicit && !step_mem->implicit) {
414 hadapt_mem->imethod = 1;
415 hadapt_mem->safety = RCONST(0.99);
416 hadapt_mem->bias = RCONST(1.2);
417 hadapt_mem->growth = RCONST(25.0);
418 hadapt_mem->k1 = RCONST(0.8);
419 hadapt_mem->k2 = RCONST(0.31);
420 hadapt_mem->etamxf = RCONST(0.3);
421
422 /* implicit */
423 } else if (step_mem->implicit && !step_mem->explicit) {
424 switch (step_mem->q) {
425 case 2: /* just use standard defaults since better ones unknown */
426 hadapt_mem->imethod = 0;
427 hadapt_mem->safety = SAFETY;
428 hadapt_mem->bias = BIAS;
429 hadapt_mem->growth = GROWTH;
430 hadapt_mem->etamxf = ETAMXF;
431 hadapt_mem->small_nef = SMALL_NEF;
432 hadapt_mem->etacf = ETACF;
433 step_mem->nlscoef = RCONST(0.001);
434 step_mem->maxcor = 5;
435 step_mem->crdown = CRDOWN;
436 step_mem->rdiv = RDIV;
437 step_mem->dgmax = DGMAX;
438 step_mem->msbp = MSBP;
439 break;
440 case 3:
441 hadapt_mem->imethod = 2;
442 hadapt_mem->safety = RCONST(0.957);
443 hadapt_mem->bias = RCONST(1.9);
444 hadapt_mem->growth = RCONST(17.6);
445 hadapt_mem->etamxf = RCONST(0.45);
446 hadapt_mem->small_nef = SMALL_NEF;
447 hadapt_mem->etacf = ETACF;
448 step_mem->nlscoef = RCONST(0.22);
449 step_mem->crdown = RCONST(0.17);
450 step_mem->rdiv = RCONST(2.3);
451 step_mem->dgmax = RCONST(0.19);
452 step_mem->msbp = 60;
453 break;
454 case 4:
455 hadapt_mem->imethod = 0;
456 hadapt_mem->safety = RCONST(0.988);
457 hadapt_mem->bias = RCONST(1.2);
458 hadapt_mem->growth = RCONST(31.5);
459 hadapt_mem->k1 = RCONST(0.535);
460 hadapt_mem->k2 = RCONST(0.209);
461 hadapt_mem->k3 = RCONST(0.148);
462 hadapt_mem->etamxf = RCONST(0.33);
463 hadapt_mem->small_nef = SMALL_NEF;
464 hadapt_mem->etacf = ETACF;
465 step_mem->nlscoef = RCONST(0.24);
466 step_mem->crdown = RCONST(0.26);
467 step_mem->rdiv = RCONST(2.3);
468 step_mem->dgmax = RCONST(0.16);
469 step_mem->msbp = 31;
470 break;
471 case 5:
472 hadapt_mem->imethod = 0;
473 hadapt_mem->safety = RCONST(0.937);
474 hadapt_mem->bias = RCONST(3.3);
475 hadapt_mem->growth = RCONST(22.0);
476 hadapt_mem->k1 = RCONST(0.56);
477 hadapt_mem->k2 = RCONST(0.338);
478 hadapt_mem->k3 = RCONST(0.14);
479 hadapt_mem->etamxf = RCONST(0.44);
480 hadapt_mem->small_nef = SMALL_NEF;
481 hadapt_mem->etacf = ETACF;
482 step_mem->nlscoef = RCONST(0.25);
483 step_mem->crdown = RCONST(0.4);
484 step_mem->rdiv = RCONST(2.3);
485 step_mem->dgmax = RCONST(0.32);
486 step_mem->msbp = 31;
487 break;
488 }
489
490 /* imex */
491 } else {
492 switch (step_mem->q) {
493 case 3:
494 hadapt_mem->imethod = 0;
495 hadapt_mem->safety = RCONST(0.965);
496 hadapt_mem->bias = RCONST(1.42);
497 hadapt_mem->growth = RCONST(28.7);
498 hadapt_mem->k1 = RCONST(0.54);
499 hadapt_mem->k2 = RCONST(0.36);
500 hadapt_mem->k3 = RCONST(0.14);
501 hadapt_mem->etamxf = RCONST(0.46);
502 hadapt_mem->small_nef = SMALL_NEF;
503 hadapt_mem->etacf = ETACF;
504 step_mem->nlscoef = RCONST(0.22);
505 step_mem->crdown = RCONST(0.17);
506 step_mem->rdiv = RCONST(2.3);
507 step_mem->dgmax = RCONST(0.19);
508 step_mem->msbp = 60;
509 break;
510 case 4:
511 hadapt_mem->imethod = 0;
512 hadapt_mem->safety = RCONST(0.97);
513 hadapt_mem->bias = RCONST(1.35);
514 hadapt_mem->growth = RCONST(25.0);
515 hadapt_mem->k1 = RCONST(0.543);
516 hadapt_mem->k2 = RCONST(0.297);
517 hadapt_mem->k3 = RCONST(0.14);
518 hadapt_mem->etamxf = RCONST(0.47);
519 hadapt_mem->small_nef = SMALL_NEF;
520 hadapt_mem->etacf = ETACF;
521 step_mem->nlscoef = RCONST(0.24);
522 step_mem->crdown = RCONST(0.26);
523 step_mem->rdiv = RCONST(2.3);
524 step_mem->dgmax = RCONST(0.16);
525 step_mem->msbp = 31;
526 break;
527 case 5:
528 hadapt_mem->imethod = 1;
529 hadapt_mem->safety = RCONST(0.993);
530 hadapt_mem->bias = RCONST(1.15);
531 hadapt_mem->growth = RCONST(28.5);
532 hadapt_mem->k1 = RCONST(0.8);
533 hadapt_mem->k2 = RCONST(0.35);
534 hadapt_mem->etamxf = RCONST(0.3);
535 hadapt_mem->small_nef = SMALL_NEF;
536 hadapt_mem->etacf = ETACF;
537 step_mem->nlscoef = RCONST(0.25);
538 step_mem->crdown = RCONST(0.4);
539 step_mem->rdiv = RCONST(2.3);
540 step_mem->dgmax = RCONST(0.32);
541 step_mem->msbp = 31;
542 break;
543 }
544
545 }
546 return(ARK_SUCCESS);
547 }
548
549
550 /*---------------------------------------------------------------
551 ARKStepSetOrder:
552
553 Specifies the method order
554
555 ** Note in documentation that this should not be called along
556 with ARKStepSetTable or ARKStepSetTableNum. This routine
557 is used to specify a desired method order using default Butcher
558 tables, whereas any user-supplied table will have their own
559 order associated with them.
560 ---------------------------------------------------------------*/
ARKStepSetOrder(void * arkode_mem,int ord)561 int ARKStepSetOrder(void *arkode_mem, int ord)
562 {
563 ARKodeMem ark_mem;
564 ARKodeARKStepMem step_mem;
565 int retval;
566
567 /* access ARKodeARKStepMem structure */
568 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetOrder",
569 &ark_mem, &step_mem);
570 if (retval != ARK_SUCCESS) return(retval);
571
572 /* set user-provided value, or default, depending on argument */
573 if (ord <= 0) {
574 step_mem->q = Q_DEFAULT;
575 } else {
576 step_mem->q = ord;
577 }
578
579 /* clear Butcher tables, since user is requesting a change in method
580 or a reset to defaults. Tables will be set in ARKInitialSetup. */
581 step_mem->stages = 0;
582 step_mem->istage = 0;
583 step_mem->p = 0;
584 ARKodeButcherTable_Free(step_mem->Be); step_mem->Be = NULL;
585 ARKodeButcherTable_Free(step_mem->Bi); step_mem->Bi = NULL;
586
587 return(ARK_SUCCESS);
588 }
589
590
591 /*---------------------------------------------------------------
592 ARKStepSetLinear:
593
594 Specifies that the implicit portion of the problem is linear,
595 and to tighten the linear solver tolerances while taking only
596 one Newton iteration. DO NOT USE IN COMBINATION WITH THE
597 FIXED-POINT SOLVER. Automatically tightens DeltaGammaMax
598 to ensure that step size changes cause Jacobian recomputation.
599
600 The argument should be 1 or 0, where 1 indicates that the
601 Jacobian of fi with respect to y depends on time, and
602 0 indicates that it is not time dependent. Alternately, when
603 using an iterative linear solver this flag denotes time
604 dependence of the preconditioner.
605 ---------------------------------------------------------------*/
ARKStepSetLinear(void * arkode_mem,int timedepend)606 int ARKStepSetLinear(void *arkode_mem, int timedepend)
607 {
608 ARKodeMem ark_mem;
609 ARKodeARKStepMem step_mem;
610 int retval;
611
612 /* access ARKodeARKStepMem structure */
613 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetLinear",
614 &ark_mem, &step_mem);
615 if (retval != ARK_SUCCESS) return(retval);
616
617 /* set parameters */
618 step_mem->linear = SUNTRUE;
619 step_mem->linear_timedep = (timedepend == 1);
620 step_mem->dgmax = RCONST(100.0)*UNIT_ROUNDOFF;
621
622 return(ARK_SUCCESS);
623 }
624
625
626 /*---------------------------------------------------------------
627 ARKStepSetNonlinear:
628
629 Specifies that the implicit portion of the problem is nonlinear.
630 Used to undo a previous call to ARKStepSetLinear. Automatically
631 loosens DeltaGammaMax back to default value.
632 ---------------------------------------------------------------*/
ARKStepSetNonlinear(void * arkode_mem)633 int ARKStepSetNonlinear(void *arkode_mem)
634 {
635 ARKodeMem ark_mem;
636 ARKodeARKStepMem step_mem;
637 int retval;
638
639 /* access ARKodeARKStepMem structure */
640 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetNonlinear",
641 &ark_mem, &step_mem);
642 if (retval != ARK_SUCCESS) return(retval);
643
644 /* set parameters */
645 step_mem->linear = SUNFALSE;
646 step_mem->linear_timedep = SUNTRUE;
647 step_mem->dgmax = DGMAX;
648
649 return(ARK_SUCCESS);
650 }
651
652
653 /*---------------------------------------------------------------
654 ARKStepSetExplicit:
655
656 Specifies that the implicit portion of the problem is disabled,
657 and to use an explicit RK method.
658 ---------------------------------------------------------------*/
ARKStepSetExplicit(void * arkode_mem)659 int ARKStepSetExplicit(void *arkode_mem)
660 {
661 ARKodeMem ark_mem;
662 ARKodeARKStepMem step_mem;
663 int retval;
664
665 /* access ARKodeARKStepMem structure */
666 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetExplicit",
667 &ark_mem, &step_mem);
668 if (retval != ARK_SUCCESS) return(retval);
669
670 /* ensure that fe is defined */
671 if (step_mem->fe == NULL) {
672 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
673 "ARKStepSetExplicit", MSG_ARK_MISSING_FE);
674 return(ARK_ILL_INPUT);
675 }
676
677 /* set the relevant parameters */
678 step_mem->explicit = SUNTRUE;
679 step_mem->implicit = SUNFALSE;
680
681 return(ARK_SUCCESS);
682 }
683
684
685 /*---------------------------------------------------------------
686 ARKStepSetImplicit:
687
688 Specifies that the explicit portion of the problem is disabled,
689 and to use an implicit RK method.
690 ---------------------------------------------------------------*/
ARKStepSetImplicit(void * arkode_mem)691 int ARKStepSetImplicit(void *arkode_mem)
692 {
693 ARKodeMem ark_mem;
694 ARKodeARKStepMem step_mem;
695 int retval;
696
697 /* access ARKodeARKStepMem structure */
698 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetImplicit",
699 &ark_mem, &step_mem);
700 if (retval != ARK_SUCCESS) return(retval);
701
702 /* ensure that fi is defined */
703 if (step_mem->fi == NULL) {
704 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
705 "ARKStepSetImplicit", MSG_ARK_MISSING_FI);
706 return(ARK_ILL_INPUT);
707 }
708
709 /* set the relevant parameters */
710 step_mem->implicit = SUNTRUE;
711 step_mem->explicit = SUNFALSE;
712
713 /* re-attach internal error weight functions if necessary */
714 if (!ark_mem->user_efun) {
715 if (ark_mem->itol == ARK_SV && ark_mem->Vabstol != NULL)
716 retval = arkSVtolerances(ark_mem, ark_mem->reltol, ark_mem->Vabstol);
717 else
718 retval = arkSStolerances(ark_mem, ark_mem->reltol, ark_mem->Sabstol);
719 if (retval != ARK_SUCCESS) return(retval);
720 }
721
722 return(ARK_SUCCESS);
723 }
724
725
726 /*---------------------------------------------------------------
727 ARKStepSetImEx:
728
729 Specifies that the specifies that problem has both implicit and
730 explicit parts, and to use an ARK method (this is the default).
731 ---------------------------------------------------------------*/
ARKStepSetImEx(void * arkode_mem)732 int ARKStepSetImEx(void *arkode_mem)
733 {
734 ARKodeMem ark_mem;
735 ARKodeARKStepMem step_mem;
736 int retval;
737
738 /* access ARKodeARKStepMem structure */
739 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetImEx",
740 &ark_mem, &step_mem);
741 if (retval != ARK_SUCCESS) return(retval);
742
743 /* ensure that fe and fi are defined */
744 if (step_mem->fe == NULL) {
745 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
746 "ARKStepSetImEx", MSG_ARK_MISSING_FE);
747 return(ARK_ILL_INPUT);
748 }
749 if (step_mem->fi == NULL) {
750 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
751 "ARKStepSetImEx", MSG_ARK_MISSING_FI);
752 return(ARK_ILL_INPUT);
753 }
754
755 /* set the relevant parameters */
756 step_mem->explicit = SUNTRUE;
757 step_mem->implicit = SUNTRUE;
758
759 /* re-attach internal error weight functions if necessary */
760 if (!ark_mem->user_efun) {
761 if (ark_mem->itol == ARK_SV && ark_mem->Vabstol != NULL)
762 retval = arkSVtolerances(ark_mem, ark_mem->reltol, ark_mem->Vabstol);
763 else
764 retval = arkSStolerances(ark_mem, ark_mem->reltol, ark_mem->Sabstol);
765 if (retval != ARK_SUCCESS) return(retval);
766 }
767
768 return(ARK_SUCCESS);
769 }
770
771
772 /*---------------------------------------------------------------
773 ARKStepSetTables:
774
775 Specifies to use customized Butcher tables for the system.
776
777 If Bi is NULL, then this sets the integrator in 'explicit' mode.
778
779 If Be is NULL, then this sets the integrator in 'implicit' mode.
780
781 Returns ARK_ILL_INPUT if both Butcher tables are not supplied.
782 ---------------------------------------------------------------*/
ARKStepSetTables(void * arkode_mem,int q,int p,ARKodeButcherTable Bi,ARKodeButcherTable Be)783 int ARKStepSetTables(void *arkode_mem, int q, int p,
784 ARKodeButcherTable Bi, ARKodeButcherTable Be)
785 {
786 int retval;
787 ARKodeMem ark_mem;
788 ARKodeARKStepMem step_mem;
789
790 /* access ARKodeARKStepMem structure */
791 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetTables",
792 &ark_mem, &step_mem);
793 if (retval != ARK_SUCCESS) return(retval);
794
795 /* check for illegal inputs */
796 if ((Bi == NULL) && (Be == NULL)) {
797 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
798 "ARKStepSetTables",
799 "At least one complete table must be supplied");
800 return(ARK_ILL_INPUT);
801 }
802
803 /* if both tables are set, check that they have the same number of stages */
804 if ((Bi != NULL) && (Be != NULL)) {
805 if (Bi->stages != Be->stages) {
806 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
807 "ARKStepSetTables",
808 "Both tables must have the same number of stages");
809 return(ARK_ILL_INPUT);
810 }
811 }
812
813 /* clear any existing parameters and Butcher tables */
814 step_mem->stages = 0;
815 step_mem->q = 0;
816 step_mem->p = 0;
817 ARKodeButcherTable_Free(step_mem->Be); step_mem->Be = NULL;
818 ARKodeButcherTable_Free(step_mem->Bi); step_mem->Bi = NULL;
819
820 /*
821 * determine mode (implicit/explicit/ImEx), and perform appropriate actions
822 */
823
824 /* explicit */
825 if (Bi == NULL) {
826
827 /* set the relevant parameters (use table q and p) */
828 step_mem->stages = Be->stages;
829 step_mem->q = Be->q;
830 step_mem->p = Be->p;
831
832 /* copy the table in step memory */
833 step_mem->Be = ARKodeButcherTable_Copy(Be);
834 if (step_mem->Be == NULL) {
835 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
836 "ARKStepSetTables", MSG_ARK_NO_MEM);
837 return(ARK_MEM_NULL);
838 }
839
840 /* set method as purely explicit */
841 retval = ARKStepSetExplicit(arkode_mem);
842 if (retval != ARK_SUCCESS) {
843 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
844 "ARKStepSetTables",
845 "Error in ARKStepSetExplicit");
846 return(retval);
847 }
848
849 /* implicit */
850 } else if (Be == NULL) {
851
852 /* set the relevant parameters (use table q and p) */
853 step_mem->stages = Bi->stages;
854 step_mem->q = Bi->q;
855 step_mem->p = Bi->p;
856
857 /* copy the table in step memory */
858 step_mem->Bi = ARKodeButcherTable_Copy(Bi);
859 if (step_mem->Bi == NULL) {
860 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
861 "ARKStepSetTables", MSG_ARK_NO_MEM);
862 return(ARK_MEM_NULL);
863 }
864
865 /* set method as purely implicit */
866 retval = ARKStepSetImplicit(arkode_mem);
867 if (retval != ARK_SUCCESS) {
868 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
869 "ARKStepSetTables",
870 "Error in ARKStepSetImplicit");
871 return(ARK_ILL_INPUT);
872 }
873
874 /* ImEx */
875 } else {
876
877 /* set the relevant parameters (use input q and p) */
878 step_mem->stages = Bi->stages;
879 step_mem->q = q;
880 step_mem->p = p;
881
882 /* copy the explicit table into step memory */
883 step_mem->Be = ARKodeButcherTable_Copy(Be);
884 if (step_mem->Be == NULL) {
885 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
886 "ARKStepSetTables", MSG_ARK_NO_MEM);
887 return(ARK_MEM_NULL);
888 }
889
890 /* copy the implicit table into step memory */
891 step_mem->Bi = ARKodeButcherTable_Copy(Bi);
892 if (step_mem->Bi == NULL) {
893 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
894 "ARKStepSetTables", MSG_ARK_NO_MEM);
895 return(ARK_MEM_NULL);
896 }
897
898 /* set method as ImEx */
899 retval = ARKStepSetImEx(arkode_mem);
900 if (retval != ARK_SUCCESS) {
901 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
902 "ARKStepSetTables",
903 "Error in ARKStepSetImEx");
904 return(ARK_ILL_INPUT);
905 }
906 }
907
908 return(ARK_SUCCESS);
909 }
910
911
912 /*---------------------------------------------------------------
913 ARKStepSetTableNum:
914
915 Specifies to use pre-existing Butcher tables for the system,
916 based on the integer flags passed to
917 ARKodeButcherTable_LoadERK() and ARKodeButcherTable_LoadDIRK()
918 within the files arkode_butcher_erk.c and arkode_butcher_dirk.c
919 (automatically calls ARKStepSetImEx).
920
921 If either argument is negative (illegal), then this disables the
922 corresponding table (e.g. itable = -1 -> explicit)
923 ---------------------------------------------------------------*/
ARKStepSetTableNum(void * arkode_mem,int itable,int etable)924 int ARKStepSetTableNum(void *arkode_mem, int itable, int etable)
925 {
926 int flag, retval;
927 ARKodeMem ark_mem;
928 ARKodeARKStepMem step_mem;
929
930 /* access ARKodeARKStepMem structure */
931 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetTableNum",
932 &ark_mem, &step_mem);
933 if (retval != ARK_SUCCESS) return(retval);
934
935 /* clear any existing parameters and Butcher tables */
936 step_mem->stages = 0;
937 step_mem->q = 0;
938 step_mem->p = 0;
939 ARKodeButcherTable_Free(step_mem->Be); step_mem->Be = NULL;
940 ARKodeButcherTable_Free(step_mem->Bi); step_mem->Bi = NULL;
941
942
943 /* determine mode (implicit/explicit/ImEx), and perform
944 appropriate actions */
945
946 /* illegal inputs */
947 if ((itable < 0) && (etable < 0)) {
948 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
949 "ARKStepSetTableNum",
950 "At least one valid table number must be supplied");
951 return(ARK_ILL_INPUT);
952
953
954 /* explicit */
955 } else if (itable < 0) {
956
957 /* check that argument specifies an explicit table */
958 if (etable<MIN_ERK_NUM || etable>MAX_ERK_NUM) {
959 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
960 "ARKStepSetTableNum",
961 "Illegal ERK table number");
962 return(ARK_ILL_INPUT);
963 }
964
965 /* fill in table based on argument */
966 step_mem->Be = ARKodeButcherTable_LoadERK(etable);
967 if (step_mem->Be == NULL) {
968 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
969 "ARKStepSetTableNum",
970 "Error setting explicit table with that index");
971 return(ARK_ILL_INPUT);
972 }
973 step_mem->stages = step_mem->Be->stages;
974 step_mem->q = step_mem->Be->q;
975 step_mem->p = step_mem->Be->p;
976
977 /* set method as purely explicit */
978 flag = ARKStepSetExplicit(arkode_mem);
979 if (flag != ARK_SUCCESS) {
980 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
981 "ARKStepSetTableNum",
982 "Error in ARKStepSetExplicit");
983 return(flag);
984 }
985
986
987 /* implicit */
988 } else if (etable < 0) {
989
990 /* check that argument specifies an implicit table */
991 if (itable<MIN_DIRK_NUM || itable>MAX_DIRK_NUM) {
992 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
993 "ARKStepSetTableNum",
994 "Illegal IRK table number");
995 return(ARK_ILL_INPUT);
996 }
997
998 /* fill in table based on argument */
999 step_mem->Bi = ARKodeButcherTable_LoadDIRK(itable);
1000 if (step_mem->Bi == NULL) {
1001 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
1002 "ARKStepSetTableNum",
1003 "Error setting table with that index");
1004 return(ARK_ILL_INPUT);
1005 }
1006 step_mem->stages = step_mem->Bi->stages;
1007 step_mem->q = step_mem->Bi->q;
1008 step_mem->p = step_mem->Bi->p;
1009
1010 /* set method as purely implicit */
1011 flag = ARKStepSetImplicit(arkode_mem);
1012 if (flag != ARK_SUCCESS) {
1013 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
1014 "ARKStepSetTableNum",
1015 "Error in ARKStepSetIxplicit");
1016 return(flag);
1017 }
1018
1019
1020 /* ImEx */
1021 } else {
1022
1023 /* ensure that tables match */
1024 if ( !((etable == ARK324L2SA_ERK_4_2_3) && (itable == ARK324L2SA_DIRK_4_2_3)) &&
1025 !((etable == ARK436L2SA_ERK_6_3_4) && (itable == ARK436L2SA_DIRK_6_3_4)) &&
1026 !((etable == ARK437L2SA_ERK_7_3_4) && (itable == ARK437L2SA_DIRK_7_3_4)) &&
1027 !((etable == ARK548L2SA_ERK_8_4_5) && (itable == ARK548L2SA_DIRK_8_4_5)) &&
1028 !((etable == ARK548L2SAb_ERK_8_4_5) && (itable == ARK548L2SAb_DIRK_8_4_5)) ) {
1029 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
1030 "ARKStepSetTableNum",
1031 "Incompatible Butcher tables for ARK method");
1032 return(ARK_ILL_INPUT);
1033 }
1034
1035 /* fill in tables based on arguments */
1036 step_mem->Bi = ARKodeButcherTable_LoadDIRK(itable);
1037 step_mem->Be = ARKodeButcherTable_LoadERK(etable);
1038 if (step_mem->Bi == NULL) {
1039 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
1040 "ARKStepSetTableNum",
1041 "Illegal IRK table number");
1042 return(ARK_ILL_INPUT);
1043 }
1044 if (step_mem->Be == NULL) {
1045 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
1046 "ARKStepSetTableNum",
1047 "Illegal ERK table number");
1048 return(ARK_ILL_INPUT);
1049 }
1050 step_mem->stages = step_mem->Bi->stages;
1051 step_mem->q = step_mem->Bi->q;
1052 step_mem->p = step_mem->Bi->p;
1053
1054 /* set method as ImEx */
1055 if (ARKStepSetImEx(arkode_mem) != ARK_SUCCESS) {
1056 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep",
1057 "ARKStepSetTableNum", MSG_ARK_MISSING_F);
1058 return(ARK_ILL_INPUT);
1059 }
1060
1061 }
1062
1063 return(ARK_SUCCESS);
1064 }
1065
1066
1067 /*---------------------------------------------------------------
1068 ARKStepSetNonlinCRDown:
1069
1070 Specifies the user-provided nonlinear convergence constant
1071 crdown. Legal values are strictly positive; illegal values
1072 imply a reset to the default.
1073 ---------------------------------------------------------------*/
ARKStepSetNonlinCRDown(void * arkode_mem,realtype crdown)1074 int ARKStepSetNonlinCRDown(void *arkode_mem, realtype crdown)
1075 {
1076 ARKodeMem ark_mem;
1077 ARKodeARKStepMem step_mem;
1078 int retval;
1079
1080 /* access ARKodeARKStepMem structure */
1081 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetNonlinCRDown",
1082 &ark_mem, &step_mem);
1083 if (retval != ARK_SUCCESS) return(retval);
1084
1085 /* if argument legal set it, otherwise set default */
1086 if (crdown <= ZERO) {
1087 step_mem->crdown = CRDOWN;
1088 } else {
1089 step_mem->crdown = crdown;
1090 }
1091
1092 return(ARK_SUCCESS);
1093 }
1094
1095
1096 /*---------------------------------------------------------------
1097 ARKStepSetNonlinRDiv:
1098
1099 Specifies the user-provided nonlinear convergence constant
1100 rdiv. Legal values are strictly positive; illegal values
1101 imply a reset to the default.
1102 ---------------------------------------------------------------*/
ARKStepSetNonlinRDiv(void * arkode_mem,realtype rdiv)1103 int ARKStepSetNonlinRDiv(void *arkode_mem, realtype rdiv)
1104 {
1105 ARKodeMem ark_mem;
1106 ARKodeARKStepMem step_mem;
1107 int retval;
1108
1109 /* access ARKodeARKStepMem structure */
1110 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetNonlinRDiv",
1111 &ark_mem, &step_mem);
1112 if (retval != ARK_SUCCESS) return(retval);
1113
1114 /* if argument legal set it, otherwise set default */
1115 if (rdiv <= ZERO) {
1116 step_mem->rdiv = RDIV;
1117 } else {
1118 step_mem->rdiv = rdiv;
1119 }
1120
1121 return(ARK_SUCCESS);
1122 }
1123
1124
1125 /*---------------------------------------------------------------
1126 ARKStepSetDeltaGammaMax:
1127
1128 Specifies the user-provided linear setup decision constant
1129 dgmax. Legal values are strictly positive; illegal values imply
1130 a reset to the default.
1131 ---------------------------------------------------------------*/
ARKStepSetDeltaGammaMax(void * arkode_mem,realtype dgmax)1132 int ARKStepSetDeltaGammaMax(void *arkode_mem, realtype dgmax)
1133 {
1134 ARKodeMem ark_mem;
1135 ARKodeARKStepMem step_mem;
1136 int retval;
1137
1138 /* access ARKodeARKStepMem structure */
1139 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetDeltaGammaMax",
1140 &ark_mem, &step_mem);
1141 if (retval != ARK_SUCCESS) return(retval);
1142
1143 /* if argument legal set it, otherwise set default */
1144 if (dgmax <= ZERO) {
1145 step_mem->dgmax = DGMAX;
1146 } else {
1147 step_mem->dgmax = dgmax;
1148 }
1149
1150 return(ARK_SUCCESS);
1151 }
1152
1153
1154 /*---------------------------------------------------------------
1155 ARKStepSetMaxStepsBetweenLSet:
1156
1157 Specifies the user-provided linear setup decision constant
1158 msbp. Positive values give the number of time steps to wait
1159 before calling lsetup; negative values imply recomputation of
1160 lsetup at each nonlinear solve; a zero value implies a reset
1161 to the default.
1162 ---------------------------------------------------------------*/
ARKStepSetMaxStepsBetweenLSet(void * arkode_mem,int msbp)1163 int ARKStepSetMaxStepsBetweenLSet(void *arkode_mem, int msbp)
1164 {
1165 ARKodeMem ark_mem;
1166 ARKodeARKStepMem step_mem;
1167 int retval;
1168
1169 /* access ARKodeARKStepMem structure */
1170 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetMaxStepsBetweenLSet",
1171 &ark_mem, &step_mem);
1172 if (retval != ARK_SUCCESS) return(retval);
1173
1174 /* if argument legal set it, otherwise set default */
1175 if (msbp == 0) {
1176 step_mem->msbp = MSBP;
1177 } else {
1178 step_mem->msbp = msbp;
1179 }
1180
1181 return(ARK_SUCCESS);
1182 }
1183
1184
1185 /*---------------------------------------------------------------
1186 ARKStepSetPredictorMethod:
1187
1188 Specifies the method to use for predicting implicit solutions.
1189 Non-default choices are {1,2,3,4}, all others will use default
1190 (trivial) predictor.
1191 ---------------------------------------------------------------*/
ARKStepSetPredictorMethod(void * arkode_mem,int pred_method)1192 int ARKStepSetPredictorMethod(void *arkode_mem, int pred_method)
1193 {
1194 ARKodeMem ark_mem;
1195 ARKodeARKStepMem step_mem;
1196 int retval;
1197
1198 /* access ARKodeARKStepMem structure */
1199 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetPredictorMethod",
1200 &ark_mem, &step_mem);
1201 if (retval != ARK_SUCCESS) return(retval);
1202
1203 /* return error if pred_method==5 and a non-NULL stage predictor function
1204 has been supplied */
1205 if ((pred_method == 5) && (step_mem->stage_predict != NULL)) {
1206 arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKode::ARKStep", "ARKStepSetPredictorMethod",
1207 "predictor 5 cannot be combined with user-supplied stage predictor");
1208 return(ARK_ILL_INPUT);
1209 }
1210
1211 /* set parameter */
1212 step_mem->predictor = pred_method;
1213
1214 return(ARK_SUCCESS);
1215 }
1216
1217
1218 /*---------------------------------------------------------------
1219 ARKStepSetMaxNonlinIters:
1220
1221 Specifies the maximum number of nonlinear iterations during
1222 one solve. A non-positive input implies a reset to the
1223 default value.
1224 ---------------------------------------------------------------*/
ARKStepSetMaxNonlinIters(void * arkode_mem,int maxcor)1225 int ARKStepSetMaxNonlinIters(void *arkode_mem, int maxcor)
1226 {
1227 ARKodeMem ark_mem;
1228 ARKodeARKStepMem step_mem;
1229 int retval;
1230
1231 /* access ARKodeARKStepMem structure */
1232 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetMaxNonlinIters",
1233 &ark_mem, &step_mem);
1234 if (retval != ARK_SUCCESS) return(retval);
1235
1236 /* Return error message if no NLS module is present */
1237 if (step_mem->NLS == NULL) {
1238 arkProcessError(ark_mem, ARK_NLS_OP_ERR, "ARKode::ARKStep",
1239 "ARKStepSetMaxNonlinIters",
1240 "No SUNNonlinearSolver object is present");
1241 return(ARK_ILL_INPUT);
1242 }
1243
1244 /* argument <= 0 sets default, otherwise set input */
1245 if (maxcor <= 0) {
1246 step_mem->maxcor = MAXCOR;
1247 } else {
1248 step_mem->maxcor = maxcor;
1249 }
1250
1251 /* send argument to NLS structure */
1252 retval = SUNNonlinSolSetMaxIters(step_mem->NLS, step_mem->maxcor);
1253 if (retval != SUN_NLS_SUCCESS) {
1254 arkProcessError(ark_mem, ARK_NLS_OP_ERR, "ARKode::ARKStep",
1255 "ARKStepSetMaxNonlinIters",
1256 "Error setting maxcor in SUNNonlinearSolver object");
1257 return(ARK_NLS_OP_ERR);
1258 }
1259
1260 return(ARK_SUCCESS);
1261 }
1262
1263
1264 /*---------------------------------------------------------------
1265 ARKStepSetNonlinConvCoef:
1266
1267 Specifies the coefficient in the nonlinear solver convergence
1268 test. A non-positive input implies a reset to the default value.
1269 ---------------------------------------------------------------*/
ARKStepSetNonlinConvCoef(void * arkode_mem,realtype nlscoef)1270 int ARKStepSetNonlinConvCoef(void *arkode_mem, realtype nlscoef)
1271 {
1272 ARKodeMem ark_mem;
1273 ARKodeARKStepMem step_mem;
1274 int retval;
1275
1276 /* access ARKodeARKStepMem structure */
1277 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepSetNonlinConvCoef",
1278 &ark_mem, &step_mem);
1279 if (retval != ARK_SUCCESS) return(retval);
1280
1281 /* argument <= 0 sets default, otherwise set input */
1282 if (nlscoef <= ZERO) {
1283 step_mem->nlscoef = NLSCOEF;
1284 } else {
1285 step_mem->nlscoef = nlscoef;
1286 }
1287
1288 return(ARK_SUCCESS);
1289 }
1290
1291
1292 /*===============================================================
1293 ARKStep optional output functions -- stepper-specific
1294 ===============================================================*/
1295
1296 /*---------------------------------------------------------------
1297 ARKStepGetCurrentGamma: Returns the current value of gamma
1298 ---------------------------------------------------------------*/
ARKStepGetCurrentGamma(void * arkode_mem,realtype * gamma)1299 int ARKStepGetCurrentGamma(void *arkode_mem, realtype *gamma)
1300 {
1301 int retval;
1302 ARKodeMem ark_mem;
1303 ARKodeARKStepMem step_mem;
1304 retval = arkStep_AccessStepMem(arkode_mem, NULL, &ark_mem, &step_mem);
1305 if (retval != ARK_SUCCESS) return(retval);
1306 *gamma = step_mem->gamma;
1307 return(retval);
1308 }
1309
1310
1311 /*---------------------------------------------------------------
1312 ARKStepGetNumRhsEvals:
1313
1314 Returns the current number of calls to fe and fi
1315 ---------------------------------------------------------------*/
ARKStepGetNumRhsEvals(void * arkode_mem,long int * fe_evals,long int * fi_evals)1316 int ARKStepGetNumRhsEvals(void *arkode_mem, long int *fe_evals,
1317 long int *fi_evals)
1318 {
1319 ARKodeMem ark_mem;
1320 ARKodeARKStepMem step_mem;
1321 int retval;
1322
1323 /* access ARKodeARKStepMem structure */
1324 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetNumRhsEvals",
1325 &ark_mem, &step_mem);
1326 if (retval != ARK_SUCCESS) return(retval);
1327
1328 /* get values from step_mem */
1329 *fe_evals = step_mem->nfe;
1330 *fi_evals = step_mem->nfi;
1331
1332 return(ARK_SUCCESS);
1333 }
1334
1335
1336 /*---------------------------------------------------------------
1337 ARKStepGetNumLinSolvSetups:
1338
1339 Returns the current number of calls to the lsetup routine
1340 ---------------------------------------------------------------*/
ARKStepGetNumLinSolvSetups(void * arkode_mem,long int * nlinsetups)1341 int ARKStepGetNumLinSolvSetups(void *arkode_mem, long int *nlinsetups)
1342 {
1343 ARKodeMem ark_mem;
1344 ARKodeARKStepMem step_mem;
1345 int retval;
1346
1347 /* access ARKodeARKStepMem structure */
1348 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetNumLinSolvSetups",
1349 &ark_mem, &step_mem);
1350 if (retval != ARK_SUCCESS) return(retval);
1351
1352 /* get value from step_mem */
1353 *nlinsetups = step_mem->nsetups;
1354
1355 return(ARK_SUCCESS);
1356 }
1357
1358
1359 /*---------------------------------------------------------------
1360 ARKStepGetCurrentButcherTables:
1361
1362 Sets pointers to the explicit and implicit Butcher tables
1363 currently in use.
1364 ---------------------------------------------------------------*/
ARKStepGetCurrentButcherTables(void * arkode_mem,ARKodeButcherTable * Bi,ARKodeButcherTable * Be)1365 int ARKStepGetCurrentButcherTables(void *arkode_mem,
1366 ARKodeButcherTable *Bi,
1367 ARKodeButcherTable *Be)
1368 {
1369 ARKodeMem ark_mem;
1370 ARKodeARKStepMem step_mem;
1371 int retval;
1372
1373 /* access ARKodeARKStepMem structure */
1374 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetCurrentButcherTables",
1375 &ark_mem, &step_mem);
1376 if (retval != ARK_SUCCESS) return(retval);
1377
1378 /* get tables from step_mem */
1379 *Bi = step_mem->Bi;
1380 *Be = step_mem->Be;
1381 return(ARK_SUCCESS);
1382 }
1383
1384
1385 /*---------------------------------------------------------------
1386 ARKStepGetEstLocalErrors: (updated to the correct vector, but
1387 need to verify that it is unchanged between filling the
1388 estimated error and the end of the time step)
1389
1390 Returns an estimate of the local error
1391 ---------------------------------------------------------------*/
ARKStepGetEstLocalErrors(void * arkode_mem,N_Vector ele)1392 int ARKStepGetEstLocalErrors(void *arkode_mem, N_Vector ele)
1393 {
1394 ARKodeMem ark_mem;
1395 ARKodeARKStepMem step_mem;
1396 int retval;
1397
1398 /* access ARKodeARKStepMem structure */
1399 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetEstLocalErrors",
1400 &ark_mem, &step_mem);
1401 if (retval != ARK_SUCCESS) return(retval);
1402
1403 /* copy vector to output */
1404 N_VScale(ONE, ark_mem->tempv1, ele);
1405
1406 return(ARK_SUCCESS);
1407 }
1408
1409
1410 /*---------------------------------------------------------------
1411 ARKStepGetTimestepperStats:
1412
1413 Returns integrator statistics
1414 ---------------------------------------------------------------*/
ARKStepGetTimestepperStats(void * arkode_mem,long int * expsteps,long int * accsteps,long int * step_attempts,long int * fe_evals,long int * fi_evals,long int * nlinsetups,long int * netfails)1415 int ARKStepGetTimestepperStats(void *arkode_mem, long int *expsteps,
1416 long int *accsteps, long int *step_attempts,
1417 long int *fe_evals, long int *fi_evals,
1418 long int *nlinsetups, long int *netfails)
1419 {
1420 ARKodeMem ark_mem;
1421 ARKodeARKStepMem step_mem;
1422 int retval;
1423
1424 /* access ARKodeARKStepMem structure */
1425 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetTimestepperStats",
1426 &ark_mem, &step_mem);
1427 if (retval != ARK_SUCCESS) return(retval);
1428
1429 /* set expsteps and accsteps from adaptivity structure */
1430 *expsteps = ark_mem->hadapt_mem->nst_exp;
1431 *accsteps = ark_mem->hadapt_mem->nst_acc;
1432
1433 /* set remaining outputs */
1434 *step_attempts = ark_mem->nst_attempts;
1435 *fe_evals = step_mem->nfe;
1436 *fi_evals = step_mem->nfi;
1437 *nlinsetups = step_mem->nsetups;
1438 *netfails = ark_mem->netf;
1439
1440 return(ARK_SUCCESS);
1441 }
1442
1443
1444 /*---------------------------------------------------------------
1445 ARKStepGetNumNonlinSolvIters:
1446
1447 Returns the current number of nonlinear solver iterations
1448 ---------------------------------------------------------------*/
ARKStepGetNumNonlinSolvIters(void * arkode_mem,long int * nniters)1449 int ARKStepGetNumNonlinSolvIters(void *arkode_mem, long int *nniters)
1450 {
1451 ARKodeMem ark_mem;
1452 ARKodeARKStepMem step_mem;
1453 int retval;
1454
1455 /* access ARKodeARKStepMem structure */
1456 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetNumNonlinSolvIters",
1457 &ark_mem, &step_mem);
1458 if (retval != ARK_SUCCESS) return(retval);
1459
1460 /* if a NLS object is present, set output from that; otherwise
1461 we took zero iterations */
1462 if (step_mem->NLS) {
1463 retval = SUNNonlinSolGetNumIters(step_mem->NLS, nniters);
1464 if (retval != SUN_NLS_SUCCESS) {
1465 arkProcessError(ark_mem, ARK_NLS_OP_ERR, "ARKode::ARKStep",
1466 "ARKStepGetNumNonlinSolvIters",
1467 "Error retrieving nniters from SUNNonlinearSolver");
1468 return(ARK_NLS_OP_ERR);
1469 }
1470 } else {
1471 *nniters = 0;
1472 }
1473
1474 return(ARK_SUCCESS);
1475 }
1476
1477
1478 /*---------------------------------------------------------------
1479 ARKStepGetNumNonlinSolvConvFails:
1480
1481 Returns the current number of nonlinear solver convergence fails
1482 ---------------------------------------------------------------*/
ARKStepGetNumNonlinSolvConvFails(void * arkode_mem,long int * nncfails)1483 int ARKStepGetNumNonlinSolvConvFails(void *arkode_mem, long int *nncfails)
1484 {
1485 ARKodeMem ark_mem;
1486 ARKodeARKStepMem step_mem;
1487 int retval;
1488
1489 /* access ARKodeARKStepMem structure */
1490 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetNumNonlinSolvConvFails",
1491 &ark_mem, &step_mem);
1492 if (retval != ARK_SUCCESS) return(retval);
1493
1494 /* set output from step_mem */
1495 *nncfails = ark_mem->ncfn;
1496
1497 return(ARK_SUCCESS);
1498 }
1499
1500
1501 /*---------------------------------------------------------------
1502 ARKStepGetNonlinSolvStats:
1503
1504 Returns nonlinear solver statistics
1505 ---------------------------------------------------------------*/
ARKStepGetNonlinSolvStats(void * arkode_mem,long int * nniters,long int * nncfails)1506 int ARKStepGetNonlinSolvStats(void *arkode_mem, long int *nniters,
1507 long int *nncfails)
1508 {
1509 ARKodeMem ark_mem;
1510 ARKodeARKStepMem step_mem;
1511 int retval;
1512
1513 /* access ARKodeARKStepMem structure */
1514 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepGetNonlinSolvStats",
1515 &ark_mem, &step_mem);
1516 if (retval != ARK_SUCCESS) return(retval);
1517
1518 /* set outputs from NLS module and step_mem structure (if present);
1519 otherwise there were zero iterations and no nonlinear failures */
1520 if (step_mem->NLS) {
1521 retval = SUNNonlinSolGetNumIters(step_mem->NLS, nniters);
1522 if (retval != SUN_NLS_SUCCESS) {
1523 arkProcessError(ark_mem, ARK_NLS_OP_ERR, "ARKode::ARKStep",
1524 "ARKStepGetNonlinSolvStats",
1525 "Error retrieving nniters from SUNNonlinearSolver");
1526 return(ARK_NLS_OP_ERR);
1527 }
1528 *nncfails = ark_mem->ncfn;
1529 } else {
1530 *nniters = 0;
1531 *nncfails = 0;
1532 }
1533
1534 return(ARK_SUCCESS);
1535 }
1536
1537
1538 /*===============================================================
1539 ARKStep parameter output
1540 ===============================================================*/
1541
1542 /*---------------------------------------------------------------
1543 ARKStepWriteParameters:
1544
1545 Outputs all solver parameters to the provided file pointer.
1546 ---------------------------------------------------------------*/
ARKStepWriteParameters(void * arkode_mem,FILE * fp)1547 int ARKStepWriteParameters(void *arkode_mem, FILE *fp)
1548 {
1549 ARKodeMem ark_mem;
1550 ARKodeARKStepMem step_mem;
1551 int flag, retval;
1552
1553 /* access ARKodeARKStepMem structure */
1554 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepWriteParameters",
1555 &ark_mem, &step_mem);
1556 if (retval != ARK_SUCCESS) return(retval);
1557
1558 /* output ARKode infrastructure parameters first */
1559 flag = arkWriteParameters(ark_mem, fp);
1560 if (flag != ARK_SUCCESS) {
1561 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
1562 "ARKStepWriteParameters",
1563 "Error writing ARKode infrastructure parameters");
1564 return(flag);
1565 }
1566
1567 /* print integrator parameters to file */
1568 fprintf(fp, "ARKStep time step module parameters:\n");
1569 fprintf(fp, " Method order %i\n",step_mem->q);
1570 if (step_mem->linear) {
1571 fprintf(fp, " Linear implicit problem");
1572 if (step_mem->linear_timedep) {
1573 fprintf(fp, " (time-dependent Jacobian)\n");
1574 } else {
1575 fprintf(fp, " (time-independent Jacobian)\n");
1576 }
1577 }
1578 if (step_mem->explicit && step_mem->implicit) {
1579 fprintf(fp, " ImEx integrator\n");
1580 } else if (step_mem->implicit) {
1581 fprintf(fp, " Implicit integrator\n");
1582 } else {
1583 fprintf(fp, " Explicit integrator\n");
1584 }
1585
1586 if (step_mem->implicit) {
1587 fprintf(fp, " Implicit predictor method = %i\n",step_mem->predictor);
1588 fprintf(fp, " Implicit solver tolerance coefficient = %"RSYM"\n",step_mem->nlscoef);
1589 fprintf(fp, " Maximum number of nonlinear corrections = %i\n",step_mem->maxcor);
1590 fprintf(fp, " Nonlinear convergence rate constant = %"RSYM"\n",step_mem->crdown);
1591 fprintf(fp, " Nonlinear divergence tolerance = %"RSYM"\n",step_mem->rdiv);
1592 fprintf(fp, " Gamma factor LSetup tolerance = %"RSYM"\n",step_mem->dgmax);
1593 fprintf(fp, " Number of steps between LSetup calls = %i\n",step_mem->msbp);
1594 }
1595 fprintf(fp, "\n");
1596
1597 return(ARK_SUCCESS);
1598 }
1599
1600
1601 /*---------------------------------------------------------------
1602 ARKStepWriteButcher:
1603
1604 Outputs Butcher tables to the provided file pointer.
1605 ---------------------------------------------------------------*/
ARKStepWriteButcher(void * arkode_mem,FILE * fp)1606 int ARKStepWriteButcher(void *arkode_mem, FILE *fp)
1607 {
1608 int retval;
1609 ARKodeMem ark_mem;
1610 ARKodeARKStepMem step_mem;
1611
1612 /* access ARKodeARKStepMem structure */
1613 retval = arkStep_AccessStepMem(arkode_mem, "ARKStepWriteButcher",
1614 &ark_mem, &step_mem);
1615 if (retval != ARK_SUCCESS) return(retval);
1616
1617 /* check that Butcher table is non-NULL (otherwise report error) */
1618 if ((step_mem->Be == NULL) && (step_mem->Bi == NULL)) {
1619 arkProcessError(ark_mem, ARK_MEM_NULL, "ARKode::ARKStep",
1620 "ARKStepWriteButcher", "Butcher table memory is NULL");
1621 return(ARK_MEM_NULL);
1622 }
1623
1624 /* print Butcher tables to file */
1625 fprintf(fp, "\nARKStep Butcher tables (stages = %i):\n", step_mem->stages);
1626 if (step_mem->explicit && (step_mem->Be != NULL)) {
1627 fprintf(fp, " Explicit Butcher table:\n");
1628 ARKodeButcherTable_Write(step_mem->Be, fp);
1629 }
1630 fprintf(fp, "\n");
1631 if (step_mem->implicit && (step_mem->Bi != NULL)) {
1632 fprintf(fp, " Implicit Butcher table:\n");
1633 ARKodeButcherTable_Write(step_mem->Bi, fp);
1634 }
1635 fprintf(fp, "\n");
1636
1637 return(ARK_SUCCESS);
1638 }
1639
1640
1641 /*---------------------------------------------------------------
1642 EOF
1643 ---------------------------------------------------------------*/
1644