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 <optional>
8 : #include <string>
9 : #include <tuple>
10 : #include <unordered_map>
11 : #include <utility>
12 : #include <vector>
13 :
14 : #include "DataStructures/ComplexDataVector.hpp"
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/SpinWeighted.hpp"
17 : #include "Evolution/Systems/Cce/Actions/WriteScriBondiQuantities.hpp"
18 : #include "Evolution/Systems/Cce/OptionTags.hpp"
19 : #include "Evolution/Systems/Cce/ScriPlusInterpolationManager.hpp"
20 : #include "Evolution/Systems/Cce/Tags.hpp"
21 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCoefficients.hpp"
22 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshTransform.hpp"
23 : #include "Parallel/AlgorithmExecution.hpp"
24 : #include "Parallel/GlobalCache.hpp"
25 : #include "Parallel/Invoke.hpp"
26 : #include "Parallel/Local.hpp"
27 : #include "Utilities/Gsl.hpp"
28 : #include "Utilities/MakeString.hpp"
29 : #include "Utilities/System/ParallelInfo.hpp"
30 : #include "Utilities/TMPL.hpp"
31 :
32 : namespace Cce {
33 : /// \cond
34 : template <typename Metavariables>
35 : struct AnalyticWorldtubeBoundary;
36 : /// \endcond
37 : namespace Actions {
38 : namespace detail {
39 : // Provide a nicer name for the output h5 files for some of the uglier
40 : // combinations we need
41 : template <typename Tag>
42 : struct ScriOutput {
43 : static std::string name() { return db::tag_name<Tag>(); }
44 : };
45 : template <typename Tag>
46 : struct ScriOutput<Tags::ScriPlus<Tag>> {
47 : static std::string name() { return pretty_type::short_name<Tag>(); }
48 : };
49 : template <>
50 : struct ScriOutput<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>> {
51 : static std::string name() { return "Psi4"; }
52 : };
53 :
54 : using weyl_correction_list =
55 : tmpl::list<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>,
56 : Tags::ScriPlus<Tags::Psi3>, Tags::ScriPlus<Tags::Psi2>,
57 : Tags::ScriPlus<Tags::Psi1>, Tags::ScriPlus<Tags::Psi0>,
58 : Tags::EthInertialRetardedTime>;
59 :
60 : void correct_weyl_scalars_for_inertial_time(
61 : gsl::not_null<Variables<weyl_correction_list>*> weyl_correction_variables);
62 : } // namespace detail
63 :
64 : /*!
65 : * \ingroup ActionsGroup
66 : * \brief Checks the interpolation managers and if they are ready, performs the
67 : * interpolation and sends the data to file.
68 : *
69 : * \details This uses the `ScriPlusInterpolationManager` to perform the
70 : * interpolations of all requested scri quantities (determined by
71 : * `scri_values_to_observe` in the metavariables), and write them to disk using
72 : * `observers::threadedActions::WriteSimpleData`. When using an analytic
73 : * worldtube solution, this action also uses the `AnalyticBoundaryDataManager`
74 : * to output the expected News value at the appropriate asymptotically inertial
75 : * time.
76 : *
77 : * \note This action also uses the `Tags::EthInertialRetardedTime`, interpolated
78 : * to the inertial frame, to perform the coordinate transformations presented in
79 : * \cite Boyle:2015nqa to the Weyl scalars after interpolation. For our
80 : * formulas, we need to adjust the signs and factors of two to be compatible
81 : * with our definitions of \f$\eth\f$ and choice of Newman-Penrose tetrad.
82 : *
83 : * The coefficient of \f$\left(\eth u^\prime\right)^2\Psi_2^{(3)}\f$ in
84 : * \f$\Psi_0^{\prime(5)}\f$ is incorrectly given as \f$3/4\f$ in
85 : * \cite Moxon2020gha, Eq. (94b), and in \cite Moxon2021gbv, Eq. (B2g) of the
86 : * arXiv version. The value \f$3/2\f$ below agrees with the BMS transformation
87 : * in \cite Boyle:2015nqa, Eqs. (17a--e).
88 : *
89 : * \f{align*}{
90 : * \Psi_0^{\prime (5)}
91 : * =& \Psi_0^{(5)} + 2 \eth u^\prime \Psi_1^{(4)}
92 : * + \frac{3}{2} \left(\eth u^\prime\right)^2 \Psi_2^{(3)}
93 : * + \frac{1}{2} \left( \eth u^\prime\right)^3 \Psi_3^{(2)}
94 : * + \frac{1}{16} \left(\eth u^\prime\right)^4 \Psi_4^{(1)}, \\
95 : * \Psi_1^{\prime (4)}
96 : * =& \Psi_1^{(4)} + \frac{3}{2} \eth u^\prime \Psi_2^{(3)}
97 : * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_3^{(2)}
98 : * + \frac{1}{8} \left(\eth u^\prime\right)^3 \Psi_4^{(1)}, \\
99 : * \Psi_2^{\prime (3)}
100 : * =& \Psi_2^{(3)}
101 : * + \eth u^\prime \Psi_3^{(2)}
102 : * + \frac{1}{4} \left(\eth u^\prime\right)^2 \Psi_4^{(1)}, \\
103 : * \Psi_3^{\prime (2)}
104 : * =& \Psi_3^{(2)} + \frac{1}{2} \eth u^{\prime} \Psi_4^{ (1)}, \\
105 : * \Psi_4^{\prime (1)}
106 : * =& \Psi_4^{(1)}.
107 : * \f}
108 : *
109 : * \note If \p WriteSynchronously is true, then a local synchronous action will
110 : * be used to write the News value rather than a threaded action.
111 : *
112 : * \ref DataBoxGroup changes:
113 : * - Adds: nothing
114 : * - Removes: nothing
115 : * - Modifies: `InterpolagionManager<ComplexDataVector, Tag>` for each `Tag` in
116 : * `Metavariables::scri_values_to_observe`
117 : */
118 : template <typename ObserverWriterComponent, typename BoundaryComponent,
119 : bool WriteSynchronously = true>
120 1 : struct ScriObserveInterpolated {
121 0 : using const_global_cache_tags = tmpl::flatten<
122 : tmpl::list<Tags::ObservationLMax,
123 : std::conditional_t<
124 : tt::is_a_v<AnalyticWorldtubeBoundary, BoundaryComponent>,
125 : tmpl::list<Tags::OutputNoninertialNews>, tmpl::list<>>>>;
126 : template <typename DbTags, typename... InboxTags, typename Metavariables,
127 : typename ArrayIndex, typename ActionList,
128 : typename ParallelComponent>
129 0 : static Parallel::iterable_action_return_t apply(
130 : db::DataBox<DbTags>& box,
131 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
132 : Parallel::GlobalCache<Metavariables>& cache,
133 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
134 : const ParallelComponent* const /*meta*/) {
135 : const size_t observation_l_max = db::get<Tags::ObservationLMax>(box);
136 : const size_t l_max = db::get<Tags::LMax>(box);
137 : const double extraction_radius =
138 : Parallel::get<Tags::ExtractionRadius>(cache);
139 : const std::string subfile_name =
140 : MakeString{} << "SpectreR" << std::setfill('0') << std::setw(4)
141 : << std::lround(extraction_radius);
142 : ComplexModalVector goldberg_modes{square(l_max + 1)};
143 :
144 : // alternative for the coordinate transformation getting scri+ values of the
145 : // weyl scalars:
146 : // need to obtain the eth of the inertial retarded time, each of the Weyl
147 : // scalars, and then we'll perform a transformation on that temporary
148 : // variables object, then output.
149 : // it won't be as general, but that's largely fine. The main frustration is
150 : // the loss of precision.
151 : Variables<detail::weyl_correction_list> corrected_scri_plus_weyl{
152 : Spectral::Swsh::number_of_swsh_collocation_points(l_max)};
153 :
154 : while (
155 : db::get<Tags::InterpolationManager<
156 : ComplexDataVector,
157 : tmpl::front<typename Metavariables::scri_values_to_observe>>>(box)
158 : .first_time_is_ready_to_interpolate()) {
159 : // first get the weyl scalars and correct them
160 : double interpolation_time = 0.0;
161 : tmpl::for_each<detail::weyl_correction_list>(
162 : [&interpolation_time, &corrected_scri_plus_weyl, &box](auto tag_v) {
163 : using tag = typename decltype(tag_v)::type;
164 : std::pair<double, ComplexDataVector> interpolation;
165 : db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
166 : [&interpolation](
167 : const gsl::not_null<
168 : ScriPlusInterpolationManager<ComplexDataVector, tag>*>
169 : interpolation_manager) {
170 : interpolation =
171 : interpolation_manager->interpolate_and_pop_first_time();
172 : },
173 : make_not_null(&box));
174 : interpolation_time = interpolation.first;
175 : get(get<tag>(corrected_scri_plus_weyl)).data() =
176 : interpolation.second;
177 : });
178 :
179 : detail::correct_weyl_scalars_for_inertial_time(
180 : make_not_null(&corrected_scri_plus_weyl));
181 :
182 : // add them to the data to write
183 : std::unordered_map<std::string, std::vector<double>> data_to_write{};
184 : tmpl::for_each<detail::weyl_correction_list>(
185 : [&data_to_write, &corrected_scri_plus_weyl, &interpolation_time,
186 : &observation_l_max, &l_max, &goldberg_modes](auto tag_v) {
187 : using tag = typename decltype(tag_v)::type;
188 : if constexpr (tmpl::list_contains_v<
189 : typename Metavariables::scri_values_to_observe,
190 : tag>) {
191 : data_to_write[detail::ScriOutput<tag>::name()] =
192 : transform_nodal_data(
193 : make_not_null(&goldberg_modes),
194 : get(get<tag>(corrected_scri_plus_weyl)).data(),
195 : interpolation_time, l_max, observation_l_max, tag{});
196 : }
197 : });
198 :
199 : // then do the interpolation and add the rest of the tags.
200 : tmpl::for_each<
201 : tmpl::list_difference<typename Metavariables::scri_values_to_observe,
202 : detail::weyl_correction_list>>(
203 : [&box, &data_to_write, &observation_l_max, &l_max,
204 : &goldberg_modes](auto tag_v) {
205 : using tag = typename decltype(tag_v)::type;
206 : std::pair<double, ComplexDataVector> interpolation;
207 : db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
208 : [&interpolation](
209 : const gsl::not_null<
210 : ScriPlusInterpolationManager<ComplexDataVector, tag>*>
211 : interpolation_manager) {
212 : interpolation =
213 : interpolation_manager->interpolate_and_pop_first_time();
214 : },
215 : make_not_null(&box));
216 : data_to_write[detail::ScriOutput<tag>::name()] =
217 : transform_nodal_data(make_not_null(&goldberg_modes),
218 : interpolation.second, interpolation.first,
219 : l_max, observation_l_max, tag{});
220 : });
221 :
222 : auto observer_proxy =
223 : Parallel::get_parallel_component<ObserverWriterComponent>(cache)[0];
224 : if constexpr (WriteSynchronously) {
225 : Parallel::local_synchronous_action<
226 : Cce::Actions::WriteScriBondiQuantities>(
227 : observer_proxy, cache, subfile_name, observation_l_max,
228 : std::move(data_to_write));
229 : } else {
230 : Parallel::threaded_action<Cce::Actions::WriteScriBondiQuantities>(
231 : observer_proxy, subfile_name, observation_l_max,
232 : std::move(data_to_write));
233 : }
234 :
235 : // output the expected news associated with the time value interpolated at
236 : // scri+.
237 : if constexpr (tmpl::list_contains_v<DbTags,
238 : Tags::AnalyticBoundaryDataManager>) {
239 : if (not db::get<Tags::OutputNoninertialNews>(box)) {
240 : db::get<Tags::AnalyticBoundaryDataManager>(box)
241 : .template write_news<ParallelComponent>(cache,
242 : interpolation_time);
243 : }
244 : }
245 : }
246 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
247 : }
248 :
249 : private:
250 : template <typename Tag>
251 0 : static std::vector<double> transform_nodal_data(
252 : const gsl::not_null<ComplexModalVector*> goldberg_mode_buffer,
253 : const ComplexDataVector& data, const double time, const size_t l_max,
254 : const size_t observation_l_max, const Tag& /*meta*/) {
255 : constexpr int Spin = Tag::type::type::spin;
256 : const SpinWeighted<ComplexDataVector, Spin> to_transform;
257 : make_const_view(make_not_null(&to_transform.data()), data, 0, data.size());
258 : SpinWeighted<ComplexModalVector, Spin> goldberg_modes;
259 : goldberg_modes.set_data_ref(goldberg_mode_buffer);
260 : Spectral::Swsh::libsharp_to_goldberg_modes(
261 : make_not_null(&goldberg_modes),
262 : Spectral::Swsh::swsh_transform(l_max, 1, to_transform), l_max);
263 :
264 : std::vector<double> modal_data(2 * square(observation_l_max + 1) + 1);
265 :
266 : modal_data[0] = time;
267 : for (size_t i = 0; i < square(observation_l_max + 1); ++i) {
268 : modal_data[2 * i + 1] = real(goldberg_modes.data()[i]);
269 : modal_data[2 * i + 2] = imag(goldberg_modes.data()[i]);
270 : }
271 :
272 : return modal_data;
273 : }
274 : };
275 : } // namespace Actions
276 : } // namespace Cce
|
