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 <string>
       9             : #include <type_traits>
      10             : #include <unordered_map>
      11             : #include <unordered_set>
      12             : 
      13             : #include "ControlSystem/CalculateMeasurementTimescales.hpp"
      14             : #include "ControlSystem/CombinedName.hpp"
      15             : #include "ControlSystem/Metafunctions.hpp"
      16             : #include "ControlSystem/Tags/SystemTags.hpp"
      17             : #include "DataStructures/DataVector.hpp"
      18             : #include "Domain/Creators/DomainCreator.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : 
      21             : namespace control_system {
      22             : /*!
      23             :  * \ingroup ControlSystemGroup
      24             :  * \brief Calculate the next expiration time for the FunctionsOfTime.
      25             :  *
      26             :  * If \p delay_update is false, this returns
      27             :  *
      28             :  * \f{equation}
      29             :  * T_\mathrm{expr}^\mathrm{FoT} = t + N * \tau_\mathrm{m}^\mathrm{new},
      30             :  * \f}
      31             :  *
      32             :  * and if it is true,
      33             :  *
      34             :  * \f{equation}
      35             :  * T_\mathrm{expr}^\mathrm{FoT} = t + \tau_\mathrm{m}^\mathrm{old}
      36             :  *      + N * \tau_\mathrm{m}^\mathrm{new},
      37             :  * \f}
      38             :  *
      39             :  * where \f$T_\mathrm{expr}^\mathrm{FoT}\f$ is the expiration time for the
      40             :  * FunctionsOfTime, \f$t\f$ is the update time,
      41             :  * \f$\tau_\mathrm{m}^\mathrm{old/new}\f$ is the measurement timescale, and
      42             :  * \f$N\f$ is the number of measurements per update.
      43             :  *
      44             :  * If \p delay_update is true, we update the functions of time
      45             :  * one (old) measurement before they actually expire.
      46             :  *
      47             :  * The choice of having the functions of time expire exactly one old measurement
      48             :  * after they are updated is arbitrary. They could expire any time between the
      49             :  * update time and one old measurement after the update. This decision was made
      50             :  * to minimize time spent waiting for the functions of time to be valid.
      51             :  *
      52             :  * Since functions of time are valid at their expiration time, we are actually
      53             :  * able to do the next measurement if the expiration time is at that
      54             :  * measurement. Thus we delay any potential waiting that may happen until the
      55             :  * subsequent measurement (and by that time, most, if not all, functions of time
      56             :  * should have been updated because an entire horizon find happened in the
      57             :  * meantime). If the expiration time was earlier than the next measurement, we'd
      58             :  * have to pause the Algorithm on the DG elements and wait until all the
      59             :  * functions of time have been updated.
      60             :  *
      61             :  * If \p delay_update is false, we update the functions of time
      62             :  * immediately, which allows for less parallelization but potentially
      63             :  * makes the control system more responsive.
      64             :  */
      65           1 : double function_of_time_expiration_time(
      66             :     double time, const DataVector& old_measurement_timescales,
      67             :     const DataVector& new_measurement_timescales, int measurements_per_update,
      68             :     bool delay_update);
      69             : 
      70             : /*!
      71             :  * \ingroup ControlSystemGroup
      72             :  * \brief Calculate the first expiration time for the FunctionsOfTime.
      73             :  *
      74             :  * The first expiration time is
      75             :  *
      76             :  * \f{equation}
      77             :  * T_\mathrm{expr}^\mathrm{FoT} = t + N * \tau_\mathrm{m},
      78             :  * \f}
      79             :  *
      80             :  * where $t$ is the initial time, $\tau_\mathrm{m}$ is the measurement
      81             :  * timescale, and \f$N\f$ is the number of measurements per update if
      82             :  * \p delay_update is true, and one less if it is false.
      83             :  */
      84           1 : double function_of_time_initial_expiration_time(
      85             :     double time, const DataVector& measurement_timescales,
      86             :     int measurements_per_update, bool delay_update);
      87             : 
      88             : /*!
      89             :  * \ingroup ControlSystemGroup
      90             :  * \brief Calculate the next expiration time for the MeasurementTimescales.
      91             :  *
      92             :  * \f{align}
      93             :  * T_\mathrm{expr}^\mathrm{m} &= T_\mathrm{expr}^\mathrm{FoT} - \frac{1}{2}
      94             :  *   \tau_\mathrm{m}^\mathrm{new} \\
      95             :  * \f}
      96             :  *
      97             :  * where \f$T_\mathrm{expr}^\mathrm{m}\f$ is the measurement expiration time,
      98             :  * \f$T_\mathrm{expr}^\mathrm{FoT}\f$ is the function of time expiration time as
      99             :  * calculate from `function_of_time_expiration_time()`, and
     100             :  * \f$\tau_\mathrm{m}^\mathrm{new}\f$ is the new measurement timescale. The
     101             :  * reason for the factor of a half is as follows:
     102             :  *
     103             :  * We update the functions of time one (old) measurement before the expiration
     104             :  * time. Based on how dense triggers are set up, which control_system::Trigger
     105             :  * is a dense trigger, you calculate the next trigger (measurement) time at the
     106             :  * current measurement time. However, at the function of time expiration time we
     107             :  * need updated damping timescales from all control systems in order to
     108             :  * calculate when the next measurement is going to be (and in turn, the next
     109             :  * measurement expiration time). Thus, at the measurement that occurs at the
     110             :  * function of time expiration time, our measurement timescales can't be valid
     111             :  * and we must wait for updated ones. To achieve this, we set the measurement
     112             :  * expiration time *before* the function of time expiration time, but *after*
     113             :  * the previous measurement (the update measurement). The factor of one half is
     114             :  * just to guarantee we are more than epsilon before the function of time
     115             :  * expiration time and more than epsilon after the update measurement.
     116             :  */
     117           1 : double measurement_expiration_time(double time,
     118             :                                    const DataVector& old_measurement_timescales,
     119             :                                    const DataVector& new_measurement_timescales,
     120             :                                    int measurements_per_update,
     121             :                                    bool delay_update);
     122             : 
     123             : /*!
     124             :  * \ingroup ControlSystemGroup
     125             :  * \brief Calculate the first expiration time for the MeasurementTimescales.
     126             :  *
     127             :  * This is related to `function_of_time_initial_expiration_time()` in
     128             :  * the same way that `measurement_expiration_time()` is related to
     129             :  * `function_of_time_expiration_time()`.
     130             :  */
     131           1 : double measurement_initial_expiration_time(
     132             :     double time, const DataVector& measurement_timescales,
     133             :     int measurements_per_update, bool delay_update);
     134             : 
     135             : /*!
     136             :  * \ingroup ControlSystemGroup
     137             :  * \brief Construct the initial expiration times for functions of time that are
     138             :  * controlled by a control system
     139             :  *
     140             :  * The expiration times are constructed using inputs from control system
     141             :  * OptionHolders as an unordered map from the name of the function of time being
     142             :  * controlled to the expiration time.
     143             :  *
     144             :  * The expiration time for each individual function of time is computed as
     145             :  * $\tau_\mathrm{exp} = \alpha_\mathrm{update} \tau_\mathrm{damp}$ where
     146             :  * $\alpha_\mathrm{update}$ is the update fraction supplied as input to the
     147             :  * Controller and $\tau_\mathrm{damp}$ is/are the damping timescales
     148             :  * supplied from the TimescaleTuner ($\tau_\mathrm{damp}$ is a DataVector
     149             :  * with as many components as the corresponding function of time, thus
     150             :  * $\tau_\mathrm{exp}$ will also be a DataVector of the same length).
     151             :  *
     152             :  * However, this expiration time calculated above is not necessarily the
     153             :  * expiration that is returned by this function. We group functions of time by
     154             :  * the `control_system::protocols::Measurement` that their corresponding control
     155             :  * systems use. This is because one measurement may be used to update many
     156             :  * functions of time. So the actual expiration time that is used for all the
     157             :  * functions of time in this group is the *minimum* of the $\tau_\mathrm{exp}$
     158             :  * of each function of time in the group.
     159             :  *
     160             :  * These groups are calculated using `control_system::system_to_combined_names`
     161             :  * where the list of control systems comes from the passed in `OptionHolder`s.
     162             :  *
     163             :  * If the control system isn't active then expiration time is
     164             :  * `std::numeric_limits<double>::infinity()`, regardless of what the groups'
     165             :  * expiration time is.
     166             :  */
     167             : template <size_t Dim, typename... OptionHolders>
     168           1 : std::unordered_map<std::string, double> initial_expiration_times(
     169             :     const double initial_time, const int measurements_per_update,
     170             :     const bool delay_update,
     171             :     const std::unique_ptr<::DomainCreator<Dim>>& domain_creator,
     172             :     const std::optional<OptionHolders>&... option_holders) {
     173             :   std::unordered_map<std::string, double> initial_expiration_times{};
     174             : 
     175             :   using control_systems = tmpl::list<typename OptionHolders::control_system...>;
     176             : 
     177             :   // First string is name of each control system. Second string is combination
     178             :   // of control systems with same measurement
     179             :   std::unordered_map<std::string, std::string> map_of_names =
     180             :       system_to_combined_names<control_systems>();
     181             : 
     182             :   std::unordered_map<std::string, double> combined_expiration_times{};
     183             :   std::unordered_set<std::string> infinite_expiration_times{};
     184             :   for (const auto& [control_system_name, combined_name] : map_of_names) {
     185             :     (void)control_system_name;
     186             :     if (combined_expiration_times.count(combined_name) != 1) {
     187             :       combined_expiration_times[combined_name] =
     188             :           std::numeric_limits<double>::infinity();
     189             :     }
     190             :   }
     191             : 
     192             :   [[maybe_unused]] const auto combine_expiration_times =
     193             :       [&initial_time, &measurements_per_update, &delay_update, &domain_creator,
     194             :        &map_of_names, &combined_expiration_times,
     195             :        &infinite_expiration_times](const auto& option_holder) {
     196             :         const std::string& control_system_name = std::decay_t<
     197             :             decltype(option_holder)>::value_type::control_system::name();
     198             :         const std::string& combined_name = map_of_names[control_system_name];
     199             : 
     200             :         // If not active, leave the expiration time as infinity
     201             :         if (not option_holder.has_value()) {
     202             :           infinite_expiration_times.insert(control_system_name);
     203             :           return;
     204             :         }
     205             : 
     206             :         auto tuner = option_holder->tuner;
     207             :         Tags::detail::initialize_tuner(make_not_null(&tuner), domain_creator,
     208             :                                        initial_time, control_system_name);
     209             : 
     210             :         const auto& controller = option_holder->controller;
     211             :         const DataVector measurement_timescales =
     212             :             calculate_measurement_timescales(controller, tuner,
     213             :                                              measurements_per_update);
     214             :         const double min_measurement_timescale = min(measurement_timescales);
     215             : 
     216             :         const double initial_expiration_time =
     217             :             function_of_time_initial_expiration_time(
     218             :                 initial_time, DataVector{1, min_measurement_timescale},
     219             :                 measurements_per_update, delay_update);
     220             : 
     221             :         combined_expiration_times[combined_name] = std::min(
     222             :             combined_expiration_times[combined_name], initial_expiration_time);
     223             :       };
     224             : 
     225             :   EXPAND_PACK_LEFT_TO_RIGHT(combine_expiration_times(option_holders));
     226             : 
     227             :   // Set all functions of time for a given measurement to the same expiration
     228             :   // time
     229             :   for (const auto& [control_system_name, combined_name] : map_of_names) {
     230             :     initial_expiration_times[control_system_name] =
     231             :         infinite_expiration_times.contains(control_system_name)
     232             :             ? std::numeric_limits<double>::infinity()
     233             :             : combined_expiration_times[combined_name];
     234             :   }
     235             : 
     236             :   return initial_expiration_times;
     237             : }
     238             : }  // namespace control_system