ObserveFields.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <functional>
8 #include <initializer_list>
9 #include <optional>
10 #include <pup.h>
11 #include <string>
12 #include <type_traits>
13 #include <unordered_set>
14 #include <utility>
15 #include <vector>
16 
18 #include "DataStructures/DataBox/TagName.hpp"
19 #include "DataStructures/DataVector.hpp"
20 #include "DataStructures/Tensor/TensorData.hpp"
22 #include "Domain/Tags.hpp"
23 #include "IO/Observer/ArrayComponentId.hpp"
24 #include "IO/Observer/ObservationId.hpp"
25 #include "IO/Observer/ObserverComponent.hpp" // IWYU pragma: keep
26 #include "IO/Observer/VolumeActions.hpp" // IWYU pragma: keep
27 #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
28 #include "Options/Auto.hpp"
29 #include "Options/Options.hpp"
30 #include "Parallel/ArrayIndex.hpp"
32 #include "Parallel/GlobalCache.hpp"
33 #include "Parallel/Invoke.hpp"
34 #include "Parallel/PupStlCpp17.hpp"
35 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
36 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
37 #include "Utilities/Algorithm.hpp"
38 #include "Utilities/Literals.hpp"
39 #include "Utilities/MakeString.hpp"
40 #include "Utilities/Numeric.hpp"
41 #include "Utilities/StdHelpers.hpp"
42 #include "Utilities/TMPL.hpp"
43 #include "Utilities/TypeTraits/IsA.hpp"
44 
45 /// \cond
46 template <size_t Dim>
47 class Mesh;
48 namespace Frame {
49 struct Inertial;
50 } // namespace Frame
51 namespace Tags {
52 struct Time;
53 } // namespace Tags
54 /// \endcond
55 
56 namespace dg {
57 namespace Events {
58 /// \cond
59 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
60  typename AnalyticSolutionTensors = tmpl::list<>,
61  typename NonSolutionTensors =
63 class ObserveFields;
64 /// \endcond
65 
66 /*!
67  * \ingroup DiscontinuousGalerkinGroup
68  * \brief %Observe volume tensor fields.
69  *
70  * A class that writes volume quantities to an h5 file during the simulation.
71  * The observed quantitites are:
72  * - `InertialCoordinates`
73  * - Tensors listed in `Tensors` template parameter
74  * - `Error(*)` = errors in `AnalyticSolutionTensors` =
75  * \f$\text{value} - \text{analytic solution}\f$
76  *
77  * The user may specify an `interpolation_mesh` to which the
78  * data is interpolated.
79  */
80 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
81  typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
82 class ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
83  tmpl::list<AnalyticSolutionTensors...>,
84  tmpl::list<NonSolutionTensors...>> : public Event {
85  private:
86  static_assert(
87  std::is_same_v<
88  tmpl::list_difference<tmpl::list<AnalyticSolutionTensors...>,
89  tmpl::list<Tensors...>>,
90  tmpl::list<>>,
91  "All AnalyticSolutionTensors must be listed in Tensors.");
93 
94  public:
95  /// The name of the subfile inside the HDF5 file
96  struct SubfileName {
97  using type = std::string;
98  static constexpr Options::String help = {
99  "The name of the subfile inside the HDF5 file without an extension and "
100  "without a preceding '/'."};
101  };
102 
103  /// \cond
104  explicit ObserveFields(CkMigrateMessage* /*unused*/) noexcept {}
105  using PUP::able::register_constructor;
106  WRAPPED_PUPable_decl_template(ObserveFields); // NOLINT
107  /// \endcond
108 
109  struct VariablesToObserve {
110  static constexpr Options::String help = "Subset of variables to observe";
112  static size_t lower_bound_on_size() noexcept { return 1; }
113  };
114 
115  struct InterpolateToMesh {
116  using type = Options::Auto<Mesh<VolumeDim>, Options::AutoLabel::None>;
117  static constexpr Options::String help =
118  "An optional mesh to which the variables are interpolated. This mesh "
119  "specifies any number of collocation points, basis, and quadrature on "
120  "which the observed quantities are evaluated. If no mesh is given, the "
121  "results will be evaluated on the mesh the simulation runs on. The "
122  "user may add several ObserveField Events e.g. with and without an "
123  "interpolating mesh to output the data both on the original mesh and "
124  "on a new mesh.";
125  };
126 
127  using options =
128  tmpl::list<SubfileName, VariablesToObserve, InterpolateToMesh>;
129  static constexpr Options::String help =
130  "Observe volume tensor fields.\n"
131  "\n"
132  "Writes volume quantities:\n"
133  " * InertialCoordinates\n"
134  " * Tensors listed in Tensors template parameter\n"
135  " * Error(*) = errors in AnalyticSolutionTensors\n"
136  " = value - analytic solution\n"
137  "\n"
138  "Warning: Currently, only one volume observation event can be\n"
139  "triggered at a given time. Causing multiple events to run at once\n"
140  "will produce unpredictable results.";
141 
142  ObserveFields() = default;
143 
144  explicit ObserveFields(const std::string& subfile_name,
145  const std::vector<std::string>& variables_to_observe,
146  std::optional<Mesh<VolumeDim>> interpolation_mesh = {},
147  const Options::Context& context = {})
148  : subfile_path_("/" + subfile_name),
149  variables_to_observe_(variables_to_observe.begin(),
150  variables_to_observe.end()),
151  interpolation_mesh_(interpolation_mesh) {
152  using ::operator<<;
153  const std::unordered_set<std::string> valid_tensors{
154  db::tag_name<Tensors>()...};
155  for (const auto& name : variables_to_observe_) {
156  if (valid_tensors.count(name) != 1) {
157  PARSE_ERROR(
158  context,
159  name << " is not an available variable. Available variables:\n"
160  << (std::vector<std::string>{db::tag_name<Tensors>()...}));
161  }
162  if (alg::count(variables_to_observe, name) != 1) {
163  PARSE_ERROR(context, name << " specified multiple times");
164  }
165  }
166  variables_to_observe_.insert(coordinates_tag::name());
167  }
168 
169  using argument_tags = tmpl::flatten<tmpl::list<
170  ObservationValueTag, domain::Tags::Mesh<VolumeDim>, coordinates_tag,
171  AnalyticSolutionTensors..., NonSolutionTensors...,
172  tmpl::conditional_t<(sizeof...(AnalyticSolutionTensors) > 0),
173  ::Tags::AnalyticSolutionsBase, tmpl::list<>>>>;
174 
175  template <typename OptionalAnalyticSolutions, typename Metavariables,
176  typename ParallelComponent>
177  void operator()(
178  const typename ObservationValueTag::type& observation_value,
179  const Mesh<VolumeDim>& mesh,
180  const tnsr::I<DataVector, VolumeDim, Frame::Inertial>&
181  inertial_coordinates,
182  const typename AnalyticSolutionTensors::
183  type&... analytic_solution_tensors,
184  const typename NonSolutionTensors::type&... non_solution_tensors,
185  const OptionalAnalyticSolutions& optional_analytic_solutions,
187  const ElementId<VolumeDim>& array_index,
188  const ParallelComponent* const component) const noexcept {
189  call_operator_impl(subfile_path_, variables_to_observe_,
190  interpolation_mesh_, observation_value, mesh,
191  inertial_coordinates, analytic_solution_tensors...,
192  non_solution_tensors..., optional_analytic_solutions,
193  cache, array_index, component);
194  }
195 
196  // This overload is called when the list of analytic-solution tensors is
197  // empty, i.e. it is clear at compile-time that no analytic solutions are
198  // available
199  template <typename Metavariables, typename ParallelComponent>
200  void operator()(
201  const typename ObservationValueTag::type& observation_value,
202  const Mesh<VolumeDim>& mesh,
203  const tnsr::I<DataVector, VolumeDim, Frame::Inertial>&
204  inertial_coordinates,
205  const typename NonSolutionTensors::type&... non_solution_tensors,
207  const ElementId<VolumeDim>& array_index,
208  const ParallelComponent* const component) const noexcept {
209  this->operator()(observation_value, mesh, inertial_coordinates,
210  non_solution_tensors..., std::nullopt, cache, array_index,
211  component);
212  }
213 
214  // We factor out the work into a static member function so it can be shared
215  // with other field observing events, like the one that deals with DG-subcell
216  // where there are two grids. This is to avoid copy-pasting all of the code.
217  template <typename OptionalAnalyticSolutions, typename Metavariables,
218  typename ParallelComponent>
219  static void call_operator_impl(
220  const std::string& subfile_path,
221  const std::unordered_set<std::string>& variables_to_observe,
222  const std::optional<Mesh<VolumeDim>>& interpolation_mesh,
223  const typename ObservationValueTag::type& observation_value,
224  const Mesh<VolumeDim>& mesh,
225  const tnsr::I<DataVector, VolumeDim, Frame::Inertial>&
226  inertial_coordinates,
227  const typename AnalyticSolutionTensors::
228  type&... analytic_solution_tensors,
229  const typename NonSolutionTensors::type&... non_solution_tensors,
230  const OptionalAnalyticSolutions& optional_analytic_solutions,
232  const ElementId<VolumeDim>& array_index,
233  const ParallelComponent* const /*meta*/) noexcept {
234  const auto analytic_solutions = [&optional_analytic_solutions]() {
235  if constexpr (tt::is_a_v<std::optional, OptionalAnalyticSolutions>) {
236  return optional_analytic_solutions.has_value()
237  ? std::make_optional(std::cref(*optional_analytic_solutions))
238  : std::nullopt;
239  } else {
240  return std::make_optional(std::cref(optional_analytic_solutions));
241  }
242  }();
243 
244  const std::string element_name =
245  MakeString{} << ElementId<VolumeDim>(array_index) << '/';
246 
247  // if no interpolation_mesh is provided, the interpolation is essentially
248  // ignored by the RegularGridInterpolant except for a single copy.
249  const intrp::RegularGrid interpolant(mesh,
250  interpolation_mesh.value_or(mesh));
251 
252  // Remove tensor types, only storing individual components.
253  std::vector<TensorComponent> components;
254  // This is larger than we need if we are only observing some
255  // tensors, but that's not a big deal and calculating the correct
256  // size is nontrivial.
257  components.reserve(alg::accumulate(
259  inertial_coordinates.size(),
260  (analytic_solutions.has_value() ? 2 : 1) *
261  AnalyticSolutionTensors::type::size()...,
262  NonSolutionTensors::type::size()...},
263  0_st));
264 
265  const auto record_tensor_components = [&components, &element_name,
266  &interpolant, &variables_to_observe](
267  const auto tensor_tag_v,
268  const auto& tensor) noexcept {
269  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
270  if (variables_to_observe.count(db::tag_name<tensor_tag>()) == 1) {
271  for (size_t i = 0; i < tensor.size(); ++i) {
272  components.emplace_back(element_name + db::tag_name<tensor_tag>() +
273  tensor.component_suffix(i),
274  interpolant.interpolate(tensor[i]));
275  }
276  }
277  };
278  record_tensor_components(tmpl::type_<coordinates_tag>{},
279  inertial_coordinates);
280  EXPAND_PACK_LEFT_TO_RIGHT(record_tensor_components(
281  tmpl::type_<AnalyticSolutionTensors>{}, analytic_solution_tensors));
282  EXPAND_PACK_LEFT_TO_RIGHT(record_tensor_components(
283  tmpl::type_<NonSolutionTensors>{}, non_solution_tensors));
284 
285  if (analytic_solutions.has_value()) {
286  const auto record_errors = [&components, &element_name,
287  &analytic_solutions, &interpolant,
288  &variables_to_observe](
289  const auto tensor_tag_v,
290  const auto& tensor) noexcept {
291  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
292  if (variables_to_observe.count(db::tag_name<tensor_tag>()) == 1) {
293  for (size_t i = 0; i < tensor.size(); ++i) {
294  DataVector error = tensor[i] - get<::Tags::Analytic<tensor_tag>>(
295  analytic_solutions->get())[i];
296  components.emplace_back(element_name + "Error(" +
297  db::tag_name<tensor_tag>() + ")" +
298  tensor.component_suffix(i),
299  interpolant.interpolate(error));
300  }
301  }
302  };
303  EXPAND_PACK_LEFT_TO_RIGHT(record_errors(
304  tmpl::type_<AnalyticSolutionTensors>{}, analytic_solution_tensors));
305 
306  (void)(record_errors); // Silence GCC warning about unused variable
307  }
308 
309  // Send data to volume observer
310  auto& local_observer =
311  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
312  cache)
313  .ckLocalBranch();
314  Parallel::simple_action<observers::Actions::ContributeVolumeData>(
315  local_observer,
316  observers::ObservationId(observation_value, subfile_path + ".vol"),
317  subfile_path,
321  std::move(components), interpolation_mesh.value_or(mesh).extents(),
322  interpolation_mesh.value_or(mesh).basis(),
323  interpolation_mesh.value_or(mesh).quadrature());
324  }
325 
326  using observation_registration_tags = tmpl::list<>;
328  get_observation_type_and_key_for_registration() const noexcept {
330  observers::ObservationKey(subfile_path_ + ".vol")};
331  }
332 
333  bool needs_evolved_variables() const noexcept override { return true; }
334 
335  // NOLINTNEXTLINE(google-runtime-references)
336  void pup(PUP::er& p) noexcept override {
337  Event::pup(p);
338  p | subfile_path_;
339  p | variables_to_observe_;
340  p | interpolation_mesh_;
341  }
342 
343  private:
344  std::string subfile_path_;
345  std::unordered_set<std::string> variables_to_observe_{};
346  std::optional<Mesh<VolumeDim>> interpolation_mesh_{};
347 };
348 
349 /// \cond
350 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
351  typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
352 PUP::able::PUP_ID
353  ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
354  tmpl::list<AnalyticSolutionTensors...>,
355  tmpl::list<NonSolutionTensors...>>::my_PUP_ID = 0; // NOLINT
356 /// \endcond
357 } // namespace Events
358 } // namespace dg
observers::ObservationId
A unique identifier for an observation representing the type of observation and the instance (e....
Definition: ObservationId.hpp:71
std::string
CharmPupable.hpp
EXPAND_PACK_LEFT_TO_RIGHT
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:565
utility
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
domain::Tags::Coordinates< VolumeDim, Frame::Inertial >
PARSE_ERROR
#define PARSE_ERROR(context, m)
Definition: Options.hpp:71
functional
unordered_set
Literals.hpp
std::pair
GlobalCache.hpp
Options.hpp
Tags::AnalyticSolutionsBase
Base tag for the analytic solution tensors. Retrieved values can be either Variables or std::optional...
Definition: Tags.hpp:106
Tags.hpp
vector
observers::ObservationKey
Used as a key in maps to keep track of how many elements have registered.
Definition: ObservationId.hpp:28
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
Options::Context
Definition: Options.hpp:41
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ApplyMassMatrix.hpp:13
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:51
ElementId.hpp
std::add_pointer_t
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
Event
Definition: Event.hpp:19
dg::Events::ObserveFields< VolumeDim, ObservationValueTag, tmpl::list< Tensors... >, tmpl::list< AnalyticSolutionTensors... >, tmpl::list< NonSolutionTensors... > >::needs_evolved_variables
bool needs_evolved_variables() const noexcept override
Whether the event uses anything depending on the evolved_variables. If this returns false,...
Definition: ObserveFields.hpp:333
cstddef
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
Options::Auto
A class indicating that a parsed value can be automatically computed instead of specified.
Definition: Auto.hpp:33
alg::count
decltype(auto) count(const Container &c, const T &value)
Convenience wrapper around std::count.
Definition: Algorithm.hpp:208
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
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:382
Time
Definition: Time.hpp:29
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
observers::TypeOfObservation::Volume
@ Volume
The sender will only perform volume observations.
Frame
Definition: IndexType.hpp:36
optional
observers::ArrayComponentId
An ID type that identifies both the parallel component and the index in the parallel component.
Definition: ArrayComponentId.hpp:27
StdHelpers.hpp
Parallel::ArrayIndex
The array index used for indexing Chare Arrays, mostly an implementation detail.
Definition: ArrayIndex.hpp:28
intrp::RegularGrid
Interpolate data from a Mesh onto a regular grid of points.
Definition: RegularGridInterpolant.hpp:45
Prefixes.hpp
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
type_traits
TMPL.hpp
alg::accumulate
decltype(auto) accumulate(const Container &c, T init)
Convenience wrapper around std::accumulate, returns std::accumulate(begin(c), end(c),...
Definition: Numeric.hpp:54
MakeString
Make a string by streaming into object.
Definition: MakeString.hpp:18
initializer_list
string