ObserveFields.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <initializer_list>
8 #include <pup.h>
9 #include <string>
10 #include <type_traits>
11 #include <unordered_set>
12 #include <utility>
13 #include <vector>
14 
16 #include "DataStructures/DataVector.hpp"
17 #include "DataStructures/Tensor/TensorData.hpp"
18 #include "Domain/ElementId.hpp"
19 #include "Domain/ElementIndex.hpp" // IWYU pragma: keep
20 #include "Domain/Tags.hpp"
21 #include "Evolution/EventsAndTriggers/Event.hpp"
22 #include "IO/Observer/ArrayComponentId.hpp"
23 #include "IO/Observer/ObservationId.hpp"
24 #include "IO/Observer/ObserverComponent.hpp" // IWYU pragma: keep
25 #include "IO/Observer/VolumeActions.hpp" // IWYU pragma: keep
26 #include "Options/Options.hpp"
27 #include "Parallel/ArrayIndex.hpp"
30 #include "Parallel/Invoke.hpp"
31 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
32 #include "Time/Time.hpp"
33 #include "Utilities/Algorithm.hpp"
34 #include "Utilities/Literals.hpp"
35 #include "Utilities/MakeString.hpp"
36 #include "Utilities/Numeric.hpp"
37 #include "Utilities/TMPL.hpp"
38 
39 /// \cond
40 template <size_t Dim>
41 class Mesh;
42 namespace Frame {
43 struct Inertial;
44 } // namespace Frame
45 namespace Tags {
46 struct Time;
47 } // namespace Tags
48 /// \endcond
49 
50 namespace dg {
51 namespace Events {
52 template <size_t VolumeDim, typename Tensors, typename AnalyticSolutionTensors,
53  typename EventRegistrars,
54  typename NonSolutionTensors =
55  tmpl::list_difference<Tensors, AnalyticSolutionTensors>>
57 
58 namespace Registrars {
59 template <size_t VolumeDim, typename Tensors,
60  typename AnalyticSolutionTensors = tmpl::list<>>
61 struct ObserveFields {
62  template <typename RegistrarList>
63  using f = Events::ObserveFields<VolumeDim, Tensors, AnalyticSolutionTensors,
64  RegistrarList>;
65 };
66 } // namespace Registrars
67 
68 template <size_t VolumeDim, typename Tensors,
69  typename AnalyticSolutionTensors = tmpl::list<>,
70  typename EventRegistrars = tmpl::list<Registrars::ObserveFields<
71  VolumeDim, Tensors, AnalyticSolutionTensors>>,
72  typename NonSolutionTensors>
73 class ObserveFields; // IWYU pragma: keep
74 
75 /*!
76  * \ingroup DiscontinuousGalerkinGroup
77  * \brief %Observe volume tensor fields.
78  *
79  * Writes volume quantities:
80  * - `InertialCoordinates`
81  * - Tensors listed in `Tensors` template parameter
82  * - `Error(*)` = errors in `AnalyticSolutionTensors` =
83  * \f$\text{value} - \text{analytic solution}\f$
84  *
85  * \warning Currently, only one volume observation event can be
86  * triggered at a given time. Causing multiple events to run at once
87  * will produce unpredictable results.
88  */
89 template <size_t VolumeDim, typename... Tensors,
90  typename... AnalyticSolutionTensors, typename EventRegistrars,
91  typename... NonSolutionTensors>
92 class ObserveFields<VolumeDim, tmpl::list<Tensors...>,
93  tmpl::list<AnalyticSolutionTensors...>, EventRegistrars,
94  tmpl::list<NonSolutionTensors...>>
95  : public Event<EventRegistrars> {
96  private:
97  static_assert(
99  tmpl::list_difference<tmpl::list<AnalyticSolutionTensors...>,
100  tmpl::list<Tensors...>>,
101  tmpl::list<>>,
102  "All AnalyticSolutionTensors must be listed in Tensors.");
104 
105  template <typename T>
106  static std::string component_suffix(const T& tensor,
107  size_t component_index) noexcept {
108  return tensor.rank() == 0
109  ? ""
110  : "_" + tensor.component_name(
111  tensor.get_tensor_index(component_index));
112  }
113 
114  public:
115  /// \cond
116  explicit ObserveFields(CkMigrateMessage* /*unused*/) noexcept {}
117  using PUP::able::register_constructor;
119  /// \endcond
120 
121  struct VariablesToObserve {
122  static constexpr OptionString help = "Subset of variables to observe";
124  static type default_value() noexcept { return {Tensors::name()...}; }
125  static size_t lower_bound_on_size() noexcept { return 1; }
126  };
127 
128  using options = tmpl::list<VariablesToObserve>;
129  static constexpr OptionString 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  explicit ObserveFields(const std::vector<std::string>& variables_to_observe =
143  VariablesToObserve::default_value(),
144  const OptionContext& context = {})
145  : variables_to_observe_(variables_to_observe.begin(),
146  variables_to_observe.end()) {
147  const std::unordered_set<std::string> valid_tensors{Tensors::name()...};
148  for (const auto& name : variables_to_observe_) {
149  if (valid_tensors.count(name) != 1) {
150  PARSE_ERROR(
151  context,
152  name << " is not an available variable. Available variables:\n"
153  << (std::vector<std::string>{Tensors::name()...}));
154  }
155  if (alg::count(variables_to_observe, name) != 1) {
156  PARSE_ERROR(context, name << " specified multiple times");
157  }
158  }
159  variables_to_observe_.insert(coordinates_tag::name());
160  }
161 
162  using argument_tags =
163  tmpl::list<::Tags::Time, ::Tags::Mesh<VolumeDim>, coordinates_tag,
164  AnalyticSolutionTensors..., NonSolutionTensors...>;
165 
166  template <typename Metavariables, typename ParallelComponent>
167  void operator()(
168  const Time& time, const Mesh<VolumeDim>& mesh,
169  const db::item_type<coordinates_tag>& inertial_coordinates,
170  const db::item_type<
171  AnalyticSolutionTensors>&... analytic_solution_tensors,
172  const db::item_type<NonSolutionTensors>&... non_solution_tensors,
174  const ElementIndex<VolumeDim>& array_index,
175  const ParallelComponent* const /*meta*/) const noexcept {
176  const std::string element_name =
177  MakeString{} << ElementId<VolumeDim>(array_index) << '/';
178 
179  // Remove tensor types, only storing individual components.
180  std::vector<TensorComponent> components;
181  // This is larger than we need if we are only observing some
182  // tensors, but that's not a big deal and calculating the correct
183  // size is nontrivial.
184  components.reserve(alg::accumulate(
186  inertial_coordinates.size(),
189  0_st));
190 
191  const auto record_tensor_components = [this, &components, &element_name](
192  const auto tensor_tag_v, const auto& tensor) noexcept {
193  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
194  if (variables_to_observe_.count(tensor_tag::name()) == 1) {
195  for (size_t i = 0; i < tensor.size(); ++i) {
196  components.emplace_back(
197  element_name + tensor_tag::name() + component_suffix(tensor, i),
198  tensor[i]);
199  }
200  }
201  };
202  record_tensor_components(tmpl::type_<coordinates_tag>{},
203  inertial_coordinates);
204  EXPAND_PACK_LEFT_TO_RIGHT(record_tensor_components(
205  tmpl::type_<AnalyticSolutionTensors>{}, analytic_solution_tensors));
206  EXPAND_PACK_LEFT_TO_RIGHT(record_tensor_components(
207  tmpl::type_<NonSolutionTensors>{}, non_solution_tensors));
208 
209  const auto record_errors =
210  [ this, &inertial_coordinates, &time, &components, &
211  element_name ](const auto tensor_tag_v, const auto& tensor,
212  const auto& local_cache) noexcept {
213  const auto analytic_solution =
214  Parallel::get<OptionTags::AnalyticSolutionBase>(local_cache)
215  .variables(inertial_coordinates, time.value(),
216  tmpl::list<AnalyticSolutionTensors...>{});
217  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
218  if (variables_to_observe_.count(tensor_tag::name()) == 1) {
219  for (size_t i = 0; i < tensor.size(); ++i) {
220  DataVector error = tensor[i] - get<tensor_tag>(analytic_solution)[i];
221  components.emplace_back(element_name + "Error(" + tensor_tag::name() +
222  ")" + component_suffix(tensor, i),
223  std::move(error));
224  }
225  }
226  };
228  record_errors(tmpl::type_<AnalyticSolutionTensors>{},
229  analytic_solution_tensors, cache));
230 
231  (void)(record_errors); // Silence GCC warning about unused variable
232 
233  // Send data to volume observer
234  auto& local_observer =
235  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
236  cache)
237  .ckLocalBranch();
238  Parallel::simple_action<observers::Actions::ContributeVolumeData>(
239  local_observer,
241  time.value(), typename Metavariables::element_observation_type{}),
242  std::string{"/element_data"},
246  std::move(components), mesh.extents());
247  }
248 
249  // NOLINTNEXTLINE(google-runtime-references)
250  void pup(PUP::er& p) noexcept override {
252  p | variables_to_observe_;
253  }
254 
255  private:
256  std::unordered_set<std::string> variables_to_observe_{};
257 };
258 
259 /// \cond
260 template <size_t VolumeDim, typename... Tensors,
261  typename... AnalyticSolutionTensors, typename EventRegistrars,
262  typename... NonSolutionTensors>
263 PUP::able::PUP_ID
264  ObserveFields<VolumeDim, tmpl::list<Tensors...>,
265  tmpl::list<AnalyticSolutionTensors...>, EventRegistrars,
266  tmpl::list<NonSolutionTensors...>>::my_PUP_ID = 0; // NOLINT
267 /// \endcond
268 } // namespace Events
269 } // namespace dg
Definition: SolvePoissonProblem.hpp:67
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:561
Definition: Digraph.hpp:11
The time in a simulation. Times can be safely compared for exact equality as long as they do not belo...
Definition: Time.hpp:31
Definition: Filtering.hpp:38
An ElementId uniquely labels an Element. It is constructed from the BlockId of the Block to which the...
Definition: ElementId.hpp:36
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
Defines classes and functions for making classes creatable from input files.
Defines macros to allow serialization of abstract template base classes.
Definition: Completion.hpp:13
Defines class ElementId.
An ID type that identifies both the parallel component and the index in the parallel component...
Definition: ArrayComponentId.hpp:27
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
Definition: ObserveFields.hpp:56
Defines Time and TimeDelta.
decltype(auto) count(const Container &c, const T &value)
Convenience wrapper around std::count.
Definition: Algorithm.hpp:166
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
#define PARSE_ERROR(context, m)
Like ERROR("\n" << (context) << m), but instead throws an exception that will be caught in a higher l...
Definition: Options.hpp:67
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
constexpr bool is_same_v
Variable template for is_same.
Definition: TypeTraits.hpp:221
Defines useful literals.
Indicates the Frame that a TensorIndexType is in.
Definition: IndexType.hpp:36
Definition: DataBoxTag.hpp:29
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
const Index< Dim > & extents() const noexcept
The number of grid points in each dimension of the grid.
Definition: Mesh.hpp:107
Stores a collection of function values.
Definition: DataVector.hpp:42
Defines class ElementIndex.
Information about the nested operations being performed by the parser, for use in printing errors...
Definition: Options.hpp:36
Wraps the template metaprogramming library used (brigand)
Base class for something that can happen during a simulation (such as an observation).
Definition: Event.hpp:30
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that is returned by the Tag. If it is a base tag then a TagList must be passed as a seco...
Definition: DataBoxTag.hpp:410
Defines tags related to domain quantities.
Definition: ObserveFields.hpp:61
A class for indexing a Charm array by Element.
Definition: ElementIndex.hpp:53
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
A type-erased identifier that combines the identifier&#39;s type and hash used to uniquely identify an ob...
Definition: ObservationId.hpp:42
Defines class template ConstGlobalCache.
The array index used for indexing Chare Arrays, mostly an implementation detail.
Definition: ArrayIndex.hpp:20