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"
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"
27 #include "Utilities/TMPL.hpp"
28 
29 namespace Cce {
30 /// \cond
31 template <typename Metavariables>
32 struct AnalyticWorldtubeBoundary;
33 /// \endcond
34 namespace Actions {
35 namespace detail {
36 // Provide a nicer name for the output h5 files for some of the uglier
37 // combinations we need
38 template <typename Tag>
39 struct ScriOutput {
40  static std::string name() noexcept { return db::tag_name<Tag>(); }
41 };
42 template <typename Tag>
43 struct ScriOutput<Tags::ScriPlus<Tag>> {
44  static std::string name() noexcept { return pretty_type::short_name<Tag>(); }
45 };
46 template <>
47 struct ScriOutput<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>> {
48  static std::string name() noexcept { return "Psi4"; }
49 };
50 
51 using weyl_correction_list =
52  tmpl::list<Tags::Du<Tags::TimeIntegral<Tags::ScriPlus<Tags::Psi4>>>,
53  Tags::ScriPlus<Tags::Psi3>, Tags::ScriPlus<Tags::Psi2>,
54  Tags::ScriPlus<Tags::Psi1>, Tags::ScriPlus<Tags::Psi0>,
55  Tags::EthInertialRetardedTime>;
56 
57 void correct_weyl_scalars_for_inertial_time(
58  gsl::not_null<Variables<weyl_correction_list>*>
59  weyl_correction_variables) noexcept;
60 } // namespace detail
61 
62 /*!
63  * \ingroup ActionsGroup
64  * \brief Checks the interpolation managers and if they are ready, performs the
65  * interpolation and sends the data to file via
66  * `observers::ThreadedActions::WriteSimpleData`.
67  *
68  * \details This uses the `ScriPlusInterpolationManager` to perform the
69  * interpolations of all requested scri quantities (determined by
70  * `scri_values_to_observe` in the metavariables), and write them to disk using
71  * `observers::threadedActions::WriteSimpleData`. When using an analytic
72  * worldtube solution, this action also uses the `AnalyticBoundaryDataManager`
73  * to output the expected News value at the appropriate asymptotically inertial
74  * time.
75  *
76  * \note This action also uses the `Tags::EthInertialRetardedTime`, interpolated
77  * to the inertial frame, to perform the coordinate transformations presented in
78  * \cite Boyle:2015nqa to the Weyl scalars after interpolation. For our
79  * formulas, we need to adjust the signs and factors of two to be compatible
80  * with our definitions of \f$\eth\f$ and choice of Newman-Penrose tetrad.
81  *
82  * \f{align*}{
83  * \Psi_0^{\prime (5)}
84  * =& \Psi_0^{(5)} + 2 \eth u^\prime \Psi_1^{(4)}
85  * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_2^{(3)}
86  * + \frac{1}{2} \left( \eth u^\prime\right)^3 \Psi_3^{(2)}
87  * + \frac{1}{16} \left(\eth u^\prime\right)^4 \Psi_4^{(1)}, \\
88  * \Psi_1^{\prime (4)}
89  * =& \Psi_1^{(4)} + \frac{3}{2} \eth u^\prime \Psi_2^{(3)}
90  * + \frac{3}{4} \left(\eth u^\prime\right)^2 \Psi_3^{(2)}
91  * + \frac{1}{8} \left(\eth u^\prime\right)^3 \Psi_4^{(1)}, \\
92  * \Psi_2^{\prime (3)}
93  * =& \Psi_2^{(3)}
94  * + \eth u^\prime \Psi_3^{(2)}
95  * + \frac{1}{4} \left(\eth u^\prime\right)^2 \Psi_4^{(1)}, \\
96  * \Psi_3^{\prime (2)}
97  * =& \Psi_3^{(2)} + \frac{1}{2} \eth u^{\prime} \Psi_4^{ (1)}, \\
98  * \Psi_4^{\prime (1)}
99  * =& \Psi_4^{(1)}.
100  * \f}
101  *
102  * \ref DataBoxGroup changes:
103  * - Adds: nothing
104  * - Removes: nothing
105  * - Modifies: `InterpolagionManager<ComplexDataVector, Tag>` for each `Tag` in
106  * `Metavariables::scri_values_to_observe`
107  */
108 template <typename ObserverWriterComponent, typename BoundaryComponent>
110  using const_global_cache_tags = tmpl::flatten<
111  tmpl::list<Tags::ObservationLMax,
113  tt::is_a_v<AnalyticWorldtubeBoundary, BoundaryComponent>,
114  tmpl::list<Tags::OutputNoninertialNews>, tmpl::list<>>>>;
115  template <typename DbTags, typename... InboxTags, typename Metavariables,
116  typename ArrayIndex, typename ActionList,
117  typename ParallelComponent>
118  static std::tuple<db::DataBox<DbTags>&&> apply(
119  db::DataBox<DbTags>& box,
120  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
122  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
123  const ParallelComponent* const /*meta*/) noexcept {
124  const size_t observation_l_max = db::get<Tags::ObservationLMax>(box);
125  const size_t l_max = db::get<Tags::LMax>(box);
126  std::vector<double> data_to_write(2 * square(observation_l_max + 1) + 1);
127  ComplexModalVector goldberg_modes{square(l_max + 1)};
128  std::vector<std::string> file_legend;
129  file_legend.reserve(2 * square(observation_l_max + 1) + 1);
130  file_legend.emplace_back("time");
131  for (int i = 0; i <= static_cast<int>(observation_l_max); ++i) {
132  for (int j = -i; j <= i; ++j) {
133  file_legend.push_back(MakeString{} << "Real Y_" << i << "," << j);
134  file_legend.push_back(MakeString{} << "Imag Y_" << i << "," << j);
135  }
136  }
137  // alternative for the coordinate transformation getting scri+ values of the
138  // weyl scalars:
139  // need to obtain the eth of the inertial retarded time, each of the Weyl
140  // scalars, and then we'll perform a transformation on that temporary
141  // variables object, then output.
142  // it won't be as general, but that's largely fine. The main frustration is
143  // the loss of precision.
144  Variables<detail::weyl_correction_list> corrected_scri_plus_weyl{
146 
147  while (
150  tmpl::front<typename Metavariables::scri_values_to_observe>>>(box)
151  .first_time_is_ready_to_interpolate()) {
152  // first get the weyl scalars and correct them
153  double interpolation_time = 0.0;
154  tmpl::for_each<detail::weyl_correction_list>([&interpolation_time,
155  &corrected_scri_plus_weyl,
156  &box](auto tag_v) noexcept {
157  using tag = typename decltype(tag_v)::type;
159  db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
160  make_not_null(&box),
161  [&interpolation](
162  const gsl::not_null<
164  interpolation_manager) {
165  interpolation =
166  interpolation_manager->interpolate_and_pop_first_time();
167  });
168  interpolation_time = interpolation.first;
169  get(get<tag>(corrected_scri_plus_weyl)).data() = interpolation.second;
170  });
171 
172  detail::correct_weyl_scalars_for_inertial_time(
173  make_not_null(&corrected_scri_plus_weyl));
174 
175  // then output each of them
176  tmpl::for_each<detail::weyl_correction_list>(
177  [&data_to_write, &corrected_scri_plus_weyl, &interpolation_time,
178  &file_legend, &observation_l_max, &l_max, &cache,
179  &goldberg_modes](auto tag_v) noexcept {
180  using tag = typename decltype(tag_v)::type;
181  if constexpr (tmpl::list_contains_v<
182  typename Metavariables::scri_values_to_observe,
183  tag>) {
184  ScriObserveInterpolated::transform_and_write<
185  tag, tag::type::type::spin, ParallelComponent>(
186  get(get<tag>(corrected_scri_plus_weyl)).data(),
187  interpolation_time, make_not_null(&goldberg_modes),
188  make_not_null(&data_to_write), file_legend, l_max,
189  observation_l_max, cache);
190  }
191  });
192 
193  // then do the interpolation and output of each of 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, &file_legend, &observation_l_max, &l_max,
198  &cache, &goldberg_modes](auto tag_v) noexcept {
199  using tag = typename decltype(tag_v)::type;
201  db::mutate<Tags::InterpolationManager<ComplexDataVector, tag>>(
202  make_not_null(&box),
203  [&interpolation](
204  const gsl::not_null<
206  interpolation_manager) {
207  interpolation =
208  interpolation_manager->interpolate_and_pop_first_time();
209  });
210  ScriObserveInterpolated::transform_and_write<
211  tag, tag::type::type::spin, ParallelComponent>(
212  interpolation.second, interpolation.first,
213  make_not_null(&goldberg_modes), make_not_null(&data_to_write),
214  file_legend, l_max, observation_l_max, cache);
215  });
216  // output the expected news associated with the time value interpolated at
217  // scri+.
218  if constexpr (tmpl::list_contains_v<DbTags,
220  if (not db::get<Tags::OutputNoninertialNews>(box)) {
221  db::get<Tags::AnalyticBoundaryDataManager>(box)
222  .template write_news<ParallelComponent>(cache,
223  interpolation_time);
224  }
225  }
226  }
227  return std::forward_as_tuple(std::move(box));
228  }
229 
230  private:
231  template <typename Tag, int Spin, typename ParallelComponent,
232  typename Metavariables>
233  static void transform_and_write(
234  const ComplexDataVector& data, const double time,
235  const gsl::not_null<ComplexModalVector*> goldberg_mode_buffer,
236  const gsl::not_null<std::vector<double>*> data_to_write_buffer,
237  const std::vector<std::string>& legend, const size_t l_max,
238  const size_t observation_l_max,
240  const SpinWeighted<ComplexDataVector, Spin> to_transform;
241  make_const_view(make_not_null(&to_transform.data()), data, 0, data.size());
243  goldberg_modes.set_data_ref(goldberg_mode_buffer);
245  make_not_null(&goldberg_modes),
246  Spectral::Swsh::swsh_transform(l_max, 1, to_transform), l_max);
247 
248  (*data_to_write_buffer)[0] = time;
249  for (size_t i = 0; i < square(observation_l_max + 1); ++i) {
250  (*data_to_write_buffer)[2 * i + 1] = real(goldberg_modes.data()[i]);
251  (*data_to_write_buffer)[2 * i + 2] = imag(goldberg_modes.data()[i]);
252  }
253  auto& my_proxy =
254  Parallel::get_parallel_component<ParallelComponent>(cache);
255  auto observer_proxy =
256  Parallel::get_parallel_component<ObserverWriterComponent>(
257  cache)[static_cast<size_t>(
258  Parallel::my_node(*my_proxy.ckLocal()))];
259  Parallel::threaded_action<observers::ThreadedActions::WriteSimpleData>(
260  observer_proxy, legend, *data_to_write_buffer,
261  "/" + detail::ScriOutput<Tag>::name());
262  }
263 };
264 } // namespace Actions
265 } // 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:660
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:367
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:786
Cce::Actions::ScriObserveInterpolated
Checks the interpolation managers and if they are ready, performs the interpolation and sends the dat...
Definition: ScriObserveInterpolated.hpp:109
SpinWeighted
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition,...
Definition: SpinWeighted.hpp:24
Cce::Tags::ObservationLMax
Definition: OptionTags.hpp:307
Parallel::my_node
int my_node(const DistribObject &distributed_object) noexcept
Index of my node.
Definition: Info.hpp:51
DataBox.hpp
cstddef
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
Spectral::Swsh::libsharp_to_goldberg_modes
void libsharp_to_goldberg_modes(gsl::not_null< SpinWeighted< ComplexModalVector, Spin > * > goldberg_modes, const SpinWeighted< ComplexModalVector, Spin > &libsharp_modes, size_t l_max) noexcept
Compute the set of Goldberg Spin-weighted spherical harmonic modes (in the convention of ) from a lib...
Cce::Tags::AnalyticBoundaryDataManager
A tag that constructs a AnalyticBoundaryDataManager from options.
Definition: OptionTags.hpp:529
ComplexModalVector
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:40
Spectral::Swsh::swsh_transform
void swsh_transform(const size_t l_max, const size_t number_of_radial_points, const gsl::not_null< SpinWeighted< ComplexModalVector, Spin > * > first_coefficient, const ModalThenNodalTypes &... coefficients_then_collocations) noexcept
Perform a forward libsharp spin-weighted spherical harmonic transform on any number of supplied SpinW...
Definition: SwshTransform.hpp:171
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:380
Cce
The set of utilities for performing Cauchy characteristic evolution and Cauchy characteristic matchin...
Definition: CharacteristicExtractFwd.hpp:6
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
std::conditional_t
ComplexDataVector
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:53
brigand::list_difference
fold< Sequence2, Sequence1, lazy::remove< _state, _element > > list_difference
Obtain the elements of Sequence1 that are not in Sequence2.
Definition: TMPL.hpp:589
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