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 : * \f{align*}{
84 : * \Psi_0^{\prime (5)}
85 : * =& \Psi_0^{(5)} + 2 \eth u^\prime \Psi_1^{(4)}
86 : * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_2^{(3)}
87 : * + \frac{1}{2} \left( \eth u^\prime\right)^3 \Psi_3^{(2)}
88 : * + \frac{1}{16} \left(\eth u^\prime\right)^4 \Psi_4^{(1)}, \\
89 : * \Psi_1^{\prime (4)}
90 : * =& \Psi_1^{(4)} + \frac{3}{2} \eth u^\prime \Psi_2^{(3)}
91 : * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_3^{(2)}
92 : * + \frac{1}{8} \left(\eth u^\prime\right)^3 \Psi_4^{(1)}, \\
93 : * \Psi_2^{\prime (3)}
94 : * =& \Psi_2^{(3)}
95 : * + \eth u^\prime \Psi_3^{(2)}
96 : * + \frac{1}{4} \left(\eth u^\prime\right)^2 \Psi_4^{(1)}, \\
97 : * \Psi_3^{\prime (2)}
98 : * =& \Psi_3^{(2)} + \frac{1}{2} \eth u^{\prime} \Psi_4^{ (1)}, \\
99 : * \Psi_4^{\prime (1)}
100 : * =& \Psi_4^{(1)}.
101 : * \f}
102 : *
103 : * \note If \p WriteSynchronously is true, then a local synchronous action will
104 : * be used to write the News value rather than a threaded action.
105 : *
106 : * \ref DataBoxGroup changes:
107 : * - Adds: nothing
108 : * - Removes: nothing
109 : * - Modifies: `InterpolagionManager<ComplexDataVector, Tag>` for each `Tag` in
110 : * `Metavariables::scri_values_to_observe`
111 : */
112 : template <typename ObserverWriterComponent, typename BoundaryComponent,
113 : bool WriteSynchronously = true>
114 1 : struct ScriObserveInterpolated {
115 0 : using const_global_cache_tags = tmpl::flatten<
116 : tmpl::list<Tags::ObservationLMax,
117 : std::conditional_t<
118 : tt::is_a_v<AnalyticWorldtubeBoundary, BoundaryComponent>,
119 : tmpl::list<Tags::OutputNoninertialNews>, tmpl::list<>>>>;
120 : template <typename DbTags, typename... InboxTags, typename Metavariables,
121 : typename ArrayIndex, typename ActionList,
122 : typename ParallelComponent>
123 0 : static Parallel::iterable_action_return_t apply(
124 : db::DataBox<DbTags>& box,
125 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
126 : Parallel::GlobalCache<Metavariables>& cache,
127 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
128 : const ParallelComponent* const /*meta*/) {
129 : const size_t observation_l_max = db::get<Tags::ObservationLMax>(box);
130 : const size_t l_max = db::get<Tags::LMax>(box);
131 : const double extraction_radius =
132 : Parallel::get<Tags::ExtractionRadius>(cache);
133 : const std::string subfile_name =
134 : MakeString{} << "SpectreR" << std::setfill('0') << std::setw(4)
135 : << std::lround(extraction_radius);
136 : ComplexModalVector goldberg_modes{square(l_max + 1)};
137 :
138 : // alternative for the coordinate transformation getting scri+ values of the
139 : // weyl scalars:
140 : // need to obtain the eth of the inertial retarded time, each of the Weyl
141 : // scalars, and then we'll perform a transformation on that temporary
142 : // variables object, then output.
143 : // it won't be as general, but that's largely fine. The main frustration is
144 : // the loss of precision.
145 : Variables<detail::weyl_correction_list> corrected_scri_plus_weyl{
146 : Spectral::Swsh::number_of_swsh_collocation_points(l_max)};
147 :
148 : while (
149 : db::get<Tags::InterpolationManager<
150 : ComplexDataVector,
151 : tmpl::front<typename Metavariables::scri_values_to_observe>>>(box)
152 : .first_time_is_ready_to_interpolate()) {
153 : // first get the weyl scalars and correct them
154 : double interpolation_time = 0.0;
155 : tmpl::for_each<detail::weyl_correction_list>(
156 : [&interpolation_time, &corrected_scri_plus_weyl, &box](auto tag_v) {
157 : using tag = typename decltype(tag_v)::type;
158 : std::pair<double, ComplexDataVector> interpolation;
159 : db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
160 : [&interpolation](
161 : const gsl::not_null<
162 : ScriPlusInterpolationManager<ComplexDataVector, tag>*>
163 : interpolation_manager) {
164 : interpolation =
165 : interpolation_manager->interpolate_and_pop_first_time();
166 : },
167 : make_not_null(&box));
168 : interpolation_time = interpolation.first;
169 : get(get<tag>(corrected_scri_plus_weyl)).data() =
170 : interpolation.second;
171 : });
172 :
173 : detail::correct_weyl_scalars_for_inertial_time(
174 : make_not_null(&corrected_scri_plus_weyl));
175 :
176 : // add them to the data to write
177 : std::unordered_map<std::string, std::vector<double>> data_to_write{};
178 : tmpl::for_each<detail::weyl_correction_list>(
179 : [&data_to_write, &corrected_scri_plus_weyl, &interpolation_time,
180 : &observation_l_max, &l_max, &goldberg_modes](auto tag_v) {
181 : using tag = typename decltype(tag_v)::type;
182 : if constexpr (tmpl::list_contains_v<
183 : typename Metavariables::scri_values_to_observe,
184 : tag>) {
185 : data_to_write[detail::ScriOutput<tag>::name()] =
186 : transform_nodal_data(
187 : make_not_null(&goldberg_modes),
188 : get(get<tag>(corrected_scri_plus_weyl)).data(),
189 : interpolation_time, l_max, observation_l_max, tag{});
190 : }
191 : });
192 :
193 : // then do the interpolation and add the rest of the tags.
194 : tmpl::for_each<
195 : tmpl::list_difference<typename Metavariables::scri_values_to_observe,
196 : detail::weyl_correction_list>>(
197 : [&box, &data_to_write, &observation_l_max, &l_max,
198 : &goldberg_modes](auto tag_v) {
199 : using tag = typename decltype(tag_v)::type;
200 : std::pair<double, ComplexDataVector> interpolation;
201 : db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
202 : [&interpolation](
203 : const gsl::not_null<
204 : ScriPlusInterpolationManager<ComplexDataVector, tag>*>
205 : interpolation_manager) {
206 : interpolation =
207 : interpolation_manager->interpolate_and_pop_first_time();
208 : },
209 : make_not_null(&box));
210 : data_to_write[detail::ScriOutput<tag>::name()] =
211 : transform_nodal_data(make_not_null(&goldberg_modes),
212 : interpolation.second, interpolation.first,
213 : l_max, observation_l_max, tag{});
214 : });
215 :
216 : auto observer_proxy =
217 : Parallel::get_parallel_component<ObserverWriterComponent>(cache)[0];
218 : if constexpr (WriteSynchronously) {
219 : Parallel::local_synchronous_action<
220 : Cce::Actions::WriteScriBondiQuantities>(
221 : observer_proxy, cache, subfile_name, observation_l_max,
222 : std::move(data_to_write));
223 : } else {
224 : Parallel::threaded_action<Cce::Actions::WriteScriBondiQuantities>(
225 : observer_proxy, subfile_name, observation_l_max,
226 : std::move(data_to_write));
227 : }
228 :
229 : // output the expected news associated with the time value interpolated at
230 : // scri+.
231 : if constexpr (tmpl::list_contains_v<DbTags,
232 : Tags::AnalyticBoundaryDataManager>) {
233 : if (not db::get<Tags::OutputNoninertialNews>(box)) {
234 : db::get<Tags::AnalyticBoundaryDataManager>(box)
235 : .template write_news<ParallelComponent>(cache,
236 : interpolation_time);
237 : }
238 : }
239 : }
240 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
241 : }
242 :
243 : private:
244 : template <typename Tag>
245 0 : static std::vector<double> transform_nodal_data(
246 : const gsl::not_null<ComplexModalVector*> goldberg_mode_buffer,
247 : const ComplexDataVector& data, const double time, const size_t l_max,
248 : const size_t observation_l_max, const Tag& /*meta*/) {
249 : constexpr int Spin = Tag::type::type::spin;
250 : const SpinWeighted<ComplexDataVector, Spin> to_transform;
251 : make_const_view(make_not_null(&to_transform.data()), data, 0, data.size());
252 : SpinWeighted<ComplexModalVector, Spin> goldberg_modes;
253 : goldberg_modes.set_data_ref(goldberg_mode_buffer);
254 : Spectral::Swsh::libsharp_to_goldberg_modes(
255 : make_not_null(&goldberg_modes),
256 : Spectral::Swsh::swsh_transform(l_max, 1, to_transform), l_max);
257 :
258 : std::vector<double> modal_data(2 * square(observation_l_max + 1) + 1);
259 :
260 : modal_data[0] = time;
261 : for (size_t i = 0; i < square(observation_l_max + 1); ++i) {
262 : modal_data[2 * i + 1] = real(goldberg_modes.data()[i]);
263 : modal_data[2 * i + 2] = imag(goldberg_modes.data()[i]);
264 : }
265 :
266 : return modal_data;
267 : }
268 : };
269 : } // namespace Actions
270 : } // namespace Cce
|
