Replace Kerr-Schild with Spherical-Kerr-Schild for FishboneMoncriefDisk

src/DataStructures/CachedTempBuffer.hpp

@@ -50,6 +50,7 @@ class CachedTempBuffer {
     }
     return get<Tag>(data_);
   }
+  size_t number_of_grid_points() const { return data_.number_of_grid_points(); }
 
  private:
   template <typename Tag>

src/DataStructures/TempBuffer.hpp

@@ -33,6 +33,8 @@ template <typename TagList>
 struct TempBuffer<TagList, true> : tuples::tagged_tuple_from_typelist<TagList> {
   explicit TempBuffer(const size_t /*size*/)
       : tuples::tagged_tuple_from_typelist<TagList>::TaggedTuple() {}
+
+  static size_t number_of_grid_points() { return 1; }
 };
 
 template <typename TagList>

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp

@@ -97,7 +97,7 @@ MagnetizedFmDisk::MagnetizedFmDisk(
       get<gr::Tags::SpatialMetric<DataVector, 3>>(unmagnetized_vars);
 
   const tnsr::I<double, 3> x_max{
-      {{fm_disk_.max_pressure_radius_, fm_disk_.bh_spin_a_, 0.0}}};
+      {{fm_disk_.max_pressure_radius_, 0.0, 0.0}}};
 
   const double b_squared_max = max(
       get(dot_product(b_field, b_field, spatial_metric)) /
@@ -137,11 +137,13 @@ tnsr::I<DataType, 3> MagnetizedFmDisk::unnormalized_magnetic_field(
   auto magnetic_field =
       make_with_value<tnsr::I<DataType, 3, Frame::NoFrame>>(x, 0.0);
 
+  auto x_ks = x;
+
   // The maximum pressure (and hence the maximum rest mass density) is located
-  // on the ring x^2 + y^2 = r_max^2 + a^2, z = 0.
+  // on the ring x^2 + y^2 = r_max^2, z = 0.
   // Note that `x` may or may not include points on this ring.
   const tnsr::I<double, 3> x_max{
-      {{fm_disk_.max_pressure_radius_, fm_disk_.bh_spin_a_, 0.0}}};
+      {{fm_disk_.max_pressure_radius_, 0.0, 0.0}}};
   const double threshold_rest_mass_density =
       threshold_density_ *
       get(get<hydro::Tags::RestMassDensity<double>>(variables(
@@ -174,13 +176,10 @@ tnsr::I<DataType, 3> MagnetizedFmDisk::unnormalized_magnetic_field(
       // (r, theta, phi) coordinates. We need to get B^r, B^theta, B^phi first,
       // then we transform to Cartesian coordinates.
       const double z_squared = square(get_element(get<2>(x), s));
-      const double a_squared = fm_disk_.bh_spin_a_ * fm_disk_.bh_spin_a_;
-
       double sin_theta_squared =
           square(get_element(get<0>(x), s)) + square(get_element(get<1>(x), s));
-      double r_squared = 0.5 * (sin_theta_squared + z_squared - a_squared);
-      r_squared += sqrt(square(r_squared) + a_squared * z_squared);
-      sin_theta_squared /= (r_squared + a_squared);
+      double r_squared = sin_theta_squared + z_squared;
+      sin_theta_squared /= r_squared;
 
       // B^i is proportional to derivatives of the magnetic potential. One has
       //
@@ -216,20 +215,53 @@ tnsr::I<DataType, 3> MagnetizedFmDisk::unnormalized_magnetic_field(
           prefactor * (mag_potential(radius - small, sin_theta_squared) -
                        mag_potential(radius + small, sin_theta_squared));
       // phi component of the field vanishes identically.
+
+      // Need x in KS coordinates instead of SKS coordinates
+      // to do transformation below. copy x into x_sks and alter at each
+      // element index s.
+      const double sks_to_ks_factor =
+          sqrt(r_squared + square(fm_disk_.bh_spin_a_)) / radius;
+      get_element(x_ks.get(0), s) = get_element(x.get(0), s) * sks_to_ks_factor;
+      get_element(x_ks.get(1), s) = get_element(x.get(1), s) * sks_to_ks_factor;
+      get_element(x_ks.get(2), s) = get_element(x.get(2), s);
+    }
+  }
+  // magnetic field in KS_coords
+  // change back to SKS_coords
+  const auto magnetic_field_ks =
+      kerr_schild_coords_.cartesian_from_spherical_ks(std::move(magnetic_field),
+                                                      x_ks);
+
+  using inv_jacobian =
+      gr::Solutions::SphericalKerrSchild::internal_tags::inv_jacobian<
+          DataType, Frame::Inertial>;
+
+  FmDisk::IntermediateVariables<DataType> vars(x);
+
+  const auto inv_jacobians =
+      get<inv_jacobian>(fm_disk_.background_spacetime_.variables(
+          x, 0.0, tmpl::list<inv_jacobian>{},
+          make_not_null(&vars.sph_kerr_schild_cache)));
+
+  auto magnetic_field_sks =
+      make_with_value<tnsr::I<DataType, 3, Frame::Inertial>>(x, 0.0);
+
+  for (size_t j = 0; j < 3; ++j) {
+    for (size_t i = 0; i < 3; ++i) {
+      magnetic_field_sks.get(j) +=
+          inv_jacobians.get(j, i) * magnetic_field_ks.get(i);
     }
   }
-  return kerr_schild_coords_.cartesian_from_spherical_ks(
-      std::move(magnetic_field), x);
+
+  return magnetic_field_sks;
 }
 
-template <typename DataType, bool NeedSpacetime>
+template <typename DataType>
 tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>
 MagnetizedFmDisk::variables(
     const tnsr::I<DataType, 3>& x,
     tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
-    const typename FmDisk::IntermediateVariables<DataType,
-                                                 NeedSpacetime>& /*vars*/,
-    const size_t /*index*/) const {
+    gsl::not_null<FmDisk::IntermediateVariables<DataType>*> /*vars*/) const {
   auto result = unnormalized_magnetic_field(x);
   for (size_t i = 0; i < 3; ++i) {
     result.get(i) *= b_field_normalization_;
@@ -252,20 +284,18 @@ bool operator!=(const MagnetizedFmDisk& lhs, const MagnetizedFmDisk& rhs) {
 }
 
 #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
-#define NEED_SPACETIME(data) BOOST_PP_TUPLE_ELEM(1, data)
 
-#define INSTANTIATE(_, data)                                               \
-  template tuples::TaggedTuple<hydro::Tags::MagneticField<DTYPE(data), 3>> \
-  MagnetizedFmDisk::variables(                                             \
-      const tnsr::I<DTYPE(data), 3>& x,                                    \
-      tmpl::list<hydro::Tags::MagneticField<DTYPE(data), 3>> /*meta*/,     \
-      const typename FmDisk::FishboneMoncriefDisk::IntermediateVariables<  \
-          DTYPE(data), NEED_SPACETIME(data)>& vars,                        \
-      const size_t) const;
+#define INSTANTIATE(_, data)                                                 \
+  template tuples::TaggedTuple<hydro::Tags::MagneticField<DTYPE(data), 3>>   \
+  MagnetizedFmDisk::variables(                                               \
+      const tnsr::I<DTYPE(data), 3>& x,                                      \
+      tmpl::list<hydro::Tags::MagneticField<DTYPE(data), 3>> /*meta*/,       \
+      gsl::not_null<                                                         \
+          FmDisk::FishboneMoncriefDisk::IntermediateVariables<DTYPE(data)>*> \
+          vars) const;
 
-GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector), (true, false))
+GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector))
 
 #undef DTYPE
-#undef NEED_SPACETIME
 #undef INSTANTIATE
 }  // namespace grmhd::AnalyticData

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp

