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"
12 #include "Evolution/Systems/GrMhd/ValenciaDivClean/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"
21 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
22 #include "PointwiseFunctions/Hydro/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 grmhd {
31 namespace ValenciaDivClean {
32 
33 namespace Actions {
34 
35 /*!
36  * \brief Temporary action for observing volume and reduction data
37  */
38 struct Observe {
39  private:
43  using reduction_data = Parallel::ReductionData<
46  l2_error_datum, l2_error_datum>;
47 
48  public:
49  struct ObserveNSlabs {
50  using type = size_t;
51  static constexpr OptionString help = {"Observe every Nth slab"};
52  };
53  struct ObserveAtT0 {
54  using type = bool;
55  static constexpr OptionString help = {"If true observe at t=0"};
56  };
57 
58  using const_global_cache_tags = tmpl::list<ObserveNSlabs, ObserveAtT0>;
59  using observed_reduction_data_tags =
60  observers::make_reduction_data_tags<tmpl::list<reduction_data>>;
61 
62  template <typename... DbTags, typename... InboxTags, typename Metavariables,
63  size_t Dim, typename ActionList, typename ParallelComponent>
64  static auto apply(db::DataBox<tmpl::list<DbTags...>>& box,
67  const ElementIndex<Dim>& array_index,
68  const ActionList /*meta*/,
69  const ParallelComponent* const /*meta*/) noexcept {
70  const auto& time_id = db::get<::Tags::TimeId>(box);
71  if (time_id.substep() != 0 or (time_id.slab_number() == 0 and
72  not Parallel::get<ObserveAtT0>(cache))) {
73  return std::forward_as_tuple(std::move(box));
74  }
75 
76  const auto& time = time_id.time();
77  const std::string element_name = MakeString{} << ElementId<Dim>(array_index)
78  << '/';
79  if (time_id.slab_number() >= 0 and time_id.time().is_at_slab_start() and
80  static_cast<size_t>(time_id.slab_number()) %
81  Parallel::get<ObserveNSlabs>(cache) ==
82  0) {
83  const auto& extents = db::get<::Tags::Mesh<Dim>>(box).extents();
84  // Retrieve the tensors and compute the solution error.
85  const auto& tilde_d = db::get<Tags::TildeD>(box);
86  const auto& tilde_tau = db::get<Tags::TildeTau>(box);
87  const auto& tilde_phi = db::get<Tags::TildePhi>(box);
88  const auto& tilde_s = db::get<Tags::TildeS<>>(box);
89  const auto& tilde_b = db::get<Tags::TildeB<>>(box);
90  const auto& rest_mass_density =
91  db::get<hydro::Tags::RestMassDensity<DataVector>>(box);
92  const auto& specific_internal_energy =
93  db::get<hydro::Tags::SpecificInternalEnergy<DataVector>>(box);
94  const auto& divergence_cleaning_field =
95  db::get<hydro::Tags::DivergenceCleaningField<DataVector>>(box);
96  const auto& pressure = db::get<hydro::Tags::Pressure<DataVector>>(box);
97  const auto& lorentz_factor =
98  db::get<hydro::Tags::LorentzFactor<DataVector>>(box);
99  const auto& specific_enthalpy =
100  db::get<hydro::Tags::SpecificEnthalpy<DataVector>>(box);
101  const auto& spatial_velocity =
102  db::get<hydro::Tags::SpatialVelocity<DataVector, Dim>>(box);
103  const auto& magnetic_field =
104  db::get<hydro::Tags::MagneticField<DataVector, Dim>>(box);
105 
106  const auto& inertial_coordinates =
107  db::get<::Tags::Coordinates<Dim, Frame::Inertial>>(box);
108 
109  // Compute the error in the solution, and generate tensor component list.
110  using PrimitiveVars =
111  typename Metavariables::analytic_variables_tags;
112  using solution_tag = OptionTags::AnalyticSolutionBase;
113  const auto exact_solution = Parallel::get<solution_tag>(cache).variables(
114  inertial_coordinates, time.value(), PrimitiveVars{});
115 
116  // Remove tensor types, only storing individual components.
117  std::vector<TensorComponent> components;
118  components.reserve(34);
119 
120  components.emplace_back(element_name + Tags::TildeD::name(),
121  tilde_d.get());
122  components.emplace_back(element_name + Tags::TildeTau::name(),
123  tilde_tau.get());
124  components.emplace_back(element_name + Tags::TildePhi::name(),
125  tilde_phi.get());
126  components.emplace_back(
128  rest_mass_density.get());
129  components.emplace_back(
130  element_name +
132  specific_internal_energy.get());
133  components.emplace_back(
134  element_name +
136  divergence_cleaning_field.get());
137  components.emplace_back(
139  pressure.get());
140  components.emplace_back(
142  lorentz_factor.get());
143  components.emplace_back(
145  specific_enthalpy.get());
146 
148  DataVector error =
149  tuples::get<hydro::Tags::RestMassDensity<DataVector>>(exact_solution)
150  .get() -
151  rest_mass_density.get();
152  const double rest_mass_density_error =
153  alg::accumulate(error, 0.0, PlusSquare{});
154  components.emplace_back(
155  element_name + "Error" +
157  error);
158  error = tuples::get<hydro::Tags::SpecificInternalEnergy<DataVector>>(
159  exact_solution)
160  .get() -
161  specific_internal_energy.get();
162  const double specific_internal_energy_error =
163  alg::accumulate(error, 0.0, PlusSquare{});
164  components.emplace_back(
165  element_name + "Error" +
167  error);
168  error =
169  tuples::get<hydro::Tags::Pressure<DataVector>>(exact_solution).get() -
170  pressure.get();
171  const double pressure_error =
172  alg::accumulate(error, 0.0, PlusSquare{});
173  components.emplace_back(
174  element_name + "Error" + hydro::Tags::Pressure<DataVector>::name(),
175  error);
176  for (size_t d = 0; d < Dim; ++d) {
177  const std::string component_suffix =
178  d == 0 ? "_x" : d == 1 ? "_y" : "_z";
179  components.emplace_back(
180  element_name + Tags::TildeS<>::name() + component_suffix,
181  tilde_s.get(d));
182  components.emplace_back(
183  element_name + Tags::TildeB<>::name() + component_suffix,
184  tilde_b.get(d));
185  components.emplace_back(
186  element_name +
188  component_suffix,
189  spatial_velocity.get(d));
190  components.emplace_back(
192  component_suffix,
193  magnetic_field.get(d));
194 
195  error = tuples::get<hydro::Tags::SpatialVelocity<DataVector, Dim>>(
196  exact_solution)
197  .get(d) -
198  spatial_velocity.get(d);
199  components.emplace_back(
200  element_name + "Error" +
202  component_suffix,
203  error);
204  error = tuples::get<hydro::Tags::MagneticField<DataVector, Dim>>(
205  exact_solution)
206  .get(d) -
207  magnetic_field.get(d);
208  components.emplace_back(
209  element_name + "Error" +
211  component_suffix,
212  error);
213  components.emplace_back(
215  component_suffix,
216  inertial_coordinates.get(d));
217  }
218 
219  // Send data to volume observer
220  auto& local_observer =
221  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
222  cache)
223  .ckLocalBranch();
224  Parallel::simple_action<observers::Actions::ContributeVolumeData>(
225  local_observer, observers::ObservationId(time),
226  std::string{"/element_data"},
230  std::move(components), extents);
231 
232  // Send data to reduction observer
233  Parallel::simple_action<observers::Actions::ContributeReductionData>(
234  local_observer, observers::ObservationId(time),
235  std::string{"/element_data"},
237  "Time", "NumberOfPoints", "RestMassDensityError",
238  "SpecificInternalEnergyError", "PressureError"},
239  reduction_data{
240  time.value(),
241  db::get<::Tags::Mesh<Dim>>(box).number_of_grid_points(),
242  rest_mass_density_error, specific_internal_energy_error,
243  pressure_error});
244  }
245  return std::forward_as_tuple(std::move(box));
246  }
247 };
248 } // namespace Actions
249 } // namespace ValenciaDivClean
250 } // namespace grmhd
Defines class tuples::TaggedTuple.
The spatial velocity .
Definition: Tags.hpp:144
The fluid pressure .
Definition: Tags.hpp:123
The specific internal energy .
Definition: Tags.hpp:176
Scalar< DataType > specific_enthalpy(const Scalar< DataType > &rest_mass_density, const Scalar< DataType > &specific_internal_energy, const Scalar< DataType > &pressure) noexcept
Computes the relativistic specific enthalpy as: where is the specific internal energy...
An ID type that identifies both the parallel component and the index in the parallel component...
Definition: ArrayComponentId.hpp:27
The magnetic field measured by an Eulerian observer, where is the normal to the spatial hypersurfac...
Definition: Tags.hpp:80
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
The divergence-cleaning field .
Definition: Tags.hpp:47
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:27
Make a string by streaming into object.
Definition: MakeString.hpp:16
The densitized momentum density .
Definition: Tags.hpp:57
Temporary action for observing volume and reduction data.
Definition: Observe.hpp:38
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
The densitized magnetic field .
Definition: Tags.hpp:64
Functional for computing + of two objects.
Definition: Functional.hpp:234
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
The Lorentz factor .
Definition: Tags.hpp:64
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.
The rest-mass density .
Definition: Tags.hpp:130
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 specific enthalpy .
Definition: Tags.hpp:169
Items related to general relativistic magnetohydrodynamics (GRMHD)
Definition: Characteristics.hpp:34
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