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 <optional>
8 : #include <pup.h>
9 : #include <string>
10 : #include <type_traits>
11 :
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/LinkedMessageId.hpp"
14 : #include "DataStructures/Tensor/IndexType.hpp"
15 : #include "DataStructures/Variables.hpp"
16 : #include "Domain/Creators/Tags/Domain.hpp"
17 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
18 : #include "Options/String.hpp"
19 : #include "Parallel/GlobalCache.hpp"
20 : #include "Parallel/Invoke.hpp"
21 : #include "Parallel/ParallelComponentHelpers.hpp"
22 : #include "ParallelAlgorithms/Actions/GetItemFromDistributedObject.hpp"
23 : #include "ParallelAlgorithms/ApparentHorizonFinder/Component.hpp"
24 : #include "ParallelAlgorithms/ApparentHorizonFinder/ComputeVarsToInterpolateToTarget.hpp"
25 : #include "ParallelAlgorithms/ApparentHorizonFinder/FindApparentHorizon.hpp"
26 : #include "ParallelAlgorithms/ApparentHorizonFinder/HorizonAliases.hpp"
27 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
28 : #include "ParallelAlgorithms/Events/Tags.hpp"
29 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
30 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
31 : #include "Utilities/Algorithm.hpp"
32 : #include "Utilities/CleanupRoutine.hpp"
33 : #include "Utilities/PrettyType.hpp"
34 : #include "Utilities/Serialization/CharmPupable.hpp"
35 : #include "Utilities/TMPL.hpp"
36 :
37 0 : namespace ah::Events {
38 : /*!
39 : * \brief Starts a horizon find by sending volume data (`ah::source_vars`) to
40 : * the horizon finder component (`ah::Component`) using the
41 : * `ah::FindApparentHorizon` simple action.
42 : *
43 : * \details Only sends data if this Element is in the
44 : * `ah::Tags::BlocksForHorizonFind` tag.
45 : *
46 : */
47 : template <typename HorizonMetavars>
48 1 : class FindApparentHorizon : public Event {
49 : public:
50 : /// \cond
51 : explicit FindApparentHorizon(CkMigrateMessage* /*unused*/) {}
52 : using PUP::able::register_constructor;
53 : WRAPPED_PUPable_decl_template(FindApparentHorizon); // NOLINT
54 : /// \endcond
55 :
56 0 : using options = tmpl::list<>;
57 0 : static constexpr Options::String help = "Start a horizon find.";
58 :
59 0 : static std::string name() { return pretty_type::name<HorizonMetavars>(); }
60 :
61 0 : FindApparentHorizon() = default;
62 :
63 : /// This constructor is not available through options
64 1 : explicit FindApparentHorizon(std::optional<std::string> dependency)
65 : : dependency_(std::move(dependency)) {}
66 :
67 0 : using compute_tags_for_observation_box =
68 : typename HorizonMetavars::compute_tags_on_element;
69 :
70 0 : using return_tags = tmpl::list<>;
71 0 : using argument_tags = tmpl::append<
72 : tmpl::list<typename HorizonMetavars::time_tag,
73 : ::Events::Tags::ObserverMesh<3>, domain::Tags::Element<3>>,
74 : ah::source_vars<3>>;
75 :
76 : template <typename Metavariables, typename ParallelComponent>
77 0 : void operator()(const LinkedMessageId<double>& time, const Mesh<3>& mesh,
78 : const Element<3>& element,
79 : const tnsr::aa<DataVector, 3>& spacetime_metric,
80 : const tnsr::aa<DataVector, 3>& pi,
81 : const tnsr::iaa<DataVector, 3>& phi,
82 : const tnsr::ijaa<DataVector, 3>& deriv_phi,
83 : Parallel::GlobalCache<Metavariables>& cache,
84 : const ElementId<3>& element_id,
85 : const ParallelComponent* const /*meta*/,
86 : const ObservationValue& /*observation_value*/) const {
87 : const auto& blocks_to_interpolate =
88 : Parallel::get<ah::Tags::BlocksForHorizonFind>(cache);
89 : ASSERT(blocks_to_interpolate.contains(name_),
90 : "Blocks to interpolate doesn't contain target " << name_);
91 : const auto& blocks_to_interpolate_for_this_target =
92 : blocks_to_interpolate.at(name_);
93 : const auto& domain = Parallel::get<domain::Tags::Domain<3>>(cache);
94 : const auto& blocks = domain.blocks();
95 : const auto& block = blocks[element_id.block_id()];
96 : const auto& block_name = block.name();
97 :
98 : // Only send data if this target needs this blocks data
99 : if (not blocks_to_interpolate_for_this_target.contains(block_name)) {
100 : return;
101 : }
102 :
103 : // Send volume data ONLY if this element intersected with the previous
104 : // horizon or if it's a neighbor of an intersecting element.
105 : // - WARNING: the algorithm WILL deadlock if the horizon has moved outside
106 : // of the elements that send data here. So we can't be too greedy with the
107 : // elements that send data. We may have to send data from corner neighbors
108 : // as well if deadlocks occur.
109 : // - WARNING: we don't currently wait for the `PreviousSurface` to be
110 : // updated by the last horizon find. This could be done by storing the
111 : // time of the last horizon find in the element DataBox and comparing it
112 : // to the time that's stored in `PreviousSurface`. However, this added
113 : // dependency would possibly introduce waiting. So far, I (NV) haven't
114 : // found that waiting for the update is necessary, meaning that whichever
115 : // intersecting element IDs are stored are close enough to the latest
116 : // state so that they cover the new apparent horizon (no deadlocks have
117 : // occurred in testing).
118 : // Alternatives: if we find that we need the up-to-date intersecting
119 : // elements we could just send the data without waiting if the
120 : // intersecting element IDs are not up-to-date, which just risks sending
121 : // more data than necessary (and therefore doing more work, potentially
122 : // undoing this performance optimization if it happens a lot).
123 : const auto& locked_previous_surface =
124 : Parallel::get<ah::Tags::PreviousSurface<HorizonMetavars>>(cache);
125 : {
126 : // Scope to lock and unlock the read lock
127 : locked_previous_surface.lock.read_lock();
128 : const CleanupRoutine unlock_read_lock = [&locked_previous_surface]() {
129 : locked_previous_surface.lock.read_unlock();
130 : };
131 : // Only skip elements if we have a previous surface
132 : if (locked_previous_surface.surface.has_value()) {
133 : const auto& previous_surface = locked_previous_surface.surface.value();
134 : // If we already found a horizon AFTER the current time, skip sending
135 : // anything (the data won't be needed anymore). Unlikely to occur.
136 : if (previous_surface.time.id >= time.id) {
137 : return;
138 : }
139 : // Check if this element or any of its neighbors overlaps an element
140 : // that intersected the horizon in the last find. Only send volume data
141 : // if so.
142 : const bool send_volume_data = alg::any_of(
143 : previous_surface.intersecting_element_ids,
144 : [&element_id,
145 : &element](const ElementId<3>& intersecting_element_id) {
146 : return overlapping(element_id, intersecting_element_id) or
147 : alg::any_of(
148 : element.neighbors(),
149 : [&intersecting_element_id](
150 : const std::pair<Direction<3>, Neighbors<3>>&
151 : direction_and_neighbors) {
152 : return alg::any_of(
153 : direction_and_neighbors.second.ids(),
154 : [&intersecting_element_id](
155 : const ElementId<3>& neighbor_id) {
156 : return overlapping(neighbor_id,
157 : intersecting_element_id);
158 : });
159 : });
160 : });
161 : if (not send_volume_data) {
162 : return;
163 : }
164 : } // if previous_surface.has_value()
165 : } // Unlock read lock here
166 :
167 : // Make Variables<ah::vars_to_interpolate_to_target> and
168 : // fill it by calling compute_vars_to_interpolate_to_target().
169 : // Then pass this variables to the FindApparentHorizon simple action.
170 : using horizon_frame = typename HorizonMetavars::frame;
171 : Variables<ah::vars_to_interpolate_to_target<3, horizon_frame>>
172 : vars_to_interpolate_to_target{mesh.number_of_grid_points()};
173 : if (block.is_time_dependent()) {
174 : if constexpr (Parallel::is_in_global_cache<
175 : Metavariables, domain::Tags::FunctionsOfTime>) {
176 : const auto& functions_of_time =
177 : Parallel::get<domain::Tags::FunctionsOfTime>(cache);
178 : ah::compute_vars_to_interpolate_to_target(
179 : make_not_null(&vars_to_interpolate_to_target), spacetime_metric, pi,
180 : phi, deriv_phi, time, domain, mesh, element_id, functions_of_time);
181 : } else {
182 : ERROR(
183 : "Block is time-dependent but FunctionsOfTime are not available "
184 : "in the global cache.");
185 : }
186 : } else {
187 : ah::compute_vars_to_interpolate_to_target(
188 : make_not_null(&vars_to_interpolate_to_target), spacetime_metric, pi,
189 : phi, deriv_phi, time, domain, mesh, element_id, {});
190 : }
191 :
192 : auto& horizon_finder_proxy = Parallel::get_parallel_component<
193 : ah::Component<Metavariables, HorizonMetavars>>(cache);
194 :
195 : Parallel::simple_action<ah::FindApparentHorizon<HorizonMetavars>>(
196 : horizon_finder_proxy, time, element_id, mesh,
197 : std::move(vars_to_interpolate_to_target), dependency_);
198 : }
199 :
200 0 : using is_ready_argument_tags = tmpl::list<>;
201 :
202 : template <typename Metavariables, typename ArrayIndex, typename Component>
203 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
204 : const ArrayIndex& /*array_index*/,
205 : const Component* const /*component*/) const {
206 : return true;
207 : }
208 :
209 1 : bool needs_evolved_variables() const override { return true; }
210 :
211 0 : void pup(PUP::er& p) override {
212 : Event::pup(p);
213 : p | dependency_;
214 : }
215 :
216 : private:
217 0 : std::optional<std::string> dependency_;
218 : // Evaluate the static name() function only once to avoid repeated allocations
219 0 : std::string name_ = name();
220 : };
221 :
222 : /// \cond
223 : // NOLINTBEGIN
224 : template <typename HorizonMetavars>
225 : PUP::able::PUP_ID FindApparentHorizon<HorizonMetavars>::my_PUP_ID = 0;
226 : // NOLINTEND
227 : /// \endcond
228 : } // namespace ah::Events
|
