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 <limits>
8 : #include <optional>
9 : #include <type_traits>
10 : #include <utility>
11 : #include <vector>
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 "Domain/Creators/Tags/Domain.hpp"
19 : #include "Domain/FunctionsOfTime/Tags.hpp"
20 : #include "Domain/Structure/ElementId.hpp"
21 : #include "IO/Logging/Verbosity.hpp"
22 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
23 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
24 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
25 : #include "Parallel/GlobalCache.hpp"
26 : #include "ParallelAlgorithms/Actions/FunctionsOfTimeAreReady.hpp"
27 : #include "ParallelAlgorithms/ApparentHorizonFinder/Callbacks/InvokeCallbacks.hpp"
28 : #include "ParallelAlgorithms/ApparentHorizonFinder/CleanUp.hpp"
29 : #include "ParallelAlgorithms/ApparentHorizonFinder/ComputeCoords.hpp"
30 : #include "ParallelAlgorithms/ApparentHorizonFinder/CurrentTime.hpp"
31 : #include "ParallelAlgorithms/ApparentHorizonFinder/Destination.hpp"
32 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
33 : #include "ParallelAlgorithms/ApparentHorizonFinder/HorizonAliases.hpp"
34 : #include "ParallelAlgorithms/ApparentHorizonFinder/InterpolateVolumeVars.hpp"
35 : #include "ParallelAlgorithms/ApparentHorizonFinder/OptionTags.hpp"
36 : #include "ParallelAlgorithms/ApparentHorizonFinder/Storage.hpp"
37 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
38 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
39 : #include "Utilities/Algorithm.hpp"
40 : #include "Utilities/Gsl.hpp"
41 : #include "Utilities/TMPL.hpp"
42 :
43 : namespace ah {
44 : namespace detail {
45 : template <typename Criterion>
46 : struct get_tags {
47 : using type = typename Criterion::compute_tags_for_observation_box;
48 : };
49 : } // namespace detail
50 :
51 : /*!
52 : * \brief Simple action run on the horizon finder by the Elements which receives
53 : * volume data, finds the apparent horizon, and calls the callbacks after the
54 : * horizon is found.
55 : *
56 : * \details First, the volume variables `ah::vars_to_interpolate_to_target` are
57 : * put into the `ah::Tags::Storage`. Then, we try to set the
58 : * `ah::Tags::CurrentTime` with the `ah::set_current_time` function. Once we
59 : * have a current time, we ensure the functions of time are up-to-date at this
60 : * current time with `ah::check_if_current_time_is_ready`. If they are not,
61 : * we re-trigger this action with `vars_have_already_been_received` set to
62 : * true and not sending any new volume data.
63 : *
64 : * Once we ensure the functions of time are ready, we compute the cartesian
65 : * coordinates for the current fast-flow iteration surface with
66 : * `ah::set_current_iteration_coords`. Once we have the coords, we interpolate
67 : * the volume date onto these coordinates. After that, we do one iteration of
68 : * the `FastFlow` algorithm to get a new surface. If we haven't converged yet,
69 : * we compute new cartesian and start another iteration. Once we have converged,
70 : * we call the callbacks with `ah::invoke_callbacks` and then clean up the
71 : * horizon finer with `ah::clean_up_horizon_finder`. Then we try to set a new
72 : * current time and the process starts over again. If `FastFlow` never
73 : * converges, we call the `horizon_find_failure_callbacks` from the
74 : * \p HorizonMetavars.
75 : */
76 : template <typename HorizonMetavars>
77 1 : struct FindApparentHorizon {
78 0 : using frame = typename HorizonMetavars::frame;
79 :
80 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
81 : typename ArrayIndex>
82 0 : static void apply(db::DataBox<DbTags>& box,
83 : Parallel::GlobalCache<Metavariables>& cache,
84 : const ArrayIndex& /*array_index*/,
85 : const LinkedMessageId<double>& incoming_time,
86 : const ElementId<3>& incoming_element_id,
87 : const ::Mesh<3>& incoming_mesh,
88 : Variables<ah::vars_to_interpolate_to_target<3, frame>>&&
89 : incoming_vars_to_interpolate,
90 : const std::optional<std::string>& dependency,
91 : const bool vars_have_already_been_received = false) {
92 : const auto& verbosity = db::get<ah::Tags::Verbosity>(box);
93 : const bool quiet_print = verbosity >= ::Verbosity::Quiet;
94 : const bool verbose_print = verbosity >= ::Verbosity::Verbose;
95 : const bool debug_print = verbosity >= ::Verbosity::Debug;
96 : const std::string& name = pretty_type::name<HorizonMetavars>();
97 :
98 : // If we've already completed this horizon find, ignore the volume data and
99 : // just return
100 : auto& completed_times = db::get_mutable_reference<ah::Tags::CompletedTimes>(
101 : make_not_null(&box));
102 : if (completed_times.contains(incoming_time)) {
103 : if (debug_print) {
104 : Parallel::printf(
105 : "%s: Received volume data at finished time %s from Element %s. "
106 : "Ignoring.\n",
107 : name, incoming_time, incoming_element_id);
108 : }
109 : return;
110 : }
111 :
112 : auto& current_time_optional =
113 : db::get_mutable_reference<ah::Tags::CurrentTime>(make_not_null(&box));
114 : auto& pending_times =
115 : db::get_mutable_reference<ah::Tags::PendingTimes>(make_not_null(&box));
116 :
117 : // Only add to pending if it isn't the current time.
118 : if (not(current_time_optional.has_value() and
119 : current_time_optional.value() == incoming_time)) {
120 : pending_times.insert(incoming_time);
121 : }
122 :
123 : auto& all_storage = db::get_mutable_reference<ah::Tags::Storage<frame>>(
124 : make_not_null(&box));
125 : auto& block_search_order =
126 : db::get_mutable_reference<ah::Tags::BlockSearchOrder>(
127 : make_not_null(&box));
128 :
129 : // Add volume variables and destination to the box if they haven't already
130 : // been received
131 : if (not vars_have_already_been_received) {
132 : auto& current_time_storage = all_storage[incoming_time];
133 : current_time_storage.all_volume_variables.emplace(
134 : incoming_element_id,
135 : ah::Storage::VolumeVariables<frame>{
136 : incoming_mesh, std::move(incoming_vars_to_interpolate)});
137 : current_time_storage.destination = HorizonMetavars::destination;
138 : }
139 :
140 : // ========================================================================
141 : // The incoming vars are moved and used in the above code. Below, we
142 : // only use quantities at the current time, except for forwarding the
143 : // incoming quantities to a callback for checking if the functions of time
144 : // are ready.
145 : // ========================================================================
146 :
147 : // Keep trying to find horizons for as long as we can.
148 : // In the first pass, interpolate only from the incoming element because
149 : // this was missing data. On subsequent passes, interpolate from all
150 : // elements. Notes:
151 : // - The `vars_have_already_been_received` flag is used only in
152 : // `check_if_current_time_is_ready` to re-trigger the AH find in a
153 : // callback but not sending any new data.
154 : // - If the incoming time is not the current time (but later), then we have
155 : // to do the full interpolation for the current time first before moving
156 : // on to the incoming time. In this case we also need to do the full
157 : // interpolation for the incoming time once we get there.
158 : bool interpolate_only_from_incoming_element =
159 : current_time_optional == incoming_time and
160 : not vars_have_already_been_received;
161 : while (true) {
162 : // Set current time using the pending times if it isn't already set
163 : set_current_time(make_not_null(¤t_time_optional),
164 : make_not_null(&pending_times), completed_times,
165 : all_storage, verbosity, name);
166 :
167 : // At this point, if we don't have a current id we can't do anything so
168 : // just exit.
169 : if (not current_time_optional.has_value()) {
170 : if (verbose_print) {
171 : Parallel::printf(
172 : "%s: No current time set. Stopping here. Incoming time: %s\n",
173 : name, incoming_time);
174 : }
175 : return;
176 : }
177 :
178 : const LinkedMessageId<double>& current_time =
179 : current_time_optional.value();
180 : ASSERT(all_storage.contains(current_time),
181 : "Current time doesn't have any data.");
182 :
183 : auto& current_time_storage = all_storage.at(current_time);
184 :
185 : // Only check if the current time is ready once
186 : if (UNLIKELY(not current_time_storage.time_is_ready)) {
187 : if (check_if_current_time_is_ready<HorizonMetavars>(
188 : current_time, cache, incoming_time, incoming_element_id,
189 : incoming_mesh, dependency)) {
190 : // Set this so we don't check it again
191 : current_time_storage.time_is_ready = true;
192 :
193 : if (verbose_print) {
194 : Parallel::printf("%s: Current time %s is ready.\n", name,
195 : current_time);
196 : }
197 : } else {
198 : if (verbose_print) {
199 : Parallel::printf("%s: Current time %s is NOT ready.\n", name,
200 : current_time);
201 : }
202 : // We need to wait for the functions of time to be updated so we can't
203 : // do anything yet.
204 : return;
205 : }
206 : }
207 :
208 : const auto& domain = Parallel::get<domain::Tags::Domain<3>>(cache);
209 : const auto& functions_of_time =
210 : Parallel::get<domain::Tags::FunctionsOfTime>(cache);
211 : const auto& options =
212 : Parallel::get<ah::Tags::ApparentHorizonOptions<HorizonMetavars>>(
213 : cache);
214 : auto& fast_flow =
215 : db::get_mutable_reference<ah::Tags::FastFlow>(make_not_null(&box));
216 :
217 : auto& all_volume_variables = current_time_storage.all_volume_variables;
218 : auto& current_iteration_storage = current_time_storage.current_iteration;
219 : auto& previous_iteration_surface =
220 : current_time_storage.previous_iteration_surface;
221 : auto& element_order = current_time_storage.element_order;
222 : auto& previous_surfaces =
223 : db::get_mutable_reference<ah::Tags::PreviousSurfaces<frame>>(
224 : make_not_null(&box));
225 :
226 : auto& current_resolution_l =
227 : db::get_mutable_reference<ah::Tags::CurrentResolutionL>(
228 : make_not_null(&box));
229 : // Set current resolution to the initial guess resolution if it is not
230 : // yet specified.
231 : if (UNLIKELY(not current_resolution_l.has_value())) {
232 : ASSERT(fast_flow.current_iteration() == 0,
233 : "Current resolution L is not set, but fast flow iteration > 0. "
234 : "Iteration 0 should set the current resolution L.");
235 : ASSERT(previous_surfaces.empty(),
236 : "Current resolution L is not set, but previous horizons have "
237 : "been found. Iteration 0 of each horizon find should set the "
238 : "current resolution L.");
239 : current_resolution_l = options.initial_guess.l_max();
240 : }
241 :
242 : // Keep finding the horizon until it's found with sufficient resolution
243 : bool rerunning_with_higher_resolution = false;
244 : while (true) {
245 : std::pair<FastFlow::Status, FastFlow::IterInfo> status_and_info;
246 :
247 : // Keep doing fast flow iterations for as long as we can until we find a
248 : // horizon
249 : while (true) {
250 : // Points haven't been set for this iteration so need to do so now
251 : if (not current_iteration_storage.block_coord_holders.has_value() or
252 : rerunning_with_higher_resolution) {
253 : const bool coords_set_successfully = set_current_iteration_coords(
254 : make_not_null(¤t_iteration_storage),
255 : make_not_null(&block_search_order), current_time, fast_flow,
256 : options.initial_guess, previous_iteration_surface,
257 : previous_surfaces, options.max_compute_coords_retries, domain,
258 : functions_of_time, current_resolution_l,
259 : rerunning_with_higher_resolution);
260 :
261 : // Some points are outside the domain and we cannot recover
262 : if (not coords_set_successfully) {
263 : // This will possibly throw an error, depending on which
264 : // callback is used in
265 : // HorizonMetavars::horizon_find_failure_callbacks
266 : // (ErrorOnFailedApparentHorizon or IgnoreFailedApparentHorizon)
267 : tmpl::for_each<
268 : typename HorizonMetavars::horizon_find_failure_callbacks>(
269 : [&]<typename Callback>(tmpl::type_<Callback> /*meta*/) {
270 : Callback::apply(box, cache,
271 : FastFlow::Status::InterpolationFailure);
272 : });
273 :
274 : // Still clean up in case we didn't error so we don't keep
275 : // accepting volume data for this time
276 : clean_up_horizon_finder(make_not_null(¤t_time_optional),
277 : make_not_null(&all_storage),
278 : make_not_null(&completed_times),
279 : make_not_null(&fast_flow));
280 :
281 : // If we didn't error, move on to the next horizon find
282 : break;
283 : }
284 :
285 : if (debug_print) {
286 : Parallel::printf(
287 : "%s: For current time %s, at iteration %zu, coords have been "
288 : "set.\n",
289 : name, current_time, fast_flow.current_iteration());
290 : }
291 : }
292 :
293 : // Interpolate to the current iteration points
294 : bool interpolation_is_complete = false;
295 : if (interpolate_only_from_incoming_element) {
296 : // Data is coming in one element at a time, so try to interpolate
297 : // only from the incoming element
298 : ASSERT(all_volume_variables.contains(incoming_element_id),
299 : "Trying to interpolate incoming data from element "
300 : << incoming_element_id
301 : << " but it's not available. This is a bug.");
302 : const bool interpolated_any_points = interpolate_volume_data(
303 : make_not_null(¤t_iteration_storage),
304 : all_volume_variables.at(incoming_element_id),
305 : incoming_element_id);
306 : // Store which elements we found points in for next iteration
307 : if (interpolated_any_points) {
308 : element_order.push_back(incoming_element_id);
309 : }
310 : interpolation_is_complete =
311 : current_iteration_storage.interpolation_is_complete();
312 : } else {
313 : // We already have some element data, so go through the elements in
314 : // the order we found points in before
315 : for (const auto& element_id : element_order) {
316 : ASSERT(all_volume_variables.contains(element_id),
317 : "Trying to interpolate data from element "
318 : << element_id
319 : << " that was found in a previous FastFlow iteration, "
320 : "but it's not available. This is a bug.");
321 : interpolate_volume_data(make_not_null(¤t_iteration_storage),
322 : all_volume_variables.at(element_id),
323 : element_id);
324 : interpolation_is_complete =
325 : current_iteration_storage.interpolation_is_complete();
326 : if (interpolation_is_complete) {
327 : break;
328 : }
329 : }
330 : // If we still need more points, try going through all the elements
331 : // that we have received data from so far
332 : if (not interpolation_is_complete) {
333 : for (const auto& [element_id, volume_vars_storage] :
334 : all_volume_variables) {
335 : if (std::find(element_order.begin(), element_order.end(),
336 : element_id) != element_order.end()) {
337 : // Already tried this element above
338 : continue;
339 : }
340 : const bool interpolated_any_points = interpolate_volume_data(
341 : make_not_null(¤t_iteration_storage),
342 : volume_vars_storage, element_id);
343 : if (interpolated_any_points) {
344 : element_order.push_back(element_id);
345 : }
346 : interpolation_is_complete =
347 : current_iteration_storage.interpolation_is_complete();
348 : if (interpolation_is_complete) {
349 : break;
350 : }
351 : }
352 : }
353 : }
354 :
355 : if (interpolation_is_complete) {
356 : // Next iteration try to interpolate from all elements
357 : interpolate_only_from_incoming_element = false;
358 : } else {
359 : // We need more points, return and wait for more volume data.
360 : // When this action is called again, we'll try to interpolate only
361 : // from the incoming element until we have enough data.
362 : if (debug_print) {
363 : Parallel::printf(
364 : "%s: For current time %s, at iteration %zu, need more volume "
365 : "data.\n",
366 : name, current_time, fast_flow.current_iteration());
367 : }
368 : return;
369 : }
370 :
371 : // Set this before we iterate the fast flow.
372 : previous_iteration_surface = current_iteration_storage.strahlkorper;
373 :
374 : // Do a single fast flow iteration.
375 : {
376 : using InverseSpatialMetric =
377 : gr::Tags::InverseSpatialMetric<DataVector, 3, frame>;
378 : using ExtrinsicCurvature =
379 : gr::Tags::ExtrinsicCurvature<DataVector, 3, frame>;
380 : using SpatialChristoffel =
381 : gr::Tags::SpatialChristoffelSecondKind<DataVector, 3, frame>;
382 :
383 : auto& interpolated_vars =
384 : current_iteration_storage.interpolated_vars;
385 : status_and_info = fast_flow.iterate_horizon_finder(
386 : make_not_null(¤t_iteration_storage.strahlkorper),
387 : get<InverseSpatialMetric>(interpolated_vars),
388 : get<ExtrinsicCurvature>(interpolated_vars),
389 : get<SpatialChristoffel>(interpolated_vars));
390 : }
391 :
392 : const auto& status = status_and_info.first;
393 : const auto& info = status_and_info.second;
394 : const auto has_converged = converged(status);
395 :
396 : if (verbose_print or (quiet_print and has_converged)) {
397 : Parallel::printf(
398 : "%s: t=%.6g: L=%zu: its=%zu: %.1e<R<%.0e, |R|=%.1g, "
399 : "|R_grid|=%.1g, %.4g<r<%.4g\n",
400 : name, current_time.id,
401 : all_storage.at(current_time)
402 : .current_iteration.strahlkorper.l_max(),
403 : info.iteration, info.min_residual, info.max_residual,
404 : info.residual_ylm, info.residual_mesh, info.r_min, info.r_max);
405 : }
406 :
407 : if (status == FastFlow::Status::SuccessfulIteration) {
408 : // Reset so we compute and interpolate to new points
409 : current_iteration_storage.reset_for_next_iteration();
410 :
411 : // Continue in the iteration while loop
412 : continue;
413 : } else if (not has_converged) {
414 : // This will possibly throw an error
415 : tmpl::for_each<
416 : typename HorizonMetavars::horizon_find_failure_callbacks>(
417 : [&](auto callback_v) {
418 : using callback = tmpl::type_from<decltype(callback_v)>;
419 : callback::apply(box, cache, status);
420 : });
421 :
422 : // Still clean up in case we didn't error so we don't keep accepting
423 : // volume data for this time
424 : clean_up_horizon_finder(make_not_null(¤t_time_optional),
425 : make_not_null(&all_storage),
426 : make_not_null(&completed_times),
427 : make_not_null(&fast_flow));
428 :
429 : // If we didn't error, move on to the next horizon find
430 : break;
431 : }
432 :
433 : // Break out of the iteration while loop. Go back into the
434 : // while loop for sufficient resolution
435 : break;
436 : } // while (true) [iterations]
437 : if (converged(status_and_info.first)) {
438 : std::vector<size_t> recommended_resolutions{};
439 : const auto& criteria = options.criteria;
440 :
441 : // Create an ObservationBox with the union of all compute tags needed
442 : // by criteria
443 : using compute_tags_for_obs_box =
444 : tmpl::remove_duplicates<tmpl::flatten<tmpl::transform<
445 : tmpl::at<
446 : typename Metavariables::factory_creation::factory_classes,
447 : ah::Criterion>,
448 : detail::get_tags<tmpl::_1>>>>;
449 : const auto observation_box =
450 : make_observation_box<compute_tags_for_obs_box>(
451 : make_not_null(&box));
452 :
453 : for (const auto& criterion : criteria) {
454 : recommended_resolutions.push_back(criterion->evaluate(
455 : observation_box, cache,
456 : current_time_storage.current_iteration.strahlkorper,
457 : status_and_info.second));
458 : }
459 : const size_t next_resolution_l =
460 : recommended_resolutions.empty()
461 : ? all_storage.at(current_time)
462 : .current_iteration.strahlkorper.l_max()
463 : : *std::max_element(recommended_resolutions.begin(),
464 : recommended_resolutions.end());
465 : if (next_resolution_l > *current_resolution_l) {
466 : current_resolution_l = next_resolution_l;
467 : rerunning_with_higher_resolution = true;
468 :
469 : // Prepare for trying the find again with higher resolution:
470 : // Reset fast flow, but don't "clean up," which also erases the
471 : // storage for the current time, resets the current time, and
472 : // updates the completed times. None of those things happen yet,
473 : // because the current time horizon finding is not yet finished.
474 : // Also reset the current iteration storage, since the next find
475 : // is effectively a new iteration.
476 : fast_flow.reset_for_next_find();
477 : current_time_storage.current_iteration.reset_for_next_iteration();
478 :
479 : continue;
480 : } else {
481 : rerunning_with_higher_resolution = false;
482 : }
483 :
484 : // We have converged to the apparent horizon. Invoke the callbacks and
485 : // clean up for the next horizon find
486 : invoke_callbacks<HorizonMetavars>(make_not_null(&box), cache,
487 : dependency, status_and_info.first);
488 :
489 : if (debug_print) {
490 : Parallel::printf(
491 : "%s: Horizon find finished at current time %s. Cleaning up and "
492 : "moving to next time\n",
493 : name, current_time);
494 : }
495 :
496 : current_resolution_l = next_resolution_l;
497 : clean_up_horizon_finder(make_not_null(¤t_time_optional),
498 : make_not_null(&all_storage),
499 : make_not_null(&completed_times),
500 : make_not_null(&fast_flow));
501 : }
502 :
503 : // We have either i) converged with sufficient resolution, or
504 : // ii) failed to converge but not thrown an error. In both cases,
505 : // continue to the next time by breaking out of the resolution
506 : // while loop, going back into the overall while loop for time.
507 : break;
508 : } // while (true) [sufficient resolution]
509 : } // while (true) [time]
510 : }
511 : };
512 : } // namespace ah
|
