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 <memory>
8 : #include <string>
9 :
10 : #include "DataStructures/SpinWeighted.hpp"
11 : #include "DataStructures/Tags.hpp"
12 : #include "DataStructures/Tensor/TypeAliases.hpp"
13 : #include "Evolution/Systems/Cce/GaugeTransformBoundaryData.hpp"
14 : #include "Evolution/Systems/Cce/Tags.hpp"
15 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCollocation.hpp"
16 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshDerivatives.hpp"
17 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshInterpolation.hpp"
18 : #include "Options/String.hpp"
19 : #include "Parallel/NodeLock.hpp"
20 : #include "Parallel/Printf/Printf.hpp"
21 : #include "Utilities/Gsl.hpp"
22 : #include "Utilities/Serialization/CharmPupable.hpp"
23 : #include "Utilities/TMPL.hpp"
24 :
25 : /// \cond
26 : class ComplexDataVector;
27 : /// \endcond
28 :
29 : namespace Cce {
30 : /// Contains utilities and \ref DataBoxGroup mutators for generating data for
31 : /// \f$J\f$ on the initial CCE hypersurface.
32 : namespace InitializeJ {
33 :
34 : namespace detail {
35 : // used to provide a default for the finalize functor in
36 : // `iteratively_adapt_angular_coordinates`
37 : struct NoOpFinalize {
38 : void operator()(
39 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& /*gauge_c*/,
40 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& /*gauge_d*/,
41 : const tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>&
42 : /*angular_cauchy_coordinates*/,
43 : const Spectral::Swsh::SwshInterpolator& /*interpolator*/) const {
44 : }
45 : };
46 :
47 : // perform an iterative solve for the set of angular coordinates. The iteration
48 : // callable `iteration_function` must have function signature:
49 : //
50 : // double iteration_function(
51 : // const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
52 : // gauge_c_step,
53 : // const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
54 : // gauge_d_step,
55 : // const Scalar<SpinWeighted<ComplexDataVector, 2>>& gauge_c,
56 : // const Scalar<SpinWeighted<ComplexDataVector, 0>>& gauge_d,
57 : // const Spectral::Swsh::SwshInterpolator& iteration_interpolator);
58 : //
59 : // but need not be a function pointer -- a callable class or lambda will also
60 : // suffice.
61 : // For each step specified by the iteration function, the coordinates are
62 : // updated via \hat \eth \delta x^i = \delta c \eth x^i|_{x^i=\hat x^i}
63 : // + \delta \bar d (\eth x^i)|_{x^i=\hat x^i}
64 : // This coordinate update is exact, and comes from expanding the chain rule to
65 : // determine Jacobian factors. However, the result is not guaranteed to
66 : // produce the desired Jacobian c and d, because \delta c and \delta d are
67 : // not necessarily consistent with the underlying coordinates.
68 : // We then update the x^i by inverting \hat \eth, which is also exact, but
69 : // assumes a no l=0 contribution to the transformation.
70 : // Depending on the choice of approximations used to specify
71 : // `iteration_function`, though, the method can be slow to converge.
72 :
73 : // However, the iterations are typically fast, and the computation is for
74 : // initial data that needs to be computed only once during a simulation, so it
75 : // is not currently an optimization priority. If this function becomes a
76 : // bottleneck, the numerical procedure of the iterative method or the choice of
77 : // approximation used for `iteration_function` should be revisited.
78 : template <typename IterationFunctor, typename FinalizeFunctor = NoOpFinalize>
79 : double iteratively_adapt_angular_coordinates(
80 : const gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
81 : const gsl::not_null<
82 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
83 : angular_cauchy_coordinates,
84 : const size_t l_max, const double tolerance, const size_t max_steps,
85 : const double error_threshold, const IterationFunctor& iteration_function,
86 : const bool require_convergence,
87 : const FinalizeFunctor finalize_function = NoOpFinalize{}) {
88 : const size_t number_of_angular_points =
89 : Spectral::Swsh::number_of_swsh_collocation_points(l_max);
90 :
91 : Spectral::Swsh::create_angular_and_cartesian_coordinates(
92 : cartesian_cauchy_coordinates, angular_cauchy_coordinates, l_max);
93 :
94 : Variables<tmpl::list<
95 : // cartesian coordinates
96 : ::Tags::TempSpinWeightedScalar<0, 0>,
97 : ::Tags::TempSpinWeightedScalar<1, 0>,
98 : ::Tags::TempSpinWeightedScalar<2, 0>,
99 : // eth of cartesian coordinates
100 : ::Tags::TempSpinWeightedScalar<3, 1>,
101 : ::Tags::TempSpinWeightedScalar<4, 1>,
102 : ::Tags::TempSpinWeightedScalar<5, 1>,
103 : // eth of gauge-transformed cartesian coordinates
104 : ::Tags::TempSpinWeightedScalar<6, 1>,
105 : ::Tags::TempSpinWeightedScalar<7, 1>,
106 : ::Tags::TempSpinWeightedScalar<8, 1>,
107 : // gauge Jacobians
108 : ::Tags::TempSpinWeightedScalar<9, 2>,
109 : ::Tags::TempSpinWeightedScalar<10, 0>,
110 : // gauge Jacobians on next iteration
111 : ::Tags::TempSpinWeightedScalar<11, 2>,
112 : ::Tags::TempSpinWeightedScalar<12, 0>,
113 : // cartesian coordinates steps
114 : ::Tags::TempSpinWeightedScalar<13, 0>,
115 : ::Tags::TempSpinWeightedScalar<14, 0>,
116 : ::Tags::TempSpinWeightedScalar<15, 0>>>
117 : computation_buffers{number_of_angular_points};
118 :
119 : auto& x = get(get<::Tags::TempSpinWeightedScalar<0, 0>>(computation_buffers));
120 : auto& y = get(get<::Tags::TempSpinWeightedScalar<1, 0>>(computation_buffers));
121 : auto& z = get(get<::Tags::TempSpinWeightedScalar<2, 0>>(computation_buffers));
122 :
123 : x.data() =
124 : std::complex<double>(1.0, 0.0) * get<0>(*cartesian_cauchy_coordinates);
125 : y.data() =
126 : std::complex<double>(1.0, 0.0) * get<1>(*cartesian_cauchy_coordinates);
127 : z.data() =
128 : std::complex<double>(1.0, 0.0) * get<2>(*cartesian_cauchy_coordinates);
129 :
130 : auto& eth_x =
131 : get(get<::Tags::TempSpinWeightedScalar<3, 1>>(computation_buffers));
132 : auto& eth_y =
133 : get(get<::Tags::TempSpinWeightedScalar<4, 1>>(computation_buffers));
134 : auto& eth_z =
135 : get(get<::Tags::TempSpinWeightedScalar<5, 1>>(computation_buffers));
136 :
137 : Spectral::Swsh::angular_derivatives<
138 : tmpl::list<Spectral::Swsh::Tags::Eth, Spectral::Swsh::Tags::Eth,
139 : Spectral::Swsh::Tags::Eth>>(l_max, 1, make_not_null(ð_x),
140 : make_not_null(ð_y),
141 : make_not_null(ð_z), x, y, z);
142 :
143 : auto& evolution_gauge_eth_x_step =
144 : get(get<::Tags::TempSpinWeightedScalar<6, 1>>(computation_buffers));
145 : auto& evolution_gauge_eth_y_step =
146 : get(get<::Tags::TempSpinWeightedScalar<7, 1>>(computation_buffers));
147 : auto& evolution_gauge_eth_z_step =
148 : get(get<::Tags::TempSpinWeightedScalar<8, 1>>(computation_buffers));
149 :
150 : auto& gauge_c =
151 : get<::Tags::TempSpinWeightedScalar<9, 2>>(computation_buffers);
152 : auto& gauge_d =
153 : get<::Tags::TempSpinWeightedScalar<10, 0>>(computation_buffers);
154 :
155 : auto& gauge_c_step =
156 : get<::Tags::TempSpinWeightedScalar<11, 2>>(computation_buffers);
157 : auto& gauge_d_step =
158 : get<::Tags::TempSpinWeightedScalar<12, 0>>(computation_buffers);
159 :
160 : auto& x_step =
161 : get(get<::Tags::TempSpinWeightedScalar<13, 0>>(computation_buffers));
162 : auto& y_step =
163 : get(get<::Tags::TempSpinWeightedScalar<14, 0>>(computation_buffers));
164 : auto& z_step =
165 : get(get<::Tags::TempSpinWeightedScalar<15, 0>>(computation_buffers));
166 :
167 : double max_error = 1.0;
168 : size_t number_of_steps = 0;
169 : Spectral::Swsh::SwshInterpolator iteration_interpolator;
170 : while (true) {
171 : GaugeUpdateAngularFromCartesian<
172 : Tags::CauchyAngularCoords,
173 : Tags::CauchyCartesianCoords>::apply(angular_cauchy_coordinates,
174 : cartesian_cauchy_coordinates);
175 :
176 : iteration_interpolator = Spectral::Swsh::SwshInterpolator{
177 : get<0>(*angular_cauchy_coordinates),
178 : get<1>(*angular_cauchy_coordinates), l_max};
179 :
180 : GaugeUpdateJacobianFromCoordinates<
181 : Tags::PartiallyFlatGaugeC, Tags::PartiallyFlatGaugeD,
182 : Tags::CauchyAngularCoords,
183 : Tags::CauchyCartesianCoords>::apply(make_not_null(&gauge_c),
184 : make_not_null(&gauge_d),
185 : angular_cauchy_coordinates,
186 : *cartesian_cauchy_coordinates,
187 : l_max);
188 :
189 : max_error = iteration_function(make_not_null(&gauge_c_step),
190 : make_not_null(&gauge_d_step), gauge_c,
191 : gauge_d, iteration_interpolator);
192 :
193 : if (max_error > error_threshold) {
194 : ERROR(
195 : "Iterative solve for surface coordinates of initial data failed. The "
196 : "strain is too large to be fully eliminated by a well-behaved "
197 : "alteration of the spherical mesh. This could be an indication that "
198 : "there is an issue with the worldtube data. If you are confident "
199 : "the worldtube data is correct, then please use an alternative "
200 : "initial data generator such as `InverseCubic`. If that fails, "
201 : "please double check that your spherical harmonic modes are decaying "
202 : "correctly with increasing (l,m).\nError: "
203 : << max_error << "\nError threshold: " << error_threshold);
204 : }
205 : ++number_of_steps;
206 : if (max_error < tolerance or number_of_steps > max_steps) {
207 : break;
208 : }
209 : // using the evolution_gauge_.._step as temporary buffers for the
210 : // interpolation results
211 : iteration_interpolator.interpolate(
212 : make_not_null(&evolution_gauge_eth_x_step), eth_x);
213 : iteration_interpolator.interpolate(
214 : make_not_null(&evolution_gauge_eth_y_step), eth_y);
215 : iteration_interpolator.interpolate(
216 : make_not_null(&evolution_gauge_eth_z_step), eth_z);
217 :
218 : evolution_gauge_eth_x_step =
219 : 0.5 * ((get(gauge_c_step)) * conj(evolution_gauge_eth_x_step) +
220 : conj((get(gauge_d_step))) * evolution_gauge_eth_x_step);
221 : evolution_gauge_eth_y_step =
222 : 0.5 * ((get(gauge_c_step)) * conj(evolution_gauge_eth_y_step) +
223 : conj((get(gauge_d_step))) * evolution_gauge_eth_y_step);
224 : evolution_gauge_eth_z_step =
225 : 0.5 * ((get(gauge_c_step)) * conj(evolution_gauge_eth_z_step) +
226 : conj((get(gauge_d_step))) * evolution_gauge_eth_z_step);
227 :
228 : Spectral::Swsh::angular_derivatives<tmpl::list<
229 : Spectral::Swsh::Tags::InverseEth, Spectral::Swsh::Tags::InverseEth,
230 : Spectral::Swsh::Tags::InverseEth>>(
231 : l_max, 1, make_not_null(&x_step), make_not_null(&y_step),
232 : make_not_null(&z_step), evolution_gauge_eth_x_step,
233 : evolution_gauge_eth_y_step, evolution_gauge_eth_z_step);
234 :
235 : get<0>(*cartesian_cauchy_coordinates) += real(x_step.data());
236 : get<1>(*cartesian_cauchy_coordinates) += real(y_step.data());
237 : get<2>(*cartesian_cauchy_coordinates) += real(z_step.data());
238 : }
239 :
240 : finalize_function(gauge_c, gauge_d, *angular_cauchy_coordinates,
241 : iteration_interpolator);
242 :
243 : if (tolerance < max_error) {
244 : if (require_convergence) {
245 : ERROR(
246 : "Initial data iterative angular solve did not reach "
247 : "target tolerance "
248 : << tolerance << ".\n"
249 : << "Exited after " << max_steps
250 : << " iterations, achieving final\n"
251 : "maximum over collocation points deviation of J from target of "
252 : << max_error);
253 : } else {
254 : Parallel::printf(
255 : "Warning: iterative angular solve did not reach "
256 : "target tolerance %e.\n"
257 : "Exited after %zu iterations, achieving final maximum over "
258 : "collocation points for deviation from target of %e\n"
259 : "Proceeding with evolution using the partial result from partial "
260 : "angular solve.\n",
261 : tolerance, max_steps, max_error);
262 : }
263 : }
264 : return max_error;
265 : }
266 :
267 : double adjust_angular_coordinates_for_j(
268 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> volume_j,
269 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
270 : gsl::not_null<
271 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
272 : angular_cauchy_coordinates,
273 : const SpinWeighted<ComplexDataVector, 2>& surface_j, size_t l_max,
274 : double tolerance, size_t max_steps, bool adjust_volume_gauge);
275 : } // namespace detail
276 :
277 : /*!
278 : * \brief Apply a radius-independent angular gauge transformation to a volume
279 : * \f$J\f$, for use with initial data generation.
280 : *
281 : * \details Performs the gauge transformation to \f$\hat J\f$,
282 : *
283 : * \f{align*}{
284 : * \hat J = \frac{1}{4 \hat{\omega}^2} \left( \bar{\hat d}^2 J(\hat x^{\hat A})
285 : * + \hat c^2 \bar J(\hat x^{\hat A})
286 : * + 2 \hat c \bar{\hat d} K(\hat x^{\hat A}) \right).
287 : * \f}
288 : *
289 : * Where \f$\hat c\f$ and \f$\hat d\f$ are the spin-weighted angular Jacobian
290 : * factors computed by `GaugeUpdateJacobianFromCoords`, and \f$\hat \omega\f$ is
291 : * the conformal factor associated with the angular coordinate transformation.
292 : * Note that the right-hand sides with explicit \f$\hat x^{\hat A}\f$ dependence
293 : * must be interpolated and that \f$K = \sqrt{1 + J \bar J}\f$.
294 : */
295 1 : struct GaugeAdjustInitialJ {
296 0 : using boundary_tags =
297 : tmpl::list<Tags::PartiallyFlatGaugeC, Tags::PartiallyFlatGaugeD,
298 : Tags::PartiallyFlatGaugeOmega, Tags::CauchyAngularCoords,
299 : Spectral::Swsh::Tags::LMax>;
300 0 : using return_tags = tmpl::list<Tags::BondiJ>;
301 0 : using argument_tags = tmpl::append<boundary_tags>;
302 :
303 0 : static void apply(
304 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> volume_j,
305 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& gauge_c,
306 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& gauge_d,
307 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& gauge_omega,
308 : const tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>&
309 : cauchy_angular_coordinates,
310 : const Spectral::Swsh::SwshInterpolator& interpolator, size_t l_max);
311 : };
312 :
313 : /// \cond
314 : struct NoIncomingRadiation;
315 : struct ZeroNonSmooth;
316 : template <bool evolve_ccm>
317 : struct InverseCubic;
318 : template <bool evolve_ccm>
319 : struct InitializeJ;
320 : struct ConformalFactor;
321 : /// \endcond
322 :
323 : /*!
324 : * \brief Abstract base class for an initial hypersurface data generator for
325 : * Cce, when the partially flat Bondi-like coordinates are evolved.
326 : *
327 : * \details The algorithm is same as `InitializeJ<false>`, but with an
328 : * additional initialization for the partially flat Bondi-like coordinates. The
329 : * functions that are required to be overriden in the derived classes are:
330 : * - `InitializeJ::get_clone()`: should return a
331 : * `std::unique_ptr<InitializeJ<true>>` with cloned state.
332 : * - `InitializeJ::operator() const`: should take as arguments, first a
333 : * set of `gsl::not_null` pointers represented by `mutate_tags`, followed by a
334 : * set of `const` references to quantities represented by `argument_tags`. \note
335 : * The `InitializeJ::operator()` should be const, and therefore not alter
336 : * the internal state of the generator. This is compatible with all known
337 : * use-cases and permits the `InitializeJ` generator to be placed in the
338 : * `GlobalCache`.
339 : */
340 : template <>
341 1 : struct InitializeJ<true> : public PUP::able {
342 0 : using boundary_tags = tmpl::list<Tags::BoundaryValue<Tags::BondiJ>,
343 : Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>,
344 : Tags::BoundaryValue<Tags::BondiR>,
345 : Tags::BoundaryValue<Tags::BondiBeta>>;
346 :
347 0 : using mutate_tags =
348 : tmpl::list<Tags::BondiJ, Tags::CauchyCartesianCoords,
349 : Tags::CauchyAngularCoords, Tags::PartiallyFlatCartesianCoords,
350 : Tags::PartiallyFlatAngularCoords>;
351 0 : using argument_tags =
352 : tmpl::push_back<boundary_tags, Tags::LMax, Tags::NumberOfRadialPoints>;
353 :
354 : // The evolution of inertial coordinates are allowed only when InverseCubic is
355 : // used
356 0 : using creatable_classes = tmpl::list<InverseCubic<true>>;
357 :
358 0 : InitializeJ() = default;
359 0 : explicit InitializeJ(CkMigrateMessage* /*msg*/) {}
360 :
361 0 : WRAPPED_PUPable_abstract(InitializeJ); // NOLINT
362 :
363 0 : virtual std::unique_ptr<InitializeJ<true>> get_clone() const = 0;
364 :
365 0 : virtual void operator()(
366 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> j,
367 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
368 : gsl::not_null<
369 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
370 : angular_cauchy_coordinates,
371 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_inertial_coordinates,
372 : gsl::not_null<
373 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
374 : angular_inertial_coordinates,
375 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_j,
376 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_dr_j,
377 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& r,
378 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& beta, size_t l_max,
379 : size_t number_of_radial_points,
380 : gsl::not_null<Parallel::NodeLock*> hdf5_lock) const = 0;
381 : };
382 :
383 : /*!
384 : * \brief Abstract base class for an initial hypersurface data generator for
385 : * Cce, when the partially flat Bondi-like coordinates are not evolved.
386 : *
387 : * \details The functions that are required to be overriden in the derived
388 : * classes are:
389 : * - `InitializeJ::get_clone()`: should return a
390 : * `std::unique_ptr<InitializeJ<false>>` with cloned state.
391 : * - `InitializeJ::operator() const`: should take as arguments, first a
392 : * set of `gsl::not_null` pointers represented by `mutate_tags`, followed by a
393 : * set of `const` references to quantities represented by `argument_tags`. \note
394 : * The `InitializeJ::operator()` should be const, and therefore not alter
395 : * the internal state of the generator. This is compatible with all known
396 : * use-cases and permits the `InitializeJ` generator to be placed in the
397 : * `GlobalCache`.
398 : */
399 : template <>
400 1 : struct InitializeJ<false> : public PUP::able {
401 0 : using boundary_tags = tmpl::list<Tags::BoundaryValue<Tags::BondiJ>,
402 : Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>,
403 : Tags::BoundaryValue<Tags::BondiR>,
404 : Tags::BoundaryValue<Tags::BondiBeta>>;
405 :
406 0 : using mutate_tags = tmpl::list<Tags::BondiJ, Tags::CauchyCartesianCoords,
407 : Tags::CauchyAngularCoords>;
408 0 : using argument_tags =
409 : tmpl::push_back<boundary_tags, Tags::LMax, Tags::NumberOfRadialPoints>;
410 :
411 0 : using creatable_classes = tmpl::list<ConformalFactor, InverseCubic<false>,
412 : NoIncomingRadiation, ZeroNonSmooth>;
413 :
414 0 : InitializeJ() = default;
415 0 : explicit InitializeJ(CkMigrateMessage* /*msg*/) {}
416 :
417 0 : WRAPPED_PUPable_abstract(InitializeJ); // NOLINT
418 :
419 0 : virtual std::unique_ptr<InitializeJ<false>> get_clone() const = 0;
420 :
421 0 : virtual void operator()(
422 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> j,
423 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
424 : gsl::not_null<
425 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
426 : angular_cauchy_coordinates,
427 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_j,
428 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_dr_j,
429 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& r,
430 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& beta, size_t l_max,
431 : size_t number_of_radial_points,
432 : gsl::not_null<Parallel::NodeLock*> hdf5_lock) const = 0;
433 : };
434 : } // namespace InitializeJ
435 : } // namespace Cce
436 :
437 : #include "Evolution/Systems/Cce/Initialize/InverseCubic.hpp"
438 : #include "Evolution/Systems/Cce/Initialize/NoIncomingRadiation.hpp"
439 : #include "Evolution/Systems/Cce/Initialize/ZeroNonSmooth.hpp"
|
