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/Tensor.hpp"
      12             : #include "DataStructures/Tensor/TypeAliases.hpp"
      13             : #include "Options/Context.hpp"
      14             : #include "Options/String.hpp"
      15             : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
      16             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      17             : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
      18             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      19             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      20             : #include "Utilities/MakeArray.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : #include "Utilities/TaggedTuple.hpp"
      23             : 
      24             : /// \cond
      25             : namespace PUP {
      26             : class er;  // IWYU pragma: keep
      27             : }  // namespace PUP
      28             : /// \endcond
      29             : 
      30             : namespace NewtonianEuler::AnalyticData {
      31             : /*!
      32             :  * \brief A cylindrical or spherical Sod explosion \cite Toro2009 \cite Sod19781
      33             :  *
      34             :  * Common initial conditions are:
      35             :  *
      36             :  * \f{align*}{
      37             :  * (\rho, v^i, p) =
      38             :  * &\left\{
      39             :  * \begin{array}{ll}
      40             :  *   (1 ,0, 1) & \mathrm{if} \; r \le 0.5 \\
      41             :  *   (0.125 ,0, 0.1) & \mathrm{if} \; r > 0.5
      42             :  * \end{array}\right.
      43             :  * \f}
      44             :  *
      45             :  * where \f$r\f$ is the cylindrical (2d) or spherical (3d) radius. This test
      46             :  * problem uses an adiabatic index of 1.4. A reference solution can be computed
      47             :  * in 1d by solving the Newtonian Euler equations in cylindrical or spherical
      48             :  * symmetry. Note that the inner and outer density and pressure, as well as
      49             :  * where the initial discontinuity is can be chosen arbitrarily.
      50             :  */
      51             : template <size_t Dim>
      52           1 : class SodExplosion : public evolution::initial_data::InitialData,
      53             :                      public MarkAsAnalyticData {
      54             :  public:
      55             :   static_assert(Dim > 1, "Sod explosion is a 2d and 3d problem.");
      56           0 :   using equation_of_state_type = EquationsOfState::IdealFluid<false>;
      57             : 
      58             :   /// Initial radius of the discontinuity
      59           1 :   struct InitialRadius {
      60           0 :     using type = double;
      61           0 :     static constexpr Options::String help = {
      62             :         "The initial radius of the discontinuity."};
      63           0 :     static type lower_bound() { return 0.0; }
      64             :   };
      65             : 
      66           0 :   struct InnerMassDensity {
      67           0 :     using type = double;
      68           0 :     static constexpr Options::String help = {"The inner mass density."};
      69           0 :     static type lower_bound() { return 0.0; }
      70             :   };
      71             : 
      72           0 :   struct InnerPressure {
      73           0 :     using type = double;
      74           0 :     static constexpr Options::String help = {"The inner pressure."};
      75           0 :     static type lower_bound() { return 0.0; }
      76             :   };
      77             : 
      78           0 :   struct OuterMassDensity {
      79           0 :     using type = double;
      80           0 :     static constexpr Options::String help = {"The outer mass density."};
      81           0 :     static type lower_bound() { return 0.0; }
      82             :   };
      83             : 
      84           0 :   struct OuterPressure {
      85           0 :     using type = double;
      86           0 :     static constexpr Options::String help = {"The outer pressure."};
      87           0 :     static type lower_bound() { return 0.0; }
      88             :   };
      89             : 
      90           0 :   using options = tmpl::list<InitialRadius, InnerMassDensity, InnerPressure,
      91             :                              OuterMassDensity, OuterPressure>;
      92             : 
      93           0 :   static constexpr Options::String help = {
      94             :       "Cylindrical or spherical Sod explosion."};
      95             : 
      96           0 :   SodExplosion() = default;
      97           0 :   SodExplosion(const SodExplosion& /*rhs*/) = default;
      98           0 :   SodExplosion& operator=(const SodExplosion& /*rhs*/) = default;
      99           0 :   SodExplosion(SodExplosion&& /*rhs*/) = default;
     100           0 :   SodExplosion& operator=(SodExplosion&& /*rhs*/) = default;
     101           0 :   ~SodExplosion() override = default;
     102             : 
     103           0 :   auto get_clone() const
     104             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     105             : 
     106             :   /// \cond
     107             :   explicit SodExplosion(CkMigrateMessage* msg);
     108             :   using PUP::able::register_constructor;
     109             :   WRAPPED_PUPable_decl_template(SodExplosion);
     110             :   /// \endcond
     111             : 
     112           0 :   SodExplosion(double initial_radius, double inner_mass_density,
     113             :                double inner_pressure, double outer_mass_density,
     114             :                double outer_pressure, const Options::Context& context = {});
     115             : 
     116             :   /// Retrieve a collection of hydrodynamic variables at position x
     117             :   template <typename... Tags>
     118           1 :   tuples::TaggedTuple<Tags...> variables(
     119             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
     120             :       tmpl::list<Tags...> /*meta*/) const {
     121             :     return {tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
     122             :   }
     123             : 
     124           0 :   const equation_of_state_type& equation_of_state() const {
     125             :     return equation_of_state_;
     126             :   }
     127             : 
     128             :   // NOLINTNEXTLINE(google-runtime-references)
     129           0 :   void pup(PUP::er& /*p*/) override;
     130             : 
     131             :  private:
     132           0 :   tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataVector>> variables(
     133             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
     134             :       tmpl::list<hydro::Tags::RestMassDensity<DataVector>> /*meta*/) const;
     135             : 
     136             :   tuples::TaggedTuple<
     137             :       hydro::Tags::SpatialVelocity<DataVector, Dim, Frame::Inertial>>
     138           0 :   variables(const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
     139             :             tmpl::list<hydro::Tags::SpatialVelocity<
     140             :                 DataVector, Dim, Frame::Inertial>> /*meta*/) const;
     141             : 
     142           0 :   tuples::TaggedTuple<hydro::Tags::Pressure<DataVector>> variables(
     143             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
     144             :       tmpl::list<hydro::Tags::Pressure<DataVector>> /*meta*/) const;
     145             : 
     146             :   tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataVector>>
     147           0 :   variables(
     148             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
     149             :       tmpl::list<hydro::Tags::SpecificInternalEnergy<DataVector>> /*meta*/)
     150             :       const;
     151             : 
     152             :   template <size_t SpatialDim>
     153             :   friend bool
     154           0 :   operator==(  // NOLINT (clang-tidy: readability-redundant-declaration)
     155             :       const SodExplosion<SpatialDim>& lhs, const SodExplosion<SpatialDim>& rhs);
     156             : 
     157           0 :   double initial_radius_ = std::numeric_limits<double>::signaling_NaN();
     158           0 :   double inner_mass_density_ = std::numeric_limits<double>::signaling_NaN();
     159           0 :   double inner_pressure_ = std::numeric_limits<double>::signaling_NaN();
     160           0 :   double outer_mass_density_ = std::numeric_limits<double>::signaling_NaN();
     161           0 :   double outer_pressure_ = std::numeric_limits<double>::signaling_NaN();
     162           0 :   equation_of_state_type equation_of_state_{};
     163             : };
     164             : 
     165             : template <size_t Dim>
     166           0 : bool operator!=(const SodExplosion<Dim>& lhs, const SodExplosion<Dim>& rhs);
     167             : }  // namespace NewtonianEuler::AnalyticData