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 <set>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <unordered_set>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataVector.hpp"
16 : #include "DataStructures/LinkedMessageId.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "Domain/BlockLogicalCoordinates.hpp"
20 : #include "Domain/Structure/ElementId.hpp"
21 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
22 : #include "NumericalAlgorithms/Strahlkorper/Strahlkorper.hpp"
23 : #include "Parallel/MultiReaderSpinlock.hpp"
24 : #include "ParallelAlgorithms/ApparentHorizonFinder/Destination.hpp"
25 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
26 : #include "ParallelAlgorithms/ApparentHorizonFinder/HorizonAliases.hpp"
27 :
28 0 : namespace ah::Storage {
29 : /*!
30 : * \brief Holds the `ah::source_vars`, mesh, and other variables on the horizon
31 : * finder for a given element.
32 : */
33 : template <typename Fr>
34 1 : struct VolumeVariables {
35 : /*!
36 : * \brief The mesh that corresponds to the volume vars.
37 : */
38 1 : Mesh<3> mesh;
39 :
40 : /*!
41 : * \brief A `Variables` of the tensors in the volume that we need to
42 : * interpolate onto the horizon.
43 : */
44 : Variables<ah::vars_to_interpolate_to_target<3, Fr>>
45 1 : vars_to_interpolate_to_target{};
46 :
47 : // NOLINTNEXTLINE(google-runtime-references)
48 0 : void pup(PUP::er& p);
49 : };
50 :
51 : template <typename Fr>
52 0 : bool operator==(const VolumeVariables<Fr>& lhs, const VolumeVariables<Fr>& rhs);
53 : template <typename Fr>
54 0 : bool operator!=(const VolumeVariables<Fr>& lhs, const VolumeVariables<Fr>& rhs);
55 :
56 : /*!
57 : * \brief Holds the `ylm::Strahlkorper` and associated quantities for a single
58 : * `FastFlow` iteration.
59 : */
60 : template <typename Fr>
61 1 : struct Iteration {
62 : /*!
63 : * \brief The surface for this iteration
64 : */
65 1 : ylm::Strahlkorper<Fr> strahlkorper;
66 : /*!
67 : * \brief Holds the list of all points (in block logical coordinates) that
68 : * need to be interpolated onto.
69 : *
70 : * \details If an element of the vector is `std::nullopt`, then that point is
71 : * outside the domain.
72 : */
73 1 : std::optional<std::vector<BlockLogicalCoords<3>>> block_coord_holders;
74 : /*!
75 : * \brief Holds the interpolated `Variables` on the points in
76 : * `block_coord_holders`.
77 : *
78 : * \details The grid points inside are indexed according to
79 : * `block_coord_holders`.
80 : */
81 1 : Variables<ah::vars_to_interpolate_to_target<3, Fr>> interpolated_vars{};
82 : /*!
83 : * \brief Keeps track of the indices in `interpolated_vars` that have
84 : * already been interpolated to.
85 : */
86 1 : std::vector<bool> indices_interpolated_to_thus_far{};
87 : /*!
88 : * \brief Holds the element IDs of all elements that intersect with the
89 : * current iteration surface. Used to determine which elements will send
90 : * data for the next horizon find.
91 : */
92 1 : std::unordered_set<ElementId<3>> intersecting_element_ids{};
93 : /*!
94 : * \brief Offsets of newly interpolated points in the overall tensor (used as
95 : * memory buffer)
96 : */
97 1 : std::vector<size_t> offsets_of_newly_interpolated_points{};
98 : /*!
99 : * \brief Logical coordinates of newly interpolated points (used as memory
100 : * buffer)
101 : *
102 : * These `std::vector`s are used to reserve memory and then append points to
103 : * them as we find them in an element. The memory is reused for each element.
104 : * Then, a non-owning DataVector is created by pointing into this memory.
105 : * That's why this is a `std::array` of `std::vector`s, not vice versa.
106 : */
107 : std::array<std::vector<double>, 3>
108 1 : x_element_logical_of_newly_interpolated_points{};
109 : /*!
110 : * \brief Buffer for newly interpolated variables (used as memory buffer)
111 : */
112 1 : std::vector<double> newly_interpolated_vars_buffer{};
113 :
114 : /*!
115 : * \brief How many times we've tried to compute the coordinates for this
116 : * iteration.
117 : */
118 1 : size_t compute_coords_retries = 0;
119 :
120 : /*!
121 : * \brief Whether all points in `interpolated_vars` have been filled.
122 : */
123 1 : bool interpolation_is_complete() const;
124 :
125 0 : void reset_for_next_iteration();
126 :
127 : // NOLINTNEXTLINE(google-runtime-references)
128 0 : void pup(PUP::er& p);
129 : };
130 :
131 : template <typename Fr>
132 0 : bool operator==(const Iteration<Fr>& lhs, const Iteration<Fr>& rhs);
133 : template <typename Fr>
134 0 : bool operator!=(const Iteration<Fr>& lhs, const Iteration<Fr>& rhs);
135 :
136 : /*!
137 : * \brief Holds all data necessary for a single horizon find.
138 : *
139 : * \details This includes volume variables which persist for the entire horizon
140 : * find, and also interpolated variables that are updated for each iteration.
141 : */
142 : template <typename Fr>
143 1 : struct SingleTimeStorage {
144 : /*!
145 : * \brief Map between `ElementId`s and the volume variables from that element.
146 : */
147 1 : std::unordered_map<ElementId<3>, VolumeVariables<Fr>> all_volume_variables;
148 : /*!
149 : * \brief Elements in which we have found points to interpolate to in previous
150 : * iterations, to try first before searching all elements.
151 : *
152 : * This is only a performance optimization. The interpolated values don't
153 : * depend on the order because points that lie directly on shared element
154 : * boundaries are always assigned to the same element (see
155 : * `element_contains`). Otherwise, the interpolated values at these points
156 : * would fluctuate with the order in which element data arrives, which can
157 : * prevent convergence (see
158 : * https://github.com/sxs-collaboration/spectre/issues/3899) and makes
159 : * simulations nondeterministic.
160 : */
161 1 : std::vector<ElementId<3>> element_order{};
162 : /*!
163 : * \brief The `Iteration` data for the current fast flow iteration.
164 : */
165 1 : Iteration<Fr> current_iteration;
166 : /*!
167 : * \brief The previous iteration surface, used if interpolation fails.
168 : */
169 1 : ylm::Strahlkorper<Fr> previous_iteration_surface;
170 : /*!
171 : * \brief The `ah::Destination` for this horizon find.
172 : */
173 1 : Destination destination{};
174 : /*!
175 : * \brief The dependency that was sent along with the volume data for this
176 : * horizon find, if any.
177 : *
178 : * \details This is copied into `ah::Tags::Dependency` right before the
179 : * callbacks for this horizon find are invoked, on success as well as on
180 : * failure, so that e.g. `ah::callbacks::SendDependencyToObserverWriter` can
181 : * tell the ObserverWriter whether to write or drop the volume data it holds
182 : * for this time.
183 : */
184 1 : std::optional<std::string> dependency = std::nullopt;
185 : /*!
186 : * \brief Whether we have checked if the functions of time are up to date for
187 : * this horizon find.
188 : */
189 1 : bool time_is_ready = false;
190 :
191 : // NOLINTNEXTLINE(google-runtime-references)
192 0 : void pup(PUP::er& p);
193 : };
194 :
195 : template <typename Fr>
196 0 : bool operator==(const SingleTimeStorage<Fr>& lhs,
197 : const SingleTimeStorage<Fr>& rhs);
198 : template <typename Fr>
199 0 : bool operator!=(const SingleTimeStorage<Fr>& lhs,
200 : const SingleTimeStorage<Fr>& rhs);
201 :
202 : /*!
203 : * \brief The time and final surface for a previous horizon find.
204 : */
205 : template <typename Fr>
206 1 : struct PreviousSurface {
207 0 : PreviousSurface() = default;
208 0 : PreviousSurface(const LinkedMessageId<double>& time_in,
209 : ylm::Strahlkorper<Fr> surface_in,
210 : std::unordered_set<ElementId<3>> intersecting_element_ids_in);
211 :
212 0 : LinkedMessageId<double> time;
213 0 : ylm::Strahlkorper<Fr> surface;
214 0 : std::unordered_set<ElementId<3>> intersecting_element_ids;
215 :
216 : // NOLINTNEXTLINE(google-runtime-references)
217 0 : void pup(PUP::er& p);
218 : };
219 :
220 : template <typename Fr>
221 0 : bool operator==(const PreviousSurface<Fr>& lhs, const PreviousSurface<Fr>& rhs);
222 : template <typename Fr>
223 0 : bool operator!=(const PreviousSurface<Fr>& lhs, const PreviousSurface<Fr>& rhs);
224 :
225 : /*!
226 : * \brief Holds a previous surface and a lock that protects it.
227 : *
228 : * \details This is used to store and update a previous surface in the global
229 : * cache and allow multiple readers (elements) to access it simultaneously.
230 : */
231 : template <typename Fr>
232 1 : struct LockedPreviousSurface {
233 0 : std::optional<PreviousSurface<Fr>> surface;
234 : // Lock is mutable so it can be retrieved from the const global cache and put
235 : // in read-lock mode by elements.
236 : // NOLINTNEXTLINE(spectre-mutable)
237 0 : mutable Parallel::MultiReaderSpinlock lock;
238 :
239 0 : LockedPreviousSurface();
240 0 : explicit LockedPreviousSurface(const PreviousSurface<Fr>& rhs);
241 0 : LockedPreviousSurface(const LockedPreviousSurface& rhs);
242 0 : LockedPreviousSurface& operator=(const LockedPreviousSurface& rhs);
243 0 : LockedPreviousSurface(LockedPreviousSurface&& rhs);
244 0 : LockedPreviousSurface& operator=(LockedPreviousSurface&& rhs);
245 0 : ~LockedPreviousSurface() = default;
246 :
247 : // NOLINTNEXTLINE(google-runtime-references)
248 0 : void pup(PUP::er& p);
249 : };
250 :
251 : template <typename Fr>
252 0 : bool operator==(const LockedPreviousSurface<Fr>& lhs,
253 : const LockedPreviousSurface<Fr>& rhs);
254 : template <typename Fr>
255 0 : bool operator!=(const LockedPreviousSurface<Fr>& lhs,
256 : const LockedPreviousSurface<Fr>& rhs);
257 :
258 : } // namespace ah::Storage
|
