Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include <array> 7 : #include <cmath> 8 : #include <cstddef> 9 : #include <optional> 10 : #include <string> 11 : 12 : #include "ControlSystem/Component.hpp" 13 : #include "ControlSystem/ControlErrors/Shape.hpp" 14 : #include "ControlSystem/Measurements/BothHorizons.hpp" 15 : #include "ControlSystem/Measurements/SingleHorizon.hpp" 16 : #include "ControlSystem/Protocols/ControlError.hpp" 17 : #include "ControlSystem/Protocols/ControlSystem.hpp" 18 : #include "ControlSystem/Protocols/Measurement.hpp" 19 : #include "ControlSystem/Tags/QueueTags.hpp" 20 : #include "ControlSystem/Tags/SystemTags.hpp" 21 : #include "ControlSystem/UpdateControlSystem.hpp" 22 : #include "DataStructures/DataBox/DataBox.hpp" 23 : #include "DataStructures/DataBox/Tag.hpp" 24 : #include "DataStructures/LinkedMessageQueue.hpp" 25 : #include "Domain/Structure/ObjectLabel.hpp" 26 : #include "NumericalAlgorithms/SphericalHarmonics/SpherepackIterator.hpp" 27 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp" 28 : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp" 29 : #include "Parallel/GlobalCache.hpp" 30 : #include "Parallel/Printf/Printf.hpp" 31 : #include "ParallelAlgorithms/Actions/UpdateMessageQueue.hpp" 32 : #include "Utilities/GetOutput.hpp" 33 : #include "Utilities/ProtocolHelpers.hpp" 34 : #include "Utilities/TMPL.hpp" 35 : 36 : /// \cond 37 : namespace Frame { 38 : struct Distorted; 39 : } // namespace Frame 40 : /// \endcond 41 : 42 : namespace control_system::Systems { 43 : /*! 44 : * \brief Controls the \link domain::CoordinateMaps::TimeDependent::Shape Shape 45 : * \endlink map 46 : * 47 : * \details Controls the functions \f$ \lambda_{lm}(t) \f$ in the \link 48 : * domain::CoordinateMaps::TimeDependent::Shape Shape \endlink map to match the 49 : * shape of the excision sphere to the shape of the horizon. 50 : * 51 : * Requirements: 52 : * - This control system requires that there be at least one excision surface in 53 : * the simulation 54 : * - Currently this control system can only be used with the \link 55 : * control_system::measurements::BothHorizons BothHorizons \endlink 56 : * measurement 57 : * - Currently this control system can only be used with the \link 58 : * control_system::ControlErrors::Shape Shape \endlink control error 59 : */ 60 : template <::domain::ObjectLabel Horizon, size_t DerivOrder, 61 : typename Measurement> 62 1 : struct Shape : tt::ConformsTo<protocols::ControlSystem> { 63 0 : static constexpr size_t deriv_order = DerivOrder; 64 : 65 0 : static std::string name() { return "Shape"s + ::domain::name(Horizon); } 66 : 67 0 : static std::optional<std::string> component_name( 68 : const size_t i, const size_t num_components) { 69 : // num_components = 2 * (l_max + 1)**2 if l_max == m_max which it is for the 70 : // shape map. This is why we can divide by 2 and take the sqrt without 71 : // worrying about odd numbers or non-perfect squares 72 : const size_t l_max = -1 + sqrt(num_components / 2); 73 : ylm::SpherepackIterator iter(l_max, l_max); 74 : const auto compact_index = iter.compact_index(i); 75 : if (compact_index.has_value()) { 76 : iter.set(compact_index.value()); 77 : const int m = iter.coefficient_array() == 78 : ylm::SpherepackIterator::CoefficientArray::a 79 : ? static_cast<int>(iter.m()) 80 : : -static_cast<int>(iter.m()); 81 : return {"l"s + get_output(iter.l()) + "m"s + get_output(m)}; 82 : } else { 83 : return std::nullopt; 84 : } 85 : } 86 : 87 0 : using measurement = Measurement; 88 : static_assert( 89 : std::is_same_v<measurement, measurements::SingleHorizon<Horizon>> or 90 : std::is_same_v<measurement, measurements::BothHorizons>, 91 : "Must use either SingleHorizon or BothHorizon measurement for Shape " 92 : "control system."); 93 : static_assert( 94 : tt::conforms_to_v<measurement, control_system::protocols::Measurement>); 95 : 96 0 : using control_error = ControlErrors::Shape<Horizon>; 97 : static_assert(tt::conforms_to_v<control_error, 98 : control_system::protocols::ControlError>); 99 : 100 : // tag goes in control component 101 0 : struct MeasurementQueue : db::SimpleTag { 102 0 : using type = 103 : LinkedMessageQueue<double, 104 : tmpl::list<QueueTags::Horizon<Frame::Distorted>>>; 105 : }; 106 : 107 0 : using simple_tags = tmpl::list<MeasurementQueue>; 108 : 109 0 : struct process_measurement { 110 : template <typename Submeasurement> 111 0 : using argument_tags = tmpl::list<ylm::Tags::Strahlkorper<Frame::Distorted>>; 112 : 113 : template <typename Metavariables> 114 0 : static void apply(typename measurements::SingleHorizon< 115 : Horizon>::Submeasurement submeasurement, 116 : const ylm::Strahlkorper<Frame::Distorted>& strahlkorper, 117 : Parallel::GlobalCache<Metavariables>& cache, 118 : const LinkedMessageId<double>& measurement_id) { 119 : auto& control_sys_proxy = Parallel::get_parallel_component< 120 : ControlComponent<Metavariables, Shape>>(cache); 121 : 122 : Parallel::simple_action<::Actions::UpdateMessageQueue< 123 : QueueTags::Horizon<Frame::Distorted>, MeasurementQueue, 124 : UpdateControlSystem<Shape>>>(control_sys_proxy, measurement_id, 125 : strahlkorper); 126 : 127 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) { 128 : Parallel::printf("%s, time = %.16f: Received measurement '%s'.\n", 129 : name(), measurement_id.id, 130 : pretty_type::name(submeasurement)); 131 : } 132 : } 133 : 134 : template <::domain::ObjectLabel MeasureHorizon, typename Metavariables> 135 0 : static void apply( 136 : measurements::BothHorizons::FindHorizon<MeasureHorizon> submeasurement, 137 : const ylm::Strahlkorper<Frame::Distorted>& strahlkorper, 138 : Parallel::GlobalCache<Metavariables>& cache, 139 : const LinkedMessageId<double>& measurement_id) { 140 : // The measurement event will call this for both horizons, but we only 141 : // need one of the horizons. So if it is called for the wrong horizon, 142 : // just do nothing. 143 : if constexpr (MeasureHorizon == Horizon) { 144 : auto& control_sys_proxy = Parallel::get_parallel_component< 145 : ControlComponent<Metavariables, Shape>>(cache); 146 : 147 : Parallel::simple_action<::Actions::UpdateMessageQueue< 148 : QueueTags::Horizon<Frame::Distorted>, MeasurementQueue, 149 : UpdateControlSystem<Shape>>>(control_sys_proxy, measurement_id, 150 : strahlkorper); 151 : 152 : if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Verbose) { 153 : Parallel::printf("%s, time = %.16f: Received measurement '%s'.\n", 154 : name(), measurement_id.id, 155 : pretty_type::name(submeasurement)); 156 : } 157 : } else { 158 : (void)submeasurement; 159 : (void)strahlkorper; 160 : (void)cache; 161 : (void)measurement_id; 162 : } 163 : } 164 : }; 165 : }; 166 : } // namespace control_system::Systems