diff --git a/src/Evolution/Actions/RunEventsAndDenseTriggers.hpp b/src/Evolution/Actions/RunEventsAndDenseTriggers.hpp index cdc300647cc1..35f04e89d492 100644 --- a/src/Evolution/Actions/RunEventsAndDenseTriggers.hpp +++ b/src/Evolution/Actions/RunEventsAndDenseTriggers.hpp @@ -207,13 +207,8 @@ struct RunEventsAndDenseTriggers { return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } - // This can only be true the first time through the loop, - // because triggers are not allowed to reschedule for the time - // they just triggered at. This check is primarily to avoid - // special-case bookkeeping for the initial simulation time. - const bool already_at_correct_time = - db::get<::Tags::Time>(box) == next_trigger; - if (not already_at_correct_time) { + // Avoid invalidating compute items unless necessary. + if (db::get<::Tags::Time>(box) != next_trigger) { time_restorer.save(); db::mutate<::Tags::Time>( make_not_null(&box), @@ -228,53 +223,46 @@ struct RunEventsAndDenseTriggers { switch (triggered) { case TriggeringState::NotReady: return {Parallel::AlgorithmExecution::Retry, std::nullopt}; - case TriggeringState::NeedsEvolvedVariables: - if (not already_at_correct_time) { - bool ready = true; - tmpl::for_each([&](auto postprocessor_v) { - using postprocessor = tmpl::type_from; - if (ready) { - if (not postprocessor::is_ready(make_not_null(&box), - make_not_null(&inboxes), cache, - array_index, component)) { - ready = false; - } + case TriggeringState::NeedsEvolvedVariables: { + bool ready = true; + tmpl::for_each([&](auto postprocessor_v) { + using postprocessor = tmpl::type_from; + if (ready) { + if (not postprocessor::is_ready(make_not_null(&box), + make_not_null(&inboxes), cache, + array_index, component)) { + ready = false; } - }); - if (not ready) { - return {Parallel::AlgorithmExecution::Retry, std::nullopt}; - } - - using history_tag = ::Tags::HistoryEvolvedVariables; - bool dense_output_succeeded = false; - // Restore the evolved variables if they have been - // previously altered so the time stepper starts in the - // correct state, then save them if they have not been - // previously altered. (So one of these two lines is a - // no-op.) - variables_restorer.restore(); - variables_restorer.save(); - db::mutate( - make_not_null(&box), - [&dense_output_succeeded, &next_trigger]( - gsl::not_null vars, - const TimeStepper& stepper, - const typename history_tag::type& history) { - dense_output_succeeded = - stepper.dense_update_u(vars, history, next_trigger); - }, - db::get<::Tags::TimeStepper<>>(box), db::get(box)); - if (not dense_output_succeeded) { - // Need to take another time step - return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } + }); + if (not ready) { + return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + } - postprocessor_restorer.save(); - tmpl::for_each([&box](auto postprocessor_v) { - using postprocessor = tmpl::type_from; - db::mutate_apply(make_not_null(&box)); - }); + using history_tag = ::Tags::HistoryEvolvedVariables; + bool dense_output_succeeded = false; + variables_restorer.save(); + db::mutate( + make_not_null(&box), + [&dense_output_succeeded, &next_trigger]( + gsl::not_null vars, + const TimeStepper& stepper, + const typename history_tag::type& history) { + dense_output_succeeded = + stepper.dense_update_u(vars, history, next_trigger); + }, + db::get<::Tags::TimeStepper<>>(box), db::get(box)); + if (not dense_output_succeeded) { + // Need to take another time step + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } + + postprocessor_restorer.save(); + tmpl::for_each([&box](auto postprocessor_v) { + using postprocessor = tmpl::type_from; + db::mutate_apply(make_not_null(&box)); + }); + } [[fallthrough]]; default: break; diff --git a/src/Time/TimeSteppers/AdamsBashforth.cpp b/src/Time/TimeSteppers/AdamsBashforth.cpp index bc0b2c6ac606..21b19bc2b90c 100644 --- a/src/Time/TimeSteppers/AdamsBashforth.cpp +++ b/src/Time/TimeSteppers/AdamsBashforth.cpp @@ -198,7 +198,9 @@ void clean_history(const MutableUntypedHistory& history) { while (history.size() > history.integration_order()) { history.pop_front(); } - history.discard_value(history.back().time_step_id); + if (history.size() > 1) { + history.discard_value(history[history.size() - 2].time_step_id); + } } } // namespace @@ -257,7 +259,6 @@ bool AdamsBashforth::update_u_impl( const gsl::not_null u, const gsl::not_null u_error, const MutableUntypedHistory& history, const TimeDelta& time_step) const { clean_history(history); - *u_error = *u; update_u_common(u, history, time_step, history.integration_order()); // the error estimate is only useful once the history has enough elements to // do more than one order of step @@ -294,6 +295,7 @@ void AdamsBashforth::update_u_common(const gsl::not_null u, get_coefficients(history_time_iterator(history_start), history_time_iterator(history.end()), time_step); + *u = *history.back().value; auto coefficient = coefficients.rbegin(); for (auto history_entry = history_start; history_entry != history.end(); @@ -364,6 +366,11 @@ void AdamsBashforth::boundary_dense_output_impl( const gsl::not_null result, const TimeSteppers::BoundaryHistoryEvaluator& coupling, const double time) const { + if ((coupling.local_end() - 1)->value() == time) { + // Nothing to do. The requested time is the start of the step, + // which is the input value of `result`. + return; + } return boundary_impl(result, coupling, ApproximateTime{time}); } diff --git a/src/Time/TimeSteppers/RungeKutta.cpp b/src/Time/TimeSteppers/RungeKutta.cpp index a9b4b47ccca9..2695c70e1ab8 100644 --- a/src/Time/TimeSteppers/RungeKutta.cpp +++ b/src/Time/TimeSteppers/RungeKutta.cpp @@ -110,9 +110,11 @@ void update_u_impl_with_tableau(const gsl::not_null u, if (history.size() > 1) { history.pop_front(); } - history.discard_value(history.back().time_step_id); + } else if (history.substeps().size() == 1) { + history.discard_value(history.front().time_step_id); } else { - history.discard_value(history.substeps().back().time_step_id); + history.discard_value( + history.substeps()[history.substeps().size() - 2].time_step_id); } ASSERT(history.size() == 1, "Have more than one step after cleanup."); @@ -125,12 +127,15 @@ void update_u_impl_with_tableau(const gsl::not_null u, if (substep == 0) { ASSERT(tableau.substep_coefficients[0].size() == 1, "First substep should use one derivative."); - *u += tableau.substep_coefficients[0][0] * dt * history.front().derivative; + *u = *history.front().value + + tableau.substep_coefficients[0][0] * dt * history.front().derivative; } else if (substep == number_of_substeps - 1) { + *u = *history.substeps().back().value; update_between_substeps(u, history, dt, tableau.substep_coefficients[substep - 1], tableau.result_coefficients); } else if (substep < number_of_substeps - 1) { + *u = *history.substeps().back().value; update_between_substeps(u, history, dt, tableau.substep_coefficients[substep - 1], tableau.substep_coefficients[substep]); @@ -192,6 +197,7 @@ bool RungeKutta::dense_update_u_impl(const gsl::not_null u, if (time == step_end) { // Special case necessary for dense output at the initial time, // before taking a step. + *u = *history.back().value; return true; } const evolution_less before{step_end > step_start}; @@ -206,6 +212,7 @@ bool RungeKutta::dense_update_u_impl(const gsl::not_null u, const auto& tableau = butcher_tableau(); + *u = *history.back().value; // The Butcher dense output coefficients are given in terms of the // start of the step, but we are passed the value at the end of the // step, so we have to undo the step first. diff --git a/src/Time/TimeSteppers/TimeStepper.hpp b/src/Time/TimeSteppers/TimeStepper.hpp index 736bb8b94f8f..33e310ad8a4d 100644 --- a/src/Time/TimeSteppers/TimeStepper.hpp +++ b/src/Time/TimeSteppers/TimeStepper.hpp @@ -72,7 +72,7 @@ class TimeStepper : public PUP::able { #undef TIME_STEPPER_DECLARE_VIRTUALS_IMPL /// \endcond - /// Add the change for the current substep to u. + /// Set \p u to the value at the end of the current substep. /// /// Derived classes must implement this as a function with signature /// @@ -90,8 +90,8 @@ class TimeStepper : public PUP::able { time_step); } - /// Add the change for the current substep to u; report the error measure when - /// available. + /// Set \p u to the value at the end of the current substep; report the error + /// measure when available. /// /// For a substep method, the error measure will only be available on full /// steps. For a multistep method, the error measure will only be available diff --git a/tests/Unit/Evolution/Actions/Test_RunEventsAndDenseTriggers.cpp b/tests/Unit/Evolution/Actions/Test_RunEventsAndDenseTriggers.cpp index 15f6a4297001..ca55758f43fa 100644 --- a/tests/Unit/Evolution/Actions/Test_RunEventsAndDenseTriggers.cpp +++ b/tests/Unit/Evolution/Actions/Test_RunEventsAndDenseTriggers.cpp @@ -110,11 +110,13 @@ class TestTrigger : public DenseTrigger { Parallel::GlobalCache& /*cache*/, const ArrayIndex& /*array_index*/, const Component* /*component*/, const double time) const { - if (time == init_time_) { + if (time == trigger_time_) { + return is_triggered_; + } else { + // First initialization call. + CHECK(time == init_time_); return false; } - CHECK(time == trigger_time_); - return is_triggered_; } using next_check_time_argument_tags = tmpl::list; @@ -123,11 +125,13 @@ class TestTrigger : public DenseTrigger { Parallel::GlobalCache& /*cache*/, const ArrayIndex& /*array_index*/, const Component* /*component*/, const double time) const { - if (time == init_time_) { + if (time == trigger_time_) { + return next_trigger_; + } else { + // First initialization call. + CHECK(time == init_time_); return trigger_time_; } - CHECK(time == trigger_time_); - return next_trigger_; } // NOLINTNEXTLINE(google-runtime-references) @@ -323,11 +327,12 @@ void test(const bool time_runs_forward) { const double done_time = (time_runs_forward ? 1.0 : -1.0) * std::numeric_limits::infinity(); const VarsType initial_vars{1, 8.0}; + const VarsType stored_vars{1, 123.0}; const DtVarsType deriv_vars{1, 1.0}; const VarsType center_vars = initial_vars + 0.5 * step_size * deriv_vars; const auto set_up_component = - [&deriv_vars, &exact_step_size, &initial_vars, &start_time, + [&deriv_vars, &exact_step_size, &initial_vars, &start_time, &stored_vars, &time_step_id]( const gsl::not_null runner, const std::vector, @@ -352,7 +357,7 @@ void test(const bool time_runs_forward) { ActionTesting::emplace_array_component( runner, ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0}, 0, time_step_id, exact_step_size, start_time, - std::optional{}, initial_vars, std::move(history), + std::optional{}, stored_vars, std::move(history), evolution::EventsAndDenseTriggers( std::move(events_and_dense_triggers)), typename evolution::dg::Tags::NeighborMesh<1>::type{}, @@ -427,11 +432,24 @@ void test(const bool time_runs_forward) { reschedule ? std::optional{done_time} : std::nullopt; set_up_component(&runner, {{step_center, true, next_check, false}}); CHECK(run_if_ready(make_not_null(&runner)) == reschedule); - TestEvent::check_calls({{step_center, initial_vars}}); + TestEvent::check_calls({{step_center, stored_vars}}); }; check_not_needed(true); check_not_needed(false); + // Variables not needed at initial time + const auto check_not_needed_initial = [&](const bool reschedule) { + MockRuntimeSystem runner{ + {std::make_unique(1)}}; + const auto next_check = + reschedule ? std::optional{done_time} : std::nullopt; + set_up_component(&runner, {{start_time, true, next_check, false}}); + CHECK(run_if_ready(make_not_null(&runner)) == reschedule); + TestEvent::check_calls({{start_time, stored_vars}}); + }; + check_not_needed_initial(true); + check_not_needed_initial(false); + // Variables needed const auto check_needed = [&](const bool reschedule) { MockRuntimeSystem runner{ @@ -444,6 +462,18 @@ void test(const bool time_runs_forward) { check_needed(true); check_needed(false); + // Variables needed at initial time + const auto check_needed_initial = [&](const bool reschedule) { + MockRuntimeSystem runner{ + {std::make_unique(1)}}; + const auto next_check = + reschedule ? std::optional{done_time} : std::nullopt; + set_up_component(&runner, {{start_time, true, next_check, true}}); + TestCase::check_dense(&runner, reschedule, {{start_time, initial_vars}}); + }; + check_needed_initial(true); + check_needed_initial(false); + // Missing dense output data const auto check_missing_dense_data = [&](const bool data_needed) { MockRuntimeSystem runner{ @@ -473,7 +503,7 @@ void test(const bool time_runs_forward) { } else { // If we don't need the data, it shouldn't matter whether it is missing. CHECK(run_if_ready(make_not_null(&runner))); - TestEvent::check_calls({{step_center, initial_vars}}); + TestEvent::check_calls({{step_center, stored_vars}}); } }; check_missing_dense_data(true); diff --git a/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.cpp b/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.cpp index b19bc3be3cfd..07d81b0c395d 100644 --- a/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.cpp +++ b/tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.cpp @@ -42,6 +42,7 @@ void take_step( ++substep) { CHECK(time_id.substep() == substep); history->insert(time_id, *y, rhs(*y, time_id.substep_time().value())); + *y = std::numeric_limits::signaling_NaN(); stepper.update_u(y, history, step_size); time_id = stepper.next_time_id(time_id, step_size); } @@ -60,6 +61,7 @@ void take_step_and_check_error( ++substep) { CHECK(time_id.substep() == substep); history->insert(time_id, *y, rhs(*y, time_id.substep_time().value())); + *y = std::numeric_limits::signaling_NaN(); bool error_updated = stepper.update_u(y, y_error, history, step_size); CAPTURE(substep); REQUIRE((substep == stepper.number_of_substeps_for_error() - 1) == @@ -454,6 +456,7 @@ void check_dense_output(const TimeStepper& stepper, auto step = step_size; for (;;) { history.insert(time_id, y, y); + y = std::numeric_limits::signaling_NaN(); if (not before((time_id.step_time() + step).value(), time)) { if (stepper.dense_update_u(make_not_null(&y), history, time)) { return y;