diff --git a/src/Evolution/Executables/Cce/CharacteristicExtractBase.hpp b/src/Evolution/Executables/Cce/CharacteristicExtractBase.hpp index 449c2c06b95a..d67f929f4dba 100644 --- a/src/Evolution/Executables/Cce/CharacteristicExtractBase.hpp +++ b/src/Evolution/Executables/Cce/CharacteristicExtractBase.hpp @@ -14,7 +14,7 @@ #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshTags.hpp" #include "Time/StepChoosers/Constant.hpp" #include "Time/StepChoosers/ErrorControl.hpp" -#include "Time/StepChoosers/Increase.hpp" +#include "Time/StepChoosers/LimitIncrease.hpp" #include "Utilities/TMPL.hpp" template @@ -112,7 +112,7 @@ struct CharacteristicExtractDefaults { tmpl::list>; using cce_step_choosers = tmpl::list, - StepChoosers::Increase, + StepChoosers::LimitIncrease, StepChoosers::ErrorControl, swsh_vars_selector>, diff --git a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp index d81df2741b05..2d7de8172ef3 100644 --- a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp +++ b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp @@ -164,13 +164,8 @@ #include "Time/Actions/UpdateU.hpp" #include "Time/ChangeSlabSize/Action.hpp" #include "Time/ChangeSlabSize/Tags.hpp" -#include "Time/StepChoosers/Cfl.hpp" -#include "Time/StepChoosers/Constant.hpp" #include "Time/StepChoosers/Factory.hpp" -#include "Time/StepChoosers/Increase.hpp" -#include "Time/StepChoosers/PreventRapidIncrease.hpp" #include "Time/StepChoosers/StepChooser.hpp" -#include "Time/StepChoosers/StepToTimes.hpp" #include "Time/Tags/StepperErrors.hpp" #include "Time/Tags/Time.hpp" #include "Time/Tags/TimeStepId.hpp" diff --git a/src/Evolution/Executables/GeneralizedHarmonic/GeneralizedHarmonicBase.hpp b/src/Evolution/Executables/GeneralizedHarmonic/GeneralizedHarmonicBase.hpp index 9e94f95d22a2..64ea0abcdd6c 100644 --- a/src/Evolution/Executables/GeneralizedHarmonic/GeneralizedHarmonicBase.hpp +++ b/src/Evolution/Executables/GeneralizedHarmonic/GeneralizedHarmonicBase.hpp @@ -131,13 +131,8 @@ #include "Time/Actions/RecordTimeStepperData.hpp" #include "Time/Actions/SelfStartActions.hpp" #include "Time/Actions/UpdateU.hpp" -#include "Time/StepChoosers/Cfl.hpp" -#include "Time/StepChoosers/Constant.hpp" #include "Time/StepChoosers/Factory.hpp" -#include "Time/StepChoosers/Increase.hpp" -#include "Time/StepChoosers/PreventRapidIncrease.hpp" #include "Time/StepChoosers/StepChooser.hpp" -#include "Time/StepChoosers/StepToTimes.hpp" #include "Time/Tags/Time.hpp" #include "Time/TimeSequence.hpp" #include "Time/TimeSteppers/Factory.hpp" diff --git a/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp b/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp index 9046d648334c..bae017435b34 100644 --- a/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp +++ b/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp @@ -221,13 +221,8 @@ #include "Time/Actions/SelfStartActions.hpp" // IWYU pragma: keep #include "Time/Actions/UpdateU.hpp" #include "Time/ChangeSlabSize/Action.hpp" -#include "Time/StepChoosers/Cfl.hpp" -#include "Time/StepChoosers/Constant.hpp" #include "Time/StepChoosers/Factory.hpp" -#include "Time/StepChoosers/Increase.hpp" -#include "Time/StepChoosers/PreventRapidIncrease.hpp" #include "Time/StepChoosers/StepChooser.hpp" -#include "Time/StepChoosers/StepToTimes.hpp" #include "Time/Tags/Time.hpp" #include "Time/Tags/TimeStepId.hpp" #include "Time/TimeSequence.hpp" diff --git a/src/Evolution/Executables/ScalarTensor/ScalarTensorBase.hpp b/src/Evolution/Executables/ScalarTensor/ScalarTensorBase.hpp index 3e6b48e2e583..a1ffec5b349d 100644 --- a/src/Evolution/Executables/ScalarTensor/ScalarTensorBase.hpp +++ b/src/Evolution/Executables/ScalarTensor/ScalarTensorBase.hpp @@ -132,13 +132,8 @@ #include "Time/Actions/RecordTimeStepperData.hpp" #include "Time/Actions/SelfStartActions.hpp" #include "Time/Actions/UpdateU.hpp" -#include "Time/StepChoosers/Cfl.hpp" -#include "Time/StepChoosers/Constant.hpp" #include "Time/StepChoosers/Factory.hpp" -#include "Time/StepChoosers/Increase.hpp" -#include "Time/StepChoosers/PreventRapidIncrease.hpp" #include "Time/StepChoosers/StepChooser.hpp" -#include "Time/StepChoosers/StepToTimes.hpp" #include "Time/Tags/Time.hpp" #include "Time/TimeSequence.hpp" #include "Time/TimeSteppers/Factory.hpp" diff --git a/src/Evolution/Initialization/Evolution.hpp b/src/Evolution/Initialization/Evolution.hpp index de6bfcdea451..6f832f623b9e 100644 --- a/src/Evolution/Initialization/Evolution.hpp +++ b/src/Evolution/Initialization/Evolution.hpp @@ -26,6 +26,7 @@ #include "Parallel/Tags/ArrayIndex.hpp" #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp" #include "Time/AdaptiveSteppingDiagnostics.hpp" +#include "Time/ChangeSlabSize/Tags.hpp" #include "Time/ChooseLtsStepSize.hpp" #include "Time/Slab.hpp" #include "Time/StepChoosers/StepChooser.hpp" @@ -114,7 +115,8 @@ struct TimeStepping { /// Tags for items in the DataBox that are mutated by the apply function using return_tags = tmpl::list<::Tags::Next<::Tags::TimeStepId>, ::Tags::TimeStep, - ::Tags::Next<::Tags::TimeStep>>; + ::Tags::Next<::Tags::TimeStep>, + ::Tags::ChangeSlabSize::SlabSizeGoal>; /// Tags for mutable DataBox items that are either default initialized or /// initialized by the apply function @@ -130,6 +132,7 @@ struct TimeStepping { static void apply(const gsl::not_null next_time_step_id, const gsl::not_null time_step, const gsl::not_null next_time_step, + const gsl::not_null slab_size_goal, const double initial_time_value, const double initial_dt_value, const double initial_slab_size, @@ -141,6 +144,8 @@ struct TimeStepping { time_runs_forward, time_stepper); *time_step = choose_lts_step_size(initial_time, initial_dt_value); *next_time_step = *time_step; + *slab_size_goal = + time_runs_forward ? initial_slab_size : -initial_slab_size; } /// Given the items fetched from a DataBox by the argument_tags, when not @@ -148,6 +153,7 @@ struct TimeStepping { static void apply(const gsl::not_null next_time_step_id, const gsl::not_null time_step, const gsl::not_null next_time_step, + const gsl::not_null slab_size_goal, const double initial_time_value, const double initial_dt_value, const double initial_slab_size, @@ -159,6 +165,8 @@ struct TimeStepping { time_runs_forward, time_stepper); *time_step = (time_runs_forward ? 1 : -1) * initial_time.slab().duration(); *next_time_step = *time_step; + *slab_size_goal = + time_runs_forward ? initial_slab_size : -initial_slab_size; } }; @@ -168,7 +176,8 @@ struct ProjectTimeStepping : tt::ConformsTo { using return_tags = tmpl::list<::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>, ::Tags::TimeStep, ::Tags::Next<::Tags::TimeStep>, ::Tags::Time, - ::Tags::AdaptiveSteppingDiagnostics>; + ::Tags::AdaptiveSteppingDiagnostics, + ::Tags::ChangeSlabSize::SlabSizeGoal>; using argument_tags = tmpl::list; static void apply( @@ -179,6 +188,7 @@ struct ProjectTimeStepping : tt::ConformsTo { const gsl::not_null /*time*/, const gsl::not_null /*adaptive_stepping_diagnostics*/, + const gsl::not_null /*slab_size_goal*/, const ElementId& /*element_id*/, const std::pair, Element>& /*old_mesh_and_element*/) { // Do not change anything for p-refinement @@ -192,6 +202,7 @@ struct ProjectTimeStepping : tt::ConformsTo { const gsl::not_null time, const gsl::not_null adaptive_stepping_diagnostics, + const gsl::not_null slab_size_goal, const ElementId& element_id, const tuples::TaggedTuple& parent_items) { *time_step_id = get<::Tags::TimeStepId>(parent_items); @@ -199,6 +210,7 @@ struct ProjectTimeStepping : tt::ConformsTo { *time_step = get<::Tags::TimeStep>(parent_items); *next_time_step = get<::Tags::Next<::Tags::TimeStep>>(parent_items); *time = get<::Tags::Time>(parent_items); + *slab_size_goal = get<::Tags::ChangeSlabSize::SlabSizeGoal>(parent_items); // Since AdaptiveSteppingDiagnostics are reduced over all elements, we // set the slab quantities to the same value over all children, and the @@ -229,6 +241,7 @@ struct ProjectTimeStepping : tt::ConformsTo { const gsl::not_null time, const gsl::not_null adaptive_stepping_diagnostics, + const gsl::not_null slab_size_goal, const ElementId& /*element_id*/, const std::unordered_map, tuples::TaggedTuple>& children_items) { @@ -248,6 +261,8 @@ struct ProjectTimeStepping : tt::ConformsTo { *time_step = get<::Tags::TimeStep>(slowest_child_items); *next_time_step = get<::Tags::Next<::Tags::TimeStep>>(slowest_child_items); *time = get<::Tags::Time>(slowest_child_items); + *slab_size_goal = + get<::Tags::ChangeSlabSize::SlabSizeGoal>(slowest_child_items); const auto& slowest_child_diagnostics = get<::Tags::AdaptiveSteppingDiagnostics>(slowest_child_items); diff --git a/src/Time/Actions/ChangeStepSize.hpp b/src/Time/Actions/ChangeStepSize.hpp index ff8982ace95c..123a30b0567b 100644 --- a/src/Time/Actions/ChangeStepSize.hpp +++ b/src/Time/Actions/ChangeStepSize.hpp @@ -15,6 +15,7 @@ #include "Time/ChooseLtsStepSize.hpp" #include "Time/Tags/AdaptiveSteppingDiagnostics.hpp" #include "Time/Tags/HistoryEvolvedVariables.hpp" +#include "Time/TimeStepRequestProcessor.hpp" #include "Time/TimeSteppers/LtsTimeStepper.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/Gsl.hpp" @@ -61,6 +62,20 @@ bool change_step_size(const gsl::not_null*> box) { const auto& time_step_id = db::get(*box); ASSERT(time_step_id.substep() == 0, "Can't change step size on a substep."); + + const auto& current_step = db::get(*box); + const double last_step_size = current_step.value(); + + TimeStepRequestProcessor step_requests(time_step_id.time_runs_forward()); + bool step_accepted = true; + for (const auto& step_chooser : step_choosers) { + const auto [step_request, step_choice_accepted] = + step_chooser->template desired_step(last_step_size, + *box); + step_requests.process(step_request); + step_accepted = step_accepted and step_choice_accepted; + } + using history_tags = ::Tags::get_all_history_tags; bool can_change_step_size = true; tmpl::for_each([&box, &can_change_step_size, &time_stepper, @@ -74,28 +89,14 @@ bool change_step_size(const gsl::not_null*> box) { time_stepper.can_change_step_size(time_step_id, history); }); if (not can_change_step_size) { + step_requests.error_on_hard_limit( + current_step.value(), + (time_step_id.step_time() + current_step).value()); return true; } - const auto& current_step = db::get(*box); - - const double last_step_size = std::abs(db::get(*box).value()); - - // The step choosers return the magnitude of the desired step, so - // we always want the minimum requirement, but we have to negate - // the final answer if time is running backwards. - double desired_step = std::numeric_limits::infinity(); - bool step_accepted = true; - for (const auto& step_chooser : step_choosers) { - const auto [step_choice, step_choice_accepted] = - step_chooser->template desired_step( - last_step_size, *box); - desired_step = std::min(desired_step, step_choice); - step_accepted = step_accepted and step_choice_accepted; - } - if (not current_step.is_positive()) { - desired_step = -desired_step; - } + const double desired_step = step_requests.step_size( + time_step_id.step_time().value(), current_step.value()); constexpr double smallest_relative_step_size = 1.0e-9; if (abs(desired_step / current_step.slab().duration().value()) < @@ -124,6 +125,9 @@ bool change_step_size(const gsl::not_null*> box) { // if step accepted, just proceed. Otherwise, change Time::Next and jump // back to the first instance of `UpdateU`. if (step_accepted) { + step_requests.error_on_hard_limit( + current_step.value(), + (time_step_id.step_time() + current_step).value()); return true; } else { db::mutate, Tags::TimeStep>( diff --git a/src/Time/CMakeLists.txt b/src/Time/CMakeLists.txt index 7b3876d3b8f0..b70d8d64d5a7 100644 --- a/src/Time/CMakeLists.txt +++ b/src/Time/CMakeLists.txt @@ -20,6 +20,8 @@ spectre_target_sources( Time.cpp TimeSequence.cpp TimeStepId.cpp + TimeStepRequest.cpp + TimeStepRequestProcessor.cpp Utilities.cpp ) @@ -42,6 +44,8 @@ spectre_target_headers( Time.hpp TimeSequence.hpp TimeStepId.hpp + TimeStepRequest.hpp + TimeStepRequestProcessor.hpp Utilities.hpp ) diff --git a/src/Time/ChangeSlabSize/Action.hpp b/src/Time/ChangeSlabSize/Action.hpp index b7f4a7ce2a11..ff35a2510a10 100644 --- a/src/Time/ChangeSlabSize/Action.hpp +++ b/src/Time/ChangeSlabSize/Action.hpp @@ -7,16 +7,18 @@ #include #include #include -#include +#include #include "DataStructures/DataBox/DataBox.hpp" #include "Parallel/AlgorithmExecution.hpp" #include "Time/ChangeSlabSize/ChangeSlabSize.hpp" #include "Time/ChangeSlabSize/Tags.hpp" #include "Time/Tags/HistoryEvolvedVariables.hpp" +#include "Time/TimeStepRequestProcessor.hpp" #include "Time/TimeSteppers/TimeStepper.hpp" #include "Utilities/Algorithm.hpp" #include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/Numeric.hpp" #include "Utilities/TMPL.hpp" /// \cond @@ -25,6 +27,7 @@ template class GlobalCache; } // namespace Parallel namespace Tags { +struct TimeStep; struct TimeStepId; template struct TimeStepper; @@ -73,78 +76,84 @@ struct ChangeSlabSize { return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } + TimeStepRequestProcessor step_requests(time_step_id.time_runs_forward()); + const auto slab_number = time_step_id.slab_number(); const auto& expected_messages_map = db::get<::Tags::ChangeSlabSize::NumberOfExpectedMessages>(box); - if (expected_messages_map.empty() or - expected_messages_map.begin()->first != slab_number) { - return {Parallel::AlgorithmExecution::Continue, std::nullopt}; - } - const size_t expected_messages = expected_messages_map.begin()->second; - ASSERT(expected_messages > 0, - "Should only create map entries when sending messages."); - - const auto& slab_size_messages = - db::get<::Tags::ChangeSlabSize::NewSlabSize>(box); - if (slab_size_messages.empty()) { - return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + if (not expected_messages_map.empty() and + expected_messages_map.begin()->first == slab_number) { + const size_t expected_messages = expected_messages_map.begin()->second; + ASSERT(expected_messages > 0, + "Should only create map entries when sending messages."); + + const auto& slab_size_messages = + db::get<::Tags::ChangeSlabSize::NewSlabSize>(box); + if (slab_size_messages.empty()) { + return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + } + const int64_t first_received_change_slab = + slab_size_messages.begin()->first; + + ASSERT(first_received_change_slab >= slab_number, + "Received data for a change at slab " << first_received_change_slab + << " but it is already slab " << slab_number); + if (first_received_change_slab != slab_number) { + return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + } + + const auto& received_changes = slab_size_messages.begin()->second; + ASSERT(expected_messages >= received_changes.size(), + "Received " << received_changes.size() + << " size change messages at slab " << slab_number + << ", but only expected " << expected_messages); + if (received_changes.size() != expected_messages) { + return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + } + + // We have all the data we need. + + step_requests = + alg::accumulate(slab_size_messages.begin()->second, step_requests); + + db::mutate<::Tags::ChangeSlabSize::NumberOfExpectedMessages, + ::Tags::ChangeSlabSize::NewSlabSize>( + [](const gsl::not_null*> expected, + const gsl::not_null< + std::map>*> + sizes) { + expected->erase(expected->begin()); + sizes->erase(sizes->begin()); + }, + make_not_null(&box)); } - const int64_t first_received_change_slab = - slab_size_messages.begin()->first; - - ASSERT(first_received_change_slab >= slab_number, - "Received data for a change at slab " << first_received_change_slab - << " but it is already slab " << slab_number); - if (first_received_change_slab != slab_number) { - return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + const double new_slab_end = step_requests.step_end( + time_step_id.step_time().value(), + db::get<::Tags::ChangeSlabSize::SlabSizeGoal>(box)); + + if (const auto new_goal = step_requests.new_step_size_goal(); + new_goal.has_value()) { + db::mutate<::Tags::ChangeSlabSize::SlabSizeGoal>( + [&](const gsl::not_null slab_size_goal) { + *slab_size_goal = *new_goal; + }, + make_not_null(&box)); } - const auto& received_changes = slab_size_messages.begin()->second; - ASSERT(expected_messages >= received_changes.size(), - "Received " << received_changes.size() - << " size change messages at slab " - << slab_number << ", but only expected " - << expected_messages); - if (received_changes.size() != expected_messages) { - return {Parallel::AlgorithmExecution::Retry, std::nullopt}; - } - - // We have all the data we need. - - const double new_slab_size = - *alg::min_element(slab_size_messages.begin()->second); - - db::mutate<::Tags::ChangeSlabSize::NumberOfExpectedMessages, - ::Tags::ChangeSlabSize::NewSlabSize>( - [](const gsl::not_null*> expected, - const gsl::not_null< - std::map>*> - sizes) { - expected->erase(expected->begin()); - sizes->erase(sizes->begin()); - }, - make_not_null(&box)); - const TimeStepper& time_stepper = db::get<::Tags::TimeStepper>(box); // Sometimes time steppers need to run with a fixed step size. // This is generally at the start of an evolution when the history // is in an unusual state. - if (not time_stepper.can_change_step_size( + if (time_stepper.can_change_step_size( time_step_id, db::get<::Tags::HistoryEvolvedVariables<>>(box))) { - return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + change_slab_size(make_not_null(&box), new_slab_end); } - const auto current_slab = time_step_id.step_time().slab(); - - const double new_slab_end = - time_step_id.time_runs_forward() - ? current_slab.start().value() + new_slab_size - : current_slab.end().value() - new_slab_size; - - change_slab_size(make_not_null(&box), new_slab_end); - + step_requests.error_on_hard_limit( + db::get<::Tags::TimeStep>(box).value(), + (time_step_id.step_time() + db::get<::Tags::TimeStep>(box)).value()); return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } }; diff --git a/src/Time/ChangeSlabSize/Event.hpp b/src/Time/ChangeSlabSize/Event.hpp index 5fa7658fb489..832fb959757f 100644 --- a/src/Time/ChangeSlabSize/Event.hpp +++ b/src/Time/ChangeSlabSize/Event.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -14,6 +15,8 @@ #include #include "DataStructures/DataBox/DataBox.hpp" +#include "Options/Context.hpp" +#include "Options/ParseError.hpp" #include "Options/String.hpp" #include "Parallel/GlobalCache.hpp" #include "Parallel/Reduction.hpp" @@ -21,7 +24,10 @@ #include "Time/ChangeSlabSize/Tags.hpp" #include "Time/StepChoosers/StepChooser.hpp" #include "Time/TimeStepId.hpp" +#include "Time/TimeStepRequest.hpp" +#include "Time/TimeStepRequestProcessor.hpp" #include "Utilities/Functional.hpp" +#include "Utilities/PrettyType.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" @@ -40,11 +46,12 @@ struct StoreNewSlabSize { static void apply(db::DataBox& box, Parallel::GlobalCache& /*cache*/, const ArrayIndex& /*array_index*/, - const int64_t slab_number, const double slab_size) { + const int64_t slab_number, + const TimeStepRequestProcessor& requests) { db::mutate<::Tags::ChangeSlabSize::NewSlabSize>( [&](const gsl::not_null< - std::map>*> - sizes) { (*sizes)[slab_number].insert(slab_size); }, + std::map>*> + sizes) { (*sizes)[slab_number].emplace_back(requests); }, make_not_null(&box)); } }; @@ -67,7 +74,7 @@ struct StoreNewSlabSize { class ChangeSlabSize : public Event { using ReductionData = Parallel::ReductionData< Parallel::ReductionDatum>, - Parallel::ReductionDatum>>; + Parallel::ReductionDatum>>; public: /// \cond @@ -97,8 +104,21 @@ class ChangeSlabSize : public Event { ChangeSlabSize() = default; ChangeSlabSize(std::vector>> step_choosers, - const uint64_t delay_change) - : step_choosers_(std::move(step_choosers)), delay_change_(delay_change) {} + const uint64_t delay_change, const Options::Context& context) + : step_choosers_(std::move(step_choosers)), delay_change_(delay_change) { + if (delay_change != 0) { + for (const auto& chooser : step_choosers_) { + if (not chooser->can_be_delayed()) { + // The runtime name might not be exactly the same as the one + // used by the factory, but hopefully it's close enough that + // the user can figure it out. + PARSE_ERROR(context, + "The " << pretty_type::get_runtime_type_name(*chooser) + << " StepChooser cannot be applied with a delay."); + } + } + } + } using compute_tags_for_observation_box = tmpl::list<>; @@ -121,16 +141,15 @@ class ChangeSlabSize : public Event { : time_step_id.slab_number() + 1; const auto slab_to_change = next_changable_slab + static_cast(delay_change_); + const auto slab_start = time_step_id.step_time(); + const auto current_slab_size = (time_step_id.time_runs_forward() ? 1 : -1) * + slab_start.slab().duration(); - double desired_slab_size = std::numeric_limits::infinity(); + TimeStepRequestProcessor step_requests(time_step_id.time_runs_forward()); bool synchronization_required = false; for (const auto& step_chooser : step_choosers_) { - desired_slab_size = std::min( - desired_slab_size, - step_chooser - ->desired_step(time_step_id.step_time().slab().duration().value(), - *box) - .first); + step_requests.process( + step_chooser->desired_step(current_slab_size.value(), *box).first); // We must synchronize if any step chooser requires it, not just // the limiting one, because choosers requiring synchronization // can be limiting on some processors and not others. @@ -151,15 +170,13 @@ class ChangeSlabSize : public Event { if (synchronization_required) { Parallel::contribute_to_reduction< ChangeSlabSize_detail::StoreNewSlabSize>( - ReductionData(slab_to_change, desired_slab_size), self_proxy, + ReductionData(slab_to_change, step_requests), self_proxy, component_proxy); } else { db::mutate<::Tags::ChangeSlabSize::NewSlabSize>( [&](const gsl::not_null< - std::map>*> - sizes) { - (*sizes)[slab_to_change].insert(desired_slab_size); - }, + std::map>*> + sizes) { (*sizes)[slab_to_change].push_back(step_requests); }, box); } } diff --git a/src/Time/ChangeSlabSize/Tags.hpp b/src/Time/ChangeSlabSize/Tags.hpp index aed2107ef59f..95cd157e543b 100644 --- a/src/Time/ChangeSlabSize/Tags.hpp +++ b/src/Time/ChangeSlabSize/Tags.hpp @@ -6,18 +6,25 @@ #include #include #include -#include +#include #include "DataStructures/DataBox/Tag.hpp" +#include "Time/TimeStepRequestProcessor.hpp" namespace Tags::ChangeSlabSize { /// Sizes requested for each slab by ChangeSlabSize events. struct NewSlabSize : db::SimpleTag { - using type = std::map>; + using type = std::map>; }; /// Number of ChangeSlabSize events changing the size at each slab. struct NumberOfExpectedMessages : db::SimpleTag { using type = std::map; }; + +/// Long-term desired slab size. Used as the default size if nothing +/// chooses a smaller one. +struct SlabSizeGoal : db::SimpleTag { + using type = double; +}; } // namespace Tags::ChangeSlabSize diff --git a/src/Time/StepChoosers/ByBlock.hpp b/src/Time/StepChoosers/ByBlock.hpp index e539445ee1f1..11c14aa3f066 100644 --- a/src/Time/StepChoosers/ByBlock.hpp +++ b/src/Time/StepChoosers/ByBlock.hpp @@ -3,6 +3,7 @@ #pragma once +#include #include #include #include // IWYU pragma: keep @@ -10,7 +11,8 @@ #include #include "Options/String.hpp" -#include "Time/StepChoosers/StepChooser.hpp" // IWYU pragma: keep +#include "Time/StepChoosers/StepChooser.hpp" +#include "Time/TimeStepRequest.hpp" #include "Utilities/ErrorHandling/Error.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" @@ -27,7 +29,7 @@ struct Element; /// \endcond namespace StepChoosers { -/// Suggests specified step sizes in each block +/// Sets a goal specified per-block. template class ByBlock : public StepChooser { public: @@ -44,25 +46,24 @@ class ByBlock : public StepChooser { "Step sizes, indexed by block number"}; }; - static constexpr Options::String help{ - "Suggests specified step sizes in each block"}; + static constexpr Options::String help{"Sets a goal specified per-block."}; using options = tmpl::list; explicit ByBlock(std::vector sizes) : sizes_(std::move(sizes)) {} using argument_tags = tmpl::list>; - std::pair operator()( - const Element& element, - const double /*last_step_magnitude*/) const { + std::pair operator()(const Element& element, + const double last_step) const { const size_t block = element.id().block_id(); if (block >= sizes_.size()) { ERROR("Step size not specified for block " << block); } - return std::make_pair(sizes_[block], true); + return {{.size_goal = std::copysign(sizes_[block], last_step)}, true}; } bool uses_local_data() const override { return true; } + bool can_be_delayed() const override { return true; } // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p) override { p | sizes_; } diff --git a/src/Time/StepChoosers/CMakeLists.txt b/src/Time/StepChoosers/CMakeLists.txt index 14281b061449..9558c65cbb9c 100644 --- a/src/Time/StepChoosers/CMakeLists.txt +++ b/src/Time/StepChoosers/CMakeLists.txt @@ -5,6 +5,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE Constant.cpp + Maximum.cpp Random.cpp StepToTimes.cpp ) @@ -19,7 +20,8 @@ spectre_target_headers( ElementSizeCfl.hpp ErrorControl.hpp Factory.hpp - Increase.hpp + LimitIncrease.hpp + Maximum.hpp PreventRapidIncrease.hpp Random.hpp StepChooser.hpp diff --git a/src/Time/StepChoosers/Cfl.hpp b/src/Time/StepChoosers/Cfl.hpp index 6c77b1e46114..40a0e1a44254 100644 --- a/src/Time/StepChoosers/Cfl.hpp +++ b/src/Time/StepChoosers/Cfl.hpp @@ -3,6 +3,7 @@ #pragma once +#include #include #include #include @@ -11,7 +12,8 @@ #include "DataStructures/DataBox/DataBox.hpp" #include "Domain/MinimumGridSpacing.hpp" #include "Options/String.hpp" -#include "Time/StepChoosers/StepChooser.hpp" // IWYU pragma: keep +#include "Time/StepChoosers/StepChooser.hpp" +#include "Time/TimeStepRequest.hpp" #include "Time/TimeSteppers/TimeStepper.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" @@ -31,7 +33,7 @@ struct MinimumGridSpacing; /// \endcond namespace StepChoosers { -/// Suggests a step size based on the CFL stability criterion. +/// Sets a goal based on the CFL stability criterion. template class Cfl : public StepChooser { public: @@ -49,7 +51,7 @@ class Cfl : public StepChooser { }; static constexpr Options::String help{ - "Suggests a step size based on the CFL stability criterion."}; + "Sets a goal based on the CFL stability criterion."}; using options = tmpl::list; explicit Cfl(const double safety_factor) : safety_factor_(safety_factor) {} @@ -63,19 +65,21 @@ class Cfl : public StepChooser { domain::Tags::MinimumGridSpacingCompute, typename System::compute_largest_characteristic_speed>; - std::pair operator()( - const double minimum_grid_spacing, - const TimeStepper& time_stepper, - const double speed, const double last_step_magnitude) const { + std::pair operator()(const double minimum_grid_spacing, + const TimeStepper& time_stepper, + const double speed, + const double last_step) const { const double time_stepper_stability_factor = time_stepper.stable_step(); const double step_size = safety_factor_ * time_stepper_stability_factor * minimum_grid_spacing / (speed * System::volume_dim); // Reject the step if the CFL condition is violated. - return std::make_pair(step_size, last_step_magnitude <= step_size); + return {{.size_goal = std::copysign(step_size, last_step)}, + abs(last_step) <= step_size}; } bool uses_local_data() const override { return true; } + bool can_be_delayed() const override { return true; } // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p) override { p | safety_factor_; } diff --git a/src/Time/StepChoosers/Constant.cpp b/src/Time/StepChoosers/Constant.cpp index 9b63ebe6da3f..2e9587cd8829 100644 --- a/src/Time/StepChoosers/Constant.cpp +++ b/src/Time/StepChoosers/Constant.cpp @@ -10,12 +10,7 @@ namespace Constant_detail { // This function lets us avoid including ParseOptions.hpp in the // header. double parse_options(const Options::Option& options) { - const auto value = options.parse_as(); - if (value <= 0.) { - PARSE_ERROR(options.context(), - "Requested step magnitude should be positive."); - } - return value; + return options.parse_as(); } } // namespace Constant_detail } // namespace StepChoosers diff --git a/src/Time/StepChoosers/Constant.hpp b/src/Time/StepChoosers/Constant.hpp index fe6b6977b499..73cd0e22a099 100644 --- a/src/Time/StepChoosers/Constant.hpp +++ b/src/Time/StepChoosers/Constant.hpp @@ -3,20 +3,21 @@ #pragma once +#include #include #include #include #include "Options/Options.hpp" #include "Options/String.hpp" -#include "Time/StepChoosers/StepChooser.hpp" // IWYU pragma: keep -#include "Utilities/ErrorHandling/Assert.hpp" +#include "Time/StepChoosers/StepChooser.hpp" +#include "Time/TimeStepRequest.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" namespace StepChoosers { -/// Suggests a constant step size. +/// Sets a constant goal. template class Constant : public StepChooser { public: @@ -27,20 +28,18 @@ class Constant : public StepChooser { WRAPPED_PUPable_decl_template(Constant); // NOLINT /// \endcond - static constexpr Options::String help{"Suggests a constant step size."}; + static constexpr Options::String help{"Sets a constant goal."}; - explicit Constant(const double value) : value_(value) { - ASSERT(value_ > 0., "Requested step magnitude should be positive."); - } + explicit Constant(const double value) : value_(value) {} using argument_tags = tmpl::list<>; - std::pair operator()( - const double /*last_step_magnitude*/) const { - return std::make_pair(value_, true); + std::pair operator()(const double last_step) const { + return {{.size_goal = std::copysign(value_, last_step)}, true}; } bool uses_local_data() const override { return false; } + bool can_be_delayed() const override { return true; } // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p) override { p | value_; } diff --git a/src/Time/StepChoosers/ElementSizeCfl.hpp b/src/Time/StepChoosers/ElementSizeCfl.hpp index ef78d3263500..d70ab4d8c143 100644 --- a/src/Time/StepChoosers/ElementSizeCfl.hpp +++ b/src/Time/StepChoosers/ElementSizeCfl.hpp @@ -3,6 +3,7 @@ #pragma once +#include #include #include #include @@ -12,6 +13,7 @@ #include "Domain/SizeOfElement.hpp" #include "Options/String.hpp" #include "Time/StepChoosers/StepChooser.hpp" +#include "Time/TimeStepRequest.hpp" #include "Time/TimeSteppers/TimeStepper.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" @@ -31,7 +33,7 @@ struct SizeOfElement; /// \endcond namespace StepChoosers { -/// Suggests a step size based on the CFL stability criterion, but uses the full +/// Sets a goal based on the CFL stability criterion, but uses the full /// size of the element as the length scale in question. /// /// This is useful as a coarse estimate for slabs, or to place a ceiling on @@ -53,7 +55,7 @@ class ElementSizeCfl : public StepChooser { }; static constexpr Options::String help{ - "Suggests a step size based on the CFL stability criterion, but in which " + "Sets a goal based on the CFL stability criterion, but in which " "the entire size of the element is used as the spacing in the " "computation. This is useful primarily for placing a ceiling on another " "dynamically-adjusted step chooser"}; @@ -70,10 +72,10 @@ class ElementSizeCfl : public StepChooser { tmpl::list, typename System::compute_largest_characteristic_speed>; - std::pair operator()( + std::pair operator()( const TimeStepper& time_stepper, const std::array& element_size, const double speed, - const double last_step_magnitude) const { + const double last_step) const { double min_size_of_element = std::numeric_limits::infinity(); for (auto face_to_face_dimension : element_size) { if (face_to_face_dimension < min_size_of_element) { @@ -84,10 +86,12 @@ class ElementSizeCfl : public StepChooser { const double step_size = safety_factor_ * time_stepper_stability_factor * min_size_of_element / (speed * Dim); // Reject the step if the CFL condition is violated. - return std::make_pair(step_size, last_step_magnitude <= step_size); + return {{.size_goal = std::copysign(step_size, last_step)}, + abs(last_step) <= step_size}; } bool uses_local_data() const override { return true; } + bool can_be_delayed() const override { return true; } // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p) override { p | safety_factor_; } diff --git a/src/Time/StepChoosers/ErrorControl.hpp b/src/Time/StepChoosers/ErrorControl.hpp index 3718fef7d61c..be7c7f7e80f3 100644 --- a/src/Time/StepChoosers/ErrorControl.hpp +++ b/src/Time/StepChoosers/ErrorControl.hpp @@ -23,6 +23,7 @@ #include "Time/Tags/IsUsingTimeSteppingErrorControl.hpp" #include "Time/Tags/StepperErrorTolerances.hpp" #include "Time/Tags/StepperErrors.hpp" +#include "Time/TimeStepRequest.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/Serialization/CharmPupable.hpp" @@ -54,8 +55,7 @@ struct ErrorControlBase : IsAnErrorControl { } // namespace ErrorControl_detail /*! - * \brief Suggests a step size based on a target absolute and/or relative error - * measure. + * \brief Sets a goal based on time-stepper truncation error. * * \details The suggested step is calculated via a simple specialization of the * scheme suggested in \cite Hairer1993. We first compute the aggregated error @@ -175,7 +175,7 @@ class ErrorControl }; static constexpr Options::String help{ - "Chooses a step based on a target relative and absolute error tolerance"}; + "Sets a goal based on time-stepper truncation error."}; using options = tmpl::list; @@ -199,16 +199,16 @@ class ErrorControl using argument_tags = tmpl::list<::Tags::StepperErrors>; - std::pair operator()( + std::pair operator()( const typename ::Tags::StepperErrors::type& errors, const double previous_step) const { - // request that the step size not be changed if there isn't a new error + // Do not request that the step size be changed if there isn't a new error // estimate if (not errors[1].has_value()) { - return std::make_pair(std::numeric_limits::infinity(), true); + return {{}, true}; } double new_step; - if (not errors[0].has_value()) { + if (not errors[0].has_value() or errors[0]->order != errors[1]->order) { new_step = errors[1]->step_size.value() * std::clamp(safety_factor_ * pow(1.0 / std::max(errors[1]->error, 1e-14), @@ -217,9 +217,6 @@ class ErrorControl } else { // From simple advice from Numerical Recipes 17.2.1 regarding a heuristic // for PI step control. - ASSERT(errors[0]->order == errors[1]->order, - "Error estimates are different orders. Need to figure our what " - "to do in this case."); const double alpha_factor = 0.7 / (errors[1]->order + 1); const double beta_factor = 0.4 / (errors[0]->order + 1); new_step = @@ -233,11 +230,12 @@ class ErrorControl // If the error is out-of-date we can still reject a step, but it // won't improve the reported error, so make sure to only do it // if we actually request a smaller step. - return {abs(new_step), - abs(new_step) >= previous_step or errors[1]->error <= 1.0}; + return {{.size_goal = new_step}, + abs(new_step) >= abs(previous_step) or errors[1]->error <= 1.0}; } bool uses_local_data() const override { return true; } + bool can_be_delayed() const override { return true; } StepperErrorTolerances tolerances() const override { return {.absolute = absolute_tolerance_, .relative = relative_tolerance_}; diff --git a/src/Time/StepChoosers/Factory.hpp b/src/Time/StepChoosers/Factory.hpp index ac5edcf1f314..25e2b979ef9a 100644 --- a/src/Time/StepChoosers/Factory.hpp +++ b/src/Time/StepChoosers/Factory.hpp @@ -9,7 +9,8 @@ #include "Time/StepChoosers/Constant.hpp" #include "Time/StepChoosers/ElementSizeCfl.hpp" #include "Time/StepChoosers/ErrorControl.hpp" -#include "Time/StepChoosers/Increase.hpp" +#include "Time/StepChoosers/LimitIncrease.hpp" +#include "Time/StepChoosers/Maximum.hpp" #include "Time/StepChoosers/PreventRapidIncrease.hpp" #include "Time/StepChoosers/StepToTimes.hpp" #include "Utilities/TMPL.hpp" @@ -30,7 +31,8 @@ using common_step_choosers = tmpl::push_back< StepChoosers::Cfl, StepChoosers::ElementSizeCfl>, tmpl::list<>>, - StepChoosers::Constant, StepChoosers::Increase>; + StepChoosers::Constant, StepChoosers::LimitIncrease, + StepChoosers::Maximum>; template using step_choosers_for_step_only = tmpl::list, diff --git a/src/Time/StepChoosers/Increase.hpp b/src/Time/StepChoosers/LimitIncrease.hpp similarity index 52% rename from src/Time/StepChoosers/Increase.hpp rename to src/Time/StepChoosers/LimitIncrease.hpp index 3ec1b4e8e6d5..8c7776f77ea7 100644 --- a/src/Time/StepChoosers/Increase.hpp +++ b/src/Time/StepChoosers/LimitIncrease.hpp @@ -9,40 +9,43 @@ #include #include "Options/String.hpp" -#include "Time/StepChoosers/StepChooser.hpp" // IWYU pragma: keep +#include "Time/StepChoosers/StepChooser.hpp" +#include "Time/TimeStepRequest.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" namespace StepChoosers { -/// Suggests increasing the step size by a constant ratio. +/// Limits step increase to a constant ratio. template -class Increase : public StepChooser { +class LimitIncrease : public StepChooser { public: /// \cond - Increase() = default; - explicit Increase(CkMigrateMessage* /*unused*/) {} + LimitIncrease() = default; + explicit LimitIncrease(CkMigrateMessage* /*unused*/) {} using PUP::able::register_constructor; - WRAPPED_PUPable_decl_template(Increase); // NOLINT + WRAPPED_PUPable_decl_template(LimitIncrease); // NOLINT /// \endcond struct Factor { using type = double; - static constexpr Options::String help{"Factor to increase by"}; + static constexpr Options::String help{"Factor to allow increase by"}; static type lower_bound() { return 1.0; } }; - static constexpr Options::String help{"Suggests a constant factor increase."}; + static constexpr Options::String help{ + "Limits step increase to a constant ratio."}; using options = tmpl::list; - explicit Increase(const double factor) : factor_(factor) {} + explicit LimitIncrease(const double factor) : factor_(factor) {} using argument_tags = tmpl::list<>; - std::pair operator()(const double last_step_magnitude) const { - return std::make_pair(last_step_magnitude * factor_, true); + std::pair operator()(const double last_step) const { + return {{.size = last_step * factor_}, true}; } bool uses_local_data() const override { return false; } + bool can_be_delayed() const override { return true; } // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p) override { p | factor_; } @@ -53,6 +56,6 @@ class Increase : public StepChooser { /// \cond template -PUP::able::PUP_ID Increase::my_PUP_ID = 0; // NOLINT +PUP::able::PUP_ID LimitIncrease::my_PUP_ID = 0; // NOLINT /// \endcond } // namespace StepChoosers diff --git a/src/Time/StepChoosers/Maximum.cpp b/src/Time/StepChoosers/Maximum.cpp new file mode 100644 index 000000000000..e799bb8c6a6c --- /dev/null +++ b/src/Time/StepChoosers/Maximum.cpp @@ -0,0 +1,14 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Time/StepChoosers/Maximum.hpp" + +#include "Options/ParseOptions.hpp" + +namespace StepChoosers::Maximum_detail { +// This function lets us avoid including ParseOptions.hpp in the +// header. +double parse_options(const Options::Option& options) { + return options.parse_as(); +} +} // namespace StepChoosers::Maximum_detail diff --git a/src/Time/StepChoosers/Maximum.hpp b/src/Time/StepChoosers/Maximum.hpp new file mode 100644 index 000000000000..ee63bf131bc0 --- /dev/null +++ b/src/Time/StepChoosers/Maximum.hpp @@ -0,0 +1,69 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include + +#include "Options/Options.hpp" +#include "Options/String.hpp" +#include "Time/StepChoosers/StepChooser.hpp" +#include "Time/TimeStepRequest.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/TMPL.hpp" + +namespace StepChoosers { + +/// Limits the step size to a constant. +template +class Maximum : public StepChooser { + public: + /// \cond + Maximum() = default; + explicit Maximum(CkMigrateMessage* /*unused*/) {} + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(Maximum); // NOLINT + /// \endcond + + static constexpr Options::String help{"Limits the step size to a constant."}; + + explicit Maximum(const double value) : value_(value) {} + + using argument_tags = tmpl::list<>; + + std::pair operator()(const double last_step) const { + return {{.size = std::copysign(value_, last_step)}, true}; + } + + bool uses_local_data() const override { return false; } + bool can_be_delayed() const override { return true; } + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& p) override { p | value_; } + + private: + double value_ = std::numeric_limits::signaling_NaN(); +}; + +/// \cond +template +PUP::able::PUP_ID Maximum::my_PUP_ID = 0; // NOLINT +/// \endcond + +namespace Maximum_detail { +double parse_options(const Options::Option& options); +} // namespace Maximum_detail +} // namespace StepChoosers + +template +struct Options::create_from_yaml> { + template + static StepChoosers::Maximum create( + const Options::Option& options) { + return StepChoosers::Maximum( + StepChoosers::Maximum_detail::parse_options(options)); + } +}; diff --git a/src/Time/StepChoosers/PreventRapidIncrease.hpp b/src/Time/StepChoosers/PreventRapidIncrease.hpp index 86c8fb28b9f7..824f8670a558 100644 --- a/src/Time/StepChoosers/PreventRapidIncrease.hpp +++ b/src/Time/StepChoosers/PreventRapidIncrease.hpp @@ -4,23 +4,25 @@ #pragma once #include // IWYU pragma: keep // for abs -#include #include #include #include "Options/String.hpp" #include "Time/StepChoosers/StepChooser.hpp" // IWYU pragma: keep #include "Time/Tags/HistoryEvolvedVariables.hpp" +#include "Time/TimeStepRequest.hpp" #include "Time/Utilities.hpp" #include "Utilities/Serialization/CharmPupable.hpp" #include "Utilities/TMPL.hpp" namespace StepChoosers { +/// Limits the time step to prevent multistep integrator instabilities. +/// /// Avoids instabilities due to rapid increases in the step size by -/// preventing the step size from increasing unless all steps in the -/// time-stepper history are the same size. If there have been recent -/// step size changes the new size bound is the size of the most -/// recent step, otherwise it is infinite (no restriction is imposed). +/// preventing the step size from increasing if any step in the +/// time-stepper history increased. If there have been recent step +/// size increases, the new size bound is the size of the most recent +/// step, otherwise no restriction is imposed. template class PreventRapidIncrease : public StepChooser { public: @@ -32,40 +34,40 @@ class PreventRapidIncrease : public StepChooser { /// \endcond static constexpr Options::String help{ - "Prevents rapid increases in time step that can cause integrator \n" - "instabilities."}; + "Limits the time step to prevent multistep integrator instabilities."}; using options = tmpl::list<>; using argument_tags = tmpl::list<::Tags::HistoryEvolvedVariables<>>; template - std::pair operator()(const History& history, - const double last_step_magnitude) const { + std::pair operator()(const History& history, + const double last_step) const { if (history.size() < 2) { - return std::make_pair(std::numeric_limits::infinity(), true); + return {{}, true}; } const double sloppiness = slab_rounding_error(history.front().time_step_id.step_time()); std::optional