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