Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <limits>
8 :
9 : #include "DataStructures/TaggedTuple.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Options/Context.hpp"
12 : #include "Options/Options.hpp"
13 : #include "Options/String.hpp"
14 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
15 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
16 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
17 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
18 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
19 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
20 : #include "Utilities/Serialization/CharmPupable.hpp"
21 : #include "Utilities/TMPL.hpp"
22 :
23 : /// \cond
24 : namespace PUP {
25 : class er;
26 : } // namespace PUP
27 : /// \endcond
28 :
29 : namespace grmhd::AnalyticData {
30 :
31 : /*!
32 : * \brief Analytic initial data for a cylindrical or spherical blast wave.
33 : *
34 : * This class implements analytic initial data for a cylindrical blast wave,
35 : * as described, e.g., in \cite Kidder2016hev Sec. 6.2.3.
36 : * A uniform magnetic field threads an ideal fluid. The solution begins with
37 : * material at fixed (typically high) density and pressure at rest inside a
38 : * cylinder of radius \f$r < r_{\rm in}\f$ and material at fixed (typically low)
39 : * density and pressure at rest in a cylindrical shell with radius
40 : * \f$r > r_{\rm out}\f$. In the region \f$ r_{\rm in} < r < r_{\rm out}\f$,
41 : * the solution transitions such that the logarithms of the density and
42 : * pressure vary linearly. E.g., if \f$\rho(r < r_{\rm in}) = \rho_{\rm in}\f$
43 : * and \f$\rho(r > r_{\rm out}) = \rho_{\rm out}\f$, then
44 : * \f[
45 : * \log \rho = [(r_{\rm in} - r) \log(\rho_{\rm out})
46 : * + (r - r_{\rm out}) \log(\rho_{\rm in})]
47 : * / (r_{\rm in} - r_{\rm out}).
48 : * \f]
49 : * Note that the cylinder's axis is the \f$z\f$ axis. To evolve this analytic
50 : * initial data, use a cubic or cylindrical domain with periodic boundary
51 : * conditions applied to the outer boundaries whose normals are parallel or
52 : * antiparallel to the z axis. In the transverse (e.g., x and y) dimensions, the
53 : * domain should be large enough that the blast wave doesn't reach the boundary
54 : * at the final time. E.g., if `InnerRadius = 0.8`, `OuterRadius = 1.0`, and
55 : * the final time is 4.0, a good domain extends from `(x,y)=(-6.0, -6.0)` to
56 : * `(x,y)=(6.0, 6.0)`.
57 : *
58 : * An analogous problem with spherical geometry has also been used
59 : * \cite CerdaDuran2007 \cite CerdaDuran2008 \cite Cipolletta2019. The magnetic
60 : * field is chosen to be in the z-direction instead of the x-direction.
61 : */
62 1 : class BlastWave : public evolution::initial_data::InitialData,
63 : public MarkAsAnalyticData,
64 : public AnalyticDataBase,
65 : public hydro::TemperatureInitialization<BlastWave> {
66 : public:
67 0 : using equation_of_state_type = EquationsOfState::IdealFluid<true>;
68 :
69 0 : enum class Geometry { Cylindrical, Spherical };
70 :
71 : /// Inside InnerRadius, density is InnerDensity.
72 1 : struct InnerRadius {
73 0 : using type = double;
74 0 : static constexpr Options::String help = {
75 : "Inside InnerRadius, density is InnerDensity."};
76 0 : static type lower_bound() { return 0.0; }
77 : };
78 : /// Outside OuterRadius, density is OuterDensity.
79 1 : struct OuterRadius {
80 0 : using type = double;
81 0 : static constexpr Options::String help = {
82 : "Outside OuterRadius, density is OuterDensity."};
83 0 : static type lower_bound() { return 0.0; }
84 : };
85 : /// Density at radii less than InnerRadius.
86 1 : struct InnerDensity {
87 0 : using type = double;
88 0 : static constexpr Options::String help = {
89 : "Density at radii less than InnerRadius."};
90 0 : static type lower_bound() { return 0.0; }
91 : };
92 : /// Density at radii greater than OuterRadius.
93 1 : struct OuterDensity {
94 0 : using type = double;
95 0 : static constexpr Options::String help = {
96 : "Density at radii greater than OuterRadius."};
97 0 : static type lower_bound() { return 0.0; }
98 : };
99 : /// Pressure at radii less than InnerRadius.
100 1 : struct InnerPressure {
101 0 : using type = double;
102 0 : static constexpr Options::String help = {
103 : "Pressure at radii less than InnerRadius."};
104 0 : static type lower_bound() { return 0.0; }
105 : };
106 : /// Pressure at radii greater than OuterRadius.
107 1 : struct OuterPressure {
108 0 : using type = double;
109 0 : static constexpr Options::String help = {
110 : "Pressure at radii greater than OuterRadius."};
111 0 : static type lower_bound() { return 0.0; }
112 : };
113 : /// The x,y,z components of the uniform magnetic field threading the matter.
114 1 : struct MagneticField {
115 0 : using type = std::array<double, 3>;
116 0 : static constexpr Options::String help = {
117 : "The x,y,z components of the uniform magnetic field."};
118 : };
119 : /// The adiabatic index of the ideal fluid.
120 1 : struct AdiabaticIndex {
121 0 : using type = double;
122 0 : static constexpr Options::String help = {
123 : "The adiabatic index of the ideal fluid."};
124 0 : static type lower_bound() { return 1.0; }
125 : };
126 : /// The geometry of the blast wave, i.e. Cylindrical or Spherical.
127 1 : struct GeometryOption {
128 0 : static std::string name() { return "Geometry"; }
129 0 : using type = Geometry;
130 0 : static constexpr Options::String help = {
131 : "The geometry of the blast wave, i.e. Cylindrical or Spherical."};
132 : };
133 :
134 0 : using options = tmpl::list<InnerRadius, OuterRadius, InnerDensity,
135 : OuterDensity, InnerPressure, OuterPressure,
136 : MagneticField, AdiabaticIndex, GeometryOption>;
137 :
138 0 : static constexpr Options::String help = {
139 : "Cylindrical or spherical blast wave analytic initial data."};
140 :
141 0 : BlastWave() = default;
142 0 : BlastWave(const BlastWave& /*rhs*/) = default;
143 0 : BlastWave& operator=(const BlastWave& /*rhs*/) = default;
144 0 : BlastWave(BlastWave&& /*rhs*/) = default;
145 0 : BlastWave& operator=(BlastWave&& /*rhs*/) = default;
146 0 : ~BlastWave() override = default;
147 :
148 0 : BlastWave(double inner_radius, double outer_radius, double inner_density,
149 : double outer_density, double inner_pressure, double outer_pressure,
150 : const std::array<double, 3>& magnetic_field, double adiabatic_index,
151 : Geometry geometry, const Options::Context& context = {});
152 :
153 0 : auto get_clone() const
154 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
155 :
156 : /// \cond
157 : explicit BlastWave(CkMigrateMessage* msg);
158 : using PUP::able::register_constructor;
159 : WRAPPED_PUPable_decl_template(BlastWave);
160 : /// \endcond
161 :
162 : /// @{
163 : /// Retrieve the GRMHD variables at a given position.
164 : template <typename DataType>
165 1 : auto variables(const tnsr::I<DataType, 3>& x,
166 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
167 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
168 :
169 : template <typename DataType>
170 1 : auto variables(const tnsr::I<DataType, 3>& x,
171 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
172 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
173 :
174 : template <typename DataType>
175 1 : auto variables(
176 : const tnsr::I<DataType, 3>& x,
177 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
178 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
179 :
180 : template <typename DataType>
181 1 : auto variables(const tnsr::I<DataType, 3>& x,
182 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
183 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
184 :
185 : template <typename DataType>
186 1 : auto variables(const tnsr::I<DataType, 3>& x,
187 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
188 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
189 :
190 : template <typename DataType>
191 1 : auto variables(const tnsr::I<DataType, 3>& x,
192 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
193 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
194 :
195 : template <typename DataType>
196 1 : auto variables(
197 : const tnsr::I<DataType, 3>& x,
198 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
199 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
200 :
201 : template <typename DataType>
202 1 : auto variables(const tnsr::I<DataType, 3>& x,
203 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
204 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
205 :
206 : template <typename DataType>
207 1 : auto variables(const tnsr::I<DataType, 3>& x,
208 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
209 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
210 :
211 : template <typename DataType>
212 1 : auto variables(const tnsr::I<DataType, 3>& x,
213 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
214 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
215 : return TemperatureInitialization::variables(
216 : x, tmpl::list<hydro::Tags::Temperature<DataType>>{});
217 : }
218 :
219 : /// @}
220 :
221 : /// Retrieve a collection of hydrodynamic variables at position x
222 : template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
223 1 : tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
224 : const tnsr::I<DataType, 3>& x,
225 : tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
226 : return {tuples::get<Tag1>(variables(x, tmpl::list<Tag1>{})),
227 : tuples::get<Tag2>(variables(x, tmpl::list<Tag2>{})),
228 : tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
229 : }
230 :
231 : /// Retrieve the metric variables
232 : template <typename DataType, typename Tag,
233 : Requires<tmpl::list_contains_v<
234 : gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
235 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
236 : tmpl::list<Tag> /*meta*/) const {
237 : constexpr double dummy_time = 0.0;
238 : return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
239 : }
240 :
241 0 : const EquationsOfState::IdealFluid<true>& equation_of_state() const {
242 : return equation_of_state_;
243 : }
244 :
245 : // NOLINTNEXTLINE(google-runtime-references)
246 0 : void pup(PUP::er& /*p*/) override;
247 :
248 : private:
249 0 : double inner_radius_ = std::numeric_limits<double>::signaling_NaN();
250 0 : double outer_radius_ = std::numeric_limits<double>::signaling_NaN();
251 0 : double inner_density_ = std::numeric_limits<double>::signaling_NaN();
252 0 : double outer_density_ = std::numeric_limits<double>::signaling_NaN();
253 0 : double inner_pressure_ = std::numeric_limits<double>::signaling_NaN();
254 0 : double outer_pressure_ = std::numeric_limits<double>::signaling_NaN();
255 0 : std::array<double, 3> magnetic_field_{
256 : {std::numeric_limits<double>::signaling_NaN(),
257 : std::numeric_limits<double>::signaling_NaN(),
258 : std::numeric_limits<double>::signaling_NaN()}};
259 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
260 0 : Geometry geometry_ = Geometry::Cylindrical;
261 0 : EquationsOfState::IdealFluid<true> equation_of_state_{};
262 0 : gr::Solutions::Minkowski<3> background_spacetime_{};
263 :
264 0 : friend bool operator==(const BlastWave& lhs, const BlastWave& rhs);
265 :
266 0 : friend bool operator!=(const BlastWave& lhs, const BlastWave& rhs);
267 : };
268 :
269 : } // namespace grmhd::AnalyticData
270 :
271 : /// \cond
272 : template <>
273 : struct Options::create_from_yaml<grmhd::AnalyticData::BlastWave::Geometry> {
274 : template <typename Metavariables>
275 : static grmhd::AnalyticData::BlastWave::Geometry create(
276 : const Options::Option& options) {
277 : return create<void>(options);
278 : }
279 : };
280 :
281 : template <>
282 : grmhd::AnalyticData::BlastWave::Geometry
283 : Options::create_from_yaml<grmhd::AnalyticData::BlastWave::Geometry>::create<
284 : void>(const Options::Option& options);
285 : /// \endcond
|
