Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <memory>
       9             : #include <string>
      10             : 
      11             : #include "DataStructures/DataVector.hpp"
      12             : #include "DataStructures/ModalVector.hpp"
      13             : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
      14             : #include "NumericalAlgorithms/Strahlkorper/InitialShape.hpp"
      15             : #include "Options/Context.hpp"
      16             : #include "Options/String.hpp"
      17             : #include "Utilities/ForceInline.hpp"
      18             : #include "Utilities/StdArrayHelpers.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : 
      21             : /// \cond
      22             : namespace PUP {
      23             : class er;
      24             : }  // namespace PUP
      25             : 
      26             : namespace Options {
      27             : template <typename... AlternativeLists>
      28             : struct Alternatives;
      29             : }  // namespace Options
      30             : /// \endcond
      31             : 
      32             : namespace ylm {
      33             : /// \ingroup SurfacesGroup
      34             : /// \brief A star-shaped surface expanded in spherical harmonics.
      35             : template <typename Frame>
      36           1 : class Strahlkorper {
      37             :  public:
      38           0 :   struct InitialL {
      39           0 :     using type = size_t;
      40           0 :     static constexpr Options::String help = {
      41             :         "Initial Strahlkorper resolution. Sets both l_max and m_max."};
      42             :   };
      43           0 :   struct InitialShapeOption {
      44           0 :     using type = std::unique_ptr<InitialShape<Frame>>;
      45           0 :     static constexpr Options::String help = {"Initial Strahlkorper shape."};
      46           0 :     static std::string name() { return "InitialShape"; }
      47             :   };
      48           0 :   using options = tmpl::list<InitialL, InitialShapeOption>;
      49             : 
      50           0 :   static constexpr Options::String help{
      51             :       "A star-shaped surface expressed as an expansion in spherical "
      52             :       "harmonics.\n"
      53             :       "A Strahlkorper can be constructed from options by specifying InitialL "
      54             :       "and an InitialShape. InitialL sets both l_max and m_max."};
      55             : 
      56             :   // Pup needs default constructor
      57           0 :   Strahlkorper() = default;
      58           0 :   Strahlkorper(const Strahlkorper&) = default;
      59           0 :   Strahlkorper(Strahlkorper&&) = default;
      60           0 :   Strahlkorper& operator=(const Strahlkorper&) = default;
      61           0 :   Strahlkorper& operator=(Strahlkorper&&) = default;
      62             : 
      63             :   /// Construct a sphere of radius `radius` with a given center.
      64           1 :   Strahlkorper(size_t l_max, size_t m_max, double radius,
      65             :                std::array<double, 3> center);
      66             : 
      67             :   /// Construct a sphere of radius `radius`, setting `m_max`=`l_max`.
      68           1 :   Strahlkorper(size_t l_max, double radius, std::array<double, 3> center)
      69             :       : Strahlkorper(l_max, l_max, radius, center) {}
      70             : 
      71             :   /// Construct from options: either from radius and center, or from file.
      72             :   /// These constructors handle both alternatives specified in the options.
      73           1 :   Strahlkorper(size_t l_max, const std::string& h5_filename,
      74             :                const std::string& subfile_name, double time,
      75             :                double time_epsilon, bool check_frame,
      76             :                const Options::Context& context);
      77             : 
      78             :   /// Construct from options from an initial shape.
      79             :   /// @{
      80           1 :   Strahlkorper(size_t initial_l, const InitialShape<Frame>& initial_shape,
      81             :                const Options::Context& context);
      82           1 :   Strahlkorper(size_t initial_l,
      83             :                const std::unique_ptr<InitialShape<Frame>>& initial_shape,
      84             :                const Options::Context& context);
      85             :   /// @}
      86             : 
      87             :   /// Construct a Strahlkorper from a DataVector containing the radius
      88             :   /// at the collocation points.
      89             :   ///
      90             :   /// \note The collocation points of the constructed Strahlkorper
      91             :   /// will not be exactly `radius_at_collocation_points`.  Instead,
      92             :   /// the constructed Strahlkorper will match the shape given by
      93             :   /// `radius_at_collocation_points` only to order (`l_max`,`m_max`).
      94             :   /// This is because the ylm::Spherepack representation of the
      95             :   /// Strahlkorper has more collocation points than spectral
      96             :   /// coefficients.  Specifically, `radius_at_collocation_points` has
      97             :   /// \f$(l_{\rm max} + 1) (2 m_{\rm max} + 1)\f$ degrees of freedom,
      98             :   /// but because there are only
      99             :   /// \f$m_{\rm max}^2+(l_{\rm max}-m_{\rm max})(2m_{\rm max}+1)\f$
     100             :   /// spectral coefficients, it is not possible to choose spectral
     101             :   /// coefficients to exactly match all points in
     102             :   /// `radius_at_collocation_points`.
     103           1 :   explicit Strahlkorper(size_t l_max, size_t m_max,
     104             :                         const DataVector& radius_at_collocation_points,
     105             :                         std::array<double, 3> center);
     106             : 
     107             :   /// \brief Construct a Strahlkorper from a `ModalVector` containing the
     108             :   /// spectral coefficients
     109             :   ///
     110             :   /// \details The spectral coefficients should be in the form defined by
     111             :   /// `ylm::Strahlkorper::coefficients() const`.
     112           1 :   explicit Strahlkorper(size_t l_max, size_t m_max,
     113             :                         const ModalVector& spectral_coefficients,
     114             :                         std::array<double, 3> center);
     115             : 
     116             :   /// Copies a Strahlkorper from another frame into this Strahlkorper.
     117             :   ///
     118             :   /// \note The `OtherFrame` is ignored.
     119             :   template <typename OtherFrame>
     120           1 :   explicit Strahlkorper(const Strahlkorper<OtherFrame>& another_strahlkorper)
     121             :       : l_max_(another_strahlkorper.l_max()),
     122             :         m_max_(another_strahlkorper.m_max()),
     123             :         ylm_(l_max_, m_max_),
     124             :         center_(another_strahlkorper.expansion_center()),
     125             :         strahlkorper_coefs_(another_strahlkorper.coefficients()) {}
     126             : 
     127             :   /// Prolong or restrict another surface to the given `l_max` and `m_max`.
     128             :   ///
     129             :   /// \note The `OtherFrame` is ignored.
     130             :   template <typename OtherFrame>
     131           1 :   Strahlkorper(size_t l_max, size_t m_max,
     132             :                const Strahlkorper<OtherFrame>& another_strahlkorper)
     133             :       : l_max_(l_max),
     134             :         m_max_(m_max),
     135             :         ylm_(l_max, m_max),
     136             :         center_(another_strahlkorper.expansion_center()),
     137             :         strahlkorper_coefs_(
     138             :             another_strahlkorper.ylm_spherepack().prolong_or_restrict(
     139             :                 another_strahlkorper.coefficients(), ylm_)) {}
     140             : 
     141             :   /// Construct a Strahlkorper from another Strahlkorper,
     142             :   /// but explicitly specifying the coefficients.
     143             :   /// Here coefficients are in the same storage scheme
     144             :   /// as the `coefficients()` member function returns.
     145             :   ///
     146             :   /// \note The `OtherFrame` is ignored.
     147             :   template <typename OtherFrame>
     148           1 :   Strahlkorper(DataVector coefs,
     149             :                const Strahlkorper<OtherFrame>& another_strahlkorper)
     150             :       : l_max_(another_strahlkorper.l_max()),
     151             :         m_max_(another_strahlkorper.m_max()),
     152             :         ylm_(another_strahlkorper.ylm_spherepack()),
     153             :         center_(another_strahlkorper.expansion_center()),
     154             :         strahlkorper_coefs_(std::move(coefs)) {
     155             :     ASSERT(
     156             :         strahlkorper_coefs_.size() ==
     157             :             another_strahlkorper.ylm_spherepack().spectral_size(),
     158             :         "Bad size " << strahlkorper_coefs_.size() << ", expected "
     159             :                     << another_strahlkorper.ylm_spherepack().spectral_size());
     160             :   }
     161             : 
     162             :   /// Serialization for Charm++
     163             :   // NOLINTNEXTLINE(google-runtime-references)
     164           1 :   void pup(PUP::er& p);
     165             : 
     166             :   /*!
     167             :    *  These coefficients are stored as SPHEREPACK coefficients.
     168             :    *  Suppose you represent a set of coefficients \f$F^{lm}\f$ in the expansion
     169             :    *  \f[
     170             :    *  f(\theta,\phi) =
     171             :    *  \sum_{l=0}^{l_{max}} \sum_{m=-l}^{l} F^{lm} Y^{lm}(\theta,\phi)
     172             :    *  \f]
     173             :    *  Here the \f$Y^{lm}(\theta,\phi)\f$ are the usual complex-valued scalar
     174             :    *  spherical harmonics, so \f$F^{lm}\f$ are also complex-valued.
     175             :    *  But here we assume that \f$f(\theta,\phi)\f$ is real, so therefore
     176             :    *  the \f$F^{lm}\f$ obey \f$F^{l-m} = (-1)^m (F^{lm})^\star\f$. So one
     177             :    *  does not need to store both real and imaginary parts for both positive
     178             :    *  and negative \f$m\f$, and the stored coefficients can all be real.
     179             :    *
     180             :    *  So the stored coefficients are:
     181             :    * \f{align}
     182             :    *  \text{coefficients()(l,m)} &= (-1)^m \sqrt{\frac{2}{\pi}}
     183             :    *     \Re(F^{lm}) \quad \text{for} \quad m\ge 0, \\
     184             :    *  \text{coefficients()(l,m)} &= (-1)^m \sqrt{\frac{2}{\pi}}
     185             :    *     \Im(F^{lm}) \quad \text{for} \quad m<0
     186             :    * \f}
     187             :    */
     188           1 :   SPECTRE_ALWAYS_INLINE const DataVector& coefficients() const {
     189             :     return strahlkorper_coefs_;
     190             :   }
     191           0 :   SPECTRE_ALWAYS_INLINE DataVector& coefficients() {
     192             :     return strahlkorper_coefs_;
     193             :   }
     194             : 
     195             :   /// Point about which the spectral basis of the Strahlkorper is expanded.
     196             :   /// The center is given in the frame in which the Strahlkorper is defined.
     197             :   /// This center must be somewhere inside the Strahlkorper, but in principle
     198             :   /// it can be anywhere.  See `physical_center()` for a different measure.
     199           1 :   SPECTRE_ALWAYS_INLINE const std::array<double, 3>& expansion_center() const {
     200             :     return center_;
     201             :   }
     202             : 
     203             :   /// Approximate physical center (determined by \f$l=1\f$ coefficients)
     204             :   /// Implementation of Eqs. (38)-(40) in \cite Hemberger2012jz
     205           1 :   std::array<double, 3> physical_center() const;
     206             : 
     207             :   /// Average radius of the surface (determined by \f$Y_{00}\f$ coefficient)
     208           1 :   double average_radius() const;
     209             : 
     210             :   /// Maximum \f$l\f$ in \f$Y_{lm}\f$ decomposition.
     211           1 :   SPECTRE_ALWAYS_INLINE size_t l_max() const { return l_max_; }
     212             : 
     213             :   /// Maximum \f$m\f$ in \f$Y_{lm}\f$ decomposition.
     214           1 :   SPECTRE_ALWAYS_INLINE size_t m_max() const { return m_max_; }
     215             : 
     216             :   /// Radius at a particular angle \f$(\theta,\phi)\f$.
     217             :   /// This is inefficient if done at multiple points many times.
     218             :   /// See ylm::Spherepack for alternative ways of computing this.
     219           1 :   double radius(double theta, double phi) const;
     220             : 
     221             :   /// Determine if a point `x` is contained inside the surface.
     222             :   /// The point must be given in Cartesian coordinates in the frame in
     223             :   /// which the Strahlkorper is defined.
     224             :   /// This is inefficient if done at multiple points many times.
     225           1 :   bool point_is_contained(const std::array<double, 3>& x) const;
     226             : 
     227           0 :   SPECTRE_ALWAYS_INLINE const ylm::Spherepack& ylm_spherepack() const {
     228             :     return ylm_;
     229             :   }
     230             : 
     231             :  private:
     232           0 :   size_t l_max_{2}, m_max_{2};
     233           0 :   ylm::Spherepack ylm_{2, 2};
     234           0 :   std::array<double, 3> center_{{0.0, 0.0, 0.0}};
     235           0 :   DataVector strahlkorper_coefs_ = DataVector(ylm_.spectral_size(), 0.0);
     236             : };
     237             : 
     238           0 : namespace OptionTags {
     239             : /// \ingroup OptionTagsGroup
     240             : /// \ingroup SurfacesGroup
     241             : /// The input file tag for a Strahlkorper.
     242             : template <typename Frame>
     243           1 : struct Strahlkorper {
     244           0 :   using type = ylm::Strahlkorper<Frame>;
     245           0 :   static constexpr Options::String help{"A star-shaped surface"};
     246             : };
     247             : }  // namespace OptionTags
     248             : 
     249             : template <typename Frame>
     250           0 : bool operator==(const Strahlkorper<Frame>& lhs,
     251             :                 const Strahlkorper<Frame>& rhs) {
     252             :   return lhs.l_max() == rhs.l_max() and lhs.m_max() == rhs.m_max() and
     253             :          lhs.expansion_center() == rhs.expansion_center() and
     254             :          lhs.coefficients() == rhs.coefficients();
     255             : }
     256             : 
     257             : template <typename Frame>
     258           0 : bool operator!=(const Strahlkorper<Frame>& lhs,
     259             :                 const Strahlkorper<Frame>& rhs) {
     260             :   return not(lhs == rhs);
     261             : }
     262             : }  // namespace ylm