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/CoordinateMaps/Tags.hpp"
23 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
24 #include "Domain/FunctionsOfTime/Tags.hpp"
26 #include "Domain/Tags.hpp"
27 #include "Evolution/DgSubcell/ActiveGrid.hpp"
28 #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
29 #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
30 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
31 #include "Evolution/TypeTraits.hpp"
32 #include "IO/Observer/ArrayComponentId.hpp"
33 #include "IO/Observer/ObservationId.hpp"
34 #include "IO/Observer/ObserverComponent.hpp" // IWYU pragma: keep
35 #include "IO/Observer/VolumeActions.hpp" // IWYU pragma: keep
36 #include "Options/Auto.hpp"
37 #include "Options/Options.hpp"
38 #include "Parallel/ArrayIndex.hpp"
40 #include "Parallel/GlobalCache.hpp"
41 #include "Parallel/Invoke.hpp"
42 #include "Parallel/PupStlCpp17.hpp"
43 #include "ParallelAlgorithms/Events/ObserveFields.hpp"
44 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
45 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
46 #include "Utilities/Algorithm.hpp"
47 #include "Utilities/Literals.hpp"
48 #include "Utilities/MakeString.hpp"
49 #include "Utilities/Numeric.hpp"
50 #include "Utilities/StdHelpers.hpp"
51 #include "Utilities/TMPL.hpp"
52 #include "Utilities/TypeTraits/IsA.hpp"
53 
54 /// \cond
55 template <size_t Dim>
56 class Mesh;
57 namespace Frame {
58 struct Inertial;
59 } // namespace Frame
60 /// \endcond
61 
62 namespace evolution::dg::subcell::Events {
63 /// \cond
64 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
65  typename AnalyticSolutionTensors = tmpl::list<>,
66  typename NonSolutionTensors =
68 class ObserveFields;
69 /// \cond
70 
71 /*!
72  * \ingroup DiscontinuousGalerkinGroup
73  * \brief %Observe volume tensor fields.
74  *
75  * A class that writes volume quantities to an h5 file during the simulation.
76  * The observed quantitites are:
77  * - `InertialCoordinates`
78  * - Tensors listed in `Tensors` template parameter
79  * - `Error(*)` = errors in `AnalyticSolutionTensors` =
80  * \f$\text{value} - \text{analytic solution}\f$
81  *
82  * The user may specify an `interpolation_mesh` to which the
83  * data is interpolated.
84  */
85 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
86  typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
87 class ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
88  tmpl::list<AnalyticSolutionTensors...>,
89  tmpl::list<NonSolutionTensors...>> : public Event {
90  private:
91  static_assert(
92  std::is_same_v<
93  tmpl::list_difference<tmpl::list<AnalyticSolutionTensors...>,
94  tmpl::list<Tensors...>>,
95  tmpl::list<>>,
96  "All AnalyticSolutionTensors must be listed in Tensors.");
97 
98  public:
99  /// The name of the subfile inside the HDF5 file
100  struct SubfileName {
101  using type = std::string;
102  static constexpr Options::String help = {
103  "The name of the subfile inside the HDF5 file without an extension and "
104  "without a preceding '/'."};
105  };
106 
107  /// \cond
108  explicit ObserveFields(CkMigrateMessage* /*unused*/) noexcept {}
109  using PUP::able::register_constructor;
110  WRAPPED_PUPable_decl_template(ObserveFields); // NOLINT
111  /// \endcond
112 
114  static constexpr Options::String help = "Subset of variables to observe";
116  static size_t lower_bound_on_size() noexcept { return 1; }
117  };
118 
120  using type = Options::Auto<Mesh<VolumeDim>, Options::AutoLabel::None>;
121  static constexpr Options::String help =
122  "An optional mesh to which the variables are interpolated. This mesh "
123  "specifies any number of collocation points, basis, and quadrature on "
124  "which the observed quantities are evaluated. If no mesh is given, the "
125  "results will be evaluated on the mesh the simulation runs on. The "
126  "user may add several ObserveField Events e.g. with and without an "
127  "interpolating mesh to output the data both on the original mesh and "
128  "on a new mesh.";
129  };
130 
131  using options =
132  tmpl::list<SubfileName, VariablesToObserve, InterpolateToMesh>;
133  static constexpr Options::String help =
134  "Observe volume tensor fields.\n"
135  "\n"
136  "Writes volume quantities:\n"
137  " * InertialCoordinates\n"
138  " * Tensors listed in Tensors template parameter\n"
139  " * Error(*) = errors in AnalyticSolutionTensors\n"
140  " = value - analytic solution\n"
141  "\n"
142  "Warning: Currently, only one volume observation event can be\n"
143  "triggered at a given time. Causing multiple events to run at once\n"
144  "will produce unpredictable results.";
145 
146  ObserveFields() = default;
147 
148  explicit ObserveFields(const std::string& subfile_name,
149  const std::vector<std::string>& variables_to_observe,
150  std::optional<Mesh<VolumeDim>> interpolation_mesh = {},
151  const Options::Context& context = {});
152 
153  using coordinates_tag =
155 
156  using argument_tags =
157  tmpl::list<ObservationValueTag, ::domain::Tags::Mesh<VolumeDim>,
162  VolumeDim, Frame::Grid, Frame::Inertial>,
163  ::domain::Tags::FunctionsOfTime, AnalyticSolutionTensors...,
164  NonSolutionTensors...>;
165 
166  template <typename Metavariables, typename ParallelComponent>
167  void operator()(
168  const double observation_value, const Mesh<VolumeDim>& dg_mesh,
169  const Mesh<VolumeDim>& subcell_mesh,
170  const subcell::ActiveGrid active_grid,
171  const tnsr::I<DataVector, VolumeDim, Frame::Grid>& dg_grid_coordinates,
172  const tnsr::I<DataVector, VolumeDim, Frame::Grid>&
173  subcell_grid_coordinates,
174  const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial,
175  VolumeDim>& grid_to_inertial_map,
176  const std::unordered_map<
177  std::string,
179  functions_of_time,
180  const typename AnalyticSolutionTensors::
181  type&... analytic_solution_tensors,
182  const typename NonSolutionTensors::type&... non_solution_tensors,
184  const ElementId<VolumeDim>& array_index,
185  const ParallelComponent* const component) const noexcept {
187  Variables<tmpl::list<::Tags::Analytic<AnalyticSolutionTensors>...>>>
188  analytic_solution_variables{};
189  const auto set_analytic_soln =
190  [&analytic_solution_variables, &cache, &observation_value](
191  const Mesh<VolumeDim>& mesh,
192  const tnsr::I<DataVector, VolumeDim, Frame::Inertial>&
193  inertial_coords) noexcept {
194  if constexpr (evolution::is_analytic_solution_v<
195  typename Metavariables::initial_data>) {
196  Variables<tmpl::list<AnalyticSolutionTensors...>> soln_vars{
197  mesh.number_of_grid_points()};
198  soln_vars.assign_subset(
199  Parallel::get<::Tags::AnalyticSolutionBase>(cache).variables(
200  inertial_coords, observation_value,
201  tmpl::list<AnalyticSolutionTensors...>{}));
202  analytic_solution_variables = std::move(soln_vars);
203  } else {
204  (void)analytic_solution_variables;
205  (void)cache;
206  (void)observation_value;
207  }
208  };
209  if (active_grid == subcell::ActiveGrid::Dg) {
210  const auto dg_inertial_coords = grid_to_inertial_map(
211  dg_grid_coordinates, observation_value, functions_of_time);
212  set_analytic_soln(dg_mesh, dg_inertial_coords);
213  ::dg::Events::ObserveFields<VolumeDim, ObservationValueTag,
214  tmpl::list<Tensors...>,
215  tmpl::list<AnalyticSolutionTensors...>,
216  tmpl::list<NonSolutionTensors...>>::
217  call_operator_impl(
218  subfile_path_, variables_to_observe_, interpolation_mesh_,
219  observation_value, dg_mesh, dg_inertial_coords,
220  analytic_solution_tensors..., non_solution_tensors...,
221  analytic_solution_variables, cache, array_index, component);
222  } else {
223  ASSERT(active_grid == subcell::ActiveGrid::Subcell,
224  "Active grid must be either Dg or Subcell");
225  const auto subcell_inertial_coords = grid_to_inertial_map(
226  subcell_grid_coordinates, observation_value, functions_of_time);
227  set_analytic_soln(subcell_mesh, subcell_inertial_coords);
228  ::dg::Events::ObserveFields<VolumeDim, ObservationValueTag,
229  tmpl::list<Tensors...>,
230  tmpl::list<AnalyticSolutionTensors...>,
231  tmpl::list<NonSolutionTensors...>>::
232  call_operator_impl(
233  subfile_path_, variables_to_observe_, interpolation_mesh_,
234  observation_value, subcell_mesh, subcell_inertial_coords,
235  analytic_solution_tensors..., non_solution_tensors...,
236  analytic_solution_variables, cache, array_index, component);
237  }
238  }
239 
240  // This overload is called when the list of analytic-solution tensors is
241  // empty, i.e. it is clear at compile-time that no analytic solutions are
242  // available
243  template <typename Metavariables, typename ParallelComponent>
244  void operator()(
245  const double observation_value, const Mesh<VolumeDim>& mesh,
246  const tnsr::I<DataVector, VolumeDim, Frame::Inertial>&
247  inertial_coordinates,
248  const typename NonSolutionTensors::type&... non_solution_tensors,
250  const ElementId<VolumeDim>& array_index,
251  const ParallelComponent* const component) const noexcept {
252  this->operator()(observation_value, mesh, inertial_coordinates,
253  non_solution_tensors..., std::nullopt, cache, array_index,
254  component);
255  }
256 
257  using observation_registration_tags = tmpl::list<>;
259  get_observation_type_and_key_for_registration() const noexcept {
261  observers::ObservationKey(subfile_path_ + ".vol")};
262  }
263 
264  // NOLINTNEXTLINE(google-runtime-references)
265  void pup(PUP::er& p) noexcept override {
266  Event::pup(p);
267  p | subfile_path_;
268  p | variables_to_observe_;
269  p | interpolation_mesh_;
270  }
271 
272  bool needs_evolved_variables() const noexcept override { return true; }
273 
274  private:
275  std::string subfile_path_;
276  std::unordered_set<std::string> variables_to_observe_{};
277  std::optional<Mesh<VolumeDim>> interpolation_mesh_{};
278 };
279 
280 /// \cond
281 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
282  typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
283 ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
284  tmpl::list<AnalyticSolutionTensors...>,
285  tmpl::list<NonSolutionTensors...>>::
286  ObserveFields(const std::string& subfile_name,
287  const std::vector<std::string>& variables_to_observe,
288  std::optional<Mesh<VolumeDim>> interpolation_mesh,
289  const Options::Context& context)
290  : subfile_path_("/" + subfile_name),
291  variables_to_observe_(variables_to_observe.begin(),
292  variables_to_observe.end()),
293  interpolation_mesh_(interpolation_mesh) {
294  using ::operator<<;
295  const std::unordered_set<std::string> valid_tensors{
296  db::tag_name<Tensors>()...};
297  for (const auto& name : variables_to_observe_) {
298  if (valid_tensors.count(name) != 1) {
299  PARSE_ERROR(
300  context,
301  name << " is not an available variable. Available variables:\n"
302  << (std::vector<std::string>{db::tag_name<Tensors>()...}));
303  }
304  if (alg::count(variables_to_observe, name) != 1) {
305  PARSE_ERROR(context, name << " specified multiple times");
306  }
307  }
308  variables_to_observe_.insert(coordinates_tag::name());
309  if (interpolation_mesh.has_value()) {
310  PARSE_ERROR(
311  context,
312  "We don't yet support interpolating to a mesh with DG-subcell because "
313  "we don't have an interpolator on the subcell grid. Trilinear "
314  "interpolation could be added to support interpolation.");
315  }
316 }
317 /// \endcond
318 
319 /// \cond
320 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
321  typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
322 PUP::able::PUP_ID
323  ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
324  tmpl::list<AnalyticSolutionTensors...>,
325  tmpl::list<NonSolutionTensors...>>::my_PUP_ID = 0; // NOLINT
326 /// \endcond
327 } // namespace evolution::dg::subcell::Events
std::string
CharmPupable.hpp
utility
Frame::Inertial
Definition: IndexType.hpp:44
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
domain::Tags::FunctionsOfTime
The FunctionsOfTime initialized from a DomainCreator or (if override_functions_of_time is true in the...
Definition: Tags.hpp:60
domain::Tags::Coordinates< VolumeDim, Frame::Inertial >
PARSE_ERROR
#define PARSE_ERROR(context, m)
Definition: Options.hpp:71
functional
evolution::dg::subcell::Events::InterpolateToMesh
Definition: ObserveFields.hpp:119
Frame::Grid
Definition: IndexType.hpp:43
unordered_set
Literals.hpp
std::pair
GlobalCache.hpp
Options.hpp
Tags.hpp
vector
observers::ObservationKey
Used as a key in maps to keep track of how many elements have registered.
Definition: ObservationId.hpp:28
CoordinateMap.hpp
Options::Context
Definition: Options.hpp:41
evolution::dg::subcell::ActiveGrid
ActiveGrid
Definition: ActiveGrid.hpp:11
evolution::dg::subcell::Events::VariablesToObserve
Definition: ObserveFields.hpp:113
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:51
ElementId.hpp
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
Event
Definition: Event.hpp:19
cstddef
Options::Auto
A class indicating that a parsed value can be automatically computed instead of specified.
Definition: Auto.hpp:33
evolution::dg::subcell::Tags::Mesh
The mesh on the subcells.
Definition: Mesh.hpp:18
alg::count
decltype(auto) count(const Container &c, const T &value)
Convenience wrapper around std::count.
Definition: Algorithm.hpp:208
domain::CoordinateMaps::Tags::CoordinateMap
Definition: Tags.hpp:30
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
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
std::begin
T begin(T... args)
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
std::end
T end(T... args)
StdHelpers.hpp
evolution::dg::subcell::Tags::Coordinates
The coordinates in a given frame.
Definition: Coordinates.hpp:27
std::unique_ptr
Prefixes.hpp
std::unordered_map
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
initializer_list
string