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