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