Allow for general magnetic field configs in MagnetizedTovStar

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp

@@ -16,27 +16,48 @@
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
 #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp"
 #include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "Utilities/Algorithm.hpp"
+#include "Utilities/CloneUniquePtrs.hpp"
 #include "Utilities/ConstantExpressions.hpp"
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/MakeWithValue.hpp"
 
 namespace grmhd::AnalyticData {
+MagnetizedTovStar::MagnetizedTovStar() = default;
+MagnetizedTovStar::MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/) = default;
+MagnetizedTovStar& MagnetizedTovStar::operator=(MagnetizedTovStar&& /*rhs*/) =
+    default;
+MagnetizedTovStar::~MagnetizedTovStar() = default;
+
 MagnetizedTovStar::MagnetizedTovStar(
     const double central_rest_mass_density,
     std::unique_ptr<MagnetizedTovStar::equation_of_state_type>
         equation_of_state,
     const RelativisticEuler::Solutions::TovCoordinates coordinate_system,
-    const size_t pressure_exponent, const double cutoff_pressure_fraction,
-    const double vector_potential_amplitude)
+    std::vector<std::unique_ptr<
+        grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
+        magnetic_fields)
     : tov_star(central_rest_mass_density, std::move(equation_of_state),
                coordinate_system),
