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/String.hpp"
11 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
12 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
14 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
16 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
17 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
18 : #include "Utilities/Serialization/CharmPupable.hpp"
19 : #include "Utilities/TMPL.hpp"
20 : #include "Utilities/TaggedTuple.hpp"
21 :
22 : /// \cond
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : namespace grmhd::AnalyticData {
29 :
30 : /*!
31 : * \brief Analytic initial data for a Kelvin-Helmholtz instability simulation.
32 : *
33 : * This is similar to the data from Section 4.7 of \cite Beckwith2011iy.
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 : * A uniform magnetic field can be added.
78 : */
79 1 : class KhInstability : public evolution::initial_data::InitialData,
80 : public MarkAsAnalyticData,
81 : public AnalyticDataBase,
82 : public hydro::TemperatureInitialization<KhInstability> {
83 : public:
84 0 : using equation_of_state_type = EquationsOfState::IdealFluid<true>;
85 :
86 : /// The adiabatic index of the fluid.
87 1 : struct AdiabaticIndex {
88 0 : using type = double;
89 0 : static constexpr Options::String help = {
90 : "The adiabatic index of the fluid."};
91 : };
92 :
93 : /// The vertical coordinate of the horizontal bimedian of the strip.
94 1 : struct StripBimedianHeight {
95 0 : using type = double;
96 0 : static constexpr Options::String help = {"The height of the strip center."};
97 : };
98 :
99 : /// The thickness of the strip.
100 1 : struct StripThickness {
101 0 : using type = double;
102 0 : static type lower_bound() { return 0.0; }
103 0 : static constexpr Options::String help = {
104 : "The thickness of the horizontal strip."};
105 : };
106 :
107 : /// The mass density in the strip
108 1 : struct StripDensity {
109 0 : using type = double;
110 0 : static type lower_bound() { return 0.0; }
111 0 : static constexpr Options::String help = {
112 : "The mass density in the horizontal strip."};
113 : };
114 :
115 : /// The velocity along \f$x\f$ in the strip
116 1 : struct StripVelocity {
117 0 : using type = double;
118 0 : static constexpr Options::String help = {
119 : "The velocity along x in the horizontal strip."};
120 : };
121 :
122 : /// The mass density outside of the strip
123 1 : struct BackgroundDensity {
124 0 : using type = double;
125 0 : static type lower_bound() { return 0.0; }
126 0 : static constexpr Options::String help = {
127 : "The mass density outside of the strip."};
128 : };
129 :
130 : /// The velocity along \f$x\f$ outside of the strip
131 1 : struct BackgroundVelocity {
132 0 : using type = double;
133 0 : static constexpr Options::String help = {
134 : "The velocity along x outside of the strip."};
135 : };
136 :
137 : /// The initial (constant) pressure of the fluid
138 1 : struct Pressure {
139 0 : using type = double;
140 0 : static type lower_bound() { return 0.0; }
141 0 : static constexpr Options::String help = {
142 : "The initial (constant) pressure."};
143 : };
144 :
145 : /// The amplitude of the perturbation
146 1 : struct PerturbAmplitude {
147 0 : using type = double;
148 0 : static constexpr Options::String help = {
149 : "The amplitude of the perturbation."};
150 : };
151 :
152 : /// The characteristic length for the width of the perturbation
153 1 : struct PerturbWidth {
154 0 : using type = double;
155 0 : static type lower_bound() { return 0.0; }
156 0 : static constexpr Options::String help = {
157 : "The characteristic length for the width of the perturbation."};
158 : };
159 :
160 : /// The uniform magnetic field
161 1 : struct MagneticField {
162 0 : using type = std::array<double, 3>;
163 0 : static constexpr Options::String help = {"The uniform magnetic field."};
164 : };
165 :
166 0 : using options = tmpl::list<AdiabaticIndex, StripBimedianHeight,
167 : StripThickness, StripDensity, StripVelocity,
168 : BackgroundDensity, BackgroundVelocity, Pressure,
169 : PerturbAmplitude, PerturbWidth, MagneticField>;
170 :
171 0 : static constexpr Options::String help = {
172 : "Initial data to simulate the magnetized KH instability."};
173 :
174 0 : KhInstability() = default;
175 0 : KhInstability(const KhInstability& /*rhs*/) = default;
176 0 : KhInstability& operator=(const KhInstability& /*rhs*/) = default;
177 0 : KhInstability(KhInstability&& /*rhs*/) = default;
178 0 : KhInstability& operator=(KhInstability&& /*rhs*/) = default;
179 0 : ~KhInstability() override = default;
180 :
181 0 : KhInstability(double adiabatic_index, double strip_bimedian_height,
182 : double strip_thickness, double strip_density,
183 : double strip_velocity, double background_density,
184 : double background_velocity, double pressure,
185 : double perturbation_amplitude, double perturbation_width,
186 : const std::array<double, 3>& magnetic_field);
187 :
188 0 : auto get_clone() const
189 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
190 :
191 : /// \cond
192 : explicit KhInstability(CkMigrateMessage* msg);
193 : using PUP::able::register_constructor;
194 : WRAPPED_PUPable_decl_template(KhInstability);
195 : /// \endcond
196 :
197 : /// @{
198 : /// Retrieve the GRMHD variables at a given position.
199 : template <typename DataType>
200 1 : auto variables(const tnsr::I<DataType, 3>& x,
201 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
202 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
203 :
204 : template <typename DataType>
205 1 : auto variables(const tnsr::I<DataType, 3>& x,
206 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
207 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
208 :
209 : template <typename DataType>
210 1 : auto variables(
211 : const tnsr::I<DataType, 3>& x,
212 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
213 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
214 :
215 : template <typename DataType>
216 1 : auto variables(const tnsr::I<DataType, 3>& x,
217 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
218 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
219 :
220 : template <typename DataType>
221 1 : auto variables(const tnsr::I<DataType, 3>& x,
222 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
223 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
224 :
225 : template <typename DataType>
226 1 : auto variables(const tnsr::I<DataType, 3>& x,
227 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
228 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
229 :
230 : template <typename DataType>
231 1 : auto variables(
232 : const tnsr::I<DataType, 3>& x,
233 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
234 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
235 :
236 : template <typename DataType>
237 1 : auto variables(const tnsr::I<DataType, 3>& x,
238 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
239 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
240 :
241 : template <typename DataType>
242 1 : auto variables(const tnsr::I<DataType, 3>& x,
243 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
244 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
245 :
246 : template <typename DataType>
247 1 : auto variables(const tnsr::I<DataType, 3>& x,
248 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
249 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
250 : return TemperatureInitialization::variables(
251 : x, tmpl::list<hydro::Tags::Temperature<DataType>>{});
252 : }
253 : /// @}
254 :
255 : /// Retrieve a collection of hydrodynamic variables at position x
256 : template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
257 1 : tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
258 : const tnsr::I<DataType, 3>& x,
259 : tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
260 : return {tuples::get<Tag1>(variables(x, tmpl::list<Tag1>{})),
261 : tuples::get<Tag2>(variables(x, tmpl::list<Tag2>{})),
262 : tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
263 : }
264 :
265 : /// Retrieve the metric variables
266 : template <typename DataType, typename Tag,
267 : Requires<tmpl::list_contains_v<
268 : gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
269 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
270 : tmpl::list<Tag> /*meta*/) const {
271 : constexpr double dummy_time = 0.0;
272 : return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
273 : }
274 :
275 0 : const EquationsOfState::IdealFluid<true>& equation_of_state() const {
276 : return equation_of_state_;
277 : }
278 :
279 : // NOLINTNEXTLINE(google-runtime-references)
280 0 : void pup(PUP::er& /*p*/) override;
281 :
282 : private:
283 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
284 0 : double strip_bimedian_height_ = std::numeric_limits<double>::signaling_NaN();
285 0 : double strip_half_thickness_ = std::numeric_limits<double>::signaling_NaN();
286 0 : double strip_density_ = std::numeric_limits<double>::signaling_NaN();
287 0 : double strip_velocity_ = std::numeric_limits<double>::signaling_NaN();
288 0 : double background_density_ = std::numeric_limits<double>::signaling_NaN();
289 0 : double background_velocity_ = std::numeric_limits<double>::signaling_NaN();
290 0 : double pressure_ = std::numeric_limits<double>::signaling_NaN();
291 0 : double perturbation_amplitude_ = std::numeric_limits<double>::signaling_NaN();
292 0 : double perturbation_width_ = std::numeric_limits<double>::signaling_NaN();
293 0 : std::array<double, 3> magnetic_field_{
294 : {std::numeric_limits<double>::signaling_NaN(),
295 : std::numeric_limits<double>::signaling_NaN(),
296 : std::numeric_limits<double>::signaling_NaN()}};
297 0 : EquationsOfState::IdealFluid<true> equation_of_state_{};
298 0 : gr::Solutions::Minkowski<3> background_spacetime_{};
299 :
300 0 : friend bool operator==(const KhInstability& lhs, const KhInstability& rhs);
301 :
302 0 : friend bool operator!=(const KhInstability& lhs, const KhInstability& rhs);
303 : };
304 : } // namespace grmhd::AnalyticData
|
