From a4e454e18bd1d708ecd59cf203bde5d7b06825c3 Mon Sep 17 00:00:00 2001 From: Jooheon Yoo Date: Thu, 5 Oct 2023 07:19:33 -0700 Subject: [PATCH] Move SphericalTorus to PointwiseFunctions/AnalyticData/GrMhd --- src/Domain/CoordinateMaps/CMakeLists.txt | 2 - .../AnalyticData/GrMhd/CMakeLists.txt | 2 + .../GrMhd/PolarMagnetizedFmDisk.cpp | 10 +- .../GrMhd/PolarMagnetizedFmDisk.hpp | 34 ++-- .../AnalyticData/GrMhd}/SphericalTorus.cpp | 147 ++++++++------- .../AnalyticData/GrMhd}/SphericalTorus.hpp | 37 ++-- .../Unit/Domain/CoordinateMaps/CMakeLists.txt | 1 - .../CoordinateMaps/Test_SphericalTorus.cpp | 125 ------------- .../AnalyticData/GrMhd/CMakeLists.txt | 1 + .../GrMhd/Test_PolarMagnetizedFmDisk.cpp | 33 ++-- .../GrMhd/Test_SphericalTorus.cpp | 171 ++++++++++++++++++ 11 files changed, 292 insertions(+), 271 deletions(-) rename src/{Domain/CoordinateMaps => PointwiseFunctions/AnalyticData/GrMhd}/SphericalTorus.cpp (71%) rename src/{Domain/CoordinateMaps => PointwiseFunctions/AnalyticData/GrMhd}/SphericalTorus.hpp (76%) delete mode 100644 tests/Unit/Domain/CoordinateMaps/Test_SphericalTorus.cpp create mode 100644 tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_SphericalTorus.cpp diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt index 3f670d8a67782..d1cb30c1ac79a 100644 --- a/src/Domain/CoordinateMaps/CMakeLists.txt +++ b/src/Domain/CoordinateMaps/CMakeLists.txt @@ -32,7 +32,6 @@ spectre_target_sources( KerrHorizonConforming.cpp Rotation.cpp SpecialMobius.cpp - SphericalTorus.cpp UniformCylindricalEndcap.cpp UniformCylindricalFlatEndcap.cpp UniformCylindricalSide.cpp @@ -73,7 +72,6 @@ spectre_target_headers( ProductMaps.tpp Rotation.hpp SpecialMobius.hpp - SphericalTorus.hpp Tags.hpp TimeDependentHelpers.hpp UniformCylindricalEndcap.hpp diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt b/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt index 7fa226f63fcba..53b1bc59d56ab 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt @@ -20,6 +20,7 @@ spectre_target_sources( PolarMagnetizedFmDisk.cpp RiemannProblem.cpp SlabJet.cpp + SphericalTorus.cpp ) spectre_target_headers( @@ -39,6 +40,7 @@ spectre_target_headers( PolarMagnetizedFmDisk.hpp RiemannProblem.hpp SlabJet.hpp + SphericalTorus.hpp ) target_link_libraries( diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.cpp index 987a9eedda7f6..740085b11c58a 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.cpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.cpp @@ -6,17 +6,12 @@ #include #include -#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( - torus_map)), - bare_torus_map_(std::move(torus_map)) {} + MagnetizedFmDisk fm_disk, grmhd::AnalyticData::SphericalTorus torus_map) + : fm_disk_(std::move(fm_disk)), torus_map_(std::move(torus_map)) {} std::unique_ptr PolarMagnetizedFmDisk::get_clone() const { @@ -29,7 +24,6 @@ PolarMagnetizedFmDisk::PolarMagnetizedFmDisk(CkMigrateMessage* msg) void PolarMagnetizedFmDisk::pup(PUP::er& p) { p | fm_disk_; p | torus_map_; - p | bare_torus_map_; } PUP::able::PUP_ID PolarMagnetizedFmDisk::my_PUP_ID = 0; // NOLINT diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp index 1ea466161f960..813cf367f5408 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp @@ -12,12 +12,11 @@ #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/AnalyticData/GrMhd/SphericalTorus.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp" #include "PointwiseFunctions/GeneralRelativity/Transform.hpp" #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" @@ -36,7 +35,7 @@ namespace grmhd::AnalyticData { /*! * \brief Magnetized fluid disk orbiting a Kerr black hole in the usual * Cartesian coordinates, but in a toroidal mesh defined from a torus map - * (see doamin::CoordinateMaps::SphericalTorus). + * (see grmhd::AnalyticData::SphericalTorus). */ class PolarMagnetizedFmDisk : public virtual evolution::initial_data::InitialData, @@ -49,7 +48,7 @@ class PolarMagnetizedFmDisk }; struct TorusParameters { - using type = domain::CoordinateMaps::SphericalTorus; + using type = grmhd::AnalyticData::SphericalTorus; static constexpr Options::String help = "Parameters for the evolution region."; }; @@ -68,7 +67,7 @@ class PolarMagnetizedFmDisk ~PolarMagnetizedFmDisk() override = default; PolarMagnetizedFmDisk(MagnetizedFmDisk fm_disk, - domain::CoordinateMaps::SphericalTorus torus_map); + grmhd::AnalyticData::SphericalTorus torus_map); auto get_clone() const -> std::unique_ptr override; @@ -91,15 +90,7 @@ class PolarMagnetizedFmDisk // 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 x_logical; - std::array 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 observation_coordinates(torus_map_(x_logical)); + const tnsr::I observation_coordinates(torus_map_(x)); using dependencies = tmpl::map< tmpl::pair::DerivShift, @@ -113,10 +104,10 @@ class PolarMagnetizedFmDisk 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 jacobian = torus_map_.jacobian(x); + const auto inv_jacobian = torus_map_.inv_jacobian(x); - const auto change_frame = [this, &inv_jacobian, &jacobian, &x_array]( + const auto change_frame = [this, &inv_jacobian, &jacobian, &x]( const auto& data, auto tag) { using Tag = decltype(tag); auto result = @@ -125,7 +116,7 @@ class PolarMagnetizedFmDisk if constexpr (std::is_same_v< Tag, gr::AnalyticSolution<3>::DerivShift>) { const auto deriv_inv_jacobian = - bare_torus_map_.derivative_of_inv_jacobian(x_array); + torus_map_.derivative_of_inv_jacobian(x); const auto& shift = get>(data); for (size_t i = 0; i < 3; ++i) { @@ -139,7 +130,7 @@ class PolarMagnetizedFmDisk } else if constexpr (std::is_same_v< Tag, gr::AnalyticSolution<3>::DerivSpatialMetric< DataType>>) { - const auto hessian = bare_torus_map_.hessian(x_array); + const auto hessian = torus_map_.hessian(x); const auto& spatial_metric = get>(data); for (size_t i = 0; i < 3; ++i) { @@ -184,10 +175,7 @@ class PolarMagnetizedFmDisk const PolarMagnetizedFmDisk& rhs); MagnetizedFmDisk fm_disk_; - domain::CoordinateMap - torus_map_; - domain::CoordinateMaps::SphericalTorus bare_torus_map_; + grmhd::AnalyticData::SphericalTorus torus_map_; }; bool operator!=(const PolarMagnetizedFmDisk& lhs, diff --git a/src/Domain/CoordinateMaps/SphericalTorus.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp similarity index 71% rename from src/Domain/CoordinateMaps/SphericalTorus.cpp rename to src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp index 8b1541cddb817..c9c999b250131 100644 --- a/src/Domain/CoordinateMaps/SphericalTorus.cpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp @@ -1,18 +1,17 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "Domain/CoordinateMaps/SphericalTorus.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp" #include #include "Options/ParseError.hpp" -#include "Utilities/DereferenceWrapper.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/MakeWithValue.hpp" -namespace domain::CoordinateMaps { +namespace grmhd::AnalyticData { SphericalTorus::SphericalTorus(const double r_min, const double r_max, const double min_polar_angle, @@ -55,34 +54,40 @@ SphericalTorus::SphericalTorus(const std::array& radial_range, // * phi : azimuthal angle // template -std::array, 3> SphericalTorus::operator()( - const std::array& source_coords) const { - const auto r = radius(source_coords[0]); - const auto theta = M_PI_2 - (pi_over_2_minus_theta_min_ * source_coords[1]); - const auto phi = M_PI * fraction_of_torus_ * source_coords[2]; - - return { - {r * sin(theta) * cos(phi), r * sin(theta) * sin(phi), r * cos(theta)}}; +tnsr::I SphericalTorus::operator()( + const tnsr::I& source_coords) const { + tnsr::I result; + const auto r = radius(get<0>(source_coords)); + const auto theta = + M_PI_2 - (pi_over_2_minus_theta_min_ * get<1>(source_coords)); + const auto phi = M_PI * fraction_of_torus_ * get<2>(source_coords); + get<0>(result) = r * sin(theta) * cos(phi); + get<1>(result) = r * sin(theta) * sin(phi); + get<2>(result) = r * cos(theta); + return result; } -std::optional> SphericalTorus::inverse( - const std::array& target_coords) const { - const double r = - std::hypot(target_coords[0], target_coords[1], target_coords[2]); - const double theta = std::atan2( - std::hypot(target_coords[0], target_coords[1]), target_coords[2]); - const double phi = std::atan2(target_coords[1], target_coords[0]); - - return {{radius_inverse(r), (M_PI_2 - theta) / pi_over_2_minus_theta_min_, - phi / (M_PI * fraction_of_torus_)}}; +tnsr::I SphericalTorus::inverse( + const tnsr::I& target_coords) const { + tnsr::I result; + const double r = std::hypot(get<0>(target_coords), get<1>(target_coords), + get<2>(target_coords)); + const double theta = + std::atan2(std::hypot(get<0>(target_coords), get<1>(target_coords)), + get<2>(target_coords)); + const double phi = std::atan2(get<1>(target_coords), get<0>(target_coords)); + get<0>(result) = radius_inverse(r); + get<1>(result) = (M_PI_2 - theta) / pi_over_2_minus_theta_min_; + get<2>(result) = phi / (M_PI * fraction_of_torus_); + return result; } template -tnsr::Ij, 3, Frame::NoFrame> -SphericalTorus::jacobian(const std::array& source_coords) const { - using UnwrappedT = tt::remove_cvref_wrap_t; - auto jacobian = make_with_value>( - dereference_wrapper(source_coords[0]), 0.0); +Jacobian SphericalTorus::jacobian( + const tnsr::I& source_coords) const { + auto jacobian = + make_with_value>( + get<0>(source_coords), 0.0); // In order to reduce number of memory allocations we use some slots of // jacobian for storing temp variables as below. @@ -90,18 +95,18 @@ SphericalTorus::jacobian(const std::array& source_coords) const { auto& sin_theta = get<2, 1>(jacobian); auto& cos_theta = get<2, 0>(jacobian); get<2, 2>(jacobian) = - M_PI_2 - (pi_over_2_minus_theta_min_ * source_coords[1]); + M_PI_2 - (pi_over_2_minus_theta_min_ * get<1>(source_coords)); sin_theta = sin(get<2, 2>(jacobian)); cos_theta = cos(get<2, 2>(jacobian)); auto& sin_phi = get<1, 1>(jacobian); auto& cos_phi = get<1, 2>(jacobian); - get<2, 2>(jacobian) = M_PI * fraction_of_torus_ * source_coords[2]; + get<2, 2>(jacobian) = M_PI * fraction_of_torus_ * get<2>(source_coords); sin_phi = sin(get<2, 2>(jacobian)); cos_phi = cos(get<2, 2>(jacobian)); auto& r = get<2, 2>(jacobian); - radius(make_not_null(&r), source_coords[0]); + radius(make_not_null(&r), get<0>(source_coords)); // Note : execution order matters here since we are overwriting each temp // variables with jacobian values corresponding to the slot. @@ -119,11 +124,11 @@ SphericalTorus::jacobian(const std::array& source_coords) const { } template -tnsr::Ij, 3, Frame::NoFrame> -SphericalTorus::inv_jacobian(const std::array& source_coords) const { - using UnwrappedT = tt::remove_cvref_wrap_t; - auto inv_jacobian = make_with_value>( - dereference_wrapper(source_coords[0]), 0.0); +InverseJacobian +SphericalTorus::inv_jacobian(const tnsr::I& source_coords) const { + auto inv_jacobian = make_with_value< + InverseJacobian>( + get<0>(source_coords), 0.0); // In order to reduce number of memory allocations we use some slots of // jacobian for storing temp variables as below. @@ -131,18 +136,18 @@ SphericalTorus::inv_jacobian(const std::array& source_coords) const { auto& sin_theta = get<1, 2>(inv_jacobian); auto& cos_theta = get<0, 2>(inv_jacobian); get<2, 2>(inv_jacobian) = - M_PI_2 - (pi_over_2_minus_theta_min_ * source_coords[1]); + M_PI_2 - (pi_over_2_minus_theta_min_ * get<1>(source_coords)); cos_theta = cos(get<2, 2>(inv_jacobian)); sin_theta = sin(get<2, 2>(inv_jacobian)); auto& sin_phi = get<1, 1>(inv_jacobian); auto& cos_phi = get<2, 1>(inv_jacobian); - get<2, 2>(inv_jacobian) = M_PI * fraction_of_torus_ * source_coords[2]; + get<2, 2>(inv_jacobian) = M_PI * fraction_of_torus_ * get<2>(source_coords); cos_phi = cos(get<2, 2>(inv_jacobian)); sin_phi = sin(get<2, 2>(inv_jacobian)); auto& r = get<2, 2>(inv_jacobian); - radius(make_not_null(&r), source_coords[0]); + radius(make_not_null(&r), get<0>(source_coords)); // Note : execution order matters here since we are overwriting each temp // variables with jacobian values corresponding to the slot. @@ -164,11 +169,10 @@ SphericalTorus::inv_jacobian(const std::array& source_coords) const { } template -tnsr::Ijj, 3, Frame::NoFrame> -SphericalTorus::hessian(const std::array& source_coords) const { - using UnwrappedT = tt::remove_cvref_wrap_t; - auto hessian = make_with_value>( - dereference_wrapper(source_coords[0]), 0.0); +tnsr::Ijj SphericalTorus::hessian( + const tnsr::I& source_coords) const { + auto hessian = make_with_value>( + get<0>(source_coords), 0.0); // In order to reduce number of memory allocations we use some slots of // hessian for storing temp variables as below. @@ -182,14 +186,14 @@ SphericalTorus::hessian(const std::array& source_coords) const { // these two slots are zeros (not used) auto& cos_theta = get<2, 1, 2>(hessian); auto& sin_theta = get<2, 2, 2>(hessian); - get<0, 0, 0>(hessian) = M_PI_2 - (theta_factor * source_coords[1]); + get<0, 0, 0>(hessian) = M_PI_2 - (theta_factor * get<1>(source_coords)); cos_theta = cos(get<0, 0, 0>(hessian)); sin_theta = sin(get<0, 0, 0>(hessian)); // these two slots are NOT zeros, but will be overwritten with hessian later auto& cos_phi = get<2, 0, 1>(hessian); auto& sin_phi = get<2, 1, 1>(hessian); - get<0, 0, 0>(hessian) = phi_factor * source_coords[2]; + get<0, 0, 0>(hessian) = phi_factor * get<2>(source_coords); cos_phi = cos(get<0, 0, 0>(hessian)); sin_phi = sin(get<0, 0, 0>(hessian)); @@ -197,7 +201,7 @@ SphericalTorus::hessian(const std::array& source_coords) const { auto& r_factor = get<0, 0, 0>(hessian); auto& r = get<1, 0, 0>(hessian); r_factor = 0.5 * (r_max_ - r_min_); - radius(make_not_null(&r), source_coords[0]); + radius(make_not_null(&r), get<0>(source_coords)); // Note : execution order matters here get<0, 0, 1>(hessian) = -r_factor * theta_factor * cos_theta * cos_phi; @@ -225,26 +229,24 @@ SphericalTorus::hessian(const std::array& source_coords) const { } template -tnsr::Ijk, 3, Frame::NoFrame> -SphericalTorus::derivative_of_inv_jacobian( - const std::array& source_coords) const { - using UnwrappedT = tt::remove_cvref_wrap_t; - auto result = make_with_value>( - dereference_wrapper(source_coords[0]), 0.0); +tnsr::Ijk SphericalTorus::derivative_of_inv_jacobian( + const tnsr::I& source_coords) const { + auto result = make_with_value>( + get<0>(source_coords), 0.0); // In order to reduce number of memory allocations we use some slots of // `result` for storing temp variables as below. auto& cos_phi = get<1, 2, 2>(result); auto& sin_phi = get<2, 2, 0>(result); - get<0, 0, 0>(result) = M_PI * fraction_of_torus_ * source_coords[2]; + get<0, 0, 0>(result) = M_PI * fraction_of_torus_ * get<2>(source_coords); cos_phi = cos(get<0, 0, 0>(result)); sin_phi = sin(get<0, 0, 0>(result)); auto& cos_theta = get<2, 2, 1>(result); auto& sin_theta = get<2, 2, 2>(result); get<0, 0, 0>(result) = - M_PI_2 - (pi_over_2_minus_theta_min_ * source_coords[1]); + M_PI_2 - (pi_over_2_minus_theta_min_ * get<1>(source_coords)); cos_theta = cos(get<0, 0, 0>(result)); sin_theta = sin(get<0, 0, 0>(result)); @@ -255,7 +257,7 @@ SphericalTorus::derivative_of_inv_jacobian( auto& r = get<0, 0, 0>(result); auto& r_factor = get<0, 1, 0>(result); - radius(make_not_null(&r), source_coords[0]); + radius(make_not_null(&r), get<0>(source_coords)); r_factor = 0.5 * (r_max_ - r_min_); get<0, 0, 1>(result) = -theta_factor / r_factor * cos_theta * cos_phi; @@ -313,18 +315,17 @@ void SphericalTorus::pup(PUP::er& p) { } template -tt::remove_cvref_wrap_t SphericalTorus::radius(const T& x) const { +T SphericalTorus::radius(const T& x) const { return 0.5 * r_min_ * (1.0 - x) + 0.5 * r_max_ * (1.0 + x); } template -void SphericalTorus::radius(const gsl::not_null*> r, - const T& x) const { +void SphericalTorus::radius(const gsl::not_null r, const T& x) const { *r = 0.5 * r_min_ * (1.0 - x) + 0.5 * r_max_ * (1.0 + x); } template -tt::remove_cvref_wrap_t SphericalTorus::radius_inverse(const T& r) const { +T SphericalTorus::radius_inverse(const T& r) const { return ((r - r_min_) - (r_max_ - r)) / (r_max_ - r_min_); } @@ -341,26 +342,24 @@ bool operator!=(const SphericalTorus& lhs, const SphericalTorus& rhs) { #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) #define INSTANTIATE(_, data) \ - template std::array, 3> \ - SphericalTorus::operator()(const std::array& source_coords) \ + template tnsr::I SphericalTorus::operator()( \ + const tnsr::I& source_coords) const; \ + template Jacobian \ + SphericalTorus::jacobian(const tnsr::I& source_coords) \ const; \ - template tnsr::Ij, 3, Frame::NoFrame> \ - SphericalTorus::jacobian(const std::array& source_coords) \ + template InverseJacobian \ + SphericalTorus::inv_jacobian(const tnsr::I& source_coords) \ const; \ - template tnsr::Ij, 3, Frame::NoFrame> \ - SphericalTorus::inv_jacobian( \ - const std::array& source_coords) const; \ - template tnsr::Ijj, 3, Frame::NoFrame> \ - SphericalTorus::hessian(const std::array& source_coords) \ - const; \ - template tnsr::Ijk, 3, Frame::NoFrame> \ + template tnsr::Ijj SphericalTorus::hessian( \ + const tnsr::I& source_coords) const; \ + template tnsr::Ijk \ SphericalTorus::derivative_of_inv_jacobian( \ - const std::array& source_coords) const; + const tnsr::I& source_coords) const; + +GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector)) -GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector, - std::reference_wrapper, - std::reference_wrapper)) #undef INSTANTIATE #undef DTYPE -} // namespace domain::CoordinateMaps +} // namespace grmhd::AnalyticData diff --git a/src/Domain/CoordinateMaps/SphericalTorus.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp similarity index 76% rename from src/Domain/CoordinateMaps/SphericalTorus.hpp rename to src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp index 7ed6933197707..d96fe811dca49 100644 --- a/src/Domain/CoordinateMaps/SphericalTorus.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp @@ -10,9 +10,9 @@ #include #include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" #include "Options/Context.hpp" #include "Options/String.hpp" -#include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp" /// \cond namespace PUP { @@ -24,10 +24,8 @@ class not_null; } // namespace gsl /// \endcond -namespace domain::CoordinateMaps { +namespace grmhd::AnalyticData { /*! - * \ingroup CoordinateMapsGroup - * * \brief Torus made by removing two polar cones from a spherical shell * * Maps source coordinates \f$(\xi, \eta, \zeta)\f$ to @@ -100,27 +98,25 @@ class SphericalTorus { SphericalTorus() = default; template - std::array, 3> operator()( - const std::array& source_coords) const; + tnsr::I operator()(const tnsr::I& source_coords) const; - std::optional> inverse( - const std::array& target_coords) const; + tnsr::I inverse(const tnsr::I& target_coords) const; template - tnsr::Ij, 3, Frame::NoFrame> jacobian( - const std::array& source_coords) const; + Jacobian jacobian( + const tnsr::I& source_coords) const; template - tnsr::Ij, 3, Frame::NoFrame> inv_jacobian( - const std::array& source_coords) const; + InverseJacobian inv_jacobian( + const tnsr::I& source_coords) const; template - tnsr::Ijj, 3, Frame::NoFrame> hessian( - const std::array& source_coords) const; + tnsr::Ijj hessian( + const tnsr::I& source_coords) const; template - tnsr::Ijk, 3, Frame::NoFrame> - derivative_of_inv_jacobian(const std::array& source_coords) const; + tnsr::Ijk derivative_of_inv_jacobian( + const tnsr::I& source_coords) const; // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p); @@ -129,14 +125,13 @@ class SphericalTorus { private: template - tt::remove_cvref_wrap_t radius(const T& x) const; + T radius(const T& x) const; template - void radius(const gsl::not_null*> r, - const T& x) const; + void radius(gsl::not_null r, const T& x) const; template - tt::remove_cvref_wrap_t radius_inverse(const T& r) const; + T radius_inverse(const T& r) const; friend bool operator==(const SphericalTorus& lhs, const SphericalTorus& rhs); @@ -148,4 +143,4 @@ class SphericalTorus { }; bool operator!=(const SphericalTorus& lhs, const SphericalTorus& rhs); -} // namespace domain::CoordinateMaps +} // namespace grmhd::AnalyticData diff --git a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt index 1b0d93306c52e..39c8e39eb1b69 100644 --- a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt +++ b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt @@ -23,7 +23,6 @@ set(LIBRARY_SOURCES Test_ProductMaps.cpp Test_Rotation.cpp Test_SpecialMobius.cpp - Test_SphericalTorus.cpp Test_Tags.cpp Test_TimeDependentHelpers.cpp Test_UniformCylindricalEndcap.cpp diff --git a/tests/Unit/Domain/CoordinateMaps/Test_SphericalTorus.cpp b/tests/Unit/Domain/CoordinateMaps/Test_SphericalTorus.cpp deleted file mode 100644 index abe422f8944da..0000000000000 --- a/tests/Unit/Domain/CoordinateMaps/Test_SphericalTorus.cpp +++ /dev/null @@ -1,125 +0,0 @@ -// Distributed under the MIT License. -// See LICENSE.txt for details. - -#include "Framework/TestingFramework.hpp" - -#include -#include - -#include "Domain/CoordinateMaps/SphericalTorus.hpp" -#include "Framework/TestHelpers.hpp" -#include "Helpers/DataStructures/MakeWithRandomValues.hpp" -#include "Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp" -#include "Utilities/StdArrayHelpers.hpp" - -SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.SphericalTorus", - "[Domain][Unit]") { - // check parse errors first - CHECK_THROWS_WITH( - ([]() { domain::CoordinateMaps::SphericalTorus(-1.0, 2.0, 0.1, 0.5); })(), - Catch::Matchers::ContainsSubstring("Minimum radius must be positive.")); - CHECK_THROWS_WITH( - ([]() { domain::CoordinateMaps::SphericalTorus(3.0, 2.0, 0.1, 0.5); })(), - Catch::Matchers::ContainsSubstring( - "Maximum radius must be greater than minimum radius.")); - CHECK_THROWS_WITH( - ([]() { domain::CoordinateMaps::SphericalTorus(1.0, 2.0, 10.0, 0.5); })(), - Catch::Matchers::ContainsSubstring( - "Minimum polar angle should be less than pi/2.")); - CHECK_THROWS_WITH( - ([]() { domain::CoordinateMaps::SphericalTorus(1.0, 2.0, -0.1, 0.5); })(), - Catch::Matchers::ContainsSubstring( - "Minimum polar angle should be positive")); - CHECK_THROWS_WITH( - ([]() { domain::CoordinateMaps::SphericalTorus(1.0, 2.0, 0.1, -0.5); })(), - Catch::Matchers::ContainsSubstring( - "Fraction of torus included must be positive.")); - CHECK_THROWS_WITH( - ([]() { domain::CoordinateMaps::SphericalTorus(1.0, 2.0, 0.1, 2.0); })(), - Catch::Matchers::ContainsSubstring( - "Fraction of torus included must be at most 1.")); - - // check constructor - CHECK(domain::CoordinateMaps::SphericalTorus(std::array{1.0, 2.0}, - 0.1, 0.5) == - domain::CoordinateMaps::SphericalTorus(1.0, 2.0, 0.1, 0.5)); - - MAKE_GENERATOR(gen); - const double r_min = std::uniform_real_distribution<>(1.0, 2.0)(gen); - const double r_max = std::uniform_real_distribution<>(3.0, 4.0)(gen); - const double phi_max = std::uniform_real_distribution<>(0.5, 1.0)(gen); - const double fraction_of_torus = - std::uniform_real_distribution<>(0.1, 1.0)(gen); - - const domain::CoordinateMaps::SphericalTorus full_torus(r_min, r_max, - phi_max); - const domain::CoordinateMaps::SphericalTorus partial_torus( - r_min, r_max, phi_max, fraction_of_torus); - - // Can't do for full_torus because it is not invertable on the boundary. - test_suite_for_map_on_unit_cube(partial_torus); - - { - const double x = std::uniform_real_distribution<>(-1.0, 1.0)(gen); - const double y = std::uniform_real_distribution<>(-1.0, 1.0)(gen); - CHECK_ITERABLE_APPROX(full_torus(std::array{x, y, -1.0}), - full_torus(std::array{x, y, 1.0})); - } - { - const double x = std::uniform_real_distribution<>(-1.0, 1.0)(gen); - const double y1 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); - const double y2 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); - const double z1 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); - const double z2 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); - CHECK(magnitude(partial_torus(std::array{x, y1, z1})) == - approx(magnitude(partial_torus(std::array{x, y2, z2})))); - } - - test_coordinate_map_implementation(full_torus); - test_coordinate_map_implementation(partial_torus); - - CHECK(full_torus == full_torus); - CHECK_FALSE(full_torus != full_torus); - CHECK(partial_torus == partial_torus); - CHECK_FALSE(partial_torus != partial_torus); - CHECK_FALSE(full_torus == partial_torus); - CHECK(full_torus != partial_torus); - - CHECK(full_torus != - domain::CoordinateMaps::SphericalTorus(r_min + 0.1, r_max, phi_max)); - CHECK(full_torus != - domain::CoordinateMaps::SphericalTorus(r_min, r_max + 0.1, phi_max)); - CHECK(full_torus != - domain::CoordinateMaps::SphericalTorus(r_min, r_max, phi_max + 0.1)); - - CHECK(not full_torus.is_identity()); - CHECK(not partial_torus.is_identity()); - - { - std::uniform_real_distribution<> dist(-1.0, 1.0); - const auto test_point = make_with_random_values>( - make_not_null(&gen), make_not_null(&dist), double{}); - - const auto analytic = partial_torus.hessian(test_point); - const auto analytic_inverse = - partial_torus.derivative_of_inv_jacobian(test_point); - for (size_t i = 0; i < 3; ++i) { - CAPTURE(i); - for (size_t j = 0; j < 3; ++j) { - CAPTURE(j); - for (size_t k = 0; k < 3; ++k) { - CAPTURE(k); - const auto numerical = numerical_derivative( - [&i, &j, &partial_torus](const std::array& x) { - return std::array{partial_torus.jacobian(x).get(i, j), - partial_torus.inv_jacobian(x).get(i, j)}; - }, - test_point, k, 1e-3); - auto deriv_approx = Approx::custom().epsilon(1.0e-9).scale(1.0); - CHECK(analytic.get(i, j, k) == deriv_approx(numerical[0])); - CHECK(analytic_inverse.get(i, j, k) == deriv_approx(numerical[1])); - } - } - } - } -} diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt index 6dfbffaf05524..8d32d5752649f 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt @@ -18,6 +18,7 @@ set(LIBRARY_SOURCES Test_PolarMagnetizedFmDisk.cpp Test_RiemannProblem.cpp Test_SlabJet.cpp + Test_SphericalTorus.cpp ) if (TARGET FUKA::Exporter) diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp index ebfa7ff259175..e6aed89e28eda 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp @@ -9,13 +9,13 @@ #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/Tensor.hpp" -#include "Domain/CoordinateMaps/SphericalTorus.hpp" #include "Framework/TestCreation.hpp" #include "Framework/TestHelpers.hpp" #include "Helpers/PointwiseFunctions/AnalyticData/TestHelpers.hpp" #include "Helpers/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/VerifyGrSolution.hpp" #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp" #include "PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp" #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" @@ -60,30 +60,29 @@ void test_create_from_options() { const auto& disk = dynamic_cast( *deserialized_option_solution); - CHECK(disk == - grmhd::AnalyticData::PolarMagnetizedFmDisk( - grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, - 0.065, 1.654, 0.42, 85.0, 6), - domain::CoordinateMaps::SphericalTorus(2.5, 3.6, 0.9, 0.7))); + CHECK(disk == grmhd::AnalyticData::PolarMagnetizedFmDisk( + grmhd::AnalyticData::MagnetizedFmDisk( + 1.3, 0.345, 6.123, 14.2, 0.065, 1.654, 0.42, 85.0, 6), + grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7))); } void test_move() { grmhd::AnalyticData::PolarMagnetizedFmDisk disk( - grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, - 0.065, 1.654, 0.42, 85.0, 6), - domain::CoordinateMaps::SphericalTorus(2.5, 3.6, 0.9, 0.7)); + grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, 0.065, + 1.654, 0.42, 85.0, 6), + grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7)); grmhd::AnalyticData::PolarMagnetizedFmDisk disk_copy( - grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, - 0.065, 1.654, 0.42, 85.0, 6), - domain::CoordinateMaps::SphericalTorus(2.5, 3.6, 0.9, 0.7)); - test_move_semantics(std::move(disk), disk_copy); + grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, 0.065, + 1.654, 0.42, 85.0, 6), + grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7)); + test_move_semantics(std::move(disk), disk_copy); } void test_serialize() { grmhd::AnalyticData::PolarMagnetizedFmDisk disk( - grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, - 0.065, 1.654, 0.42, 85.0, 6), - domain::CoordinateMaps::SphericalTorus(2.5, 3.6, 0.9, 0.7)); + grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, 0.065, + 1.654, 0.42, 85.0, 6), + grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7)); test_serialization(disk); } @@ -112,7 +111,7 @@ void test_variables(const DataType& used_for_size) { bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, polytropic_constant, polytropic_exponent, threshold_density, inverse_plasma_beta), - domain::CoordinateMaps::SphericalTorus(3.0, 20.0, 1.0, 0.3)); + grmhd::AnalyticData::SphericalTorus(3.0, 20.0, 1.0, 0.3)); const auto coords = make_with_value>(used_for_size, 0.5); diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_SphericalTorus.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_SphericalTorus.cpp new file mode 100644 index 0000000000000..0c83646c3cbd1 --- /dev/null +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_SphericalTorus.cpp @@ -0,0 +1,171 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include + +#include "DataStructures/Tensor/Tensor.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/StdArrayHelpers.hpp" + +namespace { +std::array t2a(const tnsr::I& tens) { + std::array arr{}; + for (size_t i = 0; i < 3; i++) { + gsl::at(arr,i) = tens.get(i); + } + return arr; +} + +tnsr::I a2t(const std::array& arr) { + tnsr::I tens; + for (size_t i = 0; i < 3; i++) { + tens.get(i) = gsl::at(arr,i); + } + return tens; +} +} // namespace + +SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.SphericalTorus", + "[Unit][PointwiseFunctions]") { + // check parse errors first + CHECK_THROWS_WITH( + ([]() { grmhd::AnalyticData::SphericalTorus(-1.0, 2.0, 0.1, 0.5); })(), + Catch::Matchers::ContainsSubstring("Minimum radius must be positive.")); + CHECK_THROWS_WITH( + ([]() { grmhd::AnalyticData::SphericalTorus(3.0, 2.0, 0.1, 0.5); })(), + Catch::Matchers::ContainsSubstring( + "Maximum radius must be greater than minimum radius.")); + CHECK_THROWS_WITH( + ([]() { grmhd::AnalyticData::SphericalTorus(1.0, 2.0, 10.0, 0.5); })(), + Catch::Matchers::ContainsSubstring( + "Minimum polar angle should be less than pi/2.")); + CHECK_THROWS_WITH( + ([]() { grmhd::AnalyticData::SphericalTorus(1.0, 2.0, -0.1, 0.5); })(), + Catch::Matchers::ContainsSubstring( + "Minimum polar angle should be positive")); + CHECK_THROWS_WITH( + ([]() { grmhd::AnalyticData::SphericalTorus(1.0, 2.0, 0.1, -0.5); })(), + Catch::Matchers::ContainsSubstring( + "Fraction of torus included must be positive.")); + CHECK_THROWS_WITH( + ([]() { grmhd::AnalyticData::SphericalTorus(1.0, 2.0, 0.1, 2.0); })(), + Catch::Matchers::ContainsSubstring( + "Fraction of torus included must be at most 1.")); + + // check constructor + CHECK(grmhd::AnalyticData::SphericalTorus(std::array{1.0, 2.0}, + 0.1, 0.5) == + grmhd::AnalyticData::SphericalTorus(1.0, 2.0, 0.1, 0.5)); + + MAKE_GENERATOR(gen); + const double r_min = std::uniform_real_distribution<>(1.0, 2.0)(gen); + const double r_max = std::uniform_real_distribution<>(3.0, 4.0)(gen); + const double phi_max = std::uniform_real_distribution<>(0.5, 1.0)(gen); + const double fraction_of_torus = + std::uniform_real_distribution<>(0.1, 1.0)(gen); + + const grmhd::AnalyticData::SphericalTorus full_torus(r_min, r_max, phi_max); + const grmhd::AnalyticData::SphericalTorus partial_torus(r_min, r_max, phi_max, + fraction_of_torus); + { + const double x = std::uniform_real_distribution<>(-1.0, 1.0)(gen); + const double y = std::uniform_real_distribution<>(-1.0, 1.0)(gen); + const tnsr::I test_pt1{{{x, y, -1.0}}}; + const tnsr::I test_pt2{{{x, y, 1.0}}}; + CHECK_ITERABLE_APPROX(full_torus(test_pt1), full_torus(test_pt2)); + } + + { + const double x = std::uniform_real_distribution<>(-1.0, 1.0)(gen); + const double y1 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); + const double y2 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); + const double z1 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); + const double z2 = std::uniform_real_distribution<>(-1.0, 1.0)(gen); + + tnsr::I test_coords1{{{x, y1, z1}}}; + tnsr::I test_coords2{{{x, y2, z2}}}; + test_coords1 = partial_torus(test_coords1); + test_coords2 = partial_torus(test_coords2); + + CHECK(magnitude(t2a(test_coords1)) == approx(magnitude(t2a(test_coords2)))); + } + + CHECK(full_torus == full_torus); + CHECK_FALSE(full_torus != full_torus); + CHECK(partial_torus == partial_torus); + CHECK_FALSE(partial_torus != partial_torus); + CHECK_FALSE(full_torus == partial_torus); + CHECK(full_torus != partial_torus); + + CHECK(full_torus != + grmhd::AnalyticData::SphericalTorus(r_min + 0.1, r_max, phi_max)); + CHECK(full_torus != + grmhd::AnalyticData::SphericalTorus(r_min, r_max + 0.1, phi_max)); + CHECK(full_torus != + grmhd::AnalyticData::SphericalTorus(r_min, r_max, phi_max + 0.1)); + + CHECK(not full_torus.is_identity()); + CHECK(not partial_torus.is_identity()); + + { + auto deriv_approx = Approx::custom().epsilon(1.0e-9).scale(1.0); + std::uniform_real_distribution<> dist(-1.0, 1.0); + const auto test_point = make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), double{}); + const auto test_tnsr = a2t(test_point); + const tnsr::I mapped_tnsr = partial_torus(test_tnsr); + const auto mapped_point = t2a(mapped_tnsr); + + const tnsr::I iden_tnsr = partial_torus.inverse(mapped_tnsr); + + const auto jac = partial_torus.jacobian(test_tnsr); + const auto jac_inv = partial_torus.inv_jacobian(test_tnsr); + const auto analytic = partial_torus.hessian(test_tnsr); + const auto analytic_inverse = + partial_torus.derivative_of_inv_jacobian(test_tnsr); + for (size_t i = 0; i < 3; ++i) { + CAPTURE(i); + CHECK(test_tnsr.get(i) == deriv_approx(iden_tnsr.get(i))); + for (size_t j = 0; j < 3; ++j) { + CAPTURE(j); + const auto numerical_jacobian = numerical_derivative( + [&i, &partial_torus](const std::array& x) { + const auto x_tnsr = a2t(x); + return partial_torus(x_tnsr).get(i); + }, + test_point, j, 1e-3); + + const auto numerical_inverse_jacobian = numerical_derivative( + [&i, &partial_torus](const std::array& mapped_x) { + const auto mapped_x_tnsr = a2t(mapped_x); + return partial_torus.inverse(mapped_x_tnsr).get(i); + }, + mapped_point, j, 1e-3); + + CHECK(jac.get(i, j) == deriv_approx(numerical_jacobian)); + CHECK(jac_inv.get(i, j) == deriv_approx(numerical_inverse_jacobian)); + for (size_t k = 0; k < 3; ++k) { + CAPTURE(k); + const auto numerical_hessian = numerical_derivative( + [&i, &j, &partial_torus](const std::array& y) { + const auto y_tnsr = a2t(y); + return std::array{partial_torus.jacobian(y_tnsr).get(i, j), + partial_torus.inv_jacobian(y_tnsr).get(i, j)}; + }, + test_point, k, 1e-3); + CHECK(analytic.get(i, j, k) == deriv_approx(numerical_hessian[0])); + CHECK(analytic_inverse.get(i, j, k) == + deriv_approx(numerical_hessian[1])); + } + } + } + } +}