Skip to content

Commit

Permalink
add anderson
Browse files Browse the repository at this point in the history
  • Loading branch information
@balos1
balos1 committed May 17, 2024
1 parent 285b906 commit 3d06b93
Show file tree
Hide file tree
Showing 6 changed files with 74 additions and 21 deletions.
19 changes: 14 additions & 5 deletions examples/cvode/serial/cvBrusselator.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,8 +115,8 @@ int main(int argc, char *argv[])
{
const sunrealtype T0 = SUN_RCONST(0.0); /* initial time */
const sunrealtype Tf = SUN_RCONST(10.0); /* final time */
const sunrealtype reltol = 1.0e-15; /* relative integrator tolerance */
const sunrealtype abstol = 1.0e-15; /* absolute integrator tolerance */
const sunrealtype reltol = 1.0e-6; /* relative integrator tolerance */
const sunrealtype abstol = 1.0e-10; /* absolute integrator tolerance */

int retval;
SUNContext sunctx;
Expand All @@ -140,6 +140,7 @@ int main(int argc, char *argv[])
int nonlinear_solver_type = 0;
int linear_solver_type = 0;
int reactor_type = 2;
int anderson_m = 2;

/* Parse command line arguments and setup UserData */
int argi = 0;
Expand All @@ -152,11 +153,14 @@ int main(int argc, char *argv[])
if (argc > 3) {
linear_solver_type = atoi(argv[++argi]);
}
sunrealtype dTout = SUN_RCONST(1.0);
if (argc > 4) {
anderson_m = atoi(argv[++argi]);
}
sunrealtype dTout = SUN_RCONST(1.0);
if (argc > 5) {
dTout = atof(argv[++argi]);
}
int Nt = (int) ceil(Tf/dTout);
const int Nt = (int) ceil(Tf/dTout);

/* Set the Reaction parameters according to reactor_type */
UserData udata;
Expand Down Expand Up @@ -202,7 +206,7 @@ int main(int argc, char *argv[])
if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) { return(1); }

/* Must be called before CVodeInit */
retval = CVodeSetNonlinearSolverAlgorithm(cvode_mem, nonlinear_solver_type);
retval = CVodeSetNonlinearSolverAlgorithm(cvode_mem, nonlinear_solver_type, anderson_m);
if (check_retval(&retval, "CVodeSetNonlinearSolverAlgorithm", 1)) { return(1); }

/* Call CVodeInit to initialize the integrator memory and specify the
Expand All @@ -220,6 +224,11 @@ int main(int argc, char *argv[])
retval = CVodeSStolerances(cvode_mem, reltol, abstol);
if (check_retval(&retval, "CVodeSVtolerances", 1)) { return(1); }

/* Set an initial time step size to avoid differences when changing dtout.
The switching algorithm seems to be very sensitive to small differences
in initial step for this problem. */
CVodeSetInitStep(cvode_mem, 1e-6);

