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