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/GrMhd/InitialMagneticFields/InitialMagneticField.hpp"
      12             : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
      13             : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
      14             : #include "Utilities/TMPL.hpp"
      15             : #include "Utilities/TaggedTuple.hpp"
      16             : 
      17             : /// \cond
      18             : namespace gsl {
      19             : template <typename T>
      20             : class not_null;
      21             : }  // namespace gsl
      22             : /// \endcond
      23             : 
      24             : namespace grmhd::AnalyticData::InitialMagneticFields {
      25             : 
      26             : /*!
      27             :  * \brief %Poloidal magnetic field for GRMHD initial data.
      28             :  *
      29             :  * The vector potential has the form
      30             :  *
      31             :  * \f{align*}{
      32             :  *  A_{\phi} = A_b \varpi^2 \max(p-p_{\mathrm{cut}}, 0)^{n_s} ,
      33             :  * \f}
      34             :  *
      35             :  * where \f$A_b\f$ controls the amplitude of the magnetic field,
      36             :  * \f$\varpi^2=x^2+y^2=r^2-z^2\f$ is the cylindrical radius,
      37             :  * \f$n_s\f$ controls the degree of differentiability, and
      38             :  * \f$p_{\mathrm{cut}}\f$ controls the pressure cutoff below which the magnetic
      39             :  * field is zero.
      40             :  *
      41             :  * In Cartesian coordinates the vector potential is:
      42             :  *
      43             :  * \f{align*}{
      44             :  *   A_x & = -\frac{y}{\varpi^2}A_\phi
      45             :  *            = -y A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}, \\
      46             :  *   A_y & = \frac{x}{\varpi^2}A_\phi
      47             :  *            = x A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}, \\
      48             :  *   A_z & = 0,
      49             :  * \f}
      50             :  *
      51             :  * On the region where the field is non-zero, the magnetic field is given by:
      52             :  *
      53             :  * \f{align*}{
      54             :  *   B^x & = -\frac{A_b n_s}{\sqrt{\gamma}}
      55             :  *        (p-p_{\mathrm{cut}})^{n_s-1} x \partial_z p \\
      56             :  *   B^x & = -\frac{A_b n_s}{\sqrt{\gamma}}
      57             :  *        (p-p_{\mathrm{cut}})^{n_s-1} y \partial_z p \\
      58             :  *   B^z & = \frac{A_b}{\sqrt{\gamma}}\left[
      59             :  *        2(p-p_{\mathrm{cut}})^{n_s}
      60             :  *        + n_s (p-p_{\mathrm{cut}})^{n_s-1} (x \partial_x p + y \partial_y p)
      61             :  *        \right]
      62             :  * \f}
      63             :  *
      64             :  * Taking the small-\f$r\f$ limit gives the magnetic field at the origin:
      65             :  *
      66             :  * \f{align*}{
      67             :  *   B^x&=0, \\
      68             :  *   B^y&=0, \\
      69             :  *   B^z&=\frac{A_b}{\sqrt{\gamma}}
      70             :  *        2(p-p_{\mathrm{cut}})^{n_s}.
      71             :  * \f}
      72             :  *
      73             :  * Note that the coordinates are relative to the `Center` passed in, so the
      74             :  * field can be centered about any arbitrary point. The field is also zero
      75             :  * outside of `MaxDistanceFromCenter`, so that compact support can be imposed if
      76             :  * necessary.
      77             :  *
      78             :  * \warning This assumes the magnetic field is initialized, both in size and
      79             :  * value, before being passed into the `variables` function. This is so that
      80             :  * multiple magnetic fields can be superposed. Each magnetic field
      81             :  * configuration does a `+=` to make this possible.
      82             :  */
      83           1 : class Poloidal : public InitialMagneticField {
      84             :  public:
      85           0 :   struct PressureExponent {
      86           0 :     using type = size_t;
      87           0 :     static constexpr Options::String help = {
      88             :         "The exponent n_s controlling the smoothness of the field"};
      89             :   };
      90             : 
      91           0 :   struct CutoffPressure {
      92           0 :     using type = double;
      93           0 :     static constexpr Options::String help = {
      94             :         "The pressure below which there is no magnetic field."};
      95           0 :     static type lower_bound() { return 0.0; }
      96             :   };
      97             : 
      98           0 :   struct VectorPotentialAmplitude {
      99           0 :     using type = double;
     100           0 :     static constexpr Options::String help = {
     101             :         "The amplitude A_b of the phi-component of the vector potential. This "
     102             :         "controls the magnetic field strength."};
     103           0 :     static type lower_bound() { return 0.0; }
     104             :   };
     105             : 
     106           0 :   struct Center {
     107           0 :     using type = std::array<double, 3>;
     108           0 :     static constexpr Options::String help = {
     109             :         "The center of the magnetic field."};
     110             :   };
     111             : 
     112           0 :   struct MaxDistanceFromCenter {
     113           0 :     using type = double;
     114           0 :     static constexpr Options::String help = {
     115             :         "The maximum distance from the center to compute the magnetic field. "
     116             :         "Everywhere outside the field is set to zero."};
     117           0 :     static type lower_bound() { return 0.0; }
     118             :   };
     119             : 
     120           0 :   using options =
     121             :       tmpl::list<PressureExponent, CutoffPressure, VectorPotentialAmplitude,
     122             :                  Center, MaxDistanceFromCenter>;
     123             : 
     124           0 :   static constexpr Options::String help = {"Poloidal initial magnetic field"};
     125             : 
     126           0 :   Poloidal() = default;
     127           0 :   Poloidal(const Poloidal& /*rhs*/) = default;
     128           0 :   Poloidal& operator=(const Poloidal& /*rhs*/) = default;
     129           0 :   Poloidal(Poloidal&& /*rhs*/) = default;
     130           0 :   Poloidal& operator=(Poloidal&& /*rhs*/) = default;
     131           0 :   ~Poloidal() override = default;
     132             : 
     133           0 :   Poloidal(size_t pressure_exponent, double cutoff_pressure,
     134             :            double vector_potential_amplitude, std::array<double, 3> center,
     135             :            double max_distance_from_center);
     136             : 
     137           0 :   auto get_clone() const -> std::unique_ptr<InitialMagneticField> override;
     138             : 
     139             :   /// \cond
     140             :   explicit Poloidal(CkMigrateMessage* msg);
     141             :   using PUP::able::register_constructor;
     142             :   WRAPPED_PUPable_decl_template(Poloidal);
     143             :   /// \endcond
     144             : 
     145             :   // NOLINTNEXTLINE(google-runtime-references)
     146           0 :   void pup(PUP::er& p) override;
     147             : 
     148             :   /// Retrieve magnetic fields at `(x)`
     149           1 :   void variables(gsl::not_null<tnsr::I<DataVector, 3>*> result,
     150             :                  const tnsr::I<DataVector, 3>& coords,
     151             :                  const Scalar<DataVector>& pressure,
     152             :                  const Scalar<DataVector>& sqrt_det_spatial_metric,
     153             :                  const tnsr::i<DataVector, 3>& deriv_pressure) const override;
     154             : 
     155             :   /// Retrieve magnetic fields at `(x)`
     156           1 :   void variables(gsl::not_null<tnsr::I<double, 3>*> result,
     157             :                  const tnsr::I<double, 3>& coords,
     158             :                  const Scalar<double>& pressure,
     159             :                  const Scalar<double>& sqrt_det_spatial_metric,
     160             :                  const tnsr::i<double, 3>& deriv_pressure) const override;
     161             : 
     162           0 :   bool is_equal(const InitialMagneticField& rhs) const override;
     163             : 
     164             :  private:
     165             :   template <typename DataType>
     166           0 :   void variables_impl(gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
     167             :                       const tnsr::I<DataType, 3>& coords,
     168             :                       const Scalar<DataType>& pressure,
     169             :                       const Scalar<DataType>& sqrt_det_spatial_metric,
     170             :                       const tnsr::i<DataType, 3>& deriv_pressure) const;
     171             : 
     172           0 :   size_t pressure_exponent_ = std::numeric_limits<size_t>::max();
     173           0 :   double cutoff_pressure_ = std::numeric_limits<double>::signaling_NaN();
     174           0 :   double vector_potential_amplitude_ =
     175             :       std::numeric_limits<double>::signaling_NaN();
     176           0 :   std::array<double, 3> center_{{std::numeric_limits<double>::signaling_NaN(),
     177             :                                  std::numeric_limits<double>::signaling_NaN(),
     178             :                                  std::numeric_limits<double>::signaling_NaN()}};
     179           0 :   double max_distance_from_center_ =
     180             :       std::numeric_limits<double>::signaling_NaN();
     181             : 
     182           0 :   friend bool operator==(const Poloidal& lhs, const Poloidal& rhs);
     183           0 :   friend bool operator!=(const Poloidal& lhs, const Poloidal& rhs);
     184             : };
     185             : 
     186             : }  // namespace grmhd::AnalyticData::InitialMagneticFields