-      pressure_exponent_(pressure_exponent),
-      cutoff_pressure_(cutoff_pressure_fraction *
-                       get(this->equation_of_state().pressure_from_density(
-                           Scalar<double>{central_rest_mass_density}))),
-      vector_potential_amplitude_(vector_potential_amplitude) {}
+      magnetic_fields_(std::move(magnetic_fields)) {}
+
+MagnetizedTovStar::MagnetizedTovStar(const MagnetizedTovStar& rhs)
+    : evolution::initial_data::InitialData{rhs},
+      RelativisticEuler::Solutions::TovStar(
+          static_cast<const RelativisticEuler::Solutions::TovStar&>(rhs)),
+      magnetic_fields_(clone_unique_ptrs(rhs.magnetic_fields_)) {}
+
+MagnetizedTovStar& MagnetizedTovStar::operator=(const MagnetizedTovStar& rhs) {
+  if (this == &rhs) {
+    return *this;
+  }
+  static_cast<RelativisticEuler::Solutions::TovStar&>(*this) =
+      static_cast<const RelativisticEuler::Solutions::TovStar&>(rhs);
+  magnetic_fields_ = clone_unique_ptrs(rhs.magnetic_fields_);
+  return *this;
+}
 
 std::unique_ptr<evolution::initial_data::InitialData>
 MagnetizedTovStar::get_clone() const {
@@ -47,9 +68,7 @@ MagnetizedTovStar::MagnetizedTovStar(CkMigrateMessage* msg) : tov_star(msg) {}
 
 void MagnetizedTovStar::pup(PUP::er& p) {
   tov_star::pup(p);
-  p | pressure_exponent_;
-  p | cutoff_pressure_;
-  p | vector_potential_amplitude_;
+  p | magnetic_fields_;
 }
 
 namespace magnetized_tov_detail {
@@ -58,68 +77,41 @@ void MagnetizedTovVariables<DataType, Region>::operator()(
     const gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
     const gsl::not_null<Cache*> cache,
     hydro::Tags::MagneticField<DataType, 3> /*meta*/) const {
-  const size_t num_pts = get_size(get<0>(coords));
-  const auto& pressure_profile =
-      get(cache->get_var(*this, hydro::Tags::Pressure<DataType>{}));
-  const auto& dr_pressure_profile = get(cache->get_var(
-      *this,
-      RelativisticEuler::Solutions::tov_detail::Tags::DrPressure<DataType>{}));
-  const auto& sqrt_det_spatial_metric =
-      get(cache->get_var(*this, gr::Tags::SqrtDetSpatialMetric<DataType>{}));
-  for (size_t i = 0; i < num_pts; ++i) {
-    const double pressure = get_element(pressure_profile, i);
-    if (LIKELY(get_element(radius, i) > 1.0e-16)) {
-      if (pressure < cutoff_pressure) {
-        get_element(get<0>(*magnetic_field), i) = 0.0;
-        get_element(get<1>(*magnetic_field), i) = 0.0;
-        get_element(get<2>(*magnetic_field), i) = 0.0;
-        continue;
-      }
-
-      const double x = get_element(get<0>(coords), i);
-      const double y = get_element(get<1>(coords), i);
-      const double z = get_element(get<2>(coords), i);
-      const double radius_i = get_element(radius, i);
-      const double dr_pressure = get_element(dr_pressure_profile, i);
-      const double pressure_term =
-          pow(pressure - cutoff_pressure, pressure_exponent);
-      const double deriv_pressure_term =
-          pressure_exponent *
-          pow(pressure - cutoff_pressure,
-              static_cast<int>(pressure_exponent) - 1) *
-          dr_pressure;
-
-      get_element(get<0>(*magnetic_field), i) =
-          -x * z / radius_i * deriv_pressure_term;
-
-      get_element(get<1>(*magnetic_field), i) =
-          -y * z / radius_i * deriv_pressure_term;
-
-      get_element(get<2>(*magnetic_field), i) =
-          2.0 * pressure_term +
-          (square(x) + square(y)) / radius_i * deriv_pressure_term;
-    } else {
-      get_element(get<0>(*magnetic_field), i) = 0.0;
-      get_element(get<1>(*magnetic_field), i) = 0.0;
-      get_element(get<2>(*magnetic_field), i) =
-          2.0 * pow(pressure - cutoff_pressure, pressure_exponent);
-    }
-  }
+  const Scalar<DataType>& sqrt_det_spatial_metric =
+      cache->get_var(*this, gr::Tags::SqrtDetSpatialMetric<DataType>{});
+  const Scalar<DataType>& pressure =
+      cache->get_var(*this, hydro::Tags::Pressure<DataType>{});
+  const auto& deriv_pressure =
+      cache->get_var(*this, ::Tags::deriv<hydro::Tags::Pressure<DataType>,
+                                          tmpl::size_t<3>, Frame::Inertial>{});
+  set_number_of_grid_points(magnetic_field, get_size(get<0>(coords)));
   for (size_t i = 0; i < 3; ++i) {
-    magnetic_field->get(i) *=
-        vector_potential_amplitude / sqrt_det_spatial_metric;
+    magnetic_field->get(i) = 0.0;
+  }
+
+  for (const auto& magnetic_field_compute : magnetic_fields) {
+    magnetic_field_compute->variables(magnetic_field, coords, pressure,
+                                      sqrt_det_spatial_metric, deriv_pressure);
   }
 }
 }  // namespace magnetized_tov_detail
 
 PUP::able::PUP_ID MagnetizedTovStar::my_PUP_ID = 0;
 
 bool operator==(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs) {
-  return static_cast<const typename MagnetizedTovStar::tov_star&>(lhs) ==
-             static_cast<const typename MagnetizedTovStar::tov_star&>(rhs) and
-         lhs.pressure_exponent_ == rhs.pressure_exponent_ and
-         lhs.cutoff_pressure_ == rhs.cutoff_pressure_ and
-         lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_;
+  bool equal =
+      static_cast<const typename MagnetizedTovStar::tov_star&>(lhs) ==
+          static_cast<const typename MagnetizedTovStar::tov_star&>(rhs) and
+      lhs.magnetic_fields_.size() == rhs.magnetic_fields_.size();
+  if (not equal) {
+    return false;
+  }
+  for (size_t i = 0; i < lhs.magnetic_fields_.size(); ++i) {
+    if (not lhs.magnetic_fields_[i]->is_equal(*rhs.magnetic_fields_[i])) {
+      return false;
+    }
+  }
+  return true;
 }
 
 bool operator!=(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs) {

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp

@@ -10,6 +10,8 @@
 #include "Options/String.hpp"
 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp"
 #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp"
 #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
@@ -42,20 +44,19 @@ struct MagnetizedTovVariables
   using Base::radial_solution;
   using Base::radius;
 
-  size_t pressure_exponent;
-  double cutoff_pressure;
-  double vector_potential_amplitude;
+  const std::vector<std::unique_ptr<
+      grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>&
+      magnetic_fields;
 
   MagnetizedTovVariables(
       const tnsr::I<DataType, 3>& local_x, const DataType& local_radius,
       const RelativisticEuler::Solutions::TovSolution& local_radial_solution,
       const EquationsOfState::EquationOfState<true, 1>& local_eos,
-      size_t local_pressure_exponent, double local_cutoff_pressure,
-      double local_vector_potential_amplitude)
+      const std::vector<std::unique_ptr<
+          grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>&
+          mag_fields)
       : Base(local_x, local_radius, local_radial_solution, local_eos),
-        pressure_exponent(local_pressure_exponent),
-        cutoff_pressure(local_cutoff_pressure),
-        vector_potential_amplitude(local_vector_potential_amplitude) {}
+        magnetic_fields(mag_fields) {}
 
   void operator()(
       gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
@@ -69,86 +70,11 @@ struct MagnetizedTovVariables
  * \brief Magnetized TOV star initial data, where metric terms only account for
  * the hydrodynamics not the magnetic fields.
  *
- * Superposes a poloidal magnetic field on top of a TOV solution where the
- * vector potential has the form
+ * Superposes magnetic fields on top of a TOV solution. These can be any of
+ * the classes derived from
+ * grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField
  *
- * \f{align*}{
- *  A_{\phi} = A_b \varpi^2 \max(p-p_{\mathrm{cut}}, 0)^{n_s}
- * \f}
- *
- * where \f$A_b\f$ controls the amplitude of the magnetic field,
- * \f$\varpi^2=x^2+y^2=r^2-z^2\f$ is the cylindrical radius,
- * \f$n_s\f$ controls the degree of differentiability, and
- * \f$p_{\mathrm{cut}}\f$ controls the pressure cutoff below which the magnetic
- * field is zero.
- *
- * In Cartesian coordinates the vector potential is:
- *
- * \f{align*}{
- *   A_x&=-\frac{y}{\varpi^2}A_\phi, \\
- *   A_y&=\frac{x}{\varpi^2}A_\phi, \\
- *   A_z&=0,
- * \f}
- *
- * For the poloidal field this means
- *
- * \f{align*}{
- *   A_x &= -y A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}, \\
- *   A_y &= x A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}.
- * \f}
- *
- * The magnetic field is computed from the vector potential as
- *
- * \f{align*}{
- *   B^i&=n_a\epsilon^{aijk}\partial_jA_k \\
- *      &=\frac{1}{\sqrt{\gamma}}[ijk]\partial_j A_k,
- * \f}
- *
- * where \f$[ijk]\f$ is the total anti-symmetric symbol. This means that
- *
- * \f{align*}{
- *   B^x&=\frac{1}{\sqrt{\gamma}} (\partial_y A_z - \partial_z A_y), \\
- *   B^y&=\frac{1}{\sqrt{\gamma}} (\partial_z A_x - \partial_x A_z), \\
- *   B^z&=\frac{1}{\sqrt{\gamma}} (\partial_x A_y - \partial_y A_x).
- * \f}
- *
- * Focusing on the region where the field is non-zero we have:
- *
- * \f{align*}{
- *   \partial_x A_y
- *   &= A_b(p-p_{\mathrm{cut}})^{n_s}+
- *     \frac{x^2}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p \\
- *   \partial_y A_x
- *   &= -A_b(p-p_{\mathrm{cut}})^{n_s}-
- *     \frac{y^2}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p \\
- *   \partial_z A_x
- *   &= -\frac{yz}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\
- *   \partial_z A_y
- *   &= \frac{xz}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp
- * \f}
- *
- * The magnetic field is given by:
- *
- * \f{align*}{
- *   B^x&=-\frac{1}{\sqrt{\gamma}}\frac{xz}{r}
- *        A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\
- *   B^y&=-\frac{1}{\sqrt{\gamma}}\frac{yz}{r}
- *        A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\
- *   B^z&=\frac{A_b}{\sqrt{\gamma}}\left[
- *        2(p-p_{\mathrm{cut}})^{n_s} \phantom{\frac{a}{b}}\right. \\
- *      &\left.+\frac{x^2+y^2}{r}
- *        n_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p
- *        \right]
- * \f}
- *
- * Taking the small-\f$r\f$ limit gives the \f$r=0\f$ magnetic field:
- *
- * \f{align*}{
- *   B^x&=0, \\
- *   B^y&=0, \\
- *   B^z&=\frac{A_b}{\sqrt{\gamma}}
- *        2(p-p_{\mathrm{cut}})^{n_s}.
- * \f}
+ * ### Conversion to CGS units and values for poloidal magnetic field
  *
  * While the amplitude \f$A_b\f$ is specified in the code, it is more natural
  * to work with the magnetic field strength, which is given by \f$\sqrt{b^2}\f$
