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 "Utilities/Gsl.hpp"
15 : #include "Utilities/Serialization/CharmPupable.hpp"
16 : #include "Utilities/TMPL.hpp"
17 :
18 : /// \cond
19 : class ComplexDataVector;
20 : /// \endcond
21 :
22 1 : namespace Cce::InitializeJ {
23 :
24 : /*!
25 : * \brief Initialize \f$J\f$ on the first hypersurface using a second-order
26 : * matching at the worldtube.
27 : *
28 : * \details The volume \f$J\f$ is built from the worldtube values of
29 : * \f$J\f$, \f$\partial_r J\f$, and \f$\partial_y^2 J\f$ computed from the
30 : * H hypersurface equation. The remaining angular coordinates are determined
31 : * iteratively to ensure asymptotic flatness. The angular solve can eliminate
32 : * \f$J\f$ at scri+ only through a well-behaved alteration of the spherical
33 : * mesh, so it tolerates only a small asymptotic \f$J\f$; the initialization
34 : * aborts if the asymptotic \f$J\f$ in Cauchy coordinates, or the deviation at
35 : * any iteration of the solve, exceeds `MaxAngularSolveError`. As a further
36 : * safeguard, the initialization aborts if the second radial derivative of
37 : * \f$J\f$ at scri+ of the final solution exceeds `MaxScriSecondDerivative`.
38 : */
39 1 : struct CauchySecondOrder : InitializeJ<false> {
40 0 : struct AngularCoordinateTolerance {
41 0 : using type = double;
42 0 : static std::string name() { return "AngularCoordTolerance"; }
43 0 : static constexpr Options::String help = {
44 : "Tolerance of initial angular coordinates for CCE"};
45 0 : static type lower_bound() { return 1.0e-14; }
46 0 : static type upper_bound() { return 1.0e-3; }
47 0 : static type suggested_value() { return 1.0e-12; }
48 : };
49 :
50 0 : struct MaxIterations {
51 0 : using type = size_t;
52 0 : static constexpr Options::String help = {
53 : "Number of linearized inversion iterations."};
54 0 : static type lower_bound() { return 10; }
55 0 : static type upper_bound() { return 1000; }
56 0 : static type suggested_value() { return 300; }
57 : };
58 :
59 0 : struct RequireConvergence {
60 0 : using type = bool;
61 0 : static constexpr Options::String help = {
62 : "If true, initialization will error if it hits MaxIterations"};
63 0 : static type suggested_value() { return true; }
64 : };
65 :
66 0 : struct MaxAngularSolveError {
67 0 : using type = double;
68 0 : static constexpr Options::String help = {
69 : "Largest deviation of J from zero at scri+ that the iterative angular "
70 : "solve is permitted to encounter. Initialization aborts if the "
71 : "asymptotic J in Cauchy coordinates exceeds this value before the "
72 : "solve, or if any iteration of the solve exceeds it. Raise this to "
73 : "attempt initialization from worldtube data with a larger asymptotic "
74 : "strain, at the risk of a poorly behaved angular coordinate map."};
75 0 : static type lower_bound() { return 1.0e-14; }
76 0 : static type upper_bound() { return 1.0e2; }
77 0 : static type suggested_value() { return 1.0e-1; }
78 : };
79 :
80 0 : struct MaxScriSecondDerivative {
81 0 : using type = double;
82 0 : static constexpr Options::String help = {
83 : "Abort initialization if the largest second radial derivative of J at "
84 : "scri+ of the final initial data exceeds this threshold. The "
85 : "second-order construction drives this derivative to (near) zero, so a "
86 : "large value indicates a poorly matched solution. Set to a large value "
87 : "to effectively disable the check."};
88 0 : static type lower_bound() { return 1.0e-14; }
89 0 : static type upper_bound() { return 1.0e2; }
90 0 : static type suggested_value() { return 1.0e-8; }
91 : };
92 :
93 0 : using options =
94 : tmpl::list<AngularCoordinateTolerance, MaxIterations, RequireConvergence,
95 : MaxAngularSolveError, MaxScriSecondDerivative>;
96 0 : static constexpr Options::String help = {
97 : "Second-order initial data generator for the Cauchy CCE evolution."};
98 :
99 0 : WRAPPED_PUPable_decl_template(CauchySecondOrder); // NOLINT
100 0 : explicit CauchySecondOrder(CkMigrateMessage* /*unused*/) {}
101 :
102 0 : CauchySecondOrder(double angular_coordinate_tolerance, size_t max_iterations,
103 : bool require_convergence, double max_angular_solve_error,
104 : double max_scri_second_derivative);
105 :
106 0 : CauchySecondOrder() = default;
107 :
108 0 : std::unique_ptr<InitializeJ> get_clone() const override;
109 :
110 : // Per-class tag lists. The flexible dispatch in `InitializeJ<false>` reads
111 : // these via `call_with_dynamic_type` so this generator can request more
112 : // worldtube boundary values than the simpler sibling classes do.
113 0 : using return_tags = tmpl::list<Tags::BondiJ, Tags::CauchyCartesianCoords,
114 : Tags::CauchyAngularCoords>;
115 0 : using argument_tags = tmpl::list<
116 : Tags::BoundaryValue<Tags::BondiJ>, Tags::BoundaryValue<Tags::BondiU>,
117 : Tags::BoundaryValue<Tags::BondiW>, Tags::BoundaryValue<Tags::BondiBeta>,
118 : Tags::BoundaryValue<Tags::BondiQ>,
119 : Tags::BoundaryValue<Tags::Du<Tags::BondiJ>>,
120 : Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>,
121 : Tags::BoundaryValue<Tags::Du<Tags::Dr<Tags::BondiJ>>>,
122 : Tags::BoundaryValue<Tags::Du<Tags::BondiR>>,
123 : Tags::BoundaryValue<Tags::BondiR>, Tags::LMax,
124 : Tags::NumberOfRadialPoints>;
125 :
126 0 : void operator()(
127 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> j,
128 : gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_cauchy_coordinates,
129 : gsl::not_null<
130 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*>
131 : angular_cauchy_coordinates,
132 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_j,
133 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& boundary_u,
134 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_w,
135 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_beta,
136 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& boundary_q,
137 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_du_j,
138 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_dr_j,
139 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& boundary_du_dr_j,
140 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_du_r,
141 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& r, size_t l_max,
142 : size_t number_of_radial_points,
143 : gsl::not_null<Parallel::NodeLock*> hdf5_lock) const;
144 :
145 0 : void pup(PUP::er& p) override;
146 :
147 : private:
148 0 : bool require_convergence_ = true;
149 0 : double angular_coordinate_tolerance_ =
150 : std::numeric_limits<double>::signaling_NaN();
151 0 : size_t max_iterations_ = 0;
152 0 : double max_angular_solve_error_ =
153 : std::numeric_limits<double>::signaling_NaN();
154 0 : double max_scri_second_derivative_ =
155 : std::numeric_limits<double>::signaling_NaN();
156 : };
157 : } // namespace Cce::InitializeJ
|
