ObserveVolumeIntegrals.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <pup.h>
8 #include <string>
9 #include <utility>
10 #include <vector>
11 
13 #include "DataStructures/DataVector.hpp"
14 #include "Domain/Tags.hpp"
15 #include "IO/Observer/ArrayComponentId.hpp"
16 #include "IO/Observer/Helpers.hpp"
17 #include "IO/Observer/ObservationId.hpp"
18 #include "IO/Observer/ObserverComponent.hpp"
19 #include "IO/Observer/ReductionActions.hpp"
20 #include "IO/Observer/TypeOfObservation.hpp"
22 #include "Options/Options.hpp"
23 #include "Parallel/ArrayIndex.hpp"
25 #include "Parallel/GlobalCache.hpp"
26 #include "Parallel/Invoke.hpp"
27 #include "Parallel/Reduction.hpp"
28 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
30 #include "Utilities/Functional.hpp"
31 #include "Utilities/Registration.hpp"
32 #include "Utilities/TMPL.hpp"
33 #include "Utilities/TaggedTuple.hpp"
34 
35 /// \cond
36 namespace Frame {
37 struct Inertial;
38 } // namespace Frame
39 /// \endcond
40 
41 namespace dg {
42 namespace Events {
43 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
44  typename EventRegistrars>
46 
47 namespace Registrars {
48 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors>
49 // Presence of size_t template argument requires to define this struct
50 // instead of using Registration::Registrar alias.
52  template <typename RegistrarList>
53  using f = Events::ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
54  Tensors, RegistrarList>;
55 };
56 } // namespace Registrars
57 
58 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
59  typename EventRegistrars =
61  VolumeDim, ObservationValueTag, Tensors>>>
63 
64 /*!
65  * \ingroup DiscontinuousGalerkinGroup
66  * \brief %Observe the volume integrals of the tensors over the domain.
67  *
68  * Writes reduction quantities:
69  * - `ObservationValueTag`
70  * - `Volume` = volume of the domain
71  * - `VolumeIntegral(*)` = volume integral of the tensor
72  */
73 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
74  typename EventRegistrars>
75 class ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
76  tmpl::list<Tensors...>, EventRegistrars>
77  : public Event<EventRegistrars> {
78  private:
79  using VolumeIntegralDatum =
81 
82  using ReductionData = tmpl::wrap<
83  tmpl::list<Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
86  Parallel::ReductionData>;
87 
88  public:
89  /// The name of the subfile inside the HDF5 file
90  struct SubfileName {
91  using type = std::string;
92  static constexpr Options::String help = {
93  "The name of the subfile inside the HDF5 file without an extension and "
94  "without a preceding '/'."};
95  };
96 
97  /// \cond
98  explicit ObserveVolumeIntegrals(CkMigrateMessage* /*unused*/) noexcept {}
99  using PUP::able::register_constructor;
101  /// \endcond
102 
103  using options = tmpl::list<SubfileName>;
104  static constexpr Options::String help =
105  "Observe the volume integrals of the tensors over the domain.\n"
106  "\n"
107  "Writes reduction quantities:\n"
108  " * ObservationValueTag\n"
109  " * Volume = volume of the domain\n"
110  " * VolumeIntegral(*) = volume integral of the tensor\n"
111  "\n"
112  "Warning: Currently, only one reduction observation event can be\n"
113  "triggered at a given observation value. Causing multiple events to\n"
114  "run at once will produce unpredictable results.";
115 
116  ObserveVolumeIntegrals() = default;
117  explicit ObserveVolumeIntegrals(const std::string& subfile_name) noexcept;
118 
119  using observed_reduction_data_tags =
120  observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
121 
122  using argument_tags =
123  tmpl::list<ObservationValueTag, domain::Tags::Mesh<VolumeDim>,
125  Tensors...>;
126 
127  template <typename Metavariables, typename ArrayIndex,
128  typename ParallelComponent>
129  void operator()(const typename ObservationValueTag::type& observation_value,
130  const Mesh<VolumeDim>& mesh,
131  const Scalar<DataVector>& det_inv_jacobian,
132  const typename Tensors::type&... tensors,
134  const ArrayIndex& array_index,
135  const ParallelComponent* const /*meta*/) const noexcept {
136  // Determinant of Jacobian is needed because integral is performed in
137  // logical coords.
138  const DataVector det_jacobian = 1.0 / get(det_inv_jacobian);
139  const double local_volume = definite_integral(det_jacobian, mesh);
140 
141  std::vector<double> local_volume_integrals{};
142  std::vector<std::string> reduction_names = {
143  db::tag_name<ObservationValueTag>(), "Volume"};
144  const auto record_integrals =
145  [&local_volume_integrals, &reduction_names, &det_jacobian, &mesh ](
146  const auto tensor_tag_v, const auto& tensor) noexcept {
147  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
148  for (size_t i = 0; i < tensor.size(); ++i) {
149  reduction_names.push_back("VolumeIntegral(" +
150  db::tag_name<tensor_tag>() +
151  tensor.component_suffix(i) + ")");
152  local_volume_integrals.push_back(
153  definite_integral(det_jacobian * tensor[i], mesh));
154  }
155  return 0;
156  };
158  record_integrals(tmpl::type_<Tensors>{}, tensors));
159 
160  // Send data to reduction observer
161  auto& local_observer =
162  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
163  cache)
164  .ckLocalBranch();
165  Parallel::simple_action<observers::Actions::ContributeReductionData>(
166  local_observer,
167  observers::ObservationId(observation_value, subfile_path_ + ".dat"),
170  Parallel::ArrayIndex<ArrayIndex>(array_index)},
171  subfile_path_, reduction_names,
172  ReductionData{static_cast<double>(observation_value), local_volume,
173  local_volume_integrals});
174  }
175 
176  using observation_registration_tags = tmpl::list<>;
178  get_observation_type_and_key_for_registration() const noexcept {
180  observers::ObservationKey(subfile_path_ + ".dat")};
181  }
182 
183  // NOLINTNEXTLINE(google-runtime-references)
184  void pup(PUP::er& p) override {
186  p | subfile_path_;
187  }
188 
189  private:
190  std::string subfile_path_;
191 };
192 
193 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
194  typename EventRegistrars>
195 ObserveVolumeIntegrals<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
196  EventRegistrars>::
197  ObserveVolumeIntegrals(const std::string& subfile_name) noexcept
198  : subfile_path_("/" + subfile_name) {}
199 
200 /// \cond
201 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
202  typename EventRegistrars>
203 PUP::able::PUP_ID ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
204  tmpl::list<Tensors...>,
205  EventRegistrars>::my_PUP_ID =
206  0; // NOLINT
207 /// \endcond
208 } // namespace Events
209 } // 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:563
utility
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:639
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
funcl::VectorPlus
Function for adding two std::vectors of double component-wise.
Definition: Functional.hpp:312
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:63
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
dg::Events::ObserveVolumeIntegrals
Definition: ObserveVolumeIntegrals.hpp:45
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ComputeNonconservativeBoundaryFluxes.hpp:23
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:30
cstddef
ConstantExpressions.hpp
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
definite_integral
double definite_integral(const DataVector &integrand, const Mesh< Dim > &mesh) noexcept
Compute the definite integral of a function over a manifold.
Mesh< VolumeDim >
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
dg::Events::Registrars::ObserveVolumeIntegrals
Definition: ObserveVolumeIntegrals.hpp:51
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Frame
Definition: IndexType.hpp:36
observers::TypeOfObservation::Reduction
@ Reduction
The sender will only perform reduction observations.
DefiniteIntegral.hpp
observers::ArrayComponentId
An ID type that identifies both the parallel component and the index in the parallel component.
Definition: ArrayComponentId.hpp:27
Parallel::ArrayIndex
The array index used for indexing Chare Arrays, mostly an implementation detail.
Definition: ArrayIndex.hpp:27
domain::Tags::DetInvJacobian
The determinant of the inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:200
Prefixes.hpp
TMPL.hpp
string