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 <limits>
       8             : #include <pup.h>
       9             : 
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Evolution/Systems/NewtonianEuler/Sources/NoSource.hpp"
      12             : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
      13             : #include "Options/String.hpp"
      14             : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
      15             : #include "PointwiseFunctions/AnalyticSolutions/Hydro/SmoothFlow.hpp"
      16             : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
      17             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      18             : #include "Utilities/MakeArray.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : #include "Utilities/TaggedTuple.hpp"
      21             : 
      22             : namespace NewtonianEuler::Solutions {
      23             : /*!
      24             :  * \brief Smooth density wave advecting across the domain.
      25             :  *
      26             :  * A solution with constant pressure and uniform spatial velocity provided
      27             :  * that the rest mass density satisfies the advection equation
      28             :  *
      29             :  * \f{align*}{
      30             :  * \partial_t\rho + v^i\partial_i\rho = 0,
      31             :  * \f}
      32             :  *
      33             :  * and the specific internal energy is a function of the rest mass density only,
      34             :  * \f$\epsilon = \epsilon(\rho)\f$. For testing purposes, this class implements
      35             :  * this solution for the case where \f$\rho\f$ is a sine wave. The user
      36             :  * specifies the mean flow velocity of the fluid, the wavevector of the density
      37             :  * profile, and the amplitude \f$A\f$ of the density profile. In Cartesian
      38             :  * coordinates \f$(x, y, z)\f$, and using dimensionless units, the primitive
      39             :  * variables at a given time \f$t\f$ are then
      40             :  *
      41             :  * \f{align*}{
      42             :  * \rho(\vec{x},t) &= 1 + A \sin(\vec{k}\cdot(\vec{x} - \vec{v}t)) \\
      43             :  * \vec{v}(\vec{x},t) &= [v_x, v_y, v_z]^{T},\\
      44             :  * P(\vec{x},t) &= P, \\
      45             :  * \epsilon(\vec{x}, t) &= \frac{P}{(\gamma - 1)\rho}\\
      46             :  * \f}
      47             :  *
      48             :  * where we have assumed \f$\epsilon\f$ and \f$\rho\f$ to be related through an
      49             :  * equation mathematically equivalent to the equation of state of an ideal gas,
      50             :  * where the pressure is held constant.
      51             :  */
      52             : template <size_t Dim>
      53           1 : class SmoothFlow : virtual public MarkAsAnalyticSolution,
      54             :                    private hydro::Solutions::SmoothFlow<Dim, false> {
      55           0 :   using smooth_flow = hydro::Solutions::SmoothFlow<Dim, false>;
      56             : 
      57             :  public:
      58           0 :   using options = typename smooth_flow::options;
      59             : 
      60           0 :   static constexpr Options::String help = {
      61             :       "Smooth density wave advecting across a domain."};
      62             : 
      63           0 :   SmoothFlow() = default;
      64           0 :   SmoothFlow(const SmoothFlow& /*rhs*/) = default;
      65           0 :   SmoothFlow& operator=(const SmoothFlow& /*rhs*/) = default;
      66           0 :   SmoothFlow(SmoothFlow&& /*rhs*/) = default;
      67           0 :   SmoothFlow& operator=(SmoothFlow&& /*rhs*/) = default;
      68           0 :   ~SmoothFlow() = default;
      69             : 
      70           0 :   SmoothFlow(const std::array<double, Dim>& mean_velocity,
      71             :              const std::array<double, Dim>& wavevector, double pressure,
      72             :              double adiabatic_index, double perturbation_size);
      73             : 
      74           0 :   explicit SmoothFlow(CkMigrateMessage* msg);
      75             : 
      76             :   using smooth_flow::equation_of_state;
      77             :   using typename smooth_flow::equation_of_state_type;
      78           0 :   using source_term_type = Sources::NoSource;
      79             : 
      80             :   // Overload the variables function from the base class.
      81             :   using smooth_flow::variables;
      82             : 
      83             :   template <typename DataType>
      84           0 :   tuples::TaggedTuple<Tags::MassDensity<DataType>> variables(
      85             :       const tnsr::I<DataType, Dim>& x, const double t,
      86             :       tmpl::list<Tags::MassDensity<DataType>> /*meta*/) const {
      87             :     return {tuples::get<hydro::Tags::RestMassDensity<DataType>>(
      88             :         variables(x, t, tmpl::list<hydro::Tags::RestMassDensity<DataType>>{}))};
      89             :   }
      90             : 
      91             :   template <typename DataType>
      92           0 :   tuples::TaggedTuple<Tags::Velocity<DataType, Dim>> variables(
      93             :       const tnsr::I<DataType, Dim>& x, const double t,
      94             :       tmpl::list<Tags::Velocity<DataType, Dim>> /*meta*/) const {
      95             :     return {tuples::get<hydro::Tags::SpatialVelocity<DataType, Dim>>(variables(
      96             :         x, t, tmpl::list<hydro::Tags::SpatialVelocity<DataType, Dim>>{}))};
      97             :   }
      98             : 
      99             :   template <typename DataType>
     100           0 :   tuples::TaggedTuple<Tags::Pressure<DataType>> variables(
     101             :       const tnsr::I<DataType, Dim>& x, const double t,
     102             :       tmpl::list<Tags::Pressure<DataType>> /*meta*/) const {
     103             :     return {tuples::get<hydro::Tags::Pressure<DataType>>(
     104             :         variables(x, t, tmpl::list<hydro::Tags::Pressure<DataType>>{}))};
     105             :   }
     106             : 
     107             :   template <typename DataType>
     108           0 :   tuples::TaggedTuple<Tags::SpecificInternalEnergy<DataType>> variables(
     109             :       const tnsr::I<DataType, Dim>& x, const double t,
     110             :       tmpl::list<Tags::SpecificInternalEnergy<DataType>> /*meta*/) const {
     111             :     return {
     112             :         tuples::get<hydro::Tags::SpecificInternalEnergy<DataType>>(variables(
     113             :             x, t,
     114             :             tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>>{}))};
     115             :   }
     116             : 
     117             :   /// Retrieve a collection of hydro variables at `(x, t)`
     118             :   template <typename DataType, typename... Tags>
     119           1 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, Dim>& x,
     120             :                                          const double t,
     121             :                                          tmpl::list<Tags...> /*meta*/) const {
     122             :     static_assert(sizeof...(Tags) > 1,
     123             :                   "The generic template will recurse infinitely if only one "
     124             :                   "tag is being retrieved.");
     125             :     return {tuples::get<Tags>(variables(x, t, tmpl::list<Tags>{}))...};
     126             :   }
     127             : 
     128             :   // clang-tidy: no runtime references
     129           0 :   void pup(PUP::er& /*p*/);  //  NOLINT
     130             : 
     131             :  private:
     132             :   template <size_t SpatialDim>
     133             :   friend bool
     134           0 :   operator==(  // NOLINT (clang-tidy: readability-redundant-declaration)
     135             :       const SmoothFlow<SpatialDim>& lhs, const SmoothFlow<SpatialDim>& rhs);
     136             : };
     137             : 
     138             : template <size_t Dim>
     139           0 : bool operator!=(const SmoothFlow<Dim>& lhs, const SmoothFlow<Dim>& rhs);
     140             : }  // namespace NewtonianEuler::Solutions