Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cmath>
7 : #include <cstddef>
8 : #include <iomanip>
9 : #include <iterator>
10 : #include <mutex>
11 : #include <string>
12 : #include <utility>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/TagName.hpp"
17 : #include "DataStructures/Tensor/TypeAliases.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "DataStructures/VariablesTag.hpp"
20 : #include "Evolution/Systems/Cce/BoundaryData.hpp"
21 : #include "Evolution/Systems/Cce/OptionTags.hpp"
22 : #include "Evolution/Systems/Cce/Tags.hpp"
23 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
24 : #include "IO/Observer/Actions/GetLockPointer.hpp"
25 : #include "IO/Observer/ObserverComponent.hpp"
26 : #include "IO/Observer/ReductionActions.hpp"
27 : #include "NumericalAlgorithms/Spectral/SwshCoefficients.hpp"
28 : #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
29 : #include "Parallel/GlobalCache.hpp"
30 : #include "Parallel/Invoke.hpp"
31 : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
32 : #include "ParallelAlgorithms/Interpolation/Protocols/PostInterpolationCallback.hpp"
33 : #include "ParallelAlgorithms/Interpolation/Targets/AngularOrdering.hpp"
34 : #include "ParallelAlgorithms/Interpolation/Targets/Sphere.hpp"
35 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
36 : #include "Utilities/ConstantExpressions.hpp"
37 : #include "Utilities/ErrorHandling/Assert.hpp"
38 : #include "Utilities/ErrorHandling/Error.hpp"
39 : #include "Utilities/Gsl.hpp"
40 : #include "Utilities/MakeString.hpp"
41 : #include "Utilities/PrettyType.hpp"
42 : #include "Utilities/ProtocolHelpers.hpp"
43 : #include "Utilities/TMPL.hpp"
44 :
45 1 : namespace intrp {
46 1 : namespace callbacks {
47 : /*!
48 : * \brief Post interpolation callback that dumps metric data in Bondi-Sachs form
49 : * on a number of extraction radii given by the `intrp::TargetPoints::Sphere`
50 : * target.
51 : *
52 : * To use this callback, the target must be the `intrp::TargetPoints::Sphere`
53 : * target in the inertial frame. This callback also expects that the GH source
54 : * vars on each of the target spheres are:
55 : *
56 : * - `gr::Tags::SpacetimeMetric`
57 : * - `gh::Tags::Pi`
58 : * - `gh::Tags::Phi`
59 : *
60 : * This callback will write a new `H5` file for each extraction radius in the
61 : * Sphere target. The name of this file will be a file prefix specified by the
62 : * Cce::Tags::FilePrefix prepended onto `CceRXXXX.h5` where the XXXX is the
63 : * zero-padded extraction radius rounded to the nearest integer. The quantities
64 : * that will be written are
65 : *
66 : * - `Cce::Tags::BondiBeta`
67 : * - `Cce::Tags::Dr<Cce::Tags::BondiJ>`
68 : * - `Cce::Tags::Du<Cce::Tags::BondiR>`
69 : * - `Cce::Tags::BondiH`
70 : * - `Cce::Tags::BondiJ`
71 : * - `Cce::Tags::BondiQ`
72 : * - `Cce::Tags::BondiR`
73 : * - `Cce::Tags::BondiU`
74 : * - `Cce::Tags::BondiW`
75 : *
76 : * \note For all real quantities (Beta, DuR, R, W) we omit writing the
77 : * negative m modes, and the imaginary part of the m = 0 mode.
78 : */
79 : template <typename InterpolationTargetTag>
80 1 : struct DumpBondiSachsOnWorldtube
81 : : tt::ConformsTo<intrp::protocols::PostInterpolationCallback> {
82 0 : static constexpr double fill_invalid_points_with =
83 : std::numeric_limits<double>::quiet_NaN();
84 :
85 0 : using const_global_cache_tags = tmpl::list<Cce::Tags::FilePrefix>;
86 :
87 0 : using cce_boundary_tags = Cce::Tags::characteristic_worldtube_boundary_tags<
88 : Cce::Tags::BoundaryValue>;
89 :
90 0 : using cce_tags_to_dump = db::wrap_tags_in<
91 : Cce::Tags::BoundaryValue,
92 : tmpl::list<Cce::Tags::BondiBeta, Cce::Tags::Dr<Cce::Tags::BondiJ>,
93 : Cce::Tags::Du<Cce::Tags::BondiR>, Cce::Tags::BondiH,
94 : Cce::Tags::BondiJ, Cce::Tags::BondiQ, Cce::Tags::BondiR,
95 : Cce::Tags::BondiU, Cce::Tags::BondiW>>;
96 :
97 0 : using gh_source_vars_for_cce =
98 : tmpl::list<gr::Tags::SpacetimeMetric<DataVector, 3>,
99 : gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>>;
100 :
101 0 : using gh_source_vars_from_interpolation =
102 : typename InterpolationTargetTag::vars_to_interpolate_to_target;
103 :
104 : static_assert(
105 : std::is_same_v<tmpl::list_difference<cce_tags_to_dump, cce_boundary_tags>,
106 : tmpl::list<>>,
107 : "Cce tags to dump are not in the boundary tags.");
108 :
109 : static_assert(
110 : tmpl::and_<
111 : std::is_same<tmpl::list_difference<gh_source_vars_from_interpolation,
112 : gh_source_vars_for_cce>,
113 : tmpl::list<>>,
114 : std::is_same<tmpl::list_difference<gh_source_vars_for_cce,
115 : gh_source_vars_from_interpolation>,
116 : tmpl::list<>>>::type::value,
117 : "To use DumpBondiSachsOnWorldtube, the GH source variables must be the "
118 : "spacetime metric, pi, and phi.");
119 :
120 : static_assert(
121 : std::is_same_v<typename InterpolationTargetTag::compute_target_points,
122 : intrp::TargetPoints::Sphere<InterpolationTargetTag,
123 : ::Frame::Inertial>>,
124 : "To use the DumpBondiSachsOnWorltube post interpolation callback, you "
125 : "must use the intrp::TargetPoints::Sphere target in the inertial "
126 : "frame");
127 :
128 : template <typename DbTags, typename Metavariables, typename TemporalId>
129 0 : static void apply(const db::DataBox<DbTags>& box,
130 : Parallel::GlobalCache<Metavariables>& cache,
131 : const TemporalId& temporal_id) {
132 : const auto& sphere =
133 : Parallel::get<Tags::Sphere<InterpolationTargetTag>>(cache);
134 : const auto& filename_prefix = Parallel::get<Cce::Tags::FilePrefix>(cache);
135 :
136 : if (sphere.angular_ordering != intrp::AngularOrdering::Cce) {
137 : ERROR(
138 : "To use the DumpBondiSachsOnWorldtube post interpolation callback, "
139 : "the angular ordering of the Spheres must be Cce, not "
140 : << sphere.angular_ordering);
141 : }
142 :
143 : const auto& radii = sphere.radii;
144 : const size_t l_max = sphere.l_max;
145 : const size_t num_points_single_sphere =
146 : Spectral::Swsh::number_of_swsh_collocation_points(l_max);
147 :
148 : const auto& all_gh_vars =
149 : db::get<::Tags::Variables<gh_source_vars_from_interpolation>>(box);
150 :
151 : const auto& all_spacetime_metric =
152 : get<gr::Tags::SpacetimeMetric<DataVector, 3>>(all_gh_vars);
153 : const auto& all_pi = get<gh::Tags::Pi<DataVector, 3>>(all_gh_vars);
154 : const auto& all_phi = get<gh::Tags::Phi<DataVector, 3>>(all_gh_vars);
155 :
156 : const tnsr::aa<DataVector, 3, ::Frame::Inertial> spacetime_metric;
157 : const tnsr::aa<DataVector, 3, ::Frame::Inertial> pi;
158 : const tnsr::iaa<DataVector, 3, ::Frame::Inertial> phi;
159 :
160 : // Bondi data
161 : Variables<cce_boundary_tags> bondi_boundary_data{num_points_single_sphere};
162 : ComplexModalVector goldberg_mode_buffer{square(l_max + 1)};
163 : const std::vector<std::string> all_legend = build_legend(l_max, false);
164 : const std::vector<std::string> real_legend = build_legend(l_max, true);
165 :
166 : size_t offset = 0;
167 : for (const auto& radius : radii) {
168 : // Set data references so we don't copy data unnecessarily
169 : for (size_t a = 0; a < 4; a++) {
170 : for (size_t b = 0; b < 4; b++) {
171 : make_const_view(make_not_null(&spacetime_metric.get(a, b)),
172 : all_spacetime_metric.get(a, b), offset,
173 : num_points_single_sphere);
174 : make_const_view(make_not_null(&pi.get(a, b)), all_pi.get(a, b),
175 : offset, num_points_single_sphere);
176 : for (size_t i = 0; i < 3; i++) {
177 : make_const_view(make_not_null(&phi.get(i, a, b)),
178 : all_phi.get(i, a, b), offset,
179 : num_points_single_sphere);
180 : }
181 : }
182 : }
183 :
184 : offset += num_points_single_sphere;
185 :
186 : Cce::create_bondi_boundary_data(make_not_null(&bondi_boundary_data), phi,
187 : pi, spacetime_metric, radius, l_max);
188 :
189 : const std::string filename =
190 : MakeString{} << filename_prefix << "CceR" << std::setfill('0')
191 : << std::setw(4) << std::lround(radius);
192 :
193 : tmpl::for_each<cce_tags_to_dump>(
194 : [&temporal_id, &l_max, &all_legend, &real_legend, &filename,
195 : &bondi_boundary_data, &goldberg_mode_buffer, &cache](auto tag_v) {
196 : using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
197 : constexpr int spin = tag::tag::type::type::spin;
198 : // Spin = 0 does not imply a quantity is real. However, all the tags
199 : // we want to print out that have spin = 0 happen to be real, so we
200 : // use this as an indicator.
201 : constexpr bool is_real = spin == 0;
202 :
203 : const auto& legend = is_real ? real_legend : all_legend;
204 :
205 : // `tag` is a BoundaryValue. We want the actual tag name
206 : const std::string subfile_name{
207 : "/" + replace_name(db::tag_name<typename tag::tag>())};
208 :
209 : const auto& bondi_data = get(get<tag>(bondi_boundary_data));
210 :
211 : // Convert our modal data to goldberg modes
212 : SpinWeighted<ComplexModalVector, spin> goldberg_modes;
213 : goldberg_modes.set_data_ref(make_not_null(&goldberg_mode_buffer));
214 : Spectral::Swsh::libsharp_to_goldberg_modes(
215 : make_not_null(&goldberg_modes),
216 : Spectral::Swsh::swsh_transform(l_max, 1, bondi_data), l_max);
217 :
218 : std::vector<double> data_to_write_buffer;
219 : data_to_write_buffer.reserve(number_of_components(l_max, is_real));
220 : data_to_write_buffer.emplace_back(
221 : intrp::InterpolationTarget_detail::get_temporal_id_value(
222 : temporal_id));
223 :
224 : // We loop over ell and m rather than just the total number of modes
225 : // because we don't print negative m or the imaginary part of m=0
226 : // for real quantities.
227 : for (size_t ell = 0; ell <= l_max; ell++) {
228 : for (int m = is_real ? 0 : -static_cast<int>(ell);
229 : m <= static_cast<int>(ell); m++) {
230 : const size_t goldberg_index =
231 : Spectral::Swsh::goldberg_mode_index(l_max, ell, m);
232 : data_to_write_buffer.push_back(
233 : real(goldberg_modes.data()[goldberg_index]));
234 : if (not is_real or m != 0) {
235 : data_to_write_buffer.push_back(
236 : imag(goldberg_modes.data()[goldberg_index]));
237 : }
238 : }
239 : }
240 :
241 : ASSERT(legend.size() == data_to_write_buffer.size(),
242 : "Legend (" << legend.size()
243 : << ") does not have the same number of "
244 : "components as data to write ("
245 : << data_to_write_buffer.size() << ") for tag "
246 : << db::tag_name<typename tag::tag>());
247 :
248 : // Even though no other cores should be writing to this file, we
249 : // still need to get the h5 file lock so the system hdf5 doesn't get
250 : // upset
251 : auto* hdf5_lock =
252 : Parallel::local_branch(
253 : Parallel::get_parallel_component<
254 : observers::ObserverWriter<Metavariables>>(cache))
255 : ->template local_synchronous_action<
256 : observers::Actions::GetLockPointer<
257 : observers::Tags::H5FileLock>>();
258 :
259 : {
260 : const std::lock_guard lock(*hdf5_lock);
261 : observers::ThreadedActions::ReductionActions_detail::write_data(
262 : subfile_name, observers::input_source_from_cache(cache),
263 : legend, std::make_tuple(data_to_write_buffer), filename,
264 : std::index_sequence<0>{});
265 : }
266 : });
267 : }
268 : }
269 :
270 : private:
271 : // There are exactly half the number of modes for spin = 0 quantities as their
272 : // are for spin != 0
273 0 : static size_t number_of_components(const size_t l_max, const bool is_real) {
274 : return 1 + square(l_max + 1) * (is_real ? 1 : 2);
275 : }
276 :
277 0 : static std::vector<std::string> build_legend(const size_t l_max,
278 : const bool is_real) {
279 : std::vector<std::string> legend;
280 : legend.reserve(number_of_components(l_max, is_real));
281 : legend.emplace_back("Time");
282 : for (int ell = 0; ell <= static_cast<int>(l_max); ++ell) {
283 : for (int m = is_real ? 0 : -ell; m <= ell; ++m) {
284 : legend.push_back(MakeString{} << "Re(" << ell << "," << m << ")");
285 : // For real quantities, don't include the imaginary m=0
286 : if (not is_real or m != 0) {
287 : legend.push_back(MakeString{} << "Im(" << ell << "," << m << ")");
288 : }
289 : }
290 : }
291 : return legend;
292 : }
293 :
294 : // These match names that CCE executable expects
295 0 : static std::string replace_name(const std::string& db_tag_name) {
296 : if (db_tag_name == "BondiBeta") {
297 : return "Beta";
298 : } else if (db_tag_name == "Dr(J)") {
299 : return "DrJ";
300 : } else if (db_tag_name == "Du(R)") {
301 : return "DuR";
302 : } else {
303 : return db_tag_name;
304 : }
305 : }
306 : };
307 : } // namespace callbacks
308 : } // namespace intrp
|
