Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 : #include <ostream>
8 :
9 : #include "DataStructures/TaggedTuple.hpp"
10 : #include "NumericalAlgorithms/Interpolation/CubicSpline.hpp"
11 : #include "Options/Context.hpp"
12 : #include "Options/String.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
15 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
16 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
17 :
18 : /// \cond
19 : namespace PUP {
20 : class er;
21 : } // namespace PUP
22 : namespace Tags {
23 : template <typename Tag>
24 : struct dt;
25 : } // namespace Tags
26 : /// \endcond
27 :
28 : namespace gr::Solutions {
29 :
30 : /*!
31 : * \brief Trumpet Schwarzschild solution in isotropic coordinates
32 : *
33 : * \details
34 : * This solution is a trumpet Schwarzschild black hole in the isotropic
35 : * coordinates. It is a time-independent puncture solution in 1+log slicing.
36 : * The solution cannot be written down analytically in Schwarzschild
37 : * coordinates or isotropic coordinates. Refer to \cite Hannam2008
38 : * for equations and details.
39 : *
40 : * We first set up a source grid in isotropic radial coordinate r, on which
41 : * we compute the lapse and Schwarzschild radial coordinate R. The lapse
42 : * is computed by solving eq. (54)-(56) using the toms748 algorithm, where
43 : * the integrals in eq. (54)-(56) are evaluated by a tanh_sinh integrator.
44 : *
45 : * Eq. (54) is for \f$\alpha < \alpha_s=0.1\f$:
46 : * \f{equation*}{
47 : * r(\alpha) = R(\alpha_s)^{(1/\alpha_s)} \exp \left[ -
48 : * \int_{\alpha}^{\alpha_s} \frac{1}{\bar{\alpha} R(\bar{\alpha})}
49 : * \frac{dR}{d\alpha}(\bar{\alpha}) \: d\bar{\alpha} - C_0 \right],
50 : * \f}
51 : * where eq. (55) gives
52 : * \f{equation*}{
53 : * C_0 = \int_{\alpha_s}^{1} \frac{\ln R(\alpha)}{\alpha^2} d\alpha.
54 : * \f}
55 : * Eq. (56) is for \f$\alpha > \alpha_s=0.1\f$:
56 : * \f{equation*}{
57 : * r(\alpha) = R(\alpha)^{(1/\alpha)} \exp \left[
58 : * \int_{\alpha_s}^{\alpha} \frac{\ln R(\bar{\alpha})}{\bar{\alpha}^2}
59 : * d\bar{\alpha} - C_0 \right].
60 : * \f}
61 : *
62 : * A standard Gauss quadratures will fail due to singular
63 : * behaviors of the integrands near the integration limits. The Schwarzschild R
64 : * is then computed by solving eq. (39) using again the toms748 algorithm.
65 : * \f{equation*}{
66 : * \alpha^2 = 1 - \frac{2M}{R} + \frac{C(n)^2 e^{2 \alpha / n}
67 : * }{R^4}.
68 : * \f}
69 : * Note that for a value of n near 2 (corresponding to standard 1+log slicing),
70 : * eq. (39) has two positive roots for a fixed lapse in [0., 1.).
71 : * See the following plot by Mathematica. We choose
72 : * \image html trumpet_two_branches_illustration.png
73 : * the physical root, i.e. with the correct asymptotic behaviors: R diverges as
74 : * \f$\alpha\f$ tends to 1, and R tends to the smaller solution as \f$\alpha\f$
75 : * tends to 0.
76 :
77 : * The physical root is numerically selected by finding the
78 : * critical lapse and critical Schwarzschild R, below which we tell the toms748
79 : * solver to find a root \f$R\in\f$ [min_schwarzschild_r, crit_schwarzschild_r]
80 : * or else we find a root
81 : * \f$R\in\f$ [crit_schwarzschild_r, max_schwarzschild_r]. min_schwarzschild_r
82 : * is currently selected to be 0., and max_schwarzschild_r is at least as large
83 : * as max_isotropic_r (currently 5000M), the maximum coordinate radius we
84 : * support for this initial data. In case the solver is asked to find a
85 : * Schwarzschild R greater than the max_isotropic_r, we use double the
86 : * asymptotic solution from eq. (39), i.e. \f$4/(1-\alpha^2)\f$, as solver
87 : * upper bound. This latter upper bound is necessary since the integrator
88 : * in eq. (55) needs lapse very close to 1 to converge.
89 : *
90 : * After acquiring the lapse and Schwarzschild R, we can assemble all the 3+1
91 : * quantities in the isotropic coordinates on the source grid.
92 : *
93 : * Since the user supplies grid points in the Cartesian version of the
94 : * isotropic coordinates, we compute a user grid in the isotropic radial
95 : * coordinate based on the Cartesian grid, compute the lapse and Schwarzschild
96 : * R on the source grid, interpolate to the user grid, and then assemble all
97 : * 3+1 quantities on the user grid and transform them back to the Cartesian
98 : * grid.
99 : *
100 : * To insulate our implementation from different mass parameters, we
101 : * nondimensionalize the above process using the black hole mass until
102 : * the final step of computing 3+1 quantities on the Cartesian grid,
103 : * where we restore the correct unit.
104 : *
105 : * The following are input file options that can be specified:
106 : * - Mass
107 : * - N (the parameter n in the slicing condition eq. (36))
108 : *
109 : * \note
110 : * N=2. is strongly suggested, as this gives a stationary solution
111 : * in the standard 1+log slicing. The other values of N near 2.
112 : * should work but have not been tested thoroughly. N outside of
113 : * [2., 3.] has not been tested at all.
114 : *
115 : * Some quantities very close to the puncture (<1.e-4 in isotropic radius)
116 : * may have larger truncation errors. This is expected since some quantities
117 : * such as the determiant of the spatial metric diverges at the puncture.
118 : */
119 1 : class TrumpetSchwarzschild : public MarkAsAnalyticSolution,
120 : public AnalyticSolution<3_st> {
121 : private:
122 : template <typename DataType>
123 : struct IntermediateVars;
124 :
125 : public:
126 0 : static constexpr size_t volume_dim = 3;
127 :
128 0 : struct Mass {
129 0 : using type = double;
130 0 : static constexpr Options::String help = {
131 : "Mass of the Schwarzschild black hole"};
132 0 : static type lower_bound() { return 0.1; };
133 : };
134 :
135 : // currently we have only tested around N=2; for other N eq. (39) may need to
136 : // be solved differently
137 0 : struct N {
138 0 : using type = double;
139 0 : static constexpr Options::String help = {
140 : "Parameter of the trumpet solution family. N=2. gives"
141 : "a stationary solution in standard 1+log slicing."
142 : "An N value outside [2., 3.] has not been tested."};
143 0 : static type lower_bound() { return 0.; };
144 : };
145 :
146 0 : using options = tmpl::list<Mass, N>;
147 0 : static constexpr Options::String help{
148 : "Schwarzschild solution in trumpet isotropic coordinates"};
149 :
150 0 : TrumpetSchwarzschild(double mass, double n,
151 : const Options::Context& context = {});
152 :
153 0 : TrumpetSchwarzschild() = default;
154 0 : TrumpetSchwarzschild(const TrumpetSchwarzschild& /*rhs*/) = default;
155 0 : TrumpetSchwarzschild& operator=(const TrumpetSchwarzschild& /*rhs*/) =
156 : default;
157 0 : TrumpetSchwarzschild(TrumpetSchwarzschild&& /*rhs*/) = default;
158 0 : TrumpetSchwarzschild& operator=(TrumpetSchwarzschild&& /*rhs*/) = default;
159 0 : ~TrumpetSchwarzschild() = default;
160 :
161 0 : explicit TrumpetSchwarzschild(CkMigrateMessage* /*msg*/);
162 :
163 : template <typename DataType>
164 0 : using DerivLapse = ::Tags::deriv<gr::Tags::Lapse<DataType>,
165 : tmpl::size_t<volume_dim>, Frame::Inertial>;
166 : template <typename DataType>
167 0 : using DerivShift = ::Tags::deriv<gr::Tags::Shift<DataType, volume_dim>,
168 : tmpl::size_t<volume_dim>, Frame::Inertial>;
169 : template <typename DataType>
170 0 : using DerivSpatialMetric =
171 : ::Tags::deriv<gr::Tags::SpatialMetric<DataType, volume_dim>,
172 : tmpl::size_t<volume_dim>, Frame::Inertial>;
173 :
174 : template <typename DataType>
175 0 : using tags = tmpl::list<
176 : gr::Tags::Lapse<DataType>, ::Tags::dt<gr::Tags::Lapse<DataType>>,
177 : DerivLapse<DataType>, gr::Tags::Shift<DataType, volume_dim>,
178 : ::Tags::dt<gr::Tags::Shift<DataType, volume_dim>>, DerivShift<DataType>,
179 : gr::Tags::SpatialMetric<DataType, volume_dim>,
180 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, volume_dim>>,
181 : DerivSpatialMetric<DataType>, gr::Tags::SqrtDetSpatialMetric<DataType>,
182 : gr::Tags::ExtrinsicCurvature<DataType, volume_dim>,
183 : gr::Tags::InverseSpatialMetric<DataType, volume_dim>>;
184 :
185 : // The user input Cartesian "x" is in isotropic coordinates
186 : template <typename DataType, typename... Tags>
187 0 : tuples::TaggedTuple<Tags...> variables(
188 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
189 : tmpl::list<Tags...> /*meta*/) const {
190 : const auto& vars =
191 : IntermediateVars<DataType>{mass_, n_, x, t, data_on_source_grid_};
192 : return {get<Tags>(variables(x, t, vars, tmpl::list<Tags>{}))...};
193 : }
194 :
195 : template <typename DataType, typename... Tags>
196 0 : tuples::TaggedTuple<Tags...> variables(
197 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
198 : const IntermediateVars<DataType>& vars,
199 : tmpl::list<Tags...> /*meta*/) const {
200 : static_assert(sizeof...(Tags) > 1,
201 : "Unrecognized tag requested. See the function parameters "
202 : "for the tag.");
203 : return {get<Tags>(variables(x, t, vars, tmpl::list<Tags>{}))...};
204 : }
205 :
206 : // NOLINTNEXTLINE(google-runtime-references)
207 0 : void pup(PUP::er& p);
208 :
209 0 : SPECTRE_ALWAYS_INLINE double mass() const { return mass_; }
210 0 : SPECTRE_ALWAYS_INLINE double n() const { return n_; }
211 :
212 : private:
213 : template <typename DataType>
214 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
215 : double t, const IntermediateVars<DataType>& vars,
216 : tmpl::list<gr::Tags::Lapse<DataType>> /*meta*/) const
217 : -> tuples::TaggedTuple<gr::Tags::Lapse<DataType>>;
218 :
219 : template <typename DataType>
220 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
221 : double t, const IntermediateVars<DataType>& vars,
222 : tmpl::list<::Tags::dt<gr::Tags::Lapse<DataType>>> /*meta*/)
223 : const -> tuples::TaggedTuple<::Tags::dt<gr::Tags::Lapse<DataType>>>;
224 :
225 : template <typename DataType>
226 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
227 : double t, const IntermediateVars<DataType>& vars,
228 : tmpl::list<DerivLapse<DataType>> /*meta*/) const
229 : -> tuples::TaggedTuple<DerivLapse<DataType>>;
230 :
231 : template <typename DataType>
232 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
233 : double t, const IntermediateVars<DataType>& vars,
234 : tmpl::list<gr::Tags::Shift<DataType, volume_dim>> /*meta*/)
235 : const -> tuples::TaggedTuple<gr::Tags::Shift<DataType, volume_dim>>;
236 :
237 : template <typename DataType>
238 0 : auto variables(
239 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
240 : const IntermediateVars<DataType>& vars,
241 : tmpl::list<::Tags::dt<gr::Tags::Shift<DataType, volume_dim>>> /*meta*/)
242 : const
243 : -> tuples::TaggedTuple<::Tags::dt<gr::Tags::Shift<DataType, volume_dim>>>;
244 :
245 : template <typename DataType>
246 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
247 : double t, const IntermediateVars<DataType>& vars,
248 : tmpl::list<DerivShift<DataType>> /*meta*/) const
249 : -> tuples::TaggedTuple<DerivShift<DataType>>;
250 :
251 : template <typename DataType>
252 0 : auto variables(
253 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
254 : const IntermediateVars<DataType>& vars,
255 : tmpl::list<gr::Tags::SpatialMetric<DataType, volume_dim>> /*meta*/) const
256 : -> tuples::TaggedTuple<gr::Tags::SpatialMetric<DataType, volume_dim>>;
257 :
258 : template <typename DataType>
259 0 : auto variables(
260 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
261 : const IntermediateVars<DataType>& vars,
262 : tmpl::list<
263 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, volume_dim>>> /*meta*/)
264 : const -> tuples::TaggedTuple<
265 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, volume_dim>>>;
266 :
267 : template <typename DataType>
268 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& /*x*/,
269 : double /*t*/, const IntermediateVars<DataType>& vars,
270 : tmpl::list<DerivSpatialMetric<DataType>> /*meta*/) const
271 : -> tuples::TaggedTuple<DerivSpatialMetric<DataType>>;
272 :
273 : template <typename DataType>
274 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& /*x*/,
275 : double /*t*/, const IntermediateVars<DataType>& vars,
276 : tmpl::list<gr::Tags::SqrtDetSpatialMetric<DataType>> /*meta*/)
277 : const -> tuples::TaggedTuple<gr::Tags::SqrtDetSpatialMetric<DataType>>;
278 :
279 : template <typename DataType>
280 0 : auto variables(
281 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
282 : const IntermediateVars<DataType>& vars,
283 : tmpl::list<gr::Tags::ExtrinsicCurvature<DataType, volume_dim>> /*meta*/)
284 : const -> tuples::TaggedTuple<
285 : gr::Tags::ExtrinsicCurvature<DataType, volume_dim>>;
286 :
287 : template <typename DataType>
288 0 : auto variables(
289 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
290 : const IntermediateVars<DataType>& vars,
291 : tmpl::list<gr::Tags::InverseSpatialMetric<DataType, volume_dim>> /*meta*/)
292 : const -> tuples::TaggedTuple<
293 : gr::Tags::InverseSpatialMetric<DataType, volume_dim>>;
294 :
295 : template <typename DataType>
296 0 : struct IntermediateVars {
297 0 : IntermediateVars(double mass, double n,
298 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
299 : double /*t*/,
300 : const std::array<DataVector, 2>& data_on_source_grid);
301 : // the lapse, Schwarzschild R, and isotropic r on isotropic grid points
302 0 : DataType lapse_on_user_grid{};
303 0 : DataType schwarzschild_r_on_user_grid{};
304 0 : DataType d_lapse_d_schwarzschild_r_on_user_grid{};
305 0 : DataType isotropic_r_on_user_grid{};
306 :
307 : // cached value to avoid division
308 0 : DataType one_over_schwarzschild_r_on_user_grid{};
309 0 : DataType one_over_isotropic_r_on_user_grid{};
310 0 : double one_over_mass;
311 0 : double one_over_n;
312 : };
313 :
314 0 : double mass_{std::numeric_limits<double>::signaling_NaN()};
315 0 : double n_{std::numeric_limits<double>::signaling_NaN()};
316 :
317 0 : const static tnsr::I<DataVector, 1, Frame::Inertial> source_grid_;
318 0 : std::array<DataVector, 2> data_on_source_grid_;
319 : };
320 :
321 0 : bool operator==(const TrumpetSchwarzschild& lhs,
322 : const TrumpetSchwarzschild& rhs);
323 0 : bool operator!=(const TrumpetSchwarzschild& lhs,
324 : const TrumpetSchwarzschild& rhs);
325 : } // namespace gr::Solutions
|
