SpECTRE Documentation Coverage Report
Current view: top level - ControlSystem/ControlErrors - Size.hpp Hit Total Coverage
Commit: 7f96d57c707fe8ffa5ad32d026c0a60120e43e0c Lines: 9 57 15.8 %
Date: 2024-11-04 21:21:10
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14