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/Helpers.hpp"
16 #include "IO/Observer/ObservationId.hpp"
17 #include "IO/Observer/ObserverComponent.hpp"
18 #include "IO/Observer/ReductionActions.hpp"
20 #include "Options/Options.hpp"
23 #include "Parallel/Invoke.hpp"
24 #include "Parallel/Reduction.hpp"
25 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
27 #include "Utilities/Functional.hpp"
28 #include "Utilities/TMPL.hpp"
29 #include "Utilities/TaggedTuple.hpp"
30 
31 /// \cond
32 namespace Frame {
33 struct Inertial;
34 } // namespace Frame
35 /// \endcond
36 
37 namespace dg {
38 namespace Events {
39 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
40  typename EventRegistrars>
42 
43 namespace Registrars {
44 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors>
45 // Presence of size_t template argument requires to define this struct
46 // instead of using Registration::Registrar alias.
48  template <typename RegistrarList>
49  using f = Events::ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
50  Tensors, RegistrarList>;
51 };
52 } // namespace Registrars
53 
54 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
55  typename EventRegistrars =
57  VolumeDim, ObservationValueTag, Tensors>>>
59 
60 /*!
61  * \ingroup DiscontinuousGalerkinGroup
62  * \brief %Observe the volume integrals of the tensors over the domain.
63  *
64  * Writes reduction quantities:
65  * - `ObservationValueTag`
66  * - `Volume` = volume of the domain
67  * - `VolumeIntegral(*)` = volume integral of the tensor
68  *
69  * \warning Currently, only one reduction observation event can be
70  * triggered at a given observation value. Causing multiple events to run at
71  * once will produce unpredictable results.
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  /// \cond
90  explicit ObserveVolumeIntegrals(CkMigrateMessage* /*unused*/) noexcept {}
91  using PUP::able::register_constructor;
93  /// \endcond
94 
95  using options = tmpl::list<>;
96  static constexpr OptionString help =
97  "Observe the volume integrals of the tensors over the domain.\n"
98  "\n"
99  "Writes reduction quantities:\n"
100  " * ObservationValueTag\n"
101  " * Volume = volume of the domain\n"
102  " * VolumeIntegral(*) = volume integral of the tensor\n"
103  "\n"
104  "Warning: Currently, only one reduction observation event can be\n"
105  "triggered at a given observation value. Causing multiple events to\n"
106  "run at once will produce unpredictable results.";
107 
108  ObserveVolumeIntegrals() = default;
109 
110  using observed_reduction_data_tags =
111  observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
112 
113  using argument_tags =
114  tmpl::list<ObservationValueTag, domain::Tags::Mesh<VolumeDim>,
116  Tensors...>;
117 
118  template <typename Metavariables, typename ArrayIndex,
119  typename ParallelComponent>
120  void operator()(const typename ObservationValueTag::type& observation_value,
121  const Mesh<VolumeDim>& mesh,
122  const Scalar<DataVector>& det_inv_jacobian,
123  const typename Tensors::type&... tensors,
125  const ArrayIndex& /*array_index*/,
126  const ParallelComponent* const /*meta*/) const noexcept {
127  // Determinant of Jacobian is needed because integral is performed in
128  // logical coords.
129  const DataVector det_jacobian = 1.0 / get(det_inv_jacobian);
130  const double local_volume = definite_integral(det_jacobian, mesh);
131 
132  std::vector<double> local_volume_integrals{};
133  std::vector<std::string> reduction_names = {
134  db::tag_name<ObservationValueTag>(), "Volume"};
135  const auto record_integrals =
136  [&local_volume_integrals, &reduction_names, &det_jacobian, &mesh ](
137  const auto tensor_tag_v, const auto& tensor) noexcept {
138  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
139  for (size_t i = 0; i < tensor.size(); ++i) {
140  reduction_names.push_back("VolumeIntegral(" +
141  db::tag_name<tensor_tag>() +
142  tensor.component_suffix(i) + ")");
143  local_volume_integrals.push_back(
144  definite_integral(det_jacobian * tensor[i], mesh));
145  }
146  return 0;
147  };
149  record_integrals(tmpl::type_<Tensors>{}, tensors));
150 
151  // Send data to reduction observer
152  auto& local_observer =
153  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
154  cache)
155  .ckLocalBranch();
156  Parallel::simple_action<observers::Actions::ContributeReductionData>(
157  local_observer,
159  observation_value,
160  typename Metavariables::element_observation_type{}),
161  std::string{"/volume_integrals"}, reduction_names,
162  ReductionData{static_cast<double>(observation_value), local_volume,
163  local_volume_integrals});
164  }
165 };
166 
167 /// \cond
168 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
169  typename EventRegistrars>
170 PUP::able::PUP_ID ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
171  tmpl::list<Tensors...>,
172  EventRegistrars>::my_PUP_ID =
173  0; // NOLINT
174 /// \endcond
175 } // namespace Events
176 } // namespace dg
Parallel::ConstGlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
observers::ObservationId
A type-erased identifier that combines the identifier's type and hash used to uniquely identify an ob...
Definition: ObservationId.hpp:42
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:689
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:64
Options.hpp
Tags.hpp
vector
dg::Events::ObserveVolumeIntegrals
Definition: ObserveVolumeIntegrals.hpp:41
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ComputeNonconservativeBoundaryFluxes.hpp:23
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:47
Frame
Definition: IndexType.hpp:36
DefiniteIntegral.hpp
domain::Tags::DetInvJacobian
The determinant of the inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:197
Prefixes.hpp
OptionString
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
TMPL.hpp
ConstGlobalCache.hpp
string