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 <limits>
9 : #include <optional>
10 : #include <string>
11 : #include <variant>
12 :
13 : #include "DataStructures/DataVector.hpp"
14 : #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp"
15 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
16 : #include "Domain/Structure/ObjectLabel.hpp"
17 : #include "NumericalAlgorithms/SphericalHarmonics/IO/ReadSurfaceYlm.hpp"
18 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
19 : #include "Options/Auto.hpp"
20 : #include "Options/Context.hpp"
21 : #include "Options/String.hpp"
22 : #include "Utilities/Gsl.hpp"
23 : #include "Utilities/TMPL.hpp"
24 :
25 : namespace domain::creators::time_dependent_options {
26 : /*!
27 : * \brief Mass and spin necessary for calculating the \f$ Y_{lm} \f$
28 : * coefficients of a Kerr horizon of certain Boyer-Lindquist radius for the
29 : * shape map of the Sphere domain creator.
30 : */
31 1 : struct KerrSchildFromBoyerLindquist {
32 : /// \brief The mass of the Kerr black hole.
33 1 : struct Mass {
34 0 : using type = double;
35 0 : static constexpr Options::String help = {"The mass of the Kerr BH."};
36 : };
37 : /// \brief The dimensionless spin of the Kerr black hole.
38 1 : struct Spin {
39 0 : using type = std::array<double, 3>;
40 0 : static constexpr Options::String help = {
41 : "The dim'less spin of the Kerr BH."};
42 : };
43 :
44 0 : using options = tmpl::list<Mass, Spin>;
45 :
46 0 : static constexpr Options::String help = {
47 : "Conform to an ellipsoid of constant Boyer-Lindquist radius in "
48 : "Kerr-Schild coordinates. This Boyer-Lindquist radius is chosen as the "
49 : "value of the 'InnerRadius'. To conform to the outer Kerr horizon, "
50 : "choose an 'InnerRadius' of r_+ = M + sqrt(M^2-a^2)."};
51 :
52 0 : KerrSchildFromBoyerLindquist();
53 0 : KerrSchildFromBoyerLindquist(double mass_in, std::array<double, 3> spin_in);
54 :
55 0 : double mass{std::numeric_limits<double>::signaling_NaN()};
56 0 : std::array<double, 3> spin{std::numeric_limits<double>::signaling_NaN(),
57 : std::numeric_limits<double>::signaling_NaN(),
58 : std::numeric_limits<double>::signaling_NaN()};
59 : };
60 :
61 : /// Label for shape map options
62 1 : struct Spherical {};
63 :
64 0 : struct YlmsFromFile {
65 0 : struct H5Filename {
66 0 : using type = std::string;
67 0 : static constexpr Options::String help =
68 : "Path to the data file containing the ylm coefficients and their "
69 : "derivatives.";
70 : };
71 :
72 0 : struct SubfileNames {
73 0 : using type = std::vector<std::string>;
74 0 : static constexpr Options::String help =
75 : "Subfile names for the different order derivatives of the ylm "
76 : "coefficients. You must specify the subfile name for the ylm "
77 : "coefficients themselves, and can optionally specify the subfile name "
78 : "for the first and second time derivatives as well, in that order. If "
79 : "you don't specify a derivative subfile, those coefficients will be "
80 : "defaulted to zero.";
81 0 : static size_t lower_bound_on_size() { return 1; }
82 0 : static size_t upper_bound_on_size() { return 3; }
83 : };
84 :
85 0 : struct MatchTime {
86 0 : using type = double;
87 0 : static constexpr Options::String help =
88 : "Time in the H5File to get the coefficients at. Will likely be the "
89 : "same as the initial time";
90 : };
91 :
92 0 : struct MatchTimeEpsilon {
93 0 : using type = Options::Auto<double>;
94 0 : static constexpr Options::String help =
95 : "Look for times in the H5File within this epsilon of the match time. "
96 : "This is to avoid having to know the exact time to all digits. Default "
97 : "is 1e-12.";
98 : };
99 :
100 0 : struct SetL1CoefsToZero {
101 0 : using type = bool;
102 0 : static constexpr Options::String help =
103 : "Whether to set the L=1 coefs to zero or not. This may be desirable "
104 : "because L=1 is degenerate with a translation of the BH.";
105 : };
106 :
107 0 : struct CheckFrame {
108 0 : using type = bool;
109 0 : static constexpr Options::String help =
110 : "Whether to check if the frame of the Strahlkorper in the file matches "
111 : "the Distorted frame.";
112 : };
113 :
114 0 : using options = tmpl::list<H5Filename, SubfileNames, MatchTime,
115 : MatchTimeEpsilon, SetL1CoefsToZero, CheckFrame>;
116 :
117 0 : static constexpr Options::String help = {
118 : "Read the Y_lm coefficients of a Strahlkorper from file and use those to "
119 : "initialize the coefficients of a shape map."};
120 0 : YlmsFromFile();
121 0 : YlmsFromFile(std::string h5_filename_in,
122 : std::vector<std::string> subfile_names_in, double match_time_in,
123 : std::optional<double> match_time_epsilon_in,
124 : bool set_l1_coefs_to_zero_in, bool check_frame_in = true);
125 :
126 0 : std::string h5_filename;
127 0 : std::vector<std::string> subfile_names;
128 0 : double match_time{};
129 0 : std::optional<double> match_time_epsilon;
130 0 : bool set_l1_coefs_to_zero{};
131 0 : bool check_frame{true};
132 : };
133 :
134 0 : struct YlmsFromSpEC {
135 0 : struct DatFilename {
136 0 : using type = std::string;
137 0 : static constexpr Options::String help =
138 : "Name of the SpEC dat file holding the coefficients. Note that this "
139 : "isn't a Dat file within an H5 file. This must be an actual '.dat' "
140 : "file on disk.";
141 : };
142 :
143 0 : struct MatchTime {
144 0 : using type = double;
145 0 : static constexpr Options::String help =
146 : "Time in the H5File to get the coefficients at. Will likely be the "
147 : "same as the initial time";
148 : };
149 :
150 0 : struct MatchTimeEpsilon {
151 0 : using type = Options::Auto<double>;
152 0 : static constexpr Options::String help =
153 : "Look for times in the H5File within this epsilon of the match time. "
154 : "This is to avoid having to know the exact time to all digits. Default "
155 : "is 1e-12.";
156 : };
157 :
158 0 : struct SetL1CoefsToZero {
159 0 : using type = bool;
160 0 : static constexpr Options::String help =
161 : "Whether to set the L=1 coefs to zero or not. This may be desirable "
162 : "because L=1 is degenerate with a translation of the BH.";
163 : };
164 :
165 0 : using options =
166 : tmpl::list<DatFilename, MatchTime, MatchTimeEpsilon, SetL1CoefsToZero>;
167 :
168 0 : static constexpr Options::String help = {
169 : "Read the Y_lm coefficients of a Strahlkorper from file and use those to "
170 : "initialize the coefficients of a shape map."};
171 0 : YlmsFromSpEC();
172 0 : YlmsFromSpEC(std::string dat_filename_in, double match_time_in,
173 : std::optional<double> match_time_epsilon_in,
174 : bool set_l1_coefs_to_zero_in);
175 :
176 0 : std::string dat_filename;
177 0 : double match_time{};
178 0 : std::optional<double> match_time_epsilon;
179 0 : bool set_l1_coefs_to_zero{};
180 : };
181 :
182 : namespace detail {
183 : struct TransitionEndsAtCube {
184 : using type = bool;
185 : static constexpr Options::String help = {
186 : "If 'true', the shape map transition function will be 0 at the cubical "
187 : "boundary around the object. If 'false' the transition function will "
188 : "be 0 at the outer radius of the inner sphere around the object"};
189 : };
190 : } // namespace detail
191 :
192 : /*!
193 : * \brief Specialized version of `FromVolumeFile` for the shape map
194 : *
195 : * \details This is needed because the regular `FromVolumeFile` doesn't have
196 : * options for domain settings like `TransitionEndsAtCube` or `LMax`.
197 : */
198 : template <ObjectLabel Object>
199 1 : struct FromVolumeFileShapeSize : public FromVolumeFile {
200 : public:
201 0 : struct LMax {
202 0 : using type = Options::Auto<size_t>;
203 0 : static constexpr Options::String help = {
204 : "LMax used for the number of spherical harmonic coefficients of the "
205 : "distortion map. If set to 'Auto', will use the LMax from the shape "
206 : "function of time in the volume file."};
207 : };
208 0 : using options = tmpl::push_front<FromVolumeFile::options, LMax,
209 : detail::TransitionEndsAtCube>;
210 :
211 0 : FromVolumeFileShapeSize() = default;
212 0 : FromVolumeFileShapeSize(const std::optional<size_t>& l_max_in,
213 : bool transition_ends_at_cube_in,
214 : std::string h5_filename, std::string subfile_name,
215 : const Options::Context& context = {});
216 :
217 0 : size_t l_max{};
218 0 : bool transition_ends_at_cube{};
219 :
220 : private:
221 0 : std::string h5_filename_;
222 0 : std::string subfile_name_;
223 : };
224 :
225 : /*!
226 : * \brief Class to be used as an option for initializing shape map coefficients.
227 : *
228 : * \details This class can also be used as an option tag with the \p type type
229 : * alias, `name()` function, and \p help string.
230 : *
231 : * \tparam IncludeTransitionEndsAtCube This is mainly added for the
232 : * `domain::creators::BinaryCompactObject` domain.
233 : * \tparam Object Which object that this shape map represents. Use
234 : * `domain::ObjectLabel::None` if there is only a single object in your
235 : * simulation.
236 : */
237 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
238 1 : struct ShapeMapOptions {
239 0 : using type = Options::Auto<
240 : std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
241 : FromVolumeFileShapeSize<Object>>,
242 : Options::AutoLabel::None>;
243 0 : static std::string name() { return "ShapeMap" + get_output(Object); }
244 0 : static constexpr Options::String help = {
245 : "Options for a time-dependent distortion (shape) map about the "
246 : "specified object. Specify 'None' to not use this map."};
247 :
248 0 : struct LMax {
249 0 : using type = size_t;
250 0 : static constexpr Options::String help = {
251 : "LMax used for the number of spherical harmonic coefficients of the "
252 : "distortion map."};
253 : };
254 :
255 0 : struct InitialValues {
256 0 : using type = Options::Auto<
257 : std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile, YlmsFromSpEC>,
258 : Spherical>;
259 0 : static constexpr Options::String help = {
260 : "Initial Ylm coefficients for the shape map. Specify 'Spherical' for "
261 : "all coefficients to be initialized to zero."};
262 : };
263 :
264 0 : struct SizeInitialValues {
265 0 : using type = Options::Auto<std::array<double, 3>>;
266 0 : static constexpr Options::String help = {
267 : "Initial value and two derivatives of the 00 coefficient. Specify "
268 : "'Auto' to use the 00 coefficient specified in the 'InitialValues' "
269 : "option. If you specify 'Auto', the deformed sphere will match the "
270 : "'InitialValues' surface exactly, and the original radius will only "
271 : "set the radius of the sphere in the grid frame (before deformation)."};
272 : };
273 :
274 0 : using common_options = tmpl::list<LMax, InitialValues, SizeInitialValues>;
275 :
276 0 : using options = tmpl::conditional_t<
277 : IncludeTransitionEndsAtCube,
278 : tmpl::push_back<common_options, detail::TransitionEndsAtCube>,
279 : common_options>;
280 0 : ShapeMapOptions() = default;
281 0 : ShapeMapOptions(size_t l_max_in,
282 : std::optional<std::variant<KerrSchildFromBoyerLindquist,
283 : YlmsFromFile, YlmsFromSpEC>>
284 : initial_values_in,
285 : std::optional<std::array<double, 3>> initial_size_values_in =
286 : std::nullopt,
287 : bool transition_ends_at_cube_in = false)
288 : : l_max(l_max_in),
289 : initial_values(std::move(initial_values_in)),
290 : initial_size_values(initial_size_values_in),
291 : transition_ends_at_cube(transition_ends_at_cube_in) {}
292 :
293 0 : size_t l_max{};
294 : std::optional<
295 : std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile, YlmsFromSpEC>>
296 0 : initial_values;
297 0 : std::optional<std::array<double, 3>> initial_size_values;
298 0 : bool transition_ends_at_cube{false};
299 : };
300 :
301 : /*!
302 : * \brief Helper function to get LMax from the different variants that the shape
303 : * map options could be.
304 : */
305 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
306 1 : size_t l_max_from_shape_options(
307 : const std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
308 : FromVolumeFileShapeSize<Object>>& shape_map_options);
309 :
310 : /*!
311 : * \brief Helper function to get whether the shape map transition function ends
312 : * at the cube from the different variants that the shape map options could be.
313 : */
314 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
315 1 : bool transition_ends_at_cube_from_shape_options(
316 : const std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
317 : FromVolumeFileShapeSize<Object>>& shape_map_options);
318 :
319 : /*!
320 : * \brief Helper function that takes the variant of the shape map options, and
321 : * returns the fully constructed shape and size functions of time.
322 : *
323 : * \details Even if the functions of time are read from a file, they will have a
324 : * new \p initial_time, \p shape_expiration_time, and \p size_expiration_time.
325 : * The \p deformed_radius is only used for the non-volume file variants.
326 : */
327 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
328 1 : FunctionsOfTimeMap get_shape_and_size(
329 : const std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
330 : FromVolumeFileShapeSize<Object>>& shape_map_options,
331 : double initial_time, double shape_expiration_time,
332 : double size_expiration_time, double deformed_radius);
333 : } // namespace domain::creators::time_dependent_options
|