From 14ecd14b10aa8ad724b77429d26f3e1c074a5606 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 29 Jul 2024 20:02:28 -0700 Subject: [PATCH 01/32] Remove special case of negative base in SUNRpowerR --- src/sundials/sundials_math.c | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/sundials/sundials_math.c b/src/sundials/sundials_math.c index b60e9c2540..7b05330d1c 100644 --- a/src/sundials/sundials_math.c +++ b/src/sundials/sundials_math.c @@ -22,6 +22,8 @@ #include #include +#include + sunrealtype SUNRpowerI(sunrealtype base, int exponent) { int i, expt; @@ -36,7 +38,9 @@ sunrealtype SUNRpowerI(sunrealtype base, int exponent) sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent) { - if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); } + // TODO(SBR): cleanup this and header + // if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); } + SUNAssert(base >= 0.0, "Base should be positive"); #if defined(SUNDIALS_DOUBLE_PRECISION) return (pow(base, exponent)); From dcb14bc1340af8c204f1571c55b9625786c4d574 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 29 Jul 2024 20:49:36 -0700 Subject: [PATCH 02/32] Fix assert condition --- src/sundials/sundials_math.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sundials/sundials_math.c b/src/sundials/sundials_math.c index 7b05330d1c..8e6b303290 100644 --- a/src/sundials/sundials_math.c +++ b/src/sundials/sundials_math.c @@ -40,7 +40,7 @@ sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent) { // TODO(SBR): cleanup this and header // if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); } - SUNAssert(base >= 0.0, "Base should be positive"); + SUNAssert(base > 0.0, "Base should be positive"); #if defined(SUNDIALS_DOUBLE_PRECISION) return (pow(base, exponent)); From 19f2dc957e8341c45265a7806bddc88fe2d04e91 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 29 Jul 2024 20:52:02 -0700 Subject: [PATCH 03/32] Remove magic constant --- src/arkode/arkode_adapt.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/arkode/arkode_adapt.c b/src/arkode/arkode_adapt.c index 7fef0ec2ae..6bd9aa18a1 100644 --- a/src/arkode/arkode_adapt.c +++ b/src/arkode/arkode_adapt.c @@ -131,7 +131,7 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur, "Error in explicit stability function."); return (ARK_ILL_INPUT); } - if (h_cfl <= ZERO) { h_cfl = SUN_RCONST(1.0e30) * SUNRabs(hcur); } + if (h_cfl <= ZERO) { h_cfl = INFINITY; } #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_INFO, "ARKODE::arkAdapt", From 53f3dde8403b4e4f44fac29143a707472bc04cae Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 29 Jul 2024 21:33:31 -0700 Subject: [PATCH 04/32] Remove TINY parameter --- src/arkode/arkode_adapt.c | 2 +- .../soderlind/sunadaptcontroller_soderlind.c | 11 +++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/arkode/arkode_adapt.c b/src/arkode/arkode_adapt.c index 6bd9aa18a1..a4758436eb 100644 --- a/src/arkode/arkode_adapt.c +++ b/src/arkode/arkode_adapt.c @@ -156,7 +156,7 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur, #endif /* increment the relevant step counter, set desired step */ - if (SUNRabs(h_acc) < SUNRabs(h_cfl)) { hadapt_mem->nst_acc++; } + if (SUNRabs(h_acc) <= SUNRabs(h_cfl)) { hadapt_mem->nst_acc++; } else { hadapt_mem->nst_exp++; } h_acc = int_dir * SUNMIN(SUNRabs(h_acc), SUNRabs(h_cfl)); diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index bca88252cb..b2924f12f0 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -62,7 +62,6 @@ #define DEFAULT_IMPGUS_K1 SUN_RCONST(0.98) /* Implicit Gustafsson parameters */ #define DEFAULT_IMPGUS_K2 SUN_RCONST(0.95) #define DEFAULT_BIAS SUN_RCONST(1.5) -#define TINY SUN_RCONST(1.0e-10) /* ----------------------------------------------------------------- * exported functions @@ -322,9 +321,9 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, const sunrealtype k3 = -SODERLIND_K3(C) / ord; const sunrealtype k4 = SODERLIND_K4(C); const sunrealtype k5 = SODERLIND_K5(C); - const sunrealtype e1 = SUNMAX(SODERLIND_BIAS(C) * dsm, TINY); - const sunrealtype e2 = SUNMAX(SODERLIND_EP(C), TINY); - const sunrealtype e3 = SUNMAX(SODERLIND_EPP(C), TINY); + const sunrealtype e1 = SODERLIND_BIAS(C) * dsm; + const sunrealtype e2 = SODERLIND_EP(C); + const sunrealtype e3 = SODERLIND_EPP(C); const sunrealtype hrat = h / SODERLIND_HP(C); const sunrealtype hrat2 = SODERLIND_HP(C) / SODERLIND_HPP(C); @@ -340,6 +339,10 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5); } + if (isnan(*hnew)) { + *hnew = INFINITY; + } + /* return with success */ return SUN_SUCCESS; } From e838c4ee2d6d4e76de79ddb8f7508620de0a3a06 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 29 Jul 2024 21:37:41 -0700 Subject: [PATCH 05/32] Remove TINY from imexgus controller --- .../imexgus/sunadaptcontroller_imexgus.c | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c index e4bb23e6dd..a4d28779e8 100644 --- a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c +++ b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c @@ -48,7 +48,6 @@ #define DEFAULT_K1I SUN_RCONST(0.95) #define DEFAULT_K2I SUN_RCONST(0.95) #define DEFAULT_BIAS SUN_RCONST(1.5) -#define TINY SUN_RCONST(1.0e-10) /* ----------------------------------------------------------------- * exported functions @@ -135,15 +134,15 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, /* set usable time-step adaptivity parameters -- first step */ const sunrealtype k = -SUN_RCONST(1.0) / ord; - const sunrealtype e = SUNMAX(SACIMEXGUS_BIAS(C) * dsm, TINY); + const sunrealtype e = SACIMEXGUS_BIAS(C) * dsm; /* set usable time-step adaptivity parameters -- subsequent steps */ const sunrealtype k1e = -SACIMEXGUS_K1E(C) / ord; const sunrealtype k2e = -SACIMEXGUS_K2E(C) / ord; const sunrealtype k1i = -SACIMEXGUS_K1I(C) / ord; const sunrealtype k2i = -SACIMEXGUS_K2I(C) / ord; - const sunrealtype e1 = SUNMAX(SACIMEXGUS_BIAS(C) * dsm, TINY); - const sunrealtype e2 = e1 / SUNMAX(SACIMEXGUS_EP(C), TINY); + const sunrealtype e1 = SACIMEXGUS_BIAS(C) * dsm; + const sunrealtype e2 = e1 / SACIMEXGUS_EP(C); const sunrealtype hrat = h / SACIMEXGUS_HP(C); /* compute estimated time step size, modifying the first step formula */ @@ -154,6 +153,10 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, SUNRpowerR(e1, k1e) * SUNRpowerR(e2, k2e)); } + if (isnan(*hnew)) { + *hnew = INFINITY; + } + /* return with success */ return SUN_SUCCESS; } From 9cde50af083b8d9de7d6c73c88dc6c5074e08e85 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 29 Jul 2024 21:44:20 -0700 Subject: [PATCH 06/32] Switch to assert --- src/sundials/sundials_math.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sundials/sundials_math.c b/src/sundials/sundials_math.c index 8e6b303290..8c5f5fc942 100644 --- a/src/sundials/sundials_math.c +++ b/src/sundials/sundials_math.c @@ -22,7 +22,7 @@ #include #include -#include +#include sunrealtype SUNRpowerI(sunrealtype base, int exponent) { @@ -40,7 +40,7 @@ sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent) { // TODO(SBR): cleanup this and header // if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); } - SUNAssert(base > 0.0, "Base should be positive"); + assert(base > 0.0); #if defined(SUNDIALS_DOUBLE_PRECISION) return (pow(base, exponent)); From 2c2444c9f4012790cf0abe29ed70812eb22728fa Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 29 Jul 2024 22:10:04 -0700 Subject: [PATCH 07/32] Fix sign of controller in nan case --- .../imexgus/sunadaptcontroller_imexgus.c | 7 +++++-- .../soderlind/sunadaptcontroller_soderlind.c | 7 +++++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c index a4d28779e8..8f881b49dc 100644 --- a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c +++ b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c @@ -153,8 +153,11 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, SUNRpowerR(e1, k1e) * SUNRpowerR(e2, k2e)); } - if (isnan(*hnew)) { - *hnew = INFINITY; + if (isnan(*hnew)) + { + /* hnew can be NAN if multiple e's are 0. In that case, make hnew inf with + * the same sign as h */ + *hnew = h / SUN_RCONST(0.0); } /* return with success */ diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index b2924f12f0..0ffc9313b9 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -339,8 +339,11 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5); } - if (isnan(*hnew)) { - *hnew = INFINITY; + if (isnan(*hnew)) + { + /* hnew can be NAN if multiple e's are 0. In that case, make hnew inf with + * the same sign as h */ + *hnew = h / SUN_RCONST(0.0); } /* return with success */ From 9cb7bba27aebd116eafefa2296d5aab92d356f85 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 30 Jul 2024 21:49:21 -0700 Subject: [PATCH 08/32] Start adaptivity unit test --- .../unit_tests/arkode/C_serial/CMakeLists.txt | 1 + .../arkode/C_serial/ark_test_adapt.c | 122 ++++++++++++++++++ 2 files changed, 123 insertions(+) create mode 100644 test/unit_tests/arkode/C_serial/ark_test_adapt.c diff --git a/test/unit_tests/arkode/C_serial/CMakeLists.txt b/test/unit_tests/arkode/C_serial/CMakeLists.txt index 6794a213b4..3670b76da4 100644 --- a/test/unit_tests/arkode/C_serial/CMakeLists.txt +++ b/test/unit_tests/arkode/C_serial/CMakeLists.txt @@ -16,6 +16,7 @@ # List of test tuples of the form "name\;args" set(ARKODE_unit_tests + "ark_test_adapt\;" "ark_test_arkstepsetforcing\;1 0" "ark_test_arkstepsetforcing\;1 1" "ark_test_arkstepsetforcing\;1 2" diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c new file mode 100644 index 0000000000..4e98e2fe4b --- /dev/null +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -0,0 +1,122 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): Steven B. Roberts @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2024, Lawrence Livermore National Security + * and Southern Methodist University. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------------------- + * Unit test for GetUserData functions + * ---------------------------------------------------------------------------*/ + +#include +#include + +#include "arkode/arkode_erkstep.h" +#include "nvector/nvector_serial.h" + +static void err_fn(int line, const char* func, const char* file, const char* msg, + SUNErrCode err_code, void* err_user_data, SUNContext sunctx) +{ + fprintf(stderr, "Error at line %i of %s in %s: %s\n", line, func, file, msg); + exit(err_code); +} + +static int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +{ + N_VConst(SUN_RCONST(0.0), ydot); + return 0; +} + +static int check_step(void* arkode_mem, N_Vector y, sunrealtype h_expected) +{ + static int step = 0; + step++; + + sunrealtype tret; + ARKodeEvolve(arkode_mem, SUN_RCONST(1.0), y, &tret, ARK_ONE_STEP); + + long int local_err_fails; + ARKodeGetNumErrTestFails(arkode_mem, &local_err_fails); + if (local_err_fails != 0) + { + fprintf(stderr, "Expected 0 local error failures at step %i but is %li\n", + step, local_err_fails); + } + + const N_Vector err = N_VClone(y); + ARKodeGetEstLocalErrors(arkode_mem, err); + sunrealtype err_norm = N_VMaxNorm(err); + N_VDestroy(err); + + if (err_norm != 0) + { + fprintf(stderr, "Expected local error at step %i to be 0 but is %g\n", step, + err_norm); + return 1; + } + + sunrealtype h_actual; + ARKodeGetCurrentStep(arkode_mem, &h_actual); + if (h_expected != h_actual) + { + fprintf(stderr, "Expected h at step %i to be %g but is %g\n", step, + h_expected, h_actual); + return 1; + } + + return 0; +} + +int main() +{ + SUNContext sunctx; + int retval = SUNContext_Create(SUN_COMM_NULL, &sunctx); + if (retval != 0) + { + fprintf(stderr, "SUNContext_Create returned %i\n", retval); + return 1; + } + + retval = SUNContext_PushErrHandler(sunctx, err_fn, NULL); + if (retval != 0) + { + fprintf(stderr, "SUNContext_PushErrHandler returned %i\n", retval); + return 1; + } + + const N_Vector y = N_VNew_Serial(1, sunctx); + N_VConst(SUN_RCONST(1.0), y); + + void* arkode_mem = ERKStepCreate(f, SUN_RCONST(0.0), y, sunctx); + + const sunrealtype h0 = SUN_RCONST(1.0e-4); + const sunrealtype first_growth = SUN_RCONST(1234.0); + const sunrealtype growth = SUN_RCONST(3.0); + + ARKodeSetInitStep(arkode_mem, h0); + ARKodeSetMaxFirstGrowth(arkode_mem, first_growth); + ARKodeSetMaxGrowth(arkode_mem, growth); + + sunrealtype h_expect = first_growth * h0; + + for (int i = 0; i < 4; i++) + { + retval = check_step(arkode_mem, y, h_expect); + if (retval != 0) { return retval; } + h_expect *= growth; + } + + ARKodeFree(&arkode_mem); + N_VDestroy(y); + SUNContext_Free(&sunctx); + + printf("SUCCESS\n"); + + return 0; +} From 743913a76e0eec612bced620f12f44fcdf353b7b Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 31 Jul 2024 18:13:37 -0700 Subject: [PATCH 09/32] Add const and comments to test --- .../arkode/C_serial/ark_test_adapt.c | 27 +++++++++++-------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index 4e98e2fe4b..01fa551559 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -20,25 +20,28 @@ #include "arkode/arkode_erkstep.h" #include "nvector/nvector_serial.h" -static void err_fn(int line, const char* func, const char* file, const char* msg, - SUNErrCode err_code, void* err_user_data, SUNContext sunctx) +static void err_fn(const int line, const char* const func, const char* const file, + const char* const msg, const SUNErrCode err_code, + void* const err_user_data, const SUNContext sunctx) { fprintf(stderr, "Error at line %i of %s in %s: %s\n", line, func, file, msg); exit(err_code); } -static int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +// RHS for the simple ODE y' = 0 +static int f(const sunrealtype t, const N_Vector y, const N_Vector ydot, + void* const user_data) { N_VConst(SUN_RCONST(0.0), ydot); return 0; } -static int check_step(void* arkode_mem, N_Vector y, sunrealtype h_expected) +static int check_step(void* const arkode_mem, const N_Vector y, + const sunrealtype h_expected, const int step) { - static int step = 0; - step++; - sunrealtype tret; + /* The ERK method should be able to take the maximum possible timestep for the + simple ODE y' = 0 without any rejected steps or local error */ ARKodeEvolve(arkode_mem, SUN_RCONST(1.0), y, &tret, ARK_ONE_STEP); long int local_err_fails; @@ -51,7 +54,7 @@ static int check_step(void* arkode_mem, N_Vector y, sunrealtype h_expected) const N_Vector err = N_VClone(y); ARKodeGetEstLocalErrors(arkode_mem, err); - sunrealtype err_norm = N_VMaxNorm(err); + const sunrealtype err_norm = N_VMaxNorm(err); N_VDestroy(err); if (err_norm != 0) @@ -104,10 +107,12 @@ int main() ARKodeSetMaxGrowth(arkode_mem, growth); sunrealtype h_expect = first_growth * h0; - - for (int i = 0; i < 4; i++) + /* Take several steps to see the special behavior at step one then to allow + the adaptivity controller history fill up */ + const int num_steps = 5; + for (int step = 1; step <= num_steps; step++) { - retval = check_step(arkode_mem, y, h_expect); + retval = check_step(arkode_mem, y, h_expect, step); if (retval != 0) { return retval; } h_expect *= growth; } From b75e78e9af1c37b0e4fa1b83ae7daeac3f5ff71f Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 31 Jul 2024 19:26:29 -0700 Subject: [PATCH 10/32] Simplify Soderlind controller and handle inf/nan better --- .../soderlind/sunadaptcontroller_soderlind.c | 41 ++++++++----------- 1 file changed, 17 insertions(+), 24 deletions(-) diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index 0ffc9313b9..fef0d82b6e 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -316,33 +316,26 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, const int ord = p + 1; /* set usable time-step adaptivity parameters */ - const sunrealtype k1 = -SODERLIND_K1(C) / ord; - const sunrealtype k2 = -SODERLIND_K2(C) / ord; - const sunrealtype k3 = -SODERLIND_K3(C) / ord; - const sunrealtype k4 = SODERLIND_K4(C); - const sunrealtype k5 = SODERLIND_K5(C); - const sunrealtype e1 = SODERLIND_BIAS(C) * dsm; - const sunrealtype e2 = SODERLIND_EP(C); - const sunrealtype e3 = SODERLIND_EPP(C); - const sunrealtype hrat = h / SODERLIND_HP(C); - const sunrealtype hrat2 = SODERLIND_HP(C) / SODERLIND_HPP(C); - - /* compute estimated optimal time step size */ - if (SODERLIND_FIRSTSTEPS(C) < 1) { *hnew = h * SUNRpowerR(e1, k1); } - else if (SODERLIND_FIRSTSTEPS(C) < 2) - { - *hnew = h * SUNRpowerR(e1, k1) * SUNRpowerR(e2, k2) * SUNRpowerR(hrat, k4); - } - else - { - *hnew = h * SUNRpowerR(e1, k1) * SUNRpowerR(e2, k2) * SUNRpowerR(e3, k3) * - SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5); + const sunrealtype ke[] = {-SODERLIND_K1(C) / ord, -SODERLIND_K2(C) / ord, -SODERLIND_K3(C) / ord}; + const sunrealtype kh[] = {SODERLIND_K4(C), SODERLIND_K5(C)}; + // TODO: Should e2 and e3 be set to e1 when we don't have enough history yet? + const sunrealtype e[] = {SODERLIND_BIAS(C) * dsm, SODERLIND_EP(C), SODERLIND_EPP(C)}; + const sunrealtype hrat[] = {h / SODERLIND_HP(C), SODERLIND_HP(C) / SODERLIND_HPP(C)}; + + // This assumes ke[0] is nonzero. Otherwise we could have pow(0, 0) which is undefined behavior in C + *hnew = h * SUNRpowerR(e[0], ke[0]); + for (int i = 0; i < SUNMIN(2, SODERLIND_FIRSTSTEPS(C)); i++) { + *hnew *= SUNRpowerR(hrat[i], kh[i]); + if (ke[i + 1] != 0) { + // This check avoids pow(0, 0) + *hnew *= SUNRpowerR(e[i + 1], ke[i + 1]); + } } - if (isnan(*hnew)) + if (!isfinite(*hnew)) { - /* hnew can be NAN if multiple e's are 0. In that case, make hnew inf with - * the same sign as h */ + /* hnew can be INFINITY or NAN if multiple e's are 0 or there are overflows. + * In that case, make hnew INFINITY with the same sign as h */ *hnew = h / SUN_RCONST(0.0); } From 2caec18128a0e2f92187b2de2b91061464d8fe0d Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 31 Jul 2024 20:11:37 -0700 Subject: [PATCH 11/32] Clean up controllers --- .../imexgus/sunadaptcontroller_imexgus.c | 35 ++++++++++--------- .../soderlind/sunadaptcontroller_soderlind.c | 7 +--- 2 files changed, 19 insertions(+), 23 deletions(-) diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c index 8f881b49dc..cfb73cb5a0 100644 --- a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c +++ b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c @@ -132,31 +132,32 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, /* order parameter to use */ const int ord = p + 1; - /* set usable time-step adaptivity parameters -- first step */ - const sunrealtype k = -SUN_RCONST(1.0) / ord; - const sunrealtype e = SACIMEXGUS_BIAS(C) * dsm; - - /* set usable time-step adaptivity parameters -- subsequent steps */ - const sunrealtype k1e = -SACIMEXGUS_K1E(C) / ord; - const sunrealtype k2e = -SACIMEXGUS_K2E(C) / ord; - const sunrealtype k1i = -SACIMEXGUS_K1I(C) / ord; - const sunrealtype k2i = -SACIMEXGUS_K2I(C) / ord; - const sunrealtype e1 = SACIMEXGUS_BIAS(C) * dsm; - const sunrealtype e2 = e1 / SACIMEXGUS_EP(C); - const sunrealtype hrat = h / SACIMEXGUS_HP(C); - /* compute estimated time step size, modifying the first step formula */ - if (SACIMEXGUS_FIRSTSTEP(C)) { *hnew = h * SUNRpowerR(e, k); } + if (SACIMEXGUS_FIRSTSTEP(C)) { + // TODO: could this case be handled by setting e1 = e2? + /* set usable time-step adaptivity parameters -- first step */ + const sunrealtype k = -SUN_RCONST(1.0) / ord; + const sunrealtype e = SACIMEXGUS_BIAS(C) * dsm; + *hnew = h * SUNRpowerR(e, k); + } else { + /* set usable time-step adaptivity parameters -- subsequent steps */ + const sunrealtype k1e = -SACIMEXGUS_K1E(C) / ord; + const sunrealtype k2e = -SACIMEXGUS_K2E(C) / ord; + const sunrealtype k1i = -SACIMEXGUS_K1I(C) / ord; + const sunrealtype k2i = -SACIMEXGUS_K2I(C) / ord; + const sunrealtype e1 = SACIMEXGUS_BIAS(C) * dsm; + const sunrealtype e2 = e1 / SACIMEXGUS_EP(C); + const sunrealtype hrat = h / SACIMEXGUS_HP(C); *hnew = h * SUNMIN(hrat * SUNRpowerR(e1, k1i) * SUNRpowerR(e2, k2i), SUNRpowerR(e1, k1e) * SUNRpowerR(e2, k2e)); } - if (isnan(*hnew)) + if (!isfinite(*hnew)) { - /* hnew can be NAN if multiple e's are 0. In that case, make hnew inf with - * the same sign as h */ + /* hnew can be INFINITY or NAN if multiple e's are 0 or there are overflows. + * In that case, make hnew INFINITY with the same sign as h */ *hnew = h / SUN_RCONST(0.0); } diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index fef0d82b6e..b4286ff011 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -322,14 +322,9 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, const sunrealtype e[] = {SODERLIND_BIAS(C) * dsm, SODERLIND_EP(C), SODERLIND_EPP(C)}; const sunrealtype hrat[] = {h / SODERLIND_HP(C), SODERLIND_HP(C) / SODERLIND_HPP(C)}; - // This assumes ke[0] is nonzero. Otherwise we could have pow(0, 0) which is undefined behavior in C *hnew = h * SUNRpowerR(e[0], ke[0]); for (int i = 0; i < SUNMIN(2, SODERLIND_FIRSTSTEPS(C)); i++) { - *hnew *= SUNRpowerR(hrat[i], kh[i]); - if (ke[i + 1] != 0) { - // This check avoids pow(0, 0) - *hnew *= SUNRpowerR(e[i + 1], ke[i + 1]); - } + *hnew *= SUNRpowerR(hrat[i], kh[i]) * SUNRpowerR(e[i + 1], ke[i + 1]); } if (!isfinite(*hnew)) From 3fbda5b1be8f202e64bc5d98ca439030e72b14c3 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 31 Jul 2024 20:16:36 -0700 Subject: [PATCH 12/32] Apply formatter --- .../imexgus/sunadaptcontroller_imexgus.c | 5 +++-- .../soderlind/sunadaptcontroller_soderlind.c | 12 ++++++++---- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c index cfb73cb5a0..fbbea251dc 100644 --- a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c +++ b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c @@ -133,12 +133,13 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, const int ord = p + 1; /* compute estimated time step size, modifying the first step formula */ - if (SACIMEXGUS_FIRSTSTEP(C)) { + if (SACIMEXGUS_FIRSTSTEP(C)) + { // TODO: could this case be handled by setting e1 = e2? /* set usable time-step adaptivity parameters -- first step */ const sunrealtype k = -SUN_RCONST(1.0) / ord; const sunrealtype e = SACIMEXGUS_BIAS(C) * dsm; - *hnew = h * SUNRpowerR(e, k); + *hnew = h * SUNRpowerR(e, k); } else { diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index b4286ff011..62c0e823b0 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -316,14 +316,18 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, const int ord = p + 1; /* set usable time-step adaptivity parameters */ - const sunrealtype ke[] = {-SODERLIND_K1(C) / ord, -SODERLIND_K2(C) / ord, -SODERLIND_K3(C) / ord}; + const sunrealtype ke[] = {-SODERLIND_K1(C) / ord, -SODERLIND_K2(C) / ord, + -SODERLIND_K3(C) / ord}; const sunrealtype kh[] = {SODERLIND_K4(C), SODERLIND_K5(C)}; // TODO: Should e2 and e3 be set to e1 when we don't have enough history yet? - const sunrealtype e[] = {SODERLIND_BIAS(C) * dsm, SODERLIND_EP(C), SODERLIND_EPP(C)}; - const sunrealtype hrat[] = {h / SODERLIND_HP(C), SODERLIND_HP(C) / SODERLIND_HPP(C)}; + const sunrealtype e[] = {SODERLIND_BIAS(C) * dsm, SODERLIND_EP(C), + SODERLIND_EPP(C)}; + const sunrealtype hrat[] = {h / SODERLIND_HP(C), + SODERLIND_HP(C) / SODERLIND_HPP(C)}; *hnew = h * SUNRpowerR(e[0], ke[0]); - for (int i = 0; i < SUNMIN(2, SODERLIND_FIRSTSTEPS(C)); i++) { + for (int i = 0; i < SUNMIN(2, SODERLIND_FIRSTSTEPS(C)); i++) + { *hnew *= SUNRpowerR(hrat[i], kh[i]) * SUNRpowerR(e[i + 1], ke[i + 1]); } From 299cd5cf89f80236a9b1527eb559c05b88246389 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 5 Sep 2024 21:44:42 -0700 Subject: [PATCH 13/32] Variable initialization cleanup --- .../soderlind/sunadaptcontroller_soderlind.c | 26 +++++-------------- 1 file changed, 7 insertions(+), 19 deletions(-) diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index 62c0e823b0..4b306b972c 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -75,12 +75,8 @@ SUNAdaptController SUNAdaptController_Soderlind(SUNContext sunctx) { SUNFunctionBegin(sunctx); - SUNAdaptController C; - SUNAdaptControllerContent_Soderlind content; - /* Create an empty controller object */ - C = NULL; - C = SUNAdaptController_NewEmpty(sunctx); + SUNAdaptController C = SUNAdaptController_NewEmpty(sunctx); SUNCheckLastErrNull(); /* Attach operations */ @@ -94,12 +90,8 @@ SUNAdaptController SUNAdaptController_Soderlind(SUNContext sunctx) C->ops->space = SUNAdaptController_Space_Soderlind; /* Create content */ - content = NULL; - content = (SUNAdaptControllerContent_Soderlind)malloc(sizeof *content); - SUNAssertNull(content, SUN_ERR_MALLOC_FAIL); - - /* Attach content */ - C->content = content; + C->content = (SUNAdaptControllerContent_Soderlind)malloc(sizeof *C->content); + SUNAssertNull(C->content, SUN_ERR_MALLOC_FAIL); /* Fill content with default/reset values */ SUNCheckCallNull(SUNAdaptController_SetDefaults_Soderlind(C)); @@ -168,8 +160,7 @@ SUNAdaptController SUNAdaptController_PI(SUNContext sunctx) { SUNFunctionBegin(sunctx); - SUNAdaptController C; - C = SUNAdaptController_Soderlind(sunctx); + SUNAdaptController C = SUNAdaptController_Soderlind(sunctx); SUNCheckLastErrNull(); SUNCheckCallNull( @@ -202,8 +193,7 @@ SUNAdaptController SUNAdaptController_I(SUNContext sunctx) { SUNFunctionBegin(sunctx); - SUNAdaptController C; - C = SUNAdaptController_Soderlind(sunctx); + SUNAdaptController C = SUNAdaptController_Soderlind(sunctx); SUNCheckLastErrNull(); SUNCheckCallNull(SUNAdaptController_SetParams_I(C, DEFAULT_I_K1)); @@ -234,8 +224,7 @@ SUNAdaptController SUNAdaptController_ExpGus(SUNContext sunctx) { SUNFunctionBegin(sunctx); - SUNAdaptController C; - C = SUNAdaptController_Soderlind(sunctx); + SUNAdaptController C = SUNAdaptController_Soderlind(sunctx); SUNCheckLastErrNull(); SUNCheckCallNull(SUNAdaptController_SetParams_ExpGus(C, DEFAULT_EXPGUS_K1, @@ -268,8 +257,7 @@ SUNAdaptController SUNAdaptController_ImpGus(SUNContext sunctx) { SUNFunctionBegin(sunctx); - SUNAdaptController C; - C = SUNAdaptController_Soderlind(sunctx); + SUNAdaptController C = SUNAdaptController_Soderlind(sunctx); SUNCheckLastErrNull(); SUNCheckCallNull(SUNAdaptController_SetParams_ImpGus(C, DEFAULT_IMPGUS_K1, From 405ff411ecbd6ecbce9b35a69d087f1d1d78a39d Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 5 Sep 2024 22:50:49 -0700 Subject: [PATCH 14/32] Use I controller on first 2 steps of Soderlind --- .../imexgus/sunadaptcontroller_imexgus.c | 1 - .../soderlind/sunadaptcontroller_soderlind.c | 34 +++++++++++-------- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c index fbbea251dc..5880f24a31 100644 --- a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c +++ b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c @@ -135,7 +135,6 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, /* compute estimated time step size, modifying the first step formula */ if (SACIMEXGUS_FIRSTSTEP(C)) { - // TODO: could this case be handled by setting e1 = e2? /* set usable time-step adaptivity parameters -- first step */ const sunrealtype k = -SUN_RCONST(1.0) / ord; const sunrealtype e = SACIMEXGUS_BIAS(C) * dsm; diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index 4b306b972c..1545cbab4e 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -303,20 +303,26 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, /* order parameter to use */ const int ord = p + 1; - /* set usable time-step adaptivity parameters */ - const sunrealtype ke[] = {-SODERLIND_K1(C) / ord, -SODERLIND_K2(C) / ord, - -SODERLIND_K3(C) / ord}; - const sunrealtype kh[] = {SODERLIND_K4(C), SODERLIND_K5(C)}; - // TODO: Should e2 and e3 be set to e1 when we don't have enough history yet? - const sunrealtype e[] = {SODERLIND_BIAS(C) * dsm, SODERLIND_EP(C), - SODERLIND_EPP(C)}; - const sunrealtype hrat[] = {h / SODERLIND_HP(C), - SODERLIND_HP(C) / SODERLIND_HPP(C)}; - - *hnew = h * SUNRpowerR(e[0], ke[0]); - for (int i = 0; i < SUNMIN(2, SODERLIND_FIRSTSTEPS(C)); i++) - { - *hnew *= SUNRpowerR(hrat[i], kh[i]) * SUNRpowerR(e[i + 1], ke[i + 1]); + if (SODERLIND_FIRSTSTEPS(C) > 1) { + /* After the first 2 steps, there is sufficient history */ + const sunrealtype k1 = -SODERLIND_K1(C) / ord; + const sunrealtype k2 = -SODERLIND_K2(C) / ord; + const sunrealtype k3 = -SODERLIND_K3(C) / ord; + const sunrealtype k4 = SODERLIND_K4(C); + const sunrealtype k5 = SODERLIND_K5(C); + const sunrealtype e1 = SODERLIND_BIAS(C) * dsm; + const sunrealtype e2 = SODERLIND_EP(C); + const sunrealtype e3 = SODERLIND_EPP(C); + const sunrealtype hrat = h / SODERLIND_HP(C); + const sunrealtype hrat2 = SODERLIND_HP(C) / SODERLIND_HPP(C); + + *hnew = h * SUNRpowerR(e1, k1) * SUNRpowerR(e2, k2) * SUNRpowerR(e3, k3) * + SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5); + } else { + /* Use an I controller on the first two steps */ + const sunrealtype k = -SUN_RCONST(1.0) / ord; + const sunrealtype e = SODERLIND_BIAS(C) * dsm; + *hnew = h * SUNRpowerR(e, k); } if (!isfinite(*hnew)) From bc9f312a12ca9c09b9fe00c7268116c410a4af5d Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 5 Sep 2024 22:54:04 -0700 Subject: [PATCH 15/32] Remove assert for pow testing --- src/sundials/sundials_math.c | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/sundials/sundials_math.c b/src/sundials/sundials_math.c index 8c5f5fc942..c78bb59068 100644 --- a/src/sundials/sundials_math.c +++ b/src/sundials/sundials_math.c @@ -38,10 +38,6 @@ sunrealtype SUNRpowerI(sunrealtype base, int exponent) sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent) { - // TODO(SBR): cleanup this and header - // if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); } - assert(base > 0.0); - #if defined(SUNDIALS_DOUBLE_PRECISION) return (pow(base, exponent)); #elif defined(SUNDIALS_SINGLE_PRECISION) From 54aa1d5a859bfdfa168a477ce31d9ae57f39f45a Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 5 Sep 2024 23:16:13 -0700 Subject: [PATCH 16/32] Update docs --- .../sunadaptcontroller/SUNAdaptController_Soderlind.rst | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/doc/shared/sunadaptcontroller/SUNAdaptController_Soderlind.rst b/doc/shared/sunadaptcontroller/SUNAdaptController_Soderlind.rst index 0c1b0cc620..830a7f2d17 100644 --- a/doc/shared/sunadaptcontroller/SUNAdaptController_Soderlind.rst +++ b/doc/shared/sunadaptcontroller/SUNAdaptController_Soderlind.rst @@ -27,11 +27,10 @@ and :cite:p:`Sod:06`. This controller has the form with default parameter values :math:`k_1 = 1.25`, :math:`k_2 = 0.5`, :math:`k_3 = -0.75`, :math:`k_4 = 0.25`, and :math:`k_5 = 0.75`, where -:math:`p` is the global order of the time integration method. In this estimate, -a floor of :math:`\varepsilon_* > 10^{-10}` is enforced to avoid division-by-zero -errors. During the first two steps (when :math:`\varepsilon_{n-2}`, -:math:`\varepsilon_{n-1}`, :math:`h_{n-2}`, and :math:`h_{n-2}` may be unavailable), -the corresponding terms are merely omitted during estimation of :math:`h'`. +:math:`p` is the global order of the time integration method. During the first +two steps (when :math:`\varepsilon_{n-2}`, :math:`\varepsilon_{n-1}`, +:math:`h_{n-2}`, and :math:`h_{n-2}` may be unavailable), the I controller +:math:`h' = h_n \varepsilon_1^{-1/(p+1)}` is used. The SUNAdaptController_Soderlind controller is implemented as a derived SUNAdaptController class, and defines its *content* field as: From bfdf75e0e56e60cc6721b237c7780e275c53c872 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 5 Sep 2024 23:35:46 -0700 Subject: [PATCH 17/32] Update IMEX controller docs --- doc/shared/sunadaptcontroller/SUNAdaptController_ImExGus.rst | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/doc/shared/sunadaptcontroller/SUNAdaptController_ImExGus.rst b/doc/shared/sunadaptcontroller/SUNAdaptController_ImExGus.rst index 0afdf7d786..11284031c3 100644 --- a/doc/shared/sunadaptcontroller/SUNAdaptController_ImExGus.rst +++ b/doc/shared/sunadaptcontroller/SUNAdaptController_ImExGus.rst @@ -45,9 +45,7 @@ with the form In the above formulas, the default values of :math:`k_1^E`, :math:`k_2^E`, :math:`k_1^I`, and :math:`k_2^I` are 0.367, 0.268, 0.98, and 0.95, respectively, -and :math:`p` is the global order of the time integration method. In these -estimates, a floor of :math:`\varepsilon_* > 10^{-10}` is enforced to avoid -division-by-zero errors. +and :math:`p` is the global order of the time integration method. The SUNAdaptController_ImExGus controller implements both formulas :eq:`expGusController` and :eq:`impGusController`, and sets its recommended step From dde1379cdac876a28d4ae13b880e6c0f7a8ab613 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 5 Sep 2024 23:35:54 -0700 Subject: [PATCH 18/32] Update changelog --- CHANGELOG.md | 6 ++++++ doc/shared/RecentChanges.rst | 7 +++++++ 2 files changed, 13 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3a82dd8621..641d5e6a2f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,12 @@ architecture for their system. ### Bug Fixes +Removed error floors from the SUNAdaptController implementations which could +unnecessarily limit the time size growth, particularly after the first step. + +On the first two time steps, the Soderlind controller uses an I controller +instead omitting unavailable terms. + Fixed the loading of ARKStep's default first order explicit method. Fixed a CMake bug regarding usage of missing "print_warning" macro diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index 0cd3986ee3..d77b9d1e20 100644 --- a/doc/shared/RecentChanges.rst +++ b/doc/shared/RecentChanges.rst @@ -10,6 +10,13 @@ override this value with the architecture for their system. **Bug Fixes** +Removed error floors from the SUNAdaptController implementations which could +unnecessarily limit the time size growth, particularly after the first step. + +On the first two time steps, the +:ref:`Soderlind controller ` uses an I controller +instead omitting unavailable terms. + Fixed the loading of ARKStep's default first order explicit method. Fixed a CMake bug regarding usage of missing "print_warning" macro From facfc39b935f4ac4753e01607c811c4bb9523fc2 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 6 Sep 2024 02:40:26 -0700 Subject: [PATCH 19/32] Clean up time step sign --- include/sundials/sundials_math.h | 24 +++++++++++++++++++ src/arkode/arkode_adapt.c | 18 +++++++------- .../imexgus/sunadaptcontroller_imexgus.c | 2 +- .../soderlind/sunadaptcontroller_soderlind.c | 4 ++-- 4 files changed, 35 insertions(+), 13 deletions(-) diff --git a/include/sundials/sundials_math.h b/include/sundials/sundials_math.h index 4585c7c7f4..bc6cd2f9eb 100644 --- a/include/sundials/sundials_math.h +++ b/include/sundials/sundials_math.h @@ -158,6 +158,30 @@ extern "C" { #endif #endif +/* + * ----------------------------------------------------------------- + * Function : SUNRcopysign + * ----------------------------------------------------------------- + * Usage : sunrealtype z; + * z = SUNRcopysign(x, y); + * ----------------------------------------------------------------- + * SUNRcopysign(x, y) returns x with the sign of y. + * ----------------------------------------------------------------- + */ + +#ifndef SUNRcopysign +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define SUNRcopysign(x, y) (copysign((x), (y))) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define SUNRcopysign(x, y) (copysignf((x), (y))) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define SUNRcopysign(x, y) (copysignl((x), (y))) +#else +#error \ + "SUNDIALS precision not defined, report to github.com/LLNL/sundials/issues" +#endif +#endif + /* * ----------------------------------------------------------------- * Function : SUNRpowerI diff --git a/src/arkode/arkode_adapt.c b/src/arkode/arkode_adapt.c index a4758436eb..93a09bf683 100644 --- a/src/arkode/arkode_adapt.c +++ b/src/arkode/arkode_adapt.c @@ -95,7 +95,7 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur, sunrealtype tcur, sunrealtype hcur, sunrealtype dsm) { int retval; - sunrealtype h_acc, h_cfl, int_dir; + sunrealtype h_acc, h_cfl; int controller_order; /* Request error-based step size from adaptivity controller */ @@ -120,9 +120,6 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur, return (ARK_CONTROLLER_ERR); } - /* determine direction of integration */ - int_dir = hcur / SUNRabs(hcur); - /* Call explicit stability function */ retval = hadapt_mem->expstab(ycur, tcur, &h_cfl, hadapt_mem->estab_data); if (retval != ARK_SUCCESS) @@ -141,24 +138,25 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur, /* enforce safety factors */ h_acc *= hadapt_mem->safety; - h_cfl *= hadapt_mem->cfl * int_dir; + h_cfl *= hadapt_mem->cfl; /* enforce maximum bound on time step growth */ - h_acc = int_dir * SUNMIN(SUNRabs(h_acc), SUNRabs(hadapt_mem->etamax * hcur)); + h_acc = SUNMIN(SUNRabs(h_acc), SUNRabs(hadapt_mem->etamax * hcur)); /* enforce minimum bound time step reduction */ - h_acc = int_dir * SUNMAX(SUNRabs(h_acc), SUNRabs(hadapt_mem->etamin * hcur)); + h_acc = SUNMAX(h_acc, SUNRabs(hadapt_mem->etamin * hcur)); #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_INFO, "ARKODE::arkAdapt", "new-step-after-max-min-bounds", - "h_acc = %" RSYM ", h_cfl = %" RSYM, h_acc, h_cfl); + "h_acc = %" RSYM ", h_cfl = %" RSYM, + SUNRcopysign(h_acc, hcur), h_cfl); #endif /* increment the relevant step counter, set desired step */ - if (SUNRabs(h_acc) <= SUNRabs(h_cfl)) { hadapt_mem->nst_acc++; } + if (h_acc <= h_cfl) { hadapt_mem->nst_acc++; } else { hadapt_mem->nst_exp++; } - h_acc = int_dir * SUNMIN(SUNRabs(h_acc), SUNRabs(h_cfl)); + h_acc = SUNRcopysign(SUNMIN(h_acc, h_cfl), hcur); /* enforce adaptivity bounds to retain Jacobian/preconditioner accuracy */ if (dsm <= ONE) diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c index 5880f24a31..650ffc52df 100644 --- a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c +++ b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c @@ -158,7 +158,7 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, { /* hnew can be INFINITY or NAN if multiple e's are 0 or there are overflows. * In that case, make hnew INFINITY with the same sign as h */ - *hnew = h / SUN_RCONST(0.0); + *hnew = SUNRcopysign(INFINITY, h); } /* return with success */ diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index 1545cbab4e..2f0d0e71e1 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -90,7 +90,7 @@ SUNAdaptController SUNAdaptController_Soderlind(SUNContext sunctx) C->ops->space = SUNAdaptController_Space_Soderlind; /* Create content */ - C->content = (SUNAdaptControllerContent_Soderlind)malloc(sizeof *C->content); + C->content = (SUNAdaptControllerContent_Soderlind)malloc(sizeof(struct _SUNAdaptControllerContent_Soderlind)); SUNAssertNull(C->content, SUN_ERR_MALLOC_FAIL); /* Fill content with default/reset values */ @@ -329,7 +329,7 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, { /* hnew can be INFINITY or NAN if multiple e's are 0 or there are overflows. * In that case, make hnew INFINITY with the same sign as h */ - *hnew = h / SUN_RCONST(0.0); + *hnew = SUNRcopysign(INFINITY, h); } /* return with success */ From 7badd781110b5d644c731fecf61ead24ce124f3c Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 6 Sep 2024 03:01:05 -0700 Subject: [PATCH 20/32] Additional simplifications --- src/arkode/arkode_adapt.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/arkode/arkode_adapt.c b/src/arkode/arkode_adapt.c index 93a09bf683..9f91b684de 100644 --- a/src/arkode/arkode_adapt.c +++ b/src/arkode/arkode_adapt.c @@ -155,18 +155,18 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur, /* increment the relevant step counter, set desired step */ if (h_acc <= h_cfl) { hadapt_mem->nst_acc++; } - else { hadapt_mem->nst_exp++; } - h_acc = SUNRcopysign(SUNMIN(h_acc, h_cfl), hcur); + else { hadapt_mem->nst_exp++; h_acc = h_cfl; } /* enforce adaptivity bounds to retain Jacobian/preconditioner accuracy */ if (dsm <= ONE) { - if ((SUNRabs(h_acc) > SUNRabs(hcur * hadapt_mem->lbound * ONEMSM)) && - (SUNRabs(h_acc) < SUNRabs(hcur * hadapt_mem->ubound * ONEPSM))) + if ((h_acc > SUNRabs(hcur * hadapt_mem->lbound * ONEMSM)) && + (h_acc < SUNRabs(hcur * hadapt_mem->ubound * ONEPSM))) { h_acc = hcur; } } + h_acc = SUNRcopysign(h_acc, hcur); /* set basic value of ark_eta */ ark_mem->eta = h_acc / hcur; From c472ffed440afc414b4c1dcf5f4545b3139bc313 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 6 Sep 2024 14:14:51 -0700 Subject: [PATCH 21/32] Additional testing of controller --- src/arkode/arkode_adapt.c | 6 +- .../soderlind/sunadaptcontroller_soderlind.c | 10 ++- .../arkode/C_serial/ark_test_adapt.c | 75 +++++++++++++------ 3 files changed, 63 insertions(+), 28 deletions(-) diff --git a/src/arkode/arkode_adapt.c b/src/arkode/arkode_adapt.c index 9f91b684de..88ff02fe3d 100644 --- a/src/arkode/arkode_adapt.c +++ b/src/arkode/arkode_adapt.c @@ -155,7 +155,11 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur, /* increment the relevant step counter, set desired step */ if (h_acc <= h_cfl) { hadapt_mem->nst_acc++; } - else { hadapt_mem->nst_exp++; h_acc = h_cfl; } + else + { + hadapt_mem->nst_exp++; + h_acc = h_cfl; + } /* enforce adaptivity bounds to retain Jacobian/preconditioner accuracy */ if (dsm <= ONE) diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index 2f0d0e71e1..17ef9493dc 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -90,7 +90,8 @@ SUNAdaptController SUNAdaptController_Soderlind(SUNContext sunctx) C->ops->space = SUNAdaptController_Space_Soderlind; /* Create content */ - C->content = (SUNAdaptControllerContent_Soderlind)malloc(sizeof(struct _SUNAdaptControllerContent_Soderlind)); + C->content = (SUNAdaptControllerContent_Soderlind)malloc( + sizeof(struct _SUNAdaptControllerContent_Soderlind)); SUNAssertNull(C->content, SUN_ERR_MALLOC_FAIL); /* Fill content with default/reset values */ @@ -303,7 +304,8 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, /* order parameter to use */ const int ord = p + 1; - if (SODERLIND_FIRSTSTEPS(C) > 1) { + if (SODERLIND_FIRSTSTEPS(C) > 1) + { /* After the first 2 steps, there is sufficient history */ const sunrealtype k1 = -SODERLIND_K1(C) / ord; const sunrealtype k2 = -SODERLIND_K2(C) / ord; @@ -318,7 +320,9 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, *hnew = h * SUNRpowerR(e1, k1) * SUNRpowerR(e2, k2) * SUNRpowerR(e3, k3) * SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5); - } else { + } + else + { /* Use an I controller on the first two steps */ const sunrealtype k = -SUN_RCONST(1.0) / ord; const sunrealtype e = SODERLIND_BIAS(C) * dsm; diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index 01fa551559..874378f663 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -11,7 +11,8 @@ * SPDX-License-Identifier: BSD-3-Clause * SUNDIALS Copyright End * ----------------------------------------------------------------------------- - * Unit test for GetUserData functions + * Unit test for step size growth and handling of inf/nan within a controller + * due to 0 local error. * ---------------------------------------------------------------------------*/ #include @@ -20,6 +21,7 @@ #include "arkode/arkode_erkstep.h" #include "nvector/nvector_serial.h" +/* If an error occurs, print to stderr and exit */ static void err_fn(const int line, const char* const func, const char* const file, const char* const msg, const SUNErrCode err_code, void* const err_user_data, const SUNContext sunctx) @@ -28,7 +30,7 @@ static void err_fn(const int line, const char* const func, const char* const fil exit(err_code); } -// RHS for the simple ODE y' = 0 +/* RHS for the simple ODE y'=0 */ static int f(const sunrealtype t, const N_Vector y, const N_Vector ydot, void* const user_data) { @@ -36,13 +38,16 @@ static int f(const sunrealtype t, const N_Vector y, const N_Vector ydot, return 0; } +/* Take a single time step solving y'=0 and check the max growth was attained */ static int check_step(void* const arkode_mem, const N_Vector y, const sunrealtype h_expected, const int step) { + /* Integration much farther than expected step */ + const sunrealtype tout = SUN_RCONST(100.0) * h_expected; sunrealtype tret; /* The ERK method should be able to take the maximum possible timestep for the - simple ODE y' = 0 without any rejected steps or local error */ - ARKodeEvolve(arkode_mem, SUN_RCONST(1.0), y, &tret, ARK_ONE_STEP); + without any rejected steps or local error */ + ARKodeEvolve(arkode_mem, tout, y, &tret, ARK_ONE_STEP); long int local_err_fails; ARKodeGetNumErrTestFails(arkode_mem, &local_err_fails); @@ -66,7 +71,7 @@ static int check_step(void* const arkode_mem, const N_Vector y, sunrealtype h_actual; ARKodeGetCurrentStep(arkode_mem, &h_actual); - if (h_expected != h_actual) + if (SUNRCompare(h_expected, h_actual)) { fprintf(stderr, "Expected h at step %i to be %g but is %g\n", step, h_expected, h_actual); @@ -76,6 +81,29 @@ static int check_step(void* const arkode_mem, const N_Vector y, return 0; } +/* Take several steps solving y'=0 and check the max growth was attained */ +static int check_steps(void* const arkode_mem, const N_Vector y, + const sunrealtype h0, const sunrealtype first_growth, + const sunrealtype growth) +{ + ARKodeSetInitStep(arkode_mem, h0); + ARKodeSetMaxFirstGrowth(arkode_mem, first_growth); + ARKodeSetMaxGrowth(arkode_mem, growth); + + sunrealtype h_expect = first_growth * h0; + /* Take several steps to see the special behavior at step one then to allow + the adaptivity controller history fill up */ + const int num_steps = 5; + for (int step = 1; step <= num_steps; step++) + { + const int retval = check_step_no_err(arkode_mem, y, h_expect, step); + if (retval != 0) { return retval; } + h_expect *= growth; + } + + return 0; +} + int main() { SUNContext sunctx; @@ -96,26 +124,25 @@ int main() const N_Vector y = N_VNew_Serial(1, sunctx); N_VConst(SUN_RCONST(1.0), y); + /* Forward integration from 0 */ void* arkode_mem = ERKStepCreate(f, SUN_RCONST(0.0), y, sunctx); - - const sunrealtype h0 = SUN_RCONST(1.0e-4); - const sunrealtype first_growth = SUN_RCONST(1234.0); - const sunrealtype growth = SUN_RCONST(3.0); - - ARKodeSetInitStep(arkode_mem, h0); - ARKodeSetMaxFirstGrowth(arkode_mem, first_growth); - ARKodeSetMaxGrowth(arkode_mem, growth); - - sunrealtype h_expect = first_growth * h0; - /* Take several steps to see the special behavior at step one then to allow - the adaptivity controller history fill up */ - const int num_steps = 5; - for (int step = 1; step <= num_steps; step++) - { - retval = check_step(arkode_mem, y, h_expect, step); - if (retval != 0) { return retval; } - h_expect *= growth; - } + retval = check_steps_no_err(arkode_mem, y, SUN_RCONST(1.0e-4), + SUN_RCONST(1234.0), SUN_RCONST(3.0)); + if (retval != 0) { return retval; } + + /* Backward integration from positive time */ + ERKStepReInit(arkode_mem, f, SUN_RCONST(999.0), y); + retval = check_steps_no_err(arkode_mem, y, SUN_RCONST(-1.0e-2), + SUN_RCONST(1.6), SUN_RCONST(2.3)); + if (retval != 0) { return retval; } + + /* Forward integration from a negative time */ + ERKStepReInit(arkode_mem, f, SUN_RCONST(-999.0), y); + retval = check_steps_no_err(arkode_mem, y, SUN_RCONST(20.0), + SUN_RCONST(1.0e5), SUN_RCONST(1.1e3)); + if (retval != 0) { return retval; } + + /* TODO(SBR): add additional tests for more the default controller */ ARKodeFree(&arkode_mem); N_VDestroy(y); From 2452bd691d93fa14fde2e2d3fa7249f803609f80 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 6 Sep 2024 14:20:36 -0700 Subject: [PATCH 22/32] Remove unnecessary header --- src/sundials/sundials_math.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/sundials/sundials_math.c b/src/sundials/sundials_math.c index c78bb59068..4036149cee 100644 --- a/src/sundials/sundials_math.c +++ b/src/sundials/sundials_math.c @@ -22,8 +22,6 @@ #include #include -#include - sunrealtype SUNRpowerI(sunrealtype base, int exponent) { int i, expt; From 13cce1285cfee6debba7163f3daf9c7f11b235ab Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 6 Sep 2024 14:39:50 -0700 Subject: [PATCH 23/32] Correct function names --- test/unit_tests/arkode/C_serial/ark_test_adapt.c | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index 874378f663..72c6d493f0 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -96,7 +96,7 @@ static int check_steps(void* const arkode_mem, const N_Vector y, const int num_steps = 5; for (int step = 1; step <= num_steps; step++) { - const int retval = check_step_no_err(arkode_mem, y, h_expect, step); + const int retval = check_step(arkode_mem, y, h_expect, step); if (retval != 0) { return retval; } h_expect *= growth; } @@ -126,20 +126,20 @@ int main() /* Forward integration from 0 */ void* arkode_mem = ERKStepCreate(f, SUN_RCONST(0.0), y, sunctx); - retval = check_steps_no_err(arkode_mem, y, SUN_RCONST(1.0e-4), - SUN_RCONST(1234.0), SUN_RCONST(3.0)); + retval = check_steps(arkode_mem, y, SUN_RCONST(1.0e-4), SUN_RCONST(1234.0), + SUN_RCONST(3.0)); if (retval != 0) { return retval; } /* Backward integration from positive time */ ERKStepReInit(arkode_mem, f, SUN_RCONST(999.0), y); - retval = check_steps_no_err(arkode_mem, y, SUN_RCONST(-1.0e-2), - SUN_RCONST(1.6), SUN_RCONST(2.3)); + retval = check_steps(arkode_mem, y, SUN_RCONST(-1.0e-2), SUN_RCONST(1.6), + SUN_RCONST(2.3)); if (retval != 0) { return retval; } /* Forward integration from a negative time */ ERKStepReInit(arkode_mem, f, SUN_RCONST(-999.0), y); - retval = check_steps_no_err(arkode_mem, y, SUN_RCONST(20.0), - SUN_RCONST(1.0e5), SUN_RCONST(1.1e3)); + retval = check_steps(arkode_mem, y, SUN_RCONST(20.0), SUN_RCONST(1.0e5), + SUN_RCONST(1.1e3)); if (retval != 0) { return retval; } /* TODO(SBR): add additional tests for more the default controller */ From 19503568fd4599cacb07ee017f074dc172d832e2 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 6 Sep 2024 15:05:30 -0700 Subject: [PATCH 24/32] Convert SUNRpowerR to macro --- include/sundials/sundials_math.h | 39 ++++++++++++++++++++------------ src/sundials/sundials_math.c | 14 ------------ 2 files changed, 24 insertions(+), 29 deletions(-) diff --git a/include/sundials/sundials_math.h b/include/sundials/sundials_math.h index bc6cd2f9eb..c89be0af1e 100644 --- a/include/sundials/sundials_math.h +++ b/include/sundials/sundials_math.h @@ -184,33 +184,42 @@ extern "C" { /* * ----------------------------------------------------------------- - * Function : SUNRpowerI + * Function : SUNRpowerR * ----------------------------------------------------------------- - * Usage : int exponent; - * sunrealtype base, ans; - * ans = SUNRpowerI(base,exponent); + * Usage : sunrealtype base, exponent, ans; + * ans = SUNRpowerR(base,exponent); * ----------------------------------------------------------------- - * SUNRpowerI returns the value of base^exponent, where base is of type - * sunrealtype and exponent is of type int. + * SUNRpowerR returns the value of base^exponent, where both base and + * exponent are of type sunrealtype. * ----------------------------------------------------------------- */ - -SUNDIALS_EXPORT sunrealtype SUNRpowerI(sunrealtype base, int exponent); +#ifndef SUNRpowerR +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define SUNRpowerR(base, exponent) (pow(base, exponent)) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define SUNRpowerR(base, exponent) (powf(base, exponent)) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define SUNRpowerR(base, exponent) (powl(base, exponent)) +#else +#error \ + "SUNDIALS precision not defined, report to github.com/LLNL/sundials/issues" +#endif +#endif /* * ----------------------------------------------------------------- - * Function : SUNRpowerR + * Function : SUNRpowerI * ----------------------------------------------------------------- - * Usage : sunrealtype base, exponent, ans; - * ans = SUNRpowerR(base,exponent); + * Usage : int exponent; + * sunrealtype base, ans; + * ans = SUNRpowerI(base,exponent); * ----------------------------------------------------------------- - * SUNRpowerR returns the value of base^exponent, where both base and - * exponent are of type sunrealtype. If base < ZERO, then SUNRpowerR - * returns ZERO. + * SUNRpowerI returns the value of base^exponent, where base is of type + * sunrealtype and exponent is of type int. * ----------------------------------------------------------------- */ -SUNDIALS_EXPORT sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent); +SUNDIALS_EXPORT sunrealtype SUNRpowerI(sunrealtype base, int exponent); /* * ----------------------------------------------------------------- diff --git a/src/sundials/sundials_math.c b/src/sundials/sundials_math.c index 4036149cee..4e4122c7be 100644 --- a/src/sundials/sundials_math.c +++ b/src/sundials/sundials_math.c @@ -34,20 +34,6 @@ sunrealtype SUNRpowerI(sunrealtype base, int exponent) return (prod); } -sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent) -{ -#if defined(SUNDIALS_DOUBLE_PRECISION) - return (pow(base, exponent)); -#elif defined(SUNDIALS_SINGLE_PRECISION) - return (powf(base, exponent)); -#elif defined(SUNDIALS_EXTENDED_PRECISION) - return (powl(base, exponent)); -#else -#error \ - "SUNDIALS precision not defined, report to github.com/LLNL/sundials/issues" -#endif -} - sunbooleantype SUNRCompare(sunrealtype a, sunrealtype b) { return (SUNRCompareTol(a, b, 10 * SUN_UNIT_ROUNDOFF)); From 5bd47904b0de5732cfe0fcc8b5af846c29927863 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 9 Sep 2024 21:29:50 -0700 Subject: [PATCH 25/32] Update comments based on @drreynolds suggestions --- test/unit_tests/arkode/C_serial/ark_test_adapt.c | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index 72c6d493f0..e7510bcb08 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -45,8 +45,8 @@ static int check_step(void* const arkode_mem, const N_Vector y, /* Integration much farther than expected step */ const sunrealtype tout = SUN_RCONST(100.0) * h_expected; sunrealtype tret; - /* The ERK method should be able to take the maximum possible timestep for the - without any rejected steps or local error */ + /* The ERK method should be able to take the maximum possible timestep without + any rejected steps or local error */ ARKodeEvolve(arkode_mem, tout, y, &tret, ARK_ONE_STEP); long int local_err_fails; @@ -91,8 +91,7 @@ static int check_steps(void* const arkode_mem, const N_Vector y, ARKodeSetMaxGrowth(arkode_mem, growth); sunrealtype h_expect = first_growth * h0; - /* Take several steps to see the special behavior at step one then to allow - the adaptivity controller history fill up */ + /* Take enough steps to allow the controller history to fill up */ const int num_steps = 5; for (int step = 1; step <= num_steps; step++) { @@ -142,7 +141,7 @@ int main() SUN_RCONST(1.1e3)); if (retval != 0) { return retval; } - /* TODO(SBR): add additional tests for more the default controller */ + /* TODO(SBR): test non-default controllers */ ARKodeFree(&arkode_mem); N_VDestroy(y); From 8d7a2678426ba60b8334946dedc224d0af7cee7b Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 10 Sep 2024 16:11:07 -0700 Subject: [PATCH 26/32] Fix typos in changelog --- CHANGELOG.md | 2 +- doc/shared/RecentChanges.rst | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 641d5e6a2f..6d5f294e36 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,7 +18,7 @@ Removed error floors from the SUNAdaptController implementations which could unnecessarily limit the time size growth, particularly after the first step. On the first two time steps, the Soderlind controller uses an I controller -instead omitting unavailable terms. +instead of omitting unavailable terms. Fixed the loading of ARKStep's default first order explicit method. diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index d77b9d1e20..ba3585e20c 100644 --- a/doc/shared/RecentChanges.rst +++ b/doc/shared/RecentChanges.rst @@ -15,7 +15,7 @@ unnecessarily limit the time size growth, particularly after the first step. On the first two time steps, the :ref:`Soderlind controller ` uses an I controller -instead omitting unavailable terms. +instead of omitting unavailable terms. Fixed the loading of ARKStep's default first order explicit method. From 7d5161b64033fdee8f153296040d4c49284c761e Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 10 Sep 2024 19:48:37 -0700 Subject: [PATCH 27/32] Update test/unit_tests/arkode/C_serial/ark_test_adapt.c Co-authored-by: David Gardner --- test/unit_tests/arkode/C_serial/ark_test_adapt.c | 1 + 1 file changed, 1 insertion(+) diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index e7510bcb08..b8fe34bf03 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -55,6 +55,7 @@ static int check_step(void* const arkode_mem, const N_Vector y, { fprintf(stderr, "Expected 0 local error failures at step %i but is %li\n", step, local_err_fails); + return 1; } const N_Vector err = N_VClone(y); From fcb7dc6ce2acef1bc4899b7eb552c91410d165f7 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 10 Sep 2024 19:56:14 -0700 Subject: [PATCH 28/32] Add missing SUN_RCONST --- test/unit_tests/arkode/C_serial/ark_test_adapt.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index b8fe34bf03..d3d47a7f15 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -63,7 +63,7 @@ static int check_step(void* const arkode_mem, const N_Vector y, const sunrealtype err_norm = N_VMaxNorm(err); N_VDestroy(err); - if (err_norm != 0) + if (err_norm != SUN_RCONST(0)) { fprintf(stderr, "Expected local error at step %i to be 0 but is %g\n", step, err_norm); From 1db5b33e8f411b151a44a3a504acfa9d692ddf72 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 10 Sep 2024 20:10:02 -0700 Subject: [PATCH 29/32] Update test/unit_tests/arkode/C_serial/ark_test_adapt.c Co-authored-by: David Gardner --- test/unit_tests/arkode/C_serial/ark_test_adapt.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index d3d47a7f15..deb282dcb6 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -63,7 +63,7 @@ static int check_step(void* const arkode_mem, const N_Vector y, const sunrealtype err_norm = N_VMaxNorm(err); N_VDestroy(err); - if (err_norm != SUN_RCONST(0)) + if (err_norm != SUN_RCONST(0.0)) { fprintf(stderr, "Expected local error at step %i to be 0 but is %g\n", step, err_norm); From 7669c933c7416a6202c41b4cbd97ffb3be7a4515 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 10 Sep 2024 21:00:16 -0700 Subject: [PATCH 30/32] Add more asserts to controllers --- .../imexgus/sunadaptcontroller_imexgus.c | 2 -- .../soderlind/sunadaptcontroller_soderlind.c | 1 - src/sundials/sundials_adaptcontroller.c | 7 +++++++ 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c index 650ffc52df..0bbcda4b0a 100644 --- a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c +++ b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c @@ -127,8 +127,6 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C, { SUNFunctionBegin(C->sunctx); - SUNAssert(hnew, SUN_ERR_ARG_CORRUPT); - /* order parameter to use */ const int ord = p + 1; diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c index 17ef9493dc..ebee48b849 100644 --- a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c +++ b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c @@ -299,7 +299,6 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C, sunrealtype* hnew) { SUNFunctionBegin(C->sunctx); - SUNAssert(hnew, SUN_ERR_ARG_CORRUPT); /* order parameter to use */ const int ord = p + 1; diff --git a/src/sundials/sundials_adaptcontroller.c b/src/sundials/sundials_adaptcontroller.c index 39a56480d0..925d3fae50 100644 --- a/src/sundials/sundials_adaptcontroller.c +++ b/src/sundials/sundials_adaptcontroller.c @@ -18,6 +18,7 @@ #include #include +#include #include "sundials/sundials_errors.h" @@ -131,6 +132,9 @@ SUNErrCode SUNAdaptController_EstimateStep(SUNAdaptController C, sunrealtype h, SUNErrCode ier = SUN_SUCCESS; if (C == NULL) { return SUN_ERR_ARG_CORRUPT; } SUNFunctionBegin(C->sunctx); + SUNAssert(isfinite(h), SUN_ERR_ARG_OUTOFRANGE); + SUNAssert(p > 0, SUN_ERR_ARG_OUTOFRANGE); + SUNAssert(dsm >= SUN_RCONST(0.0), SUN_ERR_ARG_OUTOFRANGE); SUNAssert(hnew, SUN_ERR_ARG_CORRUPT); *hnew = h; /* initialize output with identity */ if (C->ops->estimatestep) { ier = C->ops->estimatestep(C, h, p, dsm, hnew); } @@ -170,6 +174,7 @@ SUNErrCode SUNAdaptController_SetErrorBias(SUNAdaptController C, sunrealtype bia SUNErrCode ier = SUN_SUCCESS; if (C == NULL) { return SUN_ERR_ARG_CORRUPT; } SUNFunctionBegin(C->sunctx); + SUNAssert(bias >= SUN_RCONST(1.0), SUN_ERR_ARG_OUTOFRANGE); if (C->ops->seterrorbias) { ier = C->ops->seterrorbias(C, bias); } return (ier); } @@ -180,6 +185,8 @@ SUNErrCode SUNAdaptController_UpdateH(SUNAdaptController C, sunrealtype h, SUNErrCode ier = SUN_SUCCESS; if (C == NULL) { return SUN_ERR_ARG_CORRUPT; } SUNFunctionBegin(C->sunctx); + SUNAssert(isfinite(h), SUN_ERR_ARG_OUTOFRANGE); + SUNAssert(dsm >= SUN_RCONST(0.0), SUN_ERR_ARG_OUTOFRANGE); if (C->ops->updateh) { ier = C->ops->updateh(C, h, dsm); } return (ier); } From 3671f0e7b8e812f2e672796e6b3bbf05c06608f2 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 26 Sep 2024 11:37:41 -0700 Subject: [PATCH 31/32] Add non-default controllers in test --- src/sundials/sundials_adaptcontroller.c | 2 +- .../arkode/C_serial/ark_test_adapt.c | 20 ++++++++++++++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/sundials/sundials_adaptcontroller.c b/src/sundials/sundials_adaptcontroller.c index 925d3fae50..43430bd363 100644 --- a/src/sundials/sundials_adaptcontroller.c +++ b/src/sundials/sundials_adaptcontroller.c @@ -16,9 +16,9 @@ * operations listed in sundials_adaptcontroller.h * -----------------------------------------------------------------*/ +#include #include #include -#include #include "sundials/sundials_errors.h" diff --git a/test/unit_tests/arkode/C_serial/ark_test_adapt.c b/test/unit_tests/arkode/C_serial/ark_test_adapt.c index deb282dcb6..8885668e2d 100644 --- a/test/unit_tests/arkode/C_serial/ark_test_adapt.c +++ b/test/unit_tests/arkode/C_serial/ark_test_adapt.c @@ -142,8 +142,26 @@ int main() SUN_RCONST(1.1e3)); if (retval != 0) { return retval; } - /* TODO(SBR): test non-default controllers */ + /* Try a non-default controller */ + ERKStepReInit(arkode_mem, f, SUN_RCONST(0.0), y); + SUNAdaptController controller1 = SUNAdaptController_ExpGus(sunctx); + ARKodeSetAdaptController(arkode_mem, controller1); + retval = check_steps(arkode_mem, y, SUN_RCONST(0.1), SUN_RCONST(4.0), + SUN_RCONST(10.0)); + if (retval != 0) { return retval; } + + /* Try another non-default controller */ + ERKStepReInit(arkode_mem, f, SUN_RCONST(0.0), y); + SUNAdaptController controller2 = SUNAdaptController_Soderlind(sunctx); + SUNAdaptController_SetParams_PI(controller2, SUN_RCONST(0.123), + SUN_RCONST(-0.456)); + ARKodeSetAdaptController(arkode_mem, controller2); + retval = check_steps(arkode_mem, y, SUN_RCONST(-0.1), SUN_RCONST(4.0), + SUN_RCONST(10.0)); + if (retval != 0) { return retval; } + SUNAdaptController_Destroy(controller1); + SUNAdaptController_Destroy(controller2); ARKodeFree(&arkode_mem); N_VDestroy(y); SUNContext_Free(&sunctx); From fdb3d4a18e61b4fb3c139832b7ef01b740bf81ca Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 26 Sep 2024 13:37:41 -0700 Subject: [PATCH 32/32] Fix assertion to allow order 0 --- src/sundials/sundials_adaptcontroller.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sundials/sundials_adaptcontroller.c b/src/sundials/sundials_adaptcontroller.c index 43430bd363..124ad3f15b 100644 --- a/src/sundials/sundials_adaptcontroller.c +++ b/src/sundials/sundials_adaptcontroller.c @@ -133,7 +133,7 @@ SUNErrCode SUNAdaptController_EstimateStep(SUNAdaptController C, sunrealtype h, if (C == NULL) { return SUN_ERR_ARG_CORRUPT; } SUNFunctionBegin(C->sunctx); SUNAssert(isfinite(h), SUN_ERR_ARG_OUTOFRANGE); - SUNAssert(p > 0, SUN_ERR_ARG_OUTOFRANGE); + SUNAssert(p >= 0, SUN_ERR_ARG_OUTOFRANGE); SUNAssert(dsm >= SUN_RCONST(0.0), SUN_ERR_ARG_OUTOFRANGE); SUNAssert(hnew, SUN_ERR_ARG_CORRUPT); *hnew = h; /* initialize output with identity */