Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : // Tests for this file are combined with the ones for
5 : // ObserveTimeSeriesOnHorizon.hpp in
6 : // Test_ObserveFieldsAndTimeSeriesOnHorizon.cpp since they are so similar.
7 :
8 : #pragma once
9 :
10 : #include <cstddef>
11 : #include <string>
12 : #include <utility>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/TagName.hpp"
17 : #include "DataStructures/LinkedMessageId.hpp"
18 : #include "IO/H5/TensorData.hpp"
19 : #include "IO/Observer/ObserverComponent.hpp"
20 : #include "IO/Observer/ReductionActions.hpp"
21 : #include "IO/Observer/Tags.hpp"
22 : #include "IO/Observer/VolumeActions.hpp"
23 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
24 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
25 : #include "NumericalAlgorithms/Strahlkorper/IO/FillYlmLegendAndData.hpp"
26 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
27 : #include "NumericalAlgorithms/Strahlkorper/Tags.hpp"
28 : #include "Parallel/GlobalCache.hpp"
29 : #include "Parallel/Invoke.hpp"
30 : #include "Parallel/Local.hpp"
31 : #include "Parallel/Reduction.hpp"
32 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
33 : #include "ParallelAlgorithms/ApparentHorizonFinder/Protocols/Callback.hpp"
34 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
35 : #include "Utilities/ErrorHandling/Error.hpp"
36 : #include "Utilities/Gsl.hpp"
37 : #include "Utilities/PrettyType.hpp"
38 : #include "Utilities/ProtocolHelpers.hpp"
39 : #include "Utilities/TMPL.hpp"
40 :
41 0 : namespace ah ::callbacks {
42 : /*!
43 : * \brief A `ah::protocols::Callback` that outputs 2D "volume" data on a
44 : * surface and the surface's spherical harmonic data.
45 : *
46 : * \details The columns of spherical harmonic data written take the form
47 : *
48 : * \code
49 : * [Time, {Frame}ExpansionCenter_x, {Frame}ExpansionCenter_y,
50 : * {Frame}ExpansionCenter_z, Lmax, coef(0,0), ... coef(Lmax,Lmax)]
51 : * \endcode
52 : *
53 : * where `coef(l,m)` refers to the strahlkorper coefficients stored and defined
54 : * by `ylm::Strahlkorper::coefficients() const`, and `Frame` is
55 : * `HorizonMetavars::frame`. It is assumed that $l_{\mathrm{max}} =
56 : * m_{\mathrm{max}}$.
57 : *
58 : * \note Currently, $l_{\mathrm{max}}$ for a given surface does not change over
59 : * the course of the simulation, which means that the total number of columns of
60 : * coefficients that we need to write is also constant. The current
61 : * implementation of writing the coefficients at one time assumes
62 : * $l_{\mathrm{max}}$ of a surface remains constant. If and when in the future
63 : * functionality for an adaptive $l_{\mathrm{max}}$ is added, the implementation
64 : * for writing the coefficients will need to be updated to account for this. One
65 : * possible way to address this is to have a known maximum $l_{\mathrm{max}}$
66 : * for a given surface and write all coefficients up to that maximum
67 : * $l_{\mathrm{max}}$.
68 : */
69 : template <typename TagsToObserve, typename HorizonMetavars>
70 1 : struct ObserveFieldsOnHorizon : tt::ConformsTo<ah::protocols::Callback> {
71 0 : using Fr = typename HorizonMetavars::frame;
72 :
73 0 : using const_global_cache_tags = tmpl::list<observers::Tags::SurfaceFileName>;
74 :
75 : template <typename DbTags, typename Metavariables>
76 0 : static void apply(const db::DataBox<DbTags>& box,
77 : Parallel::GlobalCache<Metavariables>& cache,
78 : const FastFlow::Status /*status*/) {
79 : const auto& current_time = db::get<ah::Tags::CurrentTime>(box).value();
80 : const auto& strahlkorper = get<ylm::Tags::Strahlkorper<Fr>>(box);
81 : const ylm::Spherepack& ylm = strahlkorper.ylm_spherepack();
82 :
83 : // Output the inertial-frame coordinates of the Strahlkorper.
84 : // Note that these coordinates are not Spherepack-evenly-distributed over
85 : // the inertial-frame sphere (they are Spherepack-evenly-distributed over
86 : // the frame sphere).
87 : std::vector<TensorComponent> tensor_components;
88 : const auto& inertial_strahlkorper_coords =
89 : get<ylm::Tags::CartesianCoords<::Frame::Inertial>>(box);
90 : tensor_components.push_back(
91 : {"InertialCoordinates_x"s, get<0>(inertial_strahlkorper_coords)});
92 : tensor_components.push_back(
93 : {"InertialCoordinates_y"s, get<1>(inertial_strahlkorper_coords)});
94 : tensor_components.push_back(
95 : {"InertialCoordinates_z"s, get<2>(inertial_strahlkorper_coords)});
96 :
97 : tmpl::for_each<TagsToObserve>([&box, &tensor_components](auto tag_v) {
98 : using Tag = tmpl::type_from<decltype(tag_v)>;
99 : const auto tag_name = db::tag_name<Tag>();
100 : const auto& tensor = get<Tag>(box);
101 : for (size_t i = 0; i < tensor.size(); ++i) {
102 : tensor_components.emplace_back(tag_name + tensor.component_suffix(i),
103 : tensor[i]);
104 : }
105 : });
106 :
107 : const std::string& surface_name = pretty_type::name<HorizonMetavars>();
108 : const std::string subfile_path{std::string{"/"} + surface_name};
109 : const std::vector<size_t> extents_vector{
110 : {ylm.physical_extents()[0], ylm.physical_extents()[1]}};
111 : const std::vector<Spectral::Basis> bases_vector{
112 : 2, Spectral::Basis::SphericalHarmonic};
113 : const std::vector<Spectral::Quadrature> quadratures_vector{
114 : {Spectral::Quadrature::Gauss, Spectral::Quadrature::Equiangular}};
115 : const observers::ObservationId observation_id{current_time.id,
116 : subfile_path + ".vol"};
117 :
118 : auto& proxy = Parallel::get_parallel_component<
119 : observers::ObserverWriter<Metavariables>>(cache);
120 :
121 : // We call this on proxy[0] because the 0th element of a NodeGroup is
122 : // always guaranteed to be present.
123 : Parallel::threaded_action<observers::ThreadedActions::WriteVolumeData>(
124 : proxy[0], Parallel::get<observers::Tags::SurfaceFileName>(cache),
125 : subfile_path, observation_id,
126 : std::vector<ElementVolumeData>{{surface_name, tensor_components,
127 : extents_vector, bases_vector,
128 : quadratures_vector}});
129 :
130 : std::vector<std::string> ylm_legend{};
131 : std::vector<double> ylm_data{};
132 :
133 : // When option LMax is present in the global cache,
134 : // check that the Strahlkorper resolution does not exceed it.
135 : // LMax sets a maximum resolution for
136 : // the Strahlkorper resolution (when adjusted using adaptive
137 : // horizon finding) and the number of (zero padded) columns to output,
138 : // to avoid H5 errors caused by different rows (corresponding to
139 : // different times) having different numbers of columns.
140 : // If this option is not present in the cache, skip the check and just write
141 : // the modes of the Strahlkorper, assuming that the Strahlkorper
142 : // resolution never changes.
143 : size_t max_l_to_write = strahlkorper.l_max();
144 : if constexpr (Parallel::is_in_global_cache<Metavariables, ah::Tags::LMax>) {
145 : const auto& max_resolution_and_output_l =
146 : Parallel::get<ah::Tags::LMax>(cache);
147 : if (UNLIKELY(max_resolution_and_output_l < strahlkorper.l_max())) {
148 : ERROR("The option LMax ("
149 : << max_resolution_and_output_l << ") is smaller than the current "
150 : << "Strahlkorper resolution L=" << strahlkorper.l_max()
151 : << ". Increase LMax or decrease "
152 : << "Strahlkorper resolution L to avoid truncating data.");
153 : }
154 : max_l_to_write = max_resolution_and_output_l;
155 : }
156 :
157 : ylm::fill_ylm_legend_and_data(make_not_null(&ylm_legend),
158 : make_not_null(&ylm_data), strahlkorper,
159 : current_time.id, max_l_to_write);
160 :
161 : const std::string ylm_subfile_name{std::string{"/"} + surface_name +
162 : "_Ylm"};
163 :
164 : Parallel::threaded_action<
165 : observers::ThreadedActions::WriteReductionDataRow>(
166 : proxy[0], ylm_subfile_name, std::move(ylm_legend),
167 : std::make_tuple(std::move(ylm_data)));
168 : }
169 : };
170 : } // namespace ah::callbacks
|
