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