diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt index 3f670d8a6778..d1cb30c1ac79 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/Evolution/Executables/GrMhd/ValenciaDivClean/CMakeLists.txt b/src/Evolution/Executables/GrMhd/ValenciaDivClean/CMakeLists.txt index 2bd312145ada..73e0141199a4 100644 --- a/src/Evolution/Executables/GrMhd/ValenciaDivClean/CMakeLists.txt +++ b/src/Evolution/Executables/GrMhd/ValenciaDivClean/CMakeLists.txt @@ -140,6 +140,12 @@ add_grmhd_executable( "" ) +add_grmhd_executable( + PolarMagnetizedFmDisk + grmhd::AnalyticData::PolarMagnetizedFmDisk + "" + ) + add_grmhd_executable( RiemannProblem grmhd::AnalyticData::RiemannProblem diff --git a/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp b/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp index 9984ad08d574..598535424a34 100644 --- a/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp +++ b/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp @@ -143,6 +143,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" diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt b/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt index 44135fdf2e2e..53b1bc59d56a 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt @@ -17,8 +17,10 @@ spectre_target_sources( MagnetizedFmDisk.cpp MagnetizedTovStar.cpp OrszagTangVortex.cpp + PolarMagnetizedFmDisk.cpp RiemannProblem.cpp SlabJet.cpp + SphericalTorus.cpp ) spectre_target_headers( @@ -35,8 +37,10 @@ spectre_target_headers( MagnetizedFmDisk.hpp MagnetizedTovStar.hpp OrszagTangVortex.hpp + 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 new file mode 100644 index 000000000000..740085b11c58 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.cpp @@ -0,0 +1,40 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp" + +#include +#include + + +namespace grmhd::AnalyticData { + +PolarMagnetizedFmDisk::PolarMagnetizedFmDisk( + 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 { + return std::make_unique(*this); +} + +PolarMagnetizedFmDisk::PolarMagnetizedFmDisk(CkMigrateMessage* msg) + : fm_disk_(msg) {} + +void PolarMagnetizedFmDisk::pup(PUP::er& p) { + p | fm_disk_; + p | torus_map_; +} + +PUP::able::PUP_ID PolarMagnetizedFmDisk::my_PUP_ID = 0; // NOLINT + +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 diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp new file mode 100644 index 000000000000..71bbb322e5d0 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp @@ -0,0 +1,183 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "DataStructures/Tensor/EagerMath/Determinant.hpp" +#include "DataStructures/Tensor/Tensor.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" +#include "Utilities/Gsl.hpp" +#include "Utilities/Serialization/CharmPupable.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 the Kerr-Schild + * Cartesian coordinates, but in a toroidal mesh defined from a torus map + * (see grmhd::AnalyticData::SphericalTorus). + */ +class PolarMagnetizedFmDisk + : public virtual evolution::initial_data::InitialData, + public MarkAsAnalyticData, + public AnalyticDataBase { + public: + struct DiskParameters { + using type = MagnetizedFmDisk; + static constexpr Options::String help = "Parameters for the disk."; + }; + + struct TorusParameters { + using type = grmhd::AnalyticData::SphericalTorus; + static constexpr Options::String help = + "Parameters for the evolution region."; + }; + + using options = tmpl::list; + + static constexpr Options::String help = + "Magnetized Fishbone-Moncrief disk in polar coordinates."; + + PolarMagnetizedFmDisk() = default; + PolarMagnetizedFmDisk(const PolarMagnetizedFmDisk& /*rhs*/) = default; + PolarMagnetizedFmDisk& operator=(const PolarMagnetizedFmDisk& /*rhs*/) = + default; + PolarMagnetizedFmDisk(PolarMagnetizedFmDisk&& /*rhs*/) = default; + PolarMagnetizedFmDisk& operator=(PolarMagnetizedFmDisk&& /*rhs*/) = default; + ~PolarMagnetizedFmDisk() override = default; + + PolarMagnetizedFmDisk(MagnetizedFmDisk fm_disk, + grmhd::AnalyticData::SphericalTorus torus_map); + + auto get_clone() const + -> std::unique_ptr override; + + /// \cond + explicit PolarMagnetizedFmDisk(CkMigrateMessage* msg); + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(PolarMagnetizedFmDisk); + /// \endcond + + /// The grmhd variables. + /// + /// \note The functions are optimized for retrieving the hydro variables + /// before the metric variables. + template + tuples::TaggedTuple variables(const tnsr::I& x, + tmpl::list /*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. + + const tnsr::I observation_coordinates(torus_map_(x)); + + using dependencies = tmpl::map< + tmpl::pair::DerivShift, + gr::Tags::Shift>, + tmpl::pair::DerivSpatialMetric, + gr::Tags::SpatialMetric>>; + using required_tags = tmpl::remove_duplicates< + tmpl::remove...>, + tmpl::no_such_type_>>; + + auto observation_data = + fm_disk_.variables(observation_coordinates, required_tags{}); + + 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]( + const auto& data, auto tag) { + using Tag = decltype(tag); + auto result = + transform::to_different_frame(get(data), jacobian, inv_jacobian); + + if constexpr (std::is_same_v< + Tag, gr::AnalyticSolution<3>::DerivShift>) { + const auto deriv_inv_jacobian = + torus_map_.derivative_of_inv_jacobian(x); + const auto& shift = + get>(data); + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j < 3; ++j) { + for (size_t k = 0; k < 3; ++k) { + result.get(i, j) += + deriv_inv_jacobian.get(j, k, i) * shift.get(k); + } + } + } + } else if constexpr (std::is_same_v< + Tag, gr::AnalyticSolution<3>::DerivSpatialMetric< + DataType>>) { + const auto hessian = torus_map_.hessian(x); + const auto& spatial_metric = + get>(data); + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j < 3; ++j) { + for (size_t k = j; k < 3; ++k) { + for (size_t l = 0; l < 3; ++l) { + for (size_t m = 0; m < 3; ++m) { + result.get(i, j, k) += + (hessian.get(l, j, i) * jacobian.get(m, k) + + hessian.get(l, k, i) * jacobian.get(m, j)) * + spatial_metric.get(l, m); + } + } + } + } + } + } else if constexpr (std::is_same_v< + Tag, gr::Tags::SqrtDetSpatialMetric>) { + 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) override; + + private: + friend bool operator==(const PolarMagnetizedFmDisk& lhs, + const PolarMagnetizedFmDisk& rhs); + + MagnetizedFmDisk fm_disk_; + grmhd::AnalyticData::SphericalTorus torus_map_; +}; + +bool operator!=(const PolarMagnetizedFmDisk& lhs, + const PolarMagnetizedFmDisk& rhs); +} // namespace grmhd::AnalyticData 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 8b1541cddb81..c9c999b25013 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 7ed693319770..d96fe811dca4 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/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp index 52d9cce2c294..7d0cdc555263 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp +++ b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp @@ -154,7 +154,7 @@ FishboneMoncriefDisk::variables( const auto specific_enthalpy = get>( variables(x, tmpl::list>{}, vars, index)); - auto rest_mass_density = make_with_value>(x, 0.0); + auto rest_mass_density = make_with_value>(x, 1.0e-15); variables_impl(vars, [&rest_mass_density, &specific_enthalpy, this]( const size_t s, const double /*potential_at_s*/) { get_element(get(rest_mass_density), s) = diff --git a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt index 1b0d93306c52..39c8e39eb1b6 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 abe422f8944d..000000000000 --- 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/Helpers/PointwiseFunctions/AnalyticData/TestHelpers.hpp b/tests/Unit/Helpers/PointwiseFunctions/AnalyticData/TestHelpers.hpp index 24b5e24d4c07..942932de6828 100644 --- a/tests/Unit/Helpers/PointwiseFunctions/AnalyticData/TestHelpers.hpp +++ b/tests/Unit/Helpers/PointwiseFunctions/AnalyticData/TestHelpers.hpp @@ -23,9 +23,10 @@ void test_tag_retrieval(const Solution& solution, const Coords& coords, tmpl::for_each([&](const auto tag_v) { using tag = tmpl::type_from; const auto single_var = solution.variables(coords, tmpl::list{}); - CHECK(tuples::get(single_var) == tuples::get(vars_from_all_tags)); - CHECK(tuples::get(single_var) == - tuples::get(vars_from_all_tags_reversed)); + CHECK_ITERABLE_APPROX(tuples::get(single_var), + tuples::get(vars_from_all_tags)); + CHECK_ITERABLE_APPROX(tuples::get(single_var), + tuples::get(vars_from_all_tags_reversed)); }); } } // namespace AnalyticData diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt index 1fc7f5a7931c..8d32d5752649 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/CMakeLists.txt @@ -15,8 +15,10 @@ set(LIBRARY_SOURCES Test_MagnetizedFmDisk.cpp Test_MagnetizedTovStar.cpp Test_OrszagTangVortex.cpp + 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 new file mode 100644 index 000000000000..e6aed89e28ed --- /dev/null +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp @@ -0,0 +1,136 @@ +// 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 "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" +#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp" +#include "Utilities/MakeWithValue.hpp" +#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" + +namespace { + +static_assert( + not is_analytic_solution_v, + "PolarMagnetizedFmDisk should be analytic_data, and not an " + "analytic_solution"); +static_assert(is_analytic_data_v, + "PolarMagnetizedFmDisk should be analytic_data, and not an " + "analytic_solution"); + +void test_create_from_options() { + register_classes_with_charm(); + const std::unique_ptr option_solution = + TestHelpers::test_option_tag_factory_creation< + evolution::initial_data::OptionTags::InitialData, + grmhd::AnalyticData::PolarMagnetizedFmDisk>( + "PolarMagnetizedFmDisk:\n" + " DiskParameters:\n" + " BhMass: 1.3\n" + " BhDimlessSpin: 0.345\n" + " InnerEdgeRadius: 6.123\n" + " MaxPressureRadius: 14.2\n" + " PolytropicConstant: 0.065\n" + " PolytropicExponent: 1.654\n" + " ThresholdDensity: 0.42\n" + " InversePlasmaBeta: 85.0\n" + " BFieldNormGridRes: 6\n" + " TorusParameters:\n" + " RadialRange: [2.5, 3.6]\n" + " MinPolarAngle: 0.9\n" + " FractionOfTorus: 0.7\n") + ->get_clone(); + const auto deserialized_option_solution = + serialize_and_deserialize(option_solution); + 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), + 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), + 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), + 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), + grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7)); + test_serialization(disk); +} + +struct Wrapper { + template + auto variables(const tnsr::I& x, const double /*t*/, + Tags tags) const { + return disk_->variables(x, tags); + } + const grmhd::AnalyticData::PolarMagnetizedFmDisk* disk_; +}; + +template +void test_variables(const DataType& used_for_size) { + const double bh_mass = 1.12; + const double bh_dimless_spin = 0.97; + const double inner_edge_radius = 6.2; + const double max_pressure_radius = 11.6; + const double polytropic_constant = 0.034; + const double polytropic_exponent = 1.65; + const double threshold_density = 0.14; + const double inverse_plasma_beta = 0.023; + + const grmhd::AnalyticData::PolarMagnetizedFmDisk disk( + grmhd::AnalyticData::MagnetizedFmDisk( + bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, + polytropic_constant, polytropic_exponent, threshold_density, + inverse_plasma_beta), + grmhd::AnalyticData::SphericalTorus(3.0, 20.0, 1.0, 0.3)); + + const auto coords = make_with_value>(used_for_size, 0.5); + + TestHelpers::AnalyticData::test_tag_retrieval( + disk, coords, + typename grmhd::AnalyticData::PolarMagnetizedFmDisk::tags{}); + + if constexpr (std::is_same_v) { + TestHelpers::VerifyGrSolution::verify_consistency(Wrapper{&disk}, 0.0, + coords, 1.0e-2, 1.0e-10); + } +} +} // namespace + +SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.PolarMagFmDisk", + "[Unit][PointwiseFunctions]") { + test_create_from_options(); + test_serialize(); + test_move(); + test_variables(std::numeric_limits::signaling_NaN()); + test_variables(DataVector(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 000000000000..0c83646c3cbd --- /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])); + } + } + } + } +}