Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <memory>
8 : #include <optional>
9 : #include <pup.h>
10 : #include <string>
11 :
12 : #include "ControlSystem/Averager.hpp"
13 : #include "ControlSystem/ControlErrors/Size/ComovingCharSpeedDerivative.hpp"
14 : #include "ControlSystem/ControlErrors/Size/Error.hpp"
15 : #include "ControlSystem/ControlErrors/Size/Info.hpp"
16 : #include "ControlSystem/ControlErrors/Size/State.hpp"
17 : #include "ControlSystem/ControlErrors/Size/StateHistory.hpp"
18 : #include "ControlSystem/Protocols/ControlError.hpp"
19 : #include "ControlSystem/Tags/QueueTags.hpp"
20 : #include "ControlSystem/Tags/SystemTags.hpp"
21 : #include "ControlSystem/TimescaleTuner.hpp"
22 : #include "DataStructures/DataBox/DataBox.hpp"
23 : #include "DataStructures/DataBox/Prefixes.hpp"
24 : #include "DataStructures/DataVector.hpp"
25 : #include "DataStructures/TaggedTuple.hpp"
26 : #include "DataStructures/Tensor/Tensor.hpp"
27 : #include "Domain/Creators/Tags/Domain.hpp"
28 : #include "Domain/Structure/ObjectLabel.hpp"
29 : #include "IO/Logging/Verbosity.hpp"
30 : #include "IO/Observer/ReductionActions.hpp"
31 : #include "NumericalAlgorithms/Interpolation/ZeroCrossingPredictor.hpp"
32 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
33 : #include "NumericalAlgorithms/Strahlkorper/Tags.hpp"
34 : #include "Options/Auto.hpp"
35 : #include "Options/String.hpp"
36 : #include "Parallel/GlobalCache.hpp"
37 : #include "Parallel/Printf/Printf.hpp"
38 : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
39 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
40 : #include "Utilities/ConstantExpressions.hpp"
41 : #include "Utilities/Gsl.hpp"
42 : #include "Utilities/ProtocolHelpers.hpp"
43 : #include "Utilities/TMPL.hpp"
44 :
45 : /// \cond
46 : namespace domain::Tags {
47 : struct FunctionsOfTime;
48 : } // namespace domain::Tags
49 : namespace Frame {
50 : struct Grid;
51 : struct Distorted;
52 : } // namespace Frame
53 : /// \endcond
54 :
55 : namespace control_system {
56 1 : namespace size {
57 : /*!
58 : * \brief Function that computes the control error for
59 : * `control_system::size::States::DeltaR`.
60 : *
61 : * This is helpful to have calculated separately because other control errors
62 : * may make use of this quantity. The equation for the control error is given in
63 : * Eq. 96 in \cite Hemberger2012jz.
64 : *
65 : * \param horizon_00 The $l=0,m=0$ coefficient of the apparent horizon in the
66 : * distorted frame.
67 : * \param dt_horizon_00 The $l=0,m=0$ coefficient of the time derivative of the
68 : * apparent horizon in the distorted frame, where the derivative is taken in the
69 : * distorted frame as well.
70 : * \param lambda_00 The $l=0,m=0$ component of the function of time for the time
71 : * dependent map
72 : * \param dt_lambda_00 The $l=0,m=0$ component of the time derivative of the
73 : * function of time for the time dependent map
74 : * \param grid_frame_excision_sphere_radius Radius of the excision sphere in the
75 : * grid frame
76 : */
77 1 : double control_error_delta_r(const double horizon_00,
78 : const double dt_horizon_00, const double lambda_00,
79 : const double dt_lambda_00,
80 : const double grid_frame_excision_sphere_radius);
81 : } // namespace size
82 :
83 : namespace ControlErrors {
84 : /*!
85 : * \brief Control error in the for the \f$l=0\f$ component of the
86 : * `domain::CoordinateMaps::TimeDependent::Shape` map.
87 : *
88 : * \details The goal of this control error is
89 : *
90 : * 1. Keep the excision sphere inside the horizon
91 : * 2. Maintain a fixed distance between the excision surface and the horizon
92 : * surface.
93 : * 3. Prevent the characteristic field \f$ u^-_{ab} \f$ associated with the
94 : * characteristic speed \f$ v_- \f$ in `gh::characteristic_speeds` from coming
95 : * into the domain.
96 : *
97 : * For a more detailed account of how this is accomplished, see
98 : * `control_system::size::State` and `control_system::size::control_error` which
99 : * this class calls.
100 : *
101 : * This class holds a `control_system::size::Info` and three different
102 : * `intrp::ZeroCrossingPredictor`s internally which are needed to calculate the
103 : * `control_system::size::control_error`. Additionally, this class stores a
104 : * history of control errors for all `control_system::size::State`s using a
105 : * `control_system::size::StateHistory`. This is useful for when a discontinuous
106 : * change happens (switching `control_system::size::State`s) and we need to
107 : * repopulate the `Averager` with a history of the control error. It also
108 : * conforms to the `control_system::protocols::ControlError` protocol.
109 : *
110 : * In order to calculate the control error, we need the $\ell = 0, m = 0$
111 : * coefficient of the horizon and its time derivative. However, because we will
112 : * be finding the horizon fairly often, the value of the coefficient and its
113 : * derivative won't change smoothly because of horizon finder noise (different
114 : * number of iterations). But we expect these quantities to be smooth when
115 : * making state decisions. So to account for this, we use an `Averager` and a
116 : * `TimescaleTuner` to smooth out the horizon coefficient and get its
117 : * derivative. Every measurement, we update this smoothing averager with the
118 : * $\ell = 0, m = 0$ coefficient of the horizon and the current smoothing
119 : * timescale. Then, once we have enough measurements, we use the
120 : * `Averager::operator()` to get the averaged coefficient and its time
121 : * derivative. Since `Averager%s` calculate the average at an "averaged time",
122 : * we have to account for this small offset from the current time with a simple
123 : * Taylor expansion. Then we use this newly averaged and corrected coefficient
124 : * (and time derivative) in our calculation of the control error. The timescale
125 : * in the smoothing `TimescaleTuner` is then updated using the difference
126 : * between the averaged and un-averaged coefficient (and its time derivative).
127 : *
128 : * In addition to calculating the control error, if the
129 : * `control_system::Tags::WriteDataToDisk` tag inside the
130 : * `Parallel::GlobalCache` is true, then a diagnostic file named
131 : * `Diagnostics.dat` is also written to the same group that
132 : * `control_system::write_components_to_disk` would write the standard control
133 : * system output (`/ControlSystems/Size/`). The columns of this diagnostic file
134 : * are as follows (with a small explanation if the name isn't clear):
135 : *
136 : * - %Time
137 : * - ControlError
138 : * - StateNumber: Result of `control_system::size::State::number()`
139 : * - DiscontinuousChangeHasOccurred: 1.0 for true, 0.0 for false.
140 : * - FunctionOfTime
141 : * - DtFunctionOfTime
142 : * - HorizonCoef00
143 : * - AveragedDtHorizonCoef00: The averaged 00 component of the horizon
144 : * (averaging scheme detailed above.)
145 : * - RawDtHorizonCoef00: The raw 00 component of the horizon passed in to the
146 : * control error.
147 : * - SmootherTimescale: Damping timescale for the averaging of DtHorizonCoef00.
148 : * - MinDeltaR: The minimum of the `gr::surfaces::radial_distance` between the
149 : * horizon and the excision surfaces.
150 : * - MinRelativeDeltaR: MinDeltaR divided by the
151 : * `ylm::Strahlkorper::average_radius` of the horizon
152 : * - AvgDeltaR: Same as MinDeltaR except it's the average radii.
153 : * - AvgRelativeDeltaR: AvgDeltaR divided by the average radius of the horizon
154 : * - ControlErrorDeltaR: \f$ \dot{S}_{00} (\lambda_{00} -
155 : * r_{\mathrm{excision}}^{\mathrm{grid}} / Y_{00}) / S_{00} -
156 : * \dot{\lambda}_{00} \f$
157 : * - TargetCharSpeed
158 : * - MinCharSpeed
159 : * - MinComovingCharSpeed: Eq. 98 in \cite Hemberger2012jz
160 : * - CharSpeedCrossingTime: %Time at which the min char speed is predicted to
161 : * cross zero and become negative (or 0.0 if that time is in the past).
162 : * - ComovingCharSpeedCrossingTime: %Time at which the min comoving char speed
163 : * is predicted to cross zero and become negative (or 0.0 if that time is in
164 : * the past).
165 : * - DeltaRCrossingTime: %Time at which the distance between the excision and
166 : * horizon surfaces is predicted to be zero (or 0.0 if that time is in the
167 : * past).
168 : * - SuggestedTimescale: A timescale for the `TimescaleTuner` suggested by one
169 : * of the State%s (or 0.0 if no timescale was suggested)
170 : * - DampingTime
171 : */
172 : template <size_t DerivOrder, ::domain::ObjectLabel Horizon>
173 1 : struct Size : tt::ConformsTo<protocols::ControlError> {
174 0 : using object_centers = domain::object_list<Horizon>;
175 :
176 0 : struct MaxNumTimesForZeroCrossingPredictor {
177 : // Int so we get proper bounds checking
178 0 : using type = int;
179 0 : static constexpr Options::String help{
180 : "The maximum number of times used to calculate the zero crossing of "
181 : "the char speeds."};
182 0 : static int lower_bound() { return 3; }
183 : };
184 :
185 0 : struct SmoothAvgTimescaleFraction {
186 0 : using type = double;
187 0 : static constexpr Options::String help{
188 : "Average timescale fraction for smoothing horizon measurements."};
189 : };
190 :
191 0 : struct SmootherTuner {
192 0 : using type = TimescaleTuner<true>;
193 0 : static constexpr Options::String help{
194 : "TimescaleTuner for smoothing horizon measurements."};
195 : };
196 :
197 0 : struct InitialState {
198 0 : using type = std::unique_ptr<size::State>;
199 0 : static constexpr Options::String help{"Initial state to start in."};
200 : };
201 :
202 0 : struct DeltaRDriftOutwardOptions {
203 0 : using type =
204 : Options::Auto<DeltaRDriftOutwardOptions, Options::AutoLabel::None>;
205 0 : static constexpr Options::String help{
206 : "Options for State DeltaRDriftOutward. Specify 'None' to disable State "
207 : "DeltaRDriftOutward."};
208 0 : struct MaxAllowedRadialDistance {
209 0 : using type = double;
210 0 : static constexpr Options::String help{
211 : "Drift excision boundary outward if distance from horizon to "
212 : "excision exceeds this."};
213 : };
214 0 : struct OutwardDriftVelocity {
215 0 : using type = double;
216 0 : static constexpr Options::String help{
217 : "Constant drift velocity term, if triggered by "
218 : "MaxAllowedRadialDistance."};
219 : };
220 0 : struct OutwardDriftTimescale {
221 0 : using type = double;
222 0 : static constexpr Options::String help{
223 : "Denominator in non-constant drift velocity term, if triggered by "
224 : "MaxAllowedRadialDistance."};
225 : };
226 0 : using options = tmpl::list<MaxAllowedRadialDistance, OutwardDriftVelocity,
227 : OutwardDriftTimescale>;
228 0 : DeltaRDriftOutwardOptions();
229 0 : DeltaRDriftOutwardOptions(double max_allowed_radial_distance_in,
230 : double outward_drift_velocity_in,
231 : double outward_drift_timescale_in);
232 0 : void pup(PUP::er& p);
233 :
234 0 : double max_allowed_radial_distance{};
235 0 : double outward_drift_velocity{};
236 0 : double outward_drift_timescale{};
237 :
238 0 : friend bool operator==(const DeltaRDriftOutwardOptions& lhs,
239 : const DeltaRDriftOutwardOptions& rhs) {
240 : return lhs.max_allowed_radial_distance ==
241 : rhs.max_allowed_radial_distance and
242 : lhs.outward_drift_velocity == rhs.outward_drift_velocity and
243 : lhs.outward_drift_timescale == rhs.outward_drift_timescale;
244 : }
245 :
246 0 : friend bool operator!=(const DeltaRDriftOutwardOptions& lhs,
247 : const DeltaRDriftOutwardOptions& rhs) {
248 : return not(lhs == rhs);
249 : }
250 : };
251 :
252 0 : struct DeltaRDriftInwardOptions {
253 0 : using type =
254 : Options::Auto<DeltaRDriftInwardOptions, Options::AutoLabel::None>;
255 0 : static constexpr Options::String help{
256 : "Options for State DeltaRDriftInward. Specify 'None' to disable State "
257 : "DeltaRDriftInward."};
258 0 : struct MinAllowedRadialDistance {
259 0 : using type = double;
260 0 : static constexpr Options::String help{
261 : "Drift excision boundary inward if distance from horizon to "
262 : "excision is less than this."};
263 : };
264 0 : struct MinAllowedCharSpeed {
265 0 : using type = double;
266 0 : static constexpr Options::String help{
267 : "Drift excision boundary inward if min char speed is less than "
268 : "this."};
269 : };
270 0 : struct InwardDriftVelocity {
271 0 : using type = double;
272 0 : static constexpr Options::String help{
273 : "Maximum value of drift velocity term, if State DeltaRDriftInward is "
274 : "triggered by MinAllowedRadialDistance or MinAllowedCharSpeed."};
275 : };
276 0 : using options = tmpl::list<MinAllowedRadialDistance, MinAllowedCharSpeed,
277 : InwardDriftVelocity>;
278 0 : DeltaRDriftInwardOptions();
279 0 : DeltaRDriftInwardOptions(double min_allowed_radial_distance_in,
280 : double min_allowed_char_speed_in,
281 : double inward_drift_velocity_in);
282 0 : void pup(PUP::er& p);
283 :
284 0 : double min_allowed_radial_distance{};
285 0 : double min_allowed_char_speed{};
286 0 : double inward_drift_velocity{};
287 :
288 0 : friend bool operator==(const DeltaRDriftInwardOptions& lhs,
289 : const DeltaRDriftInwardOptions& rhs) {
290 : return lhs.min_allowed_radial_distance ==
291 : rhs.min_allowed_radial_distance and
292 : lhs.min_allowed_char_speed == rhs.min_allowed_char_speed and
293 : lhs.inward_drift_velocity == rhs.inward_drift_velocity;
294 : }
295 :
296 0 : friend bool operator!=(const DeltaRDriftInwardOptions& lhs,
297 : const DeltaRDriftInwardOptions& rhs) {
298 : return not(lhs == rhs);
299 : }
300 : };
301 :
302 0 : using options =
303 : tmpl::list<MaxNumTimesForZeroCrossingPredictor,
304 : SmoothAvgTimescaleFraction, SmootherTuner, InitialState,
305 : DeltaRDriftOutwardOptions, DeltaRDriftInwardOptions>;
306 0 : static constexpr Options::String help{
307 : "Computes the control error for size control. Will also write a "
308 : "diagnostics file if the control systems are allowed to write data to "
309 : "disk."};
310 :
311 0 : Size() = default;
312 0 : Size(const Size& rhs);
313 0 : Size& operator=(const Size& rhs);
314 0 : Size(Size&& /*rhs*/) = default;
315 0 : Size& operator=(Size&& /*rhs*/) = default;
316 0 : virtual ~Size() = default;
317 :
318 : /*!
319 : * \brief Initializes the `intrp::ZeroCrossingPredictor`s and the horizon
320 : * smoothing `Averager` and `TimescaleTuner`.
321 : *
322 : * \details All `intrp::ZeroCrossingPredictor`s are initialized with a minimum
323 : * number of times 3 and a maximum number of times `max_times`. The internal
324 : * `control_system::size::Info::state` is initialized to
325 : * `control_system::size::States::Initial`. The smoothing `Averager` uses the
326 : * input average timescale fraction and always smooths the "0th" deriv (aka
327 : * the horizon coefficients themselves). The input smoothing `TimescaleTuner`
328 : * is moved inside this class.
329 : */
330 1 : Size(const int max_times, const double smooth_avg_timescale_frac,
331 : TimescaleTuner<true> smoother_tuner,
332 : std::unique_ptr<size::State> initial_state,
333 : std::optional<DeltaRDriftOutwardOptions> delta_r_drift_outward_options,
334 : std::optional<DeltaRDriftInwardOptions> delta_r_drift_inward_options);
335 :
336 : /// Returns the internal `control_system::size::Info::suggested_time_scale`. A
337 : /// std::nullopt means that no timescale is suggested.
338 1 : std::optional<double> get_suggested_timescale() const;
339 :
340 : /*!
341 : * \brief Check if the `control_system::size::control_error` has decided to
342 : * switch states. Returns the internal
343 : * `control_system::size::Info::discontinuous_change_has_occurred`.
344 : */
345 1 : bool discontinuous_change_has_occurred() const;
346 :
347 : /*!
348 : * \brief Reset the internal `control_system::size::Info` using
349 : * `control_system::size::Info::reset`.
350 : */
351 1 : void reset();
352 :
353 : /*!
354 : * \brief Get a history of the control errors for the past few measurements.
355 : *
356 : * \return std::deque<std::pair<double, double>> This returns up to
357 : * `DerivOrder` entries, not including the most recent time. \see
358 : * `control_system::size::StateHistory::state_history`
359 : */
360 1 : std::deque<std::pair<double, double>> control_error_history() const;
361 :
362 0 : void pup(PUP::er& p);
363 :
364 0 : friend bool operator==(const Size& lhs, const Size& rhs) {
365 : return lhs.smoother_tuner_ == rhs.smoother_tuner_ and
366 : lhs.horizon_coef_averager_ == rhs.horizon_coef_averager_ and
367 : lhs.info_ == rhs.info_ and
368 : lhs.char_speed_predictor_ == rhs.char_speed_predictor_ and
369 : lhs.comoving_char_speed_predictor_ ==
370 : rhs.comoving_char_speed_predictor_ and
371 : lhs.delta_radius_predictor_ == rhs.delta_radius_predictor_ and
372 : lhs.drift_limit_char_speed_predictor_ ==
373 : rhs.drift_limit_char_speed_predictor_ and
374 : lhs.drift_limit_delta_radius_predictor_ ==
375 : rhs.drift_limit_delta_radius_predictor_ and
376 : lhs.state_history_ == rhs.state_history_ and
377 : lhs.legend_ == rhs.legend_ and
378 : lhs.subfile_name_ == rhs.subfile_name_ and
379 : lhs.delta_r_drift_outward_options_ ==
380 : rhs.delta_r_drift_outward_options_ and
381 : lhs.delta_r_drift_inward_options_ ==
382 : rhs.delta_r_drift_inward_options_;
383 : }
384 :
385 0 : friend bool operator!=(const Size& lhs, const Size& rhs) {
386 : return not(lhs == rhs);
387 : }
388 :
389 : /*!
390 : * \brief Actually computes the control error.
391 : *
392 : * \details The internal `control_system::size::Info::damping_time` is updated
393 : * to the minimum of the `TimescaleTuner::current_timescale()` that is passed
394 : * in. Also expects these queue tags to be in the `measurements` argument:
395 : *
396 : * - `ylm::Tags::Strahlkorper<Frame::Distorted>`
397 : * - `QueueTags::ExcisionSurface<Frame::Distorted>`
398 : * - `::Tags::dt<ylm::Tags::Strahlkorper<Frame::Distorted>>`
399 : * - `QueueTags::LapseOnExcisionSurface`
400 : * - `QueueTags::ShiftyQuantity<Frame::Distorted>`
401 : * - `QueueTags::SpatialMetricOnExcisionSurface<Frame::Distorted>`
402 : * - `QueueTags::InverseSpatialMetricOnExcisionSurface<Frame::Distorted>`
403 : *
404 : * \return DataVector should be of size 1
405 : */
406 : template <typename Metavariables, typename... TupleTags>
407 1 : DataVector operator()(const ::TimescaleTuner<false>& tuner,
408 : const Parallel::GlobalCache<Metavariables>& cache,
409 : const double time,
410 : const std::string& function_of_time_name,
411 : const tuples::TaggedTuple<TupleTags...>& measurements) {
412 : const Domain<3>& domain = get<domain::Tags::Domain<3>>(cache);
413 : const auto& excision_spheres = domain.excision_spheres();
414 : const auto& excision_sphere =
415 : excision_spheres.at("ExcisionSphere" + get_output(Horizon));
416 : const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
417 :
418 : const auto& excision_quantities =
419 : tuples::get<QueueTags::SizeExcisionQuantities<Frame::Distorted>>(
420 : measurements);
421 : const auto& horizon_quantities =
422 : tuples::get<QueueTags::SizeHorizonQuantities<Frame::Distorted>>(
423 : measurements);
424 :
425 : const double grid_frame_excision_sphere_radius = excision_sphere.radius();
426 : const ylm::Strahlkorper<Frame::Distorted>& apparent_horizon =
427 : tuples::get<ylm::Tags::Strahlkorper<Frame::Distorted>>(
428 : horizon_quantities);
429 : const ylm::Strahlkorper<Frame::Distorted>& excision_surface =
430 : tuples::get<QueueTags::ExcisionSurface<Frame::Distorted>>(
431 : excision_quantities);
432 : const ylm::Strahlkorper<Frame::Distorted>& time_deriv_apparent_horizon =
433 : tuples::get<::Tags::dt<ylm::Tags::Strahlkorper<Frame::Distorted>>>(
434 : horizon_quantities);
435 : const Scalar<DataVector>& lapse =
436 : tuples::get<QueueTags::LapseOnExcisionSurface>(excision_quantities);
437 : const tnsr::I<DataVector, 3, Frame::Distorted>& shifty_quantity =
438 : tuples::get<QueueTags::ShiftyQuantity<Frame::Distorted>>(
439 : excision_quantities);
440 : const tnsr::ii<DataVector, 3, Frame::Distorted>&
441 : spatial_metric_on_excision = tuples::get<
442 : QueueTags::SpatialMetricOnExcisionSurface<Frame::Distorted>>(
443 : excision_quantities);
444 : const tnsr::II<DataVector, 3, Frame::Distorted>&
445 : inverse_spatial_metric_on_excision = tuples::get<
446 : QueueTags::InverseSpatialMetricOnExcisionSurface<Frame::Distorted>>(
447 : excision_quantities);
448 : const tnsr::Ijj<DataVector, 3, Frame::Distorted>& spatial_christoffel =
449 : tuples::get<QueueTags::SpatialChristoffelSecondKind<Frame::Distorted>>(
450 : excision_quantities);
451 : const tnsr::i<DataVector, 3, Frame::Distorted>& deriv_lapse =
452 : tuples::get<QueueTags::DerivLapse<Frame::Distorted>>(
453 : excision_quantities);
454 : const tnsr::iJ<DataVector, 3, Frame::Distorted>& deriv_shift =
455 : tuples::get<QueueTags::DerivShift<Frame::Distorted>>(
456 : excision_quantities);
457 : const ::InverseJacobian<DataVector, 3, Frame::Grid,
458 : Frame::Distorted>& inv_jac_grid_to_distorted =
459 : tuples::get<QueueTags::InverseJacobian<Frame::Grid, Frame::Distorted>>(
460 : excision_quantities);
461 :
462 : db::mutate<gr::Tags::InverseSpatialMetric<DataVector, 3, Frame::Distorted>,
463 : ylm::Tags::Strahlkorper<Frame::Distorted>>(
464 : [&inverse_spatial_metric_on_excision, &excision_surface](
465 : const gsl::not_null<tnsr::II<DataVector, 3, Frame::Distorted>*>
466 : inv_spatial_metric,
467 : const gsl::not_null<ylm::Strahlkorper<Frame::Distorted>*>
468 : strahlkorper) {
469 : for (size_t i = 0; i < inv_spatial_metric->size(); i++) {
470 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
471 : (*inv_spatial_metric)[i].set_data_ref(const_cast<DataVector*>(
472 : &inverse_spatial_metric_on_excision[i]));
473 : }
474 : *strahlkorper = excision_surface;
475 : },
476 : make_not_null(&char_speed_box_));
477 :
478 : const double Y00 = 0.25 * M_2_SQRTPI;
479 :
480 : horizon_coef_averager_.update(time, {apparent_horizon.coefficients()[0]},
481 : smoother_tuner_.current_timescale());
482 :
483 : const std::optional<std::array<DataVector, DerivOrder + 1>>&
484 : averaged_horizon_coef_at_average_time = horizon_coef_averager_(time);
485 :
486 : // lambda_00 is the quantity of the same name in ArXiv:1211.6079,
487 : // and dt_lambda_00 is its time derivative.
488 : // This is the map parameter that maps the excision boundary in the grid
489 : // frame to the excision boundary in the distorted frame.
490 : const auto map_lambda_and_deriv =
491 : functions_of_time.at(function_of_time_name)->func_and_deriv(time);
492 : const double lambda_00 = map_lambda_and_deriv[0][0];
493 : const double dt_lambda_00 = map_lambda_and_deriv[1][0];
494 :
495 : // horizon_00 is \hat{S}_00 in ArXiv:1211.6079,
496 : // and dt_horizon_00 is its time derivative.
497 : // These are coefficients of the horizon in the distorted frame. However, we
498 : // want them averaged
499 : double horizon_00 = apparent_horizon.coefficients()[0];
500 : double dt_horizon_00 = time_deriv_apparent_horizon.coefficients()[0];
501 :
502 : // Only for the first few measurements will this not have a value
503 : if (LIKELY(averaged_horizon_coef_at_average_time.has_value())) {
504 : // We need to get the averaged time and evaluate the averaged coefs at
505 : // that time, not the time passed in
506 : const double averaged_time = horizon_coef_averager_.average_time(time);
507 :
508 : horizon_00 = averaged_horizon_coef_at_average_time.value()[0][0];
509 : dt_horizon_00 = averaged_horizon_coef_at_average_time.value()[1][0];
510 : const double d2t_horizon_00 =
511 : averaged_horizon_coef_at_average_time.value()[2][0];
512 :
513 : // Taylor expand from the averaged time t_avg to the current time t,
514 : // where time_diff = t - t_avg:
515 : // h(t) = h(t_avg) + time_diff * h'(t_avg)
516 : // + 0.5 * time_diff^2 * h''(t_avg),
517 : // h'(t) = h'(t_avg) + time_diff * h''(t_avg).
518 : const double time_diff = time - averaged_time;
519 : horizon_00 +=
520 : time_diff * dt_horizon_00 + 0.5 * square(time_diff) * d2t_horizon_00;
521 : dt_horizon_00 += time_diff * d2t_horizon_00;
522 :
523 : // The "control error" for the averaged horizon coefficients is just the
524 : // averaged coefs minus the actual coef and time derivative from
525 : // apparent_horizon and time_deriv_apparent_horizon
526 : smoother_tuner_.update_timescale(
527 : std::array{
528 : DataVector{averaged_horizon_coef_at_average_time.value()[0][0]},
529 : DataVector{averaged_horizon_coef_at_average_time.value()[1][0]}} -
530 : std::array{
531 : DataVector{apparent_horizon.coefficients()[0]},
532 : DataVector{time_deriv_apparent_horizon.coefficients()[0]}});
533 : }
534 :
535 : // This is needed because the horizon_00 (and dt) are spherepack coefs, not
536 : // spherical harmonic coefs.
537 : const double spherepack_factor = sqrt(0.5 * M_PI);
538 :
539 : // This is needed for every state
540 : const double control_error_delta_r = size::control_error_delta_r(
541 : horizon_00, dt_horizon_00, lambda_00, dt_lambda_00,
542 : grid_frame_excision_sphere_radius);
543 : const std::optional<double> control_error_delta_r_outward =
544 : delta_r_drift_outward_options_.has_value()
545 : ? std::optional<double>(control_error_delta_r -
546 : delta_r_drift_outward_options_.value()
547 : .outward_drift_velocity -
548 : (lambda_00 +
549 : spherepack_factor * horizon_00 -
550 : grid_frame_excision_sphere_radius / Y00) /
551 : delta_r_drift_outward_options_.value()
552 : .outward_drift_timescale)
553 : : std::nullopt;
554 : const std::optional<double> inward_drift_velocity =
555 : delta_r_drift_inward_options_.has_value()
556 : ? std::optional<double>(
557 : delta_r_drift_inward_options_.value().inward_drift_velocity)
558 : : std::nullopt;
559 : const std::optional<double> min_allowed_radial_distance =
560 : delta_r_drift_inward_options_.has_value()
561 : ? std::optional<double>(delta_r_drift_inward_options_.value()
562 : .min_allowed_radial_distance)
563 : : std::nullopt;
564 : const std::optional<double> min_allowed_char_speed =
565 : delta_r_drift_inward_options_.has_value()
566 : ? std::optional<double>(
567 : delta_r_drift_inward_options_.value().min_allowed_char_speed)
568 : : std::nullopt;
569 :
570 : // Currently we don't do anything with the derivative of the comoving char
571 : // speed. Eventually, we will pass it to the computation of the control
572 : // error below
573 : Scalar<DataVector> deriv_comoving_char_speed{};
574 : {
575 : const tnsr::i<DataVector, 3, Frame::Distorted>& excision_normal_one_form =
576 : db::get<ylm::Tags::NormalOneForm<Frame::Distorted>>(char_speed_box_);
577 : const DataVector& one_over_excision_normal_one_form_norm_dv =
578 : db::get<ylm::Tags::OneOverOneFormMagnitude>(char_speed_box_);
579 : Scalar<DataVector> one_over_excision_normal_one_form_norm{};
580 : get(one_over_excision_normal_one_form_norm)
581 : // NOLINTNEXTLINE
582 : .set_data_ref(const_cast<DataVector*>(
583 : &one_over_excision_normal_one_form_norm_dv));
584 :
585 : const tnsr::i<DataVector, 3, Frame::Distorted>& rhat =
586 : db::get<ylm::Tags::Rhat<Frame::Distorted>>(char_speed_box_);
587 :
588 : size::comoving_char_speed_derivative(
589 : make_not_null(&deriv_comoving_char_speed), lambda_00, dt_lambda_00,
590 : horizon_00, dt_horizon_00, grid_frame_excision_sphere_radius, rhat,
591 : excision_normal_one_form, one_over_excision_normal_one_form_norm,
592 : shifty_quantity, inverse_spatial_metric_on_excision,
593 : spatial_christoffel, deriv_lapse, deriv_shift,
594 : inv_jac_grid_to_distorted);
595 : }
596 :
597 : info_.damping_time = min(tuner.current_timescale());
598 :
599 : const size::ErrorDiagnostics error_diagnostics = size::control_error(
600 : make_not_null(&info_), make_not_null(&char_speed_predictor_),
601 : make_not_null(&comoving_char_speed_predictor_),
602 : make_not_null(&delta_radius_predictor_),
603 : make_not_null(&drift_limit_char_speed_predictor_),
604 : make_not_null(&drift_limit_delta_radius_predictor_), time,
605 : control_error_delta_r, control_error_delta_r_outward,
606 : delta_r_drift_outward_options_.has_value()
607 : ? std::optional<double>(delta_r_drift_outward_options_.value()
608 : .max_allowed_radial_distance)
609 : : std::nullopt,
610 : inward_drift_velocity, min_allowed_radial_distance,
611 : min_allowed_char_speed, horizon_00, dt_lambda_00, apparent_horizon,
612 : excision_surface, lapse, shifty_quantity, spatial_metric_on_excision,
613 : inverse_spatial_metric_on_excision, deriv_comoving_char_speed);
614 :
615 : state_history_.store(time, info_, error_diagnostics.control_error_args);
616 :
617 : if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
618 : auto& observer_writer_proxy = Parallel::get_parallel_component<
619 : observers::ObserverWriter<Metavariables>>(cache);
620 :
621 : // \Delta R = < R_ah > - < R_ex >
622 : // < R_ah > = S_00 * Y_00
623 : // < R_ex > = R_ex^grid - \lambda_00 * Y_00
624 : // < \Delta R > = \Delta R / < R_ah >
625 : const double avg_delta_r =
626 : (spherepack_factor * horizon_00 + lambda_00) * Y00 -
627 : grid_frame_excision_sphere_radius;
628 : const double avg_relative_delta_r =
629 : avg_delta_r / (spherepack_factor * horizon_00 * Y00);
630 :
631 : Parallel::threaded_action<
632 : observers::ThreadedActions::WriteReductionDataRow>(
633 : observer_writer_proxy[0], subfile_name_, legend_,
634 : std::make_tuple(
635 : time, error_diagnostics.control_error,
636 : static_cast<double>(error_diagnostics.state_number),
637 : error_diagnostics.discontinuous_change_has_occurred ? 1.0 : 0.0,
638 : lambda_00, dt_lambda_00, horizon_00, dt_horizon_00,
639 : time_deriv_apparent_horizon.coefficients()[0],
640 : smoother_tuner_.current_timescale()[0],
641 : error_diagnostics.min_delta_r,
642 : error_diagnostics.min_relative_delta_r, avg_delta_r,
643 : avg_relative_delta_r,
644 : error_diagnostics.control_error_args.control_error_delta_r,
645 : error_diagnostics.target_char_speed,
646 : error_diagnostics.control_error_args.min_char_speed,
647 : error_diagnostics.min_comoving_char_speed,
648 : error_diagnostics.char_speed_crossing_time,
649 : error_diagnostics.comoving_char_speed_crossing_time,
650 : error_diagnostics.delta_r_crossing_time,
651 : error_diagnostics.suggested_timescale,
652 : error_diagnostics.damping_timescale));
653 : }
654 :
655 : if (Parallel::get<control_system::Tags::Verbosity>(cache) >=
656 : ::Verbosity::Verbose) {
657 : Parallel::printf("%s: %s\n", function_of_time_name,
658 : error_diagnostics.update_message);
659 : }
660 :
661 : return DataVector{1, error_diagnostics.control_error};
662 : }
663 :
664 : private:
665 0 : TimescaleTuner<true> smoother_tuner_{};
666 0 : Averager<DerivOrder> horizon_coef_averager_{};
667 0 : size::Info info_{};
668 0 : intrp::ZeroCrossingPredictor char_speed_predictor_{};
669 0 : intrp::ZeroCrossingPredictor comoving_char_speed_predictor_{};
670 0 : intrp::ZeroCrossingPredictor delta_radius_predictor_{};
671 0 : intrp::ZeroCrossingPredictor drift_limit_char_speed_predictor_{};
672 0 : intrp::ZeroCrossingPredictor drift_limit_delta_radius_predictor_{};
673 0 : size::StateHistory state_history_{};
674 0 : std::vector<std::string> legend_{};
675 0 : std::string subfile_name_{};
676 0 : std::optional<DeltaRDriftOutwardOptions> delta_r_drift_outward_options_{};
677 0 : std::optional<DeltaRDriftInwardOptions> delta_r_drift_inward_options_{};
678 : db::compute_databox_type<tmpl::list<
679 : gr::Tags::InverseSpatialMetric<DataVector, 3, Frame::Distorted>,
680 : ylm::Tags::Strahlkorper<Frame::Distorted>,
681 : ylm::Tags::ThetaPhiCompute<Frame::Distorted>,
682 : ylm::Tags::InvJacobianCompute<Frame::Distorted>,
683 : ylm::Tags::RadiusCompute<Frame::Distorted>,
684 : ylm::Tags::RhatCompute<Frame::Distorted>,
685 : ylm::Tags::DxRadiusCompute<Frame::Distorted>,
686 : ylm::Tags::NormalOneFormCompute<Frame::Distorted>,
687 : ylm::Tags::OneOverOneFormMagnitudeCompute<DataVector, 3,
688 : Frame::Distorted>>>
689 0 : char_speed_box_{};
690 : };
691 : } // namespace ControlErrors
692 : } // namespace control_system
|
