Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <memory>
       8             : #include <optional>
       9             : 
      10             : #include "ControlSystem/CombinedName.hpp"
      11             : #include "ControlSystem/Metafunctions.hpp"
      12             : #include "DataStructures/DataBox/DataBox.hpp"
      13             : #include "Parallel/AlgorithmExecution.hpp"
      14             : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
      15             : #include "Parallel/ArrayCollection/PerformAlgorithmOnElement.hpp"
      16             : #include "Parallel/ArrayCollection/Tags/ElementLocations.hpp"
      17             : #include "Parallel/ArrayComponentId.hpp"
      18             : #include "Parallel/Callback.hpp"
      19             : #include "Parallel/GlobalCache.hpp"
      20             : #include "ParallelAlgorithms/Actions/GetItemFromDistributedObject.hpp"
      21             : #include "Time/ChangeSlabSize/ChangeSlabSize.hpp"
      22             : #include "Time/Tags/HistoryEvolvedVariables.hpp"
      23             : #include "Time/Tags/MinimumTimeStep.hpp"
      24             : #include "Utilities/ErrorHandling/Assert.hpp"
      25             : #include "Utilities/ErrorHandling/Error.hpp"
      26             : #include "Utilities/Gsl.hpp"
      27             : #include "Utilities/TMPL.hpp"
      28             : 
      29             : /// \cond
      30             : namespace Tags {
      31             : struct TimeStep;
      32             : struct TimeStepId;
      33             : template <typename StepperInterface>
      34             : struct TimeStepper;
      35             : }  // namespace Tags
      36             : class TimeStepper;
      37             : namespace control_system::Tags {
      38             : template <typename ControlSystems>
      39             : struct FutureMeasurements;
      40             : struct MeasurementTimescales;
      41             : }  // namespace control_system::Tags
      42             : namespace domain::Tags {
      43             : struct FunctionsOfTime;
      44             : }  // namespace domain::Tags
      45             : namespace tuples {
      46             : template <typename... Tags>
      47             : class TaggedTuple;
      48             : }  // namespace tuples
      49             : /// \endcond
      50             : 
      51             : namespace control_system::Actions {
      52             : /// \ingroup ControlSystemGroup
      53             : /// \brief Limit the step size in a GTS evolution to prevent deadlocks from
      54             : /// control system measurements.
      55             : ///
      56             : /// \details Most time steppers require evaluations of the coordinates
      57             : /// at several times during the step before they can produce dense
      58             : /// output.  If any of those evaluations require a function-of-time
      59             : /// update depending on a measurement within the step, the evolution
      60             : /// will deadlock.  This action reduces the step size if necessary to
      61             : /// prevent that from happening.
      62             : ///
      63             : /// Specifically:
      64             : /// 1. The chosen step will never be longer than the unmodified step,
      65             : ///    and will be short enough to avoid relevant function-of-time
      66             : ///    expirations.
      67             : /// 2. Given the previous, the step will cover as many control-system
      68             : ///    updates as possible.
      69             : /// 3. If the next step is likely to be limited by this action, adjust
      70             : ///    the length of the current step so that this step and the next
      71             : ///    step will be as close as possible to the same size.
      72             : template <typename ControlSystems>
      73           1 : struct LimitTimeStep {
      74             :  private:
      75           0 :   using control_system_groups =
      76             :       tmpl::transform<metafunctions::measurements_t<ControlSystems>,
      77             :                       metafunctions::control_systems_with_measurement<
      78             :                           tmpl::pin<ControlSystems>, tmpl::_1>>;
      79             : 
      80             :   template <typename Group>
      81           0 :   struct GroupExpiration {
      82           0 :     using type = double;
      83             :   };
      84             : 
      85             :  public:
      86           0 :   using const_global_cache_tags = tmpl::list<::Tags::MinimumTimeStep>;
      87             : 
      88             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
      89             :             size_t Dim, typename ActionList, typename ParallelComponent>
      90           0 :   static Parallel::iterable_action_return_t apply(
      91             :       db::DataBox<DbTagsList>& box,
      92             :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
      93             :       Parallel::GlobalCache<Metavariables>& cache,
      94             :       const ElementId<Dim>& array_index, ActionList /*meta*/,
      95             :       const ParallelComponent* const /*meta*/) {
      96             :     static_assert(not Metavariables::local_time_stepping,
      97             :                   "The control system LimitTimeStep action is only for global "
      98             :                   "time stepping.");
      99             :     const auto& time_step_id = db::get<::Tags::TimeStepId>(box);
     100             :     if (time_step_id.substep() != 0) {
     101             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     102             :     }
     103             : 
     104             :     const auto& time_stepper = db::get<::Tags::TimeStepper<TimeStepper>>(box);
     105             :     if (time_stepper.monotonic()) {
     106             :       // Monotonic steppers order operations in the same manner at the
     107             :       // control system, so they cannot introduce deadlocks.
     108             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     109             :     }
     110             : 
     111             :     auto& proxy = ::Parallel::get_parallel_component<ParallelComponent>(cache);
     112             : 
     113             :     // Minimum expiration time for any FoT in the measurement group.
     114             :     tmpl::wrap<tmpl::transform<control_system_groups,
     115             :                                tmpl::bind<GroupExpiration, tmpl::_1>>,
     116             :                tuples::TaggedTuple>
     117             :         group_expiration_times{};
     118             : 
     119             :     bool ready = true;
     120             :     // Calculate group_expiration_times
     121             :     tmpl::for_each<control_system_groups>([&](auto group_v) {
     122             :       if (not ready) {
     123             :         return;
     124             :       }
     125             :       using group = tmpl::type_from<decltype(group_v)>;
     126             : 
     127             :       auto& future_measurements =
     128             :           db::get_mutable_reference<Tags::FutureMeasurements<group>>(
     129             :               make_not_null(&box));
     130             : 
     131             :       std::optional<double> group_update = future_measurements.next_update();
     132             :       if (not group_update.has_value()) {
     133             :         Parallel::mutable_cache_item_is_ready<
     134             :             control_system::Tags::MeasurementTimescales>(
     135             :             cache,
     136             :             Parallel::make_array_component_id<ParallelComponent>(array_index),
     137             :             [&](const auto& measurement_timescales) {
     138             :               const auto& group_timescale =
     139             :                   *measurement_timescales.at(combined_name<group>());
     140             :               future_measurements.update(group_timescale);
     141             :               group_update = future_measurements.next_update();
     142             :               ready = group_update.has_value();
     143             :               if constexpr (Parallel::is_dg_element_collection_v<
     144             :                                 ParallelComponent>) {
     145             :                 // Note: The ArrayComponentId is still created with the
     146             :                 // array_index (ElementId) because we only support 1 callback
     147             :                 // per ArrayComponentId. This would mean for nodegroups we would
     148             :                 // discard a lot of callbacks that we need. Alternatively, the
     149             :                 // callback could do a broadcast to all elements on the
     150             :                 // nodegroup.
     151             :                 const auto element_location = static_cast<int>(
     152             :                     Parallel::local_synchronous_action<
     153             :                         ::Parallel::Actions::GetItemFromDistributedOject<
     154             :                             Parallel::Tags::ElementLocations<Dim>>>(proxy)
     155             :                         ->at(array_index));
     156             :                 return ready ? std::unique_ptr<Parallel::Callback>{}
     157             :                              : std::unique_ptr<Parallel::Callback>(
     158             :                                    new Parallel::ThreadedActionCallback<
     159             :                                        Parallel::Actions::
     160             :                                            PerformAlgorithmOnElement<false>,
     161             :                                        decltype(proxy[element_location]),
     162             :                                        std::decay_t<decltype(array_index)>>{
     163             :                                        proxy[element_location], array_index});
     164             :               } else {
     165             :                 return ready ? std::unique_ptr<Parallel::Callback>{}
     166             :                              : std::unique_ptr<Parallel::Callback>(
     167             :                                    new Parallel::PerformAlgorithmCallback(
     168             :                                        proxy[array_index]));
     169             :               }
     170             :             });
     171             :         if (not ready) {
     172             :           return;
     173             :         }
     174             :       }
     175             : 
     176             :       auto& group_expiration =
     177             :           get<GroupExpiration<group>>(group_expiration_times);
     178             :       group_expiration = std::numeric_limits<double>::infinity();
     179             : 
     180             :       if (*group_update == std::numeric_limits<double>::infinity()) {
     181             :         // Control measurement is not active.
     182             :         return;
     183             :       }
     184             : 
     185             :       // Calculate group_expiration
     186             :       Parallel::mutable_cache_item_is_ready<domain::Tags::FunctionsOfTime>(
     187             :           cache,
     188             :           Parallel::make_array_component_id<ParallelComponent>(array_index),
     189             :           [&](const auto& functions_of_time) {
     190             :             tmpl::for_each<group>([&](auto system) {
     191             :               using System = tmpl::type_from<decltype(system)>;
     192             :               if (not ready) {
     193             :                 return;
     194             :               }
     195             :               const auto& fot = *functions_of_time.at(System::name());
     196             :               ready = fot.time_bounds()[1] > *group_update;
     197             :               if (ready) {
     198             :                 group_expiration = std::min(
     199             :                     group_expiration, fot.expiration_after(*group_update));
     200             :               }
     201             :             });
     202             :             if constexpr (Parallel::is_dg_element_collection_v<
     203             :                               ParallelComponent>) {
     204             :               // Note: The ArrayComponentId is still created with the
     205             :               // array_index (ElementId) because we only support 1 callback
     206             :               // per ArrayComponentId. This would mean for nodegroups we would
     207             :               // discard a lot of callbacks that we need. Alternatively, the
     208             :               // callback could do a broadcast to all elements on the
     209             :               // nodegroup.
     210             :               const auto element_location = static_cast<int>(
     211             :                   Parallel::local_synchronous_action<
     212             :                       ::Parallel::Actions::GetItemFromDistributedOject<
     213             :                           Parallel::Tags::ElementLocations<Dim>>>(proxy)
     214             :                       ->at(array_index));
     215             :               return ready
     216             :                          ? std::unique_ptr<Parallel::Callback>{}
     217             :                          : std::unique_ptr<Parallel::Callback>(
     218             :                                new Parallel::ThreadedActionCallback<
     219             :                                    Parallel::Actions::PerformAlgorithmOnElement<
     220             :                                        false>,
     221             :                                    decltype(proxy[element_location]),
     222             :                                    std::decay_t<decltype(array_index)>>{
     223             :                                    proxy[element_location], array_index});
     224             :             } else {
     225             :               return ready ? std::unique_ptr<Parallel::Callback>{}
     226             :                            : std::unique_ptr<Parallel::Callback>(
     227             :                                  new Parallel::PerformAlgorithmCallback(
     228             :                                      proxy[array_index]));
     229             :             }
     230             :           });
     231             :     });
     232             :     if (not ready) {
     233             :       return {Parallel::AlgorithmExecution::Retry, std::nullopt};
     234             :     }
     235             : 
     236             :     const double orig_step_start = time_step_id.step_time().value();
     237             :     const double orig_step_end =
     238             :         (time_step_id.step_time() + db::get<::Tags::TimeStep>(box)).value();
     239             : 
     240             :     // Smallest of the current step end and the FoT expirations.  We
     241             :     // can't step any farther than this.
     242             :     const double latest_valid_step =
     243             :         tmpl::as_pack<control_system_groups>([&](auto... groups) {
     244             :           return std::min(
     245             :               {orig_step_end,
     246             :                get<GroupExpiration<tmpl::type_from<decltype(groups)>>>(
     247             :                    group_expiration_times)...});
     248             :         });
     249             : 
     250             :     if (not tmpl::as_pack<::Tags::get_all_history_tags<DbTagsList>>(
     251             :             [&]<typename... HistoryTags>(tmpl::type_<HistoryTags>... /*meta*/) {
     252             :               return (... and time_stepper.can_change_step_size(
     253             :                                   time_step_id, db::get<HistoryTags>(box)));
     254             :             })) {
     255             :       if (orig_step_end > latest_valid_step) {
     256             :         ERROR(
     257             :             "Step must be decreased to avoid control-system deadlock, but "
     258             :             "time-stepper requires a fixed step size.");
     259             :       }
     260             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     261             :     }
     262             :     ASSERT(db::get<::Tags::TimeStep>(box).fraction() == 1,
     263             :            "Trying to change GTS step, but it isn't a full slab.  Non-slab "
     264             :            "steps should only happen during self-start, but the preceding "
     265             :            "check should have ended the action if this is self-start.");
     266             : 
     267             :     // The last update that we can perform on the next step.  Don't
     268             :     // shrink the step past this time since that will force another
     269             :     // step to take the measurement.
     270             :     double last_update_time = orig_step_start;
     271             :     // Step time that produces a balanced step with the following
     272             :     // step, ignoring the restrictions on the current step.
     273             :     double preferred_step_time = orig_step_end;
     274             : 
     275             :     tmpl::for_each<control_system_groups>([&](auto group_v) {
     276             :       using group = tmpl::type_from<decltype(group_v)>;
     277             : 
     278             :       // This was used above, so it is not nullopt.
     279             :       const double group_update =
     280             :           db::get<Tags::FutureMeasurements<group>>(box).next_update().value();
     281             :       if (group_update <= latest_valid_step) {
     282             :         // We've satisfied this measurement.
     283             :         last_update_time = std::max(last_update_time, group_update);
     284             :       } else {
     285             :         // We can't make it far enough to do the final measurement.
     286             :         // Try to avoid a small step by choosing two equal-sized steps
     287             :         // to the expiration time.
     288             :         const double equal_step_time =
     289             :             0.5 * (orig_step_start +
     290             :                    get<GroupExpiration<group>>(group_expiration_times));
     291             :         preferred_step_time = std::min(preferred_step_time, equal_step_time);
     292             :       }
     293             :     });
     294             : 
     295             :     const double new_step_end =
     296             :         std::clamp(preferred_step_time, last_update_time, latest_valid_step);
     297             : 
     298             :     change_slab_size(make_not_null(&box), new_step_end);
     299             : 
     300             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     301             :   }
     302             : };
     303             : }  // namespace control_system::Actions