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 <deque>
8 : #include <optional>
9 : #include <string>
10 : #include <type_traits>
11 : #include <utility>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataVector.hpp"
15 : #include "DataStructures/LinkedMessageId.hpp"
16 : #include "DataStructures/Tensor/IndexType.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "DataStructures/VariablesTag.hpp"
19 : #include "Domain/Creators/Tags/Domain.hpp"
20 : #include "Domain/FunctionsOfTime/Tags.hpp"
21 : #include "Domain/StrahlkorperTransformations.hpp"
22 : #include "IO/Logging/Verbosity.hpp"
23 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
24 : #include "NumericalAlgorithms/Strahlkorper/StrahlkorperFunctions.hpp"
25 : #include "NumericalAlgorithms/Strahlkorper/Tags.hpp"
26 : #include "Parallel/GlobalCache.hpp"
27 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
28 : #include "ParallelAlgorithms/ApparentHorizonFinder/HorizonAliases.hpp"
29 : #include "ParallelAlgorithms/ApparentHorizonFinder/Storage.hpp"
30 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
31 : #include "Utilities/Gsl.hpp"
32 : #include "Utilities/TMPL.hpp"
33 :
34 : namespace ah {
35 :
36 : namespace detail {
37 : // Update the previous surface in the mutable global cache
38 : template <typename Fr>
39 : struct UpdatePreviousSurface {
40 : static void apply(const gsl::not_null<Storage::LockedPreviousSurface<Fr>*>
41 : locked_previous_surface,
42 : Storage::PreviousSurface<Fr> previous_surface) {
43 : auto& lock = locked_previous_surface->lock;
44 : lock.write_lock();
45 : locked_previous_surface->surface = std::move(previous_surface);
46 : lock.write_unlock();
47 : }
48 : };
49 : } // namespace detail
50 :
51 : /*!
52 : * \brief Invoke the callbacks specified in the `horizon_find_callbacks` alias
53 : * of the \p HorizonMetavars.
54 : *
55 : * \details Before invoking the callbacks, this function
56 : *
57 : * 1. Restricts the final interpolated variables from the $L_\mathrm{mesh}$ used
58 : * for the FastFlow algorithm, to the actual $L$ of the Strahlkorper.
59 : * 2. Adds the current Strahlkorper to the `ah::Tags::PreviousSurfaces`.
60 : * 3. Copies the Strahlkorper, its time derivative, and the dependency that
61 : * was received with the volume data for the current time (see
62 : * `ah::Storage::SingleTimeStorage::dependency`) into the box. Also possibly
63 : * computes the Inertial coordinates of the final Strahlkorper and stores
64 : * them in the box if the `frame` of the \p HorizonMetavars isn't the
65 : * Inertial frame.
66 : */
67 : template <typename HorizonMetavars, typename DbTags, typename Metavariables>
68 1 : void invoke_callbacks(const gsl::not_null<db::DataBox<DbTags>*> box,
69 : Parallel::GlobalCache<Metavariables>& cache,
70 : const FastFlow::Status status) {
71 : using Fr = typename HorizonMetavars::frame;
72 :
73 : const auto& current_time = db::get<ah::Tags::CurrentTime>(*box).value();
74 : const auto& all_storage = db::get<ah::Tags::Storage<Fr>>(*box);
75 : const auto& current_time_storage = all_storage.at(current_time);
76 : const auto& current_iteration = current_time_storage.current_iteration;
77 : const auto& current_strahlkorper = current_iteration.strahlkorper;
78 : auto& previous_surfaces =
79 : db::get_mutable_reference<ah::Tags::PreviousSurfaces<Fr>>(box);
80 : const auto& fast_flow = db::get<ah::Tags::FastFlow>(*box);
81 :
82 : // Update the previous surfaces. We do this before the callbacks in case any
83 : // of the callbacks need the previous strahlkorpers with the current
84 : // strahlkorper already in it.
85 : previous_surfaces.emplace_front(current_time, current_strahlkorper,
86 : current_iteration.intersecting_element_ids);
87 :
88 : // Broadcast the previous surface to the mutable global cache so that elements
89 : // know whether they need to send data for the next horizon find. We do this
90 : // early so that elements get this information as soon as possible.
91 : Parallel::mutate<Tags::PreviousSurface<HorizonMetavars>,
92 : detail::UpdatePreviousSurface<Fr>>(
93 : cache, previous_surfaces.front());
94 :
95 : // Remove old previous_strahlkorpers that are no longer relevant.
96 : const size_t num_previous_strahlkorpers = 3;
97 : while (previous_surfaces.size() > num_previous_strahlkorpers) {
98 : previous_surfaces.pop_back();
99 : }
100 :
101 : // The interpolated variables have been interpolated from the volume to
102 : // the points on the prolonged_strahlkorper, not to the points on the
103 : // actual strahlkorper. So here we do a restriction of these quantities
104 : // onto the actual strahlkorper.
105 : const auto& current_ylm = current_strahlkorper.ylm_spherepack();
106 : const size_t L_mesh = fast_flow.current_l_mesh(current_strahlkorper);
107 : const auto prolonged_strahlkorper =
108 : ylm::Strahlkorper<Fr>(L_mesh, L_mesh, current_strahlkorper);
109 : const auto& prolonged_ylm = prolonged_strahlkorper.ylm_spherepack();
110 : const auto& prolonged_interpolated_vars = current_iteration.interpolated_vars;
111 : db::mutate<::Tags::Variables<ah::vars_to_interpolate_to_target<3, Fr>>>(
112 : [&](const gsl::not_null<
113 : ::Variables<ah::vars_to_interpolate_to_target<3, Fr>>*>
114 : new_interpolated_vars) {
115 : new_interpolated_vars->initialize(current_ylm.physical_size());
116 : tmpl::for_each<ah::vars_to_interpolate_to_target<3, Fr>>(
117 : [&]<typename Tag>(tmpl::type_<Tag>) {
118 : const auto& prolonged_var = get<Tag>(prolonged_interpolated_vars);
119 : auto& new_var = get<Tag>(*new_interpolated_vars);
120 : for (size_t i = 0; i < prolonged_var.size(); i++) {
121 : new_var[i] =
122 : current_ylm.spec_to_phys(prolonged_ylm.prolong_or_restrict(
123 : prolonged_ylm.phys_to_spec(prolonged_var[i]),
124 : current_ylm));
125 : }
126 : });
127 : },
128 : box);
129 :
130 : db::mutate<ylm::Tags::Strahlkorper<Fr>, ylm::Tags::TimeDerivStrahlkorper<Fr>,
131 : ah::Tags::Dependency>(
132 : [&](const gsl::not_null<ylm::Strahlkorper<Fr>*> horizon,
133 : const gsl::not_null<ylm::Strahlkorper<Fr>*> time_deriv_of_horizon,
134 : const gsl::not_null<std::optional<std::string>*>
135 : callback_dependency) {
136 : *horizon = current_strahlkorper;
137 :
138 : // This is a hack to use ylm::time_deriv_of_strahlkorper function below
139 : std::deque<std::pair<double, ylm::Strahlkorper<Fr>>>
140 : previous_horizons{};
141 : for (const auto& previous_surface : previous_surfaces) {
142 : previous_horizons.emplace_back(previous_surface.time.id,
143 : previous_surface.surface);
144 : }
145 :
146 : *time_deriv_of_horizon = current_strahlkorper;
147 : ylm::time_deriv_of_strahlkorper(time_deriv_of_horizon,
148 : previous_horizons);
149 :
150 : *callback_dependency = current_time_storage.dependency;
151 : },
152 : box);
153 :
154 : // Put inertial coords in the box if they aren't already there
155 : if constexpr (not std::is_same_v<Fr, Frame::Inertial>) {
156 : db::mutate<ylm::Tags::CartesianCoords<Frame::Inertial>>(
157 : [&](const gsl::not_null<tnsr::I<DataVector, 3>*> inertial_coords) {
158 : const auto& domain = Parallel::get<domain::Tags::Domain<3>>(cache);
159 : const auto& functions_of_time =
160 : Parallel::get<domain::Tags::FunctionsOfTime>(cache);
161 : strahlkorper_coords_in_different_frame(
162 : inertial_coords, current_strahlkorper, domain, functions_of_time,
163 : current_time.id);
164 : },
165 : box);
166 : }
167 :
168 : // Finally call callbacks
169 : tmpl::for_each<typename HorizonMetavars::horizon_find_callbacks>(
170 : [&]<typename Callback>(tmpl::type_<Callback>) {
171 : Callback::apply(*box, cache, status);
172 : });
173 : }
174 : } // namespace ah
|
