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
|
