add PolarMagnetizedFMDisk

src/Evolution/Executables/GrMhd/ValenciaDivClean/CMakeLists.txt

@@ -130,6 +130,12 @@ add_grmhd_executable(
   "${LIBS_TO_LINK}"
   )
 
+add_grmhd_executable(
+  PolarMagnetizedFmDisk
+  grmhd::AnalyticData::PolarMagnetizedFmDisk
+  "${LIBS_TO_LINK}"
+  )
+
 add_grmhd_executable(
   RiemannProblem
   grmhd::AnalyticData::RiemannProblem

src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp

@@ -146,6 +146,7 @@
 #include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/OrszagTangVortex.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/RiemannProblem.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/SlabJet.hpp"
 #include "PointwiseFunctions/AnalyticData/Tags.hpp"

src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivCleanFwd.hpp

@@ -28,6 +28,7 @@ class MagneticFieldLoop;
 class MagneticRotor;
 class MagnetizedFmDisk;
 class MagnetizedTovStar;
+class PolarMagnetizedFmDisk;
 class OrszagTangVortex;
 class RiemannProblem;
 class SlabJet;

src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt

@@ -17,6 +17,7 @@ spectre_target_sources(
   MagnetizedFmDisk.cpp
   MagnetizedTovStar.cpp
   OrszagTangVortex.cpp
+  PolarMagnetizedFmDisk.cpp
   RiemannProblem.cpp
   SlabJet.cpp
   )
@@ -35,6 +36,7 @@ spectre_target_headers(
   MagnetizedFmDisk.hpp
   MagnetizedTovStar.hpp
   OrszagTangVortex.hpp
+  PolarMagnetizedFmDisk.hpp
   RiemannProblem.hpp
   SlabJet.hpp
   )

src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.cpp

@@ -0,0 +1,36 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp"
+
+#include <pup.h>
+#include <utility>
+
+#include "Domain/CoordinateMaps/CoordinateMap.tpp"
+
+namespace grmhd::AnalyticData {
+
+PolarMagnetizedFmDisk::PolarMagnetizedFmDisk(
+    MagnetizedFmDisk fm_disk, domain::CoordinateMaps::SphericalTorus torus_map)
+    : fm_disk_(std::move(fm_disk)),
+      torus_map_(
+          domain::make_coordinate_map<Frame::BlockLogical, Frame::Inertial>(
+              torus_map)),
+      bare_torus_map_(std::move(torus_map)) {}
+
+void PolarMagnetizedFmDisk::pup(PUP::er& p) {
+  p | fm_disk_;
+  p | torus_map_;
+  p | bare_torus_map_;
+}
+
+bool operator==(const PolarMagnetizedFmDisk& lhs,
+                const PolarMagnetizedFmDisk& rhs) {
+  return lhs.fm_disk_ == rhs.fm_disk_ and lhs.torus_map_ == rhs.torus_map_;
+}
+
+bool operator!=(const PolarMagnetizedFmDisk& lhs,
+                const PolarMagnetizedFmDisk& rhs) {
+  return not(lhs == rhs);
+}
+}  // namespace grmhd::AnalyticData

src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp

@@ -0,0 +1,166 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <iterator>
+#include <type_traits>
+
+#include "DataStructures/Tensor/EagerMath/Determinant.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/CoordinateMaps/CoordinateMap.hpp"
+#include "Domain/CoordinateMaps/SphericalTorus.hpp"
+#include "Options/Options.hpp"
+#include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Transform.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/TMPL.hpp"
+#include "Utilities/TaggedTuple.hpp"
+
+/// \cond
+namespace PUP {
+class er;
+}  // namespace PUP
+/// \endcond
+
+namespace grmhd::AnalyticData {
+/*!
+ * \brief Magnetized fluid disk orbiting a Kerr black hole in Polar coordinates.
+ * FIXME (add more)
+ */
+class PolarMagnetizedFmDisk : public MarkAsAnalyticData,
+                              public grmhd::AnalyticDataBase {
+ public:
+  struct DiskParameters {
+    using type = MagnetizedFmDisk;
+    static constexpr Options::String help = "Parameters for the disk.";
+  };
+
+  struct TorusParameters {
+    using type = domain::CoordinateMaps::SphericalTorus;
+    static constexpr Options::String help =
+        "Parameters for the evolution region.";
+  };
+
+  using options = tmpl::list<DiskParameters, TorusParameters>;
+
+  static constexpr Options::String help =
+      "Magnetized Fishbone-Moncrief disk in polar coordinates.";
+
+  PolarMagnetizedFmDisk() = default;
+
+  PolarMagnetizedFmDisk(MagnetizedFmDisk fm_disk,
+                        domain::CoordinateMaps::SphericalTorus torus_map);
+
+  /// The grmhd variables.
+  ///
+  /// \note The functions are optimized for retrieving the hydro variables
+  /// before the metric variables.
+  template <typename DataType, typename... Tags>
+  tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
+                                         tmpl::list<Tags...> /*meta*/) const {
+    // In this function, we label the coordinates this solution works
+    // in with Frame::BlockLogical, and the coordinates the wrapped
+    // solution uses Frame::Inertial.  This means the input and output
+    // have to be converted to the correct label.
+
+    // We have to copy here to match the CoordinateMap interface.
+    tnsr::I<DataType, 3, Frame::BlockLogical> x_logical;
+    std::array<DataType, 3> x_array;
+    for (size_t i = 0; i < 3; ++i) {
+      x_logical.get(i) = x.get(i);
+      gsl::at(x_array, i) = x.get(i);
+    }
+
+    const tnsr::I<DataType, 3> observation_coordinates(torus_map_(x_logical));
+
+    using dependencies = tmpl::map<
+        tmpl::pair<gr::AnalyticSolution<3>::DerivShift<DataType>,
+                   gr::Tags::Shift<DataType, 3, Frame::Inertial>>,
+        tmpl::pair<gr::AnalyticSolution<3>::DerivSpatialMetric<DataType>,
+                   gr::Tags::SpatialMetric<DataType, 3, Frame::Inertial>>>;
+    using required_tags = tmpl::remove_duplicates<
+        tmpl::remove<tmpl::list<Tags..., tmpl::at<dependencies, Tags>...>,
+                     tmpl::no_such_type_>>;
+
+    auto observation_data =
+        fm_disk_.variables(observation_coordinates, required_tags{});
+
+    const auto jacobian = torus_map_.jacobian(x_logical);
+    const auto inv_jacobian = torus_map_.inv_jacobian(x_logical);
+
+    const auto change_frame = [this, &inv_jacobian, &jacobian, &x_array](
+                                  const auto& data, auto tag) {
+      using Tag = decltype(tag);
+      auto result =
+          transform::to_different_frame(get<Tag>(data), jacobian, inv_jacobian);
+
+      if constexpr (std::is_same_v<
+                        Tag, gr::AnalyticSolution<3>::DerivShift<DataType>>) {
+        const auto deriv_inv_jacobian =
+            bare_torus_map_.derivative_of_inv_jacobian(x_array);
+        const auto& shift =
+            get<gr::Tags::Shift<DataType, 3, Frame::Inertial>>(data);
+        auto extra_term = tenex.evaluate<ti::i,ti::J>(
+          deriv_inv_jacobian(ti::J,ti::k,ti::i) * shift(ti::K));
+        auto result = tenex.evaluate<ti::i,ti::J>(result(ti::i,ti::J)+
+          extra_term(ti::i,ti::J));
+      } else if constexpr (std::is_same_v<
+                               Tag, gr::AnalyticSolution<3>::DerivSpatialMetric<
+                                        DataType>>) {
+        const auto hessian = bare_torus_map_.hessian(x_array);
+        const auto& spatial_metric =
+            get<gr::Tags::SpatialMetric<DataType, 3, Frame::Inertial>>(data);
+        auto extra_term_1 = tenex.evaluate<ti::i;ti::j,ti::k>(
+          hessian(ti::L,ti::j,ti::i) * jacobian(ti::M,ti::k) *
+          spatial_metric(ti::l,ti::m));
+        auto extra_term_2 = tenex.evaluate<ti::l,ti::k,ti::i>(
+          hessian(ti::L,ti::k,ti::i) * jacobian(ti::M,ti::j) *
+          spatial_metric(ti::l,ti::m));
+        auto result = tenex.evaluate<ti::i,ti::j,ti::k>(
+          result(ti::i,ti::j,ti::k) + extra_term_1(ti::i,ti::j,ti::k) +
+          extra_term_2(ti::i,ti::j,ti::k));
+      } else if constexpr (std::is_same_v<
+                               Tag, gr::Tags::SqrtDetSpatialMetric<DataType>>) {
+        get(result) *= abs(get(determinant(jacobian)));
+      }
+
+      typename Tag::type result_with_replaced_frame{};
+      std::copy(std::move_iterator(result.begin()),
+                std::move_iterator(result.end()),
+                result_with_replaced_frame.begin());
+      return result_with_replaced_frame;
+    };
+
+    return {change_frame(observation_data, Tags{})...};
+  }
+
+  using equation_of_state_type = MagnetizedFmDisk::equation_of_state_type;
+  const equation_of_state_type& equation_of_state() const {
+    return fm_disk_.equation_of_state();
+  }
+
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& p);
+
+ private:
+  friend bool operator==(const PolarMagnetizedFmDisk& lhs,
+                         const PolarMagnetizedFmDisk& rhs);
+
+  MagnetizedFmDisk fm_disk_;
+  domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
+                        domain::CoordinateMaps::SphericalTorus>
+      torus_map_;
+  domain::CoordinateMaps::SphericalTorus bare_torus_map_;
+};
+
+bool operator!=(const PolarMagnetizedFmDisk& lhs,
+                const PolarMagnetizedFmDisk& rhs);
+}  // namespace grmhd::AnalyticData