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