Strahlkorper.hpp
1 // 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 "ApparentHorizons/YlmSpherepack.hpp"
10 #include "DataStructures/DataVector.hpp"
11 #include "Options/Options.hpp"
13 
14 namespace PUP {
15 class er;
16 } // namespace PUP
17 
18 /// \ingroup SurfacesGroup
19 /// \brief A star-shaped surface expanded in spherical harmonics.
20 template <typename Frame>
21 class Strahlkorper {
22  public:
23  struct Lmax {
24  using type = size_t;
25  static constexpr OptionString help = {
26  "Strahlkorper is expanded in Ylms up to l=Lmax"};
27  };
28  struct Radius {
29  using type = double;
30  static constexpr OptionString help = {"Radius of spherical Strahlkorper"};
31  };
32  struct Center {
34  static constexpr OptionString help = {"Center of spherical Strahlkorper"};
35  };
36  using options = tmpl::list<Lmax, Radius, Center>;
37 
38  static constexpr OptionString help{
39  "A star-shaped surface expressed as an expansion in spherical "
40  "harmonics.\n"
41  "Currently only a spherical Strahlkorper can be constructed from\n"
42  "Options. To do this, specify parameters Center, Radius, and Lmax."};
43 
44  // Pup needs default constructor
45  Strahlkorper() = default;
46 
47  /// Construct a sphere of radius `radius` with a given center.
48  Strahlkorper(size_t l_max, size_t m_max, double radius,
49  std::array<double, 3> center) noexcept;
50 
51  /// Construct a sphere of radius `radius`, setting `m_max`=`l_max`.
52  Strahlkorper(size_t l_max, double radius,
53  std::array<double, 3> center) noexcept
54  : Strahlkorper(l_max, l_max, radius, center) {}
55 
56  /// Construct a Strahlkorper from a DataVector containing the radius
57  /// at the collocation points.
58  ///
59  /// \note The collocation points of the constructed Strahlkorper
60  /// will not be exactly `radius_at_collocation_points`. Instead,
61  /// the constructed Strahlkorper will match the shape given by
62  /// `radius_at_collocation_points` only to order (`l_max`,`m_max`).
63  /// This is because the YlmSpherepack representation of the
64  /// Strahlkorper has more collocation points than spectral
65  /// coefficients. Specifically, `radius_at_collocation_points` has
66  /// \f$(l_{\rm max} + 1) (2 m_{\rm max} + 1)\f$ degrees of freedom,
67  /// but because there are only
68  /// \f$m_{\rm max}^2+(l_{\rm max}-m_{\rm max})(2m_{\rm max}+1)\f$
69  /// spectral coefficients, it is not possible to choose spectral
70  /// coefficients to exactly match all points in
71  /// `radius_at_collocation_points`.
72  Strahlkorper(size_t l_max, size_t m_max,
73  const DataVector& radius_at_collocation_points,
74  std::array<double, 3> center) noexcept;
75 
76  /// Prolong or restrict another surface to the given `l_max` and `m_max`.
77  Strahlkorper(size_t l_max, size_t m_max,
78  const Strahlkorper& another_strahlkorper) noexcept;
79 
80  /// Construct a Strahlkorper from another Strahlkorper,
81  /// but explicitly specifying the coefficients.
82  /// Here coefficients are in the same storage scheme
83  /// as the `coefficients()` member function returns.
85  const Strahlkorper& another_strahlkorper) noexcept;
86 
87  /// Move-construct a Strahlkorper from another Strahlkorper,
88  /// explicitly specifying the coefficients.
89  Strahlkorper(DataVector coefs, Strahlkorper&& another_strahlkorper) noexcept;
90 
91  // clang-tidy: no runtime references
92  /// Serialization for Charm++
93  void pup(PUP::er& p) noexcept; // NOLINT
94 
95  /*!
96  * These coefficients are stored as SPHEREPACK coefficients.
97  * Suppose you represent a set of coefficients \f$F^{lm}\f$ in the expansion
98  * \f[
99  * f(\theta,\phi) =
100  * \sum_{l=0}^{l_{max}} \sum_{m=-l}^{l} F^{lm} Y^{lm}(\theta,\phi)
101  * \f]
102  * Here the \f$Y^{lm}(\theta,\phi)\f$ are the usual complex-valued scalar
103  * spherical harmonics, so \f$F^{lm}\f$ are also complex-valued.
104  * But here we assume that \f$f(\theta,\phi)\f$ is real, so therefore
105  * the \f$F^{lm}\f$ obey \f$F^{l-m} = (-1)^m (F^{lm})^\star\f$. So one
106  * does not need to store both real and imaginary parts for both positive
107  * and negative \f$m\f$, and the stored coefficients can all be real.
108  *
109  * So the stored coefficients are:
110  * \f{align}
111  * \text{coefficients()(l,m)} &= (-1)^m \sqrt{\frac{2}{\pi}}
112  * \Re(F^{lm}) \quad \text{for} \quad m\ge 0, \\
113  * \text{coefficients()(l,m)} &= (-1)^m \sqrt{\frac{2}{\pi}}
114  * \Im(F^{lm}) \quad \text{for} \quad m<0
115  * \f}
116  */
118  return strahlkorper_coefs_;
119  }
120  SPECTRE_ALWAYS_INLINE DataVector& coefficients() noexcept {
121  return strahlkorper_coefs_;
122  }
123 
124  /// Point about which the spectral basis of the Strahlkorper is expanded.
125  /// The center is given in the frame in which the Strahlkorper is defined.
126  /// This center must be somewhere inside the Strahlkorper, but in principle
127  /// it can be anywhere. See `physical_center()` for a different measure.
129  return center_;
130  }
131 
132  /// Approximate physical center (determined by \f$l=1\f$ coefficients)
133  /// Implementation of Eqs. (38)-(40) in \cite Hemberger2012jz
134  std::array<double, 3> physical_center() const noexcept;
135 
136  /// Average radius of the surface (determined by \f$Y_{00}\f$ coefficient)
137  double average_radius() const noexcept;
138 
139  /// Maximum \f$l\f$ in \f$Y_{lm}\f$ decomposition.
140  SPECTRE_ALWAYS_INLINE size_t l_max() const noexcept { return l_max_; }
141 
142  /// Maximum \f$m\f$ in \f$Y_{lm}\f$ decomposition.
143  SPECTRE_ALWAYS_INLINE size_t m_max() const noexcept { return m_max_; }
144 
145  /// Radius at a particular angle \f$(\theta,\phi)\f$.
146  /// This is inefficient if done at multiple points many times.
147  /// See YlmSpherepack for alternative ways of computing this.
148  double radius(double theta, double phi) const noexcept;
149 
150  /// Determine if a point `x` is contained inside the surface.
151  /// The point must be given in Cartesian coordinates in the frame in
152  /// which the Strahlkorper is defined.
153  /// This is inefficient if done at multiple points many times.
154  bool point_is_contained(const std::array<double, 3>& x) const noexcept;
155 
156  SPECTRE_ALWAYS_INLINE const YlmSpherepack& ylm_spherepack() const noexcept {
157  return ylm_;
158  }
159 
160  private:
161  size_t l_max_{2}, m_max_{2};
162  YlmSpherepack ylm_{2, 2};
163  std::array<double, 3> center_{{0.0, 0.0, 0.0}};
164  DataVector strahlkorper_coefs_ = DataVector(ylm_.spectral_size(), 0.0);
165 };
166 
167 namespace OptionTags {
168 /// \ingroup OptionTagsGroup
169 /// \ingroup SurfacesGroup
170 /// The input file tag for a Strahlkorper.
171 template <typename Frame>
172 struct Strahlkorper {
173  using type = ::Strahlkorper<Frame>;
174  static constexpr OptionString help{"A star-shaped surface"};
175 };
176 } // namespace OptionTags
177 
178 template <typename Frame>
179 bool operator==(const Strahlkorper<Frame>& lhs,
180  const Strahlkorper<Frame>& rhs) noexcept {
181  return lhs.l_max() == rhs.l_max() and lhs.m_max() == rhs.m_max() and
182  lhs.center() == rhs.center() and
183  lhs.coefficients() == rhs.coefficients();
184 }
185 
186 template <typename Frame>
187 bool operator!=(const Strahlkorper<Frame>& lhs,
188  const Strahlkorper<Frame>& rhs) noexcept {
189  return not(lhs == rhs);
190 }
Definition: Strahlkorper.hpp:14
Definition: Strahlkorper.hpp:32
Defines classes and functions for making classes creatable from input files.
size_t l_max() const noexcept
Maximum in decomposition.
Definition: Strahlkorper.hpp:140
C++ interface to SPHEREPACK.
Definition: YlmSpherepack.hpp:21
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
The input file tag for a Strahlkorper.
Definition: Strahlkorper.hpp:172
size_t m_max() const noexcept
Maximum in decomposition.
Definition: Strahlkorper.hpp:143
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
tnsr::iaa< DataType, SpatialDim, Frame > phi(const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein&#39;s equations...
Definition: ComputeGhQuantities.cpp:22
Definition: Strahlkorper.hpp:23
Definition: Strahlkorper.hpp:167
Defines macro to always inline a function.
Definition: Strahlkorper.hpp:28
A star-shaped surface expanded in spherical harmonics.
Definition: Strahlkorper.hpp:21
Stores a collection of function values.
Definition: DataVector.hpp:46
const DataVector & coefficients() const noexcept
Definition: Strahlkorper.hpp:117
const std::array< double, 3 > & center() const noexcept
Point about which the spectral basis of the Strahlkorper is expanded. The center is given in the fram...
Definition: Strahlkorper.hpp:128
Strahlkorper(size_t l_max, double radius, std::array< double, 3 > center) noexcept
Construct a sphere of radius radius, setting m_max=l_max.
Definition: Strahlkorper.hpp:52