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/TMPL.hpp"
27 
28 namespace Cce {
29 namespace Actions {
30 namespace detail {
31 // Provide a nicer name for the output h5 files for some of the uglier
32 // combinations we need
33 template <typename Tag>
34 struct ScriOutput {
35  static std::string name() noexcept { return db::tag_name<Tag>(); }
36 };
37 template <typename Tag>
38 struct ScriOutput<Tags::ScriPlus<Tag>> {
39  static std::string name() noexcept { return pretty_type::short_name<Tag>(); }
40 };
41 template <>
42 struct ScriOutput<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>> {
43  static std::string name() noexcept { return "Psi4"; }
44 };
45 
46 using weyl_correction_list =
47  tmpl::list<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>,
48  Tags::ScriPlus<Tags::Psi3>, Tags::ScriPlus<Tags::Psi2>,
49  Tags::ScriPlus<Tags::Psi1>, Tags::ScriPlus<Tags::Psi0>,
50  Tags::EthInertialRetardedTime>;
51 
52 void correct_weyl_scalars_for_inertial_time(
53  gsl::not_null<Variables<weyl_correction_list>*>
54  weyl_correction_variables) noexcept;
55 } // namespace detail
56 
57 /*!
58  * \ingroup ActionsGroup
59  * \brief Checks the interpolation managers and if they are ready, performs the
60  * interpolation and sends the data to file via
61  * `observers::ThreadedActions::WriteSimpleData`.
62  *
63  * \details This uses the `ScriPlusInterpolationManager` to perform the
64  * interpolations of all requested scri quantities (determined by
65  * `scri_values_to_observe` in the metavariables), and write them to disk using
66  * `observers::threadedActions::WriteSimpleData`.
67  *
68  * \note This action also uses the `Tags::EthInertialRetardedTime`, interpolated
69  * to the inertial frame, to perform the coordinate transformations presented in
70  * \cite Boyle:2015nqa to the Weyl scalars after interpolation. For our
71  * formulas, we need to adjust the signs and factors of two to be compatible
72  * with our definitions of \f$\eth\f$ and choice of Newman-Penrose tetrad.
73  *
74  * \f{align*}{
75  * \Psi_0^{\prime (5)}
76  * =& \Psi_0^{(5)} + 2 \eth u^\prime \Psi_1^{(4)}
77  * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_2^{(3)}
78  * + \frac{1}{2} \left( \eth u^\prime\right)^3 \Psi_3^{(2)}
79  * + \frac{1}{16} \left(\eth u^\prime\right)^4 \Psi_4^{(1)}, \\
80  * \Psi_1^{\prime (4)}
81  * =& \Psi_1^{(4)} + \frac{3}{2} \eth u^\prime \Psi_2^{(3)}
82  * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_3^{(2)}
83  * + \frac{1}{8} \left(\eth u^\prime\right)^3 \Psi_4^{(1)}, \\
84  * \Psi_2^{\prime (3)}
85  * =& \Psi_2^{(3)}
86  * + \eth u^\prime \Psi_3^{(2)}
87  * + \frac{1}{4} \left(\eth u^\prime\right)^2 \Psi_4^{(1)}, \\
88  * \Psi_3^{\prime (2)}
89  * =& \Psi_3^{(2)} + \frac{1}{2} \eth u^{\prime} \Psi_4^{ (1)}, \\
90  * \Psi_4^{\prime (1)}
91  * =& \Psi_4^{(1)}.
92  * \f}
93  *
94  * \ref DataBoxGroup changes:
95  * - Adds: nothing
96  * - Removes: nothing
97  * - Modifies: `InterpolagionManager<ComplexDataVector, Tag>` for each `Tag` in
98  * `Metavariables::scri_values_to_observe`
99  */
100 template <typename ObserverWriterComponent>
101 struct ScriObserveInterpolated {
102  using const_global_cache_tags = tmpl::list<Tags::ObservationLMax>;
103  template <typename DbTags, typename... InboxTags, typename Metavariables,
104  typename ArrayIndex, typename ActionList,
105  typename ParallelComponent>
106  static std::tuple<db::DataBox<DbTags>&&> apply(
107  db::DataBox<DbTags>& box,
108  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
109  Parallel::GlobalCache<Metavariables>& cache,
110  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
111  const ParallelComponent* const /*meta*/) noexcept {
112  const size_t observation_l_max = db::get<Tags::ObservationLMax>(box);
113  const size_t l_max = db::get<Tags::LMax>(box);
114  std::vector<double> data_to_write(2 * square(observation_l_max + 1) + 1);
115  ComplexModalVector goldberg_modes{square(l_max + 1)};
116  std::vector<std::string> file_legend;
117  file_legend.reserve(2 * square(observation_l_max + 1) + 1);
118  file_legend.emplace_back("time");
119  for (int i = 0; i <= static_cast<int>(observation_l_max); ++i) {
120  for (int j = -i; j <= i; ++j) {
121  file_legend.push_back(MakeString{} << "Real Y_" << i << "," << j);
122  file_legend.push_back(MakeString{} << "Imag Y_" << i << "," << j);
123  }
124  }
125  // alternative for the coordinate transformation getting scri+ values of the
126  // weyl scalars:
127  // need to obtain the eth of the inertial retarded time, each of the Weyl
128  // scalars, and then we'll perform a transformation on that temporary
129  // variables object, then output.
130  // it won't be as general, but that's largely fine. The main frustration is
131  // the loss of precision.
132  Variables<detail::weyl_correction_list> corrected_scri_plus_weyl{
133  Spectral::Swsh::number_of_swsh_collocation_points(l_max)};
134 
135  while (
136  db::get<Tags::InterpolationManager<
137  ComplexDataVector,
138  tmpl::front<typename Metavariables::scri_values_to_observe>>>(box)
139  .first_time_is_ready_to_interpolate()) {
140  // first get the weyl scalars and correct them
141  double interpolation_time = 0.0;
142  tmpl::for_each<detail::weyl_correction_list>([&interpolation_time,
143  &corrected_scri_plus_weyl,
144  &box](auto tag_v) noexcept {
145  using tag = typename decltype(tag_v)::type;
146  std::pair<double, ComplexDataVector> interpolation;
147  db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
148  make_not_null(&box),
149  [&interpolation](
150  const gsl::not_null<
151  ScriPlusInterpolationManager<ComplexDataVector, tag>*>
152  interpolation_manager) {
153  interpolation =
154  interpolation_manager->interpolate_and_pop_first_time();
155  });
156  interpolation_time = interpolation.first;
157  get(get<tag>(corrected_scri_plus_weyl)).data() = interpolation.second;
158  });
159 
160  detail::correct_weyl_scalars_for_inertial_time(
161  make_not_null(&corrected_scri_plus_weyl));
162 
163  // then output each of them
164  tmpl::for_each<detail::weyl_correction_list>(
165  [&data_to_write, &corrected_scri_plus_weyl, &interpolation_time,
166  &file_legend, &observation_l_max, &l_max, &cache,
167  &goldberg_modes](auto tag_v) noexcept {
168  using tag = typename decltype(tag_v)::type;
169  if constexpr (tmpl::list_contains_v<
170  typename Metavariables::scri_values_to_observe,
171  tag>) {
172  ScriObserveInterpolated::transform_and_write<
173  tag, tag::type::type::spin>(
174  get(get<tag>(corrected_scri_plus_weyl)).data(),
175  interpolation_time, make_not_null(&goldberg_modes),
176  make_not_null(&data_to_write), file_legend, l_max,
177  observation_l_max, cache);
178  }
179  });
180 
181  // then do the interpolation and output of each of the rest of the tags.
182  tmpl::for_each<
183  tmpl::list_difference<typename Metavariables::scri_values_to_observe,
184  detail::weyl_correction_list>>(
185  [&box, &data_to_write, &file_legend, &observation_l_max, &l_max,
186  &cache, &goldberg_modes](auto tag_v) noexcept {
187  using tag = typename decltype(tag_v)::type;
188  std::pair<double, ComplexDataVector> interpolation;
189  db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
190  make_not_null(&box),
191  [&interpolation](
192  const gsl::not_null<
193  ScriPlusInterpolationManager<ComplexDataVector, tag>*>
194  interpolation_manager) {
195  interpolation =
196  interpolation_manager->interpolate_and_pop_first_time();
197  });
198  ScriObserveInterpolated::transform_and_write<tag,
199  tag::type::type::spin>(
200  interpolation.second, interpolation.first,
201  make_not_null(&goldberg_modes), make_not_null(&data_to_write),
202  file_legend, l_max, observation_l_max, cache);
203  });
204  }
205  return std::forward_as_tuple(std::move(box));
206  }
207 
208  private:
209  template <typename Tag, int Spin, typename Metavariables>
210  static void transform_and_write(
211  const ComplexDataVector& data, const double time,
212  const gsl::not_null<ComplexModalVector*> goldberg_mode_buffer,
213  const gsl::not_null<std::vector<double>*> data_to_write_buffer,
214  const std::vector<std::string>& legend, const size_t l_max,
215  const size_t observation_l_max,
216  Parallel::GlobalCache<Metavariables>& cache) noexcept {
217  const SpinWeighted<ComplexDataVector, Spin> to_transform;
218  make_const_view(make_not_null(&to_transform.data()), data, 0, data.size());
219  SpinWeighted<ComplexModalVector, Spin> goldberg_modes;
220  goldberg_modes.set_data_ref(goldberg_mode_buffer);
221  Spectral::Swsh::libsharp_to_goldberg_modes(
222  make_not_null(&goldberg_modes),
223  Spectral::Swsh::swsh_transform(l_max, 1, to_transform), l_max);
224 
225  (*data_to_write_buffer)[0] = time;
226  for (size_t i = 0; i < square(observation_l_max + 1); ++i) {
227  (*data_to_write_buffer)[2 * i + 1] = real(goldberg_modes.data()[i]);
228  (*data_to_write_buffer)[2 * i + 2] = imag(goldberg_modes.data()[i]);
229  }
230  auto observer_proxy =
231  Parallel::get_parallel_component<ObserverWriterComponent>(
232  cache)[static_cast<size_t>(Parallel::my_node())];
233  Parallel::threaded_action<observers::ThreadedActions::WriteSimpleData>(
234  observer_proxy, legend, *data_to_write_buffer,
235  "/" + detail::ScriOutput<Tag>::name());
236  }
237 };
238 } // namespace Actions
239 } // 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:639
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
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:107
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:1116
Cce::Actions::ScriObserveInterpolated
Checks the interpolation managers and if they are ready, performs the interpolation and sends the dat...
Definition: ScriObserveInterpolated.hpp:101
SpinWeighted
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition,...
Definition: SpinWeighted.hpp:24
Parallel::my_node
int my_node()
Index of my node.
Definition: Info.hpp:34
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:45
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
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
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:16
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
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