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 <memory>
       9             : 
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Options/String.hpp"
      12             : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
      13             : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
      14             : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
      15             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      16             : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
      17             : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
      18             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      19             : #include "Utilities/Serialization/CharmPupable.hpp"
      20             : #include "Utilities/TMPL.hpp"
      21             : #include "Utilities/TaggedTuple.hpp"
      22             : 
      23             : /// \cond
      24             : namespace PUP {
      25             : class er;
      26             : }  // namespace PUP
      27             : /// \endcond
      28             : 
      29             : namespace grmhd::AnalyticData {
      30             : 
      31             : /*!
      32             :  * \brief Analytic initial data for a slab jet
      33             :  *
      34             :  * This test problem is described in \cite Komissarov1999, Section 7.4 and
      35             :  * Fig. 13. It involves a high Lorentz factor jet injected into an ambient fluid
      36             :  * with an initially uniform magnetic field.
      37             :  *
      38             :  * The parameters used in \cite Komissarov1999 are:
      39             :  *
      40             :  * ```yaml
      41             :  * AdiabaticIndex: 4. / 3.
      42             :  * AmbientDensity: 10.
      43             :  * AmbientPressure: 0.01
      44             :  * JetDensity: 0.1
      45             :  * JetPressure: 0.01
      46             :  * # u^i = [20, 0, 0] or W = sqrt(401)
      47             :  * JetVelocity: [0.9987523388778445, 0., 0.]
      48             :  * InletRadius: 1.
      49             :  * MagneticField: [1., 0., 0.]
      50             :  * ```
      51             :  *
      52             :  * In \cite Komissarov1999 an artificial dissipation of $\eta_u=0.2$ and
      53             :  * $\eta_b=0.15$ is used, which we don't use here (yet).
      54             :  *
      55             :  * \note The inlet is currently modeled as the region of the domain where
      56             :  * $x <= 0$ and $|y| <= \mathrm{inlet_radius}$. Since this class only sets
      57             :  * the initial data, no fluid will flow into the domain during the evolution.
      58             :  * This has to be modeled as a boundary condition and is not implemented yet.
      59             :  */
      60           1 : class SlabJet : public evolution::initial_data::InitialData,
      61             :                 public MarkAsAnalyticData,
      62             :                 public AnalyticDataBase,
      63             :                 public hydro::TemperatureInitialization<SlabJet> {
      64             :  public:
      65           0 :   using equation_of_state_type = EquationsOfState::IdealFluid<true>;
      66             : 
      67           0 :   struct AdiabaticIndex {
      68           0 :     using type = double;
      69           0 :     static constexpr Options::String help = {
      70             :         "The adiabatic index of the ideal fluid"};
      71           0 :     static double lower_bound() { return 1.; }
      72             :   };
      73           0 :   struct AmbientDensity {
      74           0 :     using type = double;
      75           0 :     static constexpr Options::String help = {
      76             :         "Fluid rest mass density outside the jet"};
      77           0 :     static double lower_bound() { return 0.; }
      78           0 :     static double suggested_value() { return 10.; }
      79             :   };
      80           0 :   struct AmbientPressure {
      81           0 :     using type = double;
      82           0 :     static constexpr Options::String help = {"Fluid pressure outside the jet"};
      83           0 :     static double lower_bound() { return 0.; }
      84           0 :     static double suggested_value() { return 0.01; }
      85             :   };
      86           0 :   struct AmbientElectronFraction {
      87           0 :     using type = double;
      88           0 :     static constexpr Options::String help = {
      89             :         "Electron fraction outside the jet"};
      90           0 :     static double lower_bound() { return 0.; }
      91           0 :     static double upper_bound() { return 1.; }
      92             :   };
      93           0 :   struct JetDensity {
      94           0 :     using type = double;
      95           0 :     static constexpr Options::String help = {
      96             :         "Fluid rest mass density of the jet inlet"};
      97           0 :     static double lower_bound() { return 0.; }
      98           0 :     static double suggested_value() { return 0.1; }
      99             :   };
     100           0 :   struct JetPressure {
     101           0 :     using type = double;
     102           0 :     static constexpr Options::String help = {"Fluid pressure of the jet inlet"};
     103           0 :     static double lower_bound() { return 0.; }
     104           0 :     static double suggested_value() { return 0.01; }
     105             :   };
     106           0 :   struct JetElectronFraction {
     107           0 :     using type = double;
     108           0 :     static constexpr Options::String help = {
     109             :         "Electron fraction of the jet inlet"};
     110           0 :     static double lower_bound() { return 0.; }
     111           0 :     static double upper_bound() { return 1.; }
     112             :   };
     113           0 :   struct JetVelocity {
     114           0 :     using type = std::array<double, 3>;
     115           0 :     static constexpr Options::String help = {
     116             :         "Fluid spatial velocity of the jet inlet"};
     117             :   };
     118           0 :   struct InletRadius {
     119           0 :     using type = double;
     120           0 :     static constexpr Options::String help = {
     121             :         "Radius of the jet inlet around y=0"};
     122           0 :     static double lower_bound() { return 0.; }
     123           0 :     static double suggested_value() { return 1.; }
     124             :   };
     125           0 :   struct MagneticField {
     126           0 :     using type = std::array<double, 3>;
     127           0 :     static constexpr Options::String help = {
     128             :         "Initially uniform magnetic field"};
     129           0 :     static std::array<double, 3> suggested_value() { return {{1., 0., 0.}}; }
     130             :   };
     131             : 
     132           0 :   using options =
     133             :       tmpl::list<AdiabaticIndex, AmbientDensity, AmbientPressure,
     134             :                  AmbientElectronFraction, JetDensity, JetPressure,
     135             :                  JetElectronFraction, JetVelocity, InletRadius, MagneticField>;
     136             : 
     137           0 :   static constexpr Options::String help = {
     138             :       "Analytic initial data for a jet test."};
     139             : 
     140           0 :   SlabJet() = default;
     141           0 :   SlabJet(const SlabJet& /*rhs*/) = default;
     142           0 :   SlabJet& operator=(const SlabJet& /*rhs*/) = default;
     143           0 :   SlabJet(SlabJet&& /*rhs*/) = default;
     144           0 :   SlabJet& operator=(SlabJet&& /*rhs*/) = default;
     145           0 :   ~SlabJet() = default;
     146             : 
     147           0 :   SlabJet(double adiabatic_index, double ambient_density,
     148             :           double ambient_pressure, double ambient_electron_fraction,
     149             :           double jet_density, double jet_pressure, double jet_electron_fraction,
     150             :           std::array<double, 3> jet_velocity, double inlet_radius,
     151             :           std::array<double, 3> magnetic_field);
     152             : 
     153           0 :   auto get_clone() const
     154             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     155             : 
     156             :   /// \cond
     157             :   explicit SlabJet(CkMigrateMessage* msg);
     158             :   using PUP::able::register_constructor;
     159             :   WRAPPED_PUPable_decl_template(SlabJet);
     160             :   /// \endcond
     161             : 
     162             :   /// @{
     163             :   /// Retrieve the GRMHD variables at a given position.
     164             :   template <typename DataType>
     165           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     166             :                  tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
     167             :       const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
     168             : 
     169             :   template <typename DataType>
     170           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     171             :                  tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
     172             :       const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
     173             : 
     174             :   template <typename DataType>
     175           1 :   auto variables(
     176             :       const tnsr::I<DataType, 3>& x,
     177             :       tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
     178             :       -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
     179             : 
     180             :   template <typename DataType>
     181           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     182             :                  tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
     183             :       -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
     184             : 
     185             :   template <typename DataType>
     186           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     187             :                  tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
     188             :       const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
     189             : 
     190             :   template <typename DataType>
     191           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     192             :                  tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
     193             :       const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
     194             : 
     195             :   template <typename DataType>
     196           1 :   auto variables(
     197             :       const tnsr::I<DataType, 3>& x,
     198             :       tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
     199             :       -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
     200             : 
     201             :   template <typename DataType>
     202           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     203             :                  tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
     204             :       const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
     205             : 
     206             :   template <typename DataType>
     207           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     208             :                  tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
     209             :       const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
     210             : 
     211             :   template <typename DataType>
     212           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     213             :                  tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
     214             :       -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
     215             :     return TemperatureInitialization::variables(
     216             :         x, tmpl::list<hydro::Tags::Temperature<DataType>>{});
     217             :   }
     218             :   /// @}
     219             : 
     220             :   /// Retrieve a collection of hydrodynamic variables at position x
     221             :   template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
     222           1 :   tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
     223             :       const tnsr::I<DataType, 3>& x,
     224             :       tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
     225             :     return {tuples::get<Tag1>(variables(x, tmpl::list<Tag1>{})),
     226             :             tuples::get<Tag2>(variables(x, tmpl::list<Tag2>{})),
     227             :             tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
     228             :   }
     229             : 
     230             :   /// Retrieve the metric variables
     231             :   template <typename DataType, typename Tag,
     232             :             Requires<tmpl::list_contains_v<
     233             :                 gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
     234           1 :   tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
     235             :                                      tmpl::list<Tag> /*meta*/) const {
     236             :     constexpr double dummy_time = 0.0;
     237             :     return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
     238             :   }
     239             : 
     240           0 :   const EquationsOfState::IdealFluid<true>& equation_of_state() const {
     241             :     return equation_of_state_;
     242             :   }
     243             : 
     244             :   // NOLINTNEXTLINE(google-runtime-references)
     245           0 :   void pup(PUP::er& /*p*/) override;
     246             : 
     247             :  private:
     248           0 :   EquationsOfState::IdealFluid<true> equation_of_state_{};
     249           0 :   gr::Solutions::Minkowski<3> background_spacetime_{};
     250             : 
     251           0 :   double ambient_density_ = std::numeric_limits<double>::signaling_NaN();
     252           0 :   double ambient_pressure_ = std::numeric_limits<double>::signaling_NaN();
     253           0 :   double ambient_electron_fraction_ =
     254             :       std::numeric_limits<double>::signaling_NaN();
     255           0 :   double jet_density_ = std::numeric_limits<double>::signaling_NaN();
     256           0 :   double jet_pressure_ = std::numeric_limits<double>::signaling_NaN();
     257           0 :   double jet_electron_fraction_ = std::numeric_limits<double>::signaling_NaN();
     258           0 :   std::array<double, 3> jet_velocity_{
     259             :       {std::numeric_limits<double>::signaling_NaN(),
     260             :        std::numeric_limits<double>::signaling_NaN(),
     261             :        std::numeric_limits<double>::signaling_NaN()}};
     262           0 :   double inlet_radius_ = std::numeric_limits<double>::signaling_NaN();
     263           0 :   std::array<double, 3> magnetic_field_{
     264             :       {std::numeric_limits<double>::signaling_NaN(),
     265             :        std::numeric_limits<double>::signaling_NaN(),
     266             :        std::numeric_limits<double>::signaling_NaN()}};
     267             : 
     268           0 :   friend bool operator==(const SlabJet& lhs, const SlabJet& rhs);
     269             : 
     270           0 :   friend bool operator!=(const SlabJet& lhs, const SlabJet& rhs);
     271             : };
     272             : 
     273             : }  // namespace grmhd::AnalyticData