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