KhInstability.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <limits>
8 
10 #include "Evolution/Systems/NewtonianEuler/Sources/NoSource.hpp"
11 #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
12 #include "Options/Options.hpp"
13 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
14 #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
15 #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
16 #include "Utilities/MakeArray.hpp"
17 #include "Utilities/TMPL.hpp"
18 #include "Utilities/TaggedTuple.hpp"
19 
20 /// \cond
21 namespace PUP {
22 class er; // IWYU pragma: keep
23 } // namespace PUP
24 /// \endcond
25 
26 namespace NewtonianEuler {
27 namespace AnalyticData {
28 
29 /*!
30  * \brief Initial data to simulate the Kelvin-Helmholtz instability.
31  *
32  * For comparison purposes, this class implements the planar shear of
33  * \cite Schaal2015, which is evolved using periodic boundary conditions.
34  * The initial state consists of a horizontal strip of mass
35  * density \f$\rho_\text{in}\f$ moving with horizontal speed
36  * \f$v_{\text{in}}\f$. The rest of the fluid possesses mass density
37  * \f$\rho_\text{out}\f$, and its horizontal velocity is \f$v_{\text{out}}\f$,
38  * both constant. Mathematically,
39  *
40  * \f{align*}
41  * \rho(x, y) =
42  * \begin{cases}
43  * \rho_\text{in}, & \left|y - y_\text{mid}\right| < b/2\\
44  * \rho_\text{out}, & \text{otherwise},
45  * \end{cases}
46  * \f}
47  *
48  * and
49  *
50  * \f{align*}
51  * v_x(x, y) =
52  * \begin{cases}
53  * v_{\text{in}}, & \left|y - y_\text{mid}\right| < b/2\\
54  * v_{\text{out}}, & \text{otherwise},
55  * \end{cases}
56  * \f}
57  *
58  * where \f$b > 0\f$ is the thickness of the strip, and \f$y = y_\text{mid}\f$
59  * is its horizontal bimedian. The initial pressure is set equal to a constant,
60  * and the system is evolved assuming an ideal fluid of known adiabatic index.
61  * Finally, in order to excite the instability, the vertical velocity is
62  * initialized to
63  *
64  * \f{align*}
65  * v_y(x, y) = A\sin(4\pi x)
66  * \left[\exp\left(-\dfrac{(y - y_\text{top})^2}{2\sigma^2}\right) +
67  * \exp\left(-\dfrac{(y - y_\text{bot})^2}{2\sigma^2}\right)\right],
68  * \f}
69  *
70  * whose net effect is to perturb the horizontal boundaries of the strip
71  * periodically along the \f$x-\f$axis. Here \f$A\f$ is the amplitude,
72  * \f$\sigma\f$ is a characteristic length for the perturbation width,
73  * and \f$y_\text{top} = y_\text{mid} + b/2\f$ and
74  * \f$y_\text{bot} = y_\text{mid} - b/2\f$ are the vertical coordinates
75  * of the top and bottom boundaries of the strip, respectively.
76  *
77  * \note The periodic form of the perturbation enforces the horizontal
78  * extent of the domain to be a multiple of 0.5. The domain chosen in
79  * \cite Schaal2015 is the square \f$[0, 1]^2\f$, which covers two wavelengths
80  * of the perturbation. In addition, the authors ran their simulations using
81  *
82  * ```
83  * AdiabaticIndex: 1.4
84  * StripBimedianHeight: 0.5
85  * StripThickness: 0.5
86  * StripDensity: 2.0
87  * StripVelocity: 0.5
88  * BackgroundDensity: 1.0
89  * BackgroundVelocity: -0.5
90  * Pressure: 2.5
91  * PerturbAmplitude: 0.1
92  * PerturbWidth: 0.035355339059327376 # 0.05/sqrt(2)
93  * ```
94  *
95  * \note This class can be used to initialize 2D and 3D data. In 3D, the strip
96  * is aligned to the \f$xy-\f$plane, and the vertical direction is taken to be
97  * the \f$z-\f$axis.
98  */
99 template <size_t Dim>
100 class KhInstability : public MarkAsAnalyticData {
101  public:
102  using equation_of_state_type = EquationsOfState::IdealFluid<false>;
103  using source_term_type = Sources::NoSource;
104 
105  /// The adiabatic index of the fluid.
106  struct AdiabaticIndex {
107  using type = double;
108  static constexpr OptionString help = {"The adiabatic index of the fluid."};
109  };
110 
111  /// The vertical coordinate of the horizontal bimedian of the strip.
113  using type = double;
114  static constexpr OptionString help = {"The height of the strip center."};
115  };
116 
117  /// The thickness of the strip.
118  struct StripThickness {
119  using type = double;
120  static type lower_bound() noexcept { return 0.0; }
121  static constexpr OptionString help = {
122  "The thickness of the horizontal strip."};
123  };
124 
125  /// The mass density in the strip
126  struct StripDensity {
127  using type = double;
128  static type lower_bound() noexcept { return 0.0; }
129  static constexpr OptionString help = {
130  "The mass density in the horizontal strip."};
131  };
132 
133  /// The velocity along \f$x\f$ in the strip
134  struct StripVelocity {
135  using type = double;
136  static constexpr OptionString help = {
137  "The velocity along x in the horizontal strip."};
138  };
139 
140  /// The mass density outside of the strip
142  using type = double;
143  static type lower_bound() noexcept { return 0.0; }
144  static constexpr OptionString help = {
145  "The mass density outside of the strip."};
146  };
147 
148  /// The velocity along \f$x\f$ outside of the strip
150  using type = double;
151  static constexpr OptionString help = {
152  "The velocity along x outside of the strip."};
153  };
154 
155  /// The initial (constant) pressure of the fluid
156  struct Pressure {
157  using type = double;
158  static type lower_bound() noexcept { return 0.0; }
159  static constexpr OptionString help = {"The initial (constant) pressure."};
160  };
161 
162  /// The amplitude of the perturbation
164  using type = double;
165  static constexpr OptionString help = {"The amplitude of the perturbation."};
166  };
167 
168  /// The characteristic length for the width of the perturbation
169  struct PerturbWidth {
170  using type = double;
171  static type lower_bound() noexcept { return 0.0; }
172  static constexpr OptionString help = {
173  "The characteristic length for the width of the perturbation."};
174  };
175 
176  using options =
180 
181  static constexpr OptionString help = {
182  "Initial data to simulate the KH instability."};
183 
184  KhInstability() = default;
185  KhInstability(const KhInstability& /*rhs*/) = delete;
186  KhInstability& operator=(const KhInstability& /*rhs*/) = delete;
187  KhInstability(KhInstability&& /*rhs*/) noexcept = default;
188  KhInstability& operator=(KhInstability&& /*rhs*/) noexcept = default;
189  ~KhInstability() = default;
190 
191  KhInstability(double adiabatic_index, double strip_bimedian_height,
192  double strip_thickness, double strip_density,
193  double strip_velocity, double background_density,
194  double background_velocity, double pressure,
195  double perturbation_amplitude, double perturbation_width);
196 
197  /// Retrieve a collection of hydrodynamic variables at position x
198  template <typename DataType, typename... Tags>
199  tuples::TaggedTuple<Tags...> variables(
200  const tnsr::I<DataType, Dim, Frame::Inertial>& x,
201  tmpl::list<Tags...> /*meta*/) const noexcept {
202  return {tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
203  }
204 
205  const EquationsOfState::IdealFluid<false>& equation_of_state() const
206  noexcept {
207  return equation_of_state_;
208  }
209 
210  // clang-tidy: no runtime references
211  void pup(PUP::er& /*p*/) noexcept; // NOLINT
212 
213  private:
214  // @{
215  /// Retrieve hydro variable at `x`
216  template <typename DataType>
217  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x,
218  tmpl::list<Tags::MassDensity<DataType>> /*meta*/
219  ) const noexcept
220  -> tuples::TaggedTuple<Tags::MassDensity<DataType>>;
221 
222  template <typename DataType>
223  auto variables(
224  const tnsr::I<DataType, Dim, Frame::Inertial>& x,
225  tmpl::list<Tags::Velocity<DataType, Dim, Frame::Inertial>> /*meta*/) const
226  noexcept
227  -> tuples::TaggedTuple<Tags::Velocity<DataType, Dim, Frame::Inertial>>;
228 
229  template <typename DataType>
230  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x,
231  tmpl::list<Tags::SpecificInternalEnergy<DataType>> /*meta*/
232  ) const noexcept
233  -> tuples::TaggedTuple<Tags::SpecificInternalEnergy<DataType>>;
234 
235  template <typename DataType>
236  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x,
237  tmpl::list<Tags::Pressure<DataType>> /*meta*/
238  ) const noexcept
239  -> tuples::TaggedTuple<Tags::Pressure<DataType>>;
240  // @}
241 
242  template <size_t SpatialDim>
243  friend bool
244  operator==( // NOLINT (clang-tidy: readability-redundant-declaration)
245  const KhInstability<SpatialDim>& lhs,
246  const KhInstability<SpatialDim>& rhs) noexcept;
247 
248  double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
249  double strip_bimedian_height_ = std::numeric_limits<double>::signaling_NaN();
250  double strip_half_thickness_ = std::numeric_limits<double>::signaling_NaN();
251  double strip_density_ = std::numeric_limits<double>::signaling_NaN();
252  double strip_velocity_ = std::numeric_limits<double>::signaling_NaN();
253  double background_density_ = std::numeric_limits<double>::signaling_NaN();
254  double background_velocity_ = std::numeric_limits<double>::signaling_NaN();
255  double pressure_ = std::numeric_limits<double>::signaling_NaN();
256  double perturbation_amplitude_ = std::numeric_limits<double>::signaling_NaN();
257  double perturbation_width_ = std::numeric_limits<double>::signaling_NaN();
258  EquationsOfState::IdealFluid<false> equation_of_state_{};
259 };
260 
261 template <size_t Dim>
262 bool operator!=(const KhInstability<Dim>& lhs,
263  const KhInstability<Dim>& rhs) noexcept;
264 
265 } // namespace AnalyticData
266 } // namespace NewtonianEuler
EquationsOfState
Contains all equations of state, including base class.
Definition: DarkEnergyFluid.hpp:26
NewtonianEuler::AnalyticData::KhInstability::StripDensity
The mass density in the strip.
Definition: KhInstability.hpp:126
std::rel_ops::operator!=
T operator!=(T... args)
Options.hpp
MakeArray.hpp
NewtonianEuler::Tags::MassDensity
The mass density of the fluid.
Definition: Tags.hpp:32
NewtonianEuler::AnalyticData::KhInstability::PerturbAmplitude
The amplitude of the perturbation.
Definition: KhInstability.hpp:163
NewtonianEuler
Items related to evolving the Newtonian Euler system.
Definition: EvolveNewtonianEulerFwd.hpp:8
NewtonianEuler::AnalyticData::KhInstability::BackgroundVelocity
The velocity along outside of the strip.
Definition: KhInstability.hpp:149
NewtonianEuler::AnalyticData::KhInstability::Pressure
The initial (constant) pressure of the fluid.
Definition: KhInstability.hpp:156
NewtonianEuler::AnalyticData::KhInstability::StripBimedianHeight
The vertical coordinate of the horizontal bimedian of the strip.
Definition: KhInstability.hpp:112
cstddef
EquationsOfState::IdealFluid< false >
NewtonianEuler::AnalyticData::KhInstability::AdiabaticIndex
The adiabatic index of the fluid.
Definition: KhInstability.hpp:106
NewtonianEuler::AnalyticData::KhInstability::variables
tuples::TaggedTuple< Tags... > variables(const tnsr::I< DataType, Dim, Frame::Inertial > &x, tmpl::list< Tags... >) const noexcept
Retrieve a collection of hydrodynamic variables at position x.
Definition: KhInstability.hpp:199
NewtonianEuler::AnalyticData::KhInstability::StripVelocity
The velocity along in the strip.
Definition: KhInstability.hpp:134
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
limits
TypeAliases.hpp
Frame
Definition: IndexType.hpp:36
NewtonianEuler::AnalyticData::KhInstability::BackgroundDensity
The mass density outside of the strip.
Definition: KhInstability.hpp:141
NewtonianEuler::AnalyticData::KhInstability::PerturbWidth
The characteristic length for the width of the perturbation.
Definition: KhInstability.hpp:169
NewtonianEuler::AnalyticData::KhInstability
Initial data to simulate the Kelvin-Helmholtz instability.
Definition: EvolveNewtonianEulerFwd.hpp:11
NewtonianEuler::AnalyticData::KhInstability::StripThickness
The thickness of the strip.
Definition: KhInstability.hpp:118
OptionString
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
TMPL.hpp