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             : 
       9             : #include "DataStructures/Tensor/TypeAliases.hpp"
      10             : #include "Options/Context.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/IdealFluid.hpp"  // IWYU pragma: keep
      16             : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
      17             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      18             : #include "Utilities/Serialization/CharmPupable.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : #include "Utilities/TaggedTuple.hpp"
      21             : 
      22             : // IWYU pragma:  no_include <pup.h>
      23             : 
      24             : /// \cond
      25             : namespace PUP {
      26             : class er;  // IWYU pragma: keep
      27             : }  // namespace PUP
      28             : /// \endcond
      29             : 
      30             : namespace grmhd::AnalyticData {
      31             : 
      32             : /*!
      33             :  * \brief Analytic initial data for an advecting magnetic field loop.
      34             :  *
      35             :  * This test, originally proposed in \cite Gardiner2005hy and presented in a
      36             :  * slightly modified form by \cite Mignone2010br, a region with annular
      37             :  * cross section with the specified `InnerRadius` and `OuterRadius` is given a
      38             :  * non-zero azimuthal magnetic field of constant magnitude `MagFieldStrength`
      39             :  * with zero magnetic field outside the loop.  Inside the `InnerRadius` the
      40             :  * magnetic field strength falls to zero quadratically. The loop is embedded in
      41             :  * an ideal fluid with the given `AdiabaticIndex`, `RestMassDensity` and
      42             :  * `Pressure` with a uniform `AdvectionVelocity`.  The magnetic field loop
      43             :  * should advect across the grid, maintaining its shape and strength, as long
      44             :  * as the magnetic pressure is negligible compared to the thermal pressure.
      45             :  *
      46             :  * This test diagnoses how well the evolution scheme preserves the no-monopole
      47             :  * condition, as well as the diffusivity of the scheme.
      48             :  *
      49             :  * The standard test setup is done on \f$x \in [-1,1]\f$, \f$y \in [-0.5,
      50             :  * 0.5]\f$, with periodic boundary conditions and with the following values
      51             :  * given for the options:
      52             :  * -  InnerRadius: 0.06
      53             :  * -  OuterRadius: 0.3
      54             :  * -  RestMassDensity: 1.0
      55             :  * -  Pressure: 1.0
      56             :  * -  AdvectionVelocity: [0.08164965809277261, 0.040824829046386304, 0.0]
      57             :  * -  MagFieldStrength: 0.001
      58             :  * -  AdiabaticIndex: 1.66666666666666667
      59             :  *
      60             :  */
      61           1 : class MagneticFieldLoop
      62             :     : public evolution::initial_data::InitialData,
      63             :       public MarkAsAnalyticData,
      64             :       public AnalyticDataBase,
      65             :       public hydro::TemperatureInitialization<MagneticFieldLoop> {
      66             :  public:
      67           0 :   using equation_of_state_type = EquationsOfState::IdealFluid<true>;
      68             : 
      69             :   /// The pressure throughout the fluid.
      70           1 :   struct Pressure {
      71           0 :     using type = double;
      72           0 :     static constexpr Options::String help = {
      73             :         "The constant pressure throughout the fluid."};
      74           0 :     static type lower_bound() { return 0.0; }
      75             :   };
      76             : 
      77             :   /// The rest mass density throughout the fluid.
      78           1 :   struct RestMassDensity {
      79           0 :     using type = double;
      80           0 :     static constexpr Options::String help = {
      81             :         "The constant density throughout the fluid."};
      82           0 :     static type lower_bound() { return 0.0; }
      83             :   };
      84             : 
      85             :   /// The adiabatic index for the ideal fluid.
      86           1 :   struct AdiabaticIndex {
      87           0 :     using type = double;
      88           0 :     static constexpr Options::String help = {
      89             :         "The adiabatic index for the ideal fluid."};
      90           0 :     static type lower_bound() { return 1.0; }
      91             :   };
      92             : 
      93             :   /// The fluid velocity.
      94           1 :   struct AdvectionVelocity {
      95           0 :     using type = std::array<double, 3>;
      96           0 :     static constexpr Options::String help = {"The advection velocity."};
      97           0 :     static type lower_bound() { return {{-1.0, -1.0, -1.0}}; }
      98           0 :     static type upper_bound() { return {{1.0, 1.0, 1.0}}; }
      99             :   };
     100             : 
     101             :   /// The strength of the magnetic field.
     102           1 :   struct MagFieldStrength {
     103           0 :     using type = double;
     104           0 :     static constexpr Options::String help = {
     105             :         "The magnitude of the magnetic field."};
     106           0 :     static type lower_bound() { return 0.0; }
     107             :   };
     108             : 
     109             :   /// The inner radius of the magnetic loop.
     110           1 :   struct InnerRadius {
     111           0 :     using type = double;
     112           0 :     static constexpr Options::String help = {
     113             :         "The inner radius of the magnetic loop."};
     114           0 :     static type lower_bound() { return 0.0; }
     115             :   };
     116             : 
     117             :   /// The outer radius of the magnetic loop.
     118           1 :   struct OuterRadius {
     119           0 :     using type = double;
     120           0 :     static constexpr Options::String help = {
     121             :         "The outer radius of the magnetic loop."};
     122           0 :     static type lower_bound() { return 0.0; }
     123             :   };
     124             : 
     125           0 :   using options =
     126             :       tmpl::list<Pressure, RestMassDensity, AdiabaticIndex, AdvectionVelocity,
     127             :                  MagFieldStrength, InnerRadius, OuterRadius>;
     128           0 :   static constexpr Options::String help = {
     129             :       "Periodic advection of a magnetic field loop in Minkowski."};
     130             : 
     131           0 :   MagneticFieldLoop() = default;
     132           0 :   MagneticFieldLoop(const MagneticFieldLoop& /*rhs*/) = default;
     133           0 :   MagneticFieldLoop& operator=(const MagneticFieldLoop& /*rhs*/) = default;
     134           0 :   MagneticFieldLoop(MagneticFieldLoop&& /*rhs*/) = default;
     135           0 :   MagneticFieldLoop& operator=(MagneticFieldLoop&& /*rhs*/) = default;
     136           0 :   ~MagneticFieldLoop() override = default;
     137             : 
     138           0 :   MagneticFieldLoop(double pressure, double rest_mass_density,
     139             :                     double adiabatic_index,
     140             :                     const std::array<double, 3>& advection_velocity,
     141             :                     double magnetic_field_magnitude, double inner_radius,
     142             :                     double outer_radius, const Options::Context& context = {});
     143             : 
     144           0 :   auto get_clone() const
     145             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     146             : 
     147             :   /// \cond
     148             :   explicit MagneticFieldLoop(CkMigrateMessage* msg);
     149             :   using PUP::able::register_constructor;
     150             :   WRAPPED_PUPable_decl_template(MagneticFieldLoop);
     151             :   /// \endcond
     152             : 
     153             :   /// @{
     154             :   /// Retrieve the GRMHD variables at a given position.
     155             :   template <typename DataType>
     156           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     157             :                  tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
     158             :       const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
     159             : 
     160             :   template <typename DataType>
     161           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     162             :                  tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
     163             :       const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
     164             : 
     165             :   template <typename DataType>
     166           1 :   auto variables(
     167             :       const tnsr::I<DataType, 3>& x,
     168             :       tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
     169             :       -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
     170             : 
     171             :   template <typename DataType>
     172           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     173             :                  tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
     174             :       -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
     175             : 
     176             :   template <typename DataType>
     177           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     178             :                  tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
     179             :       const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
     180             : 
     181             :   template <typename DataType>
     182           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     183             :                  tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
     184             :       const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
     185             : 
     186             :   template <typename DataType>
     187           1 :   auto variables(
     188             :       const tnsr::I<DataType, 3>& x,
     189             :       tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
     190             :       -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
     191             : 
     192             :   template <typename DataType>
     193           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     194             :                  tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
     195             :       const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
     196             : 
     197             :   template <typename DataType>
     198           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     199             :                  tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
     200             :       const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
     201             : 
     202             :   template <typename DataType>
     203           1 :   auto variables(const tnsr::I<DataType, 3>& x,
     204             :                  tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
     205             :       -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
     206             :     return TemperatureInitialization::variables(
     207             :         x, tmpl::list<hydro::Tags::Temperature<DataType>>{});
     208             :   }
     209             :   /// @}
     210             : 
     211             :   /// Retrieve a collection of hydrodynamic variables at position x
     212             :   template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
     213           1 :   tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
     214             :       const tnsr::I<DataType, 3>& x,
     215             :       tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
     216             :     return {tuples::get<Tag1>(variables(x, tmpl::list<Tag1>{})),
     217             :             tuples::get<Tag2>(variables(x, tmpl::list<Tag2>{})),
     218             :             tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
     219             :   }
     220             : 
     221             :   /// Retrieve the metric variables
     222             :   template <typename DataType, typename Tag,
     223             :             Requires<tmpl::list_contains_v<
     224             :                 gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
     225           1 :   tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
     226             :                                      tmpl::list<Tag> /*meta*/) const {
     227             :     constexpr double dummy_time = 0.0;
     228             :     return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
     229             :   }
     230             : 
     231           0 :   const EquationsOfState::IdealFluid<true>& equation_of_state() const {
     232             :     return equation_of_state_;
     233             :   }
     234             : 
     235             :   // NOLINTNEXTLINE(google-runtime-references)
     236           0 :   void pup(PUP::er& /*p*/) override;
     237             : 
     238             :  private:
     239           0 :   double pressure_ = std::numeric_limits<double>::signaling_NaN();
     240           0 :   double rest_mass_density_ = std::numeric_limits<double>::signaling_NaN();
     241           0 :   double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
     242           0 :   std::array<double, 3> advection_velocity_{
     243             :       {std::numeric_limits<double>::signaling_NaN(),
     244             :        std::numeric_limits<double>::signaling_NaN(),
     245             :        std::numeric_limits<double>::signaling_NaN()}};
     246           0 :   double magnetic_field_magnitude_ =
     247             :       std::numeric_limits<double>::signaling_NaN();
     248           0 :   double inner_radius_ = std::numeric_limits<double>::signaling_NaN();
     249           0 :   double outer_radius_ = std::numeric_limits<double>::signaling_NaN();
     250             : 
     251           0 :   EquationsOfState::IdealFluid<true> equation_of_state_{};
     252           0 :   gr::Solutions::Minkowski<3> background_spacetime_{};
     253             : 
     254           0 :   friend bool operator==(const MagneticFieldLoop& lhs,
     255             :                          const MagneticFieldLoop& rhs);
     256             : 
     257           0 :   friend bool operator!=(const MagneticFieldLoop& lhs,
     258             :                          const MagneticFieldLoop& rhs);
     259             : };
     260             : 
     261             : }  // namespace grmhd::AnalyticData