Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cmath>
       7             : #include <cstddef>
       8             : #include <pup.h>
       9             : #include <string>
      10             : 
      11             : #include "ControlSystem/Protocols/ControlError.hpp"
      12             : #include "ControlSystem/Tags/QueueTags.hpp"
      13             : #include "ControlSystem/Tags/SystemTags.hpp"
      14             : #include "DataStructures/DataVector.hpp"
      15             : #include "Domain/Structure/ObjectLabel.hpp"
      16             : #include "NumericalAlgorithms/SphericalHarmonics/SpherepackIterator.hpp"
      17             : #include "Options/String.hpp"
      18             : #include "Parallel/GlobalCache.hpp"
      19             : #include "Utilities/EqualWithinRoundoff.hpp"
      20             : #include "Utilities/ErrorHandling/Assert.hpp"
      21             : #include "Utilities/ForceInline.hpp"
      22             : #include "Utilities/ProtocolHelpers.hpp"
      23             : #include "Utilities/TMPL.hpp"
      24             : #include "Utilities/TaggedTuple.hpp"
      25             : 
      26             : /// \cond
      27             : namespace domain::Tags {
      28             : template <size_t Dim>
      29             : struct Domain;
      30             : struct FunctionsOfTime;
      31             : }  // namespace domain::Tags
      32             : namespace Frame {
      33             : struct Distorted;
      34             : }  // namespace Frame
      35             : /// \endcond
      36             : 
      37             : namespace control_system::ControlErrors {
      38             : namespace detail {
      39             : template <::domain::ObjectLabel Horizon>
      40             : std::string excision_sphere_name() {
      41             :   return "ExcisionSphere"s + ::domain::name(Horizon);
      42             : }
      43             : 
      44             : template <::domain::ObjectLabel Horizon>
      45             : std::string size_name() {
      46             :   return "Size"s + ::domain::name(Horizon);
      47             : }
      48             : }  // namespace detail
      49             : 
      50             : /*!
      51             :  * \brief Control error in the
      52             :  * \link domain::CoordinateMaps::TimeDependent::Shape Shape \endlink coordinate
      53             :  * map
      54             :  *
      55             :  * \details Computes the error for each \f$ l,m \f$ mode of the shape map except
      56             :  * for \f$ l=0 \f$ and \f$ l=1 \f$. This is because the \f$ l=0 \f$ mode is
      57             :  * controlled by characteristic speed control, and the \f$ l=1 \f$ mode is
      58             :  * redundant with the \link
      59             :  * domain::CoordinateMaps::TimeDependent::Translation Translation \endlink map.
      60             :  * The equation for the control error is Eq. (77) from \cite Hemberger2012jz
      61             :  * which is
      62             :  *
      63             :  * \f{align}
      64             :  * Q_{lm} &= -\frac{r_\mathrm{EB} -
      65             :  * Y_{00}\lambda_{00}(t)}{\sqrt{\frac{\pi}{2}} Y_{00}S_{00}} S_{lm} -
      66             :  * \lambda_{lm}(t), \quad l>=2 \label{eq:shape_control_error}
      67             :  * \f}
      68             :  *
      69             :  * where \f$ r_\mathrm{EB} \f$ is the radius of the excision boundary in the
      70             :  * grid frame, \f$ \lambda_{00}(t) \f$ is the size map parameter, \f$
      71             :  * \lambda_{lm}(t) \f$ for \f$ l>=2 \f$ are the shape map parameters, and \f$
      72             :  * S_{lm}\f$ are the coefficients of the harmonic expansion of the apparent
      73             :  * horizon. The coefficients \f$ \lambda_{lm}(t) \f$ (*not* including \f$ l=0
      74             :  * \f$) and \f$ S_{lm}\f$ (including \f$ l=0 \f$) are stored as the real-valued
      75             :  * coefficients \f$ a_{lm} \f$ and \f$ b_{lm} \f$ of the SPHEREPACK
      76             :  * spherical-harmonic expansion (in ylm::Spherepack). The $\lambda_{00}(t)$
      77             :  * coefficient, on the other hand, is stored as the complex coefficient \f$
      78             :  * A_{00} \f$ of the standard \f$ Y_{lm} \f$ decomposition. Because \f$ a_{00} =
      79             :  * \sqrt{\frac{2}{\pi}}A_{00} \f$, there is an extra factor of \f$
      80             :  * \sqrt{\frac{\pi}{2}} \f$ in the above formula in the denominator of the
      81             :  * fraction multiplying the \f$ S_{00}\f$ component so it is represented in the
      82             :  * \f$ Y_{lm} \f$ decomposition just like \f$ r_{EB} \f$ and \f$ \lambda_{00}
      83             :  * \f$ are). That way, we ensure the numerator and denominator are represented
      84             :  * in the same way before we take their ratio.
      85             :  *
      86             :  * Requirements:
      87             :  * - This control error requires that there be at least one excision surface in
      88             :  *   the simulation
      89             :  * - Currently this control error can only be used with the \link
      90             :  *   control_system::Systems::Shape Shape \endlink control system
      91             :  */
      92             : template <::domain::ObjectLabel Object>
      93           1 : struct Shape : tt::ConformsTo<protocols::ControlError> {
      94             :   // Shape needs an excision sphere
      95           0 :   using object_centers = domain::object_list<Object>;
      96             : 
      97           0 :   using options = tmpl::list<>;
      98           0 :   static constexpr Options::String help{
      99             :       "Computes the control error for shape control. This should not take any "
     100             :       "options."};
     101             : 
     102             :   // NOLINTNEXTLINE(readability-convert-member-functions-to-static)
     103           0 :   std::optional<double> get_suggested_timescale() const { return std::nullopt; }
     104             : 
     105           0 :   void reset() {}
     106             : 
     107           0 :   void pup(PUP::er& /*p*/) {}
     108             : 
     109             :   template <typename Metavariables, typename... TupleTags>
     110           0 :   DataVector operator()(const ::TimescaleTuner<true>& /*unused*/,
     111             :                         const Parallel::GlobalCache<Metavariables>& cache,
     112             :                         const double time,
     113             :                         const std::string& function_of_time_name,
     114             :                         const tuples::TaggedTuple<TupleTags...>& measurements) {
     115             :     const auto& domain = get<domain::Tags::Domain<3>>(cache);
     116             :     const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
     117             :     const DataVector lambda_lm_coefs =
     118             :         functions_of_time.at(function_of_time_name)->func(time)[0];
     119             :     const double lambda_00_coef =
     120             :         functions_of_time.at(detail::size_name<Object>())->func(time)[0][0];
     121             : 
     122             :     const auto& ah =
     123             :         get<control_system::QueueTags::Horizon<Frame::Distorted, Object>>(
     124             :             measurements);
     125             :     const auto& ah_coefs = ah.coefficients();
     126             : 
     127             :     ASSERT(lambda_lm_coefs.size() == ah_coefs.size(),
     128             :            "Number of coefficients for shape map '"
     129             :                << function_of_time_name << "' (" << lambda_lm_coefs.size()
     130             :                << ") does not match the number of coefficients in the AH "
     131             :                   "Strahlkorper ("
     132             :                << ah_coefs.size() << ").");
     133             : 
     134             :     const auto& excision_spheres = domain.excision_spheres();
     135             : 
     136             :     ASSERT(domain.excision_spheres().count(
     137             :                detail::excision_sphere_name<Object>()) == 1,
     138             :            "Excision sphere " << detail::excision_sphere_name<Object>()
     139             :                               << " not in the domain but is needed to "
     140             :                                  "compute Shape control error.");
     141             : 
     142             :     const double radius_excision_sphere_grid_frame =
     143             :         excision_spheres.at(detail::excision_sphere_name<Object>()).radius();
     144             : 
     145             :     const double Y00 = sqrt(0.25 / M_PI);
     146             :     ylm::SpherepackIterator iter{ah.l_max(), ah.m_max()};
     147             :     // See above docs for why we have the sqrt(pi/2) in the denominator
     148             :     const double relative_size_factor =
     149             :         (radius_excision_sphere_grid_frame / Y00 - lambda_00_coef) /
     150             :         (sqrt(0.5 * M_PI) * ah_coefs[iter.set(0, 0)()]);
     151             : 
     152             :     // The map parameters are in terms of SPHEREPACK coefficients (just like
     153             :     // strahlkorper coefficients), *not* spherical harmonic coefficients, thus
     154             :     // the control error for each l,m is in terms of SPHEREPACK coefficients
     155             :     // and no extra factors of sqrt(2/pi) are needed
     156             :     DataVector Q = -relative_size_factor * ah_coefs - lambda_lm_coefs;
     157             : 
     158             :     // Shape control is only for l > 1 so enforce that Q=0 for l=0,l=1. These
     159             :     // components of the control error won't be 0 automatically because the AH
     160             :     // can freely have nonzero l=0 and l=1 coefficients so we have to set the
     161             :     // control error components to be 0 manually.
     162             :     Q[iter.set(0, 0)()] = 0.0;
     163             :     Q[iter.set(1, -1)()] = 0.0;
     164             :     Q[iter.set(1, 0)()] = 0.0;
     165             :     Q[iter.set(1, 1)()] = 0.0;
     166             : 
     167             :     return Q;
     168             :   }
     169             : };
     170             : }  // namespace control_system::ControlErrors