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 : #include <cstddef> 8 : #include <pup.h> 9 : #include <string> 10 : 11 : #include "ControlSystem/Protocols/ControlError.hpp" 12 : #include "ControlSystem/Tags/QueueTags.hpp" 13 : #include "ControlSystem/Tags/SystemTags.hpp" 14 : #include "DataStructures/DataVector.hpp" 15 : #include "Domain/Structure/ObjectLabel.hpp" 16 : #include "NumericalAlgorithms/SphericalHarmonics/SpherepackIterator.hpp" 17 : #include "Options/String.hpp" 18 : #include "Parallel/GlobalCache.hpp" 19 : #include "Utilities/EqualWithinRoundoff.hpp" 20 : #include "Utilities/ErrorHandling/Assert.hpp" 21 : #include "Utilities/ForceInline.hpp" 22 : #include "Utilities/ProtocolHelpers.hpp" 23 : #include "Utilities/TMPL.hpp" 24 : #include "Utilities/TaggedTuple.hpp" 25 : 26 : /// \cond 27 : namespace domain::Tags { 28 : template <size_t Dim> 29 : struct Domain; 30 : struct FunctionsOfTime; 31 : } // namespace domain::Tags 32 : namespace Frame { 33 : struct Distorted; 34 : } // namespace Frame 35 : /// \endcond 36 : 37 : namespace control_system { 38 : namespace ControlErrors { 39 : namespace detail { 40 : template <::domain::ObjectLabel Horizon> 41 : std::string excision_sphere_name() { 42 : return "ExcisionSphere"s + ::domain::name(Horizon); 43 : } 44 : 45 : template <::domain::ObjectLabel Horizon> 46 : std::string size_name() { 47 : return "Size"s + ::domain::name(Horizon); 48 : } 49 : } // namespace detail 50 : 51 : /*! 52 : * \brief Control error in the 53 : * \link domain::CoordinateMaps::TimeDependent::Shape Shape \endlink coordinate 54 : * map 55 : * 56 : * \details Computes the error for each \f$ l,m \f$ mode of the shape map except 57 : * for \f$ l=0 \f$ and \f$ l=1 \f$. This is because the \f$ l=0 \f$ mode is 58 : * controlled by characteristic speed control, and the \f$ l=1 \f$ mode is 59 : * redundant with the \link 60 : * domain::CoordinateMaps::TimeDependent::Translation Translation \endlink map. 61 : * The equation for the control error is Eq. (77) from \cite Hemberger2012jz 62 : * which is 63 : * 64 : * \f{align} 65 : * Q_{lm} &= -\frac{r_\mathrm{EB} - 66 : * Y_{00}\lambda_{00}(t)}{\sqrt{\frac{\pi}{2}} Y_{00}S_{00}} S_{lm} - 67 : * \lambda_{lm}(t), \quad l>=2 \label{eq:shape_control_error} 68 : * \f} 69 : * 70 : * where \f$ r_\mathrm{EB} \f$ is the radius of the excision boundary in the 71 : * grid frame, \f$ \lambda_{00}(t) \f$ is the size map parameter, \f$ 72 : * \lambda_{lm}(t) \f$ for \f$ l>=2 \f$ are the shape map parameters, and \f$ 73 : * S_{lm}\f$ are the coefficients of the harmonic expansion of the apparent 74 : * horizon. The coefficients \f$ \lambda_{lm}(t) \f$ (*not* including \f$ l=0 75 : * \f$) and \f$ S_{lm}\f$ (including \f$ l=0 \f$) are stored as the real-valued 76 : * coefficients \f$ a_{lm} \f$ and \f$ b_{lm} \f$ of the SPHEREPACK 77 : * spherical-harmonic expansion (in ylm::Spherepack). The $\lambda_{00}(t)$ 78 : * coefficient, on the other hand, is stored as the complex coefficient \f$ 79 : * A_{00} \f$ of the standard \f$ Y_{lm} \f$ decomposition. Because \f$ a_{00} = 80 : * \sqrt{\frac{2}{\pi}}A_{00} \f$, there is an extra factor of \f$ 81 : * \sqrt{\frac{\pi}{2}} \f$ in the above formula in the denominator of the 82 : * fraction multiplying the \f$ S_{00}\f$ component so it is represented in the 83 : * \f$ Y_{lm} \f$ decomposition just like \f$ r_{EB} \f$ and \f$ \lambda_{00} 84 : * \f$ are). That way, we ensure the numerator and denominator are represented 85 : * in the same way before we take their ratio. 86 : * 87 : * Requirements: 88 : * - This control error requires that there be at least one excision surface in 89 : * the simulation 90 : * - Currently this control error can only be used with the \link 91 : * control_system::Systems::Shape Shape \endlink control system 92 : */ 93 : template <::domain::ObjectLabel Horizon> 94 1 : struct Shape : tt::ConformsTo<protocols::ControlError> { 95 : // Shape needs an excision sphere 96 0 : using object_centers = domain::object_list<Horizon>; 97 : 98 0 : using options = tmpl::list<>; 99 0 : static constexpr Options::String help{ 100 : "Computes the control error for shape control. This should not take any " 101 : "options."}; 102 : 103 0 : void pup(PUP::er& /*p*/) {} 104 : 105 : template <typename Metavariables, typename... TupleTags> 106 0 : DataVector operator()(const ::TimescaleTuner<true>& /*unused*/, 107 : const Parallel::GlobalCache<Metavariables>& cache, 108 : const double time, 109 : const std::string& function_of_time_name, 110 : const tuples::TaggedTuple<TupleTags...>& measurements) { 111 : const auto& domain = get<domain::Tags::Domain<3>>(cache); 112 : const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache); 113 : const DataVector lambda_lm_coefs = 114 : functions_of_time.at(function_of_time_name)->func(time)[0]; 115 : const double lambda_00_coef = 116 : functions_of_time.at(detail::size_name<Horizon>())->func(time)[0][0]; 117 : 118 : const auto& ah = 119 : get<control_system::QueueTags::Horizon<Frame::Distorted>>(measurements); 120 : const auto& ah_coefs = ah.coefficients(); 121 : 122 : ASSERT(lambda_lm_coefs.size() == ah_coefs.size(), 123 : "Number of coefficients for shape map '" 124 : << function_of_time_name << "' (" << lambda_lm_coefs.size() 125 : << ") does not match the number of coefficients in the AH " 126 : "Strahlkorper (" 127 : << ah_coefs.size() << ")."); 128 : 129 : const auto& excision_spheres = domain.excision_spheres(); 130 : 131 : ASSERT(domain.excision_spheres().count( 132 : detail::excision_sphere_name<Horizon>()) == 1, 133 : "Excision sphere " << detail::excision_sphere_name<Horizon>() 134 : << " not in the domain but is needed to " 135 : "compute Shape control error."); 136 : 137 : const double radius_excision_sphere_grid_frame = 138 : excision_spheres.at(detail::excision_sphere_name<Horizon>()).radius(); 139 : 140 : const double Y00 = sqrt(0.25 / M_PI); 141 : ylm::SpherepackIterator iter{ah.l_max(), ah.m_max()}; 142 : // See above docs for why we have the sqrt(pi/2) in the denominator 143 : const double relative_size_factor = 144 : (radius_excision_sphere_grid_frame / Y00 - lambda_00_coef) / 145 : (sqrt(0.5 * M_PI) * ah_coefs[iter.set(0, 0)()]); 146 : 147 : // The map parameters are in terms of SPHEREPACK coefficients (just like 148 : // strahlkorper coefficients), *not* spherical harmonic coefficients, thus 149 : // the control error for each l,m is in terms of SPHEREPACK coefficients 150 : // and no extra factors of sqrt(2/pi) are needed 151 : DataVector Q = -relative_size_factor * ah_coefs - lambda_lm_coefs; 152 : 153 : // Shape control is only for l > 1 so enforce that Q=0 for l=0,l=1. These 154 : // components of the control error won't be 0 automatically because the AH 155 : // can freely have nonzero l=0 and l=1 coefficients so we have to set the 156 : // control error components to be 0 manually. 157 : Q[iter.set(0, 0)()] = 0.0; 158 : Q[iter.set(1, -1)()] = 0.0; 159 : Q[iter.set(1, 0)()] = 0.0; 160 : Q[iter.set(1, 1)()] = 0.0; 161 : 162 : return Q; 163 : } 164 : }; 165 : } // namespace ControlErrors 166 : } // namespace control_system