diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt index 6366d317078b..340cc7d147c2 100644 --- a/src/Domain/CoordinateMaps/CMakeLists.txt +++ b/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 diff --git a/src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.cpp b/src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.cpp new file mode 100644 index 000000000000..ed8edf363c6b --- /dev/null +++ b/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 +#include +#include + +#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 +std::array, 3> +SphericalToCartesianPfaffian::operator()( + const std::array& 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> SphericalToCartesianPfaffian::inverse( + const std::array& 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 +tnsr::Ij, 3, Frame::NoFrame> +SphericalToCartesianPfaffian::jacobian( + const std::array& source_coords) const { + const auto& [r, theta, phi] = source_coords; + using DataType = tt::remove_cvref_wrap_t; + tnsr::Ij jacobian_matrix{ + make_with_value(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 +tnsr::Ij, 3, Frame::NoFrame> +SphericalToCartesianPfaffian::inv_jacobian( + const std::array& source_coords) const { + const auto& [r, theta, phi] = source_coords; + using DataType = tt::remove_cvref_wrap_t; + tnsr::Ij inv_jacobian_matrix{ + make_with_value(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, 3> \ + SphericalToCartesianPfaffian::operator()( \ + const std::array& source_coords) const; \ + template tnsr::Ij, 3, Frame::NoFrame> \ + SphericalToCartesianPfaffian::jacobian( \ + const std::array& source_coords) const; \ + template tnsr::Ij, 3, Frame::NoFrame> \ + SphericalToCartesianPfaffian::inv_jacobian( \ + const std::array& source_coords) const; + +GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE, + (double, DataVector, + std::reference_wrapper, + std::reference_wrapper)) +} // namespace domain::CoordinateMaps diff --git a/src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp b/src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp new file mode 100644 index 000000000000..77a3ae666e23 --- /dev/null +++ b/src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp @@ -0,0 +1,87 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include + +#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 + std::array, 3> operator()( + const std::array& source_coords) const; + + // NOLINTNEXTLINE(readability-convert-member-functions-to-static) + std::optional> inverse( + const std::array& target_coords) const; + + template + tnsr::Ij, 3, Frame::NoFrame> jacobian( + const std::array& source_coords) const; + + template + tnsr::Ij, 3, Frame::NoFrame> inv_jacobian( + const std::array& 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 diff --git a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt index 39c8e39eb1b6..461feb33a612 100644 --- a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt +++ b/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 diff --git a/tests/Unit/Domain/CoordinateMaps/Test_SphericalToCartesianPfaffian.cpp b/tests/Unit/Domain/CoordinateMaps/Test_SphericalToCartesianPfaffian.cpp new file mode 100644 index 000000000000..07492a9b6ea9 --- /dev/null +++ b/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 +#include +#include + +#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& source_point, + const std::array& 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 source_origin{{0.0, 0.5 * M_PI, 0.0}}; + const std::array target_origin{{0.0, 0.0, 0.0}}; + test_map_at_point(map, source_origin, target_origin); + const std::array source_north_pole{{0.5, 0.0, 0.0}}; + const std::array target_north_pole{{0.0, 0.0, 0.5}}; + test_map_at_point(map, source_north_pole, target_north_pole); + const std::array source_south_pole{{0.25, M_PI, 0.0}}; + const std::array 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 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