Observe.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 
11 #include "Domain/Tags.hpp"
13 #include "IO/Observer/Actions.hpp"
14 #include "IO/Observer/ArrayComponentId.hpp"
15 #include "IO/Observer/Helpers.hpp"
16 #include "IO/Observer/ObservationId.hpp"
17 #include "IO/Observer/ObserverComponent.hpp"
18 #include "IO/Observer/ReductionActions.hpp"
19 #include "IO/Observer/VolumeActions.hpp"
20 #include "Options/Options.hpp"
22 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
23 #include "Time/Tags.hpp"
24 #include "Utilities/Functional.hpp"
25 #include "Utilities/MakeString.hpp"
26 #include "Utilities/Numeric.hpp"
27 #include "Utilities/TMPL.hpp"
29 
30 namespace ScalarWave {
31 namespace Actions {
32 
33 /*!
34  * \brief Temporary action for observing volume and reduction data
35  *
36  * A few notes:
37  * - Writes the solution and error in \f$\Psi, \Pi\f$, and \f$\Phi_i\f$ to disk
38  * as volume data.
39  * - The RMS error of \f$\Psi\f$ and \f$\Pi\f$ are written to disk.
40  */
41 struct Observe {
42  private:
46  using reduction_data = Parallel::ReductionData<
49  l2_error_datum>;
50 
51  public:
52  struct ObserveNSlabs {
53  using type = size_t;
54  static constexpr OptionString help = {"Observe every Nth slab"};
55  };
56  struct ObserveAtT0 {
57  using type = bool;
58  static constexpr OptionString help = {"If true observe at t=0"};
59  };
60 
61  using const_global_cache_tags = tmpl::list<ObserveNSlabs, ObserveAtT0>;
62  using observed_reduction_data_tags =
63  observers::make_reduction_data_tags<tmpl::list<reduction_data>>;
64 
65  template <typename... DbTags, typename... InboxTags, typename Metavariables,
66  size_t Dim, typename ActionList, typename ParallelComponent>
67  static auto apply(db::DataBox<tmpl::list<DbTags...>>& box,
70  const ElementIndex<Dim>& array_index,
71  const ActionList /*meta*/,
72  const ParallelComponent* const /*meta*/) noexcept {
73  const auto& time_id = db::get<Tags::TimeId>(box);
74  if (time_id.substep() != 0 or (time_id.slab_number() == 0 and
75  not Parallel::get<ObserveAtT0>(cache))) {
76  return std::forward_as_tuple(std::move(box));
77  }
78 
79  const auto& time = time_id.time();
80  const std::string element_name = MakeString{} << ElementId<Dim>(array_index)
81  << '/';
82 
83  if (time_id.slab_number() >= 0 and time_id.time().is_at_slab_start() and
84  static_cast<size_t>(time_id.slab_number()) %
85  Parallel::get<ObserveNSlabs>(cache) ==
86  0) {
87  const auto& extents = db::get<Tags::Mesh<Dim>>(box).extents();
88  // Retrieve the tensors and compute the solution error.
89  const auto& psi = db::get<ScalarWave::Psi>(box);
90  const auto& pi = db::get<ScalarWave::Pi>(box);
91  const auto& phi = db::get<ScalarWave::Phi<Dim>>(box);
92  const auto& inertial_coordinates =
93  db::get<Tags::Coordinates<Dim, Frame::Inertial>>(box);
94 
95  // Compute the error in the solution, and generate tensor component list.
96  using Vars = typename Metavariables::system::variables_tag::type;
97  using solution_tag = OptionTags::AnalyticSolutionBase;
98  const auto exact_solution = Parallel::get<solution_tag>(cache).variables(
99  inertial_coordinates, time.value(), typename Vars::tags_list{});
100 
101  // Remove tensor types, only storing individual components.
102  std::vector<TensorComponent> components;
103  components.reserve(3 * Dim + 4);
104 
105  components.emplace_back(element_name + ScalarWave::Psi::name(),
106  psi.get());
108  DataVector error =
109  tuples::get<ScalarWave::Psi>(exact_solution).get() - psi.get();
110  const double psi_error = alg::accumulate(error, 0.0, PlusSquare{});
111  components.emplace_back(element_name + "Error" + ScalarWave::Psi::name(),
112  error);
113  components.emplace_back(element_name + ScalarWave::Pi::name(), pi.get());
114  error = tuples::get<ScalarWave::Pi>(exact_solution).get() - pi.get();
115  const double pi_error = alg::accumulate(error, 0.0, PlusSquare{});
116  components.emplace_back(element_name + "Error" + ScalarWave::Pi::name(),
117  error);
118  for (size_t d = 0; d < Dim; ++d) {
119  const std::string component_suffix =
120  d == 0 ? "_x" : d == 1 ? "_y" : "_z";
121  components.emplace_back(
122  element_name + ScalarWave::Phi<Dim>::name() + component_suffix,
123  phi.get(d));
124 
125  error = tuples::get<ScalarWave::Phi<Dim>>(exact_solution).get(d) -
126  phi.get(d);
127  components.emplace_back(element_name + "Error" +
129  component_suffix,
130  error);
131  components.emplace_back(
133  component_suffix,
134  inertial_coordinates.get(d));
135  }
136 
137  // Send data to volume observer
138  auto& local_observer =
139  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
140  cache)
141  .ckLocalBranch();
142  Parallel::simple_action<observers::Actions::ContributeVolumeData>(
143  local_observer, observers::ObservationId(time),
144  std::string{"/element_data"},
148  std::move(components), extents);
149 
150  // Send data to reduction observer
151  Parallel::simple_action<observers::Actions::ContributeReductionData>(
152  local_observer, observers::ObservationId(time),
153  std::string{"/element_data"},
154  std::vector<std::string>{"Time", "NumberOfPoints", "PsiError",
155  "PiError"},
156  reduction_data{time.value(),
157  db::get<Tags::Mesh<Dim>>(box).number_of_grid_points(),
158  psi_error, pi_error});
159  }
160  return std::forward_as_tuple(std::move(box));
161  }
162 };
163 } // namespace Actions
164 } // namespace ScalarWave
Defines class tuples::TaggedTuple.
Definition: Tags.hpp:29
Defines classes and functions for making classes creatable from input files.
Temporary action for observing volume and reduction data.
Definition: Observe.hpp:41
An ID type that identifies both the parallel component and the index in the parallel component...
Definition: ArrayComponentId.hpp:27
decltype(auto) accumulate(const Container &c, T init)
Convenience wrapper around std::accumulate, returns std::accumulate(begin(c), end(c), init).
Definition: Numeric.hpp:54
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
tnsr::aa< DataType, SpatialDim, Frame > pi(const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ii< DataType, SpatialDim, Frame > &dt_spatial_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the conjugate momentum of the spacetime metric .
Definition: ComputeGhQuantities.cpp:58
Make a string by streaming into object.
Definition: MakeString.hpp:16
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines DataBox tags for scalar wave system.
Items related to evolving the scalar wave equation:
Definition: Equations.cpp:20
Defines classes and functions used for manipulating DataBox&#39;s.
Functional for computing + of two objects.
Definition: Functional.hpp:234
tnsr::iaa< DataType, SpatialDim, Frame > phi(const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein&#39;s equations...
Definition: ComputeGhQuantities.cpp:22
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:76
The coordinates in a given frame.
Definition: Tags.hpp:95
Can be used to retrieve the analytic solution from the cache without having to know the template para...
Definition: Tags.hpp:12
Stores a collection of function values.
Definition: DataVector.hpp:46
Functional for computing sqrt on an object.
Definition: Functional.hpp:269
Wraps the template metaprogramming library used (brigand)
Defines tags related to domain quantities.
Definition: SolvePoissonProblem.hpp:38
A class for indexing a Charm array by Element.
Definition: ElementIndex.hpp:53
A type-erased identifier that combines the identifier&#39;s type and hash used to uniquely identify an ob...
Definition: ObservationId.hpp:43
Defines class template ConstGlobalCache.
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:64
Defines tags related to Time quantities.
Definition: ComputeTimeDerivative.hpp:28
The array index used for indexing Chare Arrays, mostly an implementation detail.
Definition: ArrayIndex.hpp:20