Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 :
9 : #include "DataStructures/CachedTempBuffer.hpp"
10 : #include "DataStructures/DataBox/Tag.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Options/String.hpp"
15 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
16 : #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Solutions.hpp"
17 : #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Tov.hpp"
18 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" // IWYU pragma: keep
19 : #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp" // IWYU pragma: keep
20 : #include "PointwiseFunctions/Hydro/Tags.hpp"
21 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
22 : #include "Utilities/Gsl.hpp"
23 : #include "Utilities/MakeWithValue.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; // IWYU pragma: keep
31 : } // namespace PUP
32 : /// \endcond
33 :
34 : namespace RelativisticEuler::Solutions {
35 : namespace tov_detail {
36 :
37 : enum class StarRegion { Center, Interior, Exterior };
38 :
39 : namespace Tags {
40 : template <typename DataType>
41 : struct MassOverRadius : db::SimpleTag {
42 : using type = Scalar<DataType>;
43 : };
44 : template <typename DataType>
45 : struct LogSpecificEnthalpy : db::SimpleTag {
46 : using type = Scalar<DataType>;
47 : };
48 : template <typename DataType>
49 : struct ConformalFactor : db::SimpleTag {
50 : using type = Scalar<DataType>;
51 : };
52 : template <typename DataType>
53 : struct DrConformalFactor : db::SimpleTag {
54 : using type = Scalar<DataType>;
55 : };
56 : template <typename DataType>
57 : struct ArealRadius : db::SimpleTag {
58 : using type = Scalar<DataType>;
59 : };
60 : template <typename DataType>
61 : struct DrArealRadius : db::SimpleTag {
62 : using type = Scalar<DataType>;
63 : };
64 : template <typename DataType>
65 : struct DrPressure : db::SimpleTag {
66 : using type = Scalar<DataType>;
67 : };
68 : template <typename DataType>
69 : struct MetricTimePotential : db::SimpleTag {
70 : using type = Scalar<DataType>;
71 : };
72 : template <typename DataType>
73 : struct DrMetricTimePotential : db::SimpleTag {
74 : using type = Scalar<DataType>;
75 : };
76 : template <typename DataType>
77 : struct MetricRadialPotential : db::SimpleTag {
78 : using type = Scalar<DataType>;
79 : };
80 : template <typename DataType>
81 : struct DrMetricRadialPotential : db::SimpleTag {
82 : using type = Scalar<DataType>;
83 : };
84 : template <typename DataType>
85 : struct MetricAngularPotential : db::SimpleTag {
86 : using type = Scalar<DataType>;
87 : };
88 : template <typename DataType>
89 : struct DrMetricAngularPotential : db::SimpleTag {
90 : using type = Scalar<DataType>;
91 : };
92 : } // namespace Tags
93 :
94 : template <typename DataType>
95 : using TovVariablesCache = cached_temp_buffer_from_typelist<tmpl::list<
96 : Tags::MassOverRadius<DataType>, Tags::LogSpecificEnthalpy<DataType>,
97 : Tags::ConformalFactor<DataType>, Tags::DrConformalFactor<DataType>,
98 : Tags::ArealRadius<DataType>, Tags::DrArealRadius<DataType>,
99 : hydro::Tags::SpecificEnthalpy<DataType>,
100 : hydro::Tags::RestMassDensity<DataType>,
101 : hydro::Tags::ElectronFraction<DataType>, hydro::Tags::Pressure<DataType>,
102 : Tags::DrPressure<DataType>,
103 : ::Tags::deriv<hydro::Tags::Pressure<DataType>, tmpl::size_t<3>,
104 : Frame::Inertial>,
105 : hydro::Tags::SpecificInternalEnergy<DataType>,
106 : hydro::Tags::Temperature<DataType>, Tags::MetricTimePotential<DataType>,
107 : Tags::DrMetricTimePotential<DataType>,
108 : Tags::MetricRadialPotential<DataType>,
109 : Tags::DrMetricRadialPotential<DataType>,
110 : Tags::MetricAngularPotential<DataType>,
111 : Tags::DrMetricAngularPotential<DataType>,
112 : hydro::Tags::LorentzFactor<DataType>,
113 : hydro::Tags::SpatialVelocity<DataType, 3>,
114 : hydro::Tags::MagneticField<DataType, 3>,
115 : hydro::Tags::DivergenceCleaningField<DataType>, gr::Tags::Lapse<DataType>,
116 : ::Tags::dt<gr::Tags::Lapse<DataType>>,
117 : ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>, Frame::Inertial>,
118 : gr::Tags::Shift<DataType, 3>, ::Tags::dt<gr::Tags::Shift<DataType, 3>>,
119 : ::Tags::deriv<gr::Tags::Shift<DataType, 3>, tmpl::size_t<3>,
120 : Frame::Inertial>,
121 : gr::Tags::SpatialMetric<DataType, 3>,
122 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, 3>>,
123 : ::Tags::deriv<gr::Tags::SpatialMetric<DataType, 3>, tmpl::size_t<3>,
124 : Frame::Inertial>,
125 : gr::Tags::SqrtDetSpatialMetric<DataType>,
126 : gr::Tags::ExtrinsicCurvature<DataType, 3>,
127 : gr::Tags::InverseSpatialMetric<DataType, 3>>>;
128 :
129 : template <typename DataType, StarRegion Region>
130 : struct TovVariables {
131 : static constexpr size_t Dim = 3;
132 : using Cache = TovVariablesCache<DataType>;
133 :
134 : TovVariables(const TovVariables&) = default;
135 : TovVariables& operator=(const TovVariables&) = default;
136 : TovVariables(TovVariables&&) = default;
137 : TovVariables& operator=(TovVariables&&) = default;
138 : virtual ~TovVariables() = default;
139 : TovVariables(
140 : const tnsr::I<DataType, 3>& local_coords, const DataType& local_radius,
141 : const RelativisticEuler::Solutions::TovSolution& local_radial_solution,
142 : const EquationsOfState::EquationOfState<true, 1>& local_eos)
143 : : coords(local_coords),
144 : radius(local_radius),
145 : radial_solution(local_radial_solution),
146 : eos(local_eos) {}
147 :
148 : const tnsr::I<DataType, 3>& coords;
149 : const DataType& radius;
150 : const RelativisticEuler::Solutions::TovSolution& radial_solution;
151 : const EquationsOfState::EquationOfState<true, 1>& eos;
152 :
153 : void operator()(gsl::not_null<Scalar<DataType>*> mass_over_radius,
154 : gsl::not_null<Cache*> cache,
155 : Tags::MassOverRadius<DataType> /*meta*/) const;
156 : void operator()(gsl::not_null<Scalar<DataType>*> log_specific_enthalpy,
157 : gsl::not_null<Cache*> cache,
158 : Tags::LogSpecificEnthalpy<DataType> /*meta*/) const;
159 : void operator()(gsl::not_null<Scalar<DataType>*> conformal_factor,
160 : gsl::not_null<Cache*> cache,
161 : Tags::ConformalFactor<DataType> /*meta*/) const;
162 : void operator()(gsl::not_null<Scalar<DataType>*> dr_conformal_factor,
163 : gsl::not_null<Cache*> cache,
164 : Tags::DrConformalFactor<DataType> /*meta*/) const;
165 : void operator()(gsl::not_null<Scalar<DataType>*> areal_radius,
166 : gsl::not_null<Cache*> cache,
167 : Tags::ArealRadius<DataType> /*meta*/) const;
168 : void operator()(gsl::not_null<Scalar<DataType>*> dr_areal_radius,
169 : gsl::not_null<Cache*> cache,
170 : Tags::DrArealRadius<DataType> /*meta*/) const;
171 : void operator()(gsl::not_null<Scalar<DataType>*> specific_enthalpy,
172 : gsl::not_null<Cache*> cache,
173 : hydro::Tags::SpecificEnthalpy<DataType> /*meta*/) const;
174 : void operator()(gsl::not_null<Scalar<DataType>*> temperature,
175 : gsl::not_null<Cache*> cache,
176 : hydro::Tags::Temperature<DataType> /*meta*/) const;
177 : void operator()(gsl::not_null<Scalar<DataType>*> rest_mass_density,
178 : gsl::not_null<Cache*> cache,
179 : hydro::Tags::RestMassDensity<DataType> /*meta*/) const;
180 : void operator()(gsl::not_null<Scalar<DataType>*> electron_fraction,
181 : gsl::not_null<Cache*> cache,
182 : hydro::Tags::ElectronFraction<DataType> /*meta*/) const;
183 : void operator()(gsl::not_null<Scalar<DataType>*> pressure,
184 : gsl::not_null<Cache*> cache,
185 : hydro::Tags::Pressure<DataType> /*meta*/) const;
186 : void operator()(gsl::not_null<Scalar<DataType>*> dr_pressure,
187 : gsl::not_null<Cache*> cache,
188 : Tags::DrPressure<DataType> /*meta*/) const;
189 : void operator()(
190 : gsl::not_null<tnsr::i<DataType, 3>*> deriv_pressure,
191 : gsl::not_null<Cache*> cache,
192 : ::Tags::deriv<hydro::Tags::Pressure<DataType>, tmpl::size_t<3>,
193 : Frame::Inertial> /*meta*/) const;
194 : void operator()(gsl::not_null<Scalar<DataType>*> specific_internal_energy,
195 : gsl::not_null<Cache*> cache,
196 : hydro::Tags::SpecificInternalEnergy<DataType> /*meta*/) const;
197 : void operator()(gsl::not_null<Scalar<DataType>*> metric_time_potential,
198 : gsl::not_null<Cache*> cache,
199 : Tags::MetricTimePotential<DataType> /*meta*/) const;
200 : void operator()(gsl::not_null<Scalar<DataType>*> dr_metric_time_potential,
201 : gsl::not_null<Cache*> cache,
202 : Tags::DrMetricTimePotential<DataType> /*meta*/) const;
203 : void operator()(gsl::not_null<Scalar<DataType>*> metric_radial_potential,
204 : gsl::not_null<Cache*> cache,
205 : Tags::MetricRadialPotential<DataType> /*meta*/) const;
206 : void operator()(gsl::not_null<Scalar<DataType>*> dr_metric_radial_potential,
207 : gsl::not_null<Cache*> cache,
208 : Tags::DrMetricRadialPotential<DataType> /*meta*/) const;
209 : void operator()(gsl::not_null<Scalar<DataType>*> metric_angular_potential,
210 : gsl::not_null<Cache*> cache,
211 : Tags::MetricAngularPotential<DataType> /*meta*/) const;
212 : void operator()(gsl::not_null<Scalar<DataType>*> dr_metric_angular_potential,
213 : gsl::not_null<Cache*> cache,
214 : Tags::DrMetricAngularPotential<DataType> /*meta*/) const;
215 : void operator()(gsl::not_null<Scalar<DataType>*> lorentz_factor,
216 : gsl::not_null<Cache*> cache,
217 : hydro::Tags::LorentzFactor<DataType> /*meta*/) const;
218 : void operator()(gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
219 : gsl::not_null<Cache*> cache,
220 : hydro::Tags::SpatialVelocity<DataType, 3> /*meta*/) const;
221 : virtual void operator()(
222 : gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
223 : gsl::not_null<Cache*> cache,
224 : hydro::Tags::MagneticField<DataType, 3> /*meta*/) const;
225 : void operator()(
226 : gsl::not_null<Scalar<DataType>*> div_cleaning_field,
227 : gsl::not_null<Cache*> cache,
228 : hydro::Tags::DivergenceCleaningField<DataType> /*meta*/) const;
229 : void operator()(gsl::not_null<Scalar<DataType>*> lapse,
230 : gsl::not_null<Cache*> cache,
231 : gr::Tags::Lapse<DataType> /*meta*/) const;
232 : void operator()(gsl::not_null<Scalar<DataType>*> dt_lapse,
233 : gsl::not_null<Cache*> cache,
234 : ::Tags::dt<gr::Tags::Lapse<DataType>> /*meta*/) const;
235 : void operator()(gsl::not_null<tnsr::i<DataType, 3>*> deriv_lapse,
236 : gsl::not_null<Cache*> cache,
237 : ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>,
238 : Frame::Inertial> /*meta*/) const;
239 : void operator()(gsl::not_null<tnsr::I<DataType, 3>*> shift,
240 : gsl::not_null<Cache*> cache,
241 : gr::Tags::Shift<DataType, 3> /*meta*/) const;
242 : void operator()(gsl::not_null<tnsr::I<DataType, 3>*> dt_shift,
243 : gsl::not_null<Cache*> cache,
244 : ::Tags::dt<gr::Tags::Shift<DataType, 3>> /*meta*/) const;
245 : void operator()(gsl::not_null<tnsr::iJ<DataType, 3>*> deriv_shift,
246 : gsl::not_null<Cache*> cache,
247 : ::Tags::deriv<gr::Tags::Shift<DataType, 3>, tmpl::size_t<3>,
248 : Frame::Inertial> /*meta*/) const;
249 : void operator()(gsl::not_null<tnsr::ii<DataType, 3>*> spatial_metric,
250 : gsl::not_null<Cache*> cache,
251 : gr::Tags::SpatialMetric<DataType, 3> /*meta*/) const;
252 : void operator()(
253 : gsl::not_null<tnsr::ii<DataType, 3>*> dt_spatial_metric,
254 : gsl::not_null<Cache*> cache,
255 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, 3>> /*meta*/) const;
256 : void operator()(
257 : gsl::not_null<tnsr::ijj<DataType, 3>*> deriv_spatial_metric,
258 : gsl::not_null<Cache*> cache,
259 : ::Tags::deriv<gr::Tags::SpatialMetric<DataType, 3>, tmpl::size_t<3>,
260 : Frame::Inertial> /*meta*/) const;
261 : void operator()(gsl::not_null<Scalar<DataType>*> sqrt_det_spatial_metric,
262 : gsl::not_null<Cache*> cache,
263 : gr::Tags::SqrtDetSpatialMetric<DataType> /*meta*/) const;
264 : void operator()(gsl::not_null<tnsr::ii<DataType, 3>*> extrinsic_curvature,
265 : gsl::not_null<Cache*> cache,
266 : gr::Tags::ExtrinsicCurvature<DataType, 3> /*meta*/) const;
267 : void operator()(gsl::not_null<tnsr::II<DataType, 3>*> inv_spatial_metric,
268 : gsl::not_null<Cache*> cache,
269 : gr::Tags::InverseSpatialMetric<DataType, 3> /*meta*/) const;
270 : };
271 : } // namespace tov_detail
272 :
273 : /*!
274 : * \brief A static spherically symmetric star
275 : *
276 : * An analytic solution for a static, spherically-symmetric star found by
277 : * solving the Tolman-Oppenheimer-Volkoff (TOV) equations. The equation of
278 : * state is assumed to be that of a polytropic fluid.
279 : *
280 : * If the spherically symmetric metric is written as
281 : *
282 : * \f[
283 : * ds^2 = - e^{2 \Phi_t} dt^2 + e^{2 \Phi_r} dr^2 + e^{2 \Phi_\Omega} r^2
284 : * d\Omega^2
285 : * \f]
286 : *
287 : * where \f$r = \delta_{mn} x^m x^n\f$ is the radial coordinate and
288 : * \f$\Phi_t\f$, \f$\Phi_r\f$, and \f$\Phi_\Omega\f$ are the metric potentials,
289 : * then the lapse, shift, and spatial metric in Cartesian coordinates are
290 : *
291 : * \f{align*}
292 : * \alpha &= e^{\Phi_t} \\
293 : * \beta^i &= 0 \\
294 : * \gamma_{ij} &= \delta_{ij} e^{2 \Phi_\Omega} + \delta_{im} \delta_{jn}
295 : * \frac{x^m x^n}{r^2} \left( e^{2 \Phi_r} - e^{2 \Phi_\Omega} \right)
296 : * \f}
297 : *
298 : * We solve the TOV equations with the method implemented in
299 : * `RelativisticEuler::Solutions::TovSolution`. It provides the areal
300 : * mass-over-radius \f$m(r)/r\f$ and the log of the specific enthalpy
301 : * \f$\log{h}\f$. In areal (Schwarzschild) coordinates the spatial metric
302 : * potentials are
303 : *
304 : * \f{align}
305 : * e^{\Phi_r} &= \left(1 - \frac{2m}{r}\right)^{-1/2} \\
306 : * e^{\Phi_\Omega} &= 1
307 : * \f}
308 : *
309 : * In isotropic coordinates the spatial metric potentials are
310 : *
311 : * \begin{equation}
312 : * e^{2\Phi_r} = e^{2\Phi_\Omega} = \psi^4
313 : * \text{,}
314 : * \end{equation}
315 : *
316 : * where $\psi = \sqrt{r / \bar{r}}$ is the conformal factor, $r$ is the areal
317 : * (Schwarzschild) radius and $\bar{r}$ is the isotropic radius. See
318 : * `RelativisticEuler::Solutions::TovSolution` for details.
319 : *
320 : * \warning Isotropic coordinates should be used because the metric derivatives
321 : * are smooth. Otherwise the grid will over-compensate with finite difference
322 : * cells.
323 : */
324 :
325 1 : class TovStar : public virtual evolution::initial_data::InitialData,
326 : public MarkAsAnalyticSolution,
327 : public AnalyticSolution<3> {
328 : public:
329 0 : using equation_of_state_type = EquationsOfState::EquationOfState<true, 1>;
330 :
331 : /// The central density of the star.
332 1 : struct CentralDensity {
333 0 : using type = double;
334 0 : static constexpr Options::String help = {
335 : "The central density of the star."};
336 0 : static type lower_bound() { return 0.; }
337 : };
338 :
339 : /// Areal (Schwarzschild) or isotropic coordinates
340 1 : struct Coordinates {
341 0 : using type = RelativisticEuler::Solutions::TovCoordinates;
342 0 : static constexpr Options::String help = {
343 : "Areal ('Schwarzschild') or 'Isotropic' coordinates."};
344 : };
345 :
346 0 : static constexpr size_t volume_dim = 3_st;
347 :
348 0 : using options =
349 : tmpl::list<CentralDensity,
350 : hydro::OptionTags::InitialDataEquationOfState<true, 1>,
351 : Coordinates>;
352 :
353 0 : static constexpr Options::String help = {
354 : "A static, spherically-symmetric star found by solving the \n"
355 : "Tolman-Oppenheimer-Volkoff (TOV) equations, with a given central \n"
356 : "density and equation of state."};
357 :
358 0 : TovStar() = default;
359 0 : TovStar(const TovStar& /*rhs*/);
360 0 : TovStar& operator=(const TovStar& /*rhs*/);
361 0 : TovStar(TovStar&& /*rhs*/) = default;
362 0 : TovStar& operator=(TovStar&& /*rhs*/) = default;
363 0 : ~TovStar() override = default;
364 0 : TovStar(double central_rest_mass_density,
365 : std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
366 : equation_of_state,
367 : const RelativisticEuler::Solutions::TovCoordinates coordinate_system =
368 : RelativisticEuler::Solutions::TovCoordinates::Schwarzschild);
369 :
370 0 : auto get_clone() const
371 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
372 :
373 : /// \cond
374 : explicit TovStar(CkMigrateMessage* msg);
375 : using PUP::able::register_constructor;
376 : WRAPPED_PUPable_decl_template(TovStar);
377 : /// \endcond
378 :
379 : /// Retrieve a collection of variables at `(x, t)`
380 : template <typename DataType, typename... Tags>
381 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
382 : const double /*t*/,
383 : tmpl::list<Tags...> /*meta*/) const {
384 : return variables_impl<tov_detail::TovVariables>(x, tmpl::list<Tags...>{});
385 : }
386 :
387 : /// NOLINTNEXTLINE(google-runtime-references)
388 1 : void pup(PUP::er& /*p*/) override;
389 :
390 0 : const EquationsOfState::EquationOfState<true, 1>& equation_of_state() const {
391 : return *equation_of_state_;
392 : }
393 :
394 : /// The radial profile of the star
395 1 : const RelativisticEuler::Solutions::TovSolution& radial_solution() const {
396 : return radial_solution_;
397 : }
398 :
399 : protected:
400 : template <template <class, tov_detail::StarRegion> class VarsComputer,
401 : typename DataType, typename... Tags, typename... VarsComputerArgs>
402 0 : tuples::TaggedTuple<Tags...> variables_impl(
403 : const tnsr::I<DataType, 3>& x, tmpl::list<Tags...> /*meta*/,
404 : VarsComputerArgs&&... vars_computer_args) const {
405 : const double outer_radius = radial_solution_.outer_radius();
406 : const double center_radius_cutoff = 1.e-30 * outer_radius;
407 : const DataType radius = get(magnitude(x));
408 : // Dispatch interior and exterior regions of the star.
409 : // - Include the equality in the conditions below for the outer radius so
410 : // the `DataVector` variants are preferred over the pointwise `double`
411 : // variant when possible.
412 : // - Order the conditions so the cheaper exterior solution is preferred over
413 : // the interior.
414 : // - A `DataVector` variant for the center of the star is not needed because
415 : // it's only a single point.
416 : if (min(radius) >= outer_radius) {
417 : // All points are outside the star. This could be replaced by a
418 : // Schwarzschild solution.
419 : using ExteriorVarsComputer =
420 : VarsComputer<DataType, tov_detail::StarRegion::Exterior>;
421 : typename ExteriorVarsComputer::Cache cache{get_size(radius)};
422 : ExteriorVarsComputer computer{
423 : x, radius, radial_solution_, *equation_of_state_,
424 : std::forward<VarsComputerArgs>(vars_computer_args)...};
425 : return {cache.get_var(computer, Tags{})...};
426 : } else if (max(radius) <= outer_radius and
427 : min(radius) > center_radius_cutoff) {
428 : // All points are in the star interior, but not at the center
429 : using InteriorVarsComputer =
430 : VarsComputer<DataType, tov_detail::StarRegion::Interior>;
431 : typename InteriorVarsComputer::Cache cache{get_size(radius)};
432 : InteriorVarsComputer computer{
433 : x, radius, radial_solution_, *equation_of_state_,
434 : std::forward<VarsComputerArgs>(vars_computer_args)...};
435 : return {cache.get_var(computer, Tags{})...};
436 : } else {
437 : // Points can be at the center, in the interior, or outside the star, so
438 : // check each point individually
439 : const size_t num_points = get_size(radius);
440 : tuples::TaggedTuple<Tags...> vars{typename Tags::type{num_points}...};
441 : const auto get_var = [&vars](const size_t point_index, auto& local_cache,
442 : const auto& local_computer, auto tag_v) {
443 : using tag = std::decay_t<decltype(tag_v)>;
444 : if constexpr (std::is_same_v<DataType, DataVector>) {
445 : using tags_double = typename VarsComputer<
446 : double, tov_detail::StarRegion::Exterior>::Cache::tags_list;
447 : using tags_dv = typename VarsComputer<
448 : DataVector, tov_detail::StarRegion::Exterior>::Cache::tags_list;
449 : using tag_double =
450 : tmpl::at<tags_double, tmpl::index_of<tags_dv, tag>>;
451 : const auto& tensor =
452 : local_cache.get_var(local_computer, tag_double{});
453 : for (size_t component = 0; component < tensor.size(); ++component) {
454 : get<tag>(vars)[component][point_index] = tensor[component];
455 : }
456 : } else {
457 : (void)point_index;
458 : get<tag>(vars) = local_cache.get_var(local_computer, tag{});
459 : }
460 : return '0';
461 : };
462 : using CenterVarsComputer =
463 : VarsComputer<double, tov_detail::StarRegion::Center>;
464 : using InteriorVarsComputer =
465 : VarsComputer<double, tov_detail::StarRegion::Interior>;
466 : using ExteriorVarsComputer =
467 : VarsComputer<double, tov_detail::StarRegion::Exterior>;
468 : for (size_t i = 0; i < num_points; ++i) {
469 : const tnsr::I<double, 3> x_i{
470 : {{get_element(get<0>(x), i), get_element(get<1>(x), i),
471 : get_element(get<2>(x), i)}}};
472 : if (get_element(radius, i) > outer_radius) {
473 : typename ExteriorVarsComputer::Cache cache{1};
474 : ExteriorVarsComputer computer{
475 : x_i, get_element(radius, i), radial_solution_,
476 : *equation_of_state_,
477 : std::forward<VarsComputerArgs>(vars_computer_args)...};
478 : expand_pack(get_var(i, cache, computer, Tags{})...);
479 : } else if (get_element(radius, i) > center_radius_cutoff) {
480 : typename InteriorVarsComputer::Cache cache{1};
481 : InteriorVarsComputer computer{
482 : x_i, get_element(radius, i), radial_solution_,
483 : *equation_of_state_,
484 : std::forward<VarsComputerArgs>(vars_computer_args)...};
485 : expand_pack(get_var(i, cache, computer, Tags{})...);
486 : } else {
487 : typename CenterVarsComputer::Cache cache{1};
488 : CenterVarsComputer computer{
489 : x_i, get_element(radius, i), radial_solution_,
490 : *equation_of_state_,
491 : std::forward<VarsComputerArgs>(vars_computer_args)...};
492 : expand_pack(get_var(i, cache, computer, Tags{})...);
493 : }
494 : }
495 : return vars;
496 : }
497 : }
498 :
499 : public:
500 : template <typename DataType>
501 0 : using tags = tmpl::list_difference<
502 : typename tov_detail::TovVariablesCache<DataType>::tags_list,
503 : tmpl::list<
504 : // Remove internal tags, which may not be available in the full domain
505 : tov_detail::Tags::MassOverRadius<DataType>,
506 : tov_detail::Tags::LogSpecificEnthalpy<DataType>,
507 : tov_detail::Tags::ConformalFactor<DataType>,
508 : tov_detail::Tags::DrConformalFactor<DataType>,
509 : tov_detail::Tags::ArealRadius<DataType>,
510 : tov_detail::Tags::DrArealRadius<DataType>,
511 : tov_detail::Tags::DrPressure<DataType>,
512 : tov_detail::Tags::MetricTimePotential<DataType>,
513 : tov_detail::Tags::DrMetricTimePotential<DataType>,
514 : tov_detail::Tags::MetricRadialPotential<DataType>,
515 : tov_detail::Tags::DrMetricRadialPotential<DataType>,
516 : tov_detail::Tags::MetricAngularPotential<DataType>,
517 : tov_detail::Tags::DrMetricAngularPotential<DataType>>>;
518 :
519 : private:
520 0 : friend bool operator==(const TovStar& lhs, const TovStar& rhs);
521 :
522 0 : double central_rest_mass_density_ =
523 : std::numeric_limits<double>::signaling_NaN();
524 0 : std::unique_ptr<equation_of_state_type> equation_of_state_;
525 0 : RelativisticEuler::Solutions::TovCoordinates coordinate_system_{};
526 0 : RelativisticEuler::Solutions::TovSolution radial_solution_{};
527 : };
528 :
529 0 : bool operator!=(const TovStar& lhs, const TovStar& rhs);
530 : } // namespace RelativisticEuler::Solutions
|
