SpECTRE Documentation Coverage Report
Current view: top level - PointwiseFunctions/AnalyticSolutions/NewtonianEuler - IsentropicVortex.hpp Hit Total Coverage
Commit: 6c78fa2058c3cb0e0713e24b0dcb28e30583111a Lines: 13 56 23.2 %
Date: 2024-03-01 16:00:07
Legend: Lines: hit not hit

          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           1 : 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

Generated by: LCOV version 1.14