@@ -70,6 +70,21 @@ namespace grmhd::AnalyticData {
  *
  * is the norm of the magnetic field in the fluid frame, with \f$v_i\f$ being
  * the spatial velocity, and \f$W\f$ the Lorentz factor.
+ *
+ * \note When using Kerr-Schild coordinates, the horizon that is at
+ * constant \f$r\f$ is not spherical, but instead spheroidal. This could make
+ * application of boundary condition and computing various fluxes
+ * across the horizon more complicated than they need to be.
+ * Thus, similar to RelativisticEuler::Solutions::FishboneMoncriefDisk
+ * we use Spherical Kerr-Schild coordinates,
+ * see gr::Solutions::SphericalKerrSchild, in which constant \f$r\f$
+ * is spherical. Because we compute variables in Kerr-Schild coordinates,
+ * there is a necessary extra step of transforming them back to
+ * Spherical Kerr-Schild coordinates.
+ *
+ * \warning Spherical Kerr-Schild coordinates and "spherical Kerr-Schild"
+ * coordinates are not same.
+ *
  */
 class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData,
                          public MarkAsAnalyticData {
@@ -153,38 +168,19 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData,
   using equation_of_state_type = typename FmDisk::equation_of_state_type;
 
   /// @{
-  /// The grmhd variables in Cartesian Kerr-Schild coordinates at `(x, t)`
-  ///
-  /// \note The functions are optimized for retrieving the hydro variables
-  /// before the metric variables.
+  /// The variables in Cartesian Kerr-Schild coordinates at `(x, t)`.
   template <typename DataType, typename... Tags>
   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
                                          tmpl::list<Tags...> /*meta*/) const {
     // Can't store IntermediateVariables as member variable because we
     // need to be threadsafe.
-    constexpr double dummy_time = 0.0;
-    typename FmDisk::IntermediateVariables<
-        DataType,
-        tmpl2::flat_any_v<(
-            std::is_same_v<Tags, hydro::Tags::SpatialVelocity<DataType, 3>> or
-            std::is_same_v<Tags, hydro::Tags::LorentzFactor<DataType>> or
-            not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tags>)...>>
-        vars(fm_disk_.bh_spin_a_, fm_disk_.background_spacetime_, x, dummy_time,
-             FmDisk::index_helper(
-                 tmpl::index_of<tmpl::list<Tags...>,
-                                hydro::Tags::SpatialVelocity<DataType, 3>>{}),
-             FmDisk::index_helper(
-                 tmpl::index_of<tmpl::list<Tags...>,
-                                hydro::Tags::LorentzFactor<DataType>>{}));
+    typename FmDisk::IntermediateVariables<DataType> vars(x);
     return {std::move(get<Tags>([this, &x, &vars]() {
       if constexpr (std::is_same_v<hydro::Tags::MagneticField<DataType, 3>,
                                    Tags>) {
-        return variables(x, tmpl::list<Tags>{}, vars,
-                         tmpl::index_of<tmpl::list<Tags...>, Tags>::value);
+        return variables(x, tmpl::list<Tags>{}, make_not_null(&vars));
       } else {
-        return fm_disk_.variables(
-            x, tmpl::list<Tags>{}, vars,
-            tmpl::index_of<tmpl::list<Tags...>, Tags>::value);
+        return fm_disk_.variables(x, tmpl::list<Tags>{}, make_not_null(&vars));
       }
     }()))...};
   }
@@ -194,20 +190,12 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData,
                                      tmpl::list<Tag> /*meta*/) const {
     // Can't store IntermediateVariables as member variable because we need to
     // be threadsafe.
-    constexpr double dummy_time = 0.0;
-    typename FmDisk::IntermediateVariables<
-        DataType,
-        std::is_same_v<Tag, hydro::Tags::SpatialVelocity<DataType, 3>> or
-            std::is_same_v<Tag, hydro::Tags::LorentzFactor<DataType>> or
-            not tmpl::list_contains_v<hydro::grmhd_tags<DataType>, Tag>>
-        intermediate_vars(fm_disk_.bh_spin_a_, fm_disk_.background_spacetime_,
-                          x, dummy_time, std::numeric_limits<size_t>::max(),
-                          std::numeric_limits<size_t>::max());
+    typename FmDisk::IntermediateVariables<DataType> vars(x);
     if constexpr (std::is_same_v<hydro::Tags::MagneticField<DataType, 3>,
                                  Tag>) {
-      return variables(x, tmpl::list<Tag>{}, intermediate_vars, 0);
+      return variables(x, tmpl::list<Tag>{}, make_not_null(&vars));
     } else {
-      return fm_disk_.variables(x, tmpl::list<Tag>{}, intermediate_vars, 0);
+      return fm_disk_.variables(x, tmpl::list<Tag>{}, make_not_null(&vars));
     }
   }
   /// @}
