SpECTRE Documentation Coverage Report
Current view: top level - ControlSystem/ControlErrors - Size.hpp Hit Total Coverage
Commit: 1afd040526f6585adf96eb1710263813d4f5ac0e Lines: 9 56 16.1 %
Date: 2024-05-23 14:31:37
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 :   static constexpr size_t expected_number_of_excisions = 1;
     170             : 
     171           0 :   using object_centers = domain::object_list<Horizon>;
     172             : 
     173           0 :   struct MaxNumTimesForZeroCrossingPredictor {
     174             :     // Int so we get proper bounds checking
     175           0 :     using type = int;
     176           0 :     static constexpr Options::String help{
     177             :         "The maximum number of times used to calculate the zero crossing of "
     178             :         "the char speeds."};
     179           0 :     static int lower_bound() { return 3; }
     180             :   };
     181             : 
     182           0 :   struct SmoothAvgTimescaleFraction {
     183           0 :     using type = double;
     184           0 :     static constexpr Options::String help{
     185             :         "Average timescale fraction for smoothing horizon measurements."};
     186             :   };
     187             : 
     188           0 :   struct SmootherTuner {
     189           0 :     using type = TimescaleTuner<true>;
     190           0 :     static constexpr Options::String help{
     191             :         "TimescaleTuner for smoothing horizon measurements."};
     192             :   };
     193             : 
     194           0 :   struct InitialState {
     195           0 :     using type = std::unique_ptr<size::State>;
     196           0 :     static constexpr Options::String help{"Initial state to start in."};
     197             :   };
     198             : 
     199           0 :   struct DeltaRDriftOutwardOptions {
     200           0 :     using type =
     201             :         Options::Auto<DeltaRDriftOutwardOptions, Options::AutoLabel::None>;
     202           0 :     static constexpr Options::String help{
     203             :         "Options for State DeltaRDriftOutward. Specify 'None' to disable State "
     204             :         "DeltaRDriftOutward."};
     205           0 :     struct MaxAllowedRadialDistance {
     206           0 :       using type = double;
     207           0 :       static constexpr Options::String help{
     208             :           "Drift excision boundary outward if distance from horizon to "
     209             :           "excision exceeds this."};
     210             :     };
     211           0 :     struct OutwardDriftVelocity {
     212           0 :       using type = double;
     213           0 :       static constexpr Options::String help{
     214             :           "Constant drift velocity term, if triggered by "
     215             :           "MaxAllowedRadialDistance."};
     216             :     };
     217           0 :     struct OutwardDriftTimescale {
     218           0 :       using type = double;
     219           0 :       static constexpr Options::String help{
     220             :           "Denominator in non-constant drift velocity term, if triggered by "
     221             :           "MaxAllowedRadialDistance."};
     222             :     };
     223           0 :     using options = tmpl::list<MaxAllowedRadialDistance, OutwardDriftVelocity,
     224             :                                OutwardDriftTimescale>;
     225           0 :     void pup(PUP::er& p) {
     226             :       p | max_allowed_radial_distance;
     227             :       p | outward_drift_velocity;
     228             :       p | outward_drift_timescale;
     229             :     }
     230             : 
     231           0 :     double max_allowed_radial_distance{};
     232           0 :     double outward_drift_velocity{};
     233           0 :     double outward_drift_timescale{};
     234             :   };
     235             : 
     236           0 :   using options = tmpl::list<MaxNumTimesForZeroCrossingPredictor,
     237             :                              SmoothAvgTimescaleFraction, SmootherTuner,
     238             :                              InitialState, DeltaRDriftOutwardOptions>;
     239           0 :   static constexpr Options::String help{
     240             :       "Computes the control error for size control. Will also write a "
     241             :       "diagnostics file if the control systems are allowed to write data to "
     242             :       "disk."};
     243             : 
     244           0 :   Size() = default;
     245             : 
     246             :   /*!
     247             :    * \brief Initializes the `intrp::ZeroCrossingPredictor`s and the horizon
     248             :    * smoothing `Averager` and `TimescaleTuner`.
     249             :    *
     250             :    * \details All `intrp::ZeroCrossingPredictor`s are initialized with a minimum
     251             :    * number of times 3 and a maximum number of times `max_times`. The internal
     252             :    * `control_system::size::Info::state` is initialized to
     253             :    * `control_system::size::States::Initial`. The smoothing `Averager` uses the
     254             :    * input average timescale fraction and always smooths the "0th" deriv (aka
     255             :    * the horizon coefficients themselves). The input smoothing `TimescaleTuner`
     256             :    * is moved inside this class.
     257             :    */
     258           1 :   Size(const int max_times, const double smooth_avg_timescale_frac,
     259             :        TimescaleTuner<true> smoother_tuner,
     260             :        std::unique_ptr<size::State> initial_state,
     261             :        std::optional<DeltaRDriftOutwardOptions> delta_r_drift_outward_options);
     262             : 
     263             :   /// Returns the internal `control_system::size::Info::suggested_time_scale`. A
     264             :   /// std::nullopt means that no timescale is suggested.
     265           1 :   const std::optional<double>& get_suggested_timescale() const;
     266             : 
     267             :   /*!
     268             :    * \brief Check if the `control_system::size::control_error` has decided to
     269             :    * switch states. Returns the internal
     270             :    * `control_system::size::Info::discontinuous_change_has_occurred`.
     271             :    */
     272           1 :   bool discontinuous_change_has_occurred() const;
     273             : 
     274             :   /*!
     275             :    * \brief Reset the internal `control_system::size::Info` using
     276             :    * `control_system::size::Info::reset`.
     277             :    */
     278           1 :   void reset();
     279             : 
     280             :   /*!
     281             :    * \brief Get a history of the control errors for the past few measurements.
     282             :    *
     283             :    * \return std::deque<std::pair<double, double>> This returns up to
     284             :    * `DerivOrder` entries, not including the most recent time. \see
     285             :    * `control_system::size::StateHistory::state_history`
     286             :    */
     287           1 :   std::deque<std::pair<double, double>> control_error_history() const;
     288             : 
     289           0 :   void pup(PUP::er& p);
     290             : 
     291             :   /*!
     292             :    * \brief Actually computes the control error.
     293             :    *
     294             :    * \details The internal `control_system::size::Info::damping_time` is updated
     295             :    * to the minimum of the `TimescaleTuner::current_timescale()` that is passed
     296             :    * in. Also expects these queue tags to be in the `measurements` argument:
     297             :    *
     298             :    * - `ylm::Tags::Strahlkorper<Frame::Distorted>`
     299             :    * - `QueueTags::ExcisionSurface<Frame::Distorted>`
     300             :    * - `::Tags::dt<ylm::Tags::Strahlkorper<Frame::Distorted>>`
     301             :    * - `QueueTags::LapseOnExcisionSurface`
     302             :    * - `QueueTags::ShiftyQuantity<Frame::Distorted>`
     303             :    * - `QueueTags::SpatialMetricOnExcisionSurface<Frame::Distorted>`
     304             :    * - `QueueTags::InverseSpatialMetricOnExcisionSurface<Frame::Distorted>`
     305             :    *
     306             :    * \return DataVector should be of size 1
     307             :    */
     308             :   template <typename Metavariables, typename... TupleTags>
     309           1 :   DataVector operator()(const ::TimescaleTuner<false>& tuner,
     310             :                         const Parallel::GlobalCache<Metavariables>& cache,
     311             :                         const double time,
     312             :                         const std::string& function_of_time_name,
     313             :                         const tuples::TaggedTuple<TupleTags...>& measurements) {
     314             :     const Domain<3>& domain = get<domain::Tags::Domain<3>>(cache);
     315             :     const auto& excision_spheres = domain.excision_spheres();
     316             :     const auto& excision_sphere =
     317             :         excision_spheres.at("ExcisionSphere" + get_output(Horizon));
     318             :     const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
     319             : 
     320             :     const auto& excision_quantities =
     321             :         tuples::get<QueueTags::SizeExcisionQuantities<Frame::Distorted>>(
     322             :             measurements);
     323             :     const auto& horizon_quantities =
     324             :         tuples::get<QueueTags::SizeHorizonQuantities<Frame::Distorted>>(
     325             :             measurements);
     326             : 
     327             :     const double grid_frame_excision_sphere_radius = excision_sphere.radius();
     328             :     const ylm::Strahlkorper<Frame::Distorted>& apparent_horizon =
     329             :         tuples::get<ylm::Tags::Strahlkorper<Frame::Distorted>>(
     330             :             horizon_quantities);
     331             :     const ylm::Strahlkorper<Frame::Distorted>& excision_surface =
     332             :         tuples::get<QueueTags::ExcisionSurface<Frame::Distorted>>(
     333             :             excision_quantities);
     334             :     const ylm::Strahlkorper<Frame::Distorted>& time_deriv_apparent_horizon =
     335             :         tuples::get<::Tags::dt<ylm::Tags::Strahlkorper<Frame::Distorted>>>(
     336             :             horizon_quantities);
     337             :     const Scalar<DataVector>& lapse =
     338             :         tuples::get<QueueTags::LapseOnExcisionSurface>(excision_quantities);
     339             :     const tnsr::I<DataVector, 3, Frame::Distorted>& shifty_quantity =
     340             :         tuples::get<QueueTags::ShiftyQuantity<Frame::Distorted>>(
     341             :             excision_quantities);
     342             :     const tnsr::ii<DataVector, 3, Frame::Distorted>&
     343             :         spatial_metric_on_excision = tuples::get<
     344             :             QueueTags::SpatialMetricOnExcisionSurface<Frame::Distorted>>(
     345             :             excision_quantities);
     346             :     const tnsr::II<DataVector, 3, Frame::Distorted>&
     347             :         inverse_spatial_metric_on_excision = tuples::get<
     348             :             QueueTags::InverseSpatialMetricOnExcisionSurface<Frame::Distorted>>(
     349             :             excision_quantities);
     350             : 
     351             :     const double Y00 = 0.25 * M_2_SQRTPI;
     352             : 
     353             :     horizon_coef_averager_.update(time, {apparent_horizon.coefficients()[0]},
     354             :                                   smoother_tuner_.current_timescale());
     355             : 
     356             :     const std::optional<std::array<DataVector, DerivOrder + 1>>&
     357             :         averaged_horizon_coef_at_average_time = horizon_coef_averager_(time);
     358             : 
     359             :     // lambda_00 is the quantity of the same name in ArXiv:1211.6079,
     360             :     // and dt_lambda_00 is its time derivative.
     361             :     // This is the map parameter that maps the excision boundary in the grid
     362             :     // frame to the excision boundary in the distorted frame.
     363             :     const auto map_lambda_and_deriv =
     364             :         functions_of_time.at(function_of_time_name)->func_and_deriv(time);
     365             :     const double lambda_00 = map_lambda_and_deriv[0][0];
     366             :     const double dt_lambda_00 = map_lambda_and_deriv[1][0];
     367             : 
     368             :     // horizon_00 is \hat{S}_00 in ArXiv:1211.6079,
     369             :     // and dt_horizon_00 is its time derivative.
     370             :     // These are coefficients of the horizon in the distorted frame. However, we
     371             :     // want them averaged
     372             :     double horizon_00 = apparent_horizon.coefficients()[0];
     373             :     double dt_horizon_00 = time_deriv_apparent_horizon.coefficients()[0];
     374             : 
     375             :     // Only for the first few measurements will this not have a value
     376             :     if (LIKELY(averaged_horizon_coef_at_average_time.has_value())) {
     377             :       // We need to get the averaged time and evaluate the averaged coefs at
     378             :       // that time, not the time passed in
     379             :       const double averaged_time = horizon_coef_averager_.average_time(time);
     380             : 
     381             :       horizon_00 = averaged_horizon_coef_at_average_time.value()[0][0];
     382             :       dt_horizon_00 = averaged_horizon_coef_at_average_time.value()[1][0];
     383             :       const double d2t_horizon_00 =
     384             :           averaged_horizon_coef_at_average_time.value()[2][0];
     385             : 
     386             :       // Must account for time offset of averaged time. Do a simple Taylor
     387             :       // expansion
     388             :       const double time_diff = time - averaged_time;
     389             :       horizon_00 += time_diff * dt_horizon_00;
     390             :       dt_horizon_00 += 0.5 * square(time_diff) * d2t_horizon_00;
     391             : 
     392             :       // The "control error" for the averaged horizon coefficients is just the
     393             :       // averaged coefs minus the actual coef and time derivative from
     394             :       // apparent_horizon and time_deriv_apparent_horizon
     395             :       smoother_tuner_.update_timescale(
     396             :           std::array{
     397             :               DataVector{averaged_horizon_coef_at_average_time.value()[0][0]},
     398             :               DataVector{averaged_horizon_coef_at_average_time.value()[1][0]}} -
     399             :           std::array{
     400             :               DataVector{apparent_horizon.coefficients()[0]},
     401             :               DataVector{time_deriv_apparent_horizon.coefficients()[0]}});
     402             :     }
     403             : 
     404             :     // This is needed because the horizon_00 (and dt) are spherepack coefs, not
     405             :     // spherical harmonic coefs.
     406             :     const double spherepack_factor = sqrt(0.5 * M_PI);
     407             : 
     408             :     // This is needed for every state
     409             :     const double control_error_delta_r = size::control_error_delta_r(
     410             :         horizon_00, dt_horizon_00, lambda_00, dt_lambda_00,
     411             :         grid_frame_excision_sphere_radius);
     412             :     const std::optional<double> control_error_delta_r_outward =
     413             :         delta_r_drift_outward_options_.has_value()
     414             :             ? std::optional<double>(control_error_delta_r -
     415             :                                     delta_r_drift_outward_options_.value()
     416             :                                         .outward_drift_velocity -
     417             :                                     (lambda_00 +
     418             :                                      spherepack_factor * horizon_00 -
     419             :                                      grid_frame_excision_sphere_radius / Y00) /
     420             :                                         delta_r_drift_outward_options_.value()
     421             :                                             .outward_drift_timescale)
     422             :             : std::nullopt;
     423             : 
     424             :     info_.damping_time = min(tuner.current_timescale());
     425             : 
     426             :     const size::ErrorDiagnostics error_diagnostics = size::control_error(
     427             :         make_not_null(&info_), make_not_null(&char_speed_predictor_),
     428             :         make_not_null(&comoving_char_speed_predictor_),
     429             :         make_not_null(&delta_radius_predictor_), time, control_error_delta_r,
     430             :         control_error_delta_r_outward,
     431             :         delta_r_drift_outward_options_.has_value()
     432             :             ? std::optional<double>(delta_r_drift_outward_options_.value()
     433             :                                         .max_allowed_radial_distance)
     434             :             : std::nullopt,
     435             :         dt_lambda_00, apparent_horizon, excision_surface, lapse,
     436             :         shifty_quantity, spatial_metric_on_excision,
     437             :         inverse_spatial_metric_on_excision);
     438             : 
     439             :     state_history_.store(time, info_, error_diagnostics.control_error_args);
     440             : 
     441             :     if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
     442             :       auto& observer_writer_proxy = Parallel::get_parallel_component<
     443             :           observers::ObserverWriter<Metavariables>>(cache);
     444             : 
     445             :       // \Delta R = < R_ah > - < R_ex >
     446             :       // < R_ah > = S_00 * Y_00
     447             :       // < R_ex > = R_ex^grid - \lambda_00 * Y_00
     448             :       // < \Delta R > = \Delta R / < R_ah >
     449             :       const double avg_delta_r =
     450             :           (spherepack_factor * horizon_00 + lambda_00) * Y00 -
     451             :           grid_frame_excision_sphere_radius;
     452             :       const double avg_relative_delta_r =
     453             :           avg_delta_r / (spherepack_factor * horizon_00 * Y00);
     454             : 
     455             :       Parallel::threaded_action<
     456             :           observers::ThreadedActions::WriteReductionDataRow>(
     457             :           observer_writer_proxy[0], subfile_name_, legend_,
     458             :           std::make_tuple(
     459             :               time, error_diagnostics.control_error,
     460             :               static_cast<double>(error_diagnostics.state_number),
     461             :               error_diagnostics.discontinuous_change_has_occurred ? 1.0 : 0.0,
     462             :               lambda_00, dt_lambda_00, horizon_00, dt_horizon_00,
     463             :               time_deriv_apparent_horizon.coefficients()[0],
     464             :               smoother_tuner_.current_timescale()[0],
     465             :               error_diagnostics.min_delta_r,
     466             :               error_diagnostics.min_relative_delta_r, avg_delta_r,
     467             :               avg_relative_delta_r,
     468             :               error_diagnostics.control_error_args.control_error_delta_r,
     469             :               error_diagnostics.target_char_speed,
     470             :               error_diagnostics.control_error_args.min_char_speed,
     471             :               error_diagnostics.min_comoving_char_speed,
     472             :               error_diagnostics.char_speed_crossing_time,
     473             :               error_diagnostics.comoving_char_speed_crossing_time,
     474             :               error_diagnostics.delta_r_crossing_time,
     475             :               error_diagnostics.suggested_timescale,
     476             :               error_diagnostics.damping_timescale));
     477             :     }
     478             : 
     479             :     if (Parallel::get<control_system::Tags::Verbosity>(cache) >=
     480             :         ::Verbosity::Verbose) {
     481             :       Parallel::printf("%s: %s\n", function_of_time_name,
     482             :                        error_diagnostics.update_message);
     483             :     }
     484             : 
     485             :     return DataVector{1, error_diagnostics.control_error};
     486             :   }
     487             : 
     488             :  private:
     489           0 :   TimescaleTuner<true> smoother_tuner_{};
     490           0 :   Averager<DerivOrder> horizon_coef_averager_{};
     491           0 :   size::Info info_{};
     492           0 :   intrp::ZeroCrossingPredictor char_speed_predictor_{};
     493           0 :   intrp::ZeroCrossingPredictor comoving_char_speed_predictor_{};
     494           0 :   intrp::ZeroCrossingPredictor delta_radius_predictor_{};
     495           0 :   size::StateHistory state_history_{};
     496           0 :   std::vector<std::string> legend_{};
     497           0 :   std::string subfile_name_{};
     498           0 :   std::optional<DeltaRDriftOutwardOptions> delta_r_drift_outward_options_{};
     499             : };
     500             : }  // namespace ControlErrors
     501             : }  // namespace control_system

Generated by: LCOV version 1.14