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