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