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