ConstantDensityStar.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <limits>
7 
8 #include "DataStructures/DataBox/Prefixes.hpp" // IWYU pragma: keep
9 #include "DataStructures/Tensor/Tensor.hpp" // IWYU pragma: keep
10 #include "Elliptic/Systems/Xcts/Tags.hpp" // IWYU pragma: keep
11 #include "Options/Options.hpp"
12 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
13 #include "Utilities/TMPL.hpp"
15 
16 // IWYU pragma: no_forward_declare Tensor
17 
18 /// \cond
19 namespace PUP {
20 class er;
21 } // namespace PUP
22 /// \endcond
23 
24 namespace Xcts {
25 namespace Solutions {
26 
27 /*!
28  * \brief A constant density star in general relativity
29  *
30  * \details This solution describes a star with constant density \f$\rho_0\f$
31  * that extends to a (conformal) radius \f$R\f$. It solves the XCTS Hamiltonian
32  * constraint that reduces to the non-linear elliptic equation
33  * \f[
34  * \Delta^2\psi+2\pi\rho\psi^5=0
35  * \f]
36  * for the conformal factor \f$\psi\f$ (see `Xcts`) under the following
37  * assumptions \cite Baumgarte2006ug :
38  *
39  * - Time-symmetry \f$K_{ij}=0\f$
40  * - Conformal flatness \f$\overline{\gamma}=\delta\f$, so \f$\Delta\f$ is the
41  * flat-space Laplacian
42  * - Spherical symmetry
43  *
44  * Imposing boundary conditions
45  * \f[
46  * \frac{\partial\psi}{\partial r}=0 \quad \text{for} \quad r=0\\
47  * \psi\rightarrow 1 \quad \text{for} \quad r\rightarrow\infty
48  * \f]
49  * and considering the energy density
50  * \f[
51  * \rho(r\leq R)=\rho_0 \quad \text{and} \quad \rho(r>R)=0
52  * \f]
53  * of the star the authors of \cite Baumgarte2006ug find the solution
54  * \f[
55  * \psi(r\leq R)=C u_\alpha(r) \quad \text{and}
56  * \quad \psi(r>R)=\frac{\beta}{r} + 1
57  * \f]
58  * with \f$C=(2\pi\rho_0/3)^{-1/4}\f$, the Sobolev functions
59  * \f[
60  * u_\alpha(r)=\sqrt{\frac{\alpha R}{r^2+(\alpha R)^2}}
61  * \f]
62  * and real parameters \f$\alpha\f$ and \f$\beta\f$ that are determined by
63  * the following relations:
64  * \f[
65  * \rho_0 R^2=\frac{3}{2\pi}f^2(\alpha) \quad \text{with}
66  * \quad f(\alpha)=\frac{\alpha^5}{(1+\alpha^2)^3} \\
67  * \frac{\beta}{R} + 1 = C u_\alpha(R)
68  * \f]
69  *
70  * This solution is described in detail in \cite Baumgarte2006ug since it
71  * exhibits the non-uniqueness properties that are typical for the XCTS system.
72  * In the simple case of the constant-density star the non-uniqueness is
73  * apparent from the function \f$f(\alpha)\f$, which has two solutions for any
74  * \f$\rho_0\f$ smaller than a critical density
75  * \f[
76  * \rho_\mathrm{crit}=\frac{3}{2\pi R^2}\frac{5^2}{6^6}
77  * \approx\frac{0.0320}{R^2} \text{,}
78  * \f]
79  * a unique solution for \f$\rho_0=\rho_\mathrm{crit}\f$ and no solutions
80  * above the critical density \cite Baumgarte2006ug . The authors identify the
81  * \f$\alpha < \alpha_\mathrm{crit}=\sqrt{5}\f$ and \f$\alpha >
82  * \alpha_\mathrm{crit}\f$ branches of solutions with the strong-field and
83  * weak-field regimes, respectively (see \cite Baumgarte2006ug for details).
84  * In this implementation we compute the weak-field solution by choosing the
85  * \f$\alpha > \alpha_\mathrm{crit}\f$ that corresponds to the density
86  * \f$\rho_0\f$ of the star. Therefore, we supply initial data
87  * \f$\psi_\mathrm{init}=1\f$ so that a nonlinear iterative numerical solver
88  * will converge to the same weak-field solution.
89  */
91  private:
92  struct Density {
93  using type = double;
94  static constexpr OptionString help{"The constant density within the star"};
95  static double lower_bound() noexcept { return 0.; }
96  };
97  struct Radius {
98  using type = double;
99  static constexpr OptionString help{"The conformal radius of the star"};
100  static double lower_bound() noexcept { return 0.; }
101  };
102 
103  public:
104  using options = tmpl::list<Density, Radius>;
105  static constexpr OptionString help{
106  "A constant density star in general relativity"};
107 
108  ConstantDensityStar() = default;
109  ConstantDensityStar(const ConstantDensityStar&) noexcept = delete;
110  ConstantDensityStar& operator=(const ConstantDensityStar&) noexcept = delete;
111  ConstantDensityStar(ConstantDensityStar&&) noexcept = default;
112  ConstantDensityStar& operator=(ConstantDensityStar&&) noexcept = default;
113  ~ConstantDensityStar() noexcept = default;
114 
115  ConstantDensityStar(double density, double radius,
116  const OptionContext& context = {});
117 
118  double density() const noexcept { return density_; }
119  double radius() const noexcept { return radius_; }
120 
121  // @{
122  /// Retrieve variable at coordinates `x`
123  template <typename DataType>
124  auto variables(
125  const tnsr::I<DataType, 3, Frame::Inertial>& x,
126  tmpl::list<Xcts::Tags::ConformalFactor<DataType>> /*meta*/) const noexcept
128 
129  template <typename DataType>
130  auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
131  tmpl::list<::Tags::Initial<
132  Xcts::Tags::ConformalFactor<DataType>>> /*meta*/) const
133  noexcept -> tuples::TaggedTuple<
135 
136  template <typename DataType>
137  auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
139  3, Frame::Inertial, DataType>>> /*meta*/) const noexcept
142 
143  template <typename DataType>
144  auto variables(
145  const tnsr::I<DataType, 3, Frame::Inertial>& x,
146  tmpl::list<
148  noexcept -> tuples::TaggedTuple<
150 
151  template <typename DataType>
152  auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
154  3, Frame::Inertial, DataType>>> /*meta*/) const noexcept
156  Xcts::Tags::ConformalFactorGradient<3, Frame::Inertial, DataType>>>;
157 
158  template <typename DataType>
159  auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
160  tmpl::list<gr::Tags::EnergyDensity<DataType>> /*meta*/) const
162  // @}
163 
164  /// Retrieve a collection of variables at coordinates `x`
165  template <typename DataType, typename... Tags>
167  const tnsr::I<DataType, 3, Frame::Inertial>& x,
168  tmpl::list<Tags...> /*meta*/) const noexcept {
169  static_assert(sizeof...(Tags) > 1,
170  "The generic template will recurse infinitely if only one "
171  "tag is being retrieved.");
172  return {tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
173  }
174 
175  // clang-tidy: no pass by reference
176  void pup(PUP::er& p) noexcept; // NOLINT
177 
178  private:
179  double density_ = std::numeric_limits<double>::signaling_NaN();
182 };
183 
184 bool operator==(const ConstantDensityStar& /*lhs*/,
185  const ConstantDensityStar& /*rhs*/) noexcept;
186 
187 bool operator!=(const ConstantDensityStar& lhs,
188  const ConstantDensityStar& rhs) noexcept;
189 
190 } // namespace Solutions
191 } // namespace Xcts
Prefix indicating a source term.
Definition: Prefixes.hpp:66
Definition: Strahlkorper.hpp:14
Defines class tuples::TaggedTuple.
A constant density star in general relativity.
Definition: ConstantDensityStar.hpp:90
T signaling_NaN(T... args)
Defines classes and functions for making classes creatable from input files.
Define prefixes for DataBox tags.
The gradient of the conformal factor .
Definition: Tags.hpp:36
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:27
tuples::TaggedTuple< Tags... > variables(const tnsr::I< DataType, 3, Frame::Inertial > &x, tmpl::list< Tags... >) const noexcept
Retrieve a collection of variables at coordinates x
Definition: ConstantDensityStar.hpp:166
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Prefix indicating the initial value of a quantity.
Definition: Prefixes.hpp:75
Definition: DataBoxTag.hpp:29
Defines classes for Tensor.
The energy density , where denotes the normal to the spatial hypersurface.
Definition: Tags.hpp:111
Information about the nested operations being performed by the parser, for use in printing errors...
Definition: Options.hpp:36
Wraps the template metaprogramming library used (brigand)
Items related to solving the Extended Conformal Thin Sandwich (XCTS) equations.
Definition: Tags.hpp:16
Definition: IndexType.hpp:44
The conformal factor that rescales the spatial metric .
Definition: Tags.hpp:24