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