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