Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cmath>
7 : #include <deque>
8 : #include <utility>
9 :
10 : #include "DataStructures/DataBox/DataBox.hpp"
11 : #include "DataStructures/VariablesTag.hpp"
12 : #include "Domain/Creators/Tags/Domain.hpp"
13 : #include "Domain/FunctionsOfTime/Tags.hpp"
14 : #include "Domain/StrahlkorperTransformations.hpp"
15 : #include "IO/Logging/Tags.hpp"
16 : #include "IO/Logging/Verbosity.hpp"
17 : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
18 : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
19 : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "Parallel/Invoke.hpp"
22 : #include "Parallel/ParallelComponentHelpers.hpp"
23 : #include "Parallel/Printf/Printf.hpp"
24 : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
25 : #include "ParallelAlgorithms/ApparentHorizonFinder/InterpolationTarget.hpp"
26 : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
27 : #include "ParallelAlgorithms/Interpolation/Actions/SendPointsToInterpolator.hpp"
28 : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
29 : #include "ParallelAlgorithms/Interpolation/Protocols/PostInterpolationCallback.hpp"
30 : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
31 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
32 : #include "Utilities/ErrorHandling/Error.hpp"
33 : #include "Utilities/Gsl.hpp"
34 : #include "Utilities/PrettyType.hpp"
35 : #include "Utilities/ProtocolHelpers.hpp"
36 : #include "Utilities/TMPL.hpp"
37 :
38 : /// \cond
39 : namespace ylm::Tags {
40 : template <typename Frame>
41 : struct Strahlkorper;
42 : } // namespace ylm::Tags
43 : namespace db {
44 : template <typename TagsList>
45 : class DataBox;
46 : } // namespace db
47 : namespace intrp {
48 : template <class Metavariables, typename InterpolationTargetTag>
49 : struct InterpolationTarget;
50 : namespace Tags {
51 : template <typename Metavariables>
52 : struct TemporalIds;
53 : } // namespace Tags
54 : } // namespace intrp
55 : /// \endcond
56 :
57 : namespace intrp::callbacks {
58 :
59 : /// \brief post interpolation callback (see
60 : /// intrp::protocols::PostInterpolationCallback) that does a FastFlow iteration
61 : /// and triggers another one until convergence.
62 : ///
63 : /// Assumes that InterpolationTargetTag contains an additional
64 : /// type alias called `post_horizon_find_callbacks`, which is a list of
65 : /// structs, each of which has a function
66 : ///
67 : /// \snippet ParallelAlgorithms/ApparentHorizonFinder/Test_ApparentHorizonFinder.cpp post_horizon_find_callback_example
68 : ///
69 : /// that is called if the FastFlow iteration has converged.
70 : /// InterpolationTargetTag also is assumed to contain an additional
71 : /// type alias called `horizon_find_failure_callbacks`, which is a list of
72 : /// structs, each of which has a function
73 : ///
74 : /// \snippet ParallelAlgorithms/ApparentHorizonFinder/Test_ApparentHorizonFinder.cpp horizon_find_failure_callbacks_example
75 : ///
76 : /// that is called if the FastFlow iteration or the interpolation has
77 : /// failed.
78 : ///
79 : /// Uses:
80 : /// - Metavariables:
81 : /// - `temporal_id`
82 : /// - DataBox:
83 : /// - `logging::Tags::Verbosity<InterpolationTargetTag>`
84 : /// - `::gr::Tags::InverseSpatialMetric<DataVector, 3, Frame>`
85 : /// - `::gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame>`
86 : /// - `::gr::Tags::SpatialChristoffelSecondKind<DataVector, 3, Frame>`
87 : /// - `::ah::Tags::FastFlow`
88 : /// - `ylm::Tags::Strahlkorper<Frame>`
89 : ///
90 : /// Modifies:
91 : /// - DataBox:
92 : /// - `::ah::Tags::FastFlow`
93 : /// - `ylm::Tags::Strahlkorper<Frame>`
94 : ///
95 : /// This can be used in InterpolationTargetTag::post_interpolation_callbacks;
96 : /// see intrp::protocols::InterpolationTargetTag for details on
97 : /// InterpolationTargetTag.
98 : ///
99 : /// ### Output
100 : ///
101 : /// Optionally, a single line of output is printed to stdout either on each
102 : /// iteration (if verbosity > Verbosity::Quiet) or on convergence
103 : /// (if verbosity > Verbosity::Silent). The output consists of the
104 : /// following labels, and values associated with each label:
105 : /// - t = time given by value of `temporal_id` argument to `apply`
106 : /// - its = current iteration number.
107 : /// - R = min and max of residual over all prolonged grid points.
108 : /// - |R| = L2 norm of residual, counting only L modes solved for.
109 : /// - |R_mesh| = L2 norm of residual over prolonged grid points.
110 : /// - r = min and max radius of trial horizon surface.
111 : ///
112 : /// #### Difference between |R| and |R_mesh|:
113 : /// The horizon is represented in a \f$Y_{lm}\f$ expansion up
114 : /// to \f$l=l_{\mathrm{surface}}\f$;
115 : /// the residual |R| represents the failure of that surface to satisfy
116 : /// the apparent horizon equation.
117 : ///
118 : /// However, at each iteration we also interpolate the horizon surface
119 : /// to a higher resolution ("prolongation"). The prolonged surface
120 : /// includes \f$Y_{lm}\f$ coefficients up to \f$l=l_{\mathrm{mesh}}\f$,
121 : /// where \f$l_{\mathrm{mesh}} > l_{\mathrm{surface}}\f$.
122 : /// The residual computed on this higher-resolution surface is |R_mesh|.
123 : ///
124 : /// As iterations proceed, |R| should decrease until it reaches
125 : /// numerical roundoff error, because we are varying all the \f$Y_{lm}\f$
126 : /// coefficients up to \f$l=l_{\mathrm{surface}}\f$
127 : /// to minimize the residual. However,
128 : /// |R_mesh| should eventually stop decreasing as iterations proceed,
129 : /// hitting a floor that represents the numerical truncation error of
130 : /// the solution.
131 : ///
132 : /// The convergence criterion looks at both |R| and |R_mesh|: Once
133 : /// |R| is small enough and |R_mesh| has stabilized, it is pointless
134 : /// to keep iterating (even though one could iterate until |R| reaches
135 : /// roundoff).
136 : ///
137 : /// Problems with convergence of the apparent horizon finder can often
138 : /// be diagnosed by looking at the behavior of |R| and |R_mesh| as a
139 : /// function of iteration.
140 : ///
141 : template <typename InterpolationTargetTag, typename Frame>
142 1 : struct FindApparentHorizon
143 : : tt::ConformsTo<intrp::protocols::PostInterpolationCallback> {
144 0 : using const_global_cache_tags =
145 : Parallel::get_const_global_cache_tags_from_actions<
146 : typename InterpolationTargetTag::post_horizon_find_callbacks>;
147 : template <typename DbTags, typename Metavariables, typename TemporalId>
148 0 : static bool apply(
149 : const gsl::not_null<db::DataBox<DbTags>*> box,
150 : const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
151 : const TemporalId& temporal_id) {
152 : bool horizon_finder_failed = false;
153 :
154 : const size_t current_iteration =
155 : get<::ah::Tags::FastFlow>(*box).current_iteration();
156 :
157 : if (current_iteration == 0) {
158 : // If we get here, we are in a new apparent horizon search, as
159 : // opposed to a subsequent iteration of the same horizon search.
160 : db::mutate<ah::Tags::PreviousIterationStrahlkorper<Frame>,
161 : ah::Tags::FailedInterpolationIterations>(
162 : [](const gsl::not_null<ylm::Strahlkorper<Frame>*>
163 : previous_fast_flow_strahlkorper,
164 : const gsl::not_null<size_t*> failed_interpolation_iterations,
165 : const ylm::Strahlkorper<Frame>& strahlkorper) {
166 : // Reset to zero. Will only be used starting from the second
167 : // iteration
168 : *failed_interpolation_iterations = 0;
169 :
170 : // First iteration the previous is just the initial guess
171 : (*previous_fast_flow_strahlkorper) = strahlkorper;
172 : },
173 : box, db::get<ylm::Tags::Strahlkorper<Frame>>(*box));
174 : }
175 :
176 : // Deal with the possibility that some of the points might be
177 : // outside of the Domain.
178 : const auto& indices_of_invalid_pts =
179 : db::get<Tags::IndicesOfInvalidInterpPoints<TemporalId>>(*box);
180 : size_t& failed_interpolation_iterations =
181 : db::get_mutable_reference<ah::Tags::FailedInterpolationIterations>(box);
182 : if (indices_of_invalid_pts.count(temporal_id) > 0 and
183 : not indices_of_invalid_pts.at(temporal_id).empty()) {
184 : ++failed_interpolation_iterations;
185 :
186 : const auto& options = Parallel::get<
187 : ah::Tags::ApparentHorizon<InterpolationTargetTag, Frame>>(*cache);
188 :
189 : // Can't recover if we've exceeded our number of attempts
190 : if (failed_interpolation_iterations <=
191 : options.max_interpolation_retries) {
192 : db::mutate<ylm::Tags::Strahlkorper<Frame>>(
193 : [&](const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
194 : const ylm::Strahlkorper<Frame>&
195 : previous_strahlkorper_iteration) {
196 : if (current_iteration == 0) {
197 : // If this is the zeroth iteration and we couldn't interpolate,
198 : // then just try increasing the size of the horizon by 50%
199 : strahlkorper->coefficients()[0] *= 1.5;
200 : } else {
201 : // Otherwise move the new trial surface halfway between the
202 : // current surface and the previous surface
203 : strahlkorper->coefficients() +=
204 : 0.5 * (previous_strahlkorper_iteration.coefficients() -
205 : strahlkorper->coefficients());
206 : }
207 : },
208 : box, db::get<ah::Tags::PreviousIterationStrahlkorper<Frame>>(*box));
209 :
210 : // Resend the new coords at the same time and iteration. Use 0 for the
211 : // array index because it's always 0 due to the targets being
212 : // singletons. No need to worry about the invalid points since those get
213 : // reset within SendPointsToInterpolator
214 : Actions::SendPointsToInterpolator<InterpolationTargetTag>::
215 : template apply<::intrp::InterpolationTarget<
216 : Metavariables, InterpolationTargetTag>>(
217 : *box, *cache, 0, temporal_id, current_iteration,
218 : failed_interpolation_iterations);
219 :
220 : // We return false because we don't want this iteration to clean
221 : // up the volume data, since we still need it.
222 : return false;
223 : } else {
224 : tmpl::for_each<
225 : typename InterpolationTargetTag::horizon_find_failure_callbacks>(
226 : [&box, &cache, &temporal_id](auto callback_v) {
227 : using callback = tmpl::type_from<decltype(callback_v)>;
228 : callback::template apply<InterpolationTargetTag>(
229 : *box, *cache, temporal_id,
230 : FastFlow::Status::InterpolationFailure);
231 : });
232 : horizon_finder_failed = true;
233 : }
234 : }
235 :
236 : // Reset now that we were able to interpolate (or if we are abandoning this
237 : // horizon find)
238 : failed_interpolation_iterations = 0;
239 :
240 : if (not horizon_finder_failed) {
241 : const auto& verbosity =
242 : db::get<logging::Tags::Verbosity<InterpolationTargetTag>>(*box);
243 : const auto& inv_g =
244 : db::get<::gr::Tags::InverseSpatialMetric<DataVector, 3, Frame>>(*box);
245 : const auto& ex_curv =
246 : db::get<::gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame>>(*box);
247 : const auto& christoffel = db::get<
248 : ::gr::Tags::SpatialChristoffelSecondKind<DataVector, 3, Frame>>(*box);
249 :
250 : std::pair<FastFlow::Status, FastFlow::IterInfo> status_and_info;
251 :
252 : // Do a FastFlow iteration.
253 : db::mutate<::ah::Tags::FastFlow, ylm::Tags::Strahlkorper<Frame>,
254 : ah::Tags::PreviousIterationStrahlkorper<Frame>>(
255 : [&inv_g, &ex_curv, &christoffel, &status_and_info](
256 : const gsl::not_null<::FastFlow*> fast_flow,
257 : const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
258 : const gsl::not_null<ylm::Strahlkorper<Frame>*>
259 : previous_fast_flow_strahlkorper) {
260 : // Set this before we iterate the fast flow. Note that this will
261 : // have to be updated once we have adaptive horizon finding
262 : (*previous_fast_flow_strahlkorper) = *strahlkorper;
263 :
264 : status_and_info = fast_flow->template iterate_horizon_finder<Frame>(
265 : strahlkorper, inv_g, ex_curv, christoffel);
266 : },
267 : box);
268 :
269 : // Determine whether we have converged, whether we need another step,
270 : // or whether we have encountered an error.
271 :
272 : const auto& status = status_and_info.first;
273 : const auto& info = status_and_info.second;
274 : const auto has_converged = converged(status);
275 :
276 : if (verbosity >= ::Verbosity::Verbose or
277 : (verbosity >= ::Verbosity::Quiet and has_converged)) {
278 : Parallel::printf(
279 : "%s: t=%.6g: its=%d: %.1e<R<%.0e, |R|=%.1g, "
280 : "|R_grid|=%.1g, %.4g<r<%.4g\n",
281 : pretty_type::name<InterpolationTargetTag>(),
282 : InterpolationTarget_detail::get_temporal_id_value(temporal_id),
283 : info.iteration, info.min_residual, info.max_residual,
284 : info.residual_ylm, info.residual_mesh, info.r_min, info.r_max);
285 : }
286 :
287 : if (status == FastFlow::Status::SuccessfulIteration) {
288 : // Do another iteration of the same horizon search.
289 : const auto& current_id =
290 : db::get<intrp::Tags::CurrentTemporalId<TemporalId>>(*box);
291 : using ::operator<<;
292 : ASSERT(current_id.has_value(),
293 : "While horizon finding, the current temporal id doesn't have a "
294 : "value, but it should.");
295 : // The iteration of these new coords is the fast flow iteration + 1
296 : // because the zeroth iteration was the initial guess. Use 0 for the
297 : // array index because it's always 0 due to the targets being singletons
298 : Actions::SendPointsToInterpolator<InterpolationTargetTag>::
299 : template apply<::intrp::InterpolationTarget<
300 : Metavariables, InterpolationTargetTag>>(
301 : *box, *cache, 0, current_id.value(), info.iteration + 1);
302 : // We return false because we don't want this iteration to clean
303 : // up the volume data, since we are using it for the next iteration
304 : // (i.e. the simple_action that we just called).
305 : return false;
306 : } else if (not has_converged) {
307 : tmpl::for_each<
308 : typename InterpolationTargetTag::horizon_find_failure_callbacks>(
309 : [&box, &cache, &temporal_id, &status](auto callback_v) {
310 : using callback = tmpl::type_from<decltype(callback_v)>;
311 : callback::template apply<InterpolationTargetTag>(
312 : *box, *cache, temporal_id, status);
313 : });
314 : horizon_finder_failed = true;
315 : }
316 : }
317 :
318 : // If it failed, don't update any variables, just reset the Strahlkorper to
319 : // it's previous value
320 : if (horizon_finder_failed) {
321 : db::mutate<ylm::Tags::Strahlkorper<Frame>>(
322 : [](const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
323 : const std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>&
324 : previous_strahlkorpers) {
325 : // Don't keep a partially-converged strahlkorper in the
326 : // DataBox. Reset to either the original initial guess or
327 : // to the last-found Strahlkorper (whichever one happens
328 : // to be in previous_strahlkorpers).
329 : *strahlkorper = previous_strahlkorpers.front().second;
330 : },
331 : box, db::get<ylm::Tags::PreviousStrahlkorpers<Frame>>(*box));
332 : } else {
333 : // The interpolated variables
334 : // Tags::Variables<InterpolationTargetTag::vars_to_interpolate_to_target>
335 : // have been interpolated from the volume to the points on the
336 : // prolonged_strahlkorper, not to the points on the actual
337 : // strahlkorper. So here we do a restriction of these
338 : // quantities onto the actual strahlkorper.
339 :
340 : // Type alias to make code more understandable.
341 : using vars_tags =
342 : typename InterpolationTargetTag::vars_to_interpolate_to_target;
343 : db::mutate_apply<
344 : tmpl::list<::Tags::Variables<vars_tags>>,
345 : tmpl::list<ylm::Tags::Strahlkorper<Frame>, ::ah::Tags::FastFlow>>(
346 : [](const gsl::not_null<Variables<vars_tags>*> vars,
347 : const ylm::Strahlkorper<Frame>& strahlkorper,
348 : const FastFlow& fast_flow) {
349 : const size_t L_mesh = fast_flow.current_l_mesh(strahlkorper);
350 : const auto prolonged_strahlkorper =
351 : ylm::Strahlkorper<Frame>(L_mesh, L_mesh, strahlkorper);
352 : auto new_vars = ::Variables<vars_tags>(
353 : strahlkorper.ylm_spherepack().physical_size());
354 :
355 : tmpl::for_each<vars_tags>([&strahlkorper, &prolonged_strahlkorper,
356 : &vars, &new_vars](auto tag_v) {
357 : using tag = typename decltype(tag_v)::type;
358 : const auto& old_var = get<tag>(*vars);
359 : auto& new_var = get<tag>(new_vars);
360 : auto old_iter = old_var.begin();
361 : auto new_iter = new_var.begin();
362 : for (; old_iter != old_var.end() and new_iter != new_var.end();
363 : ++old_iter, ++new_iter) {
364 : *new_iter = strahlkorper.ylm_spherepack().spec_to_phys(
365 : prolonged_strahlkorper.ylm_spherepack().prolong_or_restrict(
366 : prolonged_strahlkorper.ylm_spherepack().phys_to_spec(
367 : *old_iter),
368 : strahlkorper.ylm_spherepack()));
369 : }
370 : });
371 : *vars = std::move(new_vars);
372 : },
373 : box);
374 :
375 : // Compute Strahlkorper Cartesian coordinates in Inertial frame
376 : // if the current frame is not inertial.
377 : if constexpr (not std::is_same_v<Frame, ::Frame::Inertial>) {
378 : db::mutate_apply<
379 : tmpl::list<ylm::Tags::CartesianCoords<::Frame::Inertial>>,
380 : tmpl::list<ylm::Tags::Strahlkorper<Frame>,
381 : domain::Tags::Domain<Metavariables::volume_dim>>>(
382 : [&cache, &temporal_id](
383 : const gsl::not_null<tnsr::I<DataVector, 3, ::Frame::Inertial>*>
384 : inertial_strahlkorper_coords,
385 : const ylm::Strahlkorper<Frame>& strahlkorper,
386 : const Domain<Metavariables::volume_dim>& domain) {
387 : // Note that functions_of_time must already be up to
388 : // date at temporal_id because they were used in the AH
389 : // search above.
390 : const auto& functions_of_time =
391 : get<domain::Tags::FunctionsOfTime>(*cache);
392 : strahlkorper_coords_in_different_frame(
393 : inertial_strahlkorper_coords, strahlkorper, domain,
394 : functions_of_time,
395 : InterpolationTarget_detail::get_temporal_id_value(
396 : temporal_id));
397 : },
398 : box);
399 : }
400 :
401 : // Update the previous strahlkorpers. We do this before the callbacks
402 : // in case any of the callbacks need the previous strahlkorpers with the
403 : // current strahlkorper already in it.
404 : db::mutate<ylm::Tags::Strahlkorper<Frame>,
405 : ylm::Tags::PreviousStrahlkorpers<Frame>>(
406 : [&temporal_id](
407 : const gsl::not_null<ylm::Strahlkorper<Frame>*> strahlkorper,
408 : const gsl::not_null<
409 : std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>*>
410 : previous_strahlkorpers) {
411 : // This is the number of previous strahlkorpers that we
412 : // keep around.
413 : const size_t num_previous_strahlkorpers = 3;
414 :
415 : // Save a new previous_strahlkorper.
416 : previous_strahlkorpers->emplace_front(
417 : InterpolationTarget_detail::get_temporal_id_value(temporal_id),
418 : *strahlkorper);
419 :
420 : // Remove old previous_strahlkorpers that are no longer relevant.
421 : while (previous_strahlkorpers->size() >
422 : num_previous_strahlkorpers) {
423 : previous_strahlkorpers->pop_back();
424 : }
425 : },
426 : box);
427 :
428 : // Finally call callbacks
429 : tmpl::for_each<
430 : typename InterpolationTargetTag::post_horizon_find_callbacks>(
431 : [&box, &cache, &temporal_id](auto callback_v) {
432 : using callback = tmpl::type_from<decltype(callback_v)>;
433 : callback::apply(*box, *cache, temporal_id);
434 : });
435 : }
436 :
437 : // Prepare for finding horizon at a new time. Regardless of if we failed or
438 : // not, we reset fast flow.
439 : db::mutate<::ah::Tags::FastFlow>(
440 : [](const gsl::not_null<::FastFlow*> fast_flow) {
441 : fast_flow->reset_for_next_find();
442 : },
443 : box);
444 :
445 : // We return true because we are now done with all the volume data
446 : // at this temporal_id, so we want it cleaned up.
447 : return true;
448 : }
449 : };
450 : } // namespace intrp::callbacks
|
