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