Merge pull request #6413 from kidder/spherical_pfaffian_map

Add Spherical Pfaffian map

src/Domain/CoordinateMaps/CMakeLists.txt

@@ -32,6 +32,7 @@ spectre_target_sources(
   KerrHorizonConforming.cpp
   Rotation.cpp
   SpecialMobius.cpp
+  SphericalToCartesianPfaffian.cpp
   UniformCylindricalEndcap.cpp
   UniformCylindricalFlatEndcap.cpp
   UniformCylindricalSide.cpp
@@ -72,6 +73,7 @@ spectre_target_headers(
   ProductMaps.tpp
   Rotation.hpp
   SpecialMobius.hpp
+  SphericalToCartesianPfaffian.hpp
   Tags.hpp
   TimeDependentHelpers.hpp
   UniformCylindricalEndcap.hpp

src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.cpp

@@ -0,0 +1,133 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp"
+
+#include <cmath>
+#include <cstddef>
+#include <pup.h>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/DereferenceWrapper.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/MakeWithValue.hpp"
+
+namespace domain::CoordinateMaps {
+SphericalToCartesianPfaffian::SphericalToCartesianPfaffian() = default;
+SphericalToCartesianPfaffian::SphericalToCartesianPfaffian(
+    SphericalToCartesianPfaffian&&) = default;
+SphericalToCartesianPfaffian::SphericalToCartesianPfaffian(
+    const SphericalToCartesianPfaffian&) = default;
+SphericalToCartesianPfaffian& SphericalToCartesianPfaffian::operator=(
+    const SphericalToCartesianPfaffian&) = default;
+SphericalToCartesianPfaffian& SphericalToCartesianPfaffian::operator=(
+    SphericalToCartesianPfaffian&&) = default;
+
+template <typename T>
+std::array<tt::remove_cvref_wrap_t<T>, 3>
+SphericalToCartesianPfaffian::operator()(
+    const std::array<T, 3>& source_coords) const {
+  const auto& [r, theta, phi] = source_coords;
+  return {
+      {r * sin(theta) * cos(phi), r * sin(theta) * sin(phi), r * cos(theta)}};
+}
+
+// NOLINTNEXTLINE(readability-convert-member-functions-to-static)
+std::optional<std::array<double, 3>> SphericalToCartesianPfaffian::inverse(
+    const std::array<double, 3>& target_coords) const {
+  const auto& [x, y, z] = target_coords;
+  if (UNLIKELY(y == 0.0 and x == 0.0)) {
+    if (UNLIKELY(z == 0.0)) {
+      return std::array{0.0, 0.5 * M_PI, 0.0};
+    } else {
+      return std::array{std::abs(z), z > 0.0 ? 0.0 : M_PI, 0.0};
+    }
+  } else {
+    const double r = std::hypot(x, y, z);
+    return std::array{r, acos(z / r),
+                      x > 0.0 ? atan2(y, x) : atan2(y, x) + 2.0 * M_PI};
+  }
+}
+
+template <typename T>
+tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>
+SphericalToCartesianPfaffian::jacobian(
+    const std::array<T, 3>& source_coords) const {
+  const auto& [r, theta, phi] = source_coords;
+  using DataType = tt::remove_cvref_wrap_t<T>;
+  tnsr::Ij<DataType, 3, Frame::NoFrame> jacobian_matrix{
+      make_with_value<DataType>(dereference_wrapper(r), 0.0)};
+  // Pfaffian basis means phi components are 1 / sin_theta times those of a
+  // coord basis
+  const auto& cos_theta = get<2, 0>(jacobian_matrix) = cos(theta);
+  const auto& sin_theta = get<2, 1>(jacobian_matrix) = sin(theta);
+  const auto& cos_phi = get<1, 2>(jacobian_matrix) = cos(phi);
+  const auto& sin_phi = get<0, 2>(jacobian_matrix) = sin(phi);
+  get<0, 0>(jacobian_matrix) = sin_theta * cos_phi;
+  get<1, 0>(jacobian_matrix) = sin_theta * sin_phi;
+  get<0, 1>(jacobian_matrix) = r * cos_theta * cos_phi;
+  get<1, 1>(jacobian_matrix) = r * cos_theta * sin_phi;
+  get<2, 1>(jacobian_matrix) *= -r;
+  get<0, 2>(jacobian_matrix) *= -r;
+  get<1, 2>(jacobian_matrix) *= r;
+  // get<2, 2>(jacobian_matrix) is zero
+  return jacobian_matrix;
+}
+
+template <typename T>
+tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>
+SphericalToCartesianPfaffian::inv_jacobian(
+    const std::array<T, 3>& source_coords) const {
+  const auto& [r, theta, phi] = source_coords;
+  using DataType = tt::remove_cvref_wrap_t<T>;
+  tnsr::Ij<DataType, 3, Frame::NoFrame> inv_jacobian_matrix{
+      make_with_value<DataType>(dereference_wrapper(r), 0.0)};
+  // Pfaffian basis means phi components are sin_theta times those of a coord
+  // basis
+  const auto& cos_theta = get<0, 2>(inv_jacobian_matrix) = cos(theta);
+  const auto& sin_theta = get<1, 2>(inv_jacobian_matrix) = sin(theta);
+  const auto& cos_phi = get<2, 1>(inv_jacobian_matrix) = cos(phi);
+  const auto& sin_phi = get<2, 0>(inv_jacobian_matrix) = sin(phi);
+  const auto& one_over_r = get<2, 2>(inv_jacobian_matrix) = 1.0/r;
+  get<0, 0>(inv_jacobian_matrix) = sin_theta * cos_phi;
+  get<0, 1>(inv_jacobian_matrix) = sin_theta * sin_phi;
+  get<1, 0>(inv_jacobian_matrix) = cos_theta * cos_phi * one_over_r;
+  get<1, 1>(inv_jacobian_matrix) = cos_theta * sin_phi * one_over_r;
+  get<1, 2>(inv_jacobian_matrix) *= -one_over_r;
+  get<2, 0>(inv_jacobian_matrix) *= -one_over_r;
+  get<2, 1>(inv_jacobian_matrix) *= one_over_r;
+  get<2, 2>(inv_jacobian_matrix) *= 0.0;
+  return inv_jacobian_matrix;
+}
+
+void SphericalToCartesianPfaffian::pup(PUP::er& /*p*/) {}
+
+bool operator==(const SphericalToCartesianPfaffian& /*lhs*/,
+                const SphericalToCartesianPfaffian& /*rhs*/) {
+  return true;
+}
+
+bool operator!=(const SphericalToCartesianPfaffian& lhs,
+                const SphericalToCartesianPfaffian& rhs) {
+  return not(lhs == rhs);
+}
+
+#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE_DTYPE(_, data)                                           \
+  template std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>               \
+  SphericalToCartesianPfaffian::operator()(                                  \
+      const std::array<DTYPE(data), 3>& source_coords) const;                \
+  template tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame> \
+  SphericalToCartesianPfaffian::jacobian(                                    \
+      const std::array<DTYPE(data), 3>& source_coords) const;                \
+  template tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame> \
+  SphericalToCartesianPfaffian::inv_jacobian(                                \
+      const std::array<DTYPE(data), 3>& source_coords) const;
+
+GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE,
+                        (double, DataVector,
+                         std::reference_wrapper<const double>,
+                         std::reference_wrapper<const DataVector>))
+}  // namespace domain::CoordinateMaps

src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp

@@ -0,0 +1,87 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <optional>
+
+#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
+
+/// \cond
+namespace PUP {
+class er;
+}  // namespace PUP
+/// \endcond
+
+namespace domain::CoordinateMaps {
+
+/*!
+ * \ingroup CoordinateMapsGroup
+ *
+ * \brief Pfaffian transformation from spherical to Cartesian coordinates.
+ *
+ * \note This map is designed to be used together with Spherepack! Spherepack
+ * expects a Pfaffian transformation of the derivatives.
+ *
+ * \details This is a Pfaffian mapping from \f$(r,\theta,\phi) \rightarrow
+ * (x,y,z) \f$.
+ *
+ * The formula for the mapping is...
+ * \f{eqnarray*}
+ *     x &=& r \sin\theta \cos\phi \\
+ *     y &=& r \sin\theta \sin\phi \\
+ *     z &=& r \cos\theta
+ * \f}
+ *
+ * The Pfaffian basis vectors
+ * \f$ (e_{\hat r}, e_{\hat \theta}, e_{\hat \phi})\f$
+ * are related to the coordinate basis vectors
+ * \f$ (e_r, e_{\theta}, e_{\phi})\f$
+ * by...
+ * \f{eqnarray*}
+ *     e_{\hat r}      &=& e_r \\
+ *     e_{\hat \theta} &=& e_{\theta} \\
+ *     e_{\hat \phi}   &=& \frac{1}{\sin \theta} e_{\phi}
+ * \f}
+ */
+class SphericalToCartesianPfaffian {
+ public:
+  static constexpr size_t dim = 3;
+  SphericalToCartesianPfaffian();
+  ~SphericalToCartesianPfaffian() = default;
+  SphericalToCartesianPfaffian(SphericalToCartesianPfaffian&&);
+  SphericalToCartesianPfaffian(const SphericalToCartesianPfaffian&);
+  SphericalToCartesianPfaffian& operator=(const SphericalToCartesianPfaffian&);
+  SphericalToCartesianPfaffian& operator=(SphericalToCartesianPfaffian&&);
+
+  template <typename T>
+  std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
+      const std::array<T, 3>& source_coords) const;
+
+  // NOLINTNEXTLINE(readability-convert-member-functions-to-static)
+  std::optional<std::array<double, 3>> inverse(
+      const std::array<double, 3>& target_coords) const;
+
+  template <typename T>
+  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
+      const std::array<T, 3>& source_coords) const;
+
+  template <typename T>
+  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
+      const std::array<T, 3>& source_coords) const;
+
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& p);
+
+  static constexpr bool is_identity() { return false; }
+};
+
+bool operator==(const SphericalToCartesianPfaffian& lhs,
+                const SphericalToCartesianPfaffian& rhs);
+
+bool operator!=(const SphericalToCartesianPfaffian& lhs,
+                const SphericalToCartesianPfaffian& rhs);
+}  // namespace domain::CoordinateMaps

