ObserveErrorNorms.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 "Domain/Tags.hpp"
14 #include "Evolution/EventsAndTriggers/Event.hpp"
15 #include "IO/Observer/Helpers.hpp"
16 #include "IO/Observer/ObservationId.hpp"
17 #include "IO/Observer/ObserverComponent.hpp" // IWYU pragma: keep
18 #include "IO/Observer/ReductionActions.hpp" // IWYU pragma: keep
19 #include "Options/Options.hpp"
22 #include "Parallel/Invoke.hpp"
23 #include "Parallel/Reduction.hpp"
24 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
25 #include "Time/Time.hpp"
27 #include "Utilities/Functional.hpp"
28 #include "Utilities/Numeric.hpp"
29 #include "Utilities/TMPL.hpp"
31 
32 /// \cond
33 namespace Frame {
34 struct Inertial;
35 } // namespace Frame
36 namespace Tags {
37 struct Time;
38 } // namespace Tags
39 /// \endcond
40 
41 namespace dg {
42 namespace Events {
43 template <size_t VolumeDim, typename Tensors, typename EventRegistrars>
45 
46 namespace Registrars {
47 template <size_t VolumeDim, typename Tensors>
49  template <typename RegistrarList>
51 };
52 } // namespace Registrars
53 
54 template <size_t VolumeDim, typename Tensors,
55  typename EventRegistrars =
56  tmpl::list<Registrars::ObserveErrorNorms<VolumeDim, Tensors>>>
57 class ObserveErrorNorms; // IWYU pragma: keep
58 
59 /*!
60  * \ingroup DiscontinuousGalerkinGroup
61  * \brief %Observe the RMS errors in the tensors compared to their
62  * analytic solution.
63  *
64  * Writes reduction quantities:
65  * - `Time`
66  * - `NumberOfPoints` = total number of points in the domain
67  * - `Error(*)` = RMS errors in `Tensors` =
68  * \f$\operatorname{RMS}\left(\sqrt{\sum_{\text{independent components}}\left[
69  * \text{value} - \text{analytic solution}\right]^2}\right)\f$
70  * over all points
71  *
72  * \warning Currently, only one reduction observation event can be
73  * triggered at a given time. Causing multiple events to run at once
74  * will produce unpredictable results.
75  */
76 template <size_t VolumeDim, typename... Tensors, typename EventRegistrars>
77 class ObserveErrorNorms<VolumeDim, tmpl::list<Tensors...>, EventRegistrars>
78  : public Event<EventRegistrars> {
79  private:
81 
82  template <typename Tag>
83  struct LocalSquareError {
84  using type = double;
85  };
86 
90  using ReductionData = tmpl::wrap<
91  tmpl::append<
92  tmpl::list<Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
94  tmpl::filled_list<L2ErrorDatum, sizeof...(Tensors)>>,
95  Parallel::ReductionData>;
96 
97  public:
98  /// \cond
99  explicit ObserveErrorNorms(CkMigrateMessage* /*unused*/) noexcept {}
100  using PUP::able::register_constructor;
102  /// \endcond
103 
104  using options = tmpl::list<>;
105  static constexpr OptionString help =
106  "Observe the RMS errors in the tensors compared to their analytic\n"
107  "solution.\n"
108  "\n"
109  "Writes reduction quantities:\n"
110  " * Time\n"
111  " * NumberOfPoints = total number of points in the domain\n"
112  " * Error(*) = RMS errors in Tensors (see online help details)\n"
113  "\n"
114  "Warning: Currently, only one reduction observation event can be\n"
115  "triggered at a given time. Causing multiple events to run at once\n"
116  "will produce unpredictable results.";
117 
118  ObserveErrorNorms() = default;
119 
120  using observed_reduction_data_tags =
121  observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
122 
123  using argument_tags = tmpl::list<::Tags::Time, coordinates_tag, Tensors...>;
124 
125  template <typename Metavariables, typename ArrayIndex,
126  typename ParallelComponent>
127  void operator()(const Time& time,
128  const db::item_type<coordinates_tag>& inertial_coordinates,
129  const db::item_type<Tensors>&... tensors,
131  const ArrayIndex& /*array_index*/,
132  const ParallelComponent* const /*meta*/) const noexcept {
133  const auto analytic_solution =
134  Parallel::get<OptionTags::AnalyticSolutionBase>(cache).variables(
135  inertial_coordinates, time.value(), tmpl::list<Tensors...>{});
136 
137  tuples::TaggedTuple<LocalSquareError<Tensors>...> local_square_errors;
138  const auto record_errors = [&analytic_solution, &local_square_errors](
139  const auto tensor_tag_v, const auto& tensor) noexcept {
140  using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
141  double local_square_error = 0.0;
142  for (size_t i = 0; i < tensor.size(); ++i) {
143  const auto error = tensor[i] - get<tensor_tag>(analytic_solution)[i];
144  local_square_error += alg::accumulate(square(error), 0.0);
145  }
146  get<LocalSquareError<tensor_tag>>(local_square_errors) =
147  local_square_error;
148  return 0;
149  };
150  expand_pack(record_errors(tmpl::type_<Tensors>{}, tensors)...);
151  const size_t num_points = get_first_argument(tensors...).begin()->size();
152 
153  // Send data to reduction observer
154  auto& local_observer =
155  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
156  cache)
157  .ckLocalBranch();
158  Parallel::simple_action<observers::Actions::ContributeReductionData>(
159  local_observer,
162  std::string{"/element_data"},
163  std::vector<std::string>{"Time", "NumberOfPoints",
164  ("Error(" + Tensors::name() + ")")...},
165  ReductionData{
166  time.value(), num_points,
167  std::move(get<LocalSquareError<Tensors>>(local_square_errors))...});
168  }
169 };
170 
171 /// \cond
172 template <size_t VolumeDim, typename... Tensors, typename EventRegistrars>
173 PUP::able::PUP_ID ObserveErrorNorms<VolumeDim, tmpl::list<Tensors...>,
174  EventRegistrars>::my_PUP_ID = 0; // NOLINT
175 /// \endcond
176 } // namespace Events
177 } // namespace dg
Definition: SolvePoissonProblem.hpp:67
Defines class tuples::TaggedTuple.
Definition: Digraph.hpp:11
decltype(auto) constexpr get_first_argument(T &&t, Ts &&...) noexcept
Returns the first argument of a parameter pack.
Definition: TMPL.hpp:569
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
Defines classes and functions for making classes creatable from input files.
Tag for compute item for current Time (from TimeId)
Definition: Tags.hpp:50
Defines macros to allow serialization of abstract template base classes.
Definition: Completion.hpp:13
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
Defines Time and TimeDelta.
double value() const noexcept
Approximate numerical value of the Time.
Definition: Time.hpp:51
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
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:27
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Definition: ObserveErrorNorms.hpp:44
Indicates the Frame that a TensorIndexType is in.
Definition: IndexType.hpp:36
Define simple functions for constant expressions.
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
Definition: ObserveErrorNorms.hpp:48
Functional for computing sqrt on an object.
Definition: Functional.hpp:273
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: SolvePoissonProblem.hpp:37
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.
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:52
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:64
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:545