SpECTRE Documentation Coverage Report
Current view: top level - ControlSystem - UpdateControlSystem.hpp Hit Total Coverage
Commit: 697db44127ce57424c079f629b3283267ec6bd78 Lines: 1 4 25.0 %
Date: 2024-06-21 21:13:41
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cmath>
       7             : 
       8             : #include "ControlSystem/Averager.hpp"
       9             : #include "ControlSystem/CalculateMeasurementTimescales.hpp"
      10             : #include "ControlSystem/Component.hpp"
      11             : #include "ControlSystem/ControlErrors/Size/Update.hpp"
      12             : #include "ControlSystem/Controller.hpp"
      13             : #include "ControlSystem/ExpirationTimes.hpp"
      14             : #include "ControlSystem/IsSize.hpp"
      15             : #include "ControlSystem/Metafunctions.hpp"
      16             : #include "ControlSystem/Tags/IsActiveMap.hpp"
      17             : #include "ControlSystem/Tags/MeasurementTimescales.hpp"
      18             : #include "ControlSystem/Tags/SystemTags.hpp"
      19             : #include "ControlSystem/TimescaleTuner.hpp"
      20             : #include "ControlSystem/UpdateFunctionOfTime.hpp"
      21             : #include "ControlSystem/WriteData.hpp"
      22             : #include "DataStructures/DataBox/DataBox.hpp"
      23             : #include "DataStructures/DataVector.hpp"
      24             : #include "DataStructures/LinkedMessageId.hpp"
      25             : #include "Domain/FunctionsOfTime/Tags.hpp"
      26             : #include "Domain/Tags.hpp"
      27             : #include "IO/Logging/Verbosity.hpp"
      28             : #include "Parallel/GlobalCache.hpp"
      29             : #include "Parallel/Printf/Printf.hpp"
      30             : #include "Utilities/Gsl.hpp"
      31             : #include "Utilities/MakeString.hpp"
      32             : #include "Utilities/StdHelpers.hpp"
      33             : #include "Utilities/TaggedTuple.hpp"
      34             : 
      35             : /// \cond
      36             : namespace domain::Tags {
      37             : struct FunctionsOfTime;
      38             : }  // namespace domain::Tags
      39             : /// \endcond
      40             : 
      41             : namespace control_system {
      42             : /*!
      43             :  * \brief Functor for updating control systems when they are ready.
      44             :  *
      45             :  * \details The apply operator of this struct is meant to be called by the
      46             :  * UpdateMessageQueue action once an entire measurement has been made.
      47             :  *
      48             :  * Requires a few tags to be in the DataBox of the ControlComponent that this
      49             :  * is running on:
      50             :  * - \link Tags::Averager Averager \endlink
      51             :  * - \link Tags::Controller Controller \endlink
      52             :  * - \link Tags::TimescaleTuner TimescaleTuner \endlink
      53             :  * - \link Tags::ControlError ControlError \endlink
      54             :  * - \link Tags::WriteDataToDisk WriteDataToDisk \endlink
      55             :  * - \link Tags::CurrentNumberOfMeasurements CurrentNumberOfMeasurements
      56             :  *   \endlink
      57             :  *
      58             :  * And the \link control_system::Tags::MeasurementsPerUpdate \endlink must be in
      59             :  * the GlobalCache. If these tags are not present, a build error will occur.
      60             :  *
      61             :  * The algorithm to determine whether or not to update the functions of time is
      62             :  * as follows:
      63             :  * 1. Ensure this control system is active. If it isn't, end here and don't
      64             :  *    process the measurements.
      65             :  * 2. Increment the current measurement stored by
      66             :  *    `control_system::Tags::CurrentNumberOfMeasurements`. This keeps track of
      67             :  *    which measurement we are on out of
      68             :  *    `control_system::Tags::MeasurementsPerUpdate`.
      69             :  * 3. Calculate the control error and update the averager (store the current
      70             :  *    measurement). If the averager doesn't have enough information to determine
      71             :  *    the derivatives of the control error, exit now (this usually only happens
      72             :  *    for the first couple measurements of a simulation). If this is \link
      73             :  *    control_system::Systems::Size size \endlink control, an extra step happens
      74             :  *    after we calculate the control error, but before we update the averager.
      75             :  *    See `control_system::size::update_averager` for this step.
      76             :  * 4. If the `control_system::Tags::WriteDataToDisk` tag is set to `true`, write
      77             :  *    the time, function of time values, and the control error and its
      78             :  *    derivative to disk.
      79             :  * 5. Determine if we need to update. We only want to update when the current
      80             :  *    measurement is equal to the number of measurements per update. If it's not
      81             :  *    time to update, exit now. Once we determine that it is time to update, set
      82             :  *    the current measurement back to 0.
      83             :  * 6. Compute the control signal using the control error and its derivatives and
      84             :  *    update the damping timescales in the TimescaleTuner. If this is \link
      85             :  *    control_system::Systems::Size size \endlink control, there is an extra
      86             :  *    step after we update the damping timescale. See
      87             :  *    `control_system::size::update_tuner` for this step.
      88             :  * 7. Calculate the new measurement timescale based off the updated damping
      89             :  *    timescales and the number of measurements per update.
      90             :  * 8. Determine the new expiration times for the
      91             :  *    `::domain::Tags::FunctionsOfTime` and
      92             :  *    `control_system::Tags::MeasurementTimescales`. Call the
      93             :  *    `control_system::AggregateUpdate` simple action on the first control
      94             :  *    system in the `component_list` of the metavariables. This simple action
      95             :  *    will mutate the global cache tags when it has enough data.
      96             :  */
      97             : template <typename ControlSystem>
      98           1 : struct UpdateControlSystem {
      99           0 :   static constexpr size_t deriv_order = ControlSystem::deriv_order;
     100             : 
     101             :   template <typename DbTags, typename Metavariables, typename ArrayIndex,
     102             :             typename... TupleTags>
     103           0 :   static void apply(const gsl::not_null<db::DataBox<DbTags>*> box,
     104             :                     Parallel::GlobalCache<Metavariables>& cache,
     105             :                     const ArrayIndex& /*array_index*/, const double time,
     106             :                     tuples::TaggedTuple<TupleTags...> data) {
     107             :     const std::string& function_of_time_name = ControlSystem::name();
     108             : 
     109             :     // Begin step 1
     110             :     // If this control system isn't active, don't do anything
     111             :     if (not Parallel::get<control_system::Tags::IsActiveMap>(cache).at(
     112             :             function_of_time_name)) {
     113             :       return;
     114             :     }
     115             : 
     116             :     // Begin step 2
     117             :     const int measurements_per_update =
     118             :         get<control_system::Tags::MeasurementsPerUpdate>(cache);
     119             :     int& current_measurement = db::get_mutable_reference<
     120             :         control_system::Tags::CurrentNumberOfMeasurements>(box);
     121             : 
     122             :     ++current_measurement;
     123             : 
     124             :     const auto& functions_of_time =
     125             :         Parallel::get<::domain::Tags::FunctionsOfTime>(cache);
     126             :     const auto& function_of_time = functions_of_time.at(function_of_time_name);
     127             : 
     128             :     // Begin step 3
     129             :     // Get the averager, controller, tuner, and control error from the box
     130             :     auto& averager = db::get_mutable_reference<
     131             :         control_system::Tags::Averager<ControlSystem>>(box);
     132             :     auto& controller = db::get_mutable_reference<
     133             :         control_system::Tags::Controller<ControlSystem>>(box);
     134             :     auto& tuner = db::get_mutable_reference<
     135             :         control_system::Tags::TimescaleTuner<ControlSystem>>(box);
     136             :     auto& control_error = db::get_mutable_reference<
     137             :         control_system::Tags::ControlError<ControlSystem>>(box);
     138             :     const DataVector& current_timescale = tuner.current_timescale();
     139             : 
     140             :     // Compute control error
     141             :     const DataVector Q =
     142             :         control_error(tuner, cache, time, function_of_time_name, data);
     143             : 
     144             :     if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
     145             :       Parallel::printf(
     146             :           "%s, time = %.16f: current measurement = %d, control error = %s\n",
     147             :           function_of_time_name, time, current_measurement, Q);
     148             :     }
     149             : 
     150             :     if constexpr (size::is_size_v<ControlSystem>) {
     151             :       // This function must be called after we update the control error because
     152             :       // it uses the new state in its logic (to repopulate the averager).
     153             :       size::update_averager(make_not_null(&averager),
     154             :                             make_not_null(&control_error), cache, time,
     155             :                             tuner.current_timescale(), function_of_time_name,
     156             :                             current_measurement);
     157             :     }
     158             : 
     159             :     // Update the averager. We do this for every measurement because we still
     160             :     // want the averager to be up-to-date even if we aren't updating at this
     161             :     // time
     162             :     averager.update(time, Q, current_timescale);
     163             : 
     164             :     // Get the averaged values of the control error and its derivatives
     165             :     const auto& opt_avg_values = averager(time);
     166             : 
     167             :     if (not opt_avg_values.has_value()) {
     168             :       if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
     169             :         Parallel::printf(
     170             :             "%s, time = %.16f: current measurement = %d, Averager does not "
     171             :             "have enough data.\n",
     172             :             function_of_time_name, time, current_measurement);
     173             :       }
     174             :       return;
     175             :     }
     176             : 
     177             :     // Begin step 4
     178             :     // Write data for every measurement
     179             :     std::array<DataVector, 2> q_and_dtq{{Q, {Q.size(), 0.0}}};
     180             :     q_and_dtq[0] = (*opt_avg_values)[0];
     181             :     if constexpr (deriv_order > 1) {
     182             :       q_and_dtq[1] = (*opt_avg_values)[1];
     183             :     }
     184             : 
     185             :     if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
     186             :       // LCOV_EXCL_START
     187             :       write_components_to_disk<ControlSystem>(time, cache, function_of_time,
     188             :                                               q_and_dtq, current_timescale);
     189             :       // LCOV_EXCL_STOP
     190             :     }
     191             : 
     192             :     // Begin step 5
     193             :     // Check if it is time to update
     194             :     if (current_measurement != measurements_per_update) {
     195             :       if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
     196             :         Parallel::printf(
     197             :             "%s, time = %.16f: current measurement = %d is not at "
     198             :             "measurements_per_update = %d. Waiting for more data.\n",
     199             :             function_of_time_name, time, current_measurement,
     200             :             measurements_per_update);
     201             :       }
     202             :       return;
     203             :     }
     204             : 
     205             :     // Set current measurement back to 0 because we're updating now
     206             :     current_measurement = 0;
     207             : 
     208             :     // Begin step 6
     209             :     const double time_offset =
     210             :         averager.last_time_updated() - averager.average_time(time);
     211             :     const double time_offset_0th =
     212             :         averager.using_average_0th_deriv_of_q() ? time_offset : 0.0;
     213             : 
     214             :     // Calculate the control signal which will be used to update the highest
     215             :     // derivative of the FunctionOfTime
     216             :     const DataVector control_signal =
     217             :         controller(time, current_timescale, opt_avg_values.value(),
     218             :                    time_offset_0th, time_offset);
     219             : 
     220             :     tuner.update_timescale(q_and_dtq);
     221             : 
     222             :     if constexpr (size::is_size_v<ControlSystem>) {
     223             :       size::update_tuner(make_not_null(&tuner), make_not_null(&control_error),
     224             :                          cache, time, function_of_time_name);
     225             :     }
     226             : 
     227             :     // Begin step 7
     228             :     // Calculate new measurement timescales with updated damping timescales
     229             :     const DataVector new_measurement_timescale =
     230             :         calculate_measurement_timescales(controller, tuner,
     231             :                                          measurements_per_update);
     232             : 
     233             :     const auto& measurement_timescales =
     234             :         Parallel::get<Tags::MeasurementTimescales>(cache);
     235             :     const auto& system_to_combined_names =
     236             :         Parallel::get<Tags::SystemToCombinedNames>(cache);
     237             :     const auto& measurement_timescale = measurement_timescales.at(
     238             :         system_to_combined_names.at(function_of_time_name));
     239             :     const double current_fot_expiration_time =
     240             :         function_of_time->time_bounds()[1];
     241             :     const double current_measurement_expiration_time =
     242             :         measurement_timescale->time_bounds()[1];
     243             :     // This call is ok because the measurement timescales are still valid
     244             :     // because the measurement timescales expire half a measurement after this
     245             :     // time.
     246             :     const DataVector old_measurement_timescale =
     247             :         measurement_timescale->func(time)[0];
     248             : 
     249             :     // Begin step 8
     250             :     // Calculate the next expiration times for both the functions of time and
     251             :     // the measurement timescales based on the current time. Then, actually
     252             :     // update the functions of time and measurement timescales
     253             :     const double new_fot_expiration_time = function_of_time_expiration_time(
     254             :         time, old_measurement_timescale, new_measurement_timescale,
     255             :         measurements_per_update);
     256             : 
     257             :     const double new_measurement_expiration_time = measurement_expiration_time(
     258             :         time, old_measurement_timescale, new_measurement_timescale,
     259             :         measurements_per_update);
     260             : 
     261             :     if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) {
     262             :       Parallel::printf("%s, time = %.16f: Control signal = %s\n",
     263             :                        function_of_time_name, time, control_signal);
     264             :       Parallel::printf(
     265             :           "%s, time = %.16f: Updating the functions of time.\n"
     266             :           " min(old_measure_timescale) = %.16f\n"
     267             :           " min(new_measure_timescale) = %.16f\n"
     268             :           " old_measure_expr_time = %.16f\n"
     269             :           " new_measure_expr_time = %.16f\n"
     270             :           " old_function_expr_time = %.16f\n"
     271             :           " new_function_expr_time = %.16f\n",
     272             :           function_of_time_name, time, min(old_measurement_timescale),
     273             :           min(new_measurement_timescale), current_measurement_expiration_time,
     274             :           new_measurement_expiration_time, current_fot_expiration_time,
     275             :           new_fot_expiration_time);
     276             :     }
     277             : 
     278             :     using first_control_component =
     279             :         tmpl::front<metafunctions::all_control_components<Metavariables>>;
     280             : 
     281             :     auto& first_control_system_proxy =
     282             :         Parallel::get_parallel_component<first_control_component>(cache);
     283             : 
     284             :     Parallel::simple_action<AggregateUpdate<ControlSystem>>(
     285             :         first_control_system_proxy, new_measurement_timescale,
     286             :         current_measurement_expiration_time, new_measurement_expiration_time,
     287             :         control_signal, current_fot_expiration_time, new_fot_expiration_time);
     288             :   }
     289             : };
     290             : }  // namespace control_system

Generated by: LCOV version 1.14