Observe.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <string>
7 
9 #include "Domain/Tags.hpp"
11 #include "IO/Observer/Actions.hpp"
12 #include "IO/Observer/ArrayComponentId.hpp"
13 #include "IO/Observer/Helpers.hpp"
14 #include "IO/Observer/ObservationId.hpp"
15 #include "IO/Observer/ObserverComponent.hpp"
16 #include "IO/Observer/ReductionActions.hpp"
17 #include "IO/Observer/VolumeActions.hpp"
20 #include "Parallel/Reduction.hpp"
21 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
22 #include "Utilities/Functional.hpp"
23 #include "Utilities/MakeString.hpp"
24 #include "Utilities/Numeric.hpp"
25 #include "Utilities/TMPL.hpp"
27 
28 namespace Poisson {
29 namespace Actions {
30 
31 /*!
32  * \brief Action to observe volume and reduction data for the Poisson system
33  *
34  * This action observes the following:
35  * - Reduction data:
36  * - `"Iteration"`: The linear solver iteration
37  * - `"NumberOfPoints"`: The total number of grid points across all elements
38  * - `"L2Error"`: The standard L2 vector norm of the pointwise difference
39  * between the numerical and analytic solution \f$\sqrt{\sum_i \left(
40  * u^\mathrm{numerical}_i - u^\mathrm{analytic}_i\right)^2}\f$ over the grid
41  * points across all elements.
42  * - Volume data:
43  * - `Poisson::Field::name()`: The numerical solution
44  * \f$u^\mathrm{numerical}\f$
45  * - `Poisson::Field::name() + "Analytic"`: The analytic solution
46  * \f$u^\mathrm{analytic}\f$
47  * - `Poisson::Field::name() + "Error"`: The pointwise error
48  * \f$u^\mathrm{numerical} - u^\mathrm{analytic}\f$
49  * - `"InertialCoordinates_{x,y,z}"`
50  *
51  * Uses:
52  * - Metavariables:
53  * - `analytic_solution_tag`
54  * - Items required by `observers::Observer<Metavariables>`
55  * - DataBox:
56  * - `LinearSolver::Tags::IterationId`
57  * - `Tags::Mesh<Dim>`
58  * - `Poisson::Field`
59  * - `Tags::Coordinates<Dim, Frame::Inertial>`
60  *
61  * \note This action can be adjusted before compiling an executable to observe
62  * only the desired quantities.
63  */
64 struct Observe {
65  private:
66  using observed_reduction_data = Parallel::ReductionData<
72 
73  public:
74  // Compile-time interface for observers
75  using observed_reduction_data_tags =
76  observers::make_reduction_data_tags<tmpl::list<observed_reduction_data>>;
77 
78  template <typename... DbTags, typename... InboxTags, typename Metavariables,
79  size_t Dim, typename ActionList, typename ParallelComponent>
80  static auto apply(db::DataBox<tmpl::list<DbTags...>>& box,
83  const ElementIndex<Dim>& array_index,
84  const ActionList /*meta*/,
85  const ParallelComponent* const /*meta*/) noexcept {
86  const auto& iteration_id = get<LinearSolver::Tags::IterationId>(box);
87  const auto& mesh = get<Tags::Mesh<Dim>>(box);
88  const std::string element_name = MakeString{} << ElementId<Dim>(array_index)
89  << '/';
90 
91  // Retrieve the current numeric solution
92  const auto& field = get<Poisson::Field>(box);
93 
94  // Compute the analytic solution
95  const auto& inertial_coordinates =
96  db::get<::Tags::Coordinates<Dim, Frame::Inertial>>(box);
97  const auto field_analytic = get<Poisson::Field>(
98  Parallel::get<typename Metavariables::analytic_solution_tag>(cache)
99  .variables(inertial_coordinates, tmpl::list<Poisson::Field>{}));
100 
101  // Compute error between numeric and analytic solutions
102  const DataVector field_error = get(field) - get(field_analytic);
103 
104  // Compute l2 error squared over local element
105  const double local_l2_error_square = alg::accumulate(
106  field_error, 0.0, funcl::Plus<funcl::Identity, funcl::Square<>>{});
107 
108  // Collect volume data
109  // Remove tensor types, only storing individual components
110  std::vector<TensorComponent> components;
111  components.reserve(3 + Dim);
112  components.emplace_back(element_name + Poisson::Field::name(), get(field));
113  components.emplace_back(element_name + Poisson::Field::name() + "Analytic",
114  get(field_analytic));
115  components.emplace_back(element_name + Poisson::Field::name() + "Error",
116  field_error);
117  components.emplace_back(element_name + "InertialCoordinates_x",
118  get<0>(inertial_coordinates));
119  if (Dim >= 2) {
120  components.emplace_back(element_name + "InertialCoordinates_y",
121  inertial_coordinates.get(1));
122  }
123  if (Dim >= 3) {
124  components.emplace_back(element_name + "InertialCoordinates_z",
125  inertial_coordinates.get(2));
126  }
127 
128  // Send data to volume observer
129  auto& local_observer =
130  *Parallel::get_parallel_component<observers::Observer<Metavariables>>(
131  cache)
132  .ckLocalBranch();
133  Parallel::simple_action<observers::Actions::ContributeVolumeData>(
134  local_observer, observers::ObservationId(iteration_id),
135  std::string{"/element_data"},
139  std::move(components), mesh.extents());
140 
141  // Send data to reduction observer
142  Parallel::simple_action<observers::Actions::ContributeReductionData>(
143  local_observer, observers::ObservationId(iteration_id),
144  std::string{"/element_data"},
145  std::vector<std::string>{"Iteration", "NumberOfPoints", "L2Error"},
146  observed_reduction_data{iteration_id.step_number,
147  mesh.number_of_grid_points(),
148  local_l2_error_square});
149 
150  return std::forward_as_tuple(std::move(box));
151  }
152 };
153 } // namespace Actions
154 } // namespace Poisson
Defines class tuples::TaggedTuple.
Function for squaring a quantity.
Definition: Functional.hpp:299
Defines DataBox tags for the linear solver.
The identity higher order function object.
Definition: Functional.hpp:214
Items related to solving a Poisson equation .
Definition: Actions.hpp:6
An ID type that identifies both the parallel component and the index in the parallel component...
Definition: ArrayComponentId.hpp:27
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
Make a string by streaming into object.
Definition: MakeString.hpp:16
Action to observe volume and reduction data for the Poisson system.
Definition: Observe.hpp:64
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
Defines DataBox tags for the Poisson system.
Functional for computing + of two objects.
Definition: Functional.hpp:234
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:76
Stores a collection of function values.
Definition: DataVector.hpp:46
Functional for computing sqrt on an object.
Definition: Functional.hpp:269
Wraps the template metaprogramming library used (brigand)
Defines tags related to domain quantities.
Definition: SolvePoissonProblem.hpp:38
A class for indexing a Charm array by Element.
Definition: ElementIndex.hpp:53
A type-erased identifier that combines the identifier&#39;s type and hash used to uniquely identify an ob...
Definition: ObservationId.hpp:43
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
Definition: ComputeTimeDerivative.hpp:28
The array index used for indexing Chare Arrays, mostly an implementation detail.
Definition: ArrayIndex.hpp:20