@@ -161,7 +87,8 @@ struct MagnetizedTovVariables
  *   &= \sqrt{b^2} \times 8.352\times10^{19}\mathrm{G} \,.
  * \f}
  *
- * We now give values used for standard tests of magnetized stars.
+ * We now give values used for standard tests of magnetized stars with a
+ * poloidal magnetic field.
  * - \f$\rho_c(0)=1.28\times10^{-3}\f$
  * - \f$K=100\f$
  * - \f$\Gamma=2\f$
@@ -194,32 +121,14 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
   using tov_star = RelativisticEuler::Solutions::TovStar;
 
  public:
-  struct PressureExponent {
-    using type = size_t;
-    static constexpr Options::String help = {
-        "The exponent n_s controlling the smoothness of the field"};
-  };
-
-  struct VectorPotentialAmplitude {
-    using type = double;
-    static constexpr Options::String help = {
-        "The amplitude A_b of the phi-component of the vector potential. This "
-        "controls the magnetic field strength."};
-    static type lower_bound() { return 0.0; }
-  };
-
-  struct CutoffPressureFraction {
-    using type = double;
+  struct MagneticFields {
+    using type = std::vector<std::unique_ptr<
+        grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>;
     static constexpr Options::String help = {
-        "The fraction of the central pressure below which there is no magnetic "
-        "field. p_cut = Fraction * p_max."};
-    static type lower_bound() { return 0.0; }
-    static type upper_bound() { return 1.0; }
+        "Magnetic fields to superpose on the TOV solution."};
   };
 
