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 
14 #include "DataStructures/DataVector.hpp"
16 #include "Domain/ElementMap.hpp"
17 #include "Domain/Tags.hpp"
18 #include "IO/Observer/Helpers.hpp"
19 #include "IO/Observer/ObservationId.hpp"
20 #include "IO/Observer/ObserverComponent.hpp"
21 #include "IO/Observer/ReductionActions.hpp"
23 #include "Options/Options.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/TMPL.hpp"
32 #include "Utilities/TaggedTuple.hpp"
33 
34 /// \cond
35 namespace Frame {
36 struct Inertial;
37 } // namespace Frame
38 /// \endcond
39 
40 namespace dg {
41 namespace Events {
42 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
43  typename EventRegistrars>
45 
46 namespace Registrars {
47 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors>
48 // Presence of size_t template argument requires to define this struct
49 // instead of using Registration::Registrar alias.
51  template <typename RegistrarList>
52  using f = Events::ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
53  Tensors, RegistrarList>;
54 };
55 } // namespace Registrars
56 
57 template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
58  typename EventRegistrars =
60  VolumeDim, ObservationValueTag, Tensors>>>
62 
63 /*!
64  * \ingroup DiscontinuousGalerkinGroup
65  * \brief %Observe the volume integrals of the tensors over the domain.
66  *
67  * Writes reduction quantities:
68  * - `ObservationValueTag`
69  * - `Volume` = volume of the domain
70  * - `VolumeIntegral(*)` = volume integral of the tensor
71  *
72  * \warning Currently, only one reduction observation event can be
73  * triggered at a given observation value. Causing multiple events to run at
74  * once will produce unpredictable results.
75  */
76 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
77  typename EventRegistrars>
78 class ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
79  tmpl::list<Tensors...>, EventRegistrars>
80  : public Event<EventRegistrars> {
81  private:
82  using VolumeIntegralDatum =
84 
85  using ReductionData = tmpl::wrap<
86  tmpl::list<Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
89  Parallel::ReductionData>;
90 
91  public:
92  /// \cond
93  explicit ObserveVolumeIntegrals(CkMigrateMessage* /*unused*/) noexcept {}
94  using PUP::able::register_constructor;
96  /// \endcond
97 
98  using options = tmpl::list<>;
99  static constexpr OptionString help =
100  "Observe the volume integrals of the tensors over the domain.\n"
101  "\n"
102  "Writes reduction quantities:\n"
103  " * ObservationValueTag\n"
104  " * Volume = volume of the domain\n"
105  " * VolumeIntegral(*) = volume integral of the tensor\n"
106  "\n"
107  "Warning: Currently, only one reduction observation event can be\n"
108  "triggered at a given observation value. Causing multiple events to\n"
109  "run at once will produce unpredictable results.";
110 
111  ObserveVolumeIntegrals() = default;
112 
113  using observed_reduction_data_tags =
114  observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
115 
116  using argument_tags =
117  tmpl::list<ObservationValueTag, domain::Tags::Mesh<VolumeDim>,
120  Tensors...>;
121 
122  template <typename Metavariables, typename ArrayIndex,
123  typename ParallelComponent>
124  void operator()(
125  const db::const_item_type<ObservationValueTag>& observation_value,
126  const Mesh<VolumeDim>& mesh,
127  const ElementMap<VolumeDim, Frame::Inertial>& element_map,
128  const tnsr::I<DataVector, VolumeDim, Frame::Logical>& logical_coordinates,
129  const db::const_item_type<Tensors>&... tensors,
131  const ArrayIndex& /*array_index*/,
132  const ParallelComponent* const /*meta*/) const noexcept {
133  // Determinant of Jacobian is needed because integral is performed in
134  // logical coords. Currently not initialized in the Metavariables
135  // as of PR #2048 (https://github.com/sxs-collaboration/spectre/pull/2048)
136  const DataVector det_jacobian =
137  get(determinant(element_map.jacobian(logical_coordinates)));
138  const double local_volume = definite_integral(det_jacobian, mesh);
139 
140  std::vector<double> local_volume_integrals{};
141  std::vector<std::string> reduction_names = {
142  db::tag_name<ObservationValueTag>(), "Volume"};
143  const auto record_integrals =
144  [&local_volume_integrals, &reduction_names, &det_jacobian, &mesh ](
145  const auto tensor_tag_v, const auto& tensor) noexcept {
146  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
147  for (size_t i = 0; i < tensor.size(); ++i) {
148  reduction_names.push_back("VolumeIntegral(" +
149  db::tag_name<tensor_tag>() +
150  tensor.component_suffix(i) + ")");
151  local_volume_integrals.push_back(
152  definite_integral(det_jacobian * tensor[i], mesh));
153  }
154  return 0;
155  };
157  record_integrals(tmpl::type_<Tensors>{}, tensors));
158 
159  // Send data to reduction observer
160  auto& local_observer =
161  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
162  cache)
163  .ckLocalBranch();
164  Parallel::simple_action<observers::Actions::ContributeReductionData>(
165  local_observer,
167  observation_value,
168  typename Metavariables::element_observation_type{}),
169  std::string{"/volume_integrals"}, reduction_names,
170  ReductionData{static_cast<double>(observation_value), local_volume,
171  local_volume_integrals});
172  }
173 };
174 
175 /// \cond
176 template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
177  typename EventRegistrars>
178 PUP::able::PUP_ID ObserveVolumeIntegrals<VolumeDim, ObservationValueTag,
179  tmpl::list<Tensors...>,
180  EventRegistrars>::my_PUP_ID =
181  0; // NOLINT
182 /// \endcond
183 } // namespace Events
184 } // namespace dg
Definition: ObserveVolumeIntegrals.hpp:44
Defines function definite_integral.
Definition: ObserveVolumeIntegrals.hpp:50
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:562
Definition: Digraph.hpp:11
Functionality related to discontinuous Galerkin schemes.
Definition: ApplyBoundaryFluxesLocalTimeStepping.hpp:33
Defines classes and functions for making classes creatable from input files.
Define prefixes for DataBox tags.
The coordinate map from logical to grid coordinate.
Definition: Tags.hpp:113
Defines macros to allow serialization of abstract template base classes.
Definition: Completion.hpp:13
Scalar< T > determinant(const Tensor< T, Symm, index_list< Index0, Index1 >> &tensor)
Computes the determinant of a rank-2 Tensor tensor.
Definition: Determinant.hpp:138
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
The CoordinateMap for the Element from the Logical frame to the TargetFrame
Definition: ElementMap.hpp:33
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
Defines function for taking the determinant of a rank-2 tensor.
Indicates the Frame that a TensorIndexType is in.
Definition: IndexType.hpp:36
Define simple functions for constant expressions.
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:136
tnsr::I< DataVector, VolumeDim, Frame::Logical > logical_coordinates(const Mesh< VolumeDim > &mesh) noexcept
Compute the logical coordinates in an Element.
Definition: LogicalCoordinates.cpp:20
Stores a collection of function values.
Definition: DataVector.hpp:42
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
Function for adding two std::vectors of double component-wise.
Definition: Functional.hpp:312
double definite_integral(const DataVector &integrand, const Mesh< Dim > &mesh) noexcept
Compute the definite integral of a function over a manifold.
Defines tags related to domain quantities.
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 data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:64