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 : * \brief Maximal Isotropic (Horizon Penetrating) Schwarzschild coordinates
165 : *
166 : * Schwarzschild coordinates with a radial transformation such that the radius
167 : * is isotropic and the coordinates are horizon penetrating.
168 : *
169 : * These arise from first choosing a family of time-independent, maximal
170 : * slicings of the Schwarzschild spacetime and a slicing condition that give a
171 : * unique solution with a limiting surface at \f$R=3M/2\f$ and horizon at
172 : * \f$R=2M\f$ \cite Estabrook1973ue. The latter is then changed by the radial
173 : * transformation \cite Baumgarte2007ht
174 : *
175 : * \f{equation}
176 : * r = \frac{\left[2R + M + \sqrt{4R^2 + 4MR + 3M^2}\right]}{4}
177 : * \times \left[
178 : * \frac{(4 + 3\sqrt{2})(2R - 3M)}{8R + 6M + 3\sqrt{8R^2 + 8MR + 6M^2}}
179 : * \right]^{1/\sqrt{2}}
180 : * \f}
181 : *
182 : * where \f$R\f$ is the canonical
183 : * Schwarzschild radius, also referred to as "areal" radius because it is
184 : * defined such that spheres with constant \f$R\f$ have the area \f$4\pi
185 : * R^2\f$, and \f$r\f$ is the "isotropic" radius. In the isotropic
186 : * radius the Schwarzschild spatial metric is conformally flat:
187 : *
188 : * \f{equation}
189 : * \gamma_{ij}=\psi^4\eta_{ij} \quad \text{with conformal factor} \quad
190 : * \psi=\sqrt{\frac{R}{r}}.
191 : * \f}
192 : *
193 : * The lapse in the conformal radius is
194 : *
195 : * \f{equation}
196 : * \alpha = \left(1 - \frac{2M}{R} + \frac{27M^4}{16R^4} \right)^{1/2}
197 : * \f}
198 : *
199 : * and the shift is given by
200 : *
201 : * \f{equation}
202 : * \beta^r = \frac{3 \sqrt{3} M^2}{4} \frac{r}{R^3}
203 : * \f}
204 : *
205 : * The solution remains maximally sliced, i.e. \f$K=0\f$. And the horizon in
206 : * these coordinates is at \f$r\approx 0.7793271080557972 M\f$ due to the
207 : * radial transformation from \f$R=2M\f$.
208 : *
209 : */
210 : MaximalIsotropic,
211 : };
212 :
213 0 : std::ostream& operator<<(std::ostream& os, SchwarzschildCoordinates coords);
214 :
215 : } // namespace Xcts::Solutions
216 :
217 : template <>
218 0 : struct Options::create_from_yaml<Xcts::Solutions::SchwarzschildCoordinates> {
219 : template <typename Metavariables>
220 0 : static Xcts::Solutions::SchwarzschildCoordinates create(
221 : const Options::Option& options) {
222 : return create<void>(options);
223 : }
224 : };
225 :
226 : template <>
227 0 : Xcts::Solutions::SchwarzschildCoordinates
228 : Options::create_from_yaml<Xcts::Solutions::SchwarzschildCoordinates>::create<
229 : void>(const Options::Option& options);
230 :
231 : namespace Xcts::Solutions {
232 :
233 : namespace detail {
234 :
235 : struct SchwarzschildImpl {
236 : struct Mass {
237 : using type = double;
238 : static constexpr Options::String help = "Mass parameter M";
239 : };
240 :
241 : struct CoordinateSystem {
242 : static std::string name() { return "Coordinates"; }
243 : using type = SchwarzschildCoordinates;
244 : static constexpr Options::String help =
245 : "The coordinate system used to describe the solution";
246 : };
247 :
248 : using options = tmpl::list<Mass, CoordinateSystem>;
249 : static constexpr Options::String help{
250 : "Schwarzschild spacetime in general relativity"};
251 :
252 : SchwarzschildImpl() = default;
253 : SchwarzschildImpl(const SchwarzschildImpl&) = default;
254 : SchwarzschildImpl& operator=(const SchwarzschildImpl&) = default;
255 : SchwarzschildImpl(SchwarzschildImpl&&) = default;
256 : SchwarzschildImpl& operator=(SchwarzschildImpl&&) = default;
257 : ~SchwarzschildImpl() = default;
258 :
259 : explicit SchwarzschildImpl(double mass,
260 : SchwarzschildCoordinates coordinate_system);
261 :
262 : /// The mass parameter \f$M\f$.
263 : double mass() const;
264 :
265 : SchwarzschildCoordinates coordinate_system() const;
266 :
267 : /// The radius of the Schwarzschild horizon in the given coordinates.
268 : double radius_at_horizon() const;
269 :
270 : // NOLINTNEXTLINE(google-runtime-references)
271 : void pup(PUP::er& p);
272 :
273 : protected:
274 : double mass_{std::numeric_limits<double>::signaling_NaN()};
275 : SchwarzschildCoordinates coordinate_system_{};
276 : };
277 :
278 : bool operator==(const SchwarzschildImpl& lhs, const SchwarzschildImpl& rhs);
279 :
280 : bool operator!=(const SchwarzschildImpl& lhs, const SchwarzschildImpl& rhs);
281 :
282 : namespace Tags {
283 : template <typename DataType>
284 : struct Radius : db::SimpleTag {
285 : using type = Scalar<DataType>;
286 : };
287 : template <typename DataType>
288 : struct ArealRadius : db::SimpleTag {
289 : using type = Scalar<DataType>;
290 : };
291 : } // namespace Tags
292 :
293 : template <typename DataType>
294 : using SchwarzschildVariablesCache =
295 : cached_temp_buffer_from_typelist<tmpl::push_front<
296 : common_tags<DataType>, detail::Tags::Radius<DataType>,
297 : detail::Tags::ArealRadius<DataType>,
298 : ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>,
299 : Frame::Inertial>,
300 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
301 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
302 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, 3>, 0>>>;
303 :
304 : template <typename DataType>
305 : struct SchwarzschildVariables
306 : : CommonVariables<DataType, SchwarzschildVariablesCache<DataType>> {
307 : static constexpr size_t Dim = 3;
308 : static constexpr int ConformalMatterScale = 0;
309 : using Cache = SchwarzschildVariablesCache<DataType>;
310 : using Base = CommonVariables<DataType, SchwarzschildVariablesCache<DataType>>;
311 : using Base::operator();
312 :
313 : SchwarzschildVariables(
314 : std::optional<std::reference_wrapper<const Mesh<Dim>>> local_mesh,
315 : std::optional<std::reference_wrapper<const InverseJacobian<
316 : DataType, Dim, Frame::ElementLogical, Frame::Inertial>>>
317 : local_inv_jacobian,
318 : const tnsr::I<DataType, 3>& local_x, const double local_mass,
319 : const SchwarzschildCoordinates local_coordinate_system)
320 : : Base(std::move(local_mesh), std::move(local_inv_jacobian)),
321 : x(local_x),
322 : mass(local_mass),
323 : coordinate_system(local_coordinate_system) {}
324 :
325 : const tnsr::I<DataType, 3>& x;
326 : double mass;
327 : SchwarzschildCoordinates coordinate_system;
328 :
329 : void operator()(gsl::not_null<Scalar<DataType>*> radius,
330 : gsl::not_null<Cache*> cache,
331 : detail::Tags::Radius<DataType> /*meta*/) const;
332 : void operator()(gsl::not_null<Scalar<DataType>*> areal_radius,
333 : gsl::not_null<Cache*> cache,
334 : detail::Tags::ArealRadius<DataType> /*meta*/) const;
335 : void operator()(
336 : gsl::not_null<tnsr::ii<DataType, 3>*> conformal_metric,
337 : gsl::not_null<Cache*> cache,
338 : Xcts::Tags::ConformalMetric<DataType, 3, Frame::Inertial> /*meta*/)
339 : const override;
340 : void operator()(
341 : gsl::not_null<tnsr::II<DataType, 3>*> inv_conformal_metric,
342 : gsl::not_null<Cache*> cache,
343 : Xcts::Tags::InverseConformalMetric<DataType, 3, Frame::Inertial> /*meta*/)
344 : const override;
345 : void operator()(
346 : gsl::not_null<tnsr::ijj<DataType, 3>*> deriv_conformal_metric,
347 : gsl::not_null<Cache*> cache,
348 : ::Tags::deriv<Xcts::Tags::ConformalMetric<DataType, 3, Frame::Inertial>,
349 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
350 : void operator()(
351 : gsl::not_null<Scalar<DataType>*> trace_extrinsic_curvature,
352 : gsl::not_null<Cache*> cache,
353 : gr::Tags::TraceExtrinsicCurvature<DataType> /*meta*/) const override;
354 : void operator()(
355 : gsl::not_null<tnsr::i<DataType, 3>*> trace_extrinsic_curvature_gradient,
356 : gsl::not_null<Cache*> cache,
357 : ::Tags::deriv<gr::Tags::TraceExtrinsicCurvature<DataType>,
358 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
359 : void operator()(
360 : gsl::not_null<Scalar<DataType>*> dt_trace_extrinsic_curvature,
361 : gsl::not_null<Cache*> cache,
362 : ::Tags::dt<gr::Tags::TraceExtrinsicCurvature<DataType>> /*meta*/)
363 : const override;
364 : void operator()(
365 : gsl::not_null<Scalar<DataType>*> conformal_factor_minus_one,
366 : gsl::not_null<Cache*> cache,
367 : Xcts::Tags::ConformalFactorMinusOne<DataType> /*meta*/) const override;
368 : void operator()(
369 : gsl::not_null<Scalar<DataType>*> conformal_factor,
370 : gsl::not_null<Cache*> cache,
371 : Xcts::Tags::ConformalFactor<DataType> /*meta*/) const override;
372 : void operator()(
373 : gsl::not_null<tnsr::i<DataType, 3>*> conformal_factor_gradient,
374 : gsl::not_null<Cache*> cache,
375 : ::Tags::deriv<Xcts::Tags::ConformalFactorMinusOne<DataType>,
376 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
377 : void operator()(gsl::not_null<Scalar<DataType>*> lapse,
378 : gsl::not_null<Cache*> cache,
379 : gr::Tags::Lapse<DataType> /*meta*/) const override;
380 : void operator()(gsl::not_null<tnsr::i<DataType, 3>*> deriv_lapse,
381 : gsl::not_null<Cache*> cache,
382 : ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>,
383 : Frame::Inertial> /*meta*/) const;
384 : void operator()(
385 : gsl::not_null<Scalar<DataType>*> lapse_times_conformal_factor_minus_one,
386 : gsl::not_null<Cache*> cache,
387 : Xcts::Tags::LapseTimesConformalFactorMinusOne<DataType> /*meta*/)
388 : const override;
389 : void operator()(
390 : gsl::not_null<Scalar<DataType>*> lapse_times_conformal_factor,
391 : gsl::not_null<Cache*> cache,
392 : Xcts::Tags::LapseTimesConformalFactor<DataType> /*meta*/) const override;
393 : void operator()(
394 : gsl::not_null<tnsr::i<DataType, 3>*>
395 : lapse_times_conformal_factor_gradient,
396 : gsl::not_null<Cache*> cache,
397 : ::Tags::deriv<Xcts::Tags::LapseTimesConformalFactorMinusOne<DataType>,
398 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
399 : void operator()(
400 : gsl::not_null<tnsr::I<DataType, 3>*> shift_background,
401 : gsl::not_null<Cache*> cache,
402 : Xcts::Tags::ShiftBackground<DataType, 3, Frame::Inertial> /*meta*/)
403 : const override;
404 : void operator()(
405 : gsl::not_null<tnsr::iJ<DataType, 3, Frame::Inertial>*>
406 : deriv_shift_background,
407 : gsl::not_null<Cache*> cache,
408 : ::Tags::deriv<Xcts::Tags::ShiftBackground<DataType, 3, Frame::Inertial>,
409 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
410 : void operator()(gsl::not_null<tnsr::II<DataType, 3, Frame::Inertial>*>
411 : longitudinal_shift_background_minus_dt_conformal_metric,
412 : gsl::not_null<Cache*> cache,
413 : Xcts::Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
414 : DataType, 3, Frame::Inertial> /*meta*/) const override;
415 : void operator()(
416 : gsl::not_null<tnsr::I<DataType, 3>*> shift_excess,
417 : gsl::not_null<Cache*> cache,
418 : Xcts::Tags::ShiftExcess<DataType, 3, Frame::Inertial> /*meta*/)
419 : const override;
420 : void operator()(
421 : gsl::not_null<tnsr::iJ<DataType, 3>*> deriv_shift_excess,
422 : gsl::not_null<Cache*> cache,
423 : ::Tags::deriv<Xcts::Tags::ShiftExcess<DataType, 3, Frame::Inertial>,
424 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
425 : void operator()(
426 : gsl::not_null<tnsr::ii<DataType, 3>*> extrinsic_curvature,
427 : gsl::not_null<Cache*> cache,
428 : gr::Tags::ExtrinsicCurvature<DataType, 3> /*meta*/) const override;
429 : void operator()(gsl::not_null<Scalar<DataType>*> energy_density,
430 : gsl::not_null<Cache*> cache,
431 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>,
432 : ConformalMatterScale> /*meta*/) const;
433 : void operator()(gsl::not_null<Scalar<DataType>*> stress_trace,
434 : gsl::not_null<Cache*> cache,
435 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>,
436 : ConformalMatterScale> /*meta*/) const;
437 : void operator()(gsl::not_null<tnsr::I<DataType, 3>*> momentum_density,
438 : gsl::not_null<Cache*> cache,
439 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, 3>,
440 : ConformalMatterScale> /*meta*/) const;
441 : };
442 :
443 : } // namespace detail
444 :
445 : /*!
446 : * \brief Schwarzschild spacetime in general relativity
447 : *
448 : * This class implements the Schwarzschild solution with mass parameter
449 : * \f$M\f$ in various coordinate systems. See the entries of the
450 : * `Xcts::Solutions::SchwarzschildCoordinates` enum for the available coordinate
451 : * systems and for the solution variables in the respective coordinates.
452 : */
453 1 : class Schwarzschild : public elliptic::analytic_data::AnalyticSolution,
454 : public detail::SchwarzschildImpl {
455 : public:
456 0 : Schwarzschild() = default;
457 0 : Schwarzschild(const Schwarzschild&) = default;
458 0 : Schwarzschild& operator=(const Schwarzschild&) = default;
459 0 : Schwarzschild(Schwarzschild&&) = default;
460 0 : Schwarzschild& operator=(Schwarzschild&&) = default;
461 0 : ~Schwarzschild() = default;
462 :
463 : using SchwarzschildImpl::SchwarzschildImpl;
464 :
465 : /// \cond
466 : explicit Schwarzschild(CkMigrateMessage* m)
467 : : elliptic::analytic_data::AnalyticSolution(m) {}
468 : using PUP::able::register_constructor;
469 : WRAPPED_PUPable_decl_template(Schwarzschild);
470 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
471 : const override {
472 : return std::make_unique<Schwarzschild>(*this);
473 : }
474 : /// \endcond
475 :
476 : template <typename DataType, typename... RequestedTags>
477 0 : tuples::TaggedTuple<RequestedTags...> variables(
478 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
479 : tmpl::list<RequestedTags...> /*meta*/) const {
480 : return variables_impl<DataType>(x, std::nullopt, std::nullopt,
481 : tmpl::list<RequestedTags...>{});
482 : }
483 :
484 : template <typename... RequestedTags>
485 0 : tuples::TaggedTuple<RequestedTags...> variables(
486 : const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& mesh,
487 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
488 : Frame::Inertial>& inv_jacobian,
489 : tmpl::list<RequestedTags...> /*meta*/) const {
490 : return variables_impl<DataVector>(x, mesh, inv_jacobian,
491 : tmpl::list<RequestedTags...>{});
492 : }
493 :
494 0 : void pup(PUP::er& p) override {
495 : elliptic::analytic_data::AnalyticSolution::pup(p);
496 : detail::SchwarzschildImpl::pup(p);
497 : }
498 :
499 : private:
500 : template <typename DataType, typename... RequestedTags>
501 0 : tuples::TaggedTuple<RequestedTags...> variables_impl(
502 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
503 : std::optional<std::reference_wrapper<const Mesh<3>>> mesh,
504 : std::optional<std::reference_wrapper<const InverseJacobian<
505 : DataType, 3, Frame::ElementLogical, Frame::Inertial>>>
506 : inv_jacobian,
507 : tmpl::list<RequestedTags...> /*meta*/) const {
508 : using VarsComputer = detail::SchwarzschildVariables<DataType>;
509 : typename VarsComputer::Cache cache{get_size(*x.begin())};
510 : const VarsComputer computer{std::move(mesh), std::move(inv_jacobian), x,
511 : mass_, coordinate_system_};
512 : const auto get_var = [&cache, &computer, &x](auto tag_v) {
513 : using tag = std::decay_t<decltype(tag_v)>;
514 : if constexpr (tmpl::list_contains_v<hydro_tags<DataType>, tag>) {
515 : (void)cache;
516 : (void)computer;
517 : return get<tag>(Flatness{}.variables(x, tmpl::list<tag>{}));
518 : } else {
519 : (void)x;
520 : return cache.get_var(computer, tag{});
521 : }
522 : };
523 : return {get_var(RequestedTags{})...};
524 : }
525 : };
526 :
527 : } // namespace Xcts::Solutions
|