@@ -220,14 +208,11 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData,
   }
 
  private:
-  template <typename DataType, bool NeedSpacetime>
-  auto variables(
-      const tnsr::I<DataType, 3>& x,
-      tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
-      const typename FmDisk::IntermediateVariables<DataType, NeedSpacetime>&
-          vars,
-      size_t index) const
-      -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
+  template <typename DataType>
+  auto variables(const tnsr::I<DataType, 3>& x,
+                 tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
+                 gsl::not_null<FmDisk::IntermediateVariables<DataType>*> vars)
+      const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
 
   template <typename DataType>
   tnsr::I<DataType, 3> unnormalized_magnetic_field(

src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp

@@ -19,6 +19,7 @@
 #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/ForceInline.hpp"
+#include "Utilities/Gsl.hpp"
 #include "Utilities/MakeArray.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
@@ -93,7 +94,7 @@ namespace Solutions {
  * where we have defined
  *
  * \f{align}{
- *   \rho \equiv r^2 + a^2. \label{eq:rho}
+ *   \rho^2 \equiv r^2 + a^2. \label{eq:rho}
  * \f}
  *
  * Therefore, we have that \f$r\f$ satisfies the equation for a sphere in
@@ -429,20 +430,6 @@ class SphericalKerrSchild : public AnalyticSolution<3_st>,
   SphericalKerrSchild& operator=(SphericalKerrSchild&& /*rhs*/) = default;
   ~SphericalKerrSchild() = default;
 
-  template <typename DataType, typename Frame, typename... Tags>
-  tuples::TaggedTuple<Tags...> variables(
-      const tnsr::I<DataType, volume_dim, Frame>& x, double /*t*/,
-      tmpl::list<Tags...> /*meta*/) const {
-    static_assert(
-        tmpl2::flat_all_v<
-            tmpl::list_contains_v<tags<DataType, Frame>, Tags>...>,
-        "At least one of the requested tags is not supported. The requested "
-        "tags are listed as template parameters of the `variables` function.");
-    IntermediateVars<DataType, Frame> cache(get_size(*x.begin()));
-    IntermediateComputer<DataType, Frame> computer(*this, x);
-    return {cache.get_var(computer, Tags{})...};
-  }
-
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p);
 
@@ -568,6 +555,46 @@ class SphericalKerrSchild : public AnalyticSolution<3_st>,
       gr::Tags::SpatialChristoffelFirstKind<DataType, 3, Frame>,
       gr::Tags::SpatialChristoffelSecondKind<DataType, 3, Frame>>;
 
+  // forward-declaration needed.
+  template <typename DataType, typename Frame>
+  class IntermediateVars;
+
+  template <typename DataType, typename Frame = Frame::Inertial>
+  using allowed_tags =
+      tmpl::push_back<tags<DataType, Frame>,
+                      typename internal_tags::inv_jacobian<DataType, Frame>>;
+
+  template <typename DataType, typename Frame, typename... Tags>
+  tuples::TaggedTuple<Tags...> variables(
+      const tnsr::I<DataType, volume_dim, Frame>& x, double /*t*/,
+      tmpl::list<Tags...> /*meta*/) const {
+    static_assert(
+        tmpl2::flat_all_v<
+            tmpl::list_contains_v<allowed_tags<DataType, Frame>, Tags>...>,
+        "At least one of the requested tags is not supported. The requested "
+        "tags are listed as template parameters of the `variables` function.");
+    IntermediateVars<DataType, Frame> cache(get_size(*x.begin()));
+    IntermediateComputer<DataType, Frame> computer(*this, x);
+    return {cache.get_var(computer, Tags{})...};
+  }
+
+  template <typename DataType, typename Frame, typename... Tags>
+  tuples::TaggedTuple<Tags...> variables(
+      const tnsr::I<DataType, volume_dim, Frame>& x, double /*t*/,
+      tmpl::list<Tags...> /*meta*/,
+      gsl::not_null<IntermediateVars<DataType, Frame>*> cache) const {
+    static_assert(
+        tmpl2::flat_all_v<
+            tmpl::list_contains_v<allowed_tags<DataType, Frame>, Tags>...>,
+        "At least one of the requested tags is not supported. The requested "
+        "tags are listed as template parameters of the `variables` function.");
+    if (cache->number_of_grid_points() != get_size(*x.begin())) {
+      *cache = IntermediateVars<DataType, Frame>(get_size(*x.begin()));
+    }
+    IntermediateComputer<DataType, Frame> computer(*this, x);
+    return {cache->get_var(computer, Tags{})...};
+  }
+
   template <typename DataType, typename Frame = ::Frame::Inertial>
   class IntermediateComputer {
    public: