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 <string>
8 : #include <utility>
9 : #include <vector>
10 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/DataBox/TagName.hpp"
13 : #include "IO/H5/TensorData.hpp"
14 : #include "IO/Observer/ObserverComponent.hpp"
15 : #include "IO/Observer/ReductionActions.hpp"
16 : #include "IO/Observer/Tags.hpp"
17 : #include "IO/Observer/VolumeActions.hpp"
18 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
19 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
20 : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
21 : #include "NumericalAlgorithms/Strahlkorper/IO/FillYlmLegendAndData.hpp"
22 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
23 : #include "NumericalAlgorithms/Strahlkorper/Tags.hpp"
24 : #include "Parallel/GlobalCache.hpp"
25 : #include "Parallel/Invoke.hpp"
26 : #include "Parallel/Local.hpp"
27 : #include "Parallel/Reduction.hpp"
28 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
29 : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
30 : #include "ParallelAlgorithms/Interpolation/Protocols/PostInterpolationCallback.hpp"
31 : #include "Utilities/ErrorHandling/Error.hpp"
32 : #include "Utilities/Functional.hpp"
33 : #include "Utilities/Gsl.hpp"
34 : #include "Utilities/PrettyType.hpp"
35 : #include "Utilities/ProtocolHelpers.hpp"
36 : #include "Utilities/TMPL.hpp"
37 :
38 : namespace intrp {
39 : namespace callbacks {
40 : /// \brief post_interpolation_callback that outputs 2D "volume" data on a
41 : /// surface and the surface's spherical harmonic data
42 : ///
43 : /// \details
44 : /// Uses:
45 : /// - Metavariables
46 : /// - `temporal_id`
47 : /// - DataBox:
48 : /// - `TagsToObserve` (each tag must be a Scalar<DataVector>)
49 : ///
50 : /// Conforms to the intrp::protocols::PostInterpolationCallback protocol
51 : ///
52 : /// For requirements on InterpolationTargetTag, see
53 : /// intrp::protocols::InterpolationTargetTag
54 : ///
55 : /// The columns of spherical harmonic data written take the form
56 : ///
57 : /// \code
58 : /// [Time, {Frame}ExpansionCenter_x, {Frame}ExpansionCenter_y,
59 : /// {Frame}ExpansionCenter_z, Lmax, coef(0,0), ... coef(Lmax,Lmax)]
60 : /// \endcode
61 : ///
62 : /// where `coef(l,m)` refers to the strahlkorper coefficients stored and defined
63 : /// by `ylm::Strahlkorper::coefficients() const`. It is assumed that
64 : /// \f$l_{max} = m_{max}\f$.
65 : ///
66 : /// \note Currently, \f$l_{max}\f$ for a given surface does not change over the
67 : /// course of the simulation, which means that the total number of columns of
68 : /// coefficients that we need to write is also constant. The current
69 : /// implementation of writing the coefficients at one time assumes \f$l_{max}\f$
70 : /// of a surface remains constant. If and when in the future functionality for
71 : /// an adaptive \f$l_{max}\f$ is added, the implementation for writing the
72 : /// coefficients will need to be updated to account for this. One possible way
73 : /// to address this is to have a known maximum \f$l_{max}\f$ for a given surface
74 : /// and write all coefficients up to that maximum \f$l_{max}\f$.
75 : template <typename TagsToObserve, typename InterpolationTargetTag,
76 : typename HorizonFrame>
77 1 : struct ObserveSurfaceData
78 : : tt::ConformsTo<intrp::protocols::PostInterpolationCallback> {
79 0 : static constexpr double fill_invalid_points_with =
80 : std::numeric_limits<double>::quiet_NaN();
81 :
82 0 : using const_global_cache_tags = tmpl::list<observers::Tags::SurfaceFileName>;
83 :
84 : template <typename DbTags, typename Metavariables, typename TemporalId>
85 0 : static void apply(const db::DataBox<DbTags>& box,
86 : Parallel::GlobalCache<Metavariables>& cache,
87 : const TemporalId& temporal_id) {
88 : const auto& strahlkorper = get<ylm::Tags::Strahlkorper<HorizonFrame>>(box);
89 : const ylm::Spherepack& ylm = strahlkorper.ylm_spherepack();
90 :
91 : // Output the inertial-frame coordinates of the Stralhlkorper.
92 : // Note that these coordinates are not
93 : // Spherepack-evenly-distributed over the inertial-frame sphere
94 : // (they are Spherepack-evenly-distributed over the HorizonFrame
95 : // sphere).
96 : std::vector<TensorComponent> tensor_components;
97 : if constexpr (db::tag_is_retrievable_v<
98 : ylm::Tags::CartesianCoords<::Frame::Inertial>,
99 : db::DataBox<DbTags>>) {
100 : const auto& inertial_strahlkorper_coords =
101 : get<ylm::Tags::CartesianCoords<::Frame::Inertial>>(box);
102 : tensor_components.push_back(
103 : {"InertialCoordinates_x"s, get<0>(inertial_strahlkorper_coords)});
104 : tensor_components.push_back(
105 : {"InertialCoordinates_y"s, get<1>(inertial_strahlkorper_coords)});
106 : tensor_components.push_back(
107 : {"InertialCoordinates_z"s, get<2>(inertial_strahlkorper_coords)});
108 : }
109 :
110 : // Output each tag if it is a scalar. Otherwise, throw a compile-time
111 : // error. This could be generalized to handle tensors of nonzero rank by
112 : // looping over the components, so each component could be visualized
113 : // separately as a scalar. But in practice, this generalization is
114 : // probably unnecessary, because Strahlkorpers are typically only
115 : // visualized with scalar quantities (used set the color at different
116 : // points on the surface).
117 : tmpl::for_each<TagsToObserve>([&box, &tensor_components](auto tag_v) {
118 : using Tag = tmpl::type_from<decltype(tag_v)>;
119 : const auto tag_name = db::tag_name<Tag>();
120 : const auto& tensor = get<Tag>(box);
121 : for (size_t i = 0; i < tensor.size(); ++i) {
122 : tensor_components.emplace_back(tag_name + tensor.component_suffix(i),
123 : tensor[i]);
124 : }
125 : });
126 :
127 : const std::string& surface_name =
128 : pretty_type::name<InterpolationTargetTag>();
129 : const std::string subfile_path{std::string{"/"} + surface_name};
130 : const std::vector<size_t> extents_vector{
131 : {ylm.physical_extents()[0], ylm.physical_extents()[1]}};
132 : const std::vector<Spectral::Basis> bases_vector{
133 : 2, Spectral::Basis::SphericalHarmonic};
134 : const std::vector<Spectral::Quadrature> quadratures_vector{
135 : {Spectral::Quadrature::Gauss, Spectral::Quadrature::Equiangular}};
136 : const double time =
137 : InterpolationTarget_detail::get_temporal_id_value(temporal_id);
138 : const observers::ObservationId observation_id{time, subfile_path + ".vol"};
139 :
140 : auto& proxy = Parallel::get_parallel_component<
141 : observers::ObserverWriter<Metavariables>>(cache);
142 :
143 : // We call this on proxy[0] because the 0th element of a NodeGroup is
144 : // always guaranteed to be present.
145 : Parallel::threaded_action<observers::ThreadedActions::WriteVolumeData>(
146 : proxy[0], Parallel::get<observers::Tags::SurfaceFileName>(cache),
147 : subfile_path, observation_id,
148 : std::vector<ElementVolumeData>{{surface_name, tensor_components,
149 : extents_vector, bases_vector,
150 : quadratures_vector}});
151 :
152 : std::vector<std::string> ylm_legend{};
153 : std::vector<double> ylm_data{};
154 :
155 : // When option LMax is present in the global cache,
156 : // check that the Strahlkorper resolution does not exceed it.
157 : // LMax sets a maximum resolution for
158 : // the Strahlkorper resolution (when adjusted using adaptive
159 : // horizon finding) and the number of (zero padded) columns to output,
160 : // to avoid H5 errors caused by different rows (corresponding to
161 : // different times) having different numbers of columns.
162 : // If this option is not present in the cache, skip the check and just write
163 : // the modes of the Strahlkorper, assuming that the Strahlkorper
164 : // resolution never changes.
165 : size_t max_l_to_write = strahlkorper.l_max();
166 : if constexpr (Parallel::is_in_global_cache<Metavariables, ah::Tags::LMax>) {
167 : const auto& max_resolution_and_output_l =
168 : Parallel::get<ah::Tags::LMax>(cache);
169 : if (UNLIKELY(max_resolution_and_output_l < strahlkorper.l_max())) {
170 : ERROR("The option LMax ("
171 : << max_resolution_and_output_l << ") is smaller than the current "
172 : << "Strahlkorper resolution L=" << strahlkorper.l_max()
173 : << ". Increase LMax or decrease "
174 : << "Strahlkorper resolution L to avoid truncating data.");
175 : }
176 : max_l_to_write = max_resolution_and_output_l;
177 : }
178 :
179 : ylm::fill_ylm_legend_and_data(make_not_null(&ylm_legend),
180 : make_not_null(&ylm_data), strahlkorper, time,
181 : max_l_to_write);
182 :
183 : const std::string ylm_subfile_name{std::string{"/"} + surface_name +
184 : "_Ylm"};
185 :
186 : Parallel::threaded_action<
187 : observers::ThreadedActions::WriteReductionDataRow>(
188 : proxy[0], ylm_subfile_name, std::move(ylm_legend),
189 : std::make_tuple(std::move(ylm_data)));
190 : }
191 : };
192 : } // namespace callbacks
193 : } // namespace intrp
|
