Merge pull request #5948 from nilsvu/cart_to_sph

Add cartesian_to_spherical routines
sxs-collaboration · Apr 30, 2024 · 8d01120 · 8d01120
2 parents 7059ea9 + 64bf4cd
commit 8d01120
Show file tree

Hide file tree

Showing 5 changed files with 278 additions and 0 deletions.
diff --git a/src/DataStructures/Tensor/EagerMath/CMakeLists.txt b/src/DataStructures/Tensor/EagerMath/CMakeLists.txt
@@ -4,6 +4,7 @@
 spectre_target_sources(
   ${LIBRARY}
   PRIVATE
+  CartesianToSpherical.cpp
   FrameTransform.cpp
   OrthonormalOneform.cpp
   RaiseOrLowerIndex.cpp
@@ -14,6 +15,7 @@ spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
+  CartesianToSpherical.hpp
   CrossProduct.hpp
   Determinant.hpp
   DeterminantAndInverse.hpp

diff --git a/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.cpp b/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.cpp
@@ -0,0 +1,93 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "DataStructures/Tensor/EagerMath/CartesianToSpherical.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/Gsl.hpp"
+
+template <typename DataType, size_t Dim, typename CoordsFrame>
+void cartesian_to_spherical(
+    const gsl::not_null<tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>>*>
+        result,
+    const tnsr::I<DataType, Dim, CoordsFrame>& x) {
+  if constexpr (Dim == 1) {
+    get<0>(*result) = get<0>(x);
+  } else if constexpr (Dim == 2) {
+    get<0>(*result) = hypot(get<0>(x), get<1>(x));
+    get<1>(*result) = atan2(get<1>(x), get<0>(x));
+  } else if constexpr (Dim == 3) {
+    get<0>(*result) =
+        sqrt(square(get<0>(x)) + square(get<1>(x)) + square(get<2>(x)));
+    get<1>(*result) = atan2(hypot(get<0>(x), get<1>(x)), get<2>(x));
+    get<2>(*result) = atan2(get<1>(x), get<0>(x));
+  }
+}
+
+template <typename DataType, size_t Dim, typename CoordsFrame>
+tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>> cartesian_to_spherical(
+    const tnsr::I<DataType, Dim, CoordsFrame>& x) {
+  tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>> result{};
+  cartesian_to_spherical(make_not_null(&result), x);
+  return result;
+}
+
+template <typename DataType, size_t Dim, typename CoordsFrame>
+void spherical_to_cartesian(
+    const gsl::not_null<tnsr::I<DataType, Dim, CoordsFrame>*> result,
+    const tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>>& x) {
+  const auto& r = get<0>(x);
+  if constexpr (Dim == 1) {
+    get<0>(*result) = r;
+  } else if constexpr (Dim == 2) {
+    const auto& phi = get<1>(x);
+    get<0>(*result) = r * cos(phi);
+    get<1>(*result) = r * sin(phi);
+  } else if constexpr (Dim == 3) {
+    const auto& theta = get<1>(x);
+    const auto& phi = get<2>(x);
+    get<0>(*result) = r * sin(theta) * cos(phi);
+    get<1>(*result) = r * sin(theta) * sin(phi);
+    get<2>(*result) = r * cos(theta);
+  }
+}
+
+template <typename DataType, size_t Dim, typename CoordsFrame>
+tnsr::I<DataType, Dim, CoordsFrame> spherical_to_cartesian(
+    const tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>>& x) {
+  tnsr::I<DataType, Dim, CoordsFrame> result{};
+  spherical_to_cartesian(make_not_null(&result), x);
+  return result;
+}
+
+#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+#define DIM(data) BOOST_PP_TUPLE_ELEM(1, data)
+#define FRAME(data) BOOST_PP_TUPLE_ELEM(2, data)
+
+#define INSTANTIATE(_, data)                                                \
+  template void cartesian_to_spherical(                                     \
+      gsl::not_null<                                                        \
+          tnsr::I<DTYPE(data), DIM(data), Frame::Spherical<FRAME(data)>>*>  \
+          result,                                                           \
+      const tnsr::I<DTYPE(data), DIM(data), FRAME(data)>& x);               \
+  template tnsr::I<DTYPE(data), DIM(data), Frame::Spherical<FRAME(data)>>   \
+  cartesian_to_spherical(                                                   \
+      const tnsr::I<DTYPE(data), DIM(data), FRAME(data)>& x);               \
+  template void spherical_to_cartesian(                                     \
+      gsl::not_null<tnsr::I<DTYPE(data), DIM(data), FRAME(data)>*> result,  \
+      const tnsr::I<DTYPE(data), DIM(data), Frame::Spherical<FRAME(data)>>& \
+          x);                                                               \
+  template tnsr::I<DTYPE(data), DIM(data), FRAME(data)>                     \
+  spherical_to_cartesian(const tnsr::I<DTYPE(data), DIM(data),              \
+                                       Frame::Spherical<FRAME(data)>>& x);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector), (1, 2, 3),
+                        (Frame::Grid, Frame::Inertial))
+
+#undef DTYPE
+#undef DIM
+#undef FRAME
+#undef INSTANTIATE
diff --git a/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.hpp b/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.hpp
@@ -0,0 +1,39 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/Gsl.hpp"
+
+/// @{
+/*!
+ * \brief Convert between Cartesian and spherical coordinates (r, theta, phi)
+ *
+ * In 2D the order is (r, phi), and in 3D the order is (r, theta, phi), as
+ * defined by `Frame::Spherical`. The conventions for the angles are as follows:
+ *
+ * - phi is the azimuthal angle (-pi, pi], measuring the angle from the x-axis
+ * - theta is the polar angle [0, pi], measuring the angle from the z-axis
+ */
+template <typename DataType, size_t Dim, typename CoordsFrame>
+void cartesian_to_spherical(
+    gsl::not_null<tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>>*>
+        result,
+    const tnsr::I<DataType, Dim, CoordsFrame>& x);
+
+template <typename DataType, size_t Dim, typename CoordsFrame>
+tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>> cartesian_to_spherical(
+    const tnsr::I<DataType, Dim, CoordsFrame>& x);
+
+template <typename DataType, size_t Dim, typename CoordsFrame>
+void spherical_to_cartesian(
+    gsl::not_null<tnsr::I<DataType, Dim, CoordsFrame>*> result,
+    const tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>>& x);
+
+template <typename DataType, size_t Dim, typename CoordsFrame>
+tnsr::I<DataType, Dim, CoordsFrame> spherical_to_cartesian(
+    const tnsr::I<DataType, Dim, Frame::Spherical<CoordsFrame>>& x);
+/// @}
diff --git a/tests/Unit/DataStructures/Tensor/EagerMath/CMakeLists.txt b/tests/Unit/DataStructures/Tensor/EagerMath/CMakeLists.txt
@@ -4,6 +4,7 @@
 set(LIBRARY "Test_EagerMath")
 
 set(LIBRARY_SOURCES
+  Test_CartesianToSpherical.cpp
   Test_CrossProduct.cpp
   Test_Determinant.cpp
   Test_DeterminantAndInverse.cpp

diff --git a/tests/Unit/DataStructures/Tensor/EagerMath/Test_CartesianToSpherical.cpp b/tests/Unit/DataStructures/Tensor/EagerMath/Test_CartesianToSpherical.cpp
@@ -0,0 +1,143 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/EagerMath/CartesianToSpherical.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+
+SPECTRE_TEST_CASE("Unit.DataStructures.Tensor.EagerMath.CartesianToSpherical",
+                  "[DataStructures][Unit]") {
+  {
+    tnsr::I<double, 2, Frame::Grid> x{{{1., 1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(sqrt(2.0)));
+    CHECK(get<1>(x_spherical) == approx(M_PI_4));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(1.));
+    CHECK(get<1>(x_back) == approx(1.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, 1.}, DataVector{5, 0.}, DataVector{5, 0.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(1.));
+    CHECK(get<1>(x_spherical) == approx(M_PI_2));
+    CHECK(get<2>(x_spherical) == approx(0.));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(1.));
+    CHECK(get<1>(x_back) == approx(0.));
+    CHECK(get<2>(x_back) == approx(0.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, 0.}, DataVector{5, 1.}, DataVector{5, 0.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(1.));
+    CHECK(get<1>(x_spherical) == approx(M_PI_2));
+    CHECK(get<2>(x_spherical) == approx(M_PI_2));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(0.));
+    CHECK(get<1>(x_back) == approx(1.));
+    CHECK(get<2>(x_back) == approx(0.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, -1.}, DataVector{5, 0.}, DataVector{5, 0.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(1.));
+    CHECK(get<1>(x_spherical) == approx(M_PI_2));
+    CHECK(get<2>(x_spherical) == approx(M_PI));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(-1.));
+    CHECK(get<1>(x_back) == approx(0.));
+    CHECK(get<2>(x_back) == approx(0.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, 1.}, DataVector{5, 0.}, DataVector{5, 1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(sqrt(2.)));
+    CHECK(get<1>(x_spherical) == approx(M_PI_4));
+    CHECK(get<2>(x_spherical) == approx(0.));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(1.));
+    CHECK(get<1>(x_back) == approx(0.));
+    CHECK(get<2>(x_back) == approx(1.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, 1.}, DataVector{5, 0.}, DataVector{5, -1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(sqrt(2.)));
+    CHECK(get<1>(x_spherical) == approx(3. * M_PI_4));
+    CHECK(get<2>(x_spherical) == approx(0.));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(1.));
+    CHECK(get<1>(x_back) == approx(0.));
+    CHECK(get<2>(x_back) == approx(-1.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, 0.}, DataVector{5, 0.}, DataVector{5, 1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(1.));
+    CHECK(get<1>(x_spherical) == approx(0.));
+    CHECK(get<2>(x_spherical) == approx(0.));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(0.));
+    CHECK(get<1>(x_back) == approx(0.));
+    CHECK(get<2>(x_back) == approx(1.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, 0.}, DataVector{5, 0.}, DataVector{5, -1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(1.));
+    CHECK(get<1>(x_spherical) == approx(M_PI));
+    CHECK(get<2>(x_spherical) == approx(0.));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(0.));
+    CHECK(get<1>(x_back) == approx(0.));
+    CHECK(get<2>(x_back) == approx(-1.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, 1.}, DataVector{5, 1.}, DataVector{5, 1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(sqrt(3.)));
+    CHECK(get<1>(x_spherical) == approx(acos(1. / sqrt(3.))));
+    CHECK(get<2>(x_spherical) == approx(M_PI_4));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(1.));
+    CHECK(get<1>(x_back) == approx(1.));
+    CHECK(get<2>(x_back) == approx(1.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, -1.}, DataVector{5, 1.}, DataVector{5, 1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(sqrt(3.)));
+    CHECK(get<1>(x_spherical) == approx(acos(1. / sqrt(3.))));
+    CHECK(get<2>(x_spherical) == approx(3. * M_PI_4));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(-1.));
+    CHECK(get<1>(x_back) == approx(1.));
+    CHECK(get<2>(x_back) == approx(1.));
+  }
+  {
+    tnsr::I<DataVector, 3, Frame::Inertial> x{
+        {DataVector{5, -1.}, DataVector{5, -1.}, DataVector{5, 1.}}};
+    const auto x_spherical = cartesian_to_spherical(x);
+    CHECK(get<0>(x_spherical) == approx(sqrt(3.)));
+    CHECK(get<1>(x_spherical) == approx(acos(1. / sqrt(3.))));
+    CHECK(get<2>(x_spherical) == approx(-3. * M_PI_4));
+    const auto x_back = spherical_to_cartesian(x_spherical);
+    CHECK(get<0>(x_back) == approx(-1.));
+    CHECK(get<1>(x_back) == approx(-1.));
+    CHECK(get<2>(x_back) == approx(1.));
+  }
+}