Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <memory>
       7             : #include <optional>
       8             : 
       9             : #include "ControlSystem/ControlErrors/Size/State.hpp"
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Utilities/Gsl.hpp"
      12             : 
      13             : /// \cond
      14             : class DataVector;
      15             : namespace ylm {
      16             : template <typename Frame>
      17             : class Strahlkorper;
      18             : }  // namespace ylm
      19             : namespace domain::FunctionsOfTime {
      20             : class FunctionOfTime;
      21             : }  // namespace domain::FunctionsOfTime
      22             : namespace control_system::size {
      23             : class Info;
      24             : }  // namespace control_system::size
      25             : namespace intrp {
      26             : class ZeroCrossingPredictor;
      27             : }  // namespace intrp
      28             : /// \endcond
      29             : 
      30             : namespace control_system::size {
      31             : /*!
      32             :  * \brief A simple struct to hold diagnostic information about computing the
      33             :  * size control error.
      34             :  */
      35           1 : struct ErrorDiagnostics {
      36           0 :   double control_error;
      37           0 :   size_t state_number;
      38           0 :   double min_delta_r;
      39           0 :   double min_relative_delta_r;
      40           0 :   double min_comoving_char_speed;
      41           0 :   double char_speed_crossing_time;
      42           0 :   double comoving_char_speed_crossing_time;
      43           0 :   double delta_r_crossing_time;
      44           0 :   double target_char_speed;
      45           0 :   double suggested_timescale;
      46           0 :   double damping_timescale;
      47           0 :   ControlErrorArgs control_error_args;
      48           0 :   std::string update_message;
      49           0 :   bool discontinuous_change_has_occurred;
      50             : };
      51             : 
      52             : /*!
      53             :  * \brief Computes the size control error, updating the stored info.
      54             :  *
      55             :  * \tparam Frame should be ::Frame::Distorted if ::Frame::Distorted exists.
      56             :  * \param info struct containing parameters that will be used/filled. Some of
      57             :  *        the fields in info will guide the behavior of other control system
      58             :  *        components like the averager and (possibly) the time step.
      59             :  * \param predictor_char_speed ZeroCrossingPredictor for the characteristic
      60             :  *          speed.
      61             :  * \param predictor_comoving_char_speed ZeroCrossingPredictor for the
      62             :  *        comoving characteristic speed.
      63             :  * \param predictor_delta_radius ZeroCrossingPredictor for the difference
      64             :  *        in radius between the horizon and the excision boundary.
      65             :  * \param predictor_drift_limit_char_speed ZeroCrossingPredictor for the
      66             :  *        difference between the characteristic speed and
      67             :  *        min_allowed_char_speed.
      68             :  * \param predictor_drift_limit_delta_radius ZeroCrossingPredictor for the
      69             :  *        difference between delta_radius and min_allowed_radial_distance.
      70             :  * \param time the current time.
      71             :  * \param control_error_delta_r the control error for the DeltaR state. This is
      72             :  *        used in other states as well.
      73             :  * \param control_error_delta_r_outward the control error for the
      74             :  *        DeltaRDriftOutward state.  If std::nullopt, then DeltaRDriftOutward
      75             :  *        will not be used.
      76             :  * \param max_allowed_radial_distance the maximum average radial distance
      77             :  *        between the horizon and the excision boundary that is allowed without
      78             :  *        triggering the DeltaRDriftOutward state.  If std::nullopt, then
      79             :  *        DeltaRDriftOutward will not be used.
      80             :  * \param inward_drift_velocity a velocity that determines how fast the
      81             :  *        excision boundary drifts inward in state DeltaRDriftInward. If
      82             :  *        std::nullopt, then DeltaRDriftInward will not be used.
      83             :  * \param min_allowed_radial_distance the minimum average radial distance
      84             :  *        between the horizon and the excision boundary that is allowed without
      85             :  *        triggering the DeltaRDriftInward state.  If std::nullopt, then
      86             :  *        DeltaRDriftInward will not be used.
      87             :  * \param min_allowed_char_speed the minimum characteristic speed on the
      88             :  *        excision boundary that is allowed without triggering the
      89             :  *        DeltaRDriftInward state.  If std::nullopt, then
      90             :  *        DeltaRDriftInward will not be used.
      91             :  * \param horizon_00 The l=0,m=0 component of the spherepack decomposition
      92             :  *        of the apparent horizon.  This is passed separately from the
      93             :  *        full apparent_horizon below because horizon_00 is time-averaged.
      94             :  * \param dt_lambda_00 the time derivative of the map parameter lambda_00
      95             :  * \param apparent_horizon the current horizon in frame Frame.
      96             :  * \param excision_boundary a Strahlkorper representing the excision
      97             :  *        boundary in frame Frame.  Note that the excision boundary is assumed
      98             :  *        to be a sphere in the grid frame.
      99             :  * \param lapse_on_excision_boundary Lapse on the excision boundary.
     100             :  * \param frame_components_of_grid_shift The quantity
     101             :  *        \f$\beta^i \frac{\partial x^\hat{i}}{\partial x_i}\f$ (see below)
     102             :  *        evaluated on the excision boundary.  This is a tensor in frame
     103             :  *        Frame.
     104             :  * \param spatial_metric_on_excision_boundary metric in frame Frame.
     105             :  * \param inverse_spatial_metric_on_excision_boundary metric in frame Frame.
     106             :  * \param deriv_comoving_char_speed the derivative of the comoving char
     107             :  *        speed with respect to the map parameter lambda_00.
     108             :  * \return Returns an `ErrorDiagnostics` object which, in addition to the actual
     109             :  *         control error, holds a lot of diagnostic information about how the
     110             :  *         control error was calculated. This information could be used to print
     111             :  *         to a file if desired.
     112             :  *
     113             :  * The characteristic speed that is needed here is
     114             :  * \f{align}
     115             :  *     v &= -\alpha -n_i\beta^i \\
     116             :  *     v &= -\alpha -n_\hat{i}\hat{\beta}^\hat{i}
     117             :  *           - n_\hat{i}\frac{\partial x^\hat{i}}{\partial t} \\
     118             :  *     v &= -\alpha -n_\bar{i}\bar{\beta}^\bar{i}
     119             :  *           - n_\bar{i}\frac{\partial x^\bar{i}}{\partial t} \\
     120             :  *     v &= -\alpha - n_\hat{i}
     121             :  *          \frac{\partial x^\hat{i}}{\partial x^i} \beta^i,
     122             :  *  \f}
     123             :  *  where we have written many equivalent forms in terms of quantities
     124             :  *  defined in different frames.
     125             :  *
     126             :  *  Here \f$\alpha\f$ is the lapse, which is invariant under frame
     127             :  *  transformations, \f$n_i\f$, \f$n_\hat{i}\f$, and \f$n_\bar{i}\f$
     128             :  *  are the metric-normalized normal one-form to the Strahlkorper in the
     129             :  *  grid, distorted, and inertial frames, and
     130             :  *  \f$\beta^i\f$, \f$\hat{\beta}^\hat{i}\f$, and \f$\bar{\beta}^\bar{i}\f$
     131             :  *  are the shift in the grid, distorted, and inertial frames.
     132             :  *
     133             :  *  Note that we decorate the shift with hats and bars in addition to
     134             :  *  decorating its index, because the shift transforms in a non-obvious
     135             :  *  way under frame transformations so it is easy to make mistakes.
     136             :  *  To be clear, these different shifts are defined by
     137             :  * \f{align}
     138             :  *   \beta^i &= \alpha^2 g^{0i},\\
     139             :  *   \hat{\beta}^\hat{i} &= \alpha^2 g^{\hat{0}\hat{i}},\\
     140             :  *   \bar{\beta}^\bar{i} &= \alpha^2 g^{\bar{0}\bar{i}},
     141             :  * \f}
     142             :  *  where \f$g^{ab}\f$ is the spacetime metric, and they transform like
     143             :  * \f{align}
     144             :  * \hat{\beta}^\hat{i} &= \beta^i \frac{\partial x^\hat{i}}{\partial x^i}-
     145             :  *  \frac{\partial x^\hat{i}}{\partial t}.
     146             :  * \f}
     147             :  *
     148             :  * The quantity we pass as frame_components_of_grid_shift is
     149             :  * \f{align}
     150             :  * \beta^i \frac{\partial x^\hat{i}}{\partial x^i}
     151             :  * &= \hat{\beta}^\hat{i} + \frac{\partial x^\hat{i}}{\partial t} \\
     152             :  * &= \bar{\beta}^\bar{j}\frac{\partial x^\hat{i}}{\partial x^i}
     153             :  *    \frac{\partial x^i}{\partial x^\bar{j}} +
     154             :  *    \frac{\partial x^\hat{i}}{\partial x^\bar{j}}
     155             :  *    \frac{\partial x^\bar{j}}{\partial t},
     156             :  * \f}
     157             :  * where we have listed several equivalent formulas that involve quantities
     158             :  * in different frames.
     159             :  */
     160             : template <typename Frame>
     161           1 : ErrorDiagnostics control_error(
     162             :     gsl::not_null<Info*> info,
     163             :     gsl::not_null<intrp::ZeroCrossingPredictor*> predictor_char_speed,
     164             :     gsl::not_null<intrp::ZeroCrossingPredictor*>
     165             :         predictor_comoving_char_speed,
     166             :     gsl::not_null<intrp::ZeroCrossingPredictor*> predictor_delta_radius,
     167             :     gsl::not_null<intrp::ZeroCrossingPredictor*>
     168             :         predictor_drift_limit_char_speed,
     169             :     gsl::not_null<intrp::ZeroCrossingPredictor*>
     170             :         predictor_drift_limit_delta_radius,
     171             :     double time, double control_error_delta_r,
     172             :     std::optional<double> control_error_delta_r_outward,
     173             :     std::optional<double> max_allowed_radial_distance,
     174             :     std::optional<double> inward_drift_velocity,
     175             :     std::optional<double> min_allowed_radial_distance,
     176             :     std::optional<double> min_allowed_char_speed, double horizon_00,
     177             :     double dt_lambda_00, const ylm::Strahlkorper<Frame>& apparent_horizon,
     178             :     const ylm::Strahlkorper<Frame>& excision_boundary,
     179             :     const Scalar<DataVector>& lapse_on_excision_boundary,
     180             :     const tnsr::I<DataVector, 3, Frame>& frame_components_of_grid_shift,
     181             :     const tnsr::ii<DataVector, 3, Frame>& spatial_metric_on_excision_boundary,
     182             :     const tnsr::II<DataVector, 3, Frame>&
     183             :         inverse_spatial_metric_on_excision_boundary,
     184             :     const Scalar<DataVector>& deriv_comoving_char_speed);
     185             : 
     186             : }  // namespace control_system::size