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/Strahlkorper/IO/ReadSurfaceYlm.hpp"
18 : #include "NumericalAlgorithms/Strahlkorper/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 CoefficientTruncationLimit {
202 0 : using type = double;
203 0 : static constexpr Options::String help = {
204 : "Coefficients below this absolute value will be truncated from the "
205 : "Shape map. Set to 0.0 to disable truncation."};
206 0 : static constexpr type default_value = 0.0;
207 : };
208 0 : struct LMax {
209 0 : using type = Options::Auto<size_t>;
210 0 : static constexpr Options::String help = {
211 : "LMax used for the number of spherical harmonic coefficients of the "
212 : "distortion map. If set to 'Auto', will use the LMax from the shape "
213 : "function of time in the volume file."};
214 : };
215 0 : using options = tmpl::push_front<FromVolumeFile::options, LMax,
216 : CoefficientTruncationLimit,
217 : detail::TransitionEndsAtCube>;
218 :
219 0 : FromVolumeFileShapeSize() = default;
220 0 : FromVolumeFileShapeSize(const std::optional<size_t>& l_max_in,
221 : double coefficient_truncation_limit_in,
222 : bool transition_ends_at_cube_in,
223 : std::string h5_filename, std::string subfile_name,
224 : const Options::Context& context = {});
225 :
226 0 : std::optional<size_t> l_max{};
227 0 : double coefficient_truncation_limit{0.0};
228 0 : bool transition_ends_at_cube{};
229 :
230 : private:
231 0 : std::string h5_filename_;
232 0 : std::string subfile_name_;
233 : };
234 :
235 : /*!
236 : * \brief Class to be used as an option for initializing shape map coefficients.
237 : *
238 : * \details This class can also be used as an option tag with the \p type type
239 : * alias, `name()` function, and \p help string.
240 : *
241 : * \tparam IncludeTransitionEndsAtCube This is mainly added for the
242 : * `domain::creators::BinaryCompactObject` domain.
243 : * \tparam Object Which object that this shape map represents. Use
244 : * `domain::ObjectLabel::None` if there is only a single object in your
245 : * simulation.
246 : */
247 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
248 1 : struct ShapeMapOptions {
249 0 : struct CoefficientTruncationLimit {
250 0 : using type = double;
251 0 : static constexpr Options::String help = {
252 : "Coefficients below this absolute value will be truncated from the "
253 : "Shape map. Set to 0.0 to disable truncation."};
254 0 : static constexpr type default_value = 0.0;
255 : };
256 0 : using type = Options::Auto<
257 : std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
258 : FromVolumeFileShapeSize<Object>>,
259 : Options::AutoLabel::None>;
260 0 : static std::string name() { return "ShapeMap" + get_output(Object); }
261 0 : static constexpr Options::String help = {
262 : "Options for a time-dependent distortion (shape) map about the "
263 : "specified object. Specify 'None' to not use this map."};
264 :
265 0 : struct LMax {
266 0 : using type = size_t;
267 0 : static constexpr Options::String help = {
268 : "LMax used for the number of spherical harmonic coefficients of the "
269 : "distortion map."};
270 : };
271 :
272 0 : struct InitialValues {
273 0 : using type = Options::Auto<
274 : std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile, YlmsFromSpEC>,
275 : Spherical>;
276 0 : static constexpr Options::String help = {
277 : "Initial Ylm coefficients for the shape map. Specify 'Spherical' for "
278 : "all coefficients to be initialized to zero."};
279 : };
280 :
281 0 : struct SizeInitialValues {
282 0 : using type = Options::Auto<std::array<double, 3>>;
283 0 : static constexpr Options::String help = {
284 : "Initial value and two derivatives of the 00 coefficient. Specify "
285 : "'Auto' to use the 00 coefficient specified in the 'InitialValues' "
286 : "option. If you specify 'Auto', the deformed sphere will match the "
287 : "'InitialValues' surface exactly, and the original radius will only "
288 : "set the radius of the sphere in the grid frame (before deformation)."};
289 : };
290 :
291 0 : using common_options = tmpl::list<LMax, InitialValues, SizeInitialValues,
292 : CoefficientTruncationLimit>;
293 :
294 0 : using options = tmpl::conditional_t<
295 : IncludeTransitionEndsAtCube,
296 : tmpl::push_back<common_options, detail::TransitionEndsAtCube>,
297 : common_options>;
298 0 : ShapeMapOptions() = default;
299 0 : ShapeMapOptions(size_t l_max_in,
300 : std::optional<std::variant<KerrSchildFromBoyerLindquist,
301 : YlmsFromFile, YlmsFromSpEC>>
302 : initial_values_in,
303 : std::optional<std::array<double, 3>> initial_size_values_in =
304 : std::nullopt,
305 : double coefficient_truncation_limit_in = 0.0,
306 : const Options::Context& context = {});
307 0 : ShapeMapOptions(size_t l_max_in,
308 : std::optional<std::variant<KerrSchildFromBoyerLindquist,
309 : YlmsFromFile, YlmsFromSpEC>>
310 : initial_values_in,
311 : std::optional<std::array<double, 3>> initial_size_values_in,
312 : double coefficient_truncation_limit_in,
313 : bool transition_ends_at_cube_in,
314 : const Options::Context& context = {});
315 :
316 0 : size_t l_max{};
317 : std::optional<
318 : std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile, YlmsFromSpEC>>
319 0 : initial_values;
320 0 : std::optional<std::array<double, 3>> initial_size_values;
321 0 : double coefficient_truncation_limit{0.0};
322 0 : bool transition_ends_at_cube{false};
323 : };
324 :
325 : /*!
326 : * \brief Helper function to get the coefficient truncation limit from the
327 : * different variants that the shape map options could be.
328 : */
329 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
330 1 : double coefficient_truncation_limit_from_shape_options(
331 : const std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
332 : FromVolumeFileShapeSize<Object>>& shape_map_options);
333 :
334 : /*!
335 : * \brief Helper function to get whether the shape map transition function ends
336 : * at the cube from the different variants that the shape map options could be.
337 : */
338 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
339 1 : bool transition_ends_at_cube_from_shape_options(
340 : const std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
341 : FromVolumeFileShapeSize<Object>>& shape_map_options);
342 :
343 : /*!
344 : * \brief Helper function that takes the variant of the shape map options, and
345 : * returns the fully constructed shape and size functions of time.
346 : *
347 : * \details Even if the functions of time are read from a file, they will have a
348 : * new \p initial_time, \p shape_expiration_time, and \p size_expiration_time.
349 : * The \p deformed_radius is only used for the non-volume file variants.
350 : */
351 : template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
352 1 : FunctionsOfTimeMap get_shape_and_size(
353 : const std::variant<ShapeMapOptions<IncludeTransitionEndsAtCube, Object>,
354 : FromVolumeFileShapeSize<Object>>& shape_map_options,
355 : double initial_time, double shape_expiration_time,
356 : double size_expiration_time, double deformed_radius);
357 : } // namespace domain::creators::time_dependent_options
|