/* Create the SUNLinearSolver object for use by CVode */
if (linear_solver_type == 0) {
LS = SUNLinSol_SPGMR(y, SUN_PREC_NONE, 0, sunctx);
Expand Down
2 changes: 1 addition & 1 deletion include/cvode/cvode.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ SUNDIALS_EXPORT int CVodeSVtolerances(void* cvode_mem, sunrealtype reltol,
SUNDIALS_EXPORT int CVodeWFtolerances(void* cvode_mem, CVEwtFn efun);

/* Optional input functions */
SUNDIALS_EXPORT int CVodeSetNonlinearSolverAlgorithm(void* cvode_mem, int algorithm);
SUNDIALS_EXPORT int CVodeSetNonlinearSolverAlgorithm(void* cvode_mem, int algorithm, int m);
SUNDIALS_EXPORT int CVodeSetConstraints(void* cvode_mem, N_Vector constraints);
SUNDIALS_EXPORT int CVodeSetDeltaGammaMaxLSetup(void* cvode_mem,
sunrealtype dgmax_lsetup);
Expand Down
51 changes: 39 additions & 12 deletions src/cvode/cvode.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -382,6 +382,7 @@ void* CVodeCreate(int lmm, SUNContext sunctx)
/* Initialize nonlinear solver variables */
cv_mem->NLS = NULL;
cv_mem->ownNLS = SUNFALSE;
cv_mem->NLS_aavectors = 2;

/* Initialize fused operations variable */
cv_mem->cv_usefused = SUNFALSE;
Expand All @@ -390,11 +391,16 @@ void* CVodeCreate(int lmm, SUNContext sunctx)
cv_mem->cv_alpharef = ALPHAREF;
cv_mem->cv_lsodkr_strategy = 0;
cv_mem->cv_gustafsoder_strategy = 0;
cv_mem->cv_force_next_step = 0;
cv_mem->cv_stifr = INIT_STIFR;
cv_mem->cv_nst_switchtonewton_considered = 0;
cv_mem->cv_nst_switchtofixed_considered = 0;
cv_mem->cv_nst_switchtonewton_met = 0;
cv_mem->cv_nst_switchtofixed_met = 0;
cv_mem->switchtofixed_min_met = 3;
cv_mem->switchtonewton_min_met = 3;
cv_mem->switchtofixed_delay = 5;
cv_mem->switchtonewton_delay = 5;

/* Return pointer to CVODE memory block */

Expand Down Expand Up @@ -503,7 +509,7 @@ int CVodeInit(void* cvode_mem, CVRhsFn f, sunrealtype t0, N_Vector y0)
return(CV_MEM_FAIL);
}

cv_mem->NLS_fixedpoint = SUNNonlinSol_FixedPoint(y0, 0, cv_mem->cv_sunctx);
cv_mem->NLS_fixedpoint = SUNNonlinSol_FixedPoint(y0, cv_mem->NLS_aavectors, cv_mem->cv_sunctx);

/* check that nonlinear solver is non-NULL */
if (cv_mem->NLS_fixedpoint == NULL) {
Expand Down Expand Up @@ -2231,6 +2237,9 @@ static int cvHin(CVodeMem cv_mem, sunrealtype tout)

hs = hg; /* safeguard against 'uninitialized variable' warning */

SUNLogDebug(CV_LOGGER, __func__, "outer-loop", "hlb = %.16g, hub = %.16g, hg = %.16g",
hlb, hub, hg);

for (count1 = 1; count1 <= MAX_ITERS; count1++)
{
/* Attempts to estimate ydd */
Expand All @@ -2241,6 +2250,8 @@ static int cvHin(CVodeMem cv_mem, sunrealtype tout)
{
hgs = hg * sign;
retval = cvYddNorm(cv_mem, hgs, &yddnrm);
SUNLogExtraDebug(CV_LOGGER, __func__, "estimate-ydd", "hgs = %.16g, yddnrm = %.16g",
hgs, yddnrm);
/* If the RHS function failed unrecoverably, give up */
if (retval < 0) { return (CV_RHSFUNC_FAIL); }
/* If successful, we can use ydd */
Expand Down Expand Up @@ -2299,6 +2310,8 @@ static int cvHin(CVodeMem cv_mem, sunrealtype tout)
if (sign == -1) { h0 = -h0; }
cv_mem->cv_h = h0;

SUNLogDebug(CV_LOGGER, __func__, "exit", "hnew = %.16g, h = %.16g", hnew, h0);

return (CV_SUCCESS);
}

Expand Down Expand Up @@ -2528,7 +2541,7 @@ void cvChooseNlsGustafSoder1(CVodeMem cv_mem)
{
if (SUNNonlinSolGetType(cv_mem->NLS) == SUNNONLINEARSOLVER_FIXEDPOINT)
{
if (cv_mem->cv_nst_switchtonewton_considered > 5) {
if (cv_mem->cv_nst_switchtonewton_considered > cv_mem->switchtonewton_delay) {

/* TODO(CJB): should we choose zi based on the method? Soderlind and Gustafsson
suggest that 2 is fine, but a better value could be found based on the method. */
Expand All @@ -2537,9 +2550,10 @@ void cvChooseNlsGustafSoder1(CVodeMem cv_mem)
SUNLogDebug(CV_LOGGER, __func__, "maybe-switch-to-newt", "halpha = %.16g, hprime = %.16g",
cv_mem->cv_halpha, cv_mem->cv_hprime);
if (switch_to_newton) {
if (cv_mem->cv_nst_switchtonewton_met >= 3)
if (cv_mem->cv_nst_switchtonewton_met >= cv_mem->switchtonewton_min_met)
{
SUNLogDebug(CV_LOGGER, __func__, "maybe-switch-to-newton", "switch to Newton", "");
SUNLogDebug(CV_LOGGER, __func__, "switch-to-newton", "nst_switchtonewton_considered = %d",
cv_mem->cv_nst_switchtonewton_considered);
cvNlsSwitch(cv_mem, cv_mem->NLS_newton);
cv_mem->cv_nst_switchtonewton_met = 0;
cv_mem->cv_nst_switchtonewton_considered = 0;
Expand All @@ -2555,17 +2569,18 @@ void cvChooseNlsGustafSoder1(CVodeMem cv_mem)
}
else if (SUNNonlinSolGetType(cv_mem->NLS) == SUNNONLINEARSOLVER_ROOTFIND)
{
if (cv_mem->cv_nst_switchtofixed_considered > 5) {
if (cv_mem->cv_nst_switchtofixed_considered > cv_mem->switchtofixed_delay) {
SUNLogDebug(CV_LOGGER, __func__, "maybe-switch-to-fp", "stiff = %.16g", cv_mem->cv_stiff);
sunbooleantype switch_to_fixedpoint = cv_mem->cv_stiff < SUN_RCONST(2.0);
if (switch_to_fixedpoint) {
if (cv_mem->cv_nst_switchtofixed_met >= 3)
if (cv_mem->cv_nst_switchtofixed_met >= cv_mem->switchtofixed_min_met)
{
SUNLogDebug(CV_LOGGER, __func__, "maybe-switch-to-fp", "switch to fixed point", "");
cvNlsSwitch(cv_mem, cv_mem->NLS_fixedpoint);
// // Gustafsson and Soderlind suggest reducing step size by alpahref if switching to fixed point.
// cv_mem->cv_hprime = cv_mem->cv_h * cv_mem->cv_alpharef;
// cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h;
// Gustafsson and Soderlind suggest setting hprime to h*alpharef after switching to fixed-point
// However in my testing this seems to result in a lot of extra steps
cv_mem->cv_force_next_step = 0;
SUNLogDebug(CV_LOGGER, __func__, "switch-to-fixed-point", "cv_nst_switchtofixed_considered = %d, hprime = %.16g, eta = %.16g",
cv_mem->cv_nst_switchtofixed_considered, cv_mem->cv_hprime, cv_mem->cv_eta);
cv_mem->cv_nst_switchtofixed_met = 0;
cv_mem->cv_nst_switchtofixed_considered = 0;
}
Expand Down Expand Up @@ -2864,7 +2879,7 @@ static void cvPredict(CVodeMem cv_mem)

#ifdef SUNDIALS_LOGGING_EXTRA_DEBUG
SUNLogger_QueueMsg(CV_LOGGER, SUN_LOGLEVEL_DEBUG, "CVODE::cvPredict",
"return", "predictor =", "");
"return", "predictor(:) =", "");
N_VPrintFile(cv_mem->cv_zn[0], CV_LOGGER->debug_fp);
#endif
}
Expand Down Expand Up @@ -3709,10 +3724,18 @@ static void cvOkCompleteStep(CVodeMem cv_mem)
static void cvOkPrepareNextStep(CVodeMem cv_mem, sunrealtype dsm)
{
/* If etamax = 1, defer step size or order changes */
if (cv_mem->cv_etamax == ONE)
if (cv_mem->cv_force_next_step)
{
cv_mem->cv_hprime = cv_mem->cv_h * cv_mem->cv_alpharef;
cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h;
SUNLogDebug(CV_LOGGER, __func__, "force-next-step", "hprime = %.16g, eta = %.16g", cv_mem->cv_hprime, cv_mem->cv_eta);
cv_mem->cv_force_next_step = 0;
}
else if (cv_mem->cv_etamax == ONE)
{
// TODO(CJB): Should we consider using halpha here even though
// the error control heuristic wants us to defer step size or order changes?
SUNLogDebug(CV_LOGGER, __func__, "defer-step-order-change", "", "");
cv_mem->cv_qwait = SUNMAX(cv_mem->cv_qwait, 2);
cv_mem->cv_qprime = cv_mem->cv_q;
cv_mem->cv_hprime = cv_mem->cv_h;
Expand Down Expand Up @@ -3790,6 +3813,7 @@ static void cvOkSetEta(CVodeMem cv_mem)
if ((cv_mem->cv_eta > cv_mem->cv_eta_min_fx) &&
(cv_mem->cv_eta < cv_mem->cv_eta_max_fx))
{
SUNLogDebug(CV_LOGGER, __func__, "retain-step-size", "eta = %.16g, eta_min_fx = %.16g, eta_max_fx = %.16g", cv_mem->cv_eta, cv_mem->cv_eta_min_fx, cv_mem->cv_eta_max_fx);
/* Eta is within the fixed step bounds, retain step size */
cv_mem->cv_eta = ONE;
cv_mem->cv_hprime = cv_mem->cv_h;
Expand All @@ -3812,6 +3836,8 @@ static void cvOkSetEta(CVodeMem cv_mem)
}
/* Set hprime */
cv_mem->cv_hprime = cv_mem->cv_h * cv_mem->cv_eta;

/* Consider using halpha from the NLS switching algorithm instead of hprime */
cv_mem->cv_halpha = (cv_mem->cv_alpharef / cv_mem->cv_crate) * cv_mem->cv_h;
int prev_nls_iters = 0;
SUNNonlinSolGetCurIter(cv_mem->NLS, &prev_nls_iters);
Expand All @@ -3822,6 +3848,7 @@ static void cvOkSetEta(CVodeMem cv_mem)
cv_mem->cv_hprime = SUNMIN(cv_mem->cv_hprime, cv_mem->cv_halpha);
cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h;
}

if (cv_mem->cv_qprime < cv_mem->cv_q) { cv_mem->cv_nscon = 0; }
}
}
Expand Down
6 changes: 6 additions & 0 deletions src/cvode/cvode_impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -369,6 +369,7 @@ typedef struct CVodeMemRec
SUNNonlinearSolver NLS_newton;
SUNNonlinearSolver NLS_fixedpoint;
int NLS_algorithm;
int NLS_aavectors;
sunbooleantype ownNLS; /* flag indicating NLS ownership */
CVRhsFn nls_f; /* f(t,y(t)) used in the nonlinear solver */
int convfail; /* flag to indicate when a Jacobian update may
Expand Down Expand Up @@ -478,12 +479,17 @@ typedef struct CVodeMemRec
/* -----------------------
Nonlinear solver switching
------------------------ */
sunbooleantype cv_force_next_step;
int cv_lsodkr_strategy;
int cv_gustafsoder_strategy;
int cv_nst_switchtofixed_met;
int cv_nst_switchtonewton_met;
int switchtofixed_min_met;
int switchtonewton_min_met;
int cv_nst_switchtonewton_considered;
int cv_nst_switchtofixed_considered;
int switchtofixed_delay;
int switchtonewton_delay;
sunrealtype cv_halpha; /* step size dictated by nonlinear solver convergence control */
sunrealtype cv_stiff; /* stiffness estimate for nonlinear solver switching
this is ||res||/||del|| (beta_k in Gustafsson & Soderlind paper )*/
Expand Down
8 changes: 6 additions & 2 deletions src/cvode/cvode_io.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -920,7 +920,7 @@ int CVodeSetNoInactiveRootWarn(void* cvode_mem)
return (CV_SUCCESS);
}

int CVodeSetNonlinearSolverAlgorithm(void* cvode_mem, int algorithm)
int CVodeSetNonlinearSolverAlgorithm(void* cvode_mem, int algorithm, int m)
{
CVodeMem cv_mem;

Expand All @@ -933,6 +933,7 @@ int CVodeSetNonlinearSolverAlgorithm(void* cvode_mem, int algorithm)
cv_mem = (CVodeMem)cvode_mem;

cv_mem->NLS_algorithm = algorithm;
cv_mem->NLS_aavectors = m;

return CV_SUCCESS;
}
Expand Down Expand Up @@ -1645,6 +1646,7 @@ int CVodePrintAllStats(void* cvode_mem, FILE* outfile, SUNOutputFormat fmt)

cv_mem = (CVodeMem)cvode_mem;

long int total_rhs_evals = cv_mem->cv_nfe;
switch (fmt)
{
case SUN_OUTPUTFORMAT_TABLE:
Expand All @@ -1663,7 +1665,7 @@ int CVodePrintAllStats(void* cvode_mem, FILE* outfile, SUNOutputFormat fmt)
fprintf(outfile, "Stab. lim. order reductions = %ld\n", cv_mem->cv_nor);

/* function evaluations */
fprintf(outfile, "RHS fn evals = %ld\n", cv_mem->cv_nfe);
fprintf(outfile, "Integrator RHS fn evals = %ld\n", cv_mem->cv_nfe);

/* nonlinear solver stats */
fprintf(outfile, "NLS iters = %ld\n", cv_mem->cv_nni);
Expand Down Expand Up @@ -1699,6 +1701,7 @@ int CVodePrintAllStats(void* cvode_mem, FILE* outfile, SUNOutputFormat fmt)
fprintf(outfile, "Prec evals per NLS iter = %" RSYM "\n",
(sunrealtype)cvls_mem->npe / (sunrealtype)cv_mem->cv_nni);
}
total_rhs_evals += cvls_mem->nfeDQ;
}

/* rootfinding stats */
Expand All @@ -1712,6 +1715,7 @@ int CVodePrintAllStats(void* cvode_mem, FILE* outfile, SUNOutputFormat fmt)
fprintf(outfile, "Projection fails = %ld\n",
cvproj_mem->npfails);
}
fprintf(outfile, "Total RHS fn evals = %ld\n", total_rhs_evals);
break;

case SUN_OUTPUTFORMAT_CSV:
Expand Down
9 changes: 8 additions & 1 deletion src/sundials/sundials_logger_impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
#define SUNDIALS_LOGGING_INFO 3
#define SUNDIALS_LOGGING_DEBUG 4
#if SUNDIALS_LOGGING_LEVEL > SUNDIALS_LOGGING_DEBUG
#define SUNDIALS_LOGGING_EXTRA_DEBUG
#define SUNDIALS_LOGGING_EXTRA_DEBUG 5
#endif

#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_ERROR
Expand Down Expand Up @@ -60,6 +60,13 @@
#define SUNLogDebug(logger, scope, label, msg_txt, ...)
#endif

#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_EXTRA_DEBUG
#define SUNLogExtraDebug(logger, scope, label, msg_txt, ...) \
SUNLogger_QueueMsg(logger, SUN_LOGLEVEL_DEBUG, scope, label, msg_txt, __VA_ARGS__)
#else
#define SUNLogExtraDebug(logger, scope, label, msg_txt, ...)
#endif

struct SUNLogger_
{
/* MPI information */
Expand Down

0 comments on commit 3d06b93

Please sign in to comment.