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 <cstddef>
8 : #include <cstdint>
9 : #include <limits>
10 : #include <memory>
11 : #include <vector>
12 :
13 : #include "DataStructures/ComplexDataVector.hpp"
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataVector.hpp"
16 : #include "DataStructures/SpinWeighted.hpp"
17 : #include "DataStructures/Tags.hpp"
18 : #include "DataStructures/Tensor/TypeAliases.hpp"
19 : #include "DataStructures/Variables.hpp"
20 : #include "Evolution/Systems/Cce/GaugeTransformBoundaryData.hpp"
21 : #include "Evolution/Systems/Cce/Tags.hpp"
22 : #include "NumericalAlgorithms/Interpolation/SpanInterpolator.hpp"
23 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCollocation.hpp"
24 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshDerivatives.hpp"
25 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshInterpolation.hpp"
26 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshTags.hpp"
27 : #include "Options/String.hpp"
28 : #include "Parallel/NodeLock.hpp"
29 : #include "Utilities/CallWithDynamicType.hpp"
30 : #include "Utilities/Gsl.hpp"
31 : #include "Utilities/Serialization/CharmPupable.hpp"
32 : #include "Utilities/TMPL.hpp"
33 :
34 : /// \cond
35 : namespace Cce::Solutions::LinearizedBondiSachs_detail::InitializeJ {
36 : struct LinearizedBondiSachs;
37 : } // namespace Cce::Solutions::LinearizedBondiSachs_detail::InitializeJ
38 : /// \endcond
39 :
40 : namespace Cce {
41 : /// Contains utilities and \ref DataBoxGroup mutators for generating data for
42 : /// \f$J\f$ on the initial CCE hypersurface.
43 : namespace InitializeJ {
44 :
45 : namespace detail {
46 : // used to provide a default for the finalize functor in
47 : // `iteratively_adapt_angular_coordinates`
48 : struct NoOpFinalize {
49 : void operator()(
50 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& /*gauge_c*/,
51 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& /*gauge_d*/,
52 : const tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>&
53 : /*angular_cauchy_coordinates*/,
54 : const Spectral::Swsh::SwshInterpolator& /*interpolator*/) const {}
55 : };
56 :
57 : // eth of the Cartesian coordinates x^i on the identity grid. The angular solves
58 : // interpolate these onto the current map at each iteration, rather than
59 : // differentiating the displaced coordinates, which aliases and raises the
60 : // residual floor.
61 : void unit_sphere_eth_cartesian_coordinates(
62 : gsl::not_null<SpinWeighted<ComplexDataVector, 1>*> eth_x,
63 : gsl::not_null<SpinWeighted<ComplexDataVector, 1>*> eth_y,
64 : gsl::not_null<SpinWeighted<ComplexDataVector, 1>*> eth_z, size_t l_max);
65 :
66 : // Normalize the Cartesian coordinates of the current map, and update its
67 : // angular coordinates, interpolator and Jacobians to match.
68 : void update_angular_coordinates_and_jacobians(
69 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
70 : gsl::not_null<
71 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
72 : angular_cauchy_coordinates,
73 : gsl::not_null<Spectral::Swsh::SwshInterpolator*> interpolator,
74 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> gauge_c,
75 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> gauge_d,
76 : size_t l_max);
77 :
78 : // What an angular solve does when it stops above its tolerance.
79 : // `Uninitialized` catches a default-constructed value and is an error to pass.
80 : enum class UnconvergedAngularSolve : uint8_t {
81 : Uninitialized,
82 : Error,
83 : Warn,
84 : Silent
85 : };
86 :
87 : // The residual of the map an angular solve hands back, and the number of
88 : // coordinate updates it performed.
89 : struct AngularSolveResult {
90 : double max_error = std::numeric_limits<double>::signaling_NaN();
91 : size_t number_of_coordinate_updates = 0;
92 : };
93 :
94 : // Abort if `max_error` exceeds `error_threshold` or is not a number.
95 : void check_angular_solve_error_threshold(double max_error,
96 : double error_threshold);
97 :
98 : // Report an angular solve that stopped above its tolerance as an error, a
99 : // warning or not at all. The messages describe `tolerance` as a target for
100 : // `max_error`, the role it has in the default convergence test.
101 : void report_unconverged_angular_solve(UnconvergedAngularSolve if_unconverged,
102 : double tolerance,
103 : size_t number_of_iterations,
104 : double max_error);
105 :
106 : // used to provide a default for the convergence test in
107 : // `iteratively_adapt_angular_coordinates`
108 : struct ResidualBelowTolerance {
109 : bool operator()(const double max_error, const double tolerance) const {
110 : return max_error < tolerance;
111 : }
112 : };
113 :
114 : // Records the running minimum of the residuals of an angular solve, and tells
115 : // when it has improved by less than a relative `tolerance` over the last
116 : // `window` iterations. Unlike the change between iterations, this is not held
117 : // open by round-off noise on a slow descent.
118 : class ResidualPlateau {
119 : public:
120 : explicit ResidualPlateau(size_t window);
121 :
122 : void record(double residual);
123 :
124 : // False until `window + 1` residuals are recorded.
125 : bool stopped_improving(double tolerance) const;
126 :
127 : private:
128 : size_t window_;
129 : std::vector<double> smallest_residual_;
130 : };
131 :
132 : // perform an iterative solve for the set of angular coordinates. The iteration
133 : // callable `iteration_function` must have function signature:
134 : //
135 : // double iteration_function(
136 : // const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
137 : // gauge_c_step,
138 : // const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
139 : // gauge_d_step,
140 : // const Scalar<SpinWeighted<ComplexDataVector, 2>>& gauge_c,
141 : // const Scalar<SpinWeighted<ComplexDataVector, 0>>& gauge_d,
142 : // const Spectral::Swsh::SwshInterpolator& iteration_interpolator);
143 : //
144 : // but need not be a function pointer -- a callable class or lambda will also
145 : // suffice.
146 : // For each step specified by the iteration function, the coordinates are
147 : // updated via \hat \eth \delta x^i = \delta c \eth x^i|_{x^i=\hat x^i}
148 : // + \delta \bar d (\eth x^i)|_{x^i=\hat x^i}
149 : // This coordinate update is exact, and comes from expanding the chain rule to
150 : // determine Jacobian factors. However, the result is not guaranteed to
151 : // produce the desired Jacobian c and d, because \delta c and \delta d are
152 : // not necessarily consistent with the underlying coordinates.
153 : // We then update the x^i by inverting \hat \eth, which is also exact, but
154 : // assumes a no l=0 contribution to the transformation.
155 : // Depending on the choice of approximations used to specify
156 : // `iteration_function`, though, the method can be slow to converge.
157 :
158 : // However, the iterations are typically fast, and the computation is for
159 : // initial data that needs to be computed only once during a simulation, so it
160 : // is not currently an optimization priority. If this function becomes a
161 : // bottleneck, the numerical procedure of the iterative method or the choice of
162 : // approximation used for `iteration_function` should be revisited.
163 : //
164 : // The solve stops once `has_converged(max_error, tolerance)` holds (by default
165 : // `max_error < tolerance`) or after `max_steps` steps, and hands back the best
166 : // map it evaluated. `has_converged` is called once after each evaluation of
167 : // `iteration_function` in the loop, but not after the final re-evaluation of
168 : // the best map. A caller with another convergence test passes
169 : // `UnconvergedAngularSolve::Silent` and reports a failure itself, since the
170 : // other modes describe `tolerance` as a target for `max_error`. With
171 : // `initialize_coordinates = false` it continues from the map the caller
172 : // supplies.
173 : template <typename IterationFunctor, typename FinalizeFunctor = NoOpFinalize,
174 : typename ConvergenceTest = ResidualBelowTolerance>
175 : AngularSolveResult iteratively_adapt_angular_coordinates(
176 : const gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
177 : const gsl::not_null<
178 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
179 : angular_cauchy_coordinates,
180 : const size_t l_max, const double tolerance, const size_t max_steps,
181 : const double error_threshold, const IterationFunctor& iteration_function,
182 : const UnconvergedAngularSolve if_unconverged,
183 : const FinalizeFunctor finalize_function = NoOpFinalize{},
184 : const bool initialize_coordinates = true,
185 : const ConvergenceTest& has_converged = ResidualBelowTolerance{}) {
186 : const size_t number_of_angular_points =
187 : Spectral::Swsh::number_of_swsh_collocation_points(l_max);
188 :
189 : if (initialize_coordinates) {
190 : Spectral::Swsh::create_angular_and_cartesian_coordinates(
191 : cartesian_cauchy_coordinates, angular_cauchy_coordinates, l_max);
192 : }
193 :
194 : Variables<tmpl::list<
195 : // eth of cartesian coordinates
196 : ::Tags::TempSpinWeightedScalar<3, 1>,
197 : ::Tags::TempSpinWeightedScalar<4, 1>,
198 : ::Tags::TempSpinWeightedScalar<5, 1>,
199 : // eth of gauge-transformed cartesian coordinates
200 : ::Tags::TempSpinWeightedScalar<6, 1>,
201 : ::Tags::TempSpinWeightedScalar<7, 1>,
202 : ::Tags::TempSpinWeightedScalar<8, 1>,
203 : // gauge Jacobians
204 : ::Tags::TempSpinWeightedScalar<9, 2>,
205 : ::Tags::TempSpinWeightedScalar<10, 0>,
206 : // gauge Jacobians on next iteration
207 : ::Tags::TempSpinWeightedScalar<11, 2>,
208 : ::Tags::TempSpinWeightedScalar<12, 0>,
209 : // cartesian coordinates steps
210 : ::Tags::TempSpinWeightedScalar<13, 0>,
211 : ::Tags::TempSpinWeightedScalar<14, 0>,
212 : ::Tags::TempSpinWeightedScalar<15, 0>>>
213 : computation_buffers{number_of_angular_points};
214 :
215 : auto& eth_x =
216 : get(get<::Tags::TempSpinWeightedScalar<3, 1>>(computation_buffers));
217 : auto& eth_y =
218 : get(get<::Tags::TempSpinWeightedScalar<4, 1>>(computation_buffers));
219 : auto& eth_z =
220 : get(get<::Tags::TempSpinWeightedScalar<5, 1>>(computation_buffers));
221 : unit_sphere_eth_cartesian_coordinates(make_not_null(ð_x),
222 : make_not_null(ð_y),
223 : make_not_null(ð_z), l_max);
224 :
225 : auto& evolution_gauge_eth_x_step =
226 : get(get<::Tags::TempSpinWeightedScalar<6, 1>>(computation_buffers));
227 : auto& evolution_gauge_eth_y_step =
228 : get(get<::Tags::TempSpinWeightedScalar<7, 1>>(computation_buffers));
229 : auto& evolution_gauge_eth_z_step =
230 : get(get<::Tags::TempSpinWeightedScalar<8, 1>>(computation_buffers));
231 :
232 : auto& gauge_c =
233 : get<::Tags::TempSpinWeightedScalar<9, 2>>(computation_buffers);
234 : auto& gauge_d =
235 : get<::Tags::TempSpinWeightedScalar<10, 0>>(computation_buffers);
236 :
237 : auto& gauge_c_step =
238 : get<::Tags::TempSpinWeightedScalar<11, 2>>(computation_buffers);
239 : auto& gauge_d_step =
240 : get<::Tags::TempSpinWeightedScalar<12, 0>>(computation_buffers);
241 :
242 : auto& x_step =
243 : get(get<::Tags::TempSpinWeightedScalar<13, 0>>(computation_buffers));
244 : auto& y_step =
245 : get(get<::Tags::TempSpinWeightedScalar<14, 0>>(computation_buffers));
246 : auto& z_step =
247 : get(get<::Tags::TempSpinWeightedScalar<15, 0>>(computation_buffers));
248 :
249 : double max_error = 1.0;
250 : size_t number_of_steps = 0;
251 : bool converged = false;
252 : Spectral::Swsh::SwshInterpolator iteration_interpolator;
253 :
254 : const auto evaluate_current_map = [&]() {
255 : update_angular_coordinates_and_jacobians(
256 : cartesian_cauchy_coordinates, angular_cauchy_coordinates,
257 : make_not_null(&iteration_interpolator), make_not_null(&gauge_c),
258 : make_not_null(&gauge_d), l_max);
259 : return iteration_function(make_not_null(&gauge_c_step),
260 : make_not_null(&gauge_d_step), gauge_c, gauge_d,
261 : iteration_interpolator);
262 : };
263 :
264 : // The best map is copied before `evaluate_current_map` normalizes it, so that
265 : // evaluating it again on the rewind reproduces its residual bit for bit.
266 : auto unevaluated_cartesian_cauchy_coordinates = *cartesian_cauchy_coordinates;
267 : auto best_cartesian_cauchy_coordinates = *cartesian_cauchy_coordinates;
268 : double best_max_error = std::numeric_limits<double>::max();
269 :
270 : while (true) {
271 : unevaluated_cartesian_cauchy_coordinates = *cartesian_cauchy_coordinates;
272 : max_error = evaluate_current_map();
273 : if (max_error < best_max_error) {
274 : best_max_error = max_error;
275 : best_cartesian_cauchy_coordinates =
276 : unevaluated_cartesian_cauchy_coordinates;
277 : }
278 :
279 : check_angular_solve_error_threshold(max_error, error_threshold);
280 : ++number_of_steps;
281 : converged = has_converged(max_error, tolerance);
282 : if (converged or number_of_steps > max_steps) {
283 : break;
284 : }
285 : // using the evolution_gauge_.._step as temporary buffers for the
286 : // interpolation results
287 : iteration_interpolator.interpolate(
288 : make_not_null(&evolution_gauge_eth_x_step), eth_x);
289 : iteration_interpolator.interpolate(
290 : make_not_null(&evolution_gauge_eth_y_step), eth_y);
291 : iteration_interpolator.interpolate(
292 : make_not_null(&evolution_gauge_eth_z_step), eth_z);
293 :
294 : evolution_gauge_eth_x_step =
295 : 0.5 * ((get(gauge_c_step)) * conj(evolution_gauge_eth_x_step) +
296 : conj((get(gauge_d_step))) * evolution_gauge_eth_x_step);
297 : evolution_gauge_eth_y_step =
298 : 0.5 * ((get(gauge_c_step)) * conj(evolution_gauge_eth_y_step) +
299 : conj((get(gauge_d_step))) * evolution_gauge_eth_y_step);
300 : evolution_gauge_eth_z_step =
301 : 0.5 * ((get(gauge_c_step)) * conj(evolution_gauge_eth_z_step) +
302 : conj((get(gauge_d_step))) * evolution_gauge_eth_z_step);
303 :
304 : Spectral::Swsh::angular_derivatives<tmpl::list<
305 : Spectral::Swsh::Tags::InverseEth, Spectral::Swsh::Tags::InverseEth,
306 : Spectral::Swsh::Tags::InverseEth>>(
307 : l_max, 1, make_not_null(&x_step), make_not_null(&y_step),
308 : make_not_null(&z_step), evolution_gauge_eth_x_step,
309 : evolution_gauge_eth_y_step, evolution_gauge_eth_z_step);
310 :
311 : get<0>(*cartesian_cauchy_coordinates) += real(x_step.data());
312 : get<1>(*cartesian_cauchy_coordinates) += real(y_step.data());
313 : get<2>(*cartesian_cauchy_coordinates) += real(z_step.data());
314 : }
315 :
316 : if (best_max_error < max_error) {
317 : *cartesian_cauchy_coordinates = best_cartesian_cauchy_coordinates;
318 : max_error = evaluate_current_map();
319 : }
320 :
321 : finalize_function(gauge_c, gauge_d, *angular_cauchy_coordinates,
322 : iteration_interpolator);
323 :
324 : // The loop evaluates once more than it updates.
325 : const size_t number_of_coordinate_updates = number_of_steps - 1;
326 : if (not converged) {
327 : report_unconverged_angular_solve(if_unconverged, tolerance,
328 : number_of_coordinate_updates, max_error);
329 : }
330 : return {.max_error = max_error,
331 : .number_of_coordinate_updates = number_of_coordinate_updates};
332 : }
333 :
334 : // Solve for the partially flat angular coordinates with a displacement
335 : // generated by a single spin-weight-1 potential \eta,
336 : //
337 : // \delta \hat x^i = Re(\eta conj(\eth x^i)|_{\hat x}),
338 : //
339 : // where \eth x^i, the eth of the Cartesian coordinates on the identity grid, is
340 : // evaluated at the current map \hat x.
341 : // `iteratively_adapt_angular_coordinates` prescribes both Jacobian variations,
342 : // four real functions for a map that has two, so only part of each requested
343 : // step is realized. Here the variation of \hat c at the identity is \eth \eta,
344 : // so the step to the target Jacobian \hat c_* is one application of \eth^{-1},
345 : //
346 : // \eta = \eth^{-1}(\hat c_* - \hat c),
347 : //
348 : // and \hat d follows from the map. \eth^{-1} leaves the l = 1 part of \eta at
349 : // zero: \eth annihilates it, and it is the residual conformal freedom of the
350 : // problem. Since \hat c_* is the exact root of the gauge condition, a few
351 : // iterations do the work of many linearized sweeps.
352 : //
353 : // `target_function` must have the signature
354 : //
355 : // double target_function(
356 : // const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
357 : // gauge_c_target,
358 : // const Scalar<SpinWeighted<ComplexDataVector, 2>>& gauge_c,
359 : // const Scalar<SpinWeighted<ComplexDataVector, 0>>& gauge_d,
360 : // const Spectral::Swsh::SwshInterpolator& iteration_interpolator);
361 : //
362 : // setting `gauge_c_target` to \hat c_* and returning the current residual.
363 : //
364 : // The residual first falls by orders of magnitude and then creeps back up, so
365 : // the solve stops at the first iteration from the third on with
366 : // e_n > 0.9 e_{n-1}, and hands back the best map it evaluated. The returned
367 : // `number_of_coordinate_updates` counts every update the solve performed,
368 : // including those discarded by returning to the best map, so that it measures
369 : // the work done. `CauchySecondOrder` subtracts it from `J0MaxIterations` to
370 : // get the number of linearized sweeps still allowed.
371 : template <typename TargetFunctor, typename FinalizeFunctor = NoOpFinalize>
372 : AngularSolveResult adapt_angular_coordinates_via_potential(
373 : const gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
374 : const gsl::not_null<
375 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
376 : angular_cauchy_coordinates,
377 : const size_t l_max, const double tolerance, const size_t max_steps,
378 : const double error_threshold, const TargetFunctor& target_function,
379 : const UnconvergedAngularSolve if_unconverged,
380 : const FinalizeFunctor finalize_function = NoOpFinalize{},
381 : const bool initialize_coordinates = true) {
382 : constexpr double plateau_factor = 0.9;
383 : const size_t number_of_angular_points =
384 : Spectral::Swsh::number_of_swsh_collocation_points(l_max);
385 :
386 : if (initialize_coordinates) {
387 : Spectral::Swsh::create_angular_and_cartesian_coordinates(
388 : cartesian_cauchy_coordinates, angular_cauchy_coordinates, l_max);
389 : }
390 :
391 : Variables<tmpl::list<
392 : // eth of the cartesian coordinates
393 : ::Tags::TempSpinWeightedScalar<0, 1>,
394 : ::Tags::TempSpinWeightedScalar<1, 1>,
395 : ::Tags::TempSpinWeightedScalar<2, 1>,
396 : // ... interpolated onto the current map
397 : ::Tags::TempSpinWeightedScalar<3, 1>,
398 : ::Tags::TempSpinWeightedScalar<4, 1>,
399 : ::Tags::TempSpinWeightedScalar<5, 1>,
400 : // gauge Jacobians
401 : ::Tags::TempSpinWeightedScalar<6, 2>,
402 : ::Tags::TempSpinWeightedScalar<7, 0>,
403 : // the target Jacobian, then the step to it
404 : ::Tags::TempSpinWeightedScalar<8, 2>,
405 : // the potential eta
406 : ::Tags::TempSpinWeightedScalar<9, 1>>>
407 : computation_buffers{number_of_angular_points};
408 : auto& [eth_x, eth_y, eth_z, interpolated_eth_x, interpolated_eth_y,
409 : interpolated_eth_z, gauge_c, gauge_d, gauge_c_target, eta] =
410 : computation_buffers;
411 :
412 : unit_sphere_eth_cartesian_coordinates(make_not_null(&get(eth_x)),
413 : make_not_null(&get(eth_y)),
414 : make_not_null(&get(eth_z)), l_max);
415 :
416 : double max_error = 1.0;
417 : double previous_max_error = std::numeric_limits<double>::max();
418 : size_t number_of_steps = 0;
419 : Spectral::Swsh::SwshInterpolator iteration_interpolator;
420 :
421 : const auto evaluate_current_map = [&]() {
422 : update_angular_coordinates_and_jacobians(
423 : cartesian_cauchy_coordinates, angular_cauchy_coordinates,
424 : make_not_null(&iteration_interpolator), make_not_null(&gauge_c),
425 : make_not_null(&gauge_d), l_max);
426 : return target_function(make_not_null(&gauge_c_target), gauge_c, gauge_d,
427 : iteration_interpolator);
428 : };
429 :
430 : // As in `iteratively_adapt_angular_coordinates`, the best map is copied
431 : // before it is normalized, so that the rewind reproduces it bit for bit.
432 : auto unevaluated_cartesian_cauchy_coordinates = *cartesian_cauchy_coordinates;
433 : auto best_cartesian_cauchy_coordinates = *cartesian_cauchy_coordinates;
434 : double best_max_error = std::numeric_limits<double>::max();
435 :
436 : while (true) {
437 : unevaluated_cartesian_cauchy_coordinates = *cartesian_cauchy_coordinates;
438 : max_error = evaluate_current_map();
439 : if (max_error < best_max_error) {
440 : best_max_error = max_error;
441 : best_cartesian_cauchy_coordinates =
442 : unevaluated_cartesian_cauchy_coordinates;
443 : }
444 :
445 : check_angular_solve_error_threshold(max_error, error_threshold);
446 : ++number_of_steps;
447 : if (max_error < tolerance or number_of_steps > max_steps or
448 : (number_of_steps > 2 and
449 : max_error > plateau_factor * previous_max_error)) {
450 : break;
451 : }
452 : previous_max_error = max_error;
453 :
454 : // the step to the target Jacobian, then the potential that generates it
455 : get(gauge_c_target).data() -= get(gauge_c).data();
456 : Spectral::Swsh::angular_derivatives<
457 : tmpl::list<Spectral::Swsh::Tags::InverseEth>>(
458 : l_max, 1, make_not_null(&get(eta)), get(gauge_c_target));
459 :
460 : iteration_interpolator.interpolate(make_not_null(&get(interpolated_eth_x)),
461 : get(eth_x));
462 : iteration_interpolator.interpolate(make_not_null(&get(interpolated_eth_y)),
463 : get(eth_y));
464 : iteration_interpolator.interpolate(make_not_null(&get(interpolated_eth_z)),
465 : get(eth_z));
466 :
467 : // dx^i = Re(eta conj(eth x^i)) = v^A \partial_A x^i is tangent to the
468 : // sphere, since x^i eth x^i = 0; the radial part of its eth is removed when
469 : // the next iteration normalizes the coordinates.
470 : get<0>(*cartesian_cauchy_coordinates) +=
471 : real(get(eta).data() * conj(get(interpolated_eth_x).data()));
472 : get<1>(*cartesian_cauchy_coordinates) +=
473 : real(get(eta).data() * conj(get(interpolated_eth_y).data()));
474 : get<2>(*cartesian_cauchy_coordinates) +=
475 : real(get(eta).data() * conj(get(interpolated_eth_z).data()));
476 : }
477 :
478 : if (best_max_error < max_error) {
479 : *cartesian_cauchy_coordinates = best_cartesian_cauchy_coordinates;
480 : max_error = evaluate_current_map();
481 : }
482 :
483 : finalize_function(gauge_c, gauge_d, *angular_cauchy_coordinates,
484 : iteration_interpolator);
485 :
486 : // The loop evaluates once more than it updates.
487 : const size_t number_of_coordinate_updates = number_of_steps - 1;
488 : if (std::isnan(max_error) or max_error >= tolerance) {
489 : report_unconverged_angular_solve(if_unconverged, tolerance,
490 : number_of_coordinate_updates, max_error);
491 : }
492 : return {.max_error = max_error,
493 : .number_of_coordinate_updates = number_of_coordinate_updates};
494 : }
495 :
496 : double adjust_angular_coordinates_for_j(
497 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> volume_j,
498 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
499 : gsl::not_null<
500 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
501 : angular_cauchy_coordinates,
502 : const SpinWeighted<ComplexDataVector, 2>& surface_j, size_t l_max,
503 : double tolerance, size_t max_steps, bool adjust_volume_gauge);
504 : } // namespace detail
505 :
506 : /*!
507 : * \brief Apply a radius-independent angular gauge transformation to a volume
508 : * \f$J\f$, for use with initial data generation.
509 : *
510 : * \details Performs the gauge transformation to \f$\hat J\f$,
511 : *
512 : * \f{align*}{
513 : * \hat J = \frac{1}{4 \hat{\omega}^2} \left( \bar{\hat d}^2 J(\hat x^{\hat A})
514 : * + \hat c^2 \bar J(\hat x^{\hat A})
515 : * + 2 \hat c \bar{\hat d} K(\hat x^{\hat A}) \right).
516 : * \f}
517 : *
518 : * Where \f$\hat c\f$ and \f$\hat d\f$ are the spin-weighted angular Jacobian
519 : * factors computed by `GaugeUpdateJacobianFromCoords`, and \f$\hat \omega\f$ is
520 : * the conformal factor associated with the angular coordinate transformation.
521 : * Note that the right-hand sides with explicit \f$\hat x^{\hat A}\f$ dependence
522 : * must be interpolated and that \f$K = \sqrt{1 + J \bar J}\f$.
523 : */
524 1 : struct GaugeAdjustInitialJ {
525 0 : using boundary_tags =
526 : tmpl::list<Tags::PartiallyFlatGaugeC, Tags::PartiallyFlatGaugeD,
527 : Tags::PartiallyFlatGaugeOmega, Tags::CauchyAngularCoords,
528 : Spectral::Swsh::Tags::LMax>;
529 0 : using return_tags = tmpl::list<Tags::BondiJ>;
530 0 : using argument_tags = tmpl::append<boundary_tags>;
531 :
532 0 : static void apply(
533 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> volume_j,
534 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& gauge_c,
535 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& gauge_d,
536 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& gauge_omega,
537 : const tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>&
538 : cauchy_angular_coordinates,
539 : const Spectral::Swsh::SwshInterpolator& interpolator, size_t l_max);
540 : };
541 :
542 : /// \cond
543 : struct NoIncomingRadiation;
544 : struct ZeroNonSmooth;
545 : template <bool evolve_ccm>
546 : struct InverseCubic;
547 : template <bool evolve_ccm>
548 : struct InitializeJ;
549 : struct ConformalFactor;
550 : struct CauchySecondOrder;
551 : /// \endcond
552 :
553 : /*!
554 : * \brief Abstract base class for an initial hypersurface data generator for
555 : * Cce, when the partially flat Bondi-like coordinates are evolved.
556 : *
557 : * \details The algorithm is same as `InitializeJ<false>`, but with an
558 : * additional initialization for the partially flat Bondi-like coordinates. The
559 : * functions that are required to be overriden in the derived classes are:
560 : * - `InitializeJ::get_clone()`: should return a
561 : * `std::unique_ptr<InitializeJ<true>>` with cloned state.
562 : * - `InitializeJ::operator() const`: should take as arguments, first a
563 : * set of `gsl::not_null` pointers represented by `mutate_tags`, followed by a
564 : * set of `const` references to quantities represented by `argument_tags`. \note
565 : * The `InitializeJ::operator()` should be const, and therefore not alter
566 : * the internal state of the generator. This is compatible with all known
567 : * use-cases and permits the `InitializeJ` generator to be placed in the
568 : * `GlobalCache`.
569 : */
570 : template <>
571 1 : struct InitializeJ<true> : public PUP::able {
572 0 : using boundary_tags = tmpl::list<Tags::BoundaryValue<Tags::BondiJ>,
573 : Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>,
574 : Tags::BoundaryValue<Tags::BondiR>,
575 : Tags::BoundaryValue<Tags::BondiBeta>>;
576 :
577 0 : using mutate_tags =
578 : tmpl::list<Tags::BondiJ, Tags::CauchyCartesianCoords,
579 : Tags::CauchyAngularCoords, Tags::PartiallyFlatCartesianCoords,
580 : Tags::PartiallyFlatAngularCoords>;
581 0 : using return_tags = mutate_tags;
582 0 : using argument_tags =
583 : tmpl::push_back<boundary_tags, Tags::LMax, Tags::NumberOfRadialPoints>;
584 :
585 : // The evolution of inertial coordinates are allowed only when InverseCubic is
586 : // used
587 0 : using creatable_classes = tmpl::list<InverseCubic<true>>;
588 :
589 0 : InitializeJ() = default;
590 0 : explicit InitializeJ(CkMigrateMessage* /*msg*/) {}
591 :
592 0 : WRAPPED_PUPable_abstract(InitializeJ); // NOLINT
593 :
594 0 : virtual std::unique_ptr<InitializeJ<true>> get_clone() const = 0;
595 :
596 : /// \brief The interpolator this generator wants used to build the worldtube
597 : /// boundary value of \f$\partial_u \partial_r J\f$, or `nullptr` to use the
598 : /// evolution's `H5Interpolator`.
599 : ///
600 : /// \details That boundary value is consumed by the initial data alone, so a
601 : /// generator is free to choose its own time-interpolation order for it
602 : /// without affecting the evolution.
603 1 : virtual std::unique_ptr<intrp::SpanInterpolator> du_dr_j_interpolator()
604 : const {
605 : return nullptr;
606 : }
607 :
608 : // Each derived class declares its own `return_tags` and `argument_tags` and
609 : // implements a non-virtual `operator()`; the dispatch below picks the
610 : // correct dynamic type and forwards through `db::mutate_apply`.
611 : template <typename DbTags>
612 0 : void operator()(const gsl::not_null<db::DataBox<DbTags>*> box,
613 : const gsl::not_null<Parallel::NodeLock*> hdf5_lock) const {
614 : call_with_dynamic_type<void, creatable_classes>(
615 : this, [&](auto* const derived) {
616 : db::mutate_apply(*derived, box, hdf5_lock);
617 : });
618 : }
619 : };
620 :
621 : /*!
622 : * \brief Abstract base class for an initial hypersurface data generator for
623 : * Cce, when the partially flat Bondi-like coordinates are not evolved.
624 : *
625 : * \details The functions that are required to be overriden in the derived
626 : * classes are:
627 : * - `InitializeJ::get_clone()`: should return a
628 : * `std::unique_ptr<InitializeJ<false>>` with cloned state.
629 : * - `InitializeJ::operator() const`: should take as arguments, first a
630 : * set of `gsl::not_null` pointers represented by `mutate_tags`, followed by a
631 : * set of `const` references to quantities represented by `argument_tags`. \note
632 : * The `InitializeJ::operator()` should be const, and therefore not alter
633 : * the internal state of the generator. This is compatible with all known
634 : * use-cases and permits the `InitializeJ` generator to be placed in the
635 : * `GlobalCache`.
636 : */
637 : template <>
638 1 : struct InitializeJ<false> : public PUP::able {
639 : // Default boundary/mutate/argument tags used by the simple derived classes
640 : // (ConformalFactor, InverseCubic<false>, NoIncomingRadiation, ZeroNonSmooth).
641 : // A derived class can shadow these with its own list (see CauchySecondOrder),
642 : // in which case the per-class list is used during dispatch.
643 0 : using boundary_tags = tmpl::list<Tags::BoundaryValue<Tags::BondiJ>,
644 : Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>,
645 : Tags::BoundaryValue<Tags::BondiR>,
646 : Tags::BoundaryValue<Tags::BondiBeta>>;
647 :
648 0 : using mutate_tags = tmpl::list<Tags::BondiJ, Tags::CauchyCartesianCoords,
649 : Tags::CauchyAngularCoords>;
650 0 : using return_tags = mutate_tags;
651 0 : using argument_tags =
652 : tmpl::push_back<boundary_tags, Tags::LMax, Tags::NumberOfRadialPoints>;
653 :
654 0 : using creatable_classes =
655 : tmpl::list<ConformalFactor, InverseCubic<false>, NoIncomingRadiation,
656 : ZeroNonSmooth, CauchySecondOrder,
657 : ::Cce::Solutions::LinearizedBondiSachs_detail::InitializeJ::
658 : LinearizedBondiSachs>;
659 :
660 0 : InitializeJ() = default;
661 0 : explicit InitializeJ(CkMigrateMessage* /*msg*/) {}
662 :
663 0 : WRAPPED_PUPable_abstract(InitializeJ); // NOLINT
664 :
665 0 : virtual std::unique_ptr<InitializeJ<false>> get_clone() const = 0;
666 :
667 : /// \brief The interpolator this generator wants used to build the worldtube
668 : /// boundary value of \f$\partial_u \partial_r J\f$, or `nullptr` to use the
669 : /// evolution's `H5Interpolator`.
670 : ///
671 : /// \details That boundary value is consumed by the initial data alone, so a
672 : /// generator is free to choose its own time-interpolation order for it
673 : /// without affecting the evolution.
674 1 : virtual std::unique_ptr<intrp::SpanInterpolator> du_dr_j_interpolator()
675 : const {
676 : return nullptr;
677 : }
678 :
679 : // Each derived class declares its own `return_tags` and `argument_tags` and
680 : // implements a non-virtual `operator()` whose signature matches those tags.
681 : // The dispatch below picks the correct dynamic type and forwards through
682 : // `db::mutate_apply` so the per-class tag lists are honored.
683 : template <typename DbTags>
684 0 : void operator()(const gsl::not_null<db::DataBox<DbTags>*> box,
685 : const gsl::not_null<Parallel::NodeLock*> hdf5_lock) const {
686 : call_with_dynamic_type<void, creatable_classes>(
687 : this, [&](auto* const derived) {
688 : db::mutate_apply(*derived, box, hdf5_lock);
689 : });
690 : }
691 : };
692 : } // namespace InitializeJ
693 : } // namespace Cce
694 :
695 : // NOLINTBEGIN(misc-header-include-cycle)
696 : #include "Evolution/Systems/Cce/AnalyticSolutions/LinearizedBondiSachsInitializeJ.hpp"
697 : #include "Evolution/Systems/Cce/Initialize/CauchySecondOrder.hpp"
698 : #include "Evolution/Systems/Cce/Initialize/ConformalFactor.hpp"
699 : #include "Evolution/Systems/Cce/Initialize/InverseCubic.hpp"
700 : #include "Evolution/Systems/Cce/Initialize/NoIncomingRadiation.hpp"
701 : #include "Evolution/Systems/Cce/Initialize/ZeroNonSmooth.hpp"
702 : // NOLINTEND(misc-header-include-cycle)
|