tests/Unit/Domain/CoordinateMaps/CMakeLists.txt

@@ -23,6 +23,7 @@ set(LIBRARY_SOURCES
   Test_ProductMaps.cpp
   Test_Rotation.cpp
   Test_SpecialMobius.cpp
+  Test_SphericalToCartesianPfaffian.cpp
   Test_Tags.cpp
   Test_TimeDependentHelpers.cpp
   Test_UniformCylindricalEndcap.cpp

tests/Unit/Domain/CoordinateMaps/Test_SphericalToCartesianPfaffian.cpp

@@ -0,0 +1,69 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <array>
+#include <optional>
+#include <pup.h>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp"
+#include "Framework/TestHelpers.hpp"
+#include "Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp"
+#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+
+namespace domain {
+namespace {
+void test_map_at_point(const CoordinateMaps::SphericalToCartesianPfaffian& map,
+                       const std::array<double, 3>& source_point,
+                       const std::array<double, 3>& target_point) {
+  test_inverse_map(map, source_point);
+  if (source_point !=
+      std::array{0.0, 0.5 * M_PI, 0.0}) {  // inv jac singular at origin
+    test_coordinate_map_argument_types(map, source_point);
+    test_inv_jacobian(map, source_point);
+  }
+  CHECK_ITERABLE_APPROX(map(source_point), target_point);
+  CHECK_ITERABLE_APPROX(map.inverse(target_point).value(), source_point);
+}
+
+void test_map(const CoordinateMaps::SphericalToCartesianPfaffian& map) {
+  CHECK(not map.is_identity());
+  CHECK_FALSE(map != map);
+  test_serialization(map);
+  const std::array<double, 3> source_origin{{0.0, 0.5 * M_PI, 0.0}};
+  const std::array<double, 3> target_origin{{0.0, 0.0, 0.0}};
+  test_map_at_point(map, source_origin, target_origin);
+  const std::array<double, 3> source_north_pole{{0.5, 0.0, 0.0}};
+  const std::array<double, 3> target_north_pole{{0.0, 0.0, 0.5}};
+  test_map_at_point(map, source_north_pole, target_north_pole);
+  const std::array<double, 3> source_south_pole{{0.25, M_PI, 0.0}};
+  const std::array<double, 3> target_south_pole{{0.0, 0.0, -0.25}};
+  test_map_at_point(map, source_south_pole, target_south_pole);
+  const Mesh<3> mesh{
+      {5, 3, 5},
+      {Spectral::Basis::Legendre, Spectral::Basis::SphericalHarmonic,
+       Spectral::Basis::SphericalHarmonic},
+      {Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss,
+       Spectral::Quadrature::Equiangular}};
+  const auto xi = logical_coordinates(mesh);
+  const std::array<DataVector, 3> source_coords{xi[0] + 2.0, xi[1], xi[2]};
+  test_inv_jacobian(map, source_coords);
+}
+
+void test() {
+  const CoordinateMaps::SphericalToCartesianPfaffian original_map{};
+  test_map(original_map);
+  const auto serialized_map = serialize_and_deserialize(original_map);
+  test_map(serialized_map);
+}
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.SphericalToCartesianPfaffian",
+                  "[Domain][Unit]") {
+  test();
+}
+}  // namespace domain

tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp

@@ -299,7 +299,7 @@ void test_inv_jacobian(const Map& map,
       for (size_t j = 0; j < Map::dim; ++j) {
         for (size_t l = 0; l < jacobian.get(0, 0).size(); l++) {
           for (size_t k = 0; k < Map::dim; ++k) {
-            identity.get(i, j)[k] += gsl::at(jacobian.get(i, k), l) *
+            identity.get(i, j)[l] += gsl::at(jacobian.get(i, k), l) *
                                      gsl::at(inv_jacobian.get(k, j), l);
           }
         }