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