Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cmath>
8 : #include <limits>
9 : #include <tuple>
10 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "IO/Observer/ObserverComponent.hpp"
13 : #include "IO/Observer/ReductionActions.hpp"
14 : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
15 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
16 : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
17 : #include "Parallel/GlobalCache.hpp"
18 : #include "Parallel/Invoke.hpp"
19 : #include "Parallel/ParallelComponentHelpers.hpp"
20 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
21 : #include "ParallelAlgorithms/ApparentHorizonFinder/Protocols/Callback.hpp"
22 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
23 : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
24 : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
25 : #include "Utilities/Gsl.hpp"
26 : #include "Utilities/PrettyType.hpp"
27 : #include "Utilities/ProtocolHelpers.hpp"
28 :
29 : /// \cond
30 : namespace Frame {
31 : struct Grid;
32 : struct Distorted;
33 : struct Inertial;
34 : } // namespace Frame
35 : /// \endcond
36 :
37 : namespace intrp::callbacks {
38 : /*!
39 : * \brief Writes the center of an apparent horizon to disk in both the
40 : * `Frame` template parameter frame and Frame::Inertial frame. Intended to be
41 : * used in the `post_horizon_find_callbacks` list of an InterpolationTargetTag.
42 : *
43 : * The centers will be written to a subfile with the name
44 : * `/ApparentHorizons/TargetName_Centers.dat` where `TargetName` is the
45 : * pretty_type::name of the InterpolationTargetTag template parameter.
46 : *
47 : * The columns of the dat file are:
48 : * - %Time
49 : * - GridCenter_x
50 : * - GridCenter_y
51 : * - GridCenter_z
52 : * - InertialCenter_x
53 : * - InertialCenter_y
54 : * - InertialCenter_z
55 : *
56 : * The `Frame` template parameter must be either `::Frame::Grid` or
57 : * `::Frame::Distorted`. Even though the template parameter can be
58 : * `::Frame::Distorted`, we still write `GridCenter_?` because the centers of
59 : * the objects are the same in the Grid and Distorted frames.
60 : *
61 : * \note Requires ylm::Tags::Strahlkorper<Frame>
62 : * and ylm::Tags::CartesianCoords<Frame::Inertial> and
63 : * ylm::Tags::EuclideanAreaElement<Frame> to be in the DataBox
64 : * of the InterpolationTarget.
65 : */
66 : template <typename InterpolationTargetTag, typename Frame>
67 1 : struct ObserveCenters {
68 : template <typename DbTags, typename Metavariables, typename TemporalId>
69 0 : static void apply(const db::DataBox<DbTags>& box,
70 : Parallel::GlobalCache<Metavariables>& cache,
71 : const TemporalId& temporal_id) {
72 : static_assert(std::is_same_v<Frame, ::Frame::Grid> or
73 : std::is_same_v<Frame, ::Frame::Distorted>,
74 : "Frame must be either Grid or Distorted.");
75 : using HorizonTag = ylm::Tags::Strahlkorper<Frame>;
76 : using CoordsTag = ylm::Tags::CartesianCoords<::Frame::Inertial>;
77 : using EuclideanAreaElementTag = ylm::Tags::EuclideanAreaElement<Frame>;
78 : static_assert(db::tag_is_retrievable_v<HorizonTag, db::DataBox<DbTags>>,
79 : "DataBox must contain ylm::Tags::Strahlkorper<Frame>");
80 : static_assert(db::tag_is_retrievable_v<CoordsTag, db::DataBox<DbTags>>,
81 : "DataBox must contain "
82 : "ylm::Tags::CartesianCoords<Frame::Inertial>");
83 : static_assert(
84 : db::tag_is_retrievable_v<EuclideanAreaElementTag, db::DataBox<DbTags>>,
85 : "DataBox must contain ylm::Tags::EuclideanAreaElement<Frame>");
86 :
87 : // Only print the centers if we want to.
88 : if (not Parallel::get<ah::Tags::ObserveCenters>(cache)) {
89 : return;
90 : }
91 :
92 : const auto& grid_horizon = db::get<HorizonTag>(box);
93 : const std::array<double, 3> grid_center = grid_horizon.physical_center();
94 :
95 : // computes the inertial center to be the average value of the
96 : // inertial coordinates over the Strahlkorper, where the average is
97 : // computed by a surface integral.
98 : // Note that we use the Euclidean area element here, since we are trying
99 : // to find a geometric center of a surface without regard to GR.
100 : const auto& inertial_coords = db::get<CoordsTag>(box);
101 : const auto& area_element = db::get<EuclideanAreaElementTag>(box);
102 : std::array<double, 3> inertial_center =
103 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
104 : auto integrand = make_with_value<Scalar<DataVector>>(get(area_element), 0.);
105 : const double denominator = grid_horizon.ylm_spherepack().definite_integral(
106 : get(area_element).data());
107 : for (size_t i = 0; i < 3; ++i) {
108 : get(integrand) = get(area_element) * inertial_coords.get(i);
109 : gsl::at(inertial_center, i) =
110 : grid_horizon.ylm_spherepack().definite_integral(
111 : get(integrand).data()) /
112 : denominator;
113 : }
114 :
115 : // time, grid_x, grid_y, grid_z, inertial_x, inertial_y, inertial_z
116 : const auto center_tuple = std::make_tuple(
117 : intrp::InterpolationTarget_detail::get_temporal_id_value(temporal_id),
118 : grid_center[0], grid_center[1], grid_center[2], inertial_center[0],
119 : inertial_center[1], inertial_center[2]);
120 :
121 : auto& observer_writer_proxy = Parallel::get_parallel_component<
122 : observers::ObserverWriter<Metavariables>>(cache);
123 :
124 : const std::string subfile_path{"/ApparentHorizons/" +
125 : pretty_type::name<InterpolationTargetTag>() +
126 : "_Centers"};
127 :
128 : Parallel::threaded_action<
129 : observers::ThreadedActions::WriteReductionDataRow>(
130 : // Node 0 is always the writer
131 : observer_writer_proxy[0], subfile_path, legend_, center_tuple);
132 : }
133 :
134 : private:
135 0 : const static inline std::vector<std::string> legend_{
136 : {"Time", "GridCenter_x", "GridCenter_y", "GridCenter_z",
137 : "InertialCenter_x", "InertialCenter_y", "InertialCenter_z"}};
138 : };
139 : } // namespace intrp::callbacks
140 :
141 : namespace ah::callbacks {
142 : /*!
143 : * \brief Writes the center of an apparent horizon to disk in both the
144 : * `Frame` template parameter frame and Frame::Inertial frame. Intended to be
145 : * used in the `horizon_find_callbacks` list of a
146 : * `ah::protocols::HorizonMetavars`.
147 : *
148 : * The centers will be written to a subfile with the name
149 : * `/ApparentHorizons/TargetName_Centers.dat` where `TargetName` is the
150 : * pretty_type::name of the \p HorizonMetavars template parameter.
151 : *
152 : * The columns of the dat file are:
153 : * - %Time
154 : * - GridCenter_x
155 : * - GridCenter_y
156 : * - GridCenter_z
157 : * - InertialCenter_x
158 : * - InertialCenter_y
159 : * - InertialCenter_z
160 : *
161 : * The `Frame` template parameter must be either `::Frame::Grid` or
162 : * `::Frame::Distorted`. Even though the template parameter can be
163 : * `::Frame::Distorted`, we still write `GridCenter_?` because the centers of
164 : * the objects are the same in the Grid and Distorted frames.
165 : *
166 : * \note Requires `ylm::Tags::Strahlkorper<Frame>` and
167 : * `ylm::Tags::CartesianCoords<Frame::Inertial>` and
168 : * `ylm::Tags::EuclideanAreaElement<Frame>` to be in the DataBox of the horizon
169 : * finder.
170 : */
171 : template <typename HorizonMetavars>
172 1 : struct ObserveCenters : tt::ConformsTo<ah::protocols::Callback> {
173 : template <typename DbTags, typename Metavariables>
174 0 : static void apply(const db::DataBox<DbTags>& box,
175 : Parallel::GlobalCache<Metavariables>& cache,
176 : const FastFlow::Status /*status*/) {
177 : using frame = typename HorizonMetavars::frame;
178 : static_assert(std::is_same_v<frame, ::Frame::Grid> or
179 : std::is_same_v<frame, ::Frame::Distorted>,
180 : "Frame must be either Grid or Distorted.");
181 : using HorizonTag = ylm::Tags::Strahlkorper<frame>;
182 : using CoordsTag = ylm::Tags::CartesianCoords<::Frame::Inertial>;
183 : using EuclideanAreaElementTag = ylm::Tags::EuclideanAreaElement<frame>;
184 : static_assert(db::tag_is_retrievable_v<HorizonTag, db::DataBox<DbTags>>,
185 : "DataBox must contain ylm::Tags::Strahlkorper<Frame>");
186 : static_assert(db::tag_is_retrievable_v<CoordsTag, db::DataBox<DbTags>>,
187 : "DataBox must contain "
188 : "ylm::Tags::CartesianCoords<Frame::Inertial>");
189 : static_assert(
190 : db::tag_is_retrievable_v<EuclideanAreaElementTag, db::DataBox<DbTags>>,
191 : "DataBox must contain ylm::Tags::EuclideanAreaElement<Frame>");
192 :
193 : // Only print the centers if we want to.
194 : if (not Parallel::get<ah::Tags::ObserveCenters>(cache)) {
195 : return;
196 : }
197 :
198 : const auto& grid_horizon = db::get<HorizonTag>(box);
199 : const std::array<double, 3> grid_center = grid_horizon.physical_center();
200 :
201 : // computes the inertial center to be the average value of the
202 : // inertial coordinates over the Strahlkorper, where the average is
203 : // computed by a surface integral.
204 : // Note that we use the Euclidean area element here, since we are trying
205 : // to find a geometric center of a surface without regard to GR.
206 : const auto& inertial_coords = db::get<CoordsTag>(box);
207 : const auto& area_element = db::get<EuclideanAreaElementTag>(box);
208 : std::array<double, 3> inertial_center =
209 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
210 : auto integrand = make_with_value<Scalar<DataVector>>(get(area_element), 0.);
211 : const double denominator = grid_horizon.ylm_spherepack().definite_integral(
212 : get(area_element).data());
213 : for (size_t i = 0; i < 3; ++i) {
214 : get(integrand) = get(area_element) * inertial_coords.get(i);
215 : gsl::at(inertial_center, i) =
216 : grid_horizon.ylm_spherepack().definite_integral(
217 : get(integrand).data()) /
218 : denominator;
219 : }
220 :
221 : // time, grid_x, grid_y, grid_z, inertial_x, inertial_y, inertial_z
222 : const auto& current_time = db::get<ah::Tags::CurrentTime>(box).value();
223 : const auto center_tuple = std::make_tuple(
224 : current_time.id, grid_center[0], grid_center[1], grid_center[2],
225 : inertial_center[0], inertial_center[1], inertial_center[2]);
226 :
227 : auto& observer_writer_proxy = Parallel::get_parallel_component<
228 : observers::ObserverWriter<Metavariables>>(cache);
229 :
230 : const std::string subfile_path{"/ApparentHorizons/" +
231 : HorizonMetavars::name() + "_Centers"};
232 :
233 : Parallel::threaded_action<
234 : observers::ThreadedActions::WriteReductionDataRow>(
235 : // Node 0 is always the writer
236 : observer_writer_proxy[0], subfile_path, legend_, center_tuple);
237 : }
238 :
239 : private:
240 0 : const static inline std::vector<std::string> legend_{
241 : {"Time", "GridCenter_x", "GridCenter_y", "GridCenter_z",
242 : "InertialCenter_x", "InertialCenter_y", "InertialCenter_z"}};
243 : };
244 : } // namespace ah::callbacks
|
