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 <ostream>
8 :
9 : #include "NumericalAlgorithms/Interpolation/CubicSpline.hpp"
10 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
11 :
12 : /// \cond
13 : namespace Options {
14 : struct Option;
15 : template <typename T>
16 : struct create_from_yaml;
17 : } // namespace Options
18 : namespace PUP {
19 : class er;
20 : } // namespace PUP
21 : /// \endcond
22 :
23 : namespace RelativisticEuler::Solutions {
24 :
25 : /// Radial coordinate of a TOV solution.
26 1 : enum class TovCoordinates {
27 : /// Areal radial coordinate, in which the exterior TOV solution coincides with
28 : /// the Schwarzschild metric in standard Schwarzschild coordinates.
29 : Schwarzschild,
30 : /// Isotropic radial coordinate, in which the spatial metric is conformally
31 : /// flat.
32 : Isotropic
33 : };
34 :
35 0 : std::ostream& operator<<(std::ostream& os, TovCoordinates coords);
36 :
37 : } // namespace RelativisticEuler::Solutions
38 :
39 : template <>
40 0 : struct Options::create_from_yaml<RelativisticEuler::Solutions::TovCoordinates> {
41 : template <typename Metavariables>
42 0 : static RelativisticEuler::Solutions::TovCoordinates create(
43 : const Options::Option& options) {
44 : return create<void>(options);
45 : }
46 : };
47 :
48 : template <>
49 0 : RelativisticEuler::Solutions::TovCoordinates
50 : Options::create_from_yaml<RelativisticEuler::Solutions::TovCoordinates>::create<
51 : void>(const Options::Option& options);
52 :
53 : namespace RelativisticEuler::Solutions {
54 :
55 : /*!
56 : * \brief TOV solver based on Lindblom's method
57 : *
58 : * Uses Lindblom's method of integrating the TOV equations from
59 : * \cite Lindblom1998dp .
60 : *
61 : * Instead of integrating the interior mass \f$m(r)\f$ and pressure \f$p(r)\f$,
62 : * Lindblom introduces the variables \f$u=r^2\f$ and \f$v=m/r\f$. Then, the TOV
63 : * equations are integrated with the log of the specific enthalpy as the
64 : * independent variable, \f$\ln(h)\f$, from the center of the star where
65 : * \f$h(r=0) = h_c\f$ to its surface \f$h(r=R) = 1\f$. The ODEs being solved are
66 : * Eq. (A2) and (A3) in \cite Lindblom1998dp :
67 : *
68 : * \f{align}
69 : * \frac{\mathrm{d}u}{\mathrm{d}\ln{h}} &= \frac{-2u (1 - 2v)}{4\pi u p + v} \\
70 : * \frac{\mathrm{d}v}{\mathrm{d}\ln{h}} &=
71 : * -(1 - 2v) \frac{4\pi u \rho - v}{4\pi u p + v}
72 : * \f}
73 : *
74 : * Note that Lindblom's paper labels the independent variable as \f$h\f$.
75 : * However the \f$h\f$ in Lindblom's paper is **not** the specific enthalpy but
76 : * its logarithm, \f$\ln(h)\f$.
77 : *
78 : * The ODEs are solved numerically when this class is constructed, and the
79 : * quantities \f$m(r)/r\f$ and \f$\ln(h)\f$ are interpolated and exposed as
80 : * member functions. With these quantities the metric can be constructed as:
81 : *
82 : * \f{equation}
83 : * \mathrm{d}s^2 = -\alpha^2 \mathrm{d}t^2 + (1 - 2m/r)^{-1} \mathrm{d}r^2
84 : * + r^2 \mathrm{d}\Omega^2
85 : * \f}
86 : *
87 : * where the lapse is
88 : *
89 : * \f{equation}
90 : * \alpha(r < R) = \mathcal{E} / h = \alpha(r=R) / h
91 : * \text{,}
92 : * \f}
93 : *
94 : * with the conserved `injection_energy()` \f$\mathcal{E}\f$, such that the
95 : * lapse matches the exterior Schwarzschild solution:
96 : *
97 : * \f{equation}
98 : * \alpha(r \geq R) = \sqrt{1 - \frac{2M}{r}}
99 : * \f}
100 : *
101 : * \par Isotropic radial coordinate
102 : * This class also supports transforming to an isotropic radial coordinate. When
103 : * you pass `RelativisticEuler::Solutions::TovCoordinates::Isotropic` to the
104 : * constructor, an additional ODE is integrated alongside the TOV equations to
105 : * determine the conformal factor
106 : *
107 : * \f{equation}
108 : * \psi^2 = \frac{r}{\bar{r}}
109 : * \f}
110 : *
111 : * where \f$r\f$ is the areal (Schwarzschild) radius and \f$\bar{r}\f$ is the
112 : * isotropic radius. The additional ODE is:
113 : *
114 : * \f{equation}
115 : * \frac{\mathrm{d}\ln(\psi)}{\mathrm{d}\ln{h}} =
116 : * \frac{\sqrt{1 - 2v}}{1 + \sqrt{1 - 2v}} \frac{v}{4\pi u p + v}
117 : * \f}
118 : *
119 : * In isotropic coordinates, the spatial metric is conformally flat:
120 : *
121 : * \f{equation}
122 : * \mathrm{d}s^2 = -\alpha^2 \mathrm{d}t^2 + \psi^4 (\mathrm{d}\bar{r}^2 +
123 : * \bar{r}^2 \mathrm{d}\Omega^2)
124 : * \f}
125 : *
126 : * When isotropic coordinates are selected, radii returned by member functions
127 : * or taken as arguments are isotropic. An exception is `mass_over_radius()`,
128 : * which always returns the quantity \f$m / r\f$ because that is the quantity
129 : * stored internally and hence most numerically precise. See
130 : * `mass_over_radius()` for details.
131 : */
132 1 : class TovSolution {
133 : public:
134 0 : TovSolution(
135 : const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
136 : double central_mass_density,
137 : const TovCoordinates coordinate_system = TovCoordinates::Schwarzschild,
138 : double log_enthalpy_at_outer_radius = 0.0,
139 : double absolute_tolerance = 1.e-18, double relative_tolerance = 1.0e-14);
140 :
141 0 : TovSolution() = default;
142 0 : TovSolution(const TovSolution& /*rhs*/) = default;
143 0 : TovSolution& operator=(const TovSolution& /*rhs*/) = default;
144 0 : TovSolution(TovSolution&& /*rhs*/) = default;
145 0 : TovSolution& operator=(TovSolution&& /*rhs*/) = default;
146 0 : ~TovSolution() = default;
147 :
148 : /// The type of radial coordinate.
149 : ///
150 : /// \see RelativisticEuler::Solutions::TovCoordinates
151 1 : TovCoordinates coordinate_system() const { return coordinate_system_; }
152 :
153 : /// \brief The outer radius of the solution.
154 : ///
155 : /// This is the outer radius in the specified `coordinate_system()`, i.e.,
156 : /// areal or isotropic.
157 : ///
158 : /// \note This is the radius at which `log_specific_enthalpy` is equal
159 : /// to the value of `log_enthalpy_at_outer_radius` that was given when
160 : /// constructing this TovSolution
161 1 : double outer_radius() const { return outer_radius_; }
162 :
163 : /// The total mass \f$m(R)\f$, where \f$R\f$ is the outer radius.
164 1 : double total_mass() const { return total_mass_; }
165 :
166 : /*!
167 : * \brief The injection energy \f$\mathcal{E}=\alpha(r=R)=\sqrt{1-2M/R}\f$.
168 : *
169 : * The injection energy of the TOV solution is
170 : *
171 : * \f{equation}
172 : * \mathcal{E} = -h k^a u_a = h \alpha
173 : * \text{,}
174 : * \f}
175 : *
176 : * where \f$\boldsymbol{k} = \partial_t\f$ is a Killing vector of the static
177 : * solution, \f$h\f$ is the specific enthalpy, \f$u_a\f$ is the fluid
178 : * four-velocity, and \f$\alpha\f$ is the lapse (see, e.g., Eqs. (2.19) and
179 : * (4.2) in \cite Moldenhauer2014yaa). Since the TOV solution is static, the
180 : * injection energy is conserved not only along stream lines but throughout
181 : * the star,
182 : *
183 : * \f{equation}
184 : * \nabla_a \mathcal{E} = 0
185 : * \text{.}
186 : * \f}
187 : *
188 : * Therefore,
189 : *
190 : * \f{equation}
191 : * \mathcal{E} = \alpha(r=R) = \sqrt{1 - 2M/R}
192 : * \f}
193 : *
194 : * by evaluating the injection energy at the outer (areal) radius \f$R\f$,
195 : * where \f$h=1\f$ and where we match the lapse to the outer Schwarzschild
196 : * solution. The conservation also implies
197 : *
198 : * \f{equation}
199 : * \alpha = \mathcal{E} / h
200 : * \f}
201 : *
202 : * throughout the star.
203 : */
204 1 : double injection_energy() const { return injection_energy_; }
205 :
206 : /// \brief The mass inside the given radius over the areal radius,
207 : /// \f$\frac{m(r)}{r}\f$
208 : ///
209 : /// The argument to this function is the radius in the `coordinate_system()`,
210 : /// i.e., areal (Schwarzschild) or isotropic radius. The denominator \f$r\f$
211 : /// in the return value is always the areal (Schwarzschild) radius. You can
212 : /// use the conformal factor \f$\psi=\sqrt{r / \bar{r}}\f$ returned by the
213 : /// `conformal_factor()` function to obtain the mass over the isotropic
214 : /// radius, or the mass alone. The reason for this choice is that we represent
215 : /// the solution internally as the mass over the areal radius, so this is the
216 : /// most numerically precise quantity from which other quantities can be
217 : /// derived.
218 : ///
219 : /// \note `r` should be non-negative and not greater than `outer_radius()`.
220 : template <typename DataType>
221 1 : DataType mass_over_radius(const DataType& r) const;
222 :
223 : /// \brief The log of the specific enthalpy at the given radius
224 : ///
225 : /// \note `r` should be non-negative and not greater than `outer_radius()`
226 : template <typename DataType>
227 1 : DataType log_specific_enthalpy(const DataType& r) const;
228 :
229 : /// \brief The conformal factor \f$\psi=\sqrt{r / \bar{r}}\f$.
230 : ///
231 : /// The conformal factor is computed only when the `coordinate_system()` is
232 : /// `RelativisticEuler::Solution::TovCoordinates::Isotropic`. Otherwise, it is
233 : /// an error to call this function.
234 : ///
235 : /// \note `r` should be non-negative and not greater than `outer_radius()`
236 : template <typename DataType>
237 1 : DataType conformal_factor(const DataType& r) const;
238 :
239 0 : const intrp::CubicSpline& mass_over_radius_interpolant() const {
240 : return mass_over_radius_interpolant_;
241 : }
242 0 : const intrp::CubicSpline& log_specific_enthalpy_interpolant() const {
243 : return log_enthalpy_interpolant_;
244 : }
245 0 : const intrp::CubicSpline& conformal_factor_interpolant() const {
246 : return conformal_factor_interpolant_;
247 : }
248 :
249 : // NOLINTNEXTLINE(google-runtime-references)
250 0 : void pup(PUP::er& p);
251 :
252 : private:
253 : template <TovCoordinates CoordSystem>
254 0 : void integrate(
255 : const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
256 : const double central_mass_density,
257 : const double log_enthalpy_at_outer_radius,
258 : const double absolute_tolerance, const double relative_tolerance);
259 :
260 0 : TovCoordinates coordinate_system_{};
261 0 : double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
262 0 : double total_mass_{std::numeric_limits<double>::signaling_NaN()};
263 0 : double injection_energy_{std::numeric_limits<double>::signaling_NaN()};
264 0 : intrp::CubicSpline mass_over_radius_interpolant_;
265 0 : intrp::CubicSpline log_enthalpy_interpolant_;
266 0 : intrp::CubicSpline conformal_factor_interpolant_;
267 : };
268 :
269 : } // namespace RelativisticEuler::Solutions
|