-  using options =
-      tmpl::push_back<tov_star::options, PressureExponent,
-                      CutoffPressureFraction, VectorPotentialAmplitude>;
+  using options = tmpl::push_back<tov_star::options, MagneticFields>;
 
   static constexpr Options::String help = {"Magnetized TOV star."};
 
@@ -228,20 +137,21 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
   template <typename DataType>
   using tags = typename tov_star::template tags<DataType>;
 
-  MagnetizedTovStar() = default;
-  MagnetizedTovStar(const MagnetizedTovStar& /*rhs*/) = default;
-  MagnetizedTovStar& operator=(const MagnetizedTovStar& /*rhs*/) = default;
-  MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/) = default;
-  MagnetizedTovStar& operator=(MagnetizedTovStar&& /*rhs*/) = default;
-  ~MagnetizedTovStar() override = default;
+  MagnetizedTovStar();
+  MagnetizedTovStar(const MagnetizedTovStar& rhs);
+  MagnetizedTovStar& operator=(const MagnetizedTovStar& rhs);
+  MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/);
+  MagnetizedTovStar& operator=(MagnetizedTovStar&& /*rhs*/);
+  ~MagnetizedTovStar() override;
 
   MagnetizedTovStar(
       double central_rest_mass_density,
       std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
           equation_of_state,
       RelativisticEuler::Solutions::TovCoordinates coordinate_system,
-      size_t pressure_exponent, double cutoff_pressure_fraction,
-      double vector_potential_amplitude);
+      std::vector<std::unique_ptr<
+          grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
+          magnetic_fields);
 
   auto get_clone() const
       -> std::unique_ptr<evolution::initial_data::InitialData> override;
@@ -260,8 +170,7 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
                                          tmpl::list<Tags...> /*meta*/) const {
     return variables_impl<magnetized_tov_detail::MagnetizedTovVariables>(
-        x, tmpl::list<Tags...>{}, pressure_exponent_, cutoff_pressure_,
-        vector_potential_amplitude_);
+        x, tmpl::list<Tags...>{}, magnetic_fields_);
   }
 
   // NOLINTNEXTLINE(google-runtime-references)
@@ -272,10 +181,9 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
                          const MagnetizedTovStar& rhs);
 
  protected:
-  size_t pressure_exponent_ = std::numeric_limits<size_t>::max();
-  double cutoff_pressure_ = std::numeric_limits<double>::signaling_NaN();
-  double vector_potential_amplitude_ =
-      std::numeric_limits<double>::signaling_NaN();
+  std::vector<std::unique_ptr<
+      grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
+      magnetic_fields_{};
 };
 
 bool operator!=(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs);