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 <utility>
11 : #include <vector>
12 :
13 : #include "DataStructures/ComplexDataVector.hpp"
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/SpinWeighted.hpp"
16 : #include "Evolution/Systems/Cce/OptionTags.hpp"
17 : #include "Evolution/Systems/Cce/ScriPlusInterpolationManager.hpp"
18 : #include "Evolution/Systems/Cce/Tags.hpp"
19 : #include "IO/Observer/ReductionActions.hpp"
20 : #include "NumericalAlgorithms/Spectral/SwshCoefficients.hpp"
21 : #include "NumericalAlgorithms/Spectral/SwshTransform.hpp"
22 : #include "Parallel/AlgorithmExecution.hpp"
23 : #include "Parallel/GlobalCache.hpp"
24 : #include "Parallel/Invoke.hpp"
25 : #include "Parallel/Local.hpp"
26 : #include "Utilities/Gsl.hpp"
27 : #include "Utilities/MakeString.hpp"
28 : #include "Utilities/System/ParallelInfo.hpp"
29 : #include "Utilities/TMPL.hpp"
30 :
31 : namespace Cce {
32 : /// \cond
33 : template <typename Metavariables>
34 : struct AnalyticWorldtubeBoundary;
35 : /// \endcond
36 : namespace Actions {
37 : namespace detail {
38 : // Provide a nicer name for the output h5 files for some of the uglier
39 : // combinations we need
40 : template <typename Tag>
41 : struct ScriOutput {
42 : static std::string name() { return db::tag_name<Tag>(); }
43 : };
44 : template <typename Tag>
45 : struct ScriOutput<Tags::ScriPlus<Tag>> {
46 : static std::string name() { return pretty_type::short_name<Tag>(); }
47 : };
48 : template <>
49 : struct ScriOutput<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>> {
50 : static std::string name() { return "Psi4"; }
51 : };
52 :
53 : using weyl_correction_list =
54 : tmpl::list<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>,
55 : Tags::ScriPlus<Tags::Psi3>, Tags::ScriPlus<Tags::Psi2>,
56 : Tags::ScriPlus<Tags::Psi1>, Tags::ScriPlus<Tags::Psi0>,
57 : Tags::EthInertialRetardedTime>;
58 :
59 : void correct_weyl_scalars_for_inertial_time(
60 : gsl::not_null<Variables<weyl_correction_list>*> weyl_correction_variables);
61 : } // namespace detail
62 :
63 : /*!
64 : * \ingroup ActionsGroup
65 : * \brief Checks the interpolation managers and if they are ready, performs the
66 : * interpolation and sends the data to file via
67 : * `observers::ThreadedActions::WriteSimpleData`.
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 : * \ref DataBoxGroup changes:
104 : * - Adds: nothing
105 : * - Removes: nothing
106 : * - Modifies: `InterpolagionManager<ComplexDataVector, Tag>` for each `Tag` in
107 : * `Metavariables::scri_values_to_observe`
108 : */
109 : template <typename ObserverWriterComponent, typename BoundaryComponent>
110 1 : struct ScriObserveInterpolated {
111 0 : using const_global_cache_tags = tmpl::flatten<
112 : tmpl::list<Tags::ObservationLMax,
113 : std::conditional_t<
114 : tt::is_a_v<AnalyticWorldtubeBoundary, BoundaryComponent>,
115 : tmpl::list<Tags::OutputNoninertialNews>, tmpl::list<>>>>;
116 : template <typename DbTags, typename... InboxTags, typename Metavariables,
117 : typename ArrayIndex, typename ActionList,
118 : typename ParallelComponent>
119 0 : static Parallel::iterable_action_return_t apply(
120 : db::DataBox<DbTags>& box,
121 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
122 : Parallel::GlobalCache<Metavariables>& cache,
123 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
124 : const ParallelComponent* const /*meta*/) {
125 : const size_t observation_l_max = db::get<Tags::ObservationLMax>(box);
126 : const size_t l_max = db::get<Tags::LMax>(box);
127 : std::vector<double> data_to_write(2 * square(observation_l_max + 1) + 1);
128 : ComplexModalVector goldberg_modes{square(l_max + 1)};
129 : std::vector<std::string> file_legend;
130 : file_legend.reserve(2 * square(observation_l_max + 1) + 1);
131 : file_legend.emplace_back("time");
132 : for (int i = 0; i <= static_cast<int>(observation_l_max); ++i) {
133 : for (int j = -i; j <= i; ++j) {
134 : file_legend.push_back(MakeString{} << "Real Y_" << i << "," << j);
135 : file_legend.push_back(MakeString{} << "Imag Y_" << i << "," << j);
136 : }
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 : // then output each of them
177 : tmpl::for_each<detail::weyl_correction_list>(
178 : [&data_to_write, &corrected_scri_plus_weyl, &interpolation_time,
179 : &file_legend, &observation_l_max, &l_max, &cache,
180 : &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 : ScriObserveInterpolated::transform_and_write<
186 : tag, tag::type::type::spin, ParallelComponent>(
187 : get(get<tag>(corrected_scri_plus_weyl)).data(),
188 : interpolation_time, make_not_null(&goldberg_modes),
189 : make_not_null(&data_to_write), file_legend, l_max,
190 : observation_l_max, cache);
191 : }
192 : });
193 :
194 : // then do the interpolation and output of each of the rest of the tags.
195 : tmpl::for_each<
196 : tmpl::list_difference<typename Metavariables::scri_values_to_observe,
197 : detail::weyl_correction_list>>(
198 : [&box, &data_to_write, &file_legend, &observation_l_max, &l_max,
199 : &cache, &goldberg_modes](auto tag_v) {
200 : using tag = typename decltype(tag_v)::type;
201 : std::pair<double, ComplexDataVector> interpolation;
202 : db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
203 : [&interpolation](
204 : const gsl::not_null<
205 : ScriPlusInterpolationManager<ComplexDataVector, tag>*>
206 : interpolation_manager) {
207 : interpolation =
208 : interpolation_manager->interpolate_and_pop_first_time();
209 : },
210 : make_not_null(&box));
211 : ScriObserveInterpolated::transform_and_write<
212 : tag, tag::type::type::spin, ParallelComponent>(
213 : interpolation.second, interpolation.first,
214 : make_not_null(&goldberg_modes), make_not_null(&data_to_write),
215 : file_legend, l_max, observation_l_max, cache);
216 : });
217 : // output the expected news associated with the time value interpolated at
218 : // scri+.
219 : if constexpr (tmpl::list_contains_v<DbTags,
220 : Tags::AnalyticBoundaryDataManager>) {
221 : if (not db::get<Tags::OutputNoninertialNews>(box)) {
222 : db::get<Tags::AnalyticBoundaryDataManager>(box)
223 : .template write_news<ParallelComponent>(cache,
224 : interpolation_time);
225 : }
226 : }
227 : }
228 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
229 : }
230 :
231 : private:
232 : template <typename Tag, int Spin, typename ParallelComponent,
233 : typename Metavariables>
234 0 : static void transform_and_write(
235 : const ComplexDataVector& data, const double time,
236 : const gsl::not_null<ComplexModalVector*> goldberg_mode_buffer,
237 : const gsl::not_null<std::vector<double>*> data_to_write_buffer,
238 : const std::vector<std::string>& legend, const size_t l_max,
239 : const size_t observation_l_max,
240 : Parallel::GlobalCache<Metavariables>& cache) {
241 : const SpinWeighted<ComplexDataVector, Spin> to_transform;
242 : make_const_view(make_not_null(&to_transform.data()), data, 0, data.size());
243 : SpinWeighted<ComplexModalVector, Spin> goldberg_modes;
244 : goldberg_modes.set_data_ref(goldberg_mode_buffer);
245 : Spectral::Swsh::libsharp_to_goldberg_modes(
246 : make_not_null(&goldberg_modes),
247 : Spectral::Swsh::swsh_transform(l_max, 1, to_transform), l_max);
248 :
249 : (*data_to_write_buffer)[0] = time;
250 : for (size_t i = 0; i < square(observation_l_max + 1); ++i) {
251 : (*data_to_write_buffer)[2 * i + 1] = real(goldberg_modes.data()[i]);
252 : (*data_to_write_buffer)[2 * i + 2] = imag(goldberg_modes.data()[i]);
253 : }
254 : auto observer_proxy =
255 : Parallel::get_parallel_component<ObserverWriterComponent>(cache)[0];
256 : Parallel::threaded_action<
257 : observers::ThreadedActions::WriteReductionDataRow>(
258 : observer_proxy, "/Cce/" + detail::ScriOutput<Tag>::name(), legend,
259 : std::make_tuple(*data_to_write_buffer));
260 : }
261 : };
262 : } // namespace Actions
263 : } // namespace Cce
|
