IsentropicVortex.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <limits>
8 #include <pup.h>
9 
11 #include "Evolution/Systems/NewtonianEuler/Tags.hpp" // IWYU pragma: keep
12 #include "Options/Options.hpp"
13 #include "Utilities/MakeArray.hpp" // IWYU pragma: keep
14 #include "Utilities/TMPL.hpp"
16 
17 namespace NewtonianEuler {
18 namespace Solutions {
19 
20 /*!
21  * \brief Newtonian isentropic vortex in Cartesian coordinates
22  *
23  * The analytic solution to the 2-D Newtonian Euler system
24  * representing the slow advection of an incompressible, isentropic
25  * vortex \ref vortex_ref "[1]". The initial condition is the superposition of a
26  * mean uniform flow with a gaussian-profile vortex. When embedded in
27  * 3-D space, the isentropic vortex is still a solution to the corresponding 3-D
28  * system if the velocity along the third axis is a constant. In Cartesian
29  * coordinates \f$(x, y, z)\f$, and using dimensionless units, the primitive
30  * quantities at a given time \f$t\f$ are then
31  *
32  * \f{align*}
33  * \rho &= \left[1 - \dfrac{(\gamma - 1)\beta^2}{8\gamma\pi^2}\exp\left(
34  * 1 - r^2\right)\right]^{1/(\gamma - 1)}, \\
35  * v_x &= U - \dfrac{\beta\tilde y}{2\pi}\exp\left(\dfrac{1 - r^2}{2}\right),\\
36  * v_y &= V + \dfrac{\beta\tilde x}{2\pi}\exp\left(\dfrac{1 - r^2}{2}\right),\\
37  * v_z &= W,\\
38  * \epsilon &= \frac{\rho^{\gamma - 1}}{\gamma - 1},
39  * \f}
40  *
41  * with
42  *
43  * \f{align*}
44  * r^2 &= {\tilde x}^2 + {\tilde y}^2,\\
45  * \tilde x &= x - X_0 - U t,\\
46  * \tilde y &= y - Y_0 - V t,
47  * \f}
48  *
49  * where \f$(X_0, Y_0)\f$ is the position of the vortex on the \f$(x, y)\f$
50  * plane at \f$t = 0\f$, \f$(U, V, W)\f$ are the components of the mean flow
51  * velocity, \f$\beta\f$ is the vortex strength, and \f$\gamma\f$ is the
52  * adiabatic index. The pressure \f$p\f$ is then obtained from the dimensionless
53  * polytropic relation
54  *
55  * \f{align*}
56  * p = \rho^\gamma.
57  * \f}
58  *
59  * On the other hand, if the velocity along the \f$z-\f$axis is not a constant
60  * but a function of the \f$z\f$ coordinate, the resulting modified isentropic
61  * vortex is still a solution to the Newtonian Euler system, but with source
62  * terms that are proportional to \f$dv_z/dz\f$. (See
63  * NewtonianEuler::Sources::IsentropicVortexSource.) For testing purposes,
64  * we choose to write the velocity as a uniform field plus a periodic
65  * perturbation,
66  *
67  * \f{align*}
68  * v_z(z) = W + \epsilon \sin{z},
69  * \f}
70  *
71  * where \f$\epsilon\f$ is the amplitude of the perturbation. The resulting
72  * source for the Newtonian Euler system will then be proportional to
73  * \f$\epsilon \cos{z}\f$.
74  *
75  * \anchor vortex_ref [1] H.C Yee, N.D Sandham, M.J Djomehri, Low-dissipative
76  * high-order shock-capturing methods using characteristic-based filters, J.
77  * Comput. Phys. [150 (1999) 199](http://dx.doi.org/10.1006/jcph.1998.6177)
78  */
80  public:
81  /// The adiabatic index of the fluid.
82  struct AdiabaticIndex {
83  using type = double;
84  static constexpr OptionString help = {"The adiabatic index of the fluid."};
85  // Note: bounds only valid for an ideal gas.
86  static type lower_bound() noexcept { return 1.0; }
87  static type upper_bound() noexcept { return 2.0; }
88  };
89 
90  /// The position of the center of the vortex at \f$t = 0\f$
91  struct Center {
93  static constexpr OptionString help = {
94  "The coordinates of the center of the vortex at t = 0."};
95  static type default_value() noexcept { return {{0.0, 0.0, 0.0}}; }
96  };
97 
98  /// The mean flow velocity.
99  struct MeanVelocity {
100  using type = std::array<double, 3>;
101  static constexpr OptionString help = {"The mean flow velocity."};
102  };
103 
104  /// The amplitude of the perturbation generating a source term.
106  using type = double;
107  static constexpr OptionString help = {
108  "The amplitude of the perturbation producing sources."};
109  };
110 
111  /// The strength of the vortex.
112  struct Strength {
113  using type = double;
114  static constexpr OptionString help = {"The strength of the vortex."};
115  static type lower_bound() noexcept { return 0.0; }
116  };
117 
118  using options = tmpl::list<AdiabaticIndex, Center, MeanVelocity,
120  static constexpr OptionString help = {"Newtonian Isentropic Vortex."};
121 
122  IsentropicVortex() = default;
123  IsentropicVortex(const IsentropicVortex& /*rhs*/) = delete;
124  IsentropicVortex& operator=(const IsentropicVortex& /*rhs*/) = delete;
125  IsentropicVortex(IsentropicVortex&& /*rhs*/) noexcept = default;
126  IsentropicVortex& operator=(IsentropicVortex&& /*rhs*/) noexcept = default;
127  ~IsentropicVortex() = default;
128 
129  IsentropicVortex(double adiabatic_index, Center::type center,
130  MeanVelocity::type mean_velocity,
131  double perturbation_amplitude, double strength);
132 
133  explicit IsentropicVortex(CkMigrateMessage* /*unused*/) noexcept {}
134 
135  template <typename DataType>
136  using primitive_t =
137  tmpl::list<Tags::MassDensity<DataType>, Tags::Velocity<DataType, 3>,
139 
140  template <typename DataType>
141  using conservative_t = tmpl::list<Tags::MassDensity<DataType>,
144 
145  template <typename DataType>
146  Scalar<DataType> perturbation(const DataType& coord_z) const noexcept;
147 
148  template <typename DataType>
149  tuples::tagged_tuple_from_typelist<primitive_t<DataType>> primitive_variables(
150  const tnsr::I<DataType, 3>& x, double t) const noexcept;
151 
152  template <typename DataType>
153  tuples::tagged_tuple_from_typelist<conservative_t<DataType>>
154  conservative_variables(const tnsr::I<DataType, 3>& x, double t) const
155  noexcept;
156 
157  // clang-tidy: no runtime references
158  void pup(PUP::er& /*p*/) noexcept; // NOLINT
159 
160  constexpr double adiabatic_index() const noexcept { return adiabatic_index_; }
161  constexpr const Center::type& center() const noexcept { return center_; }
162  constexpr const MeanVelocity::type& mean_velocity() const noexcept {
163  return mean_velocity_;
164  }
165  constexpr double perturbation_amplitude() const noexcept {
166  return perturbation_amplitude_;
167  }
168  constexpr double strength() const noexcept { return strength_; }
169 
170  private:
171  double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
172  Center::type center_ = {{0.0, 0.0, 0.0}};
173  MeanVelocity::type mean_velocity_ =
175  double perturbation_amplitude_ = std::numeric_limits<double>::signaling_NaN();
176  double strength_ = std::numeric_limits<double>::signaling_NaN();
177 };
178 
179 inline constexpr bool operator==(const IsentropicVortex& lhs,
180  const IsentropicVortex& rhs) noexcept {
181  return lhs.adiabatic_index() == rhs.adiabatic_index() and
182  lhs.center() == rhs.center() and
183  lhs.mean_velocity() == rhs.mean_velocity() and
184  lhs.perturbation_amplitude() == rhs.perturbation_amplitude() and
185  lhs.strength() == rhs.strength();
186 }
187 
188 inline constexpr bool operator!=(const IsentropicVortex& lhs,
189  const IsentropicVortex& rhs) noexcept {
190  return not(lhs == rhs);
191 }
192 
193 } // namespace Solutions
194 } // namespace NewtonianEuler
The momentum density of the fluid.
Definition: Tags.hpp:24
The strength of the vortex.
Definition: IsentropicVortex.hpp:112
Defines class tuples::TaggedTuple.
T signaling_NaN(T... args)
The position of the center of the vortex at .
Definition: IsentropicVortex.hpp:91
Defines function make_array.
Newtonian isentropic vortex in Cartesian coordinates.
Definition: IsentropicVortex.hpp:79
The amplitude of the perturbation generating a source term.
Definition: IsentropicVortex.hpp:105
The macroscopic or flow velocity of the fluid.
Definition: Tags.hpp:38
Defines classes and functions for making classes creatable from input files.
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
The specific internal energy of the fluid.
Definition: Tags.hpp:45
The energy density of the fluid.
Definition: Tags.hpp:31
Defines a list of useful type aliases for tensors.
Wraps the template metaprogramming library used (brigand)
The adiabatic index of the fluid.
Definition: IsentropicVortex.hpp:82
The mean flow velocity.
Definition: IsentropicVortex.hpp:99
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Items related to evolving the Newtonian Euler system.
Definition: Characteristics.hpp:17