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 <memory>
8 : #include <ostream>
9 :
10 : #include "DataStructures/CachedTempBuffer.hpp"
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Elliptic/Systems/Xcts/Tags.hpp"
15 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
16 : #include "Options/String.hpp"
17 : #include "PointwiseFunctions/AnalyticSolutions/Xcts/CommonVariables.hpp"
18 : #include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp"
19 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
20 : #include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp"
21 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
22 : #include "Utilities/ContainerHelpers.hpp"
23 : #include "Utilities/Gsl.hpp"
24 : #include "Utilities/Serialization/CharmPupable.hpp"
25 : #include "Utilities/TMPL.hpp"
26 : #include "Utilities/TaggedTuple.hpp"
27 :
28 : /// \cond
29 : namespace PUP {
30 : class er;
31 : } // namespace PUP
32 : /// \endcond
33 :
34 : namespace Xcts::Solutions {
35 :
36 : /// Various coordinate systems in which to express the Schwarzschild solution
37 1 : enum class SchwarzschildCoordinates {
38 : /*!
39 : * \brief Isotropic Schwarzschild coordinates
40 : *
41 : * These arise from the canonical Schwarzschild coordinates by the radial
42 : * transformation
43 : *
44 : * \f{equation}
45 : * r = \bar{r}\left(1+\frac{M}{2\bar{r}}\right)^2
46 : * \f}
47 : *
48 : * (Eq. (1.61) in \cite BaumgarteShapiro) where \f$r\f$ is the canonical
49 : * Schwarzschild radius, also referred to as "areal" radius because it is
50 : * defined such that spheres with constant \f$r\f$ have the area \f$4\pi
51 : * r^2\f$, and \f$\bar{r}\f$ is the "isotropic" radius. In the isotropic
52 : * radius the Schwarzschild spatial metric is conformally flat:
53 : *
54 : * \f{equation}
55 : * \gamma_{ij}=\psi^4\eta_{ij} \quad \text{with conformal factor} \quad
56 : * \psi=1+\frac{M}{2\bar{r}}
57 : * \f}
58 : *
59 : * (Table 2.1 in \cite BaumgarteShapiro). The lapse in the conformal radius is
60 : *
61 : * \f{equation}
62 : * \alpha=\frac{1-M/(2\bar{r})}{1+M/(2\bar{r})}
63 : * \f}
64 : *
65 : * and the shift vanishes (\f$\beta^i=0\f$) as it does in areal Schwarzschild
66 : * coordinates. The solution also remains maximally sliced, i.e. \f$K=0\f$.
67 : *
68 : * The Schwarzschild horizon in these coordinates is at
69 : * \f$\bar{r}=\frac{M}{2}\f$ due to the radial transformation from \f$r=2M\f$.
70 : */
71 : Isotropic,
72 : /*!
73 : * \brief Painlevé-Gullstrand coordinates
74 : *
75 : * In these coordinates the spatial metric is flat and the lapse is trivial,
76 : * but contrary to (isotropic) Schwarzschild coordinates the shift is
77 : * nontrivial,
78 : *
79 : * \begin{align}
80 : * \gamma_{ij} &= \eta_{ij} \\
81 : * \alpha &= 1 \\
82 : * \beta^i &= \sqrt{\frac{2M}{r}} \frac{x^i}{r} \\
83 : * K &= \frac{3}{2}\sqrt{\frac{2M}{r^3}}
84 : * \end{align}
85 : *
86 : * (Table 2.1 in \cite BaumgarteShapiro).
87 : */
88 : PainleveGullstrand,
89 : /*!
90 : * \brief Isotropic Kerr-Schild coordinates
91 : *
92 : * Kerr-Schild coordinates with a radial transformation such that the spatial
93 : * metric is conformally flat.
94 : *
95 : * The Schwarzschild spacetime in canonical (areal) Kerr-Schild coordinates is
96 : *
97 : * \begin{align}
98 : * \gamma_{ij} &= \eta_{ij} + \frac{2M}{r}\frac{x^i x^j}{r^2} \\
99 : * \alpha &= \sqrt{1 + \frac{2M}{r}}^{-1} \\
100 : * \beta^i &= \frac{2M\alpha^2}{r} \frac{x^i}{r} \\
101 : * K &= \frac{2M\alpha^3}{r^2} \left(1 + \frac{3M}{r}\right)
102 : * \text{.}
103 : * \end{align}
104 : *
105 : * (Table 2.1 in \cite BaumgarteShapiro). Since the Schwarzschild spacetime is
106 : * spherically symmetric we can transform to a radial coordinate $\bar{r}$ in
107 : * which it is conformally flat (see, e.g., Sec. 7.4.1 in \cite Pfeiffer2005zm
108 : * for details):
109 : *
110 : * \begin{equation}
111 : * {}^{(3)}\mathrm{d}s^2
112 : * = \left(1 + \frac{2M}{r}\right)\mathrm{d}r^2 + r^2 \mathrm{d}\Omega^2
113 : * = \psi^4 \left(\mathrm{d}\bar{r}^2 +
114 : * \bar{r}^2 \mathrm{d}\Omega^2\right)
115 : * \end{equation}
116 : *
117 : * Therefore, the conformal factor is $\psi^2 = r / \bar{r}$ and
118 : *
119 : * \begin{equation}
120 : * \frac{\mathrm{d}\bar{r}}{\mathrm{d}r}
121 : * = \frac{\bar{r}}{r} \sqrt{1 + \frac{2M}{r}}
122 : * = \frac{\bar{r}}{r} \frac{1}{\alpha}
123 : * \text{,}
124 : * \end{equation}
125 : *
126 : * which has the solution
127 : *
128 : * \begin{equation}
129 : * \bar{r} = \frac{r}{4} \left(1 + \sqrt{1 + \frac{2M}{r}}\right)^2
130 : * e^{2 - 2\sqrt{1 + 2M / r}}
131 : * \end{equation}
132 : *
133 : * when we impose $\bar{r} \rightarrow r$ as $r \rightarrow \infty$. We can
134 : * invert this transformation law with a numerical root find to obtain the
135 : * areal radius $r$ for any isotropic radius $\bar{r}$.
136 : *
137 : * In the isotropic radial coordinate $\bar{r}$ the solution is then:
138 : *
139 : * \begin{align}
140 : * \gamma_{ij} &= \psi^4 \eta_{ij} \\
141 : * \psi &= \sqrt{\frac{r}{\bar{r}}}
142 : * = \frac{2e^{\sqrt{1 + 2M / r} - 1}}{1 + \sqrt{1 + 2M / r}} \\
143 : * \alpha &= \sqrt{1 + \frac{2M}{r}}^{-1} \\
144 : * \beta^i
145 : * &= \frac{\mathrm{d}\bar{r}}{\mathrm{d}r} \beta^r \frac{x^i}{\bar{r}}
146 : * = \frac{2M\alpha}{r^2} x^i \\
147 : * K &= \frac{2M\alpha^3}{r^2} \left(1 + \frac{3M}{r}\right)
148 : * \end{align}
149 : *
150 : * Here, $x^i$ are the (isotropic) Cartesian coordinates from which we compute
151 : * the isotropic radius $\bar{r}$, $r$ is the areal radius we can obtain from
152 : * the isotropic radius by a root find, and $\beta^r$ is the magnitude of the
153 : * shift in areal coordinates, as given above.
154 : *
155 : * The horizon in these coordinates is at (Eq. (7.37) in
156 : * \cite Pfeiffer2005zm):
157 : *
158 : * \begin{equation}
159 : * \bar{r}_\mathrm{AH} / M \approx 1.2727410334221052
160 : * \end{equation}
161 : */
162 : KerrSchildIsotropic,
163 : };
164 :
165 0 : std::ostream& operator<<(std::ostream& os, SchwarzschildCoordinates coords);
166 :
167 : } // namespace Xcts::Solutions
168 :
169 : template <>
170 0 : struct Options::create_from_yaml<Xcts::Solutions::SchwarzschildCoordinates> {
171 : template <typename Metavariables>
172 0 : static Xcts::Solutions::SchwarzschildCoordinates create(
173 : const Options::Option& options) {
174 : return create<void>(options);
175 : }
176 : };
177 :
178 : template <>
179 0 : Xcts::Solutions::SchwarzschildCoordinates
180 : Options::create_from_yaml<Xcts::Solutions::SchwarzschildCoordinates>::create<
181 : void>(const Options::Option& options);
182 :
183 : namespace Xcts::Solutions {
184 :
185 : namespace detail {
186 :
187 : struct SchwarzschildImpl {
188 : struct Mass {
189 : using type = double;
190 : static constexpr Options::String help = "Mass parameter M";
191 : };
192 :
193 : struct CoordinateSystem {
194 : static std::string name() { return "Coordinates"; }
195 : using type = SchwarzschildCoordinates;
196 : static constexpr Options::String help =
197 : "The coordinate system used to describe the solution";
198 : };
199 :
200 : using options = tmpl::list<Mass, CoordinateSystem>;
201 : static constexpr Options::String help{
202 : "Schwarzschild spacetime in general relativity"};
203 :
204 : SchwarzschildImpl() = default;
205 : SchwarzschildImpl(const SchwarzschildImpl&) = default;
206 : SchwarzschildImpl& operator=(const SchwarzschildImpl&) = default;
207 : SchwarzschildImpl(SchwarzschildImpl&&) = default;
208 : SchwarzschildImpl& operator=(SchwarzschildImpl&&) = default;
209 : ~SchwarzschildImpl() = default;
210 :
211 : explicit SchwarzschildImpl(double mass,
212 : SchwarzschildCoordinates coordinate_system);
213 :
214 : /// The mass parameter \f$M\f$.
215 : double mass() const;
216 :
217 : SchwarzschildCoordinates coordinate_system() const;
218 :
219 : /// The radius of the Schwarzschild horizon in the given coordinates.
220 : double radius_at_horizon() const;
221 :
222 : // NOLINTNEXTLINE(google-runtime-references)
223 : void pup(PUP::er& p);
224 :
225 : protected:
226 : double mass_{std::numeric_limits<double>::signaling_NaN()};
227 : SchwarzschildCoordinates coordinate_system_{};
228 : };
229 :
230 : bool operator==(const SchwarzschildImpl& lhs, const SchwarzschildImpl& rhs);
231 :
232 : bool operator!=(const SchwarzschildImpl& lhs, const SchwarzschildImpl& rhs);
233 :
234 : namespace Tags {
235 : template <typename DataType>
236 : struct Radius : db::SimpleTag {
237 : using type = Scalar<DataType>;
238 : };
239 : template <typename DataType>
240 : struct ArealRadius : db::SimpleTag {
241 : using type = Scalar<DataType>;
242 : };
243 : } // namespace Tags
244 :
245 : template <typename DataType>
246 : using SchwarzschildVariablesCache =
247 : cached_temp_buffer_from_typelist<tmpl::push_front<
248 : common_tags<DataType>, detail::Tags::Radius<DataType>,
249 : detail::Tags::ArealRadius<DataType>,
250 : ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>,
251 : Frame::Inertial>,
252 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
253 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
254 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, 3>, 0>>>;
255 :
256 : template <typename DataType>
257 : struct SchwarzschildVariables
258 : : CommonVariables<DataType, SchwarzschildVariablesCache<DataType>> {
259 : static constexpr size_t Dim = 3;
260 : static constexpr int ConformalMatterScale = 0;
261 : using Cache = SchwarzschildVariablesCache<DataType>;
262 : using Base = CommonVariables<DataType, SchwarzschildVariablesCache<DataType>>;
263 : using Base::operator();
264 :
265 : SchwarzschildVariables(
266 : std::optional<std::reference_wrapper<const Mesh<Dim>>> local_mesh,
267 : std::optional<std::reference_wrapper<const InverseJacobian<
268 : DataType, Dim, Frame::ElementLogical, Frame::Inertial>>>
269 : local_inv_jacobian,
270 : const tnsr::I<DataType, 3>& local_x, const double local_mass,
271 : const SchwarzschildCoordinates local_coordinate_system)
272 : : Base(std::move(local_mesh), std::move(local_inv_jacobian)),
273 : x(local_x),
274 : mass(local_mass),
275 : coordinate_system(local_coordinate_system) {}
276 :
277 : const tnsr::I<DataType, 3>& x;
278 : double mass;
279 : SchwarzschildCoordinates coordinate_system;
280 :
281 : void operator()(gsl::not_null<Scalar<DataType>*> radius,
282 : gsl::not_null<Cache*> cache,
283 : detail::Tags::Radius<DataType> /*meta*/) const;
284 : void operator()(gsl::not_null<Scalar<DataType>*> areal_radius,
285 : gsl::not_null<Cache*> cache,
286 : detail::Tags::ArealRadius<DataType> /*meta*/) const;
287 : void operator()(
288 : gsl::not_null<tnsr::ii<DataType, 3>*> conformal_metric,
289 : gsl::not_null<Cache*> cache,
290 : Xcts::Tags::ConformalMetric<DataType, 3, Frame::Inertial> /*meta*/)
291 : const override;
292 : void operator()(
293 : gsl::not_null<tnsr::II<DataType, 3>*> inv_conformal_metric,
294 : gsl::not_null<Cache*> cache,
295 : Xcts::Tags::InverseConformalMetric<DataType, 3, Frame::Inertial> /*meta*/)
296 : const override;
297 : void operator()(
298 : gsl::not_null<tnsr::ijj<DataType, 3>*> deriv_conformal_metric,
299 : gsl::not_null<Cache*> cache,
300 : ::Tags::deriv<Xcts::Tags::ConformalMetric<DataType, 3, Frame::Inertial>,
301 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
302 : void operator()(
303 : gsl::not_null<Scalar<DataType>*> trace_extrinsic_curvature,
304 : gsl::not_null<Cache*> cache,
305 : gr::Tags::TraceExtrinsicCurvature<DataType> /*meta*/) const override;
306 : void operator()(
307 : gsl::not_null<tnsr::i<DataType, 3>*> trace_extrinsic_curvature_gradient,
308 : gsl::not_null<Cache*> cache,
309 : ::Tags::deriv<gr::Tags::TraceExtrinsicCurvature<DataType>,
310 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
311 : void operator()(
312 : gsl::not_null<Scalar<DataType>*> dt_trace_extrinsic_curvature,
313 : gsl::not_null<Cache*> cache,
314 : ::Tags::dt<gr::Tags::TraceExtrinsicCurvature<DataType>> /*meta*/)
315 : const override;
316 : void operator()(
317 : gsl::not_null<Scalar<DataType>*> conformal_factor_minus_one,
318 : gsl::not_null<Cache*> cache,
319 : Xcts::Tags::ConformalFactorMinusOne<DataType> /*meta*/) const override;
320 : void operator()(
321 : gsl::not_null<Scalar<DataType>*> conformal_factor,
322 : gsl::not_null<Cache*> cache,
323 : Xcts::Tags::ConformalFactor<DataType> /*meta*/) const override;
324 : void operator()(
325 : gsl::not_null<tnsr::i<DataType, 3>*> conformal_factor_gradient,
326 : gsl::not_null<Cache*> cache,
327 : ::Tags::deriv<Xcts::Tags::ConformalFactorMinusOne<DataType>,
328 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
329 : void operator()(gsl::not_null<Scalar<DataType>*> lapse,
330 : gsl::not_null<Cache*> cache,
331 : gr::Tags::Lapse<DataType> /*meta*/) const override;
332 : void operator()(gsl::not_null<tnsr::i<DataType, 3>*> deriv_lapse,
333 : gsl::not_null<Cache*> cache,
334 : ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>,
335 : Frame::Inertial> /*meta*/) const;
336 : void operator()(
337 : gsl::not_null<Scalar<DataType>*> lapse_times_conformal_factor_minus_one,
338 : gsl::not_null<Cache*> cache,
339 : Xcts::Tags::LapseTimesConformalFactorMinusOne<DataType> /*meta*/)
340 : const override;
341 : void operator()(
342 : gsl::not_null<Scalar<DataType>*> lapse_times_conformal_factor,
343 : gsl::not_null<Cache*> cache,
344 : Xcts::Tags::LapseTimesConformalFactor<DataType> /*meta*/) const override;
345 : void operator()(
346 : gsl::not_null<tnsr::i<DataType, 3>*>
347 : lapse_times_conformal_factor_gradient,
348 : gsl::not_null<Cache*> cache,
349 : ::Tags::deriv<Xcts::Tags::LapseTimesConformalFactorMinusOne<DataType>,
350 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
351 : void operator()(
352 : gsl::not_null<tnsr::I<DataType, 3>*> shift_background,
353 : gsl::not_null<Cache*> cache,
354 : Xcts::Tags::ShiftBackground<DataType, 3, Frame::Inertial> /*meta*/)
355 : const override;
356 : void operator()(gsl::not_null<tnsr::II<DataType, 3, Frame::Inertial>*>
357 : longitudinal_shift_background_minus_dt_conformal_metric,
358 : gsl::not_null<Cache*> cache,
359 : Xcts::Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
360 : DataType, 3, Frame::Inertial> /*meta*/) const override;
361 : void operator()(
362 : gsl::not_null<tnsr::I<DataType, 3>*> shift_excess,
363 : gsl::not_null<Cache*> cache,
364 : Xcts::Tags::ShiftExcess<DataType, 3, Frame::Inertial> /*meta*/)
365 : const override;
366 : void operator()(
367 : gsl::not_null<tnsr::iJ<DataType, 3>*> deriv_shift_excess,
368 : gsl::not_null<Cache*> cache,
369 : ::Tags::deriv<Xcts::Tags::ShiftExcess<DataType, 3, Frame::Inertial>,
370 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
371 : void operator()(
372 : gsl::not_null<tnsr::ii<DataType, 3>*> extrinsic_curvature,
373 : gsl::not_null<Cache*> cache,
374 : gr::Tags::ExtrinsicCurvature<DataType, 3> /*meta*/) const override;
375 : void operator()(gsl::not_null<Scalar<DataType>*> energy_density,
376 : gsl::not_null<Cache*> cache,
377 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>,
378 : ConformalMatterScale> /*meta*/) const;
379 : void operator()(gsl::not_null<Scalar<DataType>*> stress_trace,
380 : gsl::not_null<Cache*> cache,
381 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>,
382 : ConformalMatterScale> /*meta*/) const;
383 : void operator()(gsl::not_null<tnsr::I<DataType, 3>*> momentum_density,
384 : gsl::not_null<Cache*> cache,
385 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, 3>,
386 : ConformalMatterScale> /*meta*/) const;
387 : };
388 :
389 : } // namespace detail
390 :
391 : /*!
392 : * \brief Schwarzschild spacetime in general relativity
393 : *
394 : * This class implements the Schwarzschild solution with mass parameter
395 : * \f$M\f$ in various coordinate systems. See the entries of the
396 : * `Xcts::Solutions::SchwarzschildCoordinates` enum for the available coordinate
397 : * systems and for the solution variables in the respective coordinates.
398 : */
399 1 : class Schwarzschild : public elliptic::analytic_data::AnalyticSolution,
400 : public detail::SchwarzschildImpl {
401 : public:
402 0 : Schwarzschild() = default;
403 0 : Schwarzschild(const Schwarzschild&) = default;
404 0 : Schwarzschild& operator=(const Schwarzschild&) = default;
405 0 : Schwarzschild(Schwarzschild&&) = default;
406 0 : Schwarzschild& operator=(Schwarzschild&&) = default;
407 0 : ~Schwarzschild() = default;
408 :
409 : using SchwarzschildImpl::SchwarzschildImpl;
410 :
411 : /// \cond
412 : explicit Schwarzschild(CkMigrateMessage* m)
413 : : elliptic::analytic_data::AnalyticSolution(m) {}
414 : using PUP::able::register_constructor;
415 : WRAPPED_PUPable_decl_template(Schwarzschild);
416 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
417 : const override {
418 : return std::make_unique<Schwarzschild>(*this);
419 : }
420 : /// \endcond
421 :
422 : template <typename DataType, typename... RequestedTags>
423 0 : tuples::TaggedTuple<RequestedTags...> variables(
424 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
425 : tmpl::list<RequestedTags...> /*meta*/) const {
426 : return variables_impl<DataType>(x, std::nullopt, std::nullopt,
427 : tmpl::list<RequestedTags...>{});
428 : }
429 :
430 : template <typename... RequestedTags>
431 0 : tuples::TaggedTuple<RequestedTags...> variables(
432 : const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& mesh,
433 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
434 : Frame::Inertial>& inv_jacobian,
435 : tmpl::list<RequestedTags...> /*meta*/) const {
436 : return variables_impl<DataVector>(x, mesh, inv_jacobian,
437 : tmpl::list<RequestedTags...>{});
438 : }
439 :
440 0 : void pup(PUP::er& p) override {
441 : elliptic::analytic_data::AnalyticSolution::pup(p);
442 : detail::SchwarzschildImpl::pup(p);
443 : }
444 :
445 : private:
446 : template <typename DataType, typename... RequestedTags>
447 0 : tuples::TaggedTuple<RequestedTags...> variables_impl(
448 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
449 : std::optional<std::reference_wrapper<const Mesh<3>>> mesh,
450 : std::optional<std::reference_wrapper<const InverseJacobian<
451 : DataType, 3, Frame::ElementLogical, Frame::Inertial>>>
452 : inv_jacobian,
453 : tmpl::list<RequestedTags...> /*meta*/) const {
454 : using VarsComputer = detail::SchwarzschildVariables<DataType>;
455 : typename VarsComputer::Cache cache{get_size(*x.begin())};
456 : const VarsComputer computer{std::move(mesh), std::move(inv_jacobian), x,
457 : mass_, coordinate_system_};
458 : const auto get_var = [&cache, &computer, &x](auto tag_v) {
459 : using tag = std::decay_t<decltype(tag_v)>;
460 : if constexpr (tmpl::list_contains_v<hydro_tags<DataType>, tag>) {
461 : (void)cache;
462 : (void)computer;
463 : return get<tag>(Flatness{}.variables(x, tmpl::list<tag>{}));
464 : } else {
465 : (void)x;
466 : return cache.get_var(computer, tag{});
467 : }
468 : };
469 : return {get_var(RequestedTags{})...};
470 : }
471 : };
472 :
473 : } // namespace Xcts::Solutions
|
