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