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 <pup.h>
9 :
10 : #include "DataStructures/DataBox/Prefixes.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/TaggedTuple.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Elliptic/Systems/Xcts/Tags.hpp"
15 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
16 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
17 : #include "Options/Context.hpp"
18 : #include "Options/String.hpp"
19 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
20 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
21 : #include "Utilities/Serialization/CharmPupable.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : /// \cond
25 : namespace PUP {
26 : class er;
27 : } // namespace PUP
28 : /// \endcond
29 :
30 : namespace Xcts::Solutions {
31 :
32 : /*!
33 : * \brief A constant density star in general relativity
34 : *
35 : * \details This solution describes a star with constant density \f$\rho_0\f$
36 : * that extends to a (conformal) radius \f$R\f$. It solves the XCTS Hamiltonian
37 : * constraint that reduces to the non-linear elliptic equation
38 : * \f[
39 : * \Delta^2\psi+2\pi\rho\psi^5=0
40 : * \f]
41 : * for the conformal factor \f$\psi\f$ (see `Xcts`) under the following
42 : * assumptions \cite Baumgarte2006ug :
43 : *
44 : * - Time-symmetry \f$K_{ij}=0\f$
45 : * - Conformal flatness \f$\overline{\gamma}=\delta\f$, so \f$\Delta\f$ is the
46 : * flat-space Laplacian
47 : * - Spherical symmetry
48 : *
49 : * Imposing boundary conditions
50 : * \f[
51 : * \frac{\partial\psi}{\partial r}=0 \quad \text{for} \quad r=0\\
52 : * \psi\rightarrow 1 \quad \text{for} \quad r\rightarrow\infty
53 : * \f]
54 : * and considering the energy density
55 : * \f[
56 : * \rho(r\leq R)=\rho_0 \quad \text{and} \quad \rho(r>R)=0
57 : * \f]
58 : * of the star the authors of \cite Baumgarte2006ug find the solution
59 : * \f[
60 : * \psi(r\leq R)=C u_\alpha(r) \quad \text{and}
61 : * \quad \psi(r>R)=\frac{\beta}{r} + 1
62 : * \f]
63 : * with \f$C=(2\pi\rho_0/3)^{-1/4}\f$, the Sobolev functions
64 : * \f[
65 : * u_\alpha(r)=\sqrt{\frac{\alpha R}{r^2+(\alpha R)^2}}
66 : * \f]
67 : * and real parameters \f$\alpha\f$ and \f$\beta\f$ that are determined by
68 : * the following relations:
69 : * \f[
70 : * \rho_0 R^2=\frac{3}{2\pi}f^2(\alpha) \quad \text{with}
71 : * \quad f(\alpha)=\frac{\alpha^5}{(1+\alpha^2)^3} \\
72 : * \frac{\beta}{R} + 1 = C u_\alpha(R)
73 : * \f]
74 : *
75 : * This solution is described in detail in \cite Baumgarte2006ug , and also in
76 : * Exercise 3.8 in \cite BaumgarteShapiro , since it
77 : * exhibits the non-uniqueness properties that are typical for the XCTS system.
78 : * In the simple case of the constant-density star the non-uniqueness is
79 : * apparent from the function \f$f(\alpha)\f$, which has two solutions for any
80 : * \f$\rho_0\f$ smaller than a critical density
81 : * \f[
82 : * \rho_\mathrm{crit}=\frac{3}{2\pi R^2}\frac{5^2}{6^6}
83 : * \approx\frac{0.0320}{R^2} \text{,}
84 : * \f]
85 : * a unique solution for \f$\rho_0=\rho_\mathrm{crit}\f$ and no solutions
86 : * above the critical density \cite Baumgarte2006ug . The authors identify the
87 : * \f$\alpha < \alpha_\mathrm{crit}=\sqrt{5}\f$ and \f$\alpha >
88 : * \alpha_\mathrm{crit}\f$ branches of solutions with the strong-field and
89 : * weak-field regimes, respectively (see \cite Baumgarte2006ug for details).
90 : * In this implementation we compute the weak-field solution by choosing the
91 : * \f$\alpha > \alpha_\mathrm{crit}\f$ that corresponds to the density
92 : * \f$\rho_0\f$ of the star. Therefore, we supply initial data
93 : * \f$\psi_\mathrm{init}=1\f$ so that a nonlinear iterative numerical solver
94 : * will converge to the same weak-field solution.
95 : */
96 1 : class ConstantDensityStar : public elliptic::analytic_data::AnalyticSolution {
97 : public:
98 0 : struct Density {
99 0 : using type = double;
100 0 : static constexpr Options::String help{
101 : "The constant density within the star"};
102 0 : static double lower_bound() { return 0.; }
103 : };
104 0 : struct Radius {
105 0 : using type = double;
106 0 : static constexpr Options::String help{"The conformal radius of the star"};
107 0 : static double lower_bound() { return 0.; }
108 : };
109 :
110 0 : using options = tmpl::list<Density, Radius>;
111 0 : static constexpr Options::String help{
112 : "A constant density star in general relativity"};
113 :
114 0 : ConstantDensityStar() = default;
115 0 : ConstantDensityStar(const ConstantDensityStar&) = default;
116 0 : ConstantDensityStar& operator=(const ConstantDensityStar&) = default;
117 0 : ConstantDensityStar(ConstantDensityStar&&) = default;
118 0 : ConstantDensityStar& operator=(ConstantDensityStar&&) = default;
119 0 : ~ConstantDensityStar() = default;
120 0 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
121 : const override {
122 : return std::make_unique<ConstantDensityStar>(*this);
123 : }
124 :
125 : /// \cond
126 : explicit ConstantDensityStar(CkMigrateMessage* m)
127 : : elliptic::analytic_data::AnalyticSolution(m) {}
128 : using PUP::able::register_constructor;
129 : WRAPPED_PUPable_decl_template(ConstantDensityStar);
130 : /// \endcond
131 :
132 0 : ConstantDensityStar(double density, double radius,
133 : const Options::Context& context = {});
134 :
135 0 : double density() const { return density_; }
136 0 : double radius() const { return radius_; }
137 :
138 : /// @{
139 : /// Retrieve variable at coordinates `x`
140 : template <typename DataType>
141 1 : auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
142 : tmpl::list<Xcts::Tags::ConformalFactor<DataType>> /*meta*/)
143 : const -> tuples::TaggedTuple<Xcts::Tags::ConformalFactor<DataType>>;
144 :
145 : template <typename DataType>
146 1 : auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
147 : tmpl::list<::Tags::Initial<
148 : Xcts::Tags::ConformalFactor<DataType>>> /*meta*/) const
149 : -> tuples::TaggedTuple<
150 : ::Tags::Initial<Xcts::Tags::ConformalFactor<DataType>>>;
151 :
152 : template <typename DataType>
153 1 : auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
154 : tmpl::list<::Tags::Initial<
155 : ::Tags::deriv<Xcts::Tags::ConformalFactor<DataType>,
156 : tmpl::size_t<3>, Frame::Inertial>>> /*meta*/)
157 : const -> tuples::TaggedTuple<
158 : ::Tags::Initial<::Tags::deriv<Xcts::Tags::ConformalFactor<DataType>,
159 : tmpl::size_t<3>, Frame::Inertial>>>;
160 :
161 : template <typename DataType>
162 1 : auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
163 : tmpl::list<::Tags::FixedSource<
164 : Xcts::Tags::ConformalFactor<DataType>>> /*meta*/) const
165 : -> tuples::TaggedTuple<
166 : ::Tags::FixedSource<Xcts::Tags::ConformalFactor<DataType>>>;
167 :
168 : template <typename DataType>
169 1 : auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
170 : tmpl::list<gr::Tags::EnergyDensity<DataType>> /*meta*/) const
171 : -> tuples::TaggedTuple<gr::Tags::EnergyDensity<DataType>>;
172 : /// @}
173 :
174 : /// Retrieve a collection of variables at coordinates `x`
175 : template <typename DataType, typename... Tags>
176 1 : tuples::TaggedTuple<Tags...> variables(
177 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
178 : tmpl::list<Tags...> /*meta*/) const {
179 : static_assert(sizeof...(Tags) > 1,
180 : "The generic template will recurse infinitely if only one "
181 : "tag is being retrieved.");
182 : return {tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
183 : }
184 :
185 : template <typename... Tags>
186 0 : tuples::TaggedTuple<Tags...> variables(
187 : const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& /*mesh*/,
188 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
189 : Frame::Inertial>& /*inv_jacobian*/,
190 : tmpl::list<Tags...> meta) const {
191 : return variables(x, meta);
192 : }
193 :
194 : // NOLINTNEXTLINE(google-runtime-references)
195 0 : void pup(PUP::er& p) override;
196 :
197 : private:
198 0 : double density_ = std::numeric_limits<double>::signaling_NaN();
199 0 : double radius_ = std::numeric_limits<double>::signaling_NaN();
200 0 : double alpha_ = std::numeric_limits<double>::signaling_NaN();
201 : };
202 :
203 0 : bool operator==(const ConstantDensityStar& /*lhs*/,
204 : const ConstantDensityStar& /*rhs*/);
205 :
206 0 : bool operator!=(const ConstantDensityStar& lhs, const ConstantDensityStar& rhs);
207 :
208 : } // namespace Xcts::Solutions
|
