Skip to content

Commit

Permalink
Avoid roundoff trouble in choose_lts_step_size
Browse files Browse the repository at this point in the history
  • Loading branch information
wthrowe committed Aug 26, 2024
1 parent f8b65d5 commit 3ef8cb3
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 4 deletions.
9 changes: 5 additions & 4 deletions src/Time/ChooseLtsStepSize.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,10 +18,11 @@ TimeDelta choose_lts_step_size(const Time& time, const double desired_step) {
const TimeDelta full_slab =
desired_step > 0.0 ? time.slab().duration() : -time.slab().duration();
const double desired_step_count = full_slab.value() / desired_step;
const size_t desired_step_power =
desired_step_count <= 1
? 0
: static_cast<size_t>(std::ceil(std::log2(desired_step_count)));
// log2 is a slowly increasing function, so log2(2^n + eps) may give
// n because of floating-point truncation. The inner ceil call
// avoids this problem.
const auto desired_step_power =
static_cast<size_t>(std::ceil(std::log2(std::ceil(desired_step_count))));

// Ensure we will hit the slab boundary if we continue taking
// constant-sized steps.
Expand Down
15 changes: 15 additions & 0 deletions tests/Unit/Time/Actions/Test_ChangeStepSize.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

#include "Framework/TestingFramework.hpp"

#include <limits>
#include <memory>
#include <pup.h>
#include <string>
Expand Down Expand Up @@ -207,6 +208,20 @@ SPECTRE_TEST_CASE("Unit.Time.Actions.ChangeStepSize", "[Unit][Time][Actions]") {
reject_step ? -slab.duration() / 8 : -slab.duration() / 4,
reject_step ? std::optional{std::make_unique<StepRejector>(0.5)}
: std::nullopt);

// Check for roundoff issues
check(true, std::make_unique<TimeSteppers::AdamsBashforth>(1), {},
slab.start() + slab.duration() / 4,
slab_length / 16. / (1.0 + std::numeric_limits<double>::epsilon()),
slab.duration() / 32,
reject_step ? std::optional{std::make_unique<StepRejector>(0.5)}
: std::nullopt);
check(false, std::make_unique<TimeSteppers::AdamsBashforth>(1), {},
slab.end() - slab.duration() / 4,
slab_length / 16. / (1.0 + std::numeric_limits<double>::epsilon()),
-slab.duration() / 32,
reject_step ? std::optional{std::make_unique<StepRejector>(0.5)}
: std::nullopt);
}

{
Expand Down

0 comments on commit 3ef8cb3

Please sign in to comment.