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