Schwarzschild.hpp
1 // 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/CachedTempBuffer.hpp"
11 #include "DataStructures/DataVector.hpp"
13 #include "Elliptic/Systems/Xcts/Tags.hpp"
15 #include "Options/Options.hpp"
16 #include "Options/ParseOptions.hpp"
18 #include "PointwiseFunctions/AnalyticSolutions/Xcts/AnalyticSolution.hpp"
19 #include "PointwiseFunctions/AnalyticSolutions/Xcts/CommonVariables.hpp"
20 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
21 #include "Utilities/ContainerHelpers.hpp"
22 #include "Utilities/Gsl.hpp"
23 #include "Utilities/TMPL.hpp"
24 #include "Utilities/TaggedTuple.hpp"
25 
26 /// \cond
27 namespace PUP {
28 class er;
29 } // namespace PUP
30 /// \endcond
31 
32 namespace Xcts::Solutions {
33 
34 /// Various coordinate systems in which to express the Schwarzschild solution
36  /*!
37  * \brief Isotropic Schwarzschild coordinates
38  *
39  * These arise from the canonical Schwarzschild coordinates by the radial
40  * transformation
41  *
42  * \f{equation}
43  * r = \bar{r}\left(1+\frac{M}{2\bar{r}}\right)^2
44  * \f}
45  *
46  * (Eq. (1.61) in \cite BaumgarteShapiro) where \f$r\f$ is the canonical
47  * Schwarzschild radius, also referred to as "areal" radius because it is
48  * defined such that spheres with constant \f$r\f$ have the area \f$4\pi
49  * r^2\f$, and \f$\bar{r}\f$ is the "isotropic" radius. In the isotropic
50  * radius the Schwarzschild spatial metric is conformally flat:
51  *
52  * \f{equation}
53  * \gamma_{ij}=\psi^4\eta_{ij} \quad \text{with conformal factor} \quad
54  * \psi=1+\frac{M}{2\bar{r}}
55  * \f}
56  *
57  * (Table 2.1 in \cite BaumgarteShapiro). Its lapse transforms to
58  *
59  * \f{equation}
60  * \alpha=\frac{1-M/(2\bar{r})}{1+M/(2\bar{r})}
61  * \f}
62  *
63  * and the shift vanishes (\f$\beta^i=0\f$) as it does in areal Schwarzschild
64  * coordinates. The solution also remains maximally sliced, i.e. \f$K=0\f$.
65  *
66  * The Schwarzschild horizon in these coordinates is at
67  * \f$\bar{r}=\frac{M}{2}\f$ due to the radial transformation from \f$r=2M\f$.
68  */
69  Isotropic,
70 };
71 
72 std::ostream& operator<<(std::ostream& os,
73  SchwarzschildCoordinates coords) noexcept;
74 
75 } // namespace Xcts::Solutions
76 
77 template <>
78 struct Options::create_from_yaml<Xcts::Solutions::SchwarzschildCoordinates> {
79  template <typename Metavariables>
81  const Options::Option& options) {
82  return create<void>(options);
83  }
84 };
85 
86 template <>
89  void>(const Options::Option& options);
90 
91 namespace Xcts::Solutions {
92 
93 namespace detail {
94 
95 struct SchwarzschildImpl {
96  struct Mass {
97  using type = double;
98  static constexpr Options::String help = "Mass parameter M";
99  };
100 
101  struct CoordinateSystem {
102  static std::string name() noexcept { return "Coordinates"; }
103  using type = SchwarzschildCoordinates;
104  static constexpr Options::String help =
105  "The coordinate system used to describe the solution";
106  };
107 
108  using options = tmpl::list<Mass, CoordinateSystem>;
109  static constexpr Options::String help{
110  "Schwarzschild spacetime in general relativity"};
111 
112  SchwarzschildImpl() = default;
113  SchwarzschildImpl(const SchwarzschildImpl&) noexcept = default;
114  SchwarzschildImpl& operator=(const SchwarzschildImpl&) noexcept = default;
115  SchwarzschildImpl(SchwarzschildImpl&&) noexcept = default;
116  SchwarzschildImpl& operator=(SchwarzschildImpl&&) noexcept = default;
117  ~SchwarzschildImpl() noexcept = default;
118 
119  explicit SchwarzschildImpl(
120  double mass, SchwarzschildCoordinates coordinate_system) noexcept;
121 
122  /// The mass parameter \f$M\f$.
123  double mass() const noexcept;
124 
125  SchwarzschildCoordinates coordinate_system() const noexcept;
126 
127  /// The radius of the Schwarzschild horizon in the given coordinates.
128  double radius_at_horizon() const noexcept;
129 
130  // NOLINTNEXTLINE(google-runtime-references)
131  void pup(PUP::er& p) noexcept;
132 
133  protected:
135  SchwarzschildCoordinates coordinate_system_{};
136 };
137 
138 bool operator==(const SchwarzschildImpl& lhs,
139  const SchwarzschildImpl& rhs) noexcept;
140 
141 bool operator!=(const SchwarzschildImpl& lhs,
142  const SchwarzschildImpl& rhs) noexcept;
143 
144 template <typename DataType>
145 struct SchwarzschildVariables;
146 
147 template <typename DataType>
148 using SchwarzschildVariablesCache = cached_temp_buffer_from_typelist<
149  SchwarzschildVariables<DataType>,
151  common_tags<DataType>,
152  Tags::ConformalMetric<DataType, 3, Frame::Inertial>,
154  tmpl::size_t<3>, Frame::Inertial>,
157  Tags::ConformalFactor<DataType>,
160  gr::Tags::Lapse<DataType>, Tags::LapseTimesConformalFactor<DataType>,
162  tmpl::size_t<3>, Frame::Inertial>,
163  Tags::ShiftBackground<DataType, 3, Frame::Inertial>,
164  Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
165  DataType, 3, Frame::Inertial>,
166  Tags::ShiftExcess<DataType, 3, Frame::Inertial>,
167  Tags::ShiftStrain<DataType, 3, Frame::Inertial>,
168  Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
169  Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
170  Tags::Conformal<gr::Tags::MomentumDensity<3, Frame::Inertial, DataType>,
171  0>>>;
172 
173 template <typename DataType>
174 struct SchwarzschildVariables
175  : CommonVariables<DataType, SchwarzschildVariablesCache<DataType>> {
176  static constexpr size_t Dim = 3;
177  static constexpr int ConformalMatterScale = 0;
178  using Cache = SchwarzschildVariablesCache<DataType>;
179  using CommonVariables<DataType,
180  SchwarzschildVariablesCache<DataType>>::operator();
181 
182  const tnsr::I<DataType, 3>& x;
183  double mass;
184  SchwarzschildCoordinates coordinate_system;
185 
186  void operator()(gsl::not_null<tnsr::ii<DataType, 3>*> conformal_metric,
188  Tags::ConformalMetric<DataType, 3, Frame::Inertial> /*meta*/)
189  const noexcept;
190  void operator()(
191  gsl::not_null<tnsr::II<DataType, 3>*> inv_conformal_metric,
193  Tags::InverseConformalMetric<DataType, 3, Frame::Inertial> /*meta*/)
194  const noexcept;
195  void operator()(
196  gsl::not_null<tnsr::ijj<DataType, 3>*> deriv_conformal_metric,
198  ::Tags::deriv<Tags::ConformalMetric<DataType, 3, Frame::Inertial>,
199  tmpl::size_t<3>, Frame::Inertial> /*meta*/) const noexcept;
200  void operator()(
201  gsl::not_null<Scalar<DataType>*> trace_extrinsic_curvature,
202  gsl::not_null<Cache*> cache,
203  gr::Tags::TraceExtrinsicCurvature<DataType> /*meta*/) const noexcept;
204  void operator()(
205  gsl::not_null<tnsr::i<DataType, 3>*> trace_extrinsic_curvature_gradient,
206  gsl::not_null<Cache*> cache,
207  ::Tags::deriv<gr::Tags::TraceExtrinsicCurvature<DataType>,
208  tmpl::size_t<3>, Frame::Inertial> /*meta*/) const noexcept;
209  void operator()(
210  gsl::not_null<Scalar<DataType>*> dt_trace_extrinsic_curvature,
211  gsl::not_null<Cache*> cache,
212  ::Tags::dt<gr::Tags::TraceExtrinsicCurvature<DataType>> /*meta*/)
213  const noexcept;
214  void operator()(gsl::not_null<Scalar<DataType>*> conformal_factor,
215  gsl::not_null<Cache*> cache,
216  Tags::ConformalFactor<DataType> /*meta*/) const noexcept;
217  void operator()(
218  gsl::not_null<tnsr::i<DataType, 3>*> conformal_factor_gradient,
219  gsl::not_null<Cache*> cache,
220  ::Tags::deriv<Xcts::Tags::ConformalFactor<DataType>, tmpl::size_t<3>,
221  Frame::Inertial> /*meta*/) const noexcept;
222  void operator()(gsl::not_null<Scalar<DataType>*> lapse,
223  gsl::not_null<Cache*> cache,
224  gr::Tags::Lapse<DataType> /*meta*/) const noexcept;
225  void operator()(
226  gsl::not_null<Scalar<DataType>*> lapse_times_conformal_factor,
227  gsl::not_null<Cache*> cache,
228  Tags::LapseTimesConformalFactor<DataType> /*meta*/) const noexcept;
229  void operator()(
230  gsl::not_null<tnsr::i<DataType, 3>*>
231  lapse_times_conformal_factor_gradient,
232  gsl::not_null<Cache*> cache,
233  ::Tags::deriv<Tags::LapseTimesConformalFactor<DataType>, tmpl::size_t<3>,
234  Frame::Inertial> /*meta*/) const noexcept;
235  void operator()(gsl::not_null<tnsr::I<DataType, 3>*> shift_background,
236  gsl::not_null<Cache*> cache,
237  Tags::ShiftBackground<DataType, 3, Frame::Inertial> /*meta*/)
238  const noexcept;
239  void operator()(gsl::not_null<tnsr::II<DataType, 3, Frame::Inertial>*>
240  longitudinal_shift_background_minus_dt_conformal_metric,
241  gsl::not_null<Cache*> cache,
242  Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
243  DataType, 3, Frame::Inertial> /*meta*/) const noexcept;
244  void operator()(
245  gsl::not_null<tnsr::I<DataType, 3>*> shift_excess,
246  gsl::not_null<Cache*> cache,
247  Tags::ShiftExcess<DataType, 3, Frame::Inertial> /*meta*/) const noexcept;
248  void operator()(
249  gsl::not_null<tnsr::ii<DataType, 3>*> shift_strain,
250  gsl::not_null<Cache*> cache,
251  Tags::ShiftStrain<DataType, 3, Frame::Inertial> /*meta*/) const noexcept;
252  void operator()(
253  gsl::not_null<Scalar<DataType>*> energy_density,
254  gsl::not_null<Cache*> cache,
255  Tags::Conformal<gr::Tags::EnergyDensity<DataType>,
256  ConformalMatterScale> /*meta*/) const noexcept;
257  void operator()(
258  gsl::not_null<Scalar<DataType>*> stress_trace,
259  gsl::not_null<Cache*> cache,
260  Tags::Conformal<gr::Tags::StressTrace<DataType>,
261  ConformalMatterScale> /*meta*/) const noexcept;
262  void operator()(
263  gsl::not_null<tnsr::I<DataType, 3>*> momentum_density,
264  gsl::not_null<Cache*> cache,
265  Tags::Conformal<gr::Tags::MomentumDensity<3, Frame::Inertial, DataType>,
266  ConformalMatterScale> /*meta*/) const noexcept;
267 };
268 
269 } // namespace detail
270 
271 // The following implements the registration and factory-creation mechanism
272 
273 /// \cond
274 template <typename Registrars>
275 struct Schwarzschild;
276 
277 namespace Registrars {
278 struct Schwarzschild {
279  template <typename Registrars>
280  using f = Solutions::Schwarzschild<Registrars>;
281 };
282 } // namespace Registrars
283 /// \endcond
284 
285 /*!
286  * \brief Schwarzschild spacetime in general relativity
287  *
288  * This class implements the Schwarzschild solution with mass parameter
289  * \f$M\f$ in various coordinate systems. See the entries of the
290  * `Xcts::Solutions::SchwarzschildCoordinates` enum for the available coordinate
291  * systems and for the solution variables in the respective coordinates.
292  */
293 template <typename Registrars =
294  tmpl::list<Solutions::Registrars::Schwarzschild>>
295 class Schwarzschild : public AnalyticSolution<Registrars>,
296  public detail::SchwarzschildImpl {
297  private:
299 
300  public:
301  Schwarzschild() = default;
302  Schwarzschild(const Schwarzschild&) noexcept = default;
303  Schwarzschild& operator=(const Schwarzschild&) noexcept = default;
304  Schwarzschild(Schwarzschild&&) noexcept = default;
305  Schwarzschild& operator=(Schwarzschild&&) noexcept = default;
306  ~Schwarzschild() noexcept = default;
307 
308  using SchwarzschildImpl::SchwarzschildImpl;
309 
310  /// \cond
311  explicit Schwarzschild(CkMigrateMessage* m) noexcept : Base(m) {}
312  using PUP::able::register_constructor;
314  /// \endcond
315 
316  template <typename DataType, typename... RequestedTags>
317  tuples::TaggedTuple<RequestedTags...> variables(
318  const tnsr::I<DataType, 3, Frame::Inertial>& x,
319  tmpl::list<RequestedTags...> /*meta*/) const noexcept {
320  using VarsComputer = detail::SchwarzschildVariables<DataType>;
321  typename VarsComputer::Cache cache{
322  get_size(*x.begin()),
323  VarsComputer{
324  {{std::nullopt, std::nullopt}}, x, mass_, coordinate_system_}};
325  return {cache.get_var(RequestedTags{})...};
326  }
327 
328  template <typename... RequestedTags>
329  tuples::TaggedTuple<RequestedTags...> variables(
330  const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& mesh,
331  const InverseJacobian<DataVector, 3, Frame::Logical, Frame::Inertial>&
332  inv_jacobian,
333  tmpl::list<RequestedTags...> /*meta*/) const noexcept {
334  using VarsComputer = detail::SchwarzschildVariables<DataVector>;
335  typename VarsComputer::Cache cache{
336  get_size(*x.begin()),
337  VarsComputer{{{mesh, inv_jacobian}}, x, mass_, coordinate_system_}};
338  return {cache.get_var(RequestedTags{})...};
339  }
340 
341  void pup(PUP::er& p) noexcept override {
342  Base::pup(p);
343  detail::SchwarzschildImpl::pup(p);
344  }
345 };
346 
347 /// \cond
348 template <typename Registrars>
349 PUP::able::PUP_ID Schwarzschild<Registrars>::my_PUP_ID = 0; // NOLINT
350 /// \endcond
351 
352 } // namespace Xcts::Solutions
domain::push_front
CoordinateMap< SourceFrame, TargetFrame, NewMap, Maps... > push_front(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by prepending the new map to the beginning of the old maps.
gr::Tags::TraceExtrinsicCurvature
Definition: Tags.hpp:120
Xcts
Items related to solving the Extended Conformal Thin Sandwich (XCTS) decomposition of the Einstein co...
Definition: Flatness.hpp:22
std::string
CharmPupable.hpp
Frame::Inertial
Definition: IndexType.hpp:44
std::rel_ops::operator!=
T operator!=(T... args)
Options.hpp
ParseOptions.hpp
get_size
decltype(auto) get_size(const T &t, SizeFunction size=GetContainerSize{}) noexcept
Retrieve the size of t if t.size() is a valid expression, otherwise if T is fundamental or a std::com...
Definition: ContainerHelpers.hpp:145
Options::Option
Definition: Options.hpp:108
gr::lapse
Scalar< DataType > lapse(const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute lapse from shift and spacetime metric.
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
std::ostream
Options::create_from_yaml
Definition: MinmodType.hpp:11
Xcts::Solutions::SchwarzschildCoordinates
SchwarzschildCoordinates
Various coordinate systems in which to express the Schwarzschild solution.
Definition: Schwarzschild.hpp:35
Xcts::Solutions::AnalyticSolution
Base class for analytic solutions of the XCTS equations.
Definition: AnalyticSolution.hpp:34
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
Xcts::Solutions::SchwarzschildCoordinates::Isotropic
@ Isotropic
Isotropic Schwarzschild coordinates.
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
Tags::dt
Prefix indicating a time derivative.
Definition: Prefixes.hpp:29
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:382
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
limits
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Frame
Definition: IndexType.hpp:36
gr
Definition: GaugeWave.hpp:28
Tensor.hpp
Tags::deriv
Prefix indicating spatial derivatives.
Definition: PartialDerivatives.hpp:52
PartialDerivatives.hpp
Xcts::Solutions
Analytic solutions of the XCTS equations.
Definition: AnalyticSolution.hpp:18
ostream
Prefixes.hpp
gsl
Implementations from the Guideline Support Library.
Definition: ReadSpecPiecewisePolynomial.hpp:11
gr::Tags::Lapse
Definition: Tags.hpp:52
TMPL.hpp
Xcts::Solutions::Schwarzschild
Schwarzschild spacetime in general relativity.
Definition: Schwarzschild.hpp:295
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13