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 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/DataBox/ObservationBox.hpp"
13 : #include "DataStructures/DataBox/TagName.hpp"
14 : #include "DataStructures/DataVector.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Domain/Structure/ElementId.hpp"
17 : #include "Domain/Tags.hpp"
18 : #include "IO/Observer/GetSectionObservationKey.hpp"
19 : #include "IO/Observer/Helpers.hpp"
20 : #include "IO/Observer/ObservationId.hpp"
21 : #include "IO/Observer/ObserverComponent.hpp"
22 : #include "IO/Observer/ReductionActions.hpp"
23 : #include "IO/Observer/TypeOfObservation.hpp"
24 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
25 : #include "Options/String.hpp"
26 : #include "Parallel/ArrayComponentId.hpp"
27 : #include "Parallel/ArrayIndex.hpp"
28 : #include "Parallel/GlobalCache.hpp"
29 : #include "Parallel/Invoke.hpp"
30 : #include "Parallel/Local.hpp"
31 : #include "Parallel/Reduction.hpp"
32 : #include "ParallelAlgorithms/Events/Tags.hpp"
33 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
34 : #include "Utilities/ErrorHandling/Assert.hpp"
35 : #include "Utilities/ErrorHandling/Error.hpp"
36 : #include "Utilities/Functional.hpp"
37 : #include "Utilities/Numeric.hpp"
38 : #include "Utilities/OptionalHelpers.hpp"
39 : #include "Utilities/Serialization/CharmPupable.hpp"
40 : #include "Utilities/TMPL.hpp"
41 :
42 : namespace Events {
43 : /// @{
44 : /*!
45 : * \brief Find the extremum of a Scalar<DataVector> over
46 : * all elements, as well as the value of other functions
47 : * at the location of that extremum.
48 : *
49 : *
50 : * Here is an example of an input file:
51 : *
52 : * \snippet Test_ObserveAtExtremum.cpp input_file_examples
53 : *
54 : * \par Array sections
55 : * This event supports sections (see `Parallel::Section`). Set the
56 : * `ArraySectionIdTag` template parameter to split up observations into subsets
57 : * of elements. The `observers::Tags::ObservationKey<ArraySectionIdTag>` must be
58 : * available in the DataBox. It identifies the section and is used as a suffix
59 : * for the path in the output file.
60 : */
61 : template <typename ObservableTensorTagsList,
62 : typename NonTensorComputeTagsList = tmpl::list<>,
63 : typename ArraySectionIdTag = void>
64 1 : class ObserveAtExtremum;
65 :
66 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
67 : typename ArraySectionIdTag>
68 0 : class ObserveAtExtremum<tmpl::list<ObservableTensorTags...>,
69 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag>
70 : : public Event {
71 : private:
72 : /// Reduction data will contain the time, and either the maximum
73 : /// or the minimum of a function
74 : template <typename MinMaxFunctional>
75 1 : using ReductionData = Parallel::ReductionData<
76 : // Observation value
77 : Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
78 : // Maximum of first component of a vector
79 : Parallel::ReductionDatum<std::vector<double>, MinMaxFunctional>>;
80 :
81 : /// Information about the scalar whose extremum we search for,
82 : /// the type of extremum, and the other tensors to observer at
83 : /// that extremum
84 1 : struct ObserveTensors {
85 0 : static constexpr Options::String help = {
86 : "The tensor to extremize, and other tensors to observe."};
87 :
88 0 : struct Name {
89 0 : using type = std::string;
90 0 : static constexpr Options::String help = {
91 : "The name of the scalar to extremize."};
92 : };
93 :
94 0 : struct ExtremumType {
95 0 : using type = std::string;
96 0 : static constexpr Options::String help = {
97 : "The type of extremum -- either Min or Max."};
98 : };
99 :
100 0 : struct AdditionalData {
101 0 : using type = std::vector<std::string>;
102 0 : static constexpr Options::String help = {
103 : "List of other tensors to observe at the extremum"};
104 : };
105 :
106 0 : using options = tmpl::list<Name, ExtremumType, AdditionalData>;
107 :
108 0 : ObserveTensors() = default;
109 :
110 0 : ObserveTensors(std::string in_scalar, std::string in_extremum_type,
111 : std::vector<std::string> in_additional_data,
112 : const Options::Context& context = {});
113 :
114 0 : std::string scalar_name;
115 0 : std::string extremum_type;
116 0 : std::vector<std::string> additional_data{};
117 : };
118 :
119 : public:
120 : /// The name of the subfile inside the HDF5 file
121 1 : struct SubfileName {
122 0 : using type = std::string;
123 0 : static constexpr Options::String help = {
124 : "The name of the subfile inside the HDF5 file without an extension and "
125 : "without a preceding '/'."};
126 : };
127 : /// The scalar to extremize, and other tensors to observe at extremum
128 1 : struct TensorsToObserve {
129 0 : using type = ObserveTensors;
130 0 : static constexpr Options::String help = {
131 : "Struct specifying the scalar to extremize, the type of extremum "
132 : "and other tensors to observe at that extremum."};
133 : };
134 :
135 0 : explicit ObserveAtExtremum(CkMigrateMessage* msg);
136 : using PUP::able::register_constructor;
137 0 : WRAPPED_PUPable_decl_template(ObserveAtExtremum); // NOLINT
138 :
139 0 : using options = tmpl::list<SubfileName, TensorsToObserve>;
140 :
141 0 : static constexpr Options::String help =
142 : "Observe extremum of a scalar in the DataBox.\n"
143 : "\n"
144 : "Writes reduction quantities:\n"
145 : " * Observation value (e.g. Time or IterationId)\n"
146 : " * Extremum value of the desired scalar\n"
147 : " * Additional data at extremum\n";
148 :
149 0 : ObserveAtExtremum() = default;
150 :
151 0 : ObserveAtExtremum(std::string subfile_name,
152 : ObserveTensors observe_tensors);
153 :
154 0 : using observed_reduction_data_tags = observers::make_reduction_data_tags<
155 : tmpl::list<ReductionData<funcl::Max<>>, ReductionData<funcl::Min<>>>>;
156 :
157 0 : using compute_tags_for_observation_box =
158 : tmpl::list<ObservableTensorTags..., NonTensorComputeTags...>;
159 :
160 0 : using return_tags = tmpl::list<>;
161 0 : using argument_tags = tmpl::list<::Tags::ObservationBox>;
162 :
163 : template <typename ComputeTagsList, typename DataBoxType,
164 : typename Metavariables, size_t VolumeDim,
165 : typename ParallelComponent>
166 0 : void operator()(const ObservationBox<ComputeTagsList, DataBoxType>& box,
167 : Parallel::GlobalCache<Metavariables>& cache,
168 : const ElementId<VolumeDim>& array_index,
169 : const ParallelComponent* const /*meta*/,
170 : const ObservationValue& observation_value) const;
171 :
172 0 : using observation_registration_tags = tmpl::list<::Tags::DataBox>;
173 :
174 : template <typename DbTagsList>
175 : std::optional<
176 : std::pair<observers::TypeOfObservation, observers::ObservationKey>>
177 0 : get_observation_type_and_key_for_registration(
178 : const db::DataBox<DbTagsList>& box) const {
179 : const std::optional<std::string> section_observation_key =
180 : observers::get_section_observation_key<ArraySectionIdTag>(box);
181 : if (not section_observation_key.has_value()) {
182 : return std::nullopt;
183 : }
184 : return {{observers::TypeOfObservation::Reduction,
185 : observers::ObservationKey(
186 : subfile_path_ + section_observation_key.value() + ".dat")}};
187 : }
188 :
189 0 : using is_ready_argument_tags = tmpl::list<>;
190 :
191 : template <typename Metavariables, typename ArrayIndex, typename Component>
192 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
193 : const ArrayIndex& /*array_index*/,
194 : const Component* const /*meta*/) const {
195 : return true;
196 : }
197 :
198 1 : bool needs_evolved_variables() const override { return true; }
199 :
200 : // NOLINTNEXTLINE(google-runtime-references)
201 0 : void pup(PUP::er& p) override;
202 :
203 : private:
204 0 : std::string subfile_path_;
205 0 : std::string scalar_name_;
206 0 : std::string extremum_type_;
207 0 : std::vector<std::string> additional_tensor_names_{};
208 : };
209 : /// @}
210 :
211 : /// \cond
212 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
213 : typename ArraySectionIdTag>
214 : ObserveAtExtremum<tmpl::list<ObservableTensorTags...>,
215 : tmpl::list<NonTensorComputeTags...>,
216 : ArraySectionIdTag>::ObserveAtExtremum(CkMigrateMessage* msg)
217 : : Event(msg) {}
218 :
219 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
220 : typename ArraySectionIdTag>
221 : ObserveAtExtremum<
222 : tmpl::list<ObservableTensorTags...>, tmpl::list<NonTensorComputeTags...>,
223 : ArraySectionIdTag>::ObserveAtExtremum(std::string subfile_name,
224 : ObserveTensors observe_tensors)
225 : : subfile_path_("/" + subfile_name),
226 : scalar_name_(std::move(observe_tensors.scalar_name)),
227 : extremum_type_(std::move(observe_tensors.extremum_type)),
228 : additional_tensor_names_(std::move(observe_tensors.additional_data)) {}
229 :
230 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
231 : typename ArraySectionIdTag>
232 : ObserveAtExtremum<tmpl::list<ObservableTensorTags...>,
233 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag>::
234 : ObserveTensors::ObserveTensors(std::string in_scalar,
235 : std::string in_extremum_type,
236 : std::vector<std::string> in_additional_data,
237 : const Options::Context& context)
238 : : scalar_name(std::move(in_scalar)),
239 : extremum_type(std::move(in_extremum_type)),
240 : additional_data(std::move(in_additional_data)) {
241 : if (((scalar_name != db::tag_name<ObservableTensorTags>()) and ...)) {
242 : PARSE_ERROR(
243 : context, "Tensor '"
244 : << scalar_name << "' is not known. Known tensors are: "
245 : << ((db::tag_name<ObservableTensorTags>() + ",") + ...));
246 : }
247 :
248 : tmpl::for_each<tmpl::list<ObservableTensorTags...>>(
249 : [this, &context](auto tag_v) {
250 : using tag = tmpl::type_from<decltype(tag_v)>;
251 : const std::string tensor_name = db::tag_name<tag>();
252 : if (tensor_name == scalar_name) {
253 : if constexpr (tt::is_a_v<std::optional, typename tag::type>) {
254 : if (tag::type::value_type::rank() != 0) {
255 : PARSE_ERROR(context,
256 : "ObserveAtExtremum can only observe scalars!");
257 : }
258 : } else if (tag::type::rank() != 0) {
259 : PARSE_ERROR(context, "ObserveAtExtremum can only observe scalars!");
260 : }
261 : }
262 : });
263 :
264 : if (extremum_type != "Max" and extremum_type != "Min") {
265 : PARSE_ERROR(context, "Extremum type " << extremum_type
266 : << " not recognized; use Max or Min");
267 : }
268 : for (const auto& tensor : additional_data) {
269 : if (((tensor != db::tag_name<ObservableTensorTags>()) and ...)) {
270 : PARSE_ERROR(context,
271 : "Tensor '"
272 : << tensor << "' is not known. Known tensors are: "
273 : << ((db::tag_name<ObservableTensorTags>() + ",") + ...));
274 : }
275 : }
276 : }
277 :
278 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
279 : typename ArraySectionIdTag>
280 : template <typename ComputeTagsList, typename DataBoxType,
281 : typename Metavariables, size_t VolumeDim, typename ParallelComponent>
282 : void ObserveAtExtremum<tmpl::list<ObservableTensorTags...>,
283 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag>::
284 : operator()(const ObservationBox<ComputeTagsList, DataBoxType>& box,
285 : Parallel::GlobalCache<Metavariables>& cache,
286 : const ElementId<VolumeDim>& array_index,
287 : const ParallelComponent* const /*meta*/,
288 : const ObservationValue& observation_value) const {
289 : // Skip observation on elements that are not part of a section
290 : const std::optional<std::string> section_observation_key =
291 : observers::get_section_observation_key<ArraySectionIdTag>(box);
292 : if (not section_observation_key.has_value()) {
293 : return;
294 : }
295 :
296 : using tensor_tags = tmpl::list<ObservableTensorTags...>;
297 :
298 : // Vector that will contain the local extremum, and the value
299 : // of other tensors at that extremum
300 : std::vector<double> data_to_reduce{};
301 : // Vector containing a description of the data to be reduced.
302 : std::vector<std::string> legend{observation_value.name};
303 : // Location of the local extremum
304 : size_t index_of_extremum = 0;
305 : // First, look for local extremum of desired scalar
306 : tmpl::for_each<tensor_tags>([this, &box, &data_to_reduce, &legend,
307 : &index_of_extremum](auto tag_v) {
308 : using tag = tmpl::type_from<decltype(tag_v)>;
309 : const std::string tensor_name = db::tag_name<tag>();
310 : if (tensor_name == scalar_name_) {
311 : if (UNLIKELY(not has_value(get<tag>(box)))) {
312 : ERROR("Cannot observe a norm of '"
313 : << tensor_name
314 : << "' because it is a std::optional and wasn't able to be "
315 : "computed. This can happen when you try to observe errors "
316 : "without an analytic solution.");
317 : }
318 : const auto& scalar = value(get<tag>(box));
319 : const auto components = get<1>(scalar.get_vector_of_data());
320 : if (components.size() > 1) {
321 : ERROR("Extremum should be taken on a scalar, yet we have "
322 : << components.size() << " components in tensor " << tensor_name);
323 : }
324 : for (size_t i = 1; i < components[0].size(); i++) {
325 : if ((extremum_type_ == "Max" and
326 : (components[0][i] > components[0][index_of_extremum])) or
327 : (extremum_type_ == "Min" and
328 : (components[0][i] < components[0][index_of_extremum]))) {
329 : index_of_extremum = i;
330 : }
331 : }
332 : data_to_reduce.push_back(components[0][index_of_extremum]);
333 : if (extremum_type_ == "Max") {
334 : legend.push_back("Max(" + scalar_name_ + ")");
335 : } else {
336 : legend.push_back("Min(" + scalar_name_ + ")");
337 : }
338 : }
339 : });
340 : // Now get value of additional tensors at extremum
341 : tmpl::for_each<tensor_tags>([this, &box, &data_to_reduce, &legend,
342 : &index_of_extremum](auto tag_v) {
343 : using tag = tmpl::type_from<decltype(tag_v)>;
344 : const std::string tensor_name = db::tag_name<tag>();
345 : for (size_t i = 0; i < additional_tensor_names_.size(); ++i)
346 : if (tensor_name == additional_tensor_names_[i]) {
347 : if (UNLIKELY(not has_value(get<tag>(box)))) {
348 : ERROR("Cannot observe a norm of '"
349 : << tensor_name
350 : << "' because it is a std::optional and wasn't able to be "
351 : "computed. This can happen when you try to observe errors "
352 : "without an analytic solution.");
353 : }
354 : const auto& tensor = value(get<tag>(box));
355 : const auto [component_names, components] = tensor.get_vector_of_data();
356 : for (size_t j = 0; j < components.size(); j++) {
357 : data_to_reduce.push_back(components[j][index_of_extremum]);
358 : if (components.size() > 1) {
359 : legend.push_back("At" + scalar_name_ + extremum_type_ + "(" +
360 : tensor_name + "_" + component_names[j] + ")");
361 : } else {
362 : legend.push_back("At" + scalar_name_ + extremum_type_ + "(" +
363 : tensor_name + ")");
364 : }
365 : }
366 : }
367 : });
368 :
369 : // Send data to reduction observer
370 : auto& local_observer = *Parallel::local_branch(
371 : Parallel::get_parallel_component<observers::Observer<Metavariables>>(
372 : cache));
373 : const std::string subfile_path_with_suffix =
374 : subfile_path_ + section_observation_key.value();
375 :
376 : if (extremum_type_ == "Max") {
377 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
378 : local_observer,
379 : observers::ObservationId(observation_value.value,
380 : subfile_path_with_suffix + ".dat"),
381 : Parallel::make_array_component_id<ParallelComponent>(array_index),
382 : subfile_path_with_suffix, std::move(legend),
383 : ReductionData<funcl::Max<>>{observation_value.value,
384 : std::move(data_to_reduce)});
385 : } else {
386 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
387 : local_observer,
388 : observers::ObservationId(observation_value.value,
389 : subfile_path_with_suffix + ".dat"),
390 : Parallel::make_array_component_id<ParallelComponent>(array_index),
391 : subfile_path_with_suffix, std::move(legend),
392 : ReductionData<funcl::Min<>>{observation_value.value,
393 : std::move(data_to_reduce)});
394 : }
395 : }
396 :
397 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
398 : typename ArraySectionIdTag>
399 : void ObserveAtExtremum<tmpl::list<ObservableTensorTags...>,
400 : tmpl::list<NonTensorComputeTags...>,
401 : ArraySectionIdTag>::pup(PUP::er& p) {
402 : Event::pup(p);
403 : p | subfile_path_;
404 : p | scalar_name_;
405 : p | extremum_type_;
406 : p | additional_tensor_names_;
407 : }
408 :
409 : template <typename... ObservableTensorTags, typename... NonTensorComputeTags,
410 : typename ArraySectionIdTag>
411 : PUP::able::PUP_ID ObserveAtExtremum<tmpl::list<ObservableTensorTags...>,
412 : tmpl::list<NonTensorComputeTags...>,
413 : ArraySectionIdTag>::my_PUP_ID =
414 : 0; // NOLINT
415 : /// \endcond
416 : } // namespace Events
|
