Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <pup.h>
8 : #include <string>
9 : #include <tuple>
10 : #include <type_traits>
11 : #include <unordered_set>
12 : #include <vector>
13 :
14 : #include "DataStructures/ComplexDataVector.hpp"
15 : #include "DataStructures/ComplexModalVector.hpp"
16 : #include "DataStructures/DataBox/DataBox.hpp"
17 : #include "DataStructures/DataBox/DataBoxTag.hpp"
18 : #include "DataStructures/DataBox/Prefixes.hpp"
19 : #include "Domain/Structure/ElementId.hpp"
20 : #include "Evolution/Systems/Cce/OptionTags.hpp"
21 : #include "Evolution/Systems/Cce/Tags.hpp"
22 : #include "IO/Observer/ObserverComponent.hpp"
23 : #include "IO/Observer/ReductionActions.hpp"
24 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCoefficients.hpp"
25 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCollocation.hpp"
26 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshTransform.hpp"
27 : #include "Options/Context.hpp"
28 : #include "Options/ParseError.hpp"
29 : #include "Options/String.hpp"
30 : #include "Parallel/GlobalCache.hpp"
31 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
32 : #include "Time/Tags/Time.hpp"
33 : #include "Utilities/Gsl.hpp"
34 : #include "Utilities/MakeString.hpp"
35 : #include "Utilities/TMPL.hpp"
36 : #include "Utilities/TypeTraits/IsA.hpp"
37 :
38 : /// \cond
39 : template <size_t Dim>
40 : class Mesh;
41 : namespace Frame {
42 : struct Inertial;
43 : } // namespace Frame
44 : /// \endcond
45 :
46 : namespace Cce {
47 0 : namespace Events {
48 : namespace detail {
49 : template <typename Tag>
50 : std::string name() {
51 : if constexpr (std::is_same_v<Tag, Tags::ComplexInertialRetardedTime>) {
52 : return db::tag_name<Tags::InertialRetardedTime>();
53 : } else {
54 : return db::tag_name<Tag>();
55 : }
56 : }
57 : } // namespace detail
58 :
59 : /*!
60 : * \brief Event to observe fields/variables in a characteristic evolution.
61 : *
62 : * \details Similar to `dg::Events::ObserveFields`, this event will write volume
63 : * data from the characteristic domain to disk when triggered. However, there
64 : * are several differences which are important to highlight.
65 : *
66 : * First is the fields themselves. The DG event takes the fields to observe as
67 : * template parameters because the event must work with many evolution systems.
68 : * However, since this event is specific to the characteristic evolution system,
69 : * we can hardcode the list of fields that are available to observe. The fields
70 : * available to observe are the following tags along with their first and second
71 : * `Cce::Tags::Dy` derivatives (see `Cce::Tags::Dy` for a definition of `y`):
72 : *
73 : * - `Cce::Tags::BondiBeta`
74 : * - `Cce::Tags::BondiU`
75 : * - `Cce::Tags::BondiQ`
76 : * - `Cce::Tags::BondiW`
77 : * - `Cce::Tags::BondiH` (no second derivative)
78 : * - `Cce::Tags::BondiJ`
79 : * - `Cce::Tags::Du<Cce::Tags::BondiJ>`
80 : *
81 : * Some more fields to observe are:
82 : *
83 : * - `Cce::Tags::ComplexInertialRetardedTime`
84 : * - `Cce::Tags::OneMinusY`
85 : * - `Cce::Tags::BondiR`
86 : * - `Cce::Tags::EthRDividedByR`
87 : * - `Cce::Tags::DuRDividedByR`
88 : *
89 : * The main reason that this event is separate from the DG one is because this
90 : * event writes modal data over the sphere for every radial grid point, while
91 : * the DG event writes nodal data. Every tag above is a
92 : * `Scalar<SpinWeighted<ComplexDataVector, Spin>>` for some `Spin`. While this
93 : * data itself is in nodal form, it is more convenient to transform to modal
94 : * data and decompose in spherical harmonics before writing. This means
95 : * our typical way of writing/storing volume data won't work.
96 : *
97 : * This event writes its data in the following structure in the H5 file:
98 : * `/Cce/VolumeData/TagName/CompactifiedRadius_X.dat`. Every field that is
99 : * observed will get its own subgroup called `TagName`. In this subgroup, there
100 : * will be N files corresponding to N radial grid points named
101 : * `CompactifiedRadius_X.dat` where `X` here will range from 0 to `N-1`. We call
102 : * these compactified radii because for a more "physical" radius, it goes to
103 : * infinity at future-null infinity and we can't write that in a file. Instead,
104 : * these N files will correspond to the compactified coordinate $y = 1 - 2R/r$
105 : * where $r$ is your coordinate radius and $R$ is the coordinate radius of your
106 : * worldtube. Each file will hold the modal data for that radial grid point. It
107 : * is recommended to always dump the quantity `Cce::Tags::OneMinusY` so the
108 : * values of the compactified coordinates are available as well.
109 : *
110 : * There are two notable exceptions to this format. One is
111 : * `Cce::Tags::ComplexInertialRetardedTime`. The quantity we are actually
112 : * interested in is `Cce::Tags::InertialRetardedTime` which is real and only
113 : * defined once for every direction $\theta,\phi$ (meaning it does not have
114 : * different values at the different radial grid points). However, we use
115 : * `Cce::Tags::ComplexInertialRetardedTime` because it has the same data type as
116 : * the other tags which makes the internals of the class simpler. The imaginary
117 : * part of this `ComplexDataVector` is set to zero. This quantity will be stored
118 : * in a subfile named `/Cce/VolumeData/InertialRetardedTime.dat` as a single
119 : * modal set of data so we don't repeat it N times.
120 : *
121 : * The second is `Cce::Tags::OneMinusY`. Even though this quantity is stored
122 : * as a `Scalar<SpinWeighted<ComplexDataVector, 0>>` like the others, there is
123 : * only one meaningful value per radial grid point. All angular grid points for
124 : * a given radius are set to this value, namely $1-y$. Thus we only need to
125 : * write this value once for each radial grid point. We do this in a subfile
126 : * `/Cce/VolumeData/OneMinusY.dat` where the columns are named
127 : * `CompactifiedRadius_X` corresponding to the radial subfiles written for the
128 : * spin weighted quantities above (and time as the first column).
129 : *
130 : * All data will be written into the `observers::OptionTags::ReductionFileName`
131 : * file.
132 : */
133 1 : class ObserveFields : public Event {
134 : template <typename Tag, bool IncludeSecondDeriv = true>
135 : // clang-format off
136 0 : using zero_one_two_radial_derivs = tmpl::flatten<tmpl::list<
137 : Tag,
138 : Tags::Dy<Tag>,
139 : tmpl::conditional_t<IncludeSecondDeriv,
140 : Tags::Dy<Tags::Dy<Tag>>,
141 : tmpl::list<>>>>;
142 0 : using spin_weighted_tags_to_observe = tmpl::flatten<
143 : tmpl::list<zero_one_two_radial_derivs<Tags::BondiBeta>,
144 : zero_one_two_radial_derivs<Tags::BondiU>,
145 : zero_one_two_radial_derivs<Tags::BondiQ>,
146 : zero_one_two_radial_derivs<Tags::BondiW>,
147 : zero_one_two_radial_derivs<Tags::BondiH, false>,
148 : zero_one_two_radial_derivs<Tags::BondiJ>,
149 : zero_one_two_radial_derivs<Tags::Du<Tags::BondiJ>>,
150 : Tags::BondiR,
151 : Tags::EthRDividedByR,
152 : Tags::DuRDividedByR>>;
153 : // clang-format on
154 :
155 : public:
156 0 : using available_tags_to_observe =
157 : tmpl::push_back<spin_weighted_tags_to_observe,
158 : Tags::ComplexInertialRetardedTime, Tags::OneMinusY>;
159 :
160 : /// \cond
161 : explicit ObserveFields(CkMigrateMessage* /*unused*/) {}
162 : using PUP::able::register_constructor;
163 : WRAPPED_PUPable_decl_template(ObserveFields); // NOLINT
164 : /// \endcond
165 :
166 0 : struct VariablesToObserve {
167 0 : static constexpr Options::String help = "Subset of variables to observe";
168 0 : using type = std::vector<std::string>;
169 0 : static size_t lower_bound_on_size() { return 1; }
170 : };
171 :
172 0 : using options = tmpl::list<VariablesToObserve>;
173 :
174 0 : static constexpr Options::String help =
175 : "Observe volume tensor fields on the characteristic grid. Writes volume "
176 : "quantities from the tensors listed in the 'VariablesToObserve' "
177 : "option to the `/Cce/VolumeData` subfile of the reduction h5 file.\n";
178 :
179 0 : ObserveFields() = default;
180 :
181 0 : ObserveFields(const std::vector<std::string>& variables_to_observe,
182 : const Options::Context& context = {});
183 :
184 0 : using compute_tags_for_observation_box = tmpl::list<>;
185 :
186 0 : using return_tags = tmpl::list<>;
187 0 : using argument_tags = tmpl::list<::Tags::DataBox>;
188 :
189 : template <typename DbTags, typename Metavariables, typename ArrayIndex,
190 : typename ParallelComponent>
191 0 : void operator()(const db::DataBox<DbTags>& box,
192 : Parallel::GlobalCache<Metavariables>& cache,
193 : const ArrayIndex& /*array_index*/,
194 : const ParallelComponent* const /*component*/,
195 : const ObservationValue& /*observation_value*/) const {
196 : // Number of points
197 : const size_t l_max = db::get<Tags::LMax>(box);
198 : const size_t l_max_plus_one_squared = square(l_max + 1);
199 : const size_t number_of_angular_points =
200 : Spectral::Swsh::number_of_swsh_collocation_points(l_max);
201 : const size_t number_of_radial_grid_points =
202 : db::get<Tags::NumberOfRadialPoints>(box);
203 :
204 : // Buffers/views
205 : std::vector<double> data_to_write(2 * l_max_plus_one_squared + 1);
206 : ComplexModalVector goldberg_mode_buffer{l_max_plus_one_squared};
207 : ComplexDataVector spin_weighted_data_view{};
208 :
209 : // Legend
210 : std::vector<std::string> file_legend;
211 : file_legend.reserve(2 * l_max_plus_one_squared + 1);
212 : file_legend.emplace_back("time");
213 : for (int i = 0; i <= static_cast<int>(l_max); ++i) {
214 : for (int j = -i; j <= i; ++j) {
215 : file_legend.push_back(MakeString{} << "Real Y_" << i << "," << j);
216 : file_legend.push_back(MakeString{} << "Imag Y_" << i << "," << j);
217 : }
218 : }
219 :
220 : // Time
221 : const double time = db::get<::Tags::Time>(box);
222 :
223 : // Observer writer
224 : auto observer_proxy = Parallel::get_parallel_component<
225 : ::observers::ObserverWriter<Metavariables>>(cache)[0];
226 :
227 : // Actual work to transform nodal data to modal data. Places result in
228 : // data_to_write (but starts placing data in the 1st, not 0th, element
229 : // because the 0th element is time). Also makes use of the
230 : // spin_weighted_data_view and goldberg_mode_buffer
231 : const auto transform_to_modal =
232 : [&spin_weighted_data_view, &goldberg_mode_buffer, &box, &l_max,
233 : &number_of_angular_points, &l_max_plus_one_squared,
234 : &data_to_write](auto tag_v, const auto& spin_weighted_transform,
235 : auto& goldberg_modes, const size_t radial_index) {
236 : using tag = std::decay_t<decltype(tag_v)>;
237 :
238 : // Get ComplexDataVector out of SpinWeighted out of databox tag
239 : const ComplexDataVector& tensor = get(db::get<tag>(box)).data();
240 :
241 : // Make non-owning ComplexDataVector to angular data corresponding to
242 : // this radial index
243 : // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
244 : spin_weighted_data_view.set_data_ref(
245 : const_cast<ComplexDataVector&>(tensor).data() +
246 : radial_index * number_of_angular_points,
247 : number_of_angular_points);
248 :
249 : // swsh_transform requires a SpinWeighted<ComplexDataVector>. It's
250 : // easier to make a const-view from a ComplexDataVector that is
251 : // already the proper size (spin_weighted_data_view) than it is to try
252 : // and do all the indexing into the big block of memory here in this
253 : // call. That is why we have spin_weighted_data_view above.
254 : make_const_view(make_not_null(&spin_weighted_transform.data()),
255 : spin_weighted_data_view, 0,
256 : spin_weighted_data_view.size());
257 :
258 : // libsharp_to_goldberg_modes expects
259 : // SpinWeighted<ComplexModalVector>, but we don't know the spin until
260 : // we loop over tensors, so we have goldberg_mode_buffer (a
261 : // ComplexModalVector) allocated properly above and just point into it
262 : // here.
263 : goldberg_modes.set_data_ref(make_not_null(&goldberg_mode_buffer));
264 :
265 : // Transform nodal data to modal data
266 : Spectral::Swsh::libsharp_to_goldberg_modes(
267 : make_not_null(&goldberg_modes),
268 : Spectral::Swsh::swsh_transform(l_max, 1, spin_weighted_transform),
269 : l_max);
270 :
271 : // Copy data into std::vector for writing (remember 0th component is
272 : // time and was written above).
273 : for (size_t i = 0; i < l_max_plus_one_squared; ++i) {
274 : data_to_write[2 * i + 1] = real(goldberg_modes.data()[i]);
275 : data_to_write[2 * i + 2] = imag(goldberg_modes.data()[i]);
276 : }
277 : };
278 :
279 : // The inertial retarded time is special because it's stored as a
280 : // Scalar<DataVector> because it's only real and only has one set of angular
281 : // points worth of data to write. However, all the machinery above is for a
282 : // SpinWeighted<ComplexDataVector>. Luckily there is a
283 : // ComplexInertialRetardedTime where the real part is the
284 : // InertialRetardedTime and the imaginary part is 0, so we use that instead
285 : // swapping the names and legend where necessary
286 : const std::string inertial_retarded_time_name =
287 : detail::name<Tags::ComplexInertialRetardedTime>();
288 : if (variables_to_observe_.count(inertial_retarded_time_name) == 1) {
289 : const std::string subfile_name =
290 : "/Cce/VolumeData/" + inertial_retarded_time_name;
291 :
292 : // Legend
293 : std::vector<std::string> inertial_retarded_time_legend;
294 : inertial_retarded_time_legend.reserve(l_max_plus_one_squared + 1);
295 : inertial_retarded_time_legend.emplace_back("time");
296 : for (int i = 0; i <= static_cast<int>(l_max); ++i) {
297 : for (int j = -i; j <= i; ++j) {
298 : inertial_retarded_time_legend.push_back(MakeString{} << "Y_" << i
299 : << "," << j);
300 : }
301 : }
302 :
303 : // These have to be here because of the spin template
304 : const SpinWeighted<ComplexDataVector, 0> spin_weighted_transform{};
305 : SpinWeighted<ComplexModalVector, 0> goldberg_modes{};
306 :
307 : // Actually transform the time to complex modal data. Radial index 0
308 : // because this isn't volume data. It only holds one shell of data.
309 : transform_to_modal(Tags::ComplexInertialRetardedTime{},
310 : spin_weighted_transform, goldberg_modes, 0);
311 :
312 : // Buffer to write
313 : std::vector<double> inertial_retarded_time_to_write(
314 : l_max_plus_one_squared + 1);
315 :
316 : inertial_retarded_time_to_write[0] = time;
317 : // Only copy real data
318 : for (size_t i = 0; i < l_max_plus_one_squared; ++i) {
319 : inertial_retarded_time_to_write[i + 1] = data_to_write[2 * i + 1];
320 : }
321 :
322 : // Send to observer writer
323 : Parallel::threaded_action<
324 : observers::ThreadedActions::WriteReductionDataRow>(
325 : observer_proxy, subfile_name, inertial_retarded_time_legend,
326 : std::make_tuple(std::move(inertial_retarded_time_to_write)));
327 : }
328 :
329 : // One minus y is also special because every angular grid point for a given
330 : // radius holds the same value. Thus we only need to write one double per
331 : // radial grid point corresponding to 1 - y. The subfile name is just the
332 : // name of the tag, and the column names correspond to the names of the
333 : // radial subfiles for the spin weighted quantities
334 : const std::string one_minus_y_name = detail::name<Tags::OneMinusY>();
335 : if (variables_to_observe_.count(one_minus_y_name) == 1) {
336 : const std::string subfile_name = "/Cce/VolumeData/" + one_minus_y_name;
337 : std::vector<double> one_minus_y_to_write;
338 : std::vector<std::string> one_minus_y_legend;
339 : one_minus_y_to_write.reserve(number_of_radial_grid_points + 1);
340 : one_minus_y_legend.reserve(number_of_radial_grid_points + 1);
341 : one_minus_y_to_write.emplace_back(time);
342 : one_minus_y_legend.emplace_back("time");
343 :
344 : const ComplexDataVector& one_minus_y =
345 : get(db::get<Tags::OneMinusY>(box)).data();
346 :
347 : // All nodal data for each radius are the same value so we just take the
348 : // first one
349 : for (size_t radial_index = 0; radial_index < number_of_radial_grid_points;
350 : radial_index++) {
351 : one_minus_y_to_write.emplace_back(
352 : real(one_minus_y[radial_index * number_of_angular_points]));
353 : one_minus_y_legend.emplace_back("CompactifiedRadius_" +
354 : std::to_string(radial_index));
355 : }
356 :
357 : // Send to observer writer
358 : Parallel::threaded_action<
359 : observers::ThreadedActions::WriteReductionDataRow>(
360 : observer_proxy, subfile_name, one_minus_y_legend,
361 : std::make_tuple(std::move(one_minus_y_to_write)));
362 : }
363 :
364 : // Everything needs the same time so we just write it once here. We use the
365 : // code time because the inertial retarded time is specified over the whole
366 : // sphere and is written above (as a function of the code time as well)
367 : data_to_write[0] = time;
368 :
369 : // Loop over all available spin weighted tags and check if we are observing
370 : // this tag. We just capture everything in the scope because we need a
371 : // majority of the variables anyways
372 : tmpl::for_each<spin_weighted_tags_to_observe>([&](auto tag_v) {
373 : using tag = tmpl::type_from<decltype(tag_v)>;
374 : constexpr int spin = tag::type::type::spin;
375 : const std::string name = detail::name<tag>();
376 :
377 : // If we aren't observing this tag, then skip it
378 : if (variables_to_observe_.count(name) != 1) {
379 : return;
380 : }
381 :
382 : // These have to be here because of the spin template
383 : const SpinWeighted<ComplexDataVector, spin> spin_weighted_transform{};
384 : SpinWeighted<ComplexModalVector, spin> goldberg_modes{};
385 :
386 : // If we are observing this tag, loop over all radii and write data to
387 : // separate subfiles for each radius
388 : for (size_t radial_index = 0; radial_index < number_of_radial_grid_points;
389 : radial_index++) {
390 : const std::string subfile_name = "/Cce/VolumeData/" + name +
391 : "/CompactifiedRadius_" +
392 : std::to_string(radial_index);
393 :
394 : // Actually transform the time to complex modal data.
395 : transform_to_modal(tag{}, spin_weighted_transform, goldberg_modes,
396 : radial_index);
397 :
398 : // Send to observer writer
399 : Parallel::threaded_action<
400 : observers::ThreadedActions::WriteReductionDataRow>(
401 : observer_proxy, subfile_name, file_legend,
402 : std::make_tuple(data_to_write));
403 : }
404 : });
405 : }
406 :
407 0 : using is_ready_argument_tags = tmpl::list<>;
408 :
409 : template <typename Metavariables, typename ArrayIndex, typename Component>
410 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
411 : const ArrayIndex& /*array_index*/,
412 : const Component* const /*meta*/) const {
413 : return true;
414 : }
415 :
416 1 : bool needs_evolved_variables() const override { return true; }
417 :
418 : // NOLINTNEXTLINE(google-runtime-references)
419 0 : void pup(PUP::er& p) override {
420 : Event::pup(p);
421 : p | variables_to_observe_;
422 : }
423 :
424 : private:
425 0 : std::unordered_set<std::string> variables_to_observe_{};
426 : };
427 :
428 : ObserveFields::ObserveFields(
429 : const std::vector<std::string>& variables_to_observe,
430 : const Options::Context& context)
431 : : variables_to_observe_([&context, &variables_to_observe]() {
432 : std::unordered_set<std::string> result{};
433 : for (const auto& tensor : variables_to_observe) {
434 : if (result.count(tensor) != 0) {
435 : PARSE_ERROR(
436 : context,
437 : "Listed variable '"
438 : << tensor
439 : << "' more than once in list of variables to observe.");
440 : }
441 : result.insert(tensor);
442 : }
443 : return result;
444 : }()) {
445 : std::unordered_set<std::string> valid_tensors{};
446 : tmpl::for_each<available_tags_to_observe>([&valid_tensors](auto tag_v) {
447 : using tag = tmpl::type_from<decltype(tag_v)>;
448 : valid_tensors.insert(detail::name<tag>());
449 : });
450 :
451 : for (const auto& name : variables_to_observe_) {
452 : if (valid_tensors.count(name) != 1) {
453 : PARSE_ERROR(
454 : context,
455 : name << " is not an available variable. Available variables:\n"
456 : << valid_tensors);
457 : }
458 : }
459 : }
460 :
461 : /// \cond
462 : PUP::able::PUP_ID ObserveFields::my_PUP_ID = 0; // NOLINT
463 : /// \endcond
464 : } // namespace Events
465 : } // namespace Cce
|
