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