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 <memory>
8 : #include <string>
9 : #include <tuple>
10 : #include <type_traits>
11 : #include <unordered_map>
12 : #include <unordered_set>
13 : #include <utility>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/Tag.hpp"
17 : #include "DataStructures/Tags/TempTensor.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "Domain/Creators/DomainCreator.hpp"
20 : #include "Domain/Creators/OptionTags.hpp"
21 : #include "Domain/Structure/BlockGroups.hpp"
22 : #include "IO/Logging/Tags.hpp"
23 : #include "IO/Logging/Verbosity.hpp"
24 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
25 : #include "NumericalAlgorithms/SphericalHarmonics/StrahlkorperFunctions.hpp"
26 : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
27 : #include "Options/Context.hpp"
28 : #include "Options/String.hpp"
29 : #include "Parallel/GlobalCache.hpp"
30 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
31 : #include "ParallelAlgorithms/ApparentHorizonFinder/OptionTags.hpp"
32 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
33 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
34 : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
35 : #include "ParallelAlgorithms/Interpolation/Protocols/ComputeTargetPoints.hpp"
36 : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
37 : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
38 : #include "Utilities/Gsl.hpp"
39 : #include "Utilities/PrettyType.hpp"
40 : #include "Utilities/ProtocolHelpers.hpp"
41 : #include "Utilities/Requires.hpp"
42 : #include "Utilities/TMPL.hpp"
43 : #include "Utilities/TaggedTuple.hpp"
44 :
45 : /// \cond
46 : class DataVector;
47 :
48 : namespace PUP {
49 : class er;
50 : } // namespace PUP
51 : namespace db {
52 : template <typename TagsList>
53 : class DataBox;
54 : } // namespace db
55 : namespace intrp {
56 : template <class Metavariables, typename InterpolationTargetTag>
57 : struct InterpolationTarget;
58 : namespace Tags {
59 : template <typename TemporalId>
60 : struct TemporalIds;
61 : } // namespace Tags
62 : } // namespace intrp
63 : namespace Tags {
64 : struct Verbosity;
65 : } // namespace Tags
66 : /// \endcond
67 :
68 : namespace ah {
69 0 : namespace OptionHolders {
70 : /// Options for finding an apparent horizon.
71 : template <typename Frame>
72 1 : struct ApparentHorizon {
73 : private:
74 0 : struct All {};
75 :
76 : public:
77 : /// See Strahlkorper for suboptions.
78 1 : struct InitialGuess {
79 0 : static constexpr Options::String help = {"Initial guess"};
80 0 : using type = ylm::Strahlkorper<Frame>;
81 : };
82 : /// See ::FastFlow for suboptions.
83 1 : struct FastFlow {
84 0 : static constexpr Options::String help = {"FastFlow options"};
85 0 : using type = ::FastFlow;
86 : };
87 0 : struct Verbosity {
88 0 : static constexpr Options::String help = {"Verbosity"};
89 0 : using type = ::Verbosity;
90 : };
91 0 : struct MaxInterpolationRetries {
92 0 : static constexpr Options::String help = {
93 : "Number of times to retry the interpolation where, with each retry, "
94 : "the two previous surfaces are averaged and that new surface is used."};
95 0 : using type = size_t;
96 : };
97 0 : struct BlocksForInterpolation {
98 0 : static constexpr Options::String help = {
99 : "Volume data will be sent to the interpolator for these block group "
100 : "names. Set to 'All' to send volume data from the entire domain."};
101 0 : using type = Options::Auto<std::vector<std::string>, All>;
102 : };
103 0 : using options = tmpl::list<InitialGuess, FastFlow, Verbosity,
104 : MaxInterpolationRetries, BlocksForInterpolation>;
105 0 : static constexpr Options::String help = {
106 : "Provide an initial guess for the apparent horizon surface\n"
107 : "(Strahlkorper) and apparent-horizon-finding-algorithm (FastFlow)\n"
108 : "options."};
109 :
110 0 : ApparentHorizon(
111 : ylm::Strahlkorper<Frame> initial_guess_in, ::FastFlow fast_flow_in,
112 : ::Verbosity verbosity_in, size_t max_interpolation_retries_in,
113 : std::optional<std::vector<std::string>> blocks_for_interpolation_in);
114 :
115 0 : ApparentHorizon() = default;
116 0 : ApparentHorizon(const ApparentHorizon& /*rhs*/) = default;
117 0 : ApparentHorizon& operator=(const ApparentHorizon& /*rhs*/) = delete;
118 0 : ApparentHorizon(ApparentHorizon&& /*rhs*/) = default;
119 0 : ApparentHorizon& operator=(ApparentHorizon&& /*rhs*/) = default;
120 0 : ~ApparentHorizon() = default;
121 :
122 : // NOLINTNEXTLINE(google-runtime-references)
123 0 : void pup(PUP::er& p);
124 :
125 0 : ylm::Strahlkorper<Frame> initial_guess{};
126 0 : ::FastFlow fast_flow{};
127 0 : ::Verbosity verbosity{::Verbosity::Quiet};
128 0 : size_t max_interpolation_retries{};
129 0 : std::optional<std::vector<std::string>> blocks_for_interpolation;
130 : };
131 :
132 : template <typename Frame>
133 0 : bool operator==(const ApparentHorizon<Frame>& lhs,
134 : const ApparentHorizon<Frame>& rhs);
135 : template <typename Frame>
136 0 : bool operator!=(const ApparentHorizon<Frame>& lhs,
137 : const ApparentHorizon<Frame>& rhs);
138 :
139 : } // namespace OptionHolders
140 :
141 0 : namespace OptionTags {
142 : template <typename InterpolationTargetTag, typename Frame>
143 0 : struct ApparentHorizon {
144 0 : using type = OptionHolders::ApparentHorizon<Frame>;
145 0 : static constexpr Options::String help{
146 : "Options for interpolation onto apparent horizon."};
147 0 : static std::string name() {
148 : return pretty_type::name<InterpolationTargetTag>();
149 : }
150 0 : using group = ApparentHorizonGroup;
151 : };
152 : } // namespace OptionTags
153 :
154 1 : namespace Tags {
155 : template <typename InterpolationTargetTag, typename Frame>
156 0 : struct ApparentHorizon : db::SimpleTag {
157 0 : using type = OptionHolders::ApparentHorizon<Frame>;
158 0 : using option_tags =
159 : tmpl::list<OptionTags::ApparentHorizon<InterpolationTargetTag, Frame>>;
160 :
161 0 : static constexpr bool pass_metavariables = false;
162 0 : static type create_from_options(const type& option) { return option; }
163 : };
164 :
165 : namespace detail {
166 : template <typename InterpolationTargetTags>
167 : struct get_horizon_options;
168 :
169 : template <typename... InterpolationTargetTags>
170 : struct get_horizon_options<tmpl::list<InterpolationTargetTags...>> {
171 : using type = tmpl::list<OptionTags::ApparentHorizon<
172 : InterpolationTargetTags,
173 : typename InterpolationTargetTags::compute_target_points::frame>...>;
174 : };
175 :
176 : CREATE_GET_TYPE_ALIAS_OR_DEFAULT(component_being_mocked)
177 : } // namespace detail
178 :
179 : /*!
180 : * \brief Holds a map between interpolation target tag name (aka a horizon) and
181 : * a set of block names that should be used for interpolation for that target.
182 : */
183 1 : struct BlocksForInterpolation : db::SimpleTag {
184 0 : using type = std::unordered_map<std::string, std::unordered_set<std::string>>;
185 : template <typename Metavariables>
186 0 : using option_tags = tmpl::push_front<
187 : typename detail::get_horizon_options<
188 : intrp::InterpolationTarget_detail::get_sequential_target_tags<
189 : Metavariables>>::type,
190 : ::domain::OptionTags::DomainCreator<Metavariables::volume_dim>>;
191 :
192 0 : static constexpr bool pass_metavariables = true;
193 : template <typename Metavariables, typename... HorizonOptions>
194 0 : static type create_from_options(
195 : const std::unique_ptr<::DomainCreator<Metavariables::volume_dim>>&
196 : domain_creator,
197 : const HorizonOptions&... all_horizon_options) {
198 : return create_from_options_impl<Metavariables>(
199 : domain_creator, std::forward_as_tuple(all_horizon_options...),
200 : std::make_index_sequence<sizeof...(HorizonOptions)>{});
201 : }
202 :
203 : private:
204 : // Need the names of the target tags which are in the option tags, but not the
205 : // horizon options themselves. This just expands a tuple to be able to index
206 : // the `option_tags` type alias so we can get the name of the target horizon
207 : template <typename Metavariables, typename HorizonOptionsTuple, size_t... Is>
208 0 : static type create_from_options_impl(
209 : const std::unique_ptr<::DomainCreator<Metavariables::volume_dim>>&
210 : domain_creator,
211 : const HorizonOptionsTuple& all_horizon_options,
212 : const std::index_sequence<Is...>& /*index_sequence*/
213 : ) {
214 : std::unordered_map<std::string, std::unordered_set<std::string>> result{};
215 :
216 : const auto block_names = domain_creator->block_names();
217 : const auto block_groups = domain_creator->block_groups();
218 :
219 : const auto append_to_result = [&](const std::string& name,
220 : const auto& horizon_options) {
221 : if (horizon_options.blocks_for_interpolation.has_value()) {
222 : result[name] = domain::expand_block_groups_to_block_names(
223 : horizon_options.blocks_for_interpolation.value(), block_names,
224 : block_groups);
225 : } else {
226 : // Insert all blocks
227 : result[name].insert(block_names.begin(), block_names.end());
228 : }
229 :
230 : // Needed for the expand_pack below
231 : return 0;
232 : };
233 :
234 : expand_pack(
235 : append_to_result(tmpl::at_c<option_tags<Metavariables>, Is + 1>::name(),
236 : std::get<Is>(all_horizon_options))...);
237 :
238 : return result;
239 : }
240 : };
241 : } // namespace Tags
242 :
243 0 : namespace TargetPoints {
244 : /// \brief Computes points on a trial apparent horizon`.
245 : ///
246 : /// This differs from `KerrHorizon` in the following ways:
247 : /// - It supplies points on a prolonged Strahlkorper, at a higher resolution
248 : /// than the Strahlkorper in the DataBox, as needed for horizon finding.
249 : /// - It uses a `FastFlow` in the DataBox.
250 : /// - It has different options (including those for `FastFlow`).
251 : ///
252 : /// Conforms to the intrp::protocols::ComputeTargetPoints protocol
253 : ///
254 : /// For requirements on InterpolationTargetTag, see
255 : /// intrp::protocols::InterpolationTargetTag
256 : template <typename InterpolationTargetTag, typename Frame>
257 1 : struct ApparentHorizon : tt::ConformsTo<intrp::protocols::ComputeTargetPoints> {
258 0 : using const_global_cache_tags =
259 : tmpl::list<Tags::BlocksForInterpolation,
260 : Tags::ApparentHorizon<InterpolationTargetTag, Frame>>;
261 0 : using is_sequential = std::true_type;
262 0 : using frame = Frame;
263 :
264 0 : using common_tags =
265 : tmpl::push_back<ylm::Tags::items_tags<Frame>, ::ah::Tags::FastFlow,
266 : logging::Tags::Verbosity<InterpolationTargetTag>,
267 : ylm::Tags::PreviousStrahlkorpers<Frame>,
268 : ::ah::Tags::PreviousIterationStrahlkorper<Frame>,
269 : ::ah::Tags::FailedInterpolationIterations>;
270 0 : using simple_tags = tmpl::append<
271 : common_tags,
272 : tmpl::conditional_t<
273 : std::is_same_v<Frame, ::Frame::Inertial>, tmpl::list<>,
274 : tmpl::list<ylm::Tags::CartesianCoords<::Frame::Inertial>>>>;
275 0 : using compute_tags =
276 : tmpl::append<typename ylm::Tags::compute_items_tags<Frame>,
277 : ylm::Tags::TimeDerivStrahlkorperCompute<Frame>>;
278 :
279 : template <typename DbTags, typename Metavariables>
280 0 : static void initialize(const gsl::not_null<db::DataBox<DbTags>*> box,
281 : const Parallel::GlobalCache<Metavariables>& cache) {
282 : const auto& options =
283 : Parallel::get<Tags::ApparentHorizon<InterpolationTargetTag, Frame>>(
284 : cache);
285 :
286 : // Put Strahlkorper and its ComputeItems, FastFlow, and verbosity
287 : // into a new DataBox.
288 : //
289 : // Note that if frame is not inertial,
290 : // ylm::Tags::Strahlkorper<::Frame::Inertial> is already
291 : // default initialized so there is no need to do anything special
292 : // here for ylm::Tags::Strahlkorper<::Frame::Inertial>.
293 : Initialization::mutate_assign<common_tags>(
294 : box, options.initial_guess, options.fast_flow, options.verbosity,
295 : std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>{},
296 : options.initial_guess, 0_st);
297 : }
298 :
299 : template <typename Metavariables, typename DbTags, typename TemporalId>
300 0 : static tnsr::I<DataVector, 3, Frame> points(
301 : db::DataBox<DbTags>& box, const tmpl::type_<Metavariables>& /*meta*/,
302 : const TemporalId& temporal_id) {
303 : const auto& fast_flow = db::get<::ah::Tags::FastFlow>(box);
304 :
305 : if (fast_flow.current_iteration() == 0) {
306 : // Put new initial guess into ylm::Tags::Strahlkorper<Frame>.
307 : // We need to do this now, and not at the end of the previous horizon
308 : // search, because only now do we know the temporal_id of this horizon
309 : // search.
310 : db::mutate<ylm::Tags::Strahlkorper<Frame>>(
311 : [&temporal_id](
312 : const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
313 : const std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>&
314 : previous_strahlkorpers) {
315 : // If we have zero previous_strahlkorpers, then the
316 : // initial guess is already in strahlkorper, so do
317 : // nothing.
318 : //
319 : // If we have one previous_strahlkorper, then we have had
320 : // a successful horizon find, and the initial guess for the
321 : // next horizon find is already in strahlkorper, so
322 : // again we do nothing.
323 : //
324 : // If we have 2 valid previous_strahlkorpers, then
325 : // we set the initial guess by linear extrapolation in time
326 : // using the last 2 previous_strahlkorpers.
327 : //
328 : // If we have 3 valid previous_strahlkorpers, then
329 : // we set the initial guess by quadratic extrapolation in time
330 : // using the last 3 previous_strahlkorpers.
331 : //
332 : // For extrapolation, we assume that
333 : // * Expansion center of all the Strahlkorpers are equal.
334 : // * Maximum L of all the Strahlkorpers are equal.
335 : // It is easy to relax the max L assumption once we start
336 : // adaptively changing the L of the strahlkorpers.
337 : if (LIKELY(previous_strahlkorpers.size() > 2)) {
338 : // Quadratic extrapolation
339 : const double new_time =
340 : intrp::InterpolationTarget_detail::get_temporal_id_value(
341 : temporal_id);
342 : const double dt_0 = previous_strahlkorpers[0].first - new_time;
343 : const double dt_1 = previous_strahlkorpers[1].first - new_time;
344 : const double dt_2 = previous_strahlkorpers[2].first - new_time;
345 : const double fac_0 =
346 : dt_1 * dt_2 / ((dt_1 - dt_0) * (dt_2 - dt_0));
347 : const double fac_1 =
348 : dt_0 * dt_2 / ((dt_2 - dt_1) * (dt_0 - dt_1));
349 : const double fac_2 = 1.0 - fac_0 - fac_1;
350 : strahlkorper->coefficients() =
351 : fac_0 * previous_strahlkorpers[0].second.coefficients() +
352 : fac_1 * previous_strahlkorpers[1].second.coefficients() +
353 : fac_2 * previous_strahlkorpers[2].second.coefficients();
354 : } else if (previous_strahlkorpers.size() > 1) {
355 : // Linear extrapolation
356 : const double new_time =
357 : intrp::InterpolationTarget_detail::get_temporal_id_value(
358 : temporal_id);
359 : const double dt_0 = previous_strahlkorpers[0].first - new_time;
360 : const double dt_1 = previous_strahlkorpers[1].first - new_time;
361 : const double fac_0 = dt_1 / (dt_1 - dt_0);
362 : const double fac_1 = 1.0 - fac_0;
363 : strahlkorper->coefficients() =
364 : fac_0 * previous_strahlkorpers[0].second.coefficients() +
365 : fac_1 * previous_strahlkorpers[1].second.coefficients();
366 : }
367 : },
368 : make_not_null(&box),
369 : db::get<ylm::Tags::PreviousStrahlkorpers<Frame>>(box));
370 : }
371 :
372 : const auto& strahlkorper = db::get<ylm::Tags::Strahlkorper<Frame>>(box);
373 :
374 : const size_t L_mesh = fast_flow.current_l_mesh(strahlkorper);
375 :
376 : const auto prolonged_strahlkorper =
377 : ylm::Strahlkorper<Frame>(L_mesh, L_mesh, strahlkorper);
378 :
379 : return ylm::cartesian_coords(prolonged_strahlkorper);
380 : }
381 : };
382 :
383 : } // namespace TargetPoints
384 : } // namespace ah
|
