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