Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <limits>
       8             : 
       9             : #include "DataStructures/Tensor/TypeAliases.hpp"
      10             : #include "Options/String.hpp"
      11             : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
      12             : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
      13             : #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp"
      14             : #include "PointwiseFunctions/GeneralRelativity/KerrSchildCoords.hpp"
      15             : #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
      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             : /// \cond
      23             : namespace PUP {
      24             : class er;
      25             : }  // namespace PUP
      26             : /// \endcond
      27             : 
      28             : namespace grmhd::AnalyticData {
      29             : 
      30             : /*!
      31             :  * \brief Magnetized fluid disk orbiting a Kerr black hole.
      32             :  *
      33             :  * In the context of simulating accretion disks, this class implements a widely
      34             :  * used (e.g. \cite Gammie2003, \cite Porth2016rfi, \cite White2015omx)
      35             :  * initial setup for the GRMHD variables, consisting of a Fishbone-Moncrief disk
      36             :  * \cite Fishbone1976apj (see also
      37             :  * RelativisticEuler::Solutions::FishboneMoncriefDisk),
      38             :  * threaded by a weak poloidal magnetic field. The magnetic field is constructed
      39             :  * from an axially symmetric toroidal magnetic potential which, in Kerr
      40             :  * ("spherical Kerr-Schild") coordinates, has the form
      41             :  *
      42             :  * \f{align*}
      43             :  * A_\phi(r,\theta) \propto \text{max}(\rho(r,\theta) - \rho_\text{thresh}, 0),
      44             :  * \f}
      45             :  *
      46             :  * where \f$\rho_\text{thresh}\f$ is a user-specified threshold density that
      47             :  * confines the magnetic flux to exist inside of the fluid disk only. A commonly
      48             :  * used value for this parameter is
      49             :  * \f$\rho_\text{thresh} = 0.2\rho_\text{max}\f$, where \f$\rho_\text{max}\f$
      50             :  * is the maximum value of
      51             :  * the rest mass density in the disk. Using this potential, the Eulerian
      52             :  * magnetic field takes the form
      53             :  *
      54             :  * \f{align*}
      55             :  * B^r = \frac{F_{\theta\phi}}{\sqrt{\gamma}},\quad
      56             :  * B^\theta = \frac{F_{\phi r}}{\sqrt{\gamma}},\quad B^\phi = 0,
      57             :  * \f}
      58             :  *
      59             :  * where \f$F_{ij} = \partial_i A_j - \partial_j A_i\f$ are the spatial
      60             :  * components of the Faraday tensor, and \f$\gamma\f$ is the determinant of the
      61             :  * spatial metric. The magnetic field is then normalized so that the
      62             :  * plasma-\f$\beta\f$ parameter, \f$\beta = 2p/b^2\f$, equals some value
      63             :  * specified by the user. Here, \f$p\f$ is the fluid pressure, and
      64             :  *
      65             :  * \f{align*}
      66             :  * b^2 = b^\mu b_\mu = \frac{B_iB^i}{W^2} + (B^iv_i)^2
      67             :  * \f}
      68             :  *
      69             :  * is the norm of the magnetic field in the fluid frame, with \f$v_i\f$ being
      70             :  * the spatial velocity, and \f$W\f$ the Lorentz factor.
      71             :  *
      72             :  * \note When using Kerr-Schild coordinates, the horizon that is at
      73             :  * constant \f$r\f$ is not spherical, but instead spheroidal. This could make
      74             :  * application of boundary condition and computing various fluxes
      75             :  * across the horizon more complicated than they need to be.
      76             :  * Thus, similar to RelativisticEuler::Solutions::FishboneMoncriefDisk
      77             :  * we use Spherical Kerr-Schild coordinates,
      78             :  * see gr::Solutions::SphericalKerrSchild, in which constant \f$r\f$
      79             :  * is spherical. Because we compute variables in Kerr-Schild coordinates,
      80             :  * there is a necessary extra step of transforming them back to
      81             :  * Spherical Kerr-Schild coordinates.
      82             :  *
      83             :  * \warning Spherical Kerr-Schild coordinates and "spherical Kerr-Schild"
      84             :  * coordinates are not same.
      85             :  *
      86             :  */
      87           1 : class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData,
      88             :                          public MarkAsAnalyticData {
      89             :  private:
      90           0 :   using FmDisk = RelativisticEuler::Solutions::FishboneMoncriefDisk;
      91             : 
      92             :  public:
      93             :   /// The rest mass density (in units of the maximum rest mass density in the
      94             :   /// disk) below which the matter in the disk is initially unmagetized.
      95           1 :   struct ThresholdDensity {
      96           0 :     using type = double;
      97           0 :     static constexpr Options::String help = {
      98             :         "Frac. rest mass density below which B-field vanishes."};
      99           0 :     static type lower_bound() { return 0.0; }
     100           0 :     static type upper_bound() { return 1.0; }
     101             :   };
     102             :   /// The maximum-magnetic-pressure-to-maximum-fluid-pressure ratio.
     103           1 :   struct InversePlasmaBeta {
     104           0 :     using type = double;
     105           0 :     static constexpr Options::String help = {
     106             :         "Ratio of max magnetic pressure to max fluid pressure."};
     107           0 :     static type lower_bound() { return 0.0; }
     108             :   };
     109             :   /// Grid resolution used in magnetic field normalization.
     110           1 :   struct BFieldNormGridRes {
     111           0 :     using type = size_t;
     112           0 :     static constexpr Options::String help = {
     113             :         "Grid Resolution for b-field normalization."};
     114           0 :     static type suggested_value() { return 255; }
     115           0 :     static type lower_bound() { return 4; }
     116             :   };
     117             : 
     118             :   // Unlike the other analytic data classes, we cannot get these from the
     119             :   // `AnalyticDataBase` because this case causes clang-tidy to believe that
     120             :   // there is an ambiguous inheritance problem
     121           0 :   static constexpr size_t volume_dim = 3_st;
     122             : 
     123             :   template <typename DataType>
     124           0 :   using tags =
     125             :       tmpl::push_back<typename gr::AnalyticSolution<3>::template tags<DataType>,
     126             :                       hydro::Tags::RestMassDensity<DataType>,
     127             :                       hydro::Tags::ElectronFraction<DataType>,
     128             :                       hydro::Tags::SpecificInternalEnergy<DataType>,
     129             :                       hydro::Tags::Temperature<DataType>,
     130             :                       hydro::Tags::Pressure<DataType>,
     131             :                       hydro::Tags::SpatialVelocity<DataType, 3>,
     132             :                       hydro::Tags::MagneticField<DataType, 3>,
     133             :                       hydro::Tags::DivergenceCleaningField<DataType>,
     134             :                       hydro::Tags::LorentzFactor<DataType>,
     135             :                       hydro::Tags::SpecificEnthalpy<DataType>>;
     136             : 
     137           0 :   using options = tmpl::push_back<FmDisk::options, ThresholdDensity,
     138             :                                   InversePlasmaBeta, BFieldNormGridRes>;
     139             : 
     140           0 :   static constexpr Options::String help = {
     141             :       "Magnetized Fishbone-Moncrief disk."};
     142             : 
     143           0 :   MagnetizedFmDisk() = default;
     144           0 :   MagnetizedFmDisk(const MagnetizedFmDisk& /*rhs*/) = default;
     145           0 :   MagnetizedFmDisk& operator=(const MagnetizedFmDisk& /*rhs*/) = default;
     146           0 :   MagnetizedFmDisk(MagnetizedFmDisk&& /*rhs*/) = default;
     147           0 :   MagnetizedFmDisk& operator=(MagnetizedFmDisk&& /*rhs*/) = default;
     148           0 :   ~MagnetizedFmDisk() override = default;
     149             : 
     150           0 :   MagnetizedFmDisk(
     151             :       double bh_mass, double bh_dimless_spin, double inner_edge_radius,
     152             :       double max_pressure_radius, double polytropic_constant,
     153             :       double polytropic_exponent, double noise, double threshold_density,
     154             :       double inverse_plasma_beta,
     155             :       size_t normalization_grid_res = BFieldNormGridRes::suggested_value());
     156             : 
     157           0 :   auto get_clone() const
     158             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     159             : 
     160             :   /// \cond
     161             :   explicit MagnetizedFmDisk(CkMigrateMessage* msg);
     162             :   using PUP::able::register_constructor;
     163             :   WRAPPED_PUPable_decl_template(MagnetizedFmDisk);
     164             :   /// \endcond
     165             : 
     166             :   // Overload the variables function from the base class.
     167           0 :   using equation_of_state_type = typename FmDisk::equation_of_state_type;
     168             : 
     169             :   /// @{
     170             :   /// The variables in Cartesian Kerr-Schild coordinates at `(x, t)`.
     171             :   template <typename DataType, typename... Tags>
     172           1 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
     173             :                                          tmpl::list<Tags...> /*meta*/) const {
     174             :     // Can't store IntermediateVariables as member variable because we
     175             :     // need to be threadsafe.
     176             :     typename FmDisk::IntermediateVariables<DataType> vars(x);
     177             :     return {std::move(get<Tags>([this, &x, &vars]() {
     178             :       if constexpr (std::is_same_v<hydro::Tags::MagneticField<DataType, 3>,
     179             :                                    Tags>) {
     180             :         return variables(x, tmpl::list<Tags>{}, make_not_null(&vars));
     181             :       } else {
     182             :         return fm_disk_.variables(x, tmpl::list<Tags>{}, make_not_null(&vars));
     183             :       }
     184             :     }()))...};
     185             :   }
     186             : 
     187             :   template <typename DataType, typename Tag>
     188           1 :   tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
     189             :                                      tmpl::list<Tag> /*meta*/) const {
     190             :     // Can't store IntermediateVariables as member variable because we need to
     191             :     // be threadsafe.
     192             :     typename FmDisk::IntermediateVariables<DataType> vars(x);
     193             :     if constexpr (std::is_same_v<hydro::Tags::MagneticField<DataType, 3>,
     194             :                                  Tag>) {
     195             :       return variables(x, tmpl::list<Tag>{}, make_not_null(&vars));
     196             :     } else {
     197             :       return fm_disk_.variables(x, tmpl::list<Tag>{}, make_not_null(&vars));
     198             :     }
     199             :   }
     200             :   /// @}
     201             : 
     202             :   // NOLINTNEXTLINE(google-runtime-references)
     203           0 :   void pup(PUP::er& /*p*/) override;
     204             : 
     205           0 :   const equation_of_state_type& equation_of_state() const {
     206             :     return fm_disk_.equation_of_state();
     207             :   }
     208             : 
     209             :  private:
     210             :   template <typename DataType>
     211           0 :   auto variables(const tnsr::I<DataType, 3>& x,
     212             :                  tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
     213             :                  gsl::not_null<FmDisk::IntermediateVariables<DataType>*> vars)
     214             :       const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
     215             : 
     216             :   template <typename DataType>
     217           0 :   tnsr::I<DataType, 3> unnormalized_magnetic_field(
     218             :       const tnsr::I<DataType, 3>& x) const;
     219             : 
     220           0 :   friend bool operator==(const MagnetizedFmDisk& lhs,
     221             :                          const MagnetizedFmDisk& rhs);
     222             : 
     223           0 :   RelativisticEuler::Solutions::FishboneMoncriefDisk fm_disk_{};
     224           0 :   double threshold_density_ = std::numeric_limits<double>::signaling_NaN();
     225           0 :   double inverse_plasma_beta_ = std::numeric_limits<double>::signaling_NaN();
     226           0 :   double b_field_normalization_ = std::numeric_limits<double>::signaling_NaN();
     227           0 :   size_t normalization_grid_res_ = 255;
     228           0 :   gr::KerrSchildCoords kerr_schild_coords_{};
     229             : };
     230             : 
     231           0 : bool operator!=(const MagnetizedFmDisk& lhs, const MagnetizedFmDisk& rhs);
     232             : 
     233             : }  // namespace grmhd::AnalyticData