Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 : #include <optional>
8 : #include <pup.h>
9 : #include <string>
10 : #include <utility>
11 :
12 : #include "ControlSystem/Protocols/ControlError.hpp"
13 : #include "ControlSystem/Tags/QueueTags.hpp"
14 : #include "ControlSystem/TimescaleTuner.hpp"
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataVector.hpp"
17 : #include "Domain/Creators/Tags/ObjectCenter.hpp"
18 : #include "Domain/Structure/ObjectLabel.hpp"
19 : #include "IO/Logging/Verbosity.hpp"
20 : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
21 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
22 : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
23 : #include "Options/String.hpp"
24 : #include "Parallel/GlobalCache.hpp"
25 : #include "Utilities/Algorithm.hpp"
26 : #include "Utilities/ProtocolHelpers.hpp"
27 : #include "Utilities/TMPL.hpp"
28 : #include "Utilities/TaggedTuple.hpp"
29 :
30 : /// \cond
31 : namespace domain::Tags {
32 : struct FunctionsOfTime;
33 : } // namespace domain::Tags
34 : namespace Frame {
35 : struct Grid;
36 : struct Distorted;
37 : } // namespace Frame
38 : /// \endcond
39 :
40 : namespace control_system::ControlErrors {
41 : /*!
42 : * \brief Control error in the for the
43 : * `domain::CoordinateMaps::TimeDependent::Skew` map.
44 : *
45 : * \details Computes the error in the map parameters $F_y(t)$ and $F_z(t)$ of
46 : * the \link domain::CoordinateMaps::TimeDependent::Skew Skew \endlink using
47 : * a modified version of Eq. (71) from \cite Hemberger2012jz. This modified
48 : * control error is
49 : *
50 : * \begin{equation}
51 : * Q_{F_j} = g\left(\frac{w_A \Theta_A^j + w_B \Theta_B^j}{w_A + w_B}\right)
52 : * - (1-g)F_j
53 : * \end{equation}
54 : *
55 : * where the falloff function is assumed to be $W(\vec{x}) \approx 1$ from the
56 : * \link domain::CoordinateMaps::TimeDependent::Skew Skew \endlink map,
57 : * $\Theta_H^j$ are the inclination angles between the $x$-axis and the normal
58 : * to the horizon at the intersection point $x_H^\textrm{Int}$ between the
59 : * $x$-axis and the horizon, $w_A$ and $w_B$ are averaging weights defined by
60 : *
61 : * \begin{equation}
62 : * w_H = \exp{\left(-\frac{x^0_C - x_H^\textrm{Int}}{x^0_C - C^0_H}\right)}
63 : * \end{equation}
64 : *
65 : * with $C^0_H$ being the centers of the excision boundaries, and finally with
66 : *
67 : * \begin{equation}
68 : * g = \frac{1}{2}\left(1-\tanh\left(10\frac{x_A^\textrm{Int} -
69 : * x_B^\textrm{Int}}{C^0_A - C^0_B} - 5\right)\right).
70 : * \end{equation}
71 : *
72 : * This transition function $g$ is meant to only activate skew control when the
73 : * black holes are close to merger to avoid adverse effects with junk radiation.
74 : * If $g < 0.0025$, then we turn skew control off completely and the control
75 : * error just becomes
76 : *
77 : * \begin{equation}
78 : * Q_{F_j} = -F_j.
79 : * \end{equation}
80 : *
81 : * This threshold value of $g = 0.0025$ was chosen in SpEC and seems to work
82 : * well.
83 : *
84 : * Requirements:
85 : * - This control error requires that there be exactly two objects in the
86 : * simulation
87 : * - Currently both these objects must be black holes
88 : * - Currently this control system can only be used with the \link
89 : * control_system::Systems::Skew Skew \endlink control system
90 : */
91 1 : struct Skew : tt::ConformsTo<protocols::ControlError> {
92 0 : using object_centers = domain::object_list<>;
93 :
94 0 : using options = tmpl::list<>;
95 0 : static constexpr Options::String help{
96 : "Computes the control error for skew control."};
97 :
98 : // Explicitly defined copy constructor because DataBox isn't copyable
99 0 : Skew() = default;
100 0 : Skew(Skew&& rhs) = default;
101 0 : Skew& operator=(Skew&& rhs) = default;
102 0 : Skew(const Skew&);
103 0 : Skew& operator=(const Skew&);
104 0 : ~Skew() = default;
105 :
106 : /*!
107 : * \brief Returns the internal suggested timescale. A std::nullopt means that
108 : * no timescale is suggested.
109 : */
110 1 : std::optional<double> get_suggested_timescale() const;
111 :
112 : /*!
113 : * \brief Resets the internal suggested timescale to nullopt.
114 : */
115 1 : void reset();
116 :
117 0 : void pup(PUP::er& p);
118 :
119 : template <typename Metavariables, typename... TupleTags>
120 0 : DataVector operator()(const ::TimescaleTuner<true>& /*tuner*/,
121 : const Parallel::GlobalCache<Metavariables>& cache,
122 : const double time,
123 : const std::string& function_of_time_name,
124 : const tuples::TaggedTuple<TupleTags...>& measurements) {
125 : const ylm::Strahlkorper<Frame::Distorted>& horizon_a = tuples::get<
126 : QueueTags::Horizon<Frame::Distorted, ::domain::ObjectLabel::A>>(
127 : measurements);
128 : const ylm::Strahlkorper<Frame::Distorted>& horizon_b = tuples::get<
129 : QueueTags::Horizon<Frame::Distorted, ::domain::ObjectLabel::B>>(
130 : measurements);
131 :
132 : // Copy the horizon into the box so the compute tags use the correct
133 : // strahlkorper
134 : const auto set_horizon =
135 : [](const gsl::not_null<ylm::Strahlkorper<Frame::Distorted>*>
136 : horizon_ptr,
137 : const ylm::Strahlkorper<Frame::Distorted>& horizon) {
138 : *horizon_ptr = horizon;
139 : };
140 :
141 : db::mutate<ylm::Tags::Strahlkorper<Frame::Distorted>>(
142 : set_horizon, make_not_null(&box_a_), horizon_a);
143 : db::mutate<ylm::Tags::Strahlkorper<Frame::Distorted>>(
144 : set_horizon, make_not_null(&box_b_), horizon_b);
145 :
146 : const auto& normal_one_form_a =
147 : db::get<ylm::Tags::NormalOneForm<Frame::Distorted>>(box_a_);
148 : const auto& cartesian_coords_a =
149 : db::get<ylm::Tags::CartesianCoords<Frame::Distorted>>(box_a_);
150 : const auto& center_a =
151 : Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::A>>(
152 : cache);
153 :
154 : const auto& normal_one_form_b =
155 : db::get<ylm::Tags::NormalOneForm<Frame::Distorted>>(box_b_);
156 : const auto& cartesian_coords_b =
157 : db::get<ylm::Tags::CartesianCoords<Frame::Distorted>>(box_b_);
158 : const auto& center_b =
159 : Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::B>>(
160 : cache);
161 :
162 : // The translation of the cutting plane in the grid frame is just the
163 : // average of the two centers.
164 : const double cut_x = 0.5 * (center_a[0] + center_b[0]);
165 :
166 : // For AhA since it always has x_center > cut_x, we want the point which is
167 : // closest to the cutting plane, which in terms of theta, phi = pi/2, pi.
168 : // For AhB since it always has x_center < cut_x, we want the point which is
169 : // closest to the cutting plane, which in terms of theta, phi = pi/2, 0.
170 : // To get the coord and normal one form at these theta/phi, we do an
171 : // interpolation for both Ah's
172 : const auto& ylm_a = horizon_a.ylm_spherepack();
173 : const ylm::Spherepack::InterpolationInfo<double> interpolation_info_a =
174 : ylm_a.set_up_interpolation_info(std::array{M_PI_2, M_PI});
175 : const auto& ylm_b = horizon_b.ylm_spherepack();
176 : const ylm::Spherepack::InterpolationInfo<double> interpolation_info_b =
177 : ylm_b.set_up_interpolation_info(std::array{M_PI_2, 0.0});
178 :
179 : const auto set_intersection =
180 : [](const gsl::not_null<DataVector*> inclination_angle,
181 : const gsl::not_null<std::array<double, 3>*> intersection_coord,
182 : const ylm::Spherepack& ylm,
183 : const ylm::Spherepack::InterpolationInfo<double>& interpolation_info,
184 : const auto& normal_one_form, const auto& cartesian_coords) {
185 : std::array<double, 3> intersection_normal_one_form{};
186 : for (size_t i = 0; i < 3; i++) {
187 : ylm.interpolate(make_not_null(&gsl::at(*intersection_coord, i)),
188 : make_not_null(cartesian_coords.get(i).data()),
189 : interpolation_info);
190 : ylm.interpolate(
191 : make_not_null(&gsl::at(intersection_normal_one_form, i)),
192 : make_not_null(normal_one_form.get(i).data()),
193 : interpolation_info);
194 : }
195 :
196 : for (size_t i = 0; i < 2; ++i) {
197 : (*inclination_angle)[i] =
198 : atan2(gsl::at(intersection_normal_one_form, i + 1),
199 : (intersection_normal_one_form)[0]);
200 : if ((*inclination_angle)[i] < -0.5 * M_PI) {
201 : (*inclination_angle)[i] += M_PI;
202 : } else if ((*inclination_angle)[i] > +0.5 * M_PI) {
203 : (*inclination_angle)[i] -= M_PI;
204 : }
205 : }
206 : };
207 :
208 : // These are the values we need at the theta/phi described above
209 : std::array<double, 3> intersection_coord_a{};
210 : std::array<double, 3> intersection_coord_b{};
211 :
212 : set_intersection(make_not_null(&inclination_angle_a_),
213 : make_not_null(&intersection_coord_a), ylm_a,
214 : interpolation_info_a, normal_one_form_a,
215 : cartesian_coords_a);
216 : set_intersection(make_not_null(&inclination_angle_b_),
217 : make_not_null(&intersection_coord_b), ylm_b,
218 : interpolation_info_b, normal_one_form_b,
219 : cartesian_coords_b);
220 :
221 : const double relative_delta_x =
222 : (intersection_coord_a[0] - intersection_coord_b[0]) /
223 : (center_a[0] - center_b[0]);
224 : // Hardcoded function used in SpEC
225 : const double temporal_transition_function =
226 : 0.5 * (1.0 - tanh(10.0 * relative_delta_x - 5.0));
227 :
228 : // Hardcoded value used in SpEC
229 : constexpr double activation_threshold = 0.0025;
230 :
231 : const auto& function_of_time =
232 : Parallel::get<domain::Tags::FunctionsOfTime>(cache).at(
233 : function_of_time_name);
234 :
235 : // Only activate if BHs are close enough. Otherwise control to zero
236 : DataVector func = std::move(function_of_time->func(time)[0]);
237 : DataVector& control_error = func;
238 : if (temporal_transition_function > activation_threshold) {
239 : suggested_timescale_ =
240 : std::max(std::min(abs(center_a[0]), abs(center_b[0])),
241 : std::min(abs(cut_x - intersection_coord_a[0]),
242 : abs(cut_x - intersection_coord_b[0])));
243 :
244 : const double weight_a =
245 : exp(-(cut_x - intersection_coord_a[0]) / (cut_x - center_a[0]));
246 : const double weight_b =
247 : exp(-(cut_x - intersection_coord_b[0]) / (cut_x - center_b[0]));
248 :
249 : control_error = temporal_transition_function *
250 : (weight_a * inclination_angle_a_ +
251 : weight_b * inclination_angle_b_) /
252 : (weight_a + weight_b) -
253 : func;
254 : } else {
255 : control_error *= -1.0;
256 : }
257 :
258 : return std::move(control_error);
259 : }
260 :
261 : private:
262 0 : std::optional<double> suggested_timescale_;
263 : // not pupped. These are allocated here to avoid memory allocations during
264 : // the operator() call. All their data is temporary
265 0 : DataVector inclination_angle_a_{2, 0.0};
266 0 : DataVector inclination_angle_b_{2, 0.0};
267 : db::compute_databox_type<
268 : tmpl::append<ylm::Tags::items_tags<Frame::Distorted>,
269 : ylm::Tags::compute_items_tags<Frame::Distorted>>>
270 0 : box_a_;
271 : db::compute_databox_type<
272 : tmpl::append<ylm::Tags::items_tags<Frame::Distorted>,
273 : ylm::Tags::compute_items_tags<Frame::Distorted>>>
274 0 : box_b_;
275 : };
276 : } // namespace control_system::ControlErrors
|
