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 <cstddef>
8 : #include <limits>
9 :
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Options/String.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
14 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
16 : #include "PointwiseFunctions/Hydro/Tags.hpp"
17 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
18 : #include "Utilities/MakeArray.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 : // IWYU pragma: no_forward_declare NewtonianEuler::Solutions::IsentropicVortex::IntermediateVariables
29 :
30 1 : namespace NewtonianEuler::Solutions {
31 :
32 : /*!
33 : * \brief Newtonian isentropic vortex in Cartesian coordinates
34 : *
35 : * The analytic solution to the 2-D Newtonian Euler system
36 : * representing the slow advection of an incompressible, isentropic
37 : * vortex \cite Yee1999. The initial condition is the superposition of a
38 : * mean uniform flow with a gaussian-profile vortex. When embedded in
39 : * 3-D space, the isentropic vortex is still a solution to the corresponding 3-D
40 : * system if the velocity along the third axis is a constant. In Cartesian
41 : * coordinates \f$(x, y, z)\f$, and using dimensionless units, the primitive
42 : * quantities at a given time \f$t\f$ are then
43 : *
44 : * \f{align*}
45 : * \rho &= \left[1 - \dfrac{(\gamma - 1)\beta^2}{8\gamma\pi^2}\exp\left(
46 : * 1 - r^2\right)\right]^{1/(\gamma - 1)}, \\
47 : * v_x &= U - \dfrac{\beta\tilde y}{2\pi}\exp\left(\dfrac{1 - r^2}{2}\right),\\
48 : * v_y &= V + \dfrac{\beta\tilde x}{2\pi}\exp\left(\dfrac{1 - r^2}{2}\right),\\
49 : * v_z &= W,\\
50 : * \epsilon &= \frac{\rho^{\gamma - 1}}{\gamma - 1},
51 : * \f}
52 : *
53 : * with
54 : *
55 : * \f{align*}
56 : * r^2 &= {\tilde x}^2 + {\tilde y}^2,\\
57 : * \tilde x &= x - X_0 - U t,\\
58 : * \tilde y &= y - Y_0 - V t,
59 : * \f}
60 : *
61 : * where \f$(X_0, Y_0)\f$ is the position of the vortex on the \f$(x, y)\f$
62 : * plane at \f$t = 0\f$, \f$(U, V, W)\f$ are the components of the mean flow
63 : * velocity, \f$\beta\f$ is the vortex strength, and \f$\gamma\f$ is the
64 : * adiabatic index. The pressure \f$p\f$ is then obtained from the dimensionless
65 : * polytropic relation
66 : *
67 : * \f{align*}
68 : * p = \rho^\gamma.
69 : * \f}
70 : *
71 : * On the other hand, if the velocity along the \f$z-\f$axis is not a constant
72 : * but a function of the \f$z\f$ coordinate, the resulting modified isentropic
73 : * vortex is still a solution to the Newtonian Euler system, but with source
74 : * terms that are proportional to \f$dv_z/dz\f$. (See
75 : * NewtonianEuler::Sources::VortexPerturbation.) For testing purposes,
76 : * we choose to write the velocity as a uniform field plus a periodic
77 : * perturbation,
78 : *
79 : * \f{align*}
80 : * v_z(z) = W + \epsilon \sin{z},
81 : * \f}
82 : *
83 : * where \f$\epsilon\f$ is the amplitude of the perturbation. The resulting
84 : * source for the Newtonian Euler system will then be proportional to
85 : * \f$\epsilon \cos{z}\f$.
86 : */
87 : template <size_t Dim>
88 1 : class IsentropicVortex : public evolution::initial_data::InitialData,
89 : public MarkAsAnalyticSolution {
90 : static_assert(Dim == 2 or Dim == 3,
91 : "IsentropicVortex solution works in 2 and 3 dimensions");
92 :
93 : template <typename DataType>
94 : struct IntermediateVariables;
95 :
96 : public:
97 0 : using equation_of_state_type = EquationsOfState::PolytropicFluid<false>;
98 :
99 : /// The adiabatic index of the fluid.
100 1 : struct AdiabaticIndex {
101 0 : using type = double;
102 0 : static constexpr Options::String help = {
103 : "The adiabatic index of the fluid."};
104 : };
105 :
106 : /// The position of the center of the vortex at \f$t = 0\f$
107 1 : struct Center {
108 0 : using type = std::array<double, Dim>;
109 0 : static constexpr Options::String help = {
110 : "The coordinates of the center of the vortex at t = 0."};
111 : };
112 :
113 : /// The mean flow velocity.
114 1 : struct MeanVelocity {
115 0 : using type = std::array<double, Dim>;
116 0 : static constexpr Options::String help = {"The mean flow velocity."};
117 : };
118 :
119 : /// The amplitude of the perturbation generating a source term.
120 1 : struct PerturbAmplitude {
121 0 : using type = double;
122 0 : static constexpr Options::String help = {
123 : "The amplitude of the perturbation producing sources."};
124 : };
125 :
126 : /// The strength of the vortex.
127 1 : struct Strength {
128 0 : using type = double;
129 0 : static constexpr Options::String help = {"The strength of the vortex."};
130 0 : static type lower_bound() { return 0.0; }
131 : };
132 :
133 0 : using options = tmpl::conditional_t<
134 : Dim == 3,
135 : tmpl::list<AdiabaticIndex, Center, MeanVelocity, Strength,
136 : PerturbAmplitude>,
137 : tmpl::list<AdiabaticIndex, Center, MeanVelocity, Strength>>;
138 :
139 0 : static constexpr Options::String help = {
140 : "Newtonian Isentropic Vortex. Works in 2 and 3 dimensions."};
141 :
142 0 : IsentropicVortex() = default;
143 0 : IsentropicVortex(const IsentropicVortex& /*rhs*/) = default;
144 0 : IsentropicVortex& operator=(const IsentropicVortex& /*rhs*/) = default;
145 0 : IsentropicVortex(IsentropicVortex&& /*rhs*/) = default;
146 0 : IsentropicVortex& operator=(IsentropicVortex&& /*rhs*/) = default;
147 0 : ~IsentropicVortex() override = default;
148 :
149 0 : auto get_clone() const
150 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
151 :
152 : /// \cond
153 : explicit IsentropicVortex(CkMigrateMessage* msg);
154 : using PUP::able::register_constructor;
155 : WRAPPED_PUPable_decl_template(IsentropicVortex);
156 : /// \endcond
157 :
158 0 : IsentropicVortex(double adiabatic_index,
159 : const std::array<double, Dim>& center,
160 : const std::array<double, Dim>& mean_velocity,
161 : double strength, double perturbation_amplitude = 0.0);
162 :
163 : /// Retrieve a collection of hydrodynamic variables at position x and time t
164 : template <typename DataType, typename... Tags>
165 1 : tuples::TaggedTuple<Tags...> variables(
166 : const tnsr::I<DataType, Dim, Frame::Inertial>& x,
167 : const double t, // NOLINT
168 : tmpl::list<Tags...> /*meta*/) const {
169 : static_assert(sizeof...(Tags) > 1,
170 : "The generic template will recurse infinitely if only one "
171 : "tag is being retrieved.");
172 : const IntermediateVariables<DataType> vars(x, t, center_, mean_velocity_,
173 : strength_);
174 : return {tuples::get<Tags>(variables(tmpl::list<Tags>{}, vars))...};
175 : }
176 :
177 : /// Function of `z` coordinate to compute the perturbation generating
178 : /// a source term. Public so the corresponding source class can also use it.
179 : template <typename DataType>
180 1 : DataType perturbation_profile(const DataType& z) const;
181 :
182 : template <typename DataType>
183 0 : DataType deriv_of_perturbation_profile(const DataType& z) const;
184 :
185 : // To be used by VortexPerturbation source term
186 0 : double perturbation_amplitude() const { return perturbation_amplitude_; }
187 :
188 0 : const EquationsOfState::PolytropicFluid<false>& equation_of_state() const {
189 : return equation_of_state_;
190 : }
191 :
192 : // NOLINTNEXTLINE(google-runtime-references)
193 0 : void pup(PUP::er& /*p*/) override;
194 :
195 : private:
196 : /// @{
197 : /// Retrieve hydro variable at `(x, t)`
198 : template <typename DataType>
199 1 : auto variables(tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/,
200 : const IntermediateVariables<DataType>& vars) const
201 : -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
202 :
203 : template <typename DataType>
204 1 : auto variables(
205 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, Dim>> /*meta*/,
206 : const IntermediateVariables<DataType>& vars) const
207 : -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, Dim>>;
208 :
209 : template <typename DataType>
210 1 : auto variables(
211 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/,
212 : const IntermediateVariables<DataType>& vars) const
213 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
214 :
215 : template <typename DataType>
216 1 : auto variables(tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
217 : const IntermediateVariables<DataType>& vars) const
218 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
219 : /// @}
220 :
221 : // Intermediate variables needed to compute the primitives
222 : template <typename DataType>
223 0 : struct IntermediateVariables {
224 0 : IntermediateVariables(const tnsr::I<DataType, Dim, Frame::Inertial>& x,
225 : double t, const std::array<double, Dim>& center,
226 : const std::array<double, Dim>& mean_velocity,
227 : double strength);
228 0 : DataType x_tilde{};
229 0 : DataType y_tilde{};
230 0 : DataType profile{};
231 : // (3D only) z-coordinate to compute perturbation term
232 0 : DataType z_coord{};
233 : };
234 :
235 : template <size_t SpatialDim>
236 : friend bool
237 0 : operator==( // NOLINT (clang-tidy: readability-redundant-declaration)
238 : const IsentropicVortex<SpatialDim>& lhs,
239 : const IsentropicVortex<SpatialDim>& rhs);
240 :
241 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
242 0 : std::array<double, Dim> center_ =
243 : make_array<Dim>(std::numeric_limits<double>::signaling_NaN());
244 0 : std::array<double, Dim> mean_velocity_ =
245 : make_array<Dim>(std::numeric_limits<double>::signaling_NaN());
246 0 : double perturbation_amplitude_ = 0.0;
247 0 : double strength_ = std::numeric_limits<double>::signaling_NaN();
248 :
249 : // This is an ideal gas undergoing an isentropic process,
250 : // so the relation between the pressure and the mass density is polytropic,
251 : // where the polytropic exponent corresponds to the adiabatic index.
252 0 : EquationsOfState::PolytropicFluid<false> equation_of_state_{};
253 : };
254 :
255 : template <size_t Dim>
256 0 : bool operator!=(const IsentropicVortex<Dim>& lhs,
257 : const IsentropicVortex<Dim>& rhs);
258 :
259 : } // namespace NewtonianEuler::Solutions
|
