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