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 <limits>
8 : #include <memory>
9 : #include <string>
10 :
11 : #include "DataStructures/SpinWeighted.hpp"
12 : #include "DataStructures/Tensor/TypeAliases.hpp"
13 : #include "Evolution/Systems/Cce/Initialize/InitializeJ.hpp"
14 : #include "NumericalAlgorithms/Interpolation/SpanInterpolator.hpp"
15 : #include "Utilities/Gsl.hpp"
16 : #include "Utilities/Serialization/CharmPupable.hpp"
17 : #include "Utilities/TMPL.hpp"
18 :
19 : /// \cond
20 : class ComplexDataVector;
21 : /// \endcond
22 :
23 1 : namespace Cce::InitializeJ {
24 :
25 : namespace CauchySecondOrder_detail {
26 :
27 : /// \brief The dependence of the coefficients of the Cauchy-gauge ansatz on
28 : /// the free parameter \f$x = \tilde J^{(2)}\f$.
29 : ///
30 : /// \details The matching conditions at the worldtube are linear, so each
31 : /// coefficient \f$\tilde J^{(n)}\f$ is an affine function of \f$x\f$ with slope
32 : /// \f$\alpha^{(n)} = (-1)^n 2^{3-n} / (n! (3-n)!)\f$.
33 : /// @{
34 : constexpr double j0_x_coefficient = 4.0 / 3.0;
35 : constexpr double j1_x_coefficient = -2.0;
36 : constexpr double j3_x_coefficient = -1.0 / 6.0;
37 : /// @}
38 :
39 : /*!
40 : * \brief Compute the coefficients \f$\tilde J^{(0)}\f$, \f$\tilde J^{(1)}\f$
41 : * and \f$\tilde J^{(3)}\f$ of the Cauchy-gauge ansatz for a given
42 : * \f$x = \tilde J^{(2)}\f$.
43 : *
44 : * \details The ansatz is
45 : *
46 : * \f{align*}{
47 : * \tilde J(\tilde y; x) = \tilde J^{(0)}(x)
48 : * + \tilde J^{(1)}(x)\,(1 - \tilde y) + x\,(1 - \tilde y)^2
49 : * + \tilde J^{(3)}(x)\,(1 - \tilde y)^3 ,
50 : * \f}
51 : *
52 : * and its coefficients are fixed by matching the worldtube
53 : * (\f$\tilde y = -1\f$) values \f$\tilde J|_\Gamma\f$,
54 : * \f$\partial_{\tilde y} \tilde J|_\Gamma = (R/2)\, \partial_r J|_\Gamma\f$ and
55 : * \f$\partial_{\tilde y}^2 \tilde J|_\Gamma\f$, where \f$R\f$ is the worldtube
56 : * radius. They are affine in \f$x\f$, with slopes `j0_x_coefficient`,
57 : * `j1_x_coefficient` and `j3_x_coefficient`.
58 : */
59 : void radial_ansatz_coefficients(gsl::not_null<ComplexDataVector*> j0,
60 : gsl::not_null<ComplexDataVector*> j1,
61 : gsl::not_null<ComplexDataVector*> j3,
62 : const ComplexDataVector& j2,
63 : const ComplexDataVector& boundary_j,
64 : const ComplexDataVector& boundary_dr_j,
65 : const ComplexDataVector& boundary_dy_dy_j,
66 : const ComplexDataVector& boundary_r);
67 :
68 : /*!
69 : * \brief Determine \f$x = \tilde J^{(2)}\f$ of the Cauchy-gauge ansatz by
70 : * fixed-point iteration.
71 : *
72 : * \details Matching the worldtube value of \f$\tilde J\f$ and its first two
73 : * radial derivatives leaves \f$x\f$ to be determined by the remaining
74 : * partially flat condition \f$\breve J^{(2)} = 0\f$. Written in Cauchy-gauge
75 : * quantities, it is the nonlinear equation
76 : *
77 : * \f{align*}{
78 : * x = \Phi(x) \equiv \frac{1}{2} \tilde J^{(0)}
79 : * \left[\big|\tilde J^{(1)}\big|^2 - \big(\tilde K^{(1)}\big)^2\right],
80 : * \quad
81 : * \tilde K^{(1)} = \frac{\Re\big(\tilde J^{(0)} \bar{\tilde J}^{(1)}\big)}
82 : * {\sqrt{1 + |\tilde J^{(0)}|^2}},
83 : * \f}
84 : *
85 : * where \f$\tilde J^{(0)}\f$ and \f$\tilde J^{(1)}\f$ are the affine functions
86 : * of \f$x\f$ fixed by the matching conditions, whose \f$x\f$-independent parts
87 : * `j0_at_zero` and `j1_at_zero` are supplied by the caller. With
88 : *
89 : * \f{align*}{
90 : * \varrho = \max\left(\big\lVert \tilde J^{(0)}(0) \big\rVert_\infty,
91 : * \big\lVert \tilde J^{(1)}(0) \big\rVert_\infty\right),
92 : * \f}
93 : *
94 : * \f$\Phi\f$ maps the disc \f$|x| \le \varrho\f$ into itself and is a
95 : * contraction there whenever \f$\varrho \le 7\times 10^{-2}\f$, so the
96 : * iteration from \f$x = 0\f$ converges to the unique root in the disc.
97 : * Worldtube data in the radiation zone satisfy this bound comfortably, and one
98 : * or two iterations suffice for them.
99 : *
100 : * Iteration stops once the largest change of \f$x\f$ over the angular
101 : * collocation points falls below `tolerance`, or after `max_iterations`
102 : * iterations. `max_iterations` of zero performs no iterations and leaves
103 : * \f$x = 0\f$.
104 : *
105 : * \returns the number of iterations performed; `final_step` is set to the
106 : * largest change of \f$x\f$ on the last of them (infinite if none were taken).
107 : * If \f$x\f$ becomes non-finite, because the inputs contain non-finite values
108 : * or the iteration diverges, the iteration stops and `final_step` is set to
109 : * NaN.
110 : */
111 : size_t solve_asymptotic_j2(gsl::not_null<ComplexDataVector*> j2,
112 : gsl::not_null<double*> final_step,
113 : const ComplexDataVector& j0_at_zero,
114 : const ComplexDataVector& j1_at_zero,
115 : double tolerance, size_t max_iterations);
116 :
117 : /*!
118 : * \brief The sentence of the `MaxPartiallyFlatJ0` error that says whether
119 : * raising `J0MaxIterations` can help.
120 : *
121 : * \details Empty if no linearized sweeps ran. Otherwise it says whether the
122 : * sweeps stopped because \f$\max|\breve J^{(0)}|\f$ had stopped improving
123 : * (`stopped_improving`), or because the angular solve reached
124 : * `j0_max_iterations` iterations, either while \f$\max|\breve J^{(0)}|\f$
125 : * was still
126 : * improving or before the `CauchySecondOrder::j0_plateau_sweeps` sweeps that
127 : * telling the two apart takes. If `j0_max_iterations` is already at the upper
128 : * bound of the `J0MaxIterations` option, it says so instead of suggesting a
129 : * larger value.
130 : */
131 : std::string j0_max_iterations_hint(size_t number_of_linearized_sweeps,
132 : bool stopped_improving,
133 : size_t j0_max_iterations);
134 : } // namespace CauchySecondOrder_detail
135 :
136 : /*!
137 : * \brief Initialize \f$J\f$ on the first hypersurface by second-order Cauchy
138 : * matching at the worldtube.
139 : *
140 : * \details Here \f$\tilde J\f$ denotes \f$J\f$ in the Cauchy coordinates and
141 : * \f$\breve J\f$ its counterpart in the partially flat coordinates, and
142 : * \f$\tilde J^{(n)}\f$ and \f$\breve J^{(n)}\f$ are their coefficients of
143 : * \f$(1 - \tilde y)^n\f$, where \f$\tilde y = 1 - 2R/r\f$ is the compactified
144 : * radial coordinate and \f$R\f$ the worldtube radius. The initial data must
145 : * satisfy the partially flat conditions \f$\breve J^{(0)} = 0\f$ and
146 : * \f$\breve J^{(2)} = 0\f$. They are constructed in three steps:
147 : *
148 : * 1. Match the initial \f$\tilde J\f$ to the worldtube value and its first two
149 : * radial derivatives in the Cauchy gauge, using the ansatz of
150 : * `CauchySecondOrder_detail::radial_ansatz_coefficients`, and solve the
151 : * fixed-point problem that the remaining condition \f$\breve J^{(2)} = 0\f$
152 : * poses for \f$x = \tilde J^{(2)}\f$
153 : * (`CauchySecondOrder_detail::solve_asymptotic_j2`). The worldtube
154 : * \f$\partial_{\tilde y}^2 \tilde J|_\Gamma\f$ is obtained by inverting the
155 : * \f$H\f$-hypersurface equation at the worldtube
156 : * (`CauchySecondOrder_detail::compute_dy_dy_j`).
157 : * 2. Use the resulting \f$\tilde J^{(0)}\f$ to solve the Beltrami problem
158 : * \f$\breve J^{(0)} = 0\f$ for the initial angular map, as described below.
159 : * 3. Transform the whole profile into the partially flat gauge with this map
160 : * (`GaugeAdjustInitialJ`).
161 : *
162 : * In step 2, \f$\breve J^{(0)} = 0\f$ specifies the target angular Jacobian
163 : *
164 : * \f{align*}{
165 : * \hat a_\star = -\frac{\bar{\hat b}_\star\, \tilde J^{(0)}}
166 : * {1 + \tilde K^{(0)}}, \qquad
167 : * \tilde K^{(0)} = \sqrt{1 + \big|\tilde J^{(0)}\big|^2}
168 : * \f}
169 : *
170 : * at every point of the sphere, where \f$\hat a\f$ and \f$\hat b\f$ are the
171 : * spin-weighted Jacobian factors of the angular map (\f$\hat c\f$ and
172 : * \f$\hat d\f$ in `GaugeUpdateJacobianFromCoordinates`). The map is recovered
173 : * iteratively, in up to two stages that together take at most
174 : * `J0MaxIterations` iterations. A potential-based solve
175 : * (`detail::adapt_angular_coordinates_via_potential`) generates each
176 : * displacement of the Cartesian coordinates \f$\hat x^i\f$ of the map from a
177 : * single spin-weight-1 potential \f$\eta\f$,
178 : *
179 : * \f{align*}{
180 : * \delta \hat x^i = \mathrm{Re}\big(\eta\, \overline{\eth \hat x^i}\big),
181 : * \qquad \eta = \eth^{-1}\big(\hat a_\star - \hat a\big),
182 : * \f}
183 : *
184 : * and runs to its minimum residual, which takes a few iterations. If that
185 : * brings \f$\max|\breve J^{(0)}|\f$ within `MaxPartiallyFlatJ0`, its map is the
186 : * solution. Otherwise linearized sweeps
187 : * (`detail::iteratively_adapt_angular_coordinates`) continue from it until the
188 : * smallest \f$\max|\breve J^{(0)}|\f$ reached has improved by less than 1% over
189 : * the last 50 sweeps, or until the two stages together have taken
190 : * `J0MaxIterations` iterations. Both stages keep the best map they evaluated.
191 : *
192 : * Initialization aborts if \f$\max|\tilde J^{(0)}(0)|\f$, or
193 : * \f$\max|\breve J^{(0)}|\f$ at any iteration of step 2, exceeds
194 : * `MaxCauchyJ0`, or if the fixed-point iteration of step 1 does not reach
195 : * `J2Tolerance` within `J2MaxIterations` iterations. After step 3 it prints a
196 : * summary of steps 1 and 2 and of \f$\max|\breve J^{(0)}|\f$ and
197 : * \f$\max|\breve J^{(2)}|\f$, and aborts if either exceeds its bound,
198 : * `MaxPartiallyFlatJ0` or `MaxPartiallyFlatJ2`.
199 : *
200 : * The worldtube \f$\partial_{\tilde u}\partial_r J|_\Gamma\f$ that enters the
201 : * \f$H\f$-hypersurface equation is obtained by differentiating an interpolant
202 : * of the worldtube \f$\partial_r J|_\Gamma\f$ in time with
203 : * `DuDrJInterpolator`, which this generator supplies to the worldtube data
204 : * manager. It is separate from the evolution's `H5Interpolator` because the
205 : * match wants a low interpolation order for that derivative while the
206 : * evolution wants a high one for the values it interpolates every step.
207 : */
208 1 : struct CauchySecondOrder : InitializeJ<false> {
209 : /// The largest number of iterations of the potential solve. It reaches its
210 : /// minimum within a few, so this cap only bounds its cost on data it cannot
211 : /// solve.
212 1 : static constexpr size_t max_potential_iterations = 20;
213 : /// The linearized sweeps continue while the smallest
214 : /// \f$\max|\breve J^{(0)}|\f$ reached improves by at least the fraction
215 : /// `j0_plateau_improvement` over `j0_plateau_sweeps` sweeps.
216 : /// @{
217 1 : static constexpr size_t j0_plateau_sweeps = 50;
218 1 : static constexpr double j0_plateau_improvement = 1.0e-2;
219 : /// @}
220 :
221 0 : struct MaxPartiallyFlatJ0 {
222 0 : using type = double;
223 0 : static constexpr Options::String help = {
224 : "Largest allowed max|J^(0)| = max|J| at scri+ in the partially flat "
225 : "gauge. The potential solve that runs first is enough if it reaches "
226 : "this value; otherwise linearized sweeps follow until max|J^(0)| "
227 : "stops improving, and initialization aborts if it is still above this "
228 : "value."};
229 0 : static type lower_bound() { return 1.0e-16; }
230 0 : static type upper_bound() { return 1.0; }
231 0 : static type suggested_value() { return 5.0e-12; }
232 : };
233 :
234 0 : struct J0MaxIterations {
235 0 : using type = size_t;
236 0 : static constexpr Options::String help = {
237 : "Largest number of angular coordinate iterations, shared by the "
238 : "potential solve and the linearized sweeps that follow it if it does "
239 : "not reach MaxPartiallyFlatJ0."};
240 0 : static type lower_bound() { return 10; }
241 0 : static type upper_bound() { return 2000; }
242 0 : static type suggested_value() { return 1500; }
243 : };
244 :
245 0 : struct J2Tolerance {
246 0 : using type = double;
247 0 : static constexpr Options::String help = {
248 : "Convergence tolerance of the fixed-point iteration that determines "
249 : "the second radial expansion coefficient J^(2) of the Cauchy-gauge "
250 : "ansatz from the partially flat gauge condition. The iteration stops "
251 : "once the largest change of J^(2) over the collocation points falls "
252 : "below this value; the comparison is absolute, as J is dimensionless "
253 : "and of the size of the strain."};
254 0 : static type lower_bound() { return 1.0e-16; }
255 0 : static type upper_bound() { return 1.0e-3; }
256 0 : static type suggested_value() { return 1.0e-14; }
257 : };
258 :
259 0 : struct J2MaxIterations {
260 0 : using type = size_t;
261 0 : static constexpr Options::String help = {
262 : "Largest number of fixed-point iterations used to determine the "
263 : "Cauchy-gauge J^(2). For worldtube data in the radiation zone one or "
264 : "two iterations suffice. "
265 : "If this is nonzero and the iteration has not converged after this "
266 : "many iterations, initialization aborts. "
267 : "Zero skips the solve and leaves J^(2) = 0, which does not satisfy "
268 : "the partially flat gauge condition."};
269 0 : static type upper_bound() { return 100; }
270 0 : static type suggested_value() { return 10; }
271 : };
272 :
273 0 : struct MaxCauchyJ0 {
274 0 : using type = double;
275 0 : static constexpr Options::String help = {
276 : "Largest allowed max|J^(0)| in the Cauchy gauge, where J^(0) is J at "
277 : "scri+. The same bound also limits the partially flat J^(0) during "
278 : "the angular iterations as a divergence guard. Exceeding either "
279 : "bound aborts initialization. MaxPartiallyFlatJ0 bounds the partially "
280 : "flat J^(0) after the angular solve. Larger bounds allow larger "
281 : "asymptotic strains but risk a poorly behaved angular coordinate map."};
282 0 : static type lower_bound() { return 1.0e-14; }
283 0 : static type upper_bound() { return 1.0e2; }
284 0 : static type suggested_value() { return 5.0e-2; }
285 : };
286 :
287 0 : struct MaxPartiallyFlatJ2 {
288 0 : using type = double;
289 0 : static constexpr Options::String help = {
290 : "Largest allowed max|J^(2)| in the partially flat gauge after the "
291 : "gauge transformation, where J^(2) = (1/2) Dy^2 J at scri+ is the "
292 : "coefficient of (1 - y)^2. The fixed-point solve for the Cauchy-gauge "
293 : "J^(2) makes it vanish up to discretization error, so initialization "
294 : "aborts above "
295 : "this value, which usually means the angular or radial resolution is "
296 : "too low for the worldtube data."};
297 0 : static type lower_bound() { return 1.0e-16; }
298 0 : static type upper_bound() { return 1.0; }
299 0 : static type suggested_value() { return 1.0e-12; }
300 : };
301 :
302 0 : struct DuDrJInterpolator {
303 0 : using type = std::unique_ptr<intrp::SpanInterpolator>;
304 0 : static constexpr Options::String help = {
305 : "Interpolator used to time-differentiate the worldtube Dr(J) into the "
306 : "Du(Dr(J)) that this generator matches to. That boundary value feeds "
307 : "the initial data alone, so this is independent of `H5Interpolator`: "
308 : "the evolution wants a high interpolation order, while the "
309 : "second-order match wants a low one (a barycentric order of 2 to 4), "
310 : "which keeps the high-frequency content of the worldtube out of the "
311 : "initial data."};
312 : };
313 :
314 0 : using options =
315 : tmpl::list<MaxPartiallyFlatJ0, J0MaxIterations, MaxCauchyJ0, J2Tolerance,
316 : J2MaxIterations, MaxPartiallyFlatJ2, DuDrJInterpolator>;
317 0 : static constexpr Options::String help = {
318 : "Second-order initial data generator for the Cauchy CCE evolution."};
319 :
320 0 : WRAPPED_PUPable_decl_template(CauchySecondOrder); // NOLINT
321 0 : explicit CauchySecondOrder(CkMigrateMessage* /*unused*/) {}
322 :
323 0 : CauchySecondOrder(
324 : double max_partially_flat_j0, size_t j0_max_iterations,
325 : double max_cauchy_j0, double j2_tolerance, size_t j2_max_iterations,
326 : double max_partially_flat_j2,
327 : std::unique_ptr<intrp::SpanInterpolator> du_dr_j_interpolator);
328 :
329 0 : CauchySecondOrder() = default;
330 :
331 0 : std::unique_ptr<InitializeJ> get_clone() const override;
332 :
333 1 : std::unique_ptr<intrp::SpanInterpolator> du_dr_j_interpolator()
334 : const override;
335 :
336 : // Per-class tag lists. The flexible dispatch in `InitializeJ<false>` reads
337 : // these via `call_with_dynamic_type` so this generator can request more
338 : // worldtube boundary values than the simpler sibling classes do.
339 0 : using return_tags = tmpl::list<Tags::BondiJ, Tags::CauchyCartesianCoords,
340 : Tags::CauchyAngularCoords>;
341 0 : using argument_tags = tmpl::list<
342 : Tags::BoundaryValue<Tags::BondiJ>, Tags::BoundaryValue<Tags::BondiU>,
343 : Tags::BoundaryValue<Tags::BondiW>, Tags::BoundaryValue<Tags::BondiBeta>,
344 : Tags::BoundaryValue<Tags::BondiQ>,
345 : Tags::BoundaryValue<Tags::Du<Tags::BondiJ>>,
346 : Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>,
347 : Tags::BoundaryValue<Tags::Du<Tags::Dr<Tags::BondiJ>>>,
348 : Tags::BoundaryValue<Tags::Du<Tags::BondiR>>,
349 : Tags::BoundaryValue<Tags::BondiR>, Tags::LMax,
350 : Tags::NumberOfRadialPoints>;
351 :
352 0 : void operator()(
353 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> j,
354 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
355 : gsl::not_null<
356 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
357 : angular_cauchy_coordinates,
358 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_j,
359 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& boundary_u,
360 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_w,
361 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_beta,
362 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& boundary_q,
363 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_du_j,
364 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_dr_j,
365 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_du_dr_j,
366 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_du_r,
367 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& r, size_t l_max,
368 : size_t number_of_radial_points,
369 : gsl::not_null<Parallel::NodeLock*> hdf5_lock) const;
370 :
371 0 : void pup(PUP::er& p) override;
372 :
373 : private:
374 0 : double max_partially_flat_j0_ = std::numeric_limits<double>::signaling_NaN();
375 0 : size_t j0_max_iterations_ = 0;
376 0 : double max_cauchy_j0_ = std::numeric_limits<double>::signaling_NaN();
377 0 : double j2_tolerance_ = std::numeric_limits<double>::signaling_NaN();
378 0 : size_t j2_max_iterations_ = 0;
379 0 : double max_partially_flat_j2_ = std::numeric_limits<double>::signaling_NaN();
380 0 : std::unique_ptr<intrp::SpanInterpolator> du_dr_j_interpolator_;
381 : };
382 : } // namespace Cce::InitializeJ
|
