Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cmath>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <pup.h>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <utility>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/ObservationBox.hpp"
17 : #include "DataStructures/DataBox/TagName.hpp"
18 : #include "DataStructures/DataVector.hpp"
19 : #include "DataStructures/Tensor/Tensor.hpp"
20 : #include "Domain/Structure/ElementId.hpp"
21 : #include "Domain/Tags.hpp"
22 : #include "IO/Observer/GetSectionObservationKey.hpp"
23 : #include "IO/Observer/Helpers.hpp"
24 : #include "IO/Observer/ObservationId.hpp"
25 : #include "IO/Observer/ObserverComponent.hpp"
26 : #include "IO/Observer/ReductionActions.hpp"
27 : #include "IO/Observer/TypeOfObservation.hpp"
28 : #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp"
29 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
30 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
31 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
32 : #include "Options/String.hpp"
33 : #include "Parallel/ArrayIndex.hpp"
34 : #include "Parallel/GlobalCache.hpp"
35 : #include "Parallel/Invoke.hpp"
36 : #include "Parallel/Local.hpp"
37 : #include "Parallel/Reduction.hpp"
38 : #include "Parallel/TypeTraits.hpp"
39 : #include "ParallelAlgorithms/Events/Tags.hpp"
40 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
41 : #include "Utilities/ErrorHandling/Assert.hpp"
42 : #include "Utilities/ErrorHandling/Error.hpp"
43 : #include "Utilities/Functional.hpp"
44 : #include "Utilities/OptionalHelpers.hpp"
45 : #include "Utilities/PrettyType.hpp"
46 : #include "Utilities/Serialization/CharmPupable.hpp"
47 : #include "Utilities/TMPL.hpp"
48 :
49 : namespace Events {
50 : /// @{
51 : /*!
52 : * \brief Compute norms of tensors in the DataBox and write them to disk.
53 : *
54 : * The L1 norm is computed as the mean absolute value, so
55 : *
56 : * \f{align*}{
57 : * L_1(u)=\frac{1}{N}\sum_{i=0}^{N-1} |u_i|
58 : * \f}
59 : *
60 : * where \f$N\f$ is the number of grid points.
61 : *
62 : * The L2 norm is computed as the RMS, so
63 : *
64 : * \f{align*}{
65 : * L_2(u)=\sqrt{\frac{1}{N}\sum_{i=0}^{N-1} u_i^2}
66 : * \f}
67 : *
68 : * The norm can be taken for each individual component, or summed over
69 : * components. For the max/min it is then the max/min over all components, while
70 : * for the L1 norm we have (for a 3d vector, 2d and 1d are similar)
71 : *
72 : * \f{align*}{
73 : * L_1(v^k)=\frac{1}{N}\sum_{i=0}^{N-1} \left[|v^x_i| + |v^y_i|
74 : * + |v^z_i|\right]
75 : * \f}
76 : *
77 : * and for the L2 norm
78 : *
79 : * \f{align*}{
80 : * L_2(v^k)=\sqrt{\frac{1}{N}\sum_{i=0}^{N-1} \left[(v^x_i)^2 + (v^y_i)^2
81 : * + (v^z_i)^2\right]}
82 : * \f}
83 : *
84 : * The L1 integral norm is:
85 : *
86 : * \begin{equation}
87 : * L_{1,\mathrm{int}}(v^k) = \frac{1}{V}\int_\Omega \left[
88 : * |v^x_i| + |v^y_i| + |v^z_i|\right] \mathrm{d}V
89 : * \end{equation}
90 : *
91 : * The L2 integral norm is:
92 : *
93 : * \begin{equation}
94 : * L_{2,\mathrm{int}}(v^k) = \sqrt{\frac{1}{V}\int_\Omega \left[
95 : * (v^x_i)^2 + (v^y_i)^2 + (v^z_i)^2\right] \mathrm{d}V}
96 : * \end{equation}
97 : *
98 : * where $V=\int_\Omega$ is the volume of the entire domain in inertial
99 : * coordinates.
100 : *
101 : * VolumeIntegral only computes the volume integral without any normalization.
102 : *
103 : * Here is an example of an input file:
104 : *
105 : * \snippet Test_ObserveNorms.cpp input_file_examples
106 : *
107 : * \note The `NonTensorComputeTags` are intended to be used for `Variables`
108 : * compute tags like `Tags::DerivCompute`
109 : *
110 : * \par Array sections
111 : * This event supports sections (see `Parallel::Section`). Set the
112 : * `ArraySectionIdTag` template parameter to split up observations into subsets
113 : * of elements. The `observers::Tags::ObservationKey<ArraySectionIdTag>` must be
114 : * available in the DataBox. It identifies the section and is used as a suffix
115 : * for the path in the output file.
116 : *
117 : * \par Option name
118 : * The `OptionName` template parameter is used to give the event a name in the
119 : * input file. If it is not specified, the name defaults to "ObserveNorms". If
120 : * you have multiple `ObserveNorms` events in the input file, you must specify a
121 : * unique name for each one. This can happen, for example, if you want to
122 : * observe norms the full domain and also over a section of the domain.
123 : */
124 : template <typename ObservableTensorTagsList,
125 : typename NonTensorComputeTagsList = tmpl::list<>,
126 : typename ArraySectionIdTag = void, typename OptionName = void>
127 1 : class ObserveNorms;
128 :
129 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
130 : typename ArraySectionIdTag, typename OptionName>
131 0 : class ObserveNorms<tmpl::list<ObservableTensorTags...>,
132 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag,
133 : OptionName> : public Event {
134 : private:
135 0 : struct ObserveTensor {
136 0 : static constexpr Options::String help = {
137 : "The tensor to reduce, and how to reduce it."};
138 :
139 0 : struct Name {
140 0 : using type = std::string;
141 0 : static constexpr Options::String help = {
142 : "The name of the tensor to observe."};
143 : };
144 0 : struct NormType {
145 0 : using type = std::string;
146 0 : static constexpr Options::String help = {
147 : "The type of norm to use. Must be one of Max, Min, L1Norm, "
148 : "L1IntegralNorm, L2Norm, L2IntegralNorm, or VolumeIntegral."};
149 : };
150 0 : struct Components {
151 0 : using type = std::string;
152 0 : static constexpr Options::String help = {
153 : "How to handle tensor components. Must be Individual or Sum."};
154 : };
155 :
156 0 : using options = tmpl::list<Name, NormType, Components>;
157 :
158 0 : ObserveTensor() = default;
159 :
160 0 : ObserveTensor(std::string in_tensor, std::string in_norm_type,
161 : std::string in_components,
162 : const Options::Context& context = {});
163 :
164 0 : std::string tensor{};
165 0 : std::string norm_type{};
166 0 : std::string components{};
167 : };
168 :
169 0 : using ReductionData = Parallel::ReductionData<
170 : // Observation value
171 : Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
172 : // Number of grid points
173 : Parallel::ReductionDatum<size_t, funcl::Plus<>>,
174 : // Total volume
175 : Parallel::ReductionDatum<double, funcl::Plus<>>,
176 : // Max
177 : Parallel::ReductionDatum<std::vector<double>,
178 : funcl::ElementWise<funcl::Max<>>>,
179 : // Min
180 : Parallel::ReductionDatum<std::vector<double>,
181 : funcl::ElementWise<funcl::Min<>>>,
182 : // L1Norm
183 : Parallel::ReductionDatum<
184 : std::vector<double>, funcl::ElementWise<funcl::Plus<>>,
185 : funcl::ElementWise<funcl::Divides<>>, std::index_sequence<1>>,
186 : // L1IntegralNorm
187 : Parallel::ReductionDatum<
188 : std::vector<double>, funcl::ElementWise<funcl::Plus<>>,
189 : funcl::ElementWise<funcl::Divides<>>, std::index_sequence<2>>,
190 : // L2Norm
191 : Parallel::ReductionDatum<
192 : std::vector<double>, funcl::ElementWise<funcl::Plus<>>,
193 : funcl::ElementWise<funcl::Sqrt<funcl::Divides<>>>,
194 : std::index_sequence<1>>,
195 : // L2IntegralNorm
196 : Parallel::ReductionDatum<
197 : std::vector<double>, funcl::ElementWise<funcl::Plus<>>,
198 : funcl::ElementWise<funcl::Sqrt<funcl::Divides<>>>,
199 : std::index_sequence<2>>,
200 : // VolumeIntegral
201 : Parallel::ReductionDatum<std::vector<double>,
202 : funcl::ElementWise<funcl::Plus<>>>>;
203 :
204 : public:
205 0 : static std::string name() {
206 : if constexpr (std::is_same_v<OptionName, void>) {
207 : return "ObserveNorms";
208 : } else {
209 : return pretty_type::name<OptionName>();
210 : }
211 : }
212 :
213 : /// The name of the subfile inside the HDF5 file
214 1 : struct SubfileName {
215 0 : using type = std::string;
216 0 : static constexpr Options::String help = {
217 : "The name of the subfile inside the HDF5 file without an extension and "
218 : "without a preceding '/'."};
219 : };
220 : /// The tensor to observe and how to do the reduction
221 1 : struct TensorsToObserve {
222 0 : using type = std::vector<ObserveTensor>;
223 0 : static constexpr Options::String help = {
224 : "List specifying each tensor to observe and how it is reduced."};
225 : };
226 :
227 0 : explicit ObserveNorms(CkMigrateMessage* msg);
228 : using PUP::able::register_constructor;
229 0 : WRAPPED_PUPable_decl_template(ObserveNorms); // NOLINT
230 :
231 0 : using options = tmpl::list<SubfileName, TensorsToObserve>;
232 :
233 0 : static constexpr Options::String help =
234 : "Observe norms of tensors in the DataBox.\n"
235 : "\n"
236 : "You can choose the norm type for each observation. Note that 'L1Norm'\n"
237 : "(mean absolute value) and 'L2Norm' (root mean square) emphasize "
238 : "regions\n"
239 : "of the domain with many grid points, whereas 'L1IntegralNorm' and\n"
240 : "'L2IntegralNorm' emphasize regions of the domain with large volume.\n"
241 : "Choose wisely! When in doubt, try the 'L2Norm' first.\n"
242 : "\n"
243 : "Writes reduction quantities:\n"
244 : " * Observation value (e.g. Time or IterationId)\n"
245 : " * NumberOfPoints = total number of points in the domain\n"
246 : " * Volume = total volume of the domain in inertial coordinates\n"
247 : " * Max values\n"
248 : " * Min values\n"
249 : " * L1-norm values\n"
250 : " * L1 integral norm values\n"
251 : " * L2-norm values\n"
252 : " * L2 integral norm values\n"
253 : " * Volume integral values\n";
254 :
255 0 : ObserveNorms() = default;
256 :
257 0 : ObserveNorms(const std::string& subfile_name,
258 : const std::vector<ObserveTensor>& observe_tensors);
259 :
260 0 : using observed_reduction_data_tags =
261 : observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
262 :
263 0 : using compute_tags_for_observation_box =
264 : tmpl::list<ObservableTensorTags..., NonTensorComputeTags...>;
265 :
266 0 : using return_tags = tmpl::list<>;
267 0 : using argument_tags = tmpl::list<::Tags::ObservationBox>;
268 :
269 : template <typename TensorToObserveTag, typename ComputeTagsList,
270 : typename DataBoxType, size_t Dim>
271 0 : void observe_norms_impl(
272 : gsl::not_null<
273 : std::unordered_map<std::string, std::pair<std::vector<double>,
274 : std::vector<std::string>>>*>
275 : norm_values_and_names,
276 : const ObservationBox<ComputeTagsList, DataBoxType>& box,
277 : const Mesh<Dim>& mesh, const DataVector& det_jacobian,
278 : size_t number_of_points) const;
279 :
280 : template <typename ComputeTagsList, typename DataBoxType,
281 : typename Metavariables, size_t VolumeDim,
282 : typename ParallelComponent>
283 0 : void operator()(const ObservationBox<ComputeTagsList, DataBoxType>& box,
284 : Parallel::GlobalCache<Metavariables>& cache,
285 : const ElementId<VolumeDim>& array_index,
286 : const ParallelComponent* const /*meta*/,
287 : const ObservationValue& observation_value) const;
288 :
289 0 : using observation_registration_tags = tmpl::list<::Tags::DataBox>;
290 :
291 : template <typename DbTagsList>
292 : std::optional<
293 : std::pair<observers::TypeOfObservation, observers::ObservationKey>>
294 0 : get_observation_type_and_key_for_registration(
295 : const db::DataBox<DbTagsList>& box) const {
296 : const std::optional<std::string> section_observation_key =
297 : observers::get_section_observation_key<ArraySectionIdTag>(box);
298 : if (not section_observation_key.has_value()) {
299 : return std::nullopt;
300 : }
301 : return {{observers::TypeOfObservation::Reduction,
302 : observers::ObservationKey(
303 : subfile_path_ + section_observation_key.value() + ".dat")}};
304 : }
305 :
306 0 : using is_ready_argument_tags = tmpl::list<>;
307 :
308 : template <typename Metavariables, typename ArrayIndex, typename Component>
309 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
310 : const ArrayIndex& /*array_index*/,
311 : const Component* const /*meta*/) const {
312 : return true;
313 : }
314 :
315 1 : bool needs_evolved_variables() const override { return true; }
316 :
317 : // NOLINTNEXTLINE(google-runtime-references)
318 0 : void pup(PUP::er& p) override;
319 :
320 : private:
321 0 : std::string subfile_path_;
322 0 : std::vector<std::string> tensor_names_{};
323 0 : std::vector<std::string> tensor_norm_types_{};
324 0 : std::vector<std::string> tensor_components_{};
325 : };
326 : /// @}
327 :
328 : /// \cond
329 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
330 : typename ArraySectionIdTag, typename OptionName>
331 : ObserveNorms<tmpl::list<ObservableTensorTags...>,
332 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag,
333 : OptionName>::ObserveNorms(CkMigrateMessage* msg)
334 : : Event(msg) {}
335 :
336 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
337 : typename ArraySectionIdTag, typename OptionName>
338 : ObserveNorms<tmpl::list<ObservableTensorTags...>,
339 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag,
340 : OptionName>::ObserveNorms(const std::string& subfile_name,
341 : const std::vector<ObserveTensor>&
342 : observe_tensors)
343 : : subfile_path_("/" + subfile_name) {
344 : tensor_names_.reserve(observe_tensors.size());
345 : tensor_norm_types_.reserve(observe_tensors.size());
346 : tensor_components_.reserve(observe_tensors.size());
347 : for (const auto& observe_tensor : observe_tensors) {
348 : tensor_names_.push_back(observe_tensor.tensor);
349 : tensor_norm_types_.push_back(observe_tensor.norm_type);
350 : tensor_components_.push_back(observe_tensor.components);
351 : }
352 : }
353 :
354 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
355 : typename ArraySectionIdTag, typename OptionName>
356 : ObserveNorms<
357 : tmpl::list<ObservableTensorTags...>, tmpl::list<NonTensorComputeTags...>,
358 : ArraySectionIdTag,
359 : OptionName>::ObserveTensor::ObserveTensor(std::string in_tensor,
360 : std::string in_norm_type,
361 : std::string in_components,
362 : const Options::Context& context)
363 : : tensor(std::move(in_tensor)),
364 : norm_type(std::move(in_norm_type)),
365 : components(std::move(in_components)) {
366 : if (((tensor != db::tag_name<ObservableTensorTags>()) and ...)) {
367 : PARSE_ERROR(
368 : context, "Tensor '"
369 : << tensor << "' is not known. Known tensors are: "
370 : << ((db::tag_name<ObservableTensorTags>() + ",") + ...));
371 : }
372 : if (norm_type != "Max" and norm_type != "Min" and norm_type != "L1Norm" and
373 : norm_type != "L1IntegralNorm" and norm_type != "L2Norm" and
374 : norm_type != "L2IntegralNorm" and norm_type != "VolumeIntegral") {
375 : PARSE_ERROR(
376 : context,
377 : "NormType must be one of Max, Min, L1Norm, L1IntegralNorm, L2Norm, "
378 : "L2IntegralNorm, or VolumeIntegral not "
379 : << norm_type);
380 : }
381 : if (components != "Individual" and components != "Sum") {
382 : PARSE_ERROR(context,
383 : "Components must be Individual or Sum, not " << components);
384 : }
385 : }
386 :
387 : // implementation of ObserveNorms::operator() factored out to save on compile
388 : // time and compile memory
389 : namespace ObserveNorms_impl {
390 : void check_norm_is_observable(const std::string& tensor_name,
391 : bool tag_has_value);
392 :
393 : template <size_t Dim>
394 : void fill_norm_values_and_names(
395 : gsl::not_null<std::unordered_map<
396 : std::string, std::pair<std::vector<double>, std::vector<std::string>>>*>
397 : norm_values_and_names,
398 : const std::pair<std::vector<std::string>, std::vector<DataVector>>&
399 : names_and_components,
400 : const Mesh<Dim>& mesh, const DataVector& det_jacobian,
401 : const std::string& tensor_name, const std::string& tensor_norm_type,
402 : const std::string& tensor_component, size_t number_of_points);
403 :
404 : // Expand complex data into real and imaginary parts, or just forward real data
405 : std::pair<std::vector<std::string>, std::vector<DataVector>>
406 : split_complex_vector_of_data(
407 : std::pair<std::vector<std::string>, std::vector<DataVector>>&&
408 : names_and_components);
409 : std::pair<std::vector<std::string>, std::vector<DataVector>>
410 : split_complex_vector_of_data(
411 : const std::pair<std::vector<std::string>, std::vector<ComplexDataVector>>&
412 : names_and_components);
413 : } // namespace ObserveNorms_impl
414 :
415 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
416 : typename ArraySectionIdTag, typename OptionName>
417 : template <typename TensorToObserveTag, typename ComputeTagsList,
418 : typename DataBoxType, size_t Dim>
419 : void ObserveNorms<tmpl::list<ObservableTensorTags...>,
420 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag,
421 : OptionName>::
422 : observe_norms_impl(
423 : const gsl::not_null<std::unordered_map<
424 : std::string,
425 : std::pair<std::vector<double>, std::vector<std::string>>>*>
426 : norm_values_and_names,
427 : const ObservationBox<ComputeTagsList, DataBoxType>& box,
428 : const Mesh<Dim>& mesh, const DataVector& det_jacobian,
429 : const size_t number_of_points) const {
430 : const std::string tensor_name = db::tag_name<TensorToObserveTag>();
431 : for (size_t i = 0; i < tensor_names_.size(); ++i) {
432 : if (tensor_name == tensor_names_[i]) {
433 : ObserveNorms_impl::check_norm_is_observable(
434 : tensor_name, has_value(get<TensorToObserveTag>(box)));
435 : ObserveNorms_impl::fill_norm_values_and_names(
436 : norm_values_and_names,
437 : ObserveNorms_impl::split_complex_vector_of_data(
438 : value(get<TensorToObserveTag>(box)).get_vector_of_data()),
439 : mesh, det_jacobian, tensor_name, tensor_norm_types_[i],
440 : tensor_components_[i], number_of_points);
441 : }
442 : }
443 : }
444 :
445 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
446 : typename ArraySectionIdTag, typename OptionName>
447 : template <typename ComputeTagsList, typename DataBoxType,
448 : typename Metavariables, size_t VolumeDim, typename ParallelComponent>
449 : void ObserveNorms<tmpl::list<ObservableTensorTags...>,
450 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag,
451 : OptionName>::
452 : operator()(const ObservationBox<ComputeTagsList, DataBoxType>& box,
453 : Parallel::GlobalCache<Metavariables>& cache,
454 : const ElementId<VolumeDim>& array_index,
455 : const ParallelComponent* const /*meta*/,
456 : const ObservationValue& observation_value) const {
457 : // Skip observation on elements that are not part of a section
458 : const std::optional<std::string> section_observation_key =
459 : observers::get_section_observation_key<ArraySectionIdTag>(box);
460 : if (not section_observation_key.has_value()) {
461 : return;
462 : }
463 :
464 : const auto& mesh = get<::Events::Tags::ObserverMesh<VolumeDim>>(box);
465 : const auto det_jacobian = [&box, &mesh]() -> DataVector {
466 : if constexpr (VolumeDim > 1 and
467 : db::tag_is_retrievable_v<::Events::Tags::ObserverCoordinates<
468 : VolumeDim, Frame::Inertial>,
469 : std::decay_t<decltype(box)>>) {
470 : if (mesh.basis(0) == Spectral::Basis::ZernikeB2) {
471 : // disk (2D) or cylinder (3D), need rho cartesian to polar jacobian
472 : const auto& inertial_coords = get<
473 : ::Events::Tags::ObserverCoordinates<VolumeDim, Frame::Inertial>>(
474 : box);
475 : const DataVector radius = sqrt(square(get<0>(inertial_coords)) +
476 : square(get<1>(inertial_coords)));
477 : return radius / get(get<::Events::Tags::ObserverDetInvJacobian<
478 : Frame::ElementLogical, Frame::Inertial>>(box));
479 : }
480 : if constexpr (VolumeDim == 3) {
481 : if (mesh.basis(0) == Spectral::Basis::ZernikeB3) {
482 : // ball, needs r^2 cartesian to spherical jacobian
483 : const auto& inertial_coords = get<
484 : ::Events::Tags::ObserverCoordinates<VolumeDim, Frame::Inertial>>(
485 : box);
486 : const DataVector r_squared = square(get<0>(inertial_coords)) +
487 : square(get<1>(inertial_coords)) +
488 : square(get<2>(inertial_coords));
489 : return r_squared /
490 : get(get<::Events::Tags::ObserverDetInvJacobian<
491 : Frame::ElementLogical, Frame::Inertial>>(box));
492 : } else if (mesh.basis(2) == Spectral::Basis::Cartoon) {
493 : if (mesh.quadrature(2) == Spectral::Quadrature::SphericalSymmetry) {
494 : // Spherical Symmetry, needs x^2 cartesian to spherical jacobian
495 : return square(get<0>(get<::Events::Tags::ObserverCoordinates<
496 : VolumeDim, Frame::Inertial>>(box))) /
497 : get(get<::Events::Tags::ObserverDetInvJacobian<
498 : Frame::ElementLogical, Frame::Inertial>>(box));
499 : } else {
500 : // Axial Symmetry, needs x cartesian to cylindrical jacobian
501 : ASSERT(mesh.quadrature(2) == Spectral::Quadrature::AxialSymmetry,
502 : "Unexpected quadrature " << mesh.quadrature(2)
503 : << " (expected AxialSymmetry)");
504 : return get<0>(get<::Events::Tags::ObserverCoordinates<
505 : VolumeDim, Frame::Inertial>>(box)) /
506 : get(get<::Events::Tags::ObserverDetInvJacobian<
507 : Frame::ElementLogical, Frame::Inertial>>(box));
508 : }
509 : }
510 : }
511 : }
512 : (void)mesh;
513 : return 1. / get(get<::Events::Tags::ObserverDetInvJacobian<
514 : Frame::ElementLogical, Frame::Inertial>>(box));
515 : }();
516 : const size_t number_of_points = mesh.number_of_grid_points();
517 : const double local_volume = definite_integral(det_jacobian, mesh);
518 :
519 : std::unordered_map<std::string,
520 : std::pair<std::vector<double>, std::vector<std::string>>>
521 : norm_values_and_names{};
522 : // Loop over ObservableTensorTags and see if it was requested to be observed.
523 : // This approach allows us to delay evaluating any compute tags until they're
524 : // actually needed for observing.
525 : (observe_norms_impl<ObservableTensorTags>(
526 : make_not_null(&norm_values_and_names), box, mesh, det_jacobian,
527 : number_of_points),
528 : ...);
529 :
530 : // Concatenate the legend info together.
531 : std::vector<std::string> legend{observation_value.name, "NumberOfPoints",
532 : "Volume"};
533 : legend.insert(legend.end(), norm_values_and_names["Max"].second.begin(),
534 : norm_values_and_names["Max"].second.end());
535 : legend.insert(legend.end(), norm_values_and_names["Min"].second.begin(),
536 : norm_values_and_names["Min"].second.end());
537 : legend.insert(legend.end(), norm_values_and_names["L1Norm"].second.begin(),
538 : norm_values_and_names["L1Norm"].second.end());
539 : legend.insert(legend.end(),
540 : norm_values_and_names["L1IntegralNorm"].second.begin(),
541 : norm_values_and_names["L1IntegralNorm"].second.end());
542 : legend.insert(legend.end(), norm_values_and_names["L2Norm"].second.begin(),
543 : norm_values_and_names["L2Norm"].second.end());
544 : legend.insert(legend.end(),
545 : norm_values_and_names["L2IntegralNorm"].second.begin(),
546 : norm_values_and_names["L2IntegralNorm"].second.end());
547 : legend.insert(legend.end(),
548 : norm_values_and_names["VolumeIntegral"].second.begin(),
549 : norm_values_and_names["VolumeIntegral"].second.end());
550 :
551 : const std::string subfile_path_with_suffix =
552 : subfile_path_ + section_observation_key.value();
553 : // Send data to reduction observer
554 : auto& local_observer = *Parallel::local_branch(
555 : Parallel::get_parallel_component<
556 : tmpl::conditional_t<Parallel::is_nodegroup_v<ParallelComponent>,
557 : observers::ObserverWriter<Metavariables>,
558 : observers::Observer<Metavariables>>>(cache));
559 : observers::ObservationId observation_id{observation_value.value,
560 : subfile_path_with_suffix + ".dat"};
561 : Parallel::ArrayComponentId array_component_id{
562 : std::add_pointer_t<ParallelComponent>{nullptr},
563 : Parallel::ArrayIndex<ElementId<VolumeDim>>(array_index)};
564 : ReductionData reduction_data{
565 : observation_value.value,
566 : number_of_points,
567 : local_volume,
568 : std::move(norm_values_and_names["Max"].first),
569 : std::move(norm_values_and_names["Min"].first),
570 : std::move(norm_values_and_names["L1Norm"].first),
571 : std::move(norm_values_and_names["L1IntegralNorm"].first),
572 : std::move(norm_values_and_names["L2Norm"].first),
573 : std::move(norm_values_and_names["L2IntegralNorm"].first),
574 : std::move(norm_values_and_names["VolumeIntegral"].first)};
575 :
576 : if constexpr (Parallel::is_nodegroup_v<ParallelComponent>) {
577 : Parallel::threaded_action<
578 : observers::ThreadedActions::CollectReductionDataOnNode>(
579 : local_observer, std::move(observation_id),
580 : std::move(array_component_id), subfile_path_with_suffix,
581 : std::move(legend), std::move(reduction_data));
582 : } else {
583 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
584 : local_observer, std::move(observation_id),
585 : std::move(array_component_id), subfile_path_with_suffix,
586 : std::move(legend), std::move(reduction_data));
587 : }
588 : }
589 :
590 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
591 : typename ArraySectionIdTag, typename OptionName>
592 : void ObserveNorms<tmpl::list<ObservableTensorTags...>,
593 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag,
594 : OptionName>::pup(PUP::er& p) {
595 : Event::pup(p);
596 : p | subfile_path_;
597 : p | tensor_names_;
598 : p | tensor_norm_types_;
599 : p | tensor_components_;
600 : }
601 :
602 : // NOLINTBEGIN(cppcoreguidelines-avoid-non-const-global-variables)
603 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
604 : typename ArraySectionIdTag, typename OptionName>
605 : PUP::able::PUP_ID ObserveNorms<tmpl::list<ObservableTensorTags...>,
606 : tmpl::list<NonTensorComputeTags...>,
607 : ArraySectionIdTag, OptionName>::my_PUP_ID = 0;
608 : // NOLINTEND(cppcoreguidelines-avoid-non-const-global-variables)
609 : /// \endcond
610 : } // namespace Events
|
