Make magnetic field configuration runtime choice

src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/CMakeLists.txt

@@ -12,6 +12,7 @@ spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
+  Factory.hpp
   InitialMagneticField.hpp
   Poloidal.hpp
   Toroidal.hpp

src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp

@@ -0,0 +1,8 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Toroidal.hpp"

src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp

@@ -6,7 +6,22 @@
 #include <memory>
 #include <pup.h>
 
+#include "DataStructures/Tensor/TypeAliases.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+class DataVector;
+namespace grmhd::AnalyticData::InitialMagneticFields {
+class Poloidal;
+class Toroidal;
+}  // namespace grmhd::AnalyticData::InitialMagneticFields
+
+namespace gsl {
+template <typename T>
+class not_null;
+}  // namespace gsl
+/// \endcond
 
 namespace grmhd::AnalyticData {
 
@@ -32,21 +47,47 @@ namespace grmhd::AnalyticData {
  *   B^z & = \frac{1}{\sqrt{\gamma}} (\partial_x A_y - \partial_y A_x).
  * \f}
  *
+ *
+ * \warning The magnetic field classes assume the magnetic field is initialized,
+ * both in size and value, before being passed into the `variables` function.
+ * This is so that multiple magnetic fields can be superposed. Each magnetic
+ * field configuration does a `+=` to make this possible.
  */
 namespace InitialMagneticFields {
 
 /*!
  * \brief The abstract base class for initial magnetic field configurations.
+ *
+ * \warning This assumes the magnetic field is initialized, both in size and
+ * value, before being passed into the `variables` function. This is so that
+ * multiple magnetic fields can be superposed. Each magnetic field
+ * configuration does a `+=` to make this possible.
  */
 class InitialMagneticField : public PUP::able {
  protected:
   InitialMagneticField() = default;
 
  public:
+  using creatable_classes = tmpl::list<Poloidal, Toroidal>;
+
   ~InitialMagneticField() override = default;
 
   virtual auto get_clone() const -> std::unique_ptr<InitialMagneticField> = 0;
 
+  virtual void variables(
+      gsl::not_null<tnsr::I<DataVector, 3>*> result,
+      const tnsr::I<DataVector, 3>& coords, const Scalar<DataVector>& pressure,
+      const Scalar<DataVector>& sqrt_det_spatial_metric,
+      const tnsr::i<DataVector, 3>& deriv_pressure) const = 0;
+
+  virtual void variables(gsl::not_null<tnsr::I<double, 3>*> result,
+                         const tnsr::I<double, 3>& coords,
+                         const Scalar<double>& pressure,
+                         const Scalar<double>& sqrt_det_spatial_metric,
+                         const tnsr::i<double, 3>& deriv_pressure) const = 0;
+
+  virtual bool is_equal(const InitialMagneticField& rhs) const = 0;
+
   /// \cond
   explicit InitialMagneticField(CkMigrateMessage* msg) : PUP::able(msg) {}
   WRAPPED_PUPable_abstract(InitialMagneticField);

src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.cpp

@@ -24,72 +24,103 @@ void Poloidal::pup(PUP::er& p) {
   p | pressure_exponent_;
   p | cutoff_pressure_;
   p | vector_potential_amplitude_;
+  p | center_;
+  p | max_distance_from_center_;
 }
 
 // NOLINTNEXTLINE
 PUP::able::PUP_ID Poloidal::my_PUP_ID = 0;
 
 Poloidal::Poloidal(const size_t pressure_exponent, const double cutoff_pressure,
-                   const double vector_potential_amplitude)
+                   const double vector_potential_amplitude,
+                   const std::array<double, 3> center,
+                   const double max_distance_from_center)
     : pressure_exponent_(pressure_exponent),
       cutoff_pressure_(cutoff_pressure),
-      vector_potential_amplitude_(vector_potential_amplitude) {}
+      vector_potential_amplitude_(vector_potential_amplitude),
+      center_(center),
+      max_distance_from_center_(max_distance_from_center) {}
+
+void Poloidal::variables(const gsl::not_null<tnsr::I<DataVector, 3>*> result,
+                         const tnsr::I<DataVector, 3>& coords,
+                         const Scalar<DataVector>& pressure,
+                         const Scalar<DataVector>& sqrt_det_spatial_metric,
+                         const tnsr::i<DataVector, 3>& deriv_pressure) const {
+  ASSERT(result->get(0).size() == get(pressure).size(),
+         "Result must be of size " << get(pressure).size() << " but got "
+                                   << result->get(0).size());
+  variables_impl(result, coords, pressure, sqrt_det_spatial_metric,
+                 deriv_pressure);
+}
+
+void Poloidal::variables(const gsl::not_null<tnsr::I<double, 3>*> result,
+                         const tnsr::I<double, 3>& coords,
+                         const Scalar<double>& pressure,
+                         const Scalar<double>& sqrt_det_spatial_metric,
+                         const tnsr::i<double, 3>& deriv_pressure) const {
+  variables_impl(result, coords, pressure, sqrt_det_spatial_metric,
+                 deriv_pressure);
+}
+
+bool Poloidal::is_equal(const InitialMagneticField& rhs) const {
+  const auto& derived_ptr = dynamic_cast<const Poloidal* const>(&rhs);
+  return derived_ptr != nullptr and *derived_ptr == *this;
+}
 
 template <typename DataType>
-tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>
-Poloidal::variables(const tnsr::I<DataType, 3>& coords,
-                    const Scalar<DataType>& pressure,
-                    const Scalar<DataType>& sqrt_det_spatial_metric,
-                    const tnsr::i<DataType, 3>& dcoords_pressure) const {
-  auto magnetic_field = make_with_value<tnsr::I<DataType, 3>>(coords, 0.0);
+void Poloidal::variables_impl(
+    const gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
+    const tnsr::I<DataType, 3>& coords, const Scalar<DataType>& pressure,
+    const Scalar<DataType>& sqrt_det_spatial_metric,
+    const tnsr::i<DataType, 3>& deriv_pressure) const {
   const size_t num_pts = get_size(get(pressure));
 
   for (size_t i = 0; i < num_pts; ++i) {
     const double pressure_i = get_element(get(pressure), i);
-    if (pressure_i < cutoff_pressure_) {
-      get_element(magnetic_field.get(0), i) = 0.0;
-      get_element(magnetic_field.get(1), i) = 0.0;
-      get_element(magnetic_field.get(2), i) = 0.0;
+    const double x = get_element(coords.get(0), i) - center_[0];
+    const double y = get_element(coords.get(1), i) - center_[1];
+    const double z = get_element(coords.get(2), i) - center_[2];
+    const double radius = sqrt(x * x + y * y + z * z);
+    if (pressure_i < cutoff_pressure_ or radius > max_distance_from_center_) {
       continue;
     }
 
     // (p - p_c)^{n_s}
-    const double pressure_term =
-        pow(pressure_i - cutoff_pressure_, pressure_exponent_);
+    const double pressure_term = pow(pressure_i - cutoff_pressure_,
+                                     static_cast<int>(pressure_exponent_));
     // n_s * (p - p_c)^{n_s-1}
     const double n_times_pressure_to_n_minus_1 =
-        pressure_exponent_ * pow(pressure_i - cutoff_pressure_,
-                                 static_cast<int>(pressure_exponent_) - 1);
+        static_cast<double>(pressure_exponent_) *
+        pow(pressure_i - cutoff_pressure_,
+            static_cast<int>(pressure_exponent_) - 1);
 
-    const double x = get_element(coords.get(0), i);
-    const double y = get_element(coords.get(1), i);
-
-    const auto& dp_dx = get_element(dcoords_pressure.get(0), i);
-    const auto& dp_dy = get_element(dcoords_pressure.get(1), i);
-    const auto& dp_dz = get_element(dcoords_pressure.get(2), i);
+    const auto& dp_dx = get_element(deriv_pressure.get(0), i);
+    const auto& dp_dy = get_element(deriv_pressure.get(1), i);
+    const auto& dp_dz = get_element(deriv_pressure.get(2), i);
 
     // Assign Bx, By, Bz
-    get_element(magnetic_field.get(0), i) =
-        -n_times_pressure_to_n_minus_1 * x * dp_dz;
-    get_element(magnetic_field.get(1), i) =
-        -n_times_pressure_to_n_minus_1 * y * dp_dz;
-    get_element(magnetic_field.get(2), i) =
-        2.0 * pressure_term +
-        n_times_pressure_to_n_minus_1 * (x * dp_dx + y * dp_dy);
-  }
-
-  for (size_t d = 0; d < 3; ++d) {
-    magnetic_field.get(d) *=
-        vector_potential_amplitude_ / get(sqrt_det_spatial_metric);
+    get_element(magnetic_field->get(0), i) +=
+        vector_potential_amplitude_ *
+        (-n_times_pressure_to_n_minus_1 * x * dp_dz) /
+        get_element(get(sqrt_det_spatial_metric), i);
+    get_element(magnetic_field->get(1), i) +=
+        vector_potential_amplitude_ *
+        (-n_times_pressure_to_n_minus_1 * y * dp_dz) /
+        get_element(get(sqrt_det_spatial_metric), i);
+    get_element(magnetic_field->get(2), i) +=
+        vector_potential_amplitude_ *
+        (2.0 * pressure_term +
+         n_times_pressure_to_n_minus_1 * (x * dp_dx + y * dp_dy)) /
+        get_element(get(sqrt_det_spatial_metric), i);
   }
-
-  return {std::move(magnetic_field)};
 }
 
 bool operator==(const Poloidal& lhs, const Poloidal& rhs) {
   return lhs.pressure_exponent_ == rhs.pressure_exponent_ and
          lhs.cutoff_pressure_ == rhs.cutoff_pressure_ and
-         lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_;
+         lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_ and
+         lhs.center_ == rhs.center_ and
+         lhs.max_distance_from_center_ == rhs.max_distance_from_center_;
 }
 
 bool operator!=(const Poloidal& lhs, const Poloidal& rhs) {
@@ -98,13 +129,13 @@ bool operator!=(const Poloidal& lhs, const Poloidal& rhs) {
 
 #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
 
-#define INSTANTIATE(_, data)                                               \
-  template tuples::TaggedTuple<hydro::Tags::MagneticField<DTYPE(data), 3>> \
-  Poloidal::variables<DTYPE(data)>(                                        \
-      const tnsr::I<DTYPE(data), 3>& coords,                               \
-      const Scalar<DTYPE(data)>& pressure,                                 \
-      const Scalar<DTYPE(data)>& sqrt_det_spatial_metric,                  \
-      const tnsr::i<DTYPE(data), 3>& dcoords_pressure) const;
+#define INSTANTIATE(_, data)                                  \
+  template void Poloidal::variables_impl<DTYPE(data)>(        \
+      gsl::not_null<tnsr::I<DTYPE(data), 3>*> magnetic_field, \
+      const tnsr::I<DTYPE(data), 3>& coords,                  \
+      const Scalar<DTYPE(data)>& pressure,                    \
+      const Scalar<DTYPE(data)>& sqrt_det_spatial_metric,     \
+      const tnsr::i<DTYPE(data), 3>& deriv_pressure) const;
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector))
 

src/PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Poloidal.hpp

@@ -14,6 +14,13 @@
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
 
+/// \cond
+namespace gsl {
+template <typename T>
+class not_null;
+}  // namespace gsl
+/// \endcond
+
 namespace grmhd::AnalyticData::InitialMagneticFields {
 
 /*!
@@ -63,6 +70,15 @@ namespace grmhd::AnalyticData::InitialMagneticFields {
  *        2(p-p_{\mathrm{cut}})^{n_s}.
  * \f}
  *
+ * Note that the coordinates are relative to the `Center` passed in, so the
+ * field can be centered about any arbitrary point. The field is also zero
+ * outside of `MaxDistanceFromCenter`, so that compact support can be imposed if
+ * necessary.
+ *
+ * \warning This assumes the magnetic field is initialized, both in size and
+ * value, before being passed into the `variables` function. This is so that
+ * multiple magnetic fields can be superposed. Each magnetic field
+ * configuration does a `+=` to make this possible.
  */
 class Poloidal : public InitialMagneticField {
  public:
@@ -87,8 +103,23 @@ class Poloidal : public InitialMagneticField {
     static type lower_bound() { return 0.0; }
   };
 
+  struct Center {
+    using type = std::array<double, 3>;
+    static constexpr Options::String help = {
+        "The center of the magnetic field."};
+  };
+
+  struct MaxDistanceFromCenter {
+    using type = double;
+    static constexpr Options::String help = {
+        "The maximum distance from the center to compute the magnetic field. "
+        "Everywhere outside the field is set to zero."};
+    static type lower_bound() { return 0.0; }
+  };
+
   using options =
-      tmpl::list<PressureExponent, CutoffPressure, VectorPotentialAmplitude>;
+      tmpl::list<PressureExponent, CutoffPressure, VectorPotentialAmplitude,
+                 Center, MaxDistanceFromCenter>;
 
   static constexpr Options::String help = {"Poloidal initial magnetic field"};
 
@@ -100,7 +131,8 @@ class Poloidal : public InitialMagneticField {
   ~Poloidal() override = default;
 
   Poloidal(size_t pressure_exponent, double cutoff_pressure,
-           double vector_potential_amplitude);
+           double vector_potential_amplitude, std::array<double, 3> center,
+           double max_distance_from_center);
 
   auto get_clone() const -> std::unique_ptr<InitialMagneticField> override;
 
@@ -114,18 +146,38 @@ class Poloidal : public InitialMagneticField {
   void pup(PUP::er& p) override;
 
   /// Retrieve magnetic fields at `(x)`
-  template <typename DataType>
-  auto variables(const tnsr::I<DataType, 3>& coords,
-                 const Scalar<DataType>& pressure,
-                 const Scalar<DataType>& sqrt_det_spatial_metric,
-                 const tnsr::i<DataType, 3>& dcoords_pressure) const
-      -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
+  void variables(gsl::not_null<tnsr::I<DataVector, 3>*> result,
+                 const tnsr::I<DataVector, 3>& coords,
+                 const Scalar<DataVector>& pressure,
+                 const Scalar<DataVector>& sqrt_det_spatial_metric,
+                 const tnsr::i<DataVector, 3>& deriv_pressure) const override;
+
+  /// Retrieve magnetic fields at `(x)`
+  void variables(gsl::not_null<tnsr::I<double, 3>*> result,
+                 const tnsr::I<double, 3>& coords,
+                 const Scalar<double>& pressure,
+                 const Scalar<double>& sqrt_det_spatial_metric,
+                 const tnsr::i<double, 3>& deriv_pressure) const override;
+
+  bool is_equal(const InitialMagneticField& rhs) const override;
 
  private:
+  template <typename DataType>
+  void variables_impl(gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
+                      const tnsr::I<DataType, 3>& coords,
+                      const Scalar<DataType>& pressure,
+                      const Scalar<DataType>& sqrt_det_spatial_metric,
+                      const tnsr::i<DataType, 3>& deriv_pressure) const;
+
   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::array<double, 3> center_{{std::numeric_limits<double>::signaling_NaN(),
+                                 std::numeric_limits<double>::signaling_NaN(),
+                                 std::numeric_limits<double>::signaling_NaN()}};
+  double max_distance_from_center_ =
+      std::numeric_limits<double>::signaling_NaN();
 
   friend bool operator==(const Poloidal& lhs, const Poloidal& rhs);
   friend bool operator!=(const Poloidal& lhs, const Poloidal& rhs);