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