ScriObserveInterpolated.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <string>
8 #include <tuple>
9 #include <utility>
10 #include <vector>
11 
12 #include "DataStructures/ComplexDataVector.hpp"
13 #include "DataStructures/DataBox/DataBox.hpp"
14 #include "DataStructures/SpinWeighted.hpp"
15 #include "Evolution/Systems/Cce/OptionTags.hpp"
16 #include "Evolution/Systems/Cce/ScriPlusInterpolationManager.hpp"
17 #include "Evolution/Systems/Cce/Tags.hpp"
18 #include "IO/Observer/WriteSimpleData.hpp"
19 #include "NumericalAlgorithms/Spectral/SwshCoefficients.hpp"
20 #include "NumericalAlgorithms/Spectral/SwshTransform.hpp"
21 #include "Parallel/GlobalCache.hpp"
22 #include "Parallel/Invoke.hpp"
23 #include "Time/Tags.hpp"
24 #include "Utilities/Gsl.hpp"
25 #include "Utilities/MakeString.hpp"
26 #include "Utilities/System/ParallelInfo.hpp"
27 #include "Utilities/TMPL.hpp"
28 
29 namespace Cce {
30 namespace Actions {
31 namespace detail {
32 // Provide a nicer name for the output h5 files for some of the uglier
33 // combinations we need
34 template <typename Tag>
35 struct ScriOutput {
36  static std::string name() noexcept { return db::tag_name<Tag>(); }
37 };
38 template <typename Tag>
39 struct ScriOutput<Tags::ScriPlus<Tag>> {
40  static std::string name() noexcept { return pretty_type::short_name<Tag>(); }
41 };
42 template <>
43 struct ScriOutput<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>> {
44  static std::string name() noexcept { return "Psi4"; }
45 };
46 
47 using weyl_correction_list =
48  tmpl::list<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>,
49  Tags::ScriPlus<Tags::Psi3>, Tags::ScriPlus<Tags::Psi2>,
50  Tags::ScriPlus<Tags::Psi1>, Tags::ScriPlus<Tags::Psi0>,
51  Tags::EthInertialRetardedTime>;
52 
53 void correct_weyl_scalars_for_inertial_time(
54  gsl::not_null<Variables<weyl_correction_list>*>
55  weyl_correction_variables) noexcept;
56 } // namespace detail
57 
58 /*!
59  * \ingroup ActionsGroup
60  * \brief Checks the interpolation managers and if they are ready, performs the
61  * interpolation and sends the data to file via
62  * `observers::ThreadedActions::WriteSimpleData`.
63  *
64  * \details This uses the `ScriPlusInterpolationManager` to perform the
65  * interpolations of all requested scri quantities (determined by
66  * `scri_values_to_observe` in the metavariables), and write them to disk using
67  * `observers::threadedActions::WriteSimpleData`.
68  *
69  * \note This action also uses the `Tags::EthInertialRetardedTime`, interpolated
70  * to the inertial frame, to perform the coordinate transformations presented in
71  * \cite Boyle:2015nqa to the Weyl scalars after interpolation. For our
72  * formulas, we need to adjust the signs and factors of two to be compatible
73  * with our definitions of \f$\eth\f$ and choice of Newman-Penrose tetrad.
74  *
75  * \f{align*}{
76  * \Psi_0^{\prime (5)}
77  * =& \Psi_0^{(5)} + 2 \eth u^\prime \Psi_1^{(4)}
78  * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_2^{(3)}
79  * + \frac{1}{2} \left( \eth u^\prime\right)^3 \Psi_3^{(2)}
80  * + \frac{1}{16} \left(\eth u^\prime\right)^4 \Psi_4^{(1)}, \\
81  * \Psi_1^{\prime (4)}
82  * =& \Psi_1^{(4)} + \frac{3}{2} \eth u^\prime \Psi_2^{(3)}
83  * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_3^{(2)}
84  * + \frac{1}{8} \left(\eth u^\prime\right)^3 \Psi_4^{(1)}, \\
85  * \Psi_2^{\prime (3)}
86  * =& \Psi_2^{(3)}
87  * + \eth u^\prime \Psi_3^{(2)}
88  * + \frac{1}{4} \left(\eth u^\prime\right)^2 \Psi_4^{(1)}, \\
89  * \Psi_3^{\prime (2)}
90  * =& \Psi_3^{(2)} + \frac{1}{2} \eth u^{\prime} \Psi_4^{ (1)}, \\
91  * \Psi_4^{\prime (1)}
92  * =& \Psi_4^{(1)}.
93  * \f}
94  *
95  * \ref DataBoxGroup changes:
96  * - Adds: nothing
97  * - Removes: nothing
98  * - Modifies: `InterpolagionManager<ComplexDataVector, Tag>` for each `Tag` in
99  * `Metavariables::scri_values_to_observe`
100  */
101 template <typename ObserverWriterComponent>
102 struct ScriObserveInterpolated {
103  using const_global_cache_tags = tmpl::list<Tags::ObservationLMax>;
104  template <typename DbTags, typename... InboxTags, typename Metavariables,
105  typename ArrayIndex, typename ActionList,
106  typename ParallelComponent>
107  static std::tuple<db::DataBox<DbTags>&&> apply(
108  db::DataBox<DbTags>& box,
109  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
110  Parallel::GlobalCache<Metavariables>& cache,
111  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
112  const ParallelComponent* const /*meta*/) noexcept {
113  const size_t observation_l_max = db::get<Tags::ObservationLMax>(box);
114  const size_t l_max = db::get<Tags::LMax>(box);
115  std::vector<double> data_to_write(2 * square(observation_l_max + 1) + 1);
116  ComplexModalVector goldberg_modes{square(l_max + 1)};
117  std::vector<std::string> file_legend;
118  file_legend.reserve(2 * square(observation_l_max + 1) + 1);
119  file_legend.emplace_back("time");
120  for (int i = 0; i <= static_cast<int>(observation_l_max); ++i) {
121  for (int j = -i; j <= i; ++j) {
122  file_legend.push_back(MakeString{} << "Real Y_" << i << "," << j);
123  file_legend.push_back(MakeString{} << "Imag Y_" << i << "," << j);
124  }
125  }
126  // alternative for the coordinate transformation getting scri+ values of the
127  // weyl scalars:
128  // need to obtain the eth of the inertial retarded time, each of the Weyl
129  // scalars, and then we'll perform a transformation on that temporary
130  // variables object, then output.
131  // it won't be as general, but that's largely fine. The main frustration is
132  // the loss of precision.
133  Variables<detail::weyl_correction_list> corrected_scri_plus_weyl{
134  Spectral::Swsh::number_of_swsh_collocation_points(l_max)};
135 
136  while (
137  db::get<Tags::InterpolationManager<
138  ComplexDataVector,
139  tmpl::front<typename Metavariables::scri_values_to_observe>>>(box)
140  .first_time_is_ready_to_interpolate()) {
141  // first get the weyl scalars and correct them
142  double interpolation_time = 0.0;
143  tmpl::for_each<detail::weyl_correction_list>([&interpolation_time,
144  &corrected_scri_plus_weyl,
145  &box](auto tag_v) noexcept {
146  using tag = typename decltype(tag_v)::type;
147  std::pair<double, ComplexDataVector> interpolation;
148  db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
149  make_not_null(&box),
150  [&interpolation](
151  const gsl::not_null<
152  ScriPlusInterpolationManager<ComplexDataVector, tag>*>
153  interpolation_manager) {
154  interpolation =
155  interpolation_manager->interpolate_and_pop_first_time();
156  });
157  interpolation_time = interpolation.first;
158  get(get<tag>(corrected_scri_plus_weyl)).data() = interpolation.second;
159  });
160 
161  detail::correct_weyl_scalars_for_inertial_time(
162  make_not_null(&corrected_scri_plus_weyl));
163 
164  // then output each of them
165  tmpl::for_each<detail::weyl_correction_list>(
166  [&data_to_write, &corrected_scri_plus_weyl, &interpolation_time,
167  &file_legend, &observation_l_max, &l_max, &cache,
168  &goldberg_modes](auto tag_v) noexcept {
169  using tag = typename decltype(tag_v)::type;
170  if constexpr (tmpl::list_contains_v<
171  typename Metavariables::scri_values_to_observe,
172  tag>) {
173  ScriObserveInterpolated::transform_and_write<
174  tag, tag::type::type::spin>(
175  get(get<tag>(corrected_scri_plus_weyl)).data(),
176  interpolation_time, make_not_null(&goldberg_modes),
177  make_not_null(&data_to_write), file_legend, l_max,
178  observation_l_max, cache);
179  }
180  });
181 
182  // then do the interpolation and output of each of the rest of the tags.
183  tmpl::for_each<
184  tmpl::list_difference<typename Metavariables::scri_values_to_observe,
185  detail::weyl_correction_list>>(
186  [&box, &data_to_write, &file_legend, &observation_l_max, &l_max,
187  &cache, &goldberg_modes](auto tag_v) noexcept {
188  using tag = typename decltype(tag_v)::type;
189  std::pair<double, ComplexDataVector> interpolation;
190  db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
191  make_not_null(&box),
192  [&interpolation](
193  const gsl::not_null<
194  ScriPlusInterpolationManager<ComplexDataVector, tag>*>
195  interpolation_manager) {
196  interpolation =
197  interpolation_manager->interpolate_and_pop_first_time();
198  });
199  ScriObserveInterpolated::transform_and_write<tag,
200  tag::type::type::spin>(
201  interpolation.second, interpolation.first,
202  make_not_null(&goldberg_modes), make_not_null(&data_to_write),
203  file_legend, l_max, observation_l_max, cache);
204  });
205  }
206  return std::forward_as_tuple(std::move(box));
207  }
208 
209  private:
210  template <typename Tag, int Spin, typename Metavariables>
211  static void transform_and_write(
212  const ComplexDataVector& data, const double time,
213  const gsl::not_null<ComplexModalVector*> goldberg_mode_buffer,
214  const gsl::not_null<std::vector<double>*> data_to_write_buffer,
215  const std::vector<std::string>& legend, const size_t l_max,
216  const size_t observation_l_max,
217  Parallel::GlobalCache<Metavariables>& cache) noexcept {
218  const SpinWeighted<ComplexDataVector, Spin> to_transform;
219  make_const_view(make_not_null(&to_transform.data()), data, 0, data.size());
220  SpinWeighted<ComplexModalVector, Spin> goldberg_modes;
221  goldberg_modes.set_data_ref(goldberg_mode_buffer);
222  Spectral::Swsh::libsharp_to_goldberg_modes(
223  make_not_null(&goldberg_modes),
224  Spectral::Swsh::swsh_transform(l_max, 1, to_transform), l_max);
225 
226  (*data_to_write_buffer)[0] = time;
227  for (size_t i = 0; i < square(observation_l_max + 1); ++i) {
228  (*data_to_write_buffer)[2 * i + 1] = real(goldberg_modes.data()[i]);
229  (*data_to_write_buffer)[2 * i + 2] = imag(goldberg_modes.data()[i]);
230  }
231  auto observer_proxy =
232  Parallel::get_parallel_component<ObserverWriterComponent>(
233  cache)[static_cast<size_t>(sys::my_node())];
234  Parallel::threaded_action<observers::ThreadedActions::WriteSimpleData>(
235  observer_proxy, legend, *data_to_write_buffer,
236  "/" + detail::ScriOutput<Tag>::name());
237  }
238 };
239 } // namespace Actions
240 } // namespace Cce
std::string
utility
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
std::pair
GlobalCache.hpp
vector
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
Cce::Tags::InterpolationManager
Definition: Tags.hpp:374
Spectral::Swsh::swsh_transform
SpinWeighted< ComplexModalVector, Spin > swsh_transform(const size_t l_max, const size_t number_of_radial_points, const SpinWeighted< ComplexDataVector, Spin > &collocation) noexcept
Perform a forward libsharp spin-weighted spherical harmonic transform on a single supplied SpinWeight...
Definition: SwshTransform.cpp:106
ParallelInfo.hpp
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:969
Cce::Actions::ScriObserveInterpolated
Checks the interpolation managers and if they are ready, performs the interpolation and sends the dat...
Definition: ScriObserveInterpolated.hpp:102
SpinWeighted
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition,...
Definition: SpinWeighted.hpp:24
DataBox.hpp
cstddef
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
ComplexModalVector
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:39
Cce
The set of utilities for performing Cauchy characteristic evolution and Cauchy characteristic matchin...
Definition: BoundaryComputeAndSendToEvolution.hpp:28
Gsl.hpp
Cce::ScriPlusInterpolationManager
Stores necessary data and interpolates on to new time points at scri+.
Definition: ScriPlusInterpolationManager.hpp:48
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
ComplexDataVector
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:47
make_const_view
void make_const_view(const gsl::not_null< const SpinWeightedType * > view, const SpinWeightedType &spin_weighted, const size_t offset, const size_t extent) noexcept
Definition: SpinWeighted.hpp:544
TMPL.hpp
MakeString
Make a string by streaming into object.
Definition: MakeString.hpp:18
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
Spectral::Swsh::number_of_swsh_collocation_points
constexpr size_t number_of_swsh_collocation_points(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonic collocation value...
Definition: SwshCollocation.hpp:25
string
© Copyright 2017 - 2020 SXS Collaboration, MIT License