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