diff --git a/src/DataStructures/CachedTempBuffer.hpp b/src/DataStructures/CachedTempBuffer.hpp index 5574344e40a3a..dcc437bca4b48 100644 --- a/src/DataStructures/CachedTempBuffer.hpp +++ b/src/DataStructures/CachedTempBuffer.hpp @@ -50,6 +50,7 @@ class CachedTempBuffer { } return get(data_); } + size_t number_of_grid_points() const { return data_.number_of_grid_points(); } private: template diff --git a/src/DataStructures/TempBuffer.hpp b/src/DataStructures/TempBuffer.hpp index 72d21a7df757b..7b99121809387 100644 --- a/src/DataStructures/TempBuffer.hpp +++ b/src/DataStructures/TempBuffer.hpp @@ -33,6 +33,8 @@ template struct TempBuffer : tuples::tagged_tuple_from_typelist { explicit TempBuffer(const size_t /*size*/) : tuples::tagged_tuple_from_typelist::TaggedTuple() {} + + static size_t number_of_grid_points() { return 1; } }; template diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp index 2b8bb09b6b461..72ef6fa2bb4a3 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp @@ -97,7 +97,7 @@ MagnetizedFmDisk::MagnetizedFmDisk( get>(unmagnetized_vars); const tnsr::I x_max{ - {{fm_disk_.max_pressure_radius_, fm_disk_.bh_spin_a_, 0.0}}}; + {{fm_disk_.max_pressure_radius_, 0.0, 0.0}}}; const double b_squared_max = max( get(dot_product(b_field, b_field, spatial_metric)) / @@ -137,11 +137,13 @@ tnsr::I MagnetizedFmDisk::unnormalized_magnetic_field( auto magnetic_field = make_with_value>(x, 0.0); + auto x_ks = x; + // The maximum pressure (and hence the maximum rest mass density) is located - // on the ring x^2 + y^2 = r_max^2 + a^2, z = 0. + // on the ring x^2 + y^2 = r_max^2, z = 0. // Note that `x` may or may not include points on this ring. const tnsr::I x_max{ - {{fm_disk_.max_pressure_radius_, fm_disk_.bh_spin_a_, 0.0}}}; + {{fm_disk_.max_pressure_radius_, 0.0, 0.0}}}; const double threshold_rest_mass_density = threshold_density_ * get(get>(variables( @@ -174,13 +176,10 @@ tnsr::I MagnetizedFmDisk::unnormalized_magnetic_field( // (r, theta, phi) coordinates. We need to get B^r, B^theta, B^phi first, // then we transform to Cartesian coordinates. const double z_squared = square(get_element(get<2>(x), s)); - const double a_squared = fm_disk_.bh_spin_a_ * fm_disk_.bh_spin_a_; - double sin_theta_squared = square(get_element(get<0>(x), s)) + square(get_element(get<1>(x), s)); - double r_squared = 0.5 * (sin_theta_squared + z_squared - a_squared); - r_squared += sqrt(square(r_squared) + a_squared * z_squared); - sin_theta_squared /= (r_squared + a_squared); + double r_squared = sin_theta_squared + z_squared; + sin_theta_squared /= r_squared; // B^i is proportional to derivatives of the magnetic potential. One has // @@ -216,20 +215,54 @@ tnsr::I MagnetizedFmDisk::unnormalized_magnetic_field( prefactor * (mag_potential(radius - small, sin_theta_squared) - mag_potential(radius + small, sin_theta_squared)); // phi component of the field vanishes identically. + + // Need x in KS coordinates instead of SKS coordinates + // to do transformation below. copy x into x_sks and alter at each + // element index s. + const double sks_to_ks_factor = + sqrt(r_squared + square(fm_disk_.bh_spin_a_)) / radius; + get_element(x_ks.get(0), s) = get_element(x.get(0), s) * sks_to_ks_factor; + get_element(x_ks.get(1), s) = get_element(x.get(1), s) * sks_to_ks_factor; + get_element(x_ks.get(2), s) = get_element(x.get(2), s); } } - return kerr_schild_coords_.cartesian_from_spherical_ks( - std::move(magnetic_field), x); + // magnetic field in KS_coords + // change back to SKS_coords + const auto magnetic_field_ks = + kerr_schild_coords_.cartesian_from_spherical_ks(std::move(magnetic_field), + x_ks); + + using inv_jacobian = + gr::Solutions::SphericalKerrSchild::internal_tags::inv_jacobian< + DataType, Frame::Inertial>; + + FmDisk::IntermediateVariables vars(x); + + const auto inv_jacobians = + get(fm_disk_.background_spacetime_.variables( + x, 0.0, tmpl::list{}, + make_not_null(&vars.sph_kerr_schild_cache))); + + auto magnetic_field_sks = + make_with_value>(x, 0.0); + for (size_t s = 0; s < get_size(get<0>(x)); ++s) { + for (size_t j = 0; j < 3; ++j) { + for (size_t i = 0; i < 3; ++i) { + get_element(magnetic_field_sks.get(j), s) += + get_element(inv_jacobians.get(j, i), s) * + get_element(magnetic_field_ks.get(i), s); + } + } + } + return magnetic_field_sks; } -template +template tuples::TaggedTuple> MagnetizedFmDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const typename FmDisk::IntermediateVariables& /*vars*/, - const size_t /*index*/) const { + gsl::not_null*> /*vars*/) const { auto result = unnormalized_magnetic_field(x); for (size_t i = 0; i < 3; ++i) { result.get(i) *= b_field_normalization_; @@ -252,20 +285,18 @@ bool operator!=(const MagnetizedFmDisk& lhs, const MagnetizedFmDisk& rhs) { } #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define NEED_SPACETIME(data) BOOST_PP_TUPLE_ELEM(1, data) -#define INSTANTIATE(_, data) \ - template tuples::TaggedTuple> \ - MagnetizedFmDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const typename FmDisk::FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; +#define INSTANTIATE(_, data) \ + template tuples::TaggedTuple> \ + MagnetizedFmDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null< \ + FmDisk::FishboneMoncriefDisk::IntermediateVariables*> \ + vars) const; -GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector), (true, false)) +GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector)) #undef DTYPE -#undef NEED_SPACETIME #undef INSTANTIATE } // namespace grmhd::AnalyticData diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp index 819ed2a119d71..62b56e66a0bc1 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp @@ -153,38 +153,19 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData, using equation_of_state_type = typename FmDisk::equation_of_state_type; /// @{ - /// The grmhd variables in Cartesian Kerr-Schild coordinates at `(x, t)` - /// - /// \note The functions are optimized for retrieving the hydro variables - /// before the metric variables. + /// The variables in Cartesian Kerr-Schild coordinates at `(x, t)`. template tuples::TaggedTuple variables(const tnsr::I& x, tmpl::list /*meta*/) const { // Can't store IntermediateVariables as member variable because we // need to be threadsafe. - constexpr double dummy_time = 0.0; - typename FmDisk::IntermediateVariables< - DataType, - tmpl2::flat_any_v<( - std::is_same_v> or - std::is_same_v> or - not tmpl::list_contains_v, Tags>)...>> - vars(fm_disk_.bh_spin_a_, fm_disk_.background_spacetime_, x, dummy_time, - FmDisk::index_helper( - tmpl::index_of, - hydro::Tags::SpatialVelocity>{}), - FmDisk::index_helper( - tmpl::index_of, - hydro::Tags::LorentzFactor>{})); + typename FmDisk::IntermediateVariables vars(x); return {std::move(get([this, &x, &vars]() { if constexpr (std::is_same_v, Tags>) { - return variables(x, tmpl::list{}, vars, - tmpl::index_of, Tags>::value); + return variables(x, tmpl::list{}, make_not_null(&vars)); } else { - return fm_disk_.variables( - x, tmpl::list{}, vars, - tmpl::index_of, Tags>::value); + return fm_disk_.variables(x, tmpl::list{}, make_not_null(&vars)); } }()))...}; } @@ -194,20 +175,12 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData, tmpl::list /*meta*/) const { // Can't store IntermediateVariables as member variable because we need to // be threadsafe. - constexpr double dummy_time = 0.0; - typename FmDisk::IntermediateVariables< - DataType, - std::is_same_v> or - std::is_same_v> or - not tmpl::list_contains_v, Tag>> - intermediate_vars(fm_disk_.bh_spin_a_, fm_disk_.background_spacetime_, - x, dummy_time, std::numeric_limits::max(), - std::numeric_limits::max()); + typename FmDisk::IntermediateVariables vars(x); if constexpr (std::is_same_v, Tag>) { - return variables(x, tmpl::list{}, intermediate_vars, 0); + return variables(x, tmpl::list{}, make_not_null(&vars)); } else { - return fm_disk_.variables(x, tmpl::list{}, intermediate_vars, 0); + return fm_disk_.variables(x, tmpl::list{}, make_not_null(&vars)); } } /// @} @@ -220,14 +193,11 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData, } private: - template - auto variables( - const tnsr::I& x, - tmpl::list> /*meta*/, - const typename FmDisk::IntermediateVariables& - vars, - size_t index) const - -> tuples::TaggedTuple>; + template + auto variables(const tnsr::I& x, + tmpl::list> /*meta*/, + gsl::not_null*> vars) + const -> tuples::TaggedTuple>; template tnsr::I unnormalized_magnetic_field( diff --git a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp index bb055182ffeec..8ea1e142f67b6 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp @@ -19,6 +19,7 @@ #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp" #include "Utilities/ContainerHelpers.hpp" #include "Utilities/ForceInline.hpp" +#include "Utilities/Gsl.hpp" #include "Utilities/MakeArray.hpp" #include "Utilities/TMPL.hpp" #include "Utilities/TaggedTuple.hpp" @@ -429,20 +430,6 @@ class SphericalKerrSchild : public AnalyticSolution<3_st>, SphericalKerrSchild& operator=(SphericalKerrSchild&& /*rhs*/) = default; ~SphericalKerrSchild() = default; - template - tuples::TaggedTuple variables( - const tnsr::I& x, double /*t*/, - tmpl::list /*meta*/) const { - static_assert( - tmpl2::flat_all_v< - tmpl::list_contains_v, Tags>...>, - "At least one of the requested tags is not supported. The requested " - "tags are listed as template parameters of the `variables` function."); - IntermediateVars cache(get_size(*x.begin())); - IntermediateComputer computer(*this, x); - return {cache.get_var(computer, Tags{})...}; - } - // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p); @@ -568,6 +555,47 @@ class SphericalKerrSchild : public AnalyticSolution<3_st>, gr::Tags::SpatialChristoffelFirstKind, gr::Tags::SpatialChristoffelSecondKind>; + // forward-declaration needed. + template + class IntermediateVars; + + template + using allowed_tags = + tmpl::push_back, + typename internal_tags::inv_jacobian>; + + // Overload and call the second one. + template + tuples::TaggedTuple variables( + const tnsr::I& x, double /*t*/, + tmpl::list /*meta*/) const { + static_assert( + tmpl2::flat_all_v< + tmpl::list_contains_v, Tags>...>, + "At least one of the requested tags is not supported. The requested " + "tags are listed as template parameters of the `variables` function."); + IntermediateVars cache(get_size(*x.begin())); + IntermediateComputer computer(*this, x); + return {cache.get_var(computer, Tags{})...}; + } + + template + tuples::TaggedTuple variables( + const tnsr::I& x, double /*t*/, + tmpl::list /*meta*/, + gsl::not_null*> cache) const { + static_assert( + tmpl2::flat_all_v< + tmpl::list_contains_v, Tags>...>, + "At least one of the requested tags is not supported. The requested " + "tags are listed as template parameters of the `variables` function."); + if (cache->number_of_grid_points() != get_size(*x.begin())) { + *cache = IntermediateVars(get_size(*x.begin())); + } + IntermediateComputer computer(*this, x); + return {cache->get_var(computer, Tags{})...}; + } + template class IntermediateComputer { public: diff --git a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp index 7d0cdc5552639..4cc5eb195c107 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp +++ b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp @@ -11,7 +11,7 @@ #include "DataStructures/DataVector.hpp" // IWYU pragma: keep #include "DataStructures/Tensor/EagerMath/DotProduct.hpp" #include "DataStructures/Tensor/Tensor.hpp" -#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" #include "Utilities/ConstantExpressions.hpp" @@ -75,13 +75,13 @@ void FishboneMoncriefDisk::pup(PUP::er& p) { p | equation_of_state_; p | background_spacetime_; } - +// Sigma in Fishbone&Moncrief eqn (3.5) template DataType FishboneMoncriefDisk::sigma(const DataType& r_sqrd, const DataType& sin_theta_sqrd) const { return r_sqrd + square(bh_spin_a_) * (1.0 - sin_theta_sqrd); } - +// Inverse of A in Fishbone&Moncrief eqn (3.5) template DataType FishboneMoncriefDisk::inv_ucase_a(const DataType& r_sqrd, const DataType& sin_theta_sqrd, @@ -123,38 +123,24 @@ DataType FishboneMoncriefDisk::potential(const DataType& r_sqrd, log(sqrt(four_velocity_t_sqrd(r_sqrd, sin_theta_sqrd))); } -template -FishboneMoncriefDisk::IntermediateVariables:: - IntermediateVariables(const double bh_spin_a, - const gr::Solutions::KerrSchild& background_spacetime, - const tnsr::I& x, const double t, - size_t in_spatial_velocity_index, - size_t in_lorentz_factor_index) - : spatial_velocity_index(in_spatial_velocity_index), - lorentz_factor_index(in_lorentz_factor_index) { - const double a_squared = bh_spin_a * bh_spin_a; +template +FishboneMoncriefDisk::IntermediateVariables::IntermediateVariables( + const tnsr::I& x) { sin_theta_squared = square(get<0>(x)) + square(get<1>(x)); - r_squared = 0.5 * (sin_theta_squared + square(get<2>(x)) - a_squared); - r_squared += sqrt(square(r_squared) + a_squared * square(get<2>(x))); - sin_theta_squared /= (r_squared + a_squared); - - if (NeedSpacetime) { - kerr_schild_soln = background_spacetime.variables( - x, t, gr::Solutions::KerrSchild::tags{}); - } + r_squared = sin_theta_squared + square(get<2>(x)); + sin_theta_squared /= r_squared; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - const size_t index) const { - const auto specific_enthalpy = get>( - variables(x, tmpl::list>{}, vars, - index)); - auto rest_mass_density = make_with_value>(x, 1.0e-15); + gsl::not_null*> vars) const { + const auto specific_enthalpy = + get>(variables( + x, tmpl::list>{}, vars)); + auto rest_mass_density = make_with_value>(x, 0.0); 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) = @@ -164,13 +150,12 @@ FishboneMoncriefDisk::variables( return {std::move(rest_mass_density)}; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& /* vars */, - const size_t /* index */) const { + gsl::not_null*> /* vars */) const { // Need to add EoS call to get correct electron fraction // when using tables @@ -179,13 +164,12 @@ FishboneMoncriefDisk::variables( return {std::move(ye)}; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - const size_t /*index*/) const { + gsl::not_null*> vars) const { const double inner_edge_potential = potential(square(inner_edge_radius_), 1.0); auto specific_enthalpy = make_with_value>(x, 1.0); @@ -197,16 +181,14 @@ FishboneMoncriefDisk::variables( return {std::move(specific_enthalpy)}; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - const size_t index) const { + gsl::not_null*> vars) const { const auto rest_mass_density = get>( - variables(x, tmpl::list>{}, vars, - index)); + variables(x, tmpl::list>{}, vars)); auto pressure = make_with_value>(x, 0.0); variables_impl(vars, [&pressure, &rest_mass_density, this]( const size_t s, const double /*potential_at_s*/) { @@ -217,16 +199,14 @@ FishboneMoncriefDisk::variables( return {std::move(pressure)}; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - const size_t index) const { + gsl::not_null*> vars) const { const auto rest_mass_density = get>( - variables(x, tmpl::list>{}, vars, - index)); + variables(x, tmpl::list>{}, vars)); auto specific_internal_energy = make_with_value>(x, 0.0); variables_impl(vars, [&specific_internal_energy, &rest_mass_density, this]( const size_t s, const double /*potential_at_s*/) { @@ -237,16 +217,14 @@ FishboneMoncriefDisk::variables( return {std::move(specific_internal_energy)}; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - const size_t index) const { + gsl::not_null*> vars) const { const auto rest_mass_density = get>( - variables(x, tmpl::list>{}, vars, - index)); + variables(x, tmpl::list>{}, vars)); auto temperature = make_with_value>(x, 0.0); variables_impl(vars, [&temperature, &rest_mass_density, this]( @@ -263,26 +241,57 @@ tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - const size_t /*index*/) const { + gsl::not_null*> vars) const { auto spatial_velocity = make_with_value>(x, 0.0); variables_impl(vars, [&spatial_velocity, &vars, &x, this]( const size_t s, const double /*potential_at_s*/) { - const double ang_velocity = angular_velocity( - get_element(vars.r_squared, s), get_element(vars.sin_theta_squared, s)); - - auto transport_velocity = make_array<3>(0.0); - transport_velocity[0] -= ang_velocity * get_element(x.get(1), s); - transport_velocity[1] += ang_velocity * get_element(x.get(0), s); - + const double ang_velocity = + angular_velocity(get_element(vars->r_squared, s), + get_element(vars->sin_theta_squared, s)); + // We first compute the transport velocity in Kerr-Schild coordinates + // and then transform this vector back to Spherical Kerr-Schild coordinates. + auto transport_velocity_ks = make_array<3>(0.0); + auto transport_velocity_sks = make_array<3>(0.0); + + const double sks_to_ks_factor = + sqrt(square(bh_spin_a_) + get_element(vars->r_squared, s)) / + sqrt(get_element(vars->r_squared, s)); + transport_velocity_ks[0] -= + ang_velocity * get_element(x.get(1), s) * sks_to_ks_factor; + transport_velocity_ks[1] += + ang_velocity * get_element(x.get(0), s) * sks_to_ks_factor; + + using inv_jacobian = + gr::Solutions::SphericalKerrSchild::internal_tags::inv_jacobian< + DataType, Frame::Inertial>; + + const auto inv_jacobians = + get(background_spacetime_.variables( + x, 0.0, tmpl::list{}, + make_not_null(&vars->sph_kerr_schild_cache))); + + for (size_t j = 0; j < 3; ++j) { + for (size_t i = 0; i < 3; ++i) { + gsl::at(transport_velocity_sks, j) += + get_element(inv_jacobians.get(j, i), s) * + gsl::at(transport_velocity_ks, i); + } + } for (size_t i = 0; i < 3; ++i) { get_element(spatial_velocity.get(i), s) = - (gsl::at(transport_velocity, i) + + (gsl::at(transport_velocity_sks, i) + get_element( - get>(vars.kerr_schild_soln).get(i), + get>( + background_spacetime_.variables( + x, 0.0, tmpl::list>{}, + make_not_null(&vars->sph_kerr_schild_cache))) + .get(i), s)) / - get_element( - get(get>(vars.kerr_schild_soln)), s); + get_element(get(get>( + background_spacetime_.variables( + x, 0.0, tmpl::list>{}, + make_not_null(&vars->sph_kerr_schild_cache)))), + s); } }); return {std::move(spatial_velocity)}; @@ -293,44 +302,45 @@ tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - const size_t index) const { - const auto spatial_velocity = get>( - variables(x, tmpl::list>{}, - vars, index)); + gsl::not_null*> vars) const { + const auto spatial_velocity = + get>(variables( + x, tmpl::list>{}, vars)); Scalar lorentz_factor{ 1.0 / - sqrt(1.0 - get(dot_product(spatial_velocity, spatial_velocity, - get>( - vars.kerr_schild_soln))))}; + sqrt(1.0 - get(dot_product( + spatial_velocity, spatial_velocity, + get>( + background_spacetime_.variables( + x, 0.0, + tmpl::list>{}, + make_not_null(&vars->sph_kerr_schild_cache))))))}; return {std::move(lorentz_factor)}; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& /*vars*/, - const size_t /*index*/) const { + gsl::not_null*> /*vars*/) const { return {make_with_value>(x, 0.0)}; } -template +template tuples::TaggedTuple> FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& /*vars*/, - const size_t /*index*/) const { + gsl::not_null*> /*vars*/) const { return {make_with_value>(x, 0.0)}; } -template +template void FishboneMoncriefDisk::variables_impl( - const IntermediateVariables& vars, Func f) const { - const DataType& r_squared = vars.r_squared; - const DataType& sin_theta_squared = vars.sin_theta_squared; + gsl::not_null*> vars, Func f) const { + const DataType& r_squared = vars->r_squared; + const DataType& sin_theta_squared = vars->sin_theta_squared; const double inner_edge_potential = potential(square(inner_edge_radius_), 1.0); @@ -368,95 +378,77 @@ bool operator!=(const FishboneMoncriefDisk& lhs, } #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define NEED_SPACETIME(data) BOOST_PP_TUPLE_ELEM(1, data) - -#define INSTANTIATE(_, data) \ - template class FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>; \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; \ - template tuples::TaggedTuple< \ - hydro::Tags::SpecificInternalEnergy> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; \ - template tuples::TaggedTuple< \ - hydro::Tags::DivergenceCleaningField> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables< \ - DTYPE(data), NEED_SPACETIME(data)>& vars, \ - const size_t) const; - -GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector), (true, false)) -#undef INSTANTIATE -#define INSTANTIATE(_, data) \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables& \ - vars, \ - const size_t) const; \ - template tuples::TaggedTuple> \ - FishboneMoncriefDisk::variables( \ - const tnsr::I& x, \ - tmpl::list> /*meta*/, \ - const FishboneMoncriefDisk::IntermediateVariables& \ - vars, \ - const size_t) const; \ - template DTYPE(data) FishboneMoncriefDisk::potential( \ +#define INSTANTIATE(_, data) \ + template class FishboneMoncriefDisk::IntermediateVariables; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple< \ + hydro::Tags::SpecificInternalEnergy> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple< \ + hydro::Tags::DivergenceCleaningField> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template tuples::TaggedTuple> \ + FishboneMoncriefDisk::variables( \ + const tnsr::I& x, \ + tmpl::list> /*meta*/, \ + gsl::not_null*> \ + vars) const; \ + template DTYPE(data) FishboneMoncriefDisk::potential( \ const DTYPE(data) & r_sqrd, const DTYPE(data) & sin_theta_sqrd) const; GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector)) #undef DTYPE -#undef NEED_SPACETIME #undef INSTANTIATE } // namespace RelativisticEuler::Solutions diff --git a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp index 28bee599c04a4..e178c73e112fc 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp @@ -9,7 +9,7 @@ #include "DataStructures/Tensor/TypeAliases.hpp" #include "Options/String.hpp" #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp" -#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp" #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Solutions.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp" // IWYU pragma: keep #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -51,10 +51,7 @@ namespace Solutions { * where \f$u^\mu\f$ is the 4-velocity. Then, all the fluid quantities are * assumed to share the same symmetries as those of the background spacetime, * namely they are stationary (independent of \f$t\f$), and axially symmetric - * (independent of \f$\phi\f$). Note that \f$r\f$ and \f$\theta\f$ can also be - * interpreted as Kerr (a.k.a. "spherical" Kerr-Schild) coordinates - * (see gr::KerrSchildCoords), and that the symmetries of the equilibrium ensure - * that \f$u^\mu\f$ as defined above is also the 4-velocity in Kerr coordinates. + * (independent of \f$\phi\f$). * * Self-gravity is neglected, so that the fluid * variables are determined as functions of the metric. Following the treatment @@ -161,7 +158,7 @@ class FishboneMoncriefDisk public AnalyticSolution<3>, public hydro::TemperatureInitialization { protected: - template + template struct IntermediateVariables; public: @@ -258,64 +255,32 @@ class FishboneMoncriefDisk DataType potential(const DataType& r_sqrd, const DataType& sin_theta_sqrd) const; - SPECTRE_ALWAYS_INLINE static size_t index_helper( - tmpl::no_such_type_ /*meta*/) { - return std::numeric_limits::max(); - } - template - SPECTRE_ALWAYS_INLINE static size_t index_helper(T /*meta*/) { - return T::value; - } - template - using tags = tmpl::append, - typename gr::Solutions::KerrSchild::tags>; + using tags = + tmpl::append, + typename gr::Solutions::SphericalKerrSchild::tags>; /// @{ - /// The fluid variables in Cartesian Kerr-Schild coordinates at `(x, t)` - /// - /// \note The functions are optimized for retrieving the hydro variables - /// before the metric variables. + /// The variables in Cartesian Spherical-Kerr-Schild coordinates at `(x, t)` template - tuples::TaggedTuple variables( - const tnsr::I& x, - const double t, // NOLINT(readability-avoid-const-params-in-decls) - tmpl::list /*meta*/) const { + tuples::TaggedTuple variables(const tnsr::I& x, + const double /*t*/, + tmpl::list /*meta*/) const { // Can't store IntermediateVariables as member variable because we need to // be threadsafe. - IntermediateVariables< - DataType, - tmpl2::flat_any_v<( - std::is_same_v> or - std::is_same_v> or - not tmpl::list_contains_v, Tags>)...>> - vars(bh_spin_a_, background_spacetime_, x, t, - index_helper( - tmpl::index_of, - hydro::Tags::SpatialVelocity>{}), - index_helper( - tmpl::index_of, - hydro::Tags::LorentzFactor>{})); - return {std::move(get( - variables(x, tmpl::list{}, vars, - tmpl::index_of, Tags>::value)))...}; + IntermediateVariables vars(x); + return {std::move( + get(variables(x, tmpl::list{}, make_not_null(&vars))))...}; } template tuples::TaggedTuple variables(const tnsr::I& x, - const double t, // NOLINT + const double /*t*/, tmpl::list /*meta*/) const { // Can't store IntermediateVariables as member variable because we need to // be threadsafe. - IntermediateVariables< - DataType, - std::is_same_v> or - std::is_same_v> or - not tmpl::list_contains_v, Tag>> - intermediate_vars(bh_spin_a_, background_spacetime_, x, t, - std::numeric_limits::max(), - std::numeric_limits::max()); - return variables(x, tmpl::list{}, intermediate_vars, 0); + IntermediateVariables vars(x); + return variables(x, tmpl::list{}, make_not_null(&vars)); } /// @} @@ -327,124 +292,102 @@ class FishboneMoncriefDisk } protected: - template + template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; - template + template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; - template + template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; - template + template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; - template + template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; - template + template auto variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; - template + template auto variables(const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; - template + template auto variables( const tnsr::I& x, tmpl::list> /*meta*/, - const IntermediateVariables& vars, - size_t index) const + gsl::not_null*> vars) const -> tuples::TaggedTuple>; // Grab the metric variables template , - hydro::Tags::SpecificEnthalpy>, + hydro::Tags::SpecificEnthalpy, + hydro::Tags::SpatialVelocity, + hydro::Tags::LorentzFactor>, Tag>> = nullptr> tuples::TaggedTuple variables( - const tnsr::I& /*x*/, tmpl::list /*meta*/, - IntermediateVariables& vars, const size_t index) const { - // The only hydro variables that use the GR solution are spatial velocity - // and Lorentz factor. If those have already been set (their index in the - // typelist is lower than our index) then we can safely std::move the GR - // solution out, assuming nobody gets the same variable twice (e.g. trying - // to retrieve the lapse twice from the solution in the same call to the - // function), but then TaggedTuple would explode horribly, so nothing to - // worry about. - // Analytic solutions are used in Dirichlet boundary conditions and thus are - // called quite frequently, making some optimization worthwhile. - if (index > vars.spatial_velocity_index and - index > vars.lorentz_factor_index) { - return {std::move(get(vars.kerr_schild_soln))}; - } - return {get(vars.kerr_schild_soln)}; + const tnsr::I& x, tmpl::list /*meta*/, + gsl::not_null*> vars) const { + return {get(background_spacetime_.variables( + x, 0.0, tmpl::list{}, + make_not_null(&vars->sph_kerr_schild_cache)))}; } - template - void variables_impl( - const IntermediateVariables& vars, Func f) const; + template + void variables_impl(gsl::not_null*> vars, + Func f) const; // Intermediate variables needed to set several of the Fishbone-Moncrief // solution's variables. - template + + template struct IntermediateVariables { - IntermediateVariables(double bh_spin_a, - const gr::Solutions::KerrSchild& background_spacetime, - const tnsr::I& x, double t, - size_t in_spatial_velocity_index, - size_t in_lorentz_factor_index); + explicit IntermediateVariables(const tnsr::I& x); DataType r_squared{}; DataType sin_theta_squared{}; - tuples::tagged_tuple_from_typelist< - typename gr::Solutions::KerrSchild::tags> - kerr_schild_soln{}; - size_t spatial_velocity_index; - size_t lorentz_factor_index; + gr::Solutions::SphericalKerrSchild::IntermediateVars + sph_kerr_schild_cache = + gr::Solutions::SphericalKerrSchild::IntermediateVars< + DataType, Frame::Inertial>(0); }; friend bool operator==(const FishboneMoncriefDisk& lhs, @@ -459,7 +402,7 @@ class FishboneMoncriefDisk double polytropic_exponent_ = std::numeric_limits::signaling_NaN(); double angular_momentum_ = std::numeric_limits::signaling_NaN(); EquationsOfState::PolytropicFluid equation_of_state_{}; - gr::Solutions::KerrSchild background_spacetime_{}; + gr::Solutions::SphericalKerrSchild background_spacetime_{}; }; bool operator!=(const FishboneMoncriefDisk& lhs, diff --git a/src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt b/src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt index 00abdb311c26b..35b196a04ee69 100644 --- a/src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt +++ b/src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt @@ -26,6 +26,7 @@ spectre_target_sources( SpacetimeNormalOneForm.cpp SpacetimeNormalVector.cpp SpatialMetric.cpp + SphericalKerrSchildCoords.cpp TimeDerivativeOfSpacetimeMetric.cpp TimeDerivativeOfSpatialMetric.cpp WeylElectric.cpp @@ -57,6 +58,7 @@ spectre_target_headers( SpacetimeNormalOneForm.hpp SpacetimeNormalVector.hpp SpatialMetric.hpp + SphericalKerrSchildCoords.hpp Tags.hpp TagsDeclarations.hpp TimeDerivativeOfSpacetimeMetric.hpp diff --git a/src/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.cpp b/src/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.cpp new file mode 100644 index 0000000000000..f20b83e036a93 --- /dev/null +++ b/src/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.cpp @@ -0,0 +1,122 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.hpp" + +#include +#include +#include +#include + +#include "DataStructures/Tensor/EagerMath/DotProduct.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/ConstantExpressions.hpp" +#include "Utilities/ContainerHelpers.hpp" +#include "Utilities/EqualWithinRoundoff.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/MakeWithValue.hpp" + +namespace gr { + +SphericalKerrSchildCoords::SphericalKerrSchildCoords( + const double bh_mass, const double bh_dimless_spin) + : spin_a_(bh_mass * bh_dimless_spin) { + ASSERT(bh_mass > 0, + "The mass must be positive. The value given was " << bh_mass << "."); + ASSERT(bh_dimless_spin > -1.0 and bh_dimless_spin < 1.0, + "The dimensionless spin must be in the range (-1, 1). The value given " + "was " + << bh_dimless_spin << "."); +} + +void SphericalKerrSchildCoords::pup(PUP::er& p) { p | spin_a_; } + +template +tnsr::Ij +SphericalKerrSchildCoords::jacobian_matrix(const DataType& x, const DataType& y, + const DataType& z) const { + auto result = make_with_value>(x, 0.0); + DataType r_squared = square(x)+square(y)+square(z); + DataType r = sqrt(r_squared); + DataType prefactor = sqrt(1 - (square(z) / r_squared)); + + get<0, 0>(result) = x / r; + get<0, 1>(result) = x * z / (r * prefactor); + get<0, 2>(result) = -y; + get<1, 0>(result) = y / r; + get<1, 1>(result) = y * z / (r * prefactor); + get<1, 2>(result) = x; + get<2, 0>(result) = z / r; + get<2, 1>(result) = -r * prefactor; + // get<2, 2>(result) vanishes identically + + return result; +} + +template +tnsr::I +SphericalKerrSchildCoords::cartesian_from_spherical_ks( + const tnsr::I& spatial_vector, + const tnsr::I& cartesian_coords) const { + auto result = make_with_value>( + cartesian_coords, 0.0); + + for (size_t s = 0; s < get_size(get<0>(cartesian_coords)); ++s) { + const double& x = get_element(get<0>(cartesian_coords), s); + const double& y = get_element(get<1>(cartesian_coords), s); + const double& z = get_element(get<2>(cartesian_coords), s); + + // For a point off the z-axis, transforming with Jacobian is well defined. + if (LIKELY(not(equal_within_roundoff(x, 0.0) and + equal_within_roundoff(y, 0.0)))) { + const auto jacobian = jacobian_matrix(x, y, z); + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j < 3; ++j) { + get_element(result.get(i), s) += + jacobian.get(i, j) * get_element(spatial_vector.get(j), s); + } + } + } else { + // On the z-axis, the Jacobian provides the multivalued transformation + // A well defined transformation is possible if v^\theta vanishes at the + // z-axis, in which case the transformed vector points along the z-axis. + ASSERT(equal_within_roundoff(get_element(get<1>(spatial_vector), s), 0.0), + "The input vector must have a vanishing theta component on the " + "z-axis in order to perform a single-valued transformation. The " + "vector passed has v^theta = " + << get_element(get<1>(spatial_vector), s) << " at z = " << z + << " on the z-axis."); + get_element(get<2>(result), s) = + ((z > 0.0) ? get_element(get<0>(spatial_vector), s) + : -1.0 * get_element(get<0>(spatial_vector), s)); + } + } + return result; +} + +bool operator==(const SphericalKerrSchildCoords& lhs, + const SphericalKerrSchildCoords& rhs) { + return lhs.spin_a_ == rhs.spin_a_; +} + +bool operator!=(const SphericalKerrSchildCoords& lhs, + const SphericalKerrSchildCoords& rhs) { + return not(lhs == rhs); +} + +#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) + +#define INSTANTIATE(_, data) \ + template tnsr::I \ + SphericalKerrSchildCoords::cartesian_from_spherical_ks( \ + const tnsr::I& spatial_vector, \ + const tnsr::I& cartesian_coords) const; + +GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector)) + +#undef DTYPE +#undef INSTANTIATE + +} // namespace gr diff --git a/src/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.hpp b/src/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.hpp new file mode 100644 index 0000000000000..71d924498d101 --- /dev/null +++ b/src/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.hpp @@ -0,0 +1,66 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/Tensor/TypeAliases.hpp" + +// IWYU pragma: no_include + +/// \cond +namespace PUP { +class er; // IWYU pragma: keep +} // namespace PUP +/// \endcond + +namespace gr { +/*! + * \brief Contains helper functions for transforming tensors for Kerr + * spacetime in Spherical-Kerr-Schild coordinates. + */ +class SphericalKerrSchildCoords { + public: + SphericalKerrSchildCoords() = default; + SphericalKerrSchildCoords(const SphericalKerrSchildCoords& /*rhs*/) = default; + SphericalKerrSchildCoords& operator=( + const SphericalKerrSchildCoords& /*rhs*/) = default; + SphericalKerrSchildCoords(SphericalKerrSchildCoords&& /*rhs*/) = default; + SphericalKerrSchildCoords& operator=(SphericalKerrSchildCoords&& /*rhs*/) = + default; + ~SphericalKerrSchildCoords() = default; + + SphericalKerrSchildCoords(double bh_mass, double bh_dimless_spin); + + // Transforms a spatial vector from spherical + // coordinates to Cartesian for Spherical-Kerr-Schild coordinates. + // If applied on points on the z-axis, the vector to transform + // must have a vanishing \f$v^\vartheta\f$ in order + // for the transformation to be single-valued. + template + tnsr::I cartesian_from_spherical_ks( + const tnsr::I& spatial_vector, + const tnsr::I& cartesian_coords) const; + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& /*p*/); + + private: + friend bool operator==(const SphericalKerrSchildCoords& lhs, + const SphericalKerrSchildCoords& rhs); + + // The spatial components of the Jacobian of the transformation from + // spherical coordinates (r, theta, phi) to Cartesian coordinates (x, y, z) + // for Spherical-Kerr-Schild. + template + tnsr::Ij jacobian_matrix( + const DataType& x, const DataType& y, const DataType& z) const; + + double spin_a_ = std::numeric_limits::signaling_NaN(); +}; + +bool operator!=(const SphericalKerrSchildCoords& lhs, + const SphericalKerrSchildCoords& rhs); + +} // namespace gr diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py index 2797c3e0e6529..993cd4fcc1520 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py @@ -153,7 +153,7 @@ def magnetic_field( polytropic_constant, polytropic_exponent, ) - x_max = np.array([r_max, spin_a, 0.0]) + x_max = np.array([r_max, 0.0, 0.0]) threshold_rho = threshold_density * ( fm_disk.rest_mass_density( x_max, @@ -166,13 +166,15 @@ def magnetic_field( polytropic_exponent, ) ) + # FIXME: temporary x in Kerr Schild coordinates. + x_ks = [0.0, 0.0, 0.0] + if rho > threshold_rho: small = 0.0001 * bh_mass a_squared = spin_a**2 sin_theta_squared = x[0] ** 2 + x[1] ** 2 - r_squared = 0.5 * (sin_theta_squared + x[2] ** 2 - a_squared) - r_squared += np.sqrt(r_squared**2 + a_squared * x[2] ** 2) - sin_theta_squared /= r_squared + a_squared + r_squared = sin_theta_squared + x[2] ** 2 + sin_theta_squared /= r_squared r = np.sqrt(r_squared) sin_theta = np.sqrt(sin_theta_squared) @@ -233,13 +235,21 @@ def magnetic_field( threshold_rho, ) ) * prefactor - - # This normalization is specific for the default normalization resolution - # grid and for the specific member variables of the disk in test_variables - normalization = 1.7162625566578704 - return normalization * ks_coords.cartesian_from_spherical_ks( - result, x, bh_mass, bh_dimless_spin - ) + sks_to_ks_factor = np.sqrt(r_squared + a_squared) / np.sqrt(r_squared) + x_ks[0] = x[0] * sks_to_ks_factor + x_ks[1] = x[1] * sks_to_ks_factor + x_ks[2] = x[2] + # This normalization is specific for the lower (than default) + # normalization resolution for a faster testing. + # normalization = 1.71896376101037762218 # 51 before the change + normalization = 1.3733086119266160185503622 # 51 after the change + # normalization = 1.3686080052126354811292686 # 255 after the change. + result = normalization * ks_coords.cartesian_from_spherical_ks( + result, x_ks, bh_mass, bh_dimless_spin + ) + inv_jac = fm_disk.inverse_jacobian_matrix(x, spin_a) + result = np.matmul(inv_jac.T, result) + return result def divergence_cleaning_field( diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp index f105d9e5eba6a..af2b5e5e5771f 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp @@ -19,7 +19,7 @@ #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp" #include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp" #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp" -#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" @@ -84,7 +84,7 @@ void test_create_from_options() { " PolytropicExponent: 1.654\n" " ThresholdDensity: 0.42\n" " InversePlasmaBeta: 85.0\n" - " BFieldNormGridRes: 6") + " BFieldNormGridRes: 4") ->get_clone(); const auto deserialized_option_solution = serialize_and_deserialize(option_solution); @@ -92,7 +92,7 @@ void test_create_from_options() { *deserialized_option_solution); CHECK(disk == grmhd::AnalyticData::MagnetizedFmDisk( - 1.3, 0.345, 6.123, 14.2, 0.065, 1.654, 0.42, 85.0, 6)); + 1.3, 0.345, 6.123, 14.2, 0.065, 1.654, 0.42, 85.0, 4)); } void test_move() { @@ -111,7 +111,8 @@ void test_serialize() { template void test_variables(const DataType& used_for_size) { - const double bh_mass = 1.12; + // const double bh_mass = 1.12; + const double bh_mass = 1.; const double bh_dimless_spin = 0.97; const double inner_edge_radius = 6.2; const double max_pressure_radius = 11.6; @@ -119,11 +120,12 @@ void test_variables(const DataType& used_for_size) { const double polytropic_exponent = 1.65; const double threshold_density = 0.14; const double inverse_plasma_beta = 0.023; - + const size_t b_field_normalization = 51; // Using lower than default + // resolution for faster testing. MagnetizedFmDiskProxy disk(bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, polytropic_constant, polytropic_exponent, threshold_density, - inverse_plasma_beta); + inverse_plasma_beta, b_field_normalization); const auto member_variables = std::make_tuple( bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, polytropic_constant, polytropic_exponent, threshold_density, @@ -137,24 +139,23 @@ void test_variables(const DataType& used_for_size) { "divergence_cleaning_field"}, {{{-20., 20.}}}, member_variables, used_for_size, 1.0e-8); - // Test a few of the GR components to make sure that the implementation - // correctly forwards to the background solution of the base class. + // // Test a few of the GR components to make sure that the implementation + // // correctly forwards to the background solution of the base class. const auto coords = make_with_value>(used_for_size, 1.0); - const gr::Solutions::KerrSchild ks_soln{ + const gr::Solutions::SphericalKerrSchild sks_soln{ bh_mass, {{0.0, 0.0, bh_dimless_spin}}, {{0.0, 0.0, 0.0}}}; - CHECK_ITERABLE_APPROX( - get>(ks_soln.variables( - coords, 0.0, gr::Solutions::KerrSchild::tags{})), - get>( - disk.variables(coords, tmpl::list>{}))); - CHECK_ITERABLE_APPROX( - get>(ks_soln.variables( - coords, 0.0, gr::Solutions::KerrSchild::tags{})), - get>(disk.variables( - coords, tmpl::list>{}))); + const auto [expected_lapse, expected_sqrt_det_gamma] = sks_soln.variables( + coords, 0.0, + tmpl::list, + gr::Tags::SqrtDetSpatialMetric>{}); + const auto [lapse, sqrt_det_gamma] = disk.variables( + coords, tmpl::list, + gr::Tags::SqrtDetSpatialMetric>{}); + CHECK_ITERABLE_APPROX(expected_lapse, lapse); + CHECK_ITERABLE_APPROX(expected_sqrt_det_gamma, sqrt_det_gamma); const auto expected_spatial_metric = - get>(ks_soln.variables( - coords, 0.0, gr::Solutions::KerrSchild::tags{})); + get>(sks_soln.variables( + coords, 0.0, gr::Solutions::SphericalKerrSchild::tags{})); const auto spatial_metric = get>(disk.variables( coords, tmpl::list>{})); @@ -172,18 +173,16 @@ void test_variables(const DataType& used_for_size) { {"rest_mass_density", "spatial_velocity", "specific_internal_energy", "pressure", "lorentz_factor", "specific_enthalpy"}, {{{-20., 20.}}}, member_variables, used_for_size, 1.0e-8); - const auto magnetic_field = - get>(another_disk.variables( - coords, tmpl::list>{})); + const auto [magnetic_field, div_clean_field] = + another_disk.variables( + coords, tmpl::list, + hydro::Tags::DivergenceCleaningField>{}); const auto expected_magnetic_field = make_with_value>(used_for_size, 0.0); + const auto expected_div_clean_field = + make_with_value>(used_for_size, 0.0); CHECK_ITERABLE_APPROX(magnetic_field, expected_magnetic_field); - CHECK_ITERABLE_APPROX( - get>( - another_disk.variables( - coords, - tmpl::list>{})), - make_with_value>(used_for_size, 0.0)); + CHECK_ITERABLE_APPROX(div_clean_field, expected_div_clean_field); } } // namespace @@ -199,93 +198,95 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.MagFmDisk", test_variables(std::numeric_limits::signaling_NaN()); test_variables(DataVector(5)); -#ifdef SPECTRE_DEBUG - CHECK_THROWS_WITH( - []() { - grmhd::AnalyticData::MagnetizedFmDisk disk( - 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, -0.2, 0.023, 4); - }(), - Catch::Matchers::ContainsSubstring( - "The threshold density must be in the range (0, 1)")); - CHECK_THROWS_WITH( - []() { - grmhd::AnalyticData::MagnetizedFmDisk disk( - 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 1.45, 0.023, 4); - }(), - Catch::Matchers::ContainsSubstring( - "The threshold density must be in the range (0, 1)")); - CHECK_THROWS_WITH( - []() { - grmhd::AnalyticData::MagnetizedFmDisk disk( - 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 0.2, -0.153, 4); - }(), - Catch::Matchers::ContainsSubstring( - "The inverse plasma beta must be non-negative.")); - CHECK_THROWS_WITH( - []() { - grmhd::AnalyticData::MagnetizedFmDisk disk(0.7, 0.61, 5.4, 9.182, - 11.123, 1.44, 0.2, 0.153, 2); - }(), - Catch::Matchers::ContainsSubstring( - "The grid resolution used in the magnetic field " - "normalization must be at least 4 points.")); - CHECK_THROWS_WITH( - []() { - grmhd::AnalyticData::MagnetizedFmDisk disk(0.7, 0.61, 5.4, 9.182, - 11.123, 1.44, 0.2, 0.023, 5); - }(), - Catch::Matchers::ContainsSubstring("Max b squared is zero.")); -#endif - CHECK_THROWS_WITH( - TestHelpers::test_creation( - "BhMass: 13.45\n" - "BhDimlessSpin: 0.45\n" - "InnerEdgeRadius: 6.1\n" - "MaxPressureRadius: 7.6\n" - "PolytropicConstant: 2.42\n" - "PolytropicExponent: 1.33\n" - "ThresholdDensity: -0.01\n" - "InversePlasmaBeta: 0.016\n" - "BFieldNormGridRes: 4"), - Catch::Matchers::ContainsSubstring( - "Value -0.01 is below the lower bound of 0")); - CHECK_THROWS_WITH( - TestHelpers::test_creation( - "BhMass: 1.5\n" - "BhDimlessSpin: 0.94\n" - "InnerEdgeRadius: 6.4\n" - "MaxPressureRadius: 8.2\n" - "PolytropicConstant: 41.1\n" - "PolytropicExponent: 1.8\n" - "ThresholdDensity: 4.1\n" - "InversePlasmaBeta: 0.03\n" - "BFieldNormGridRes: 4"), - Catch::Matchers::ContainsSubstring( - "Value 4.1 is above the upper bound of 1")); - CHECK_THROWS_WITH( - TestHelpers::test_creation( - "BhMass: 1.4\n" - "BhDimlessSpin: 0.91\n" - "InnerEdgeRadius: 6.5\n" - "MaxPressureRadius: 7.8\n" - "PolytropicConstant: 13.5\n" - "PolytropicExponent: 1.54\n" - "ThresholdDensity: 0.22\n" - "InversePlasmaBeta: -0.03\n" - "BFieldNormGridRes: 4"), - Catch::Matchers::ContainsSubstring( - "Value -0.03 is below the lower bound of 0")); - CHECK_THROWS_WITH( - TestHelpers::test_creation( - "BhMass: 1.4\n" - "BhDimlessSpin: 0.91\n" - "InnerEdgeRadius: 6.5\n" - "MaxPressureRadius: 7.8\n" - "PolytropicConstant: 13.5\n" - "PolytropicExponent: 1.54\n" - "ThresholdDensity: 0.22\n" - "InversePlasmaBeta: 0.03\n" - "BFieldNormGridRes: 2"), - Catch::Matchers::ContainsSubstring( - "Value 2 is below the lower bound of 4")); + #ifdef SPECTRE_DEBUG + CHECK_THROWS_WITH( + []() { + grmhd::AnalyticData::MagnetizedFmDisk disk( + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, -0.2, 0.023, 4); + }(), + Catch::Matchers::ContainsSubstring( + "The threshold density must be in the range (0, 1)")); + CHECK_THROWS_WITH( + []() { + grmhd::AnalyticData::MagnetizedFmDisk disk( + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 1.45, 0.023, 4); + }(), + Catch::Matchers::ContainsSubstring( + "The threshold density must be in the range (0, 1)")); + CHECK_THROWS_WITH( + []() { + grmhd::AnalyticData::MagnetizedFmDisk disk( + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 0.2, -0.153, 4); + }(), + Catch::Matchers::ContainsSubstring( + "The inverse plasma beta must be non-negative.")); + CHECK_THROWS_WITH( + []() { + grmhd::AnalyticData::MagnetizedFmDisk disk(0.7, 0.61, 5.4, 9.182, + 11.123, 1.44, 0.2, + 0.153, 2); + }(), + Catch::Matchers::ContainsSubstring( + "The grid resolution used in the magnetic field " + "normalization must be at least 4 points.")); + CHECK_THROWS_WITH( + []() { + grmhd::AnalyticData::MagnetizedFmDisk disk(0.7, 0.61, 5.4, 9.182, + 11.123, 1.44, 0.2, + 0.023, 4); + }(), + Catch::Matchers::ContainsSubstring("Max b squared is zero.")); + #endif + CHECK_THROWS_WITH( + TestHelpers::test_creation( + "BhMass: 13.45\n" + "BhDimlessSpin: 0.45\n" + "InnerEdgeRadius: 6.1\n" + "MaxPressureRadius: 7.6\n" + "PolytropicConstant: 2.42\n" + "PolytropicExponent: 1.33\n" + "ThresholdDensity: -0.01\n" + "InversePlasmaBeta: 0.016\n" + "BFieldNormGridRes: 4"), + Catch::Matchers::ContainsSubstring( + "Value -0.01 is below the lower bound of 0")); + CHECK_THROWS_WITH( + TestHelpers::test_creation( + "BhMass: 1.5\n" + "BhDimlessSpin: 0.94\n" + "InnerEdgeRadius: 6.4\n" + "MaxPressureRadius: 8.2\n" + "PolytropicConstant: 41.1\n" + "PolytropicExponent: 1.8\n" + "ThresholdDensity: 4.1\n" + "InversePlasmaBeta: 0.03\n" + "BFieldNormGridRes: 4"), + Catch::Matchers::ContainsSubstring( + "Value 4.1 is above the upper bound of 1")); + CHECK_THROWS_WITH( + TestHelpers::test_creation( + "BhMass: 1.4\n" + "BhDimlessSpin: 0.91\n" + "InnerEdgeRadius: 6.5\n" + "MaxPressureRadius: 7.8\n" + "PolytropicConstant: 13.5\n" + "PolytropicExponent: 1.54\n" + "ThresholdDensity: 0.22\n" + "InversePlasmaBeta: -0.03\n" + "BFieldNormGridRes: 4"), + Catch::Matchers::ContainsSubstring( + "Value -0.03 is below the lower bound of 0")); + CHECK_THROWS_WITH( + TestHelpers::test_creation( + "BhMass: 1.4\n" + "BhDimlessSpin: 0.91\n" + "InnerEdgeRadius: 6.5\n" + "MaxPressureRadius: 7.8\n" + "PolytropicConstant: 13.5\n" + "PolytropicExponent: 1.54\n" + "ThresholdDensity: 0.22\n" + "InversePlasmaBeta: 0.03\n" + "BFieldNormGridRes: 2"), + Catch::Matchers::ContainsSubstring( + "Value 2 is below the lower bound of 4")); } diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp index f58968ac30225..dc4fa4b198f0a 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp @@ -49,7 +49,7 @@ void test_create_from_options() { " PolytropicExponent: 1.654\n" " ThresholdDensity: 0.42\n" " InversePlasmaBeta: 85.0\n" - " BFieldNormGridRes: 6\n" + " BFieldNormGridRes: 4\n" " TorusParameters:\n" " RadialRange: [2.5, 3.6]\n" " MinPolarAngle: 0.9\n" @@ -62,18 +62,18 @@ void test_create_from_options() { *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), + 1.3, 0.345, 6.123, 14.2, 0.065, 1.654, 0.42, 85.0, 4), 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), + 1.654, 0.42, 85.0, 4), 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), + 1.654, 0.42, 85.0, 4), grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7)); test_move_semantics(std::move(disk), disk_copy); } @@ -81,7 +81,7 @@ void test_move() { 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), + 1.654, 0.42, 85.0, 4), grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7)); test_serialization(disk); } @@ -105,12 +105,13 @@ void test_variables(const DataType& used_for_size) { const double polytropic_exponent = 1.65; const double threshold_density = 0.14; const double inverse_plasma_beta = 0.023; + const size_t b_field_normalization = 51; 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), + inverse_plasma_beta, b_field_normalization), 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/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py index 685d58a929eca..df414b9d05909 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py @@ -8,6 +8,46 @@ def delta(r_sqrd, m, a): return r_sqrd - 2.0 * m * np.sqrt(r_sqrd) + a**2 +def transformation_matrix(x, a): + # Coordinate transformation matrix for KS to Spherical KS. + # Returns the P^j_ibar from equation 10 of Spherical-KS documentation. + # Note, this assuemes the black hole spin is along z-direction. + r_sqrd = boyer_lindquist_r_sqrd(x) + r = np.sqrt(r_sqrd) + rho = np.sqrt(r_sqrd + a**2) + P_ij = np.diag([rho / r, rho / r, 1.0]) + return P_ij + + +def jacobian_matrix(x, a): + # Jacbian matrix for coordinate transformation from KS to Spherical KS. + # Returns T^i_jbar from equation 16 of Spherical-KS documentation. + # Note, this assuemes the black hole spin is along z-direction. + P_ij = transformation_matrix(x, a) + + r_sqrd = boyer_lindquist_r_sqrd(x) + r = np.sqrt(r_sqrd) + rho = np.sqrt(r_sqrd + a**2) + + F_ij = (-1.0 / (rho * r**3)) * np.diag([a**2, a**2, 0.0]) + x_vec = np.array(x).T + + T_ij = P_ij + np.outer(np.matmul(F_ij, x_vec), x_vec.T) + + return T_ij.T + + +def inverse_jacobian_matrix(x, a): + # Inverse Jacbian matrix for coordinate transformation + # from KS to Spherical KS. + # Returns S^i_jbar from equation 17 of Spherical-KS documentation. + # Note, this assuemes the black hole spin is along z-direction. + T_ij = jacobian_matrix(x, a) + S_ij = np.linalg.inv(T_ij) + + return S_ij.T + + def sigma(r_sqrd, sin_theta_sqrd, a): return r_sqrd + (1.0 - sin_theta_sqrd) * a**2 @@ -37,9 +77,9 @@ def boyer_lindquist_gff(r_sqrd, sin_theta_sqrd, m, a): ) -def boyer_lindquist_r_sqrd(x, a): - half_diff = 0.5 * (x[0] ** 2 + x[1] ** 2 + x[2] ** 2 - a**2) - return half_diff + np.sqrt(half_diff**2 + a**2 * x[2] ** 2) +def boyer_lindquist_r_sqrd(x): + # Here, we assume x is in Spherical-KS coordinates + return x[0] ** 2 + x[1] ** 2 + x[2] ** 2 def boyer_lindquist_sin_theta_sqrd(z_sqrd, r_sqrd): @@ -47,17 +87,50 @@ def boyer_lindquist_sin_theta_sqrd(z_sqrd, r_sqrd): def kerr_schild_h(x, m, a): - r_sqrd = boyer_lindquist_r_sqrd(x, a) + r_sqrd = boyer_lindquist_r_sqrd(x) return m * r_sqrd * np.sqrt(r_sqrd) / (r_sqrd**2 + a**2 * x[2] ** 2) def kerr_schild_spatial_null_form(x, m, a): - r_sqrd = boyer_lindquist_r_sqrd(x, a) + r_sqrd = boyer_lindquist_r_sqrd(x) r = np.sqrt(r_sqrd) - denom = 1.0 / (r_sqrd + a**2) - return np.array( - [(r * x[0] + a * x[1]) * denom, (r * x[1] - a * x[0]) * denom, x[2] / r] - ) + rho = np.sqrt(r_sqrd + a**2) + denom = 1.0 / rho**2 + # Again, we assume Spherical-KS coordinates + # xbar/r = x/rho, ybar/r = y/rho + # where rho^2 = r^2+a^2 and xbar and ybar are Spherical-KS + # Thus, we need to stick in the converseion factor of rho/r. + conv_fac = rho / r + l_vec = np.array( + [ + (r * x[0] * conv_fac + a * x[1] * conv_fac) * denom, + (r * x[1] * conv_fac - a * x[0] * conv_fac) * denom, + x[2] / r, + ] + ).T + jac = jacobian_matrix(x, a) + return (np.matmul(jac, l_vec.T)).T + + +def kerr_schild_spatial_null_vec(x, m, a): + r_sqrd = boyer_lindquist_r_sqrd(x) + r = np.sqrt(r_sqrd) + rho = np.sqrt(r_sqrd + a**2) + denom = 1.0 / rho**2 + # Again, we assume Spherical-KS coordinates + # xbar/r = x/rho, ybar/r = y/rho + # where rho^2 = r^2+a^2 and xbar and ybar are Spherical-KS + # Thus, we need to stick in the converseion factor of rho/r. + conv_fac = rho / r + l_vec = np.array( + [ + (r * x[0] * conv_fac + a * x[1] * conv_fac) * denom, + (r * x[1] * conv_fac - a * x[0] * conv_fac) * denom, + x[2] / r, + ] + ).T + inv_jac = inverse_jacobian_matrix(x, a) + return (np.matmul(inv_jac, l_vec)).T def kerr_schild_lapse(x, m, a): @@ -75,13 +148,17 @@ def kerr_schild_shift(x, m, a): * kerr_schild_h(x, m, a) * null_vector_0 * kerr_schild_lapse(x, m, a) ** 2 - ) * kerr_schild_spatial_null_form(x, m, a) + ) * kerr_schild_spatial_null_vec(x, m, a) def kerr_schild_spatial_metric(x, m, a): prefactor = 2.0 * kerr_schild_h(x, m, a) null_form = kerr_schild_spatial_null_form(x, m, a) - return np.identity(x.size) + prefactor * np.outer(null_form, null_form) + T_ij = jacobian_matrix(x, a) + result = np.matmul(T_ij, T_ij.T) + prefactor * np.outer( + null_form.T, null_form.T + ) + return result def angular_momentum(m, a, rmax): @@ -142,7 +219,7 @@ def specific_enthalpy( bh_spin_a = bh_mass * bh_dimless_a l = angular_momentum(bh_mass, bh_spin_a, bh_mass * dimless_r_max) Win = potential(l, r_in * r_in, 1.0, bh_mass, bh_spin_a) - r_sqrd = boyer_lindquist_r_sqrd(x, bh_spin_a) + r_sqrd = boyer_lindquist_r_sqrd(x) sin_theta_sqrd = boyer_lindquist_sin_theta_sqrd(x[2] * x[2], r_sqrd) result = 1.0 if np.sqrt(r_sqrd * sin_theta_sqrd) >= r_in: @@ -251,18 +328,26 @@ def spatial_velocity( bh_spin_a = bh_mass * bh_dimless_a l = angular_momentum(bh_mass, bh_spin_a, bh_mass * dimless_r_max) Win = potential(l, r_in * r_in, 1.0, bh_mass, bh_spin_a) - r_sqrd = boyer_lindquist_r_sqrd(x, bh_spin_a) + r_sqrd = boyer_lindquist_r_sqrd(x) sin_theta_sqrd = boyer_lindquist_sin_theta_sqrd(x[2] * x[2], r_sqrd) + sks_to_ks_factor = np.sqrt(r_sqrd + bh_spin_a**2) / np.sqrt(r_sqrd) result = np.array([0.0, 0.0, 0.0]) if np.sqrt(r_sqrd * sin_theta_sqrd) >= r_in: W = potential(l, r_sqrd, sin_theta_sqrd, bh_mass, bh_spin_a) if W < Win: - result += ( - np.array([-x[1], x[0], 0.0]) + transport_velocity_ks = ( + sks_to_ks_factor + * np.array([-x[1], x[0], 0.0]) * angular_velocity( l, r_sqrd, sin_theta_sqrd, bh_mass, bh_spin_a ) + ) + transport_velocity_sks = np.matmul( + inverse_jacobian_matrix(x, bh_spin_a), transport_velocity_ks.T + ) + result += ( + transport_velocity_sks + kerr_schild_shift(x, bh_mass, bh_spin_a) ) / kerr_schild_lapse(x, bh_mass, bh_spin_a) return result diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp index a0749b237b8f7..30490094ffbdd 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp @@ -24,7 +24,7 @@ #include "NumericalAlgorithms/Spectral/Basis.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "NumericalAlgorithms/Spectral/Quadrature.hpp" -#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp" #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -101,7 +101,7 @@ void test_serialize() { template void test_variables(const DataType& used_for_size) { - const double bh_mass = 1.34; + const double bh_mass = 1.0; const double bh_dimless_spin = 0.94432; const double inner_edge_radius = 6.0; const double max_pressure_radius = 12.0; @@ -133,22 +133,22 @@ void test_variables(const DataType& used_for_size) { // Test a few of the GR components to make sure that the implementation // correctly forwards to the background solution. Not meant to be extensive. const auto coords = make_with_value>(used_for_size, 1.0); - const gr::Solutions::KerrSchild ks_soln{ + const gr::Solutions::SphericalKerrSchild sks_soln{ bh_mass, {{0.0, 0.0, bh_dimless_spin}}, {{0.0, 0.0, 0.0}}}; CHECK_ITERABLE_APPROX( - get>(ks_soln.variables( - coords, 0.0, gr::Solutions::KerrSchild::tags{})), + get>(sks_soln.variables( + coords, 0.0, gr::Solutions::SphericalKerrSchild::tags{})), get>(disk.variables( coords, 0.0, tmpl::list>{}))); CHECK_ITERABLE_APPROX( - get>(ks_soln.variables( - coords, 0.0, gr::Solutions::KerrSchild::tags{})), + get>(sks_soln.variables( + coords, 0.0, gr::Solutions::SphericalKerrSchild::tags{})), get>(disk.variables( coords, 0.0, tmpl::list>{}))); const auto expected_spatial_metric = - get>(ks_soln.variables( - coords, 0.0, gr::Solutions::KerrSchild::tags{})); + get>(sks_soln.variables( + coords, 0.0, gr::Solutions::SphericalKerrSchild::tags{})); const auto spatial_metric = get>(disk.variables( coords, 0.0, tmpl::list>{})); diff --git a/tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt b/tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt index e88223e37a5b1..070d8c2b2c8f0 100644 --- a/tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt +++ b/tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt @@ -13,6 +13,7 @@ set(LIBRARY_SOURCES Test_Psi4.cpp Test_Ricci.cpp Test_SpacetimeDerivativeOfGothG.cpp + Test_SphericalKerrSchildCoords.cpp Test_Tags.cpp Test_WeylElectric.cpp Test_WeylMagnetic.cpp diff --git a/tests/Unit/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.py b/tests/Unit/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.py new file mode 100644 index 0000000000000..58fac2e79a84a --- /dev/null +++ b/tests/Unit/PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.py @@ -0,0 +1,26 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import numpy as np + + +def jacobian(coords): + result = np.zeros((3, 3)) + x, y, z = coords + r_squared = x**2 + y**2 + z**2 + r = np.sqrt(r_squared) + prefactor = np.sqrt(1 - ((z**2) / r_squared)) + + result[0, 0] = x / r + result[0, 1] = x * z / (r * prefactor) + result[0, 2] = -y + result[1, 0] = y / r + result[1, 1] = y * z / (r * prefactor) + result[1, 2] = x + result[2, 0] = z / r + result[2, 1] = -r * prefactor + return result + + +def cartesian_from_spherical_ks(vector, cartesian_coords): + return np.einsum("ij,j", jacobian(cartesian_coords), vector) diff --git a/tests/Unit/PointwiseFunctions/GeneralRelativity/Test_SphericalKerrSchildCoords.cpp b/tests/Unit/PointwiseFunctions/GeneralRelativity/Test_SphericalKerrSchildCoords.cpp new file mode 100644 index 0000000000000..3fff3b01c4e6b --- /dev/null +++ b/tests/Unit/PointwiseFunctions/GeneralRelativity/Test_SphericalKerrSchildCoords.cpp @@ -0,0 +1,169 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" // IWYU pragma: keep +#include "Framework/CheckWithRandomValues.hpp" +#include "Framework/Pypp.hpp" +#include "Framework/SetupLocalPythonEnvironment.hpp" +#include "Framework/TestHelpers.hpp" +#include "PointwiseFunctions/GeneralRelativity/SphericalKerrSchildCoords.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/MakeWithValue.hpp" + +namespace { + +template +void test_kerr_schild_coords(const DataType& used_for_size) { + const double bh_mass = 2.54; + const double bh_dimless_spin = 0.9375; + + gr::SphericalKerrSchildCoords sph_kerr_schild_coords(bh_mass, + bh_dimless_spin); + + pypp::check_with_random_values<1>( + &gr::SphericalKerrSchildCoords::cartesian_from_spherical_ks, + sph_kerr_schild_coords, "SphericalKerrSchildCoords", + "cartesian_from_spherical_ks", {{{-10.0, 10.0}}}, {}, used_for_size); +} + +template +void test_coord_transformation_on_xy_plane(const DataType& used_for_size) { + const double bh_mass = 0.34; + const double bh_dimless_spin = 0.6332; + const double spin_a = bh_mass * bh_dimless_spin; + + MAKE_GENERATOR(generator); + + // get random point outside of the ring singularity x^2 + y^2 = a^2 + // to avoid possible FPE when computing the Jacobian. + std::uniform_real_distribution<> distribution_mod(spin_a, 10); + const double random_mod = distribution_mod(generator); + std::uniform_real_distribution<> distribution_angle(0, 2.0 * M_PI); + const double random_angle = distribution_angle(generator); + auto random_point_on_xy_plane = + make_with_value>(used_for_size, 0.0); + get<0>(random_point_on_xy_plane) = random_mod * cos(random_angle); + get<1>(random_point_on_xy_plane) = random_mod * sin(random_angle); + + std::uniform_real_distribution<> distribution_vector(-3.0, 3.0); + auto random_vector_to_transform = + make_with_value>(used_for_size, 0.0); + for (size_t i = 0; i < 3; ++i) { + random_vector_to_transform.get(i) = distribution_vector(generator); + } + + CHECK_ITERABLE_APPROX( + (gr::SphericalKerrSchildCoords{bh_mass, bh_dimless_spin} + .cartesian_from_spherical_ks(random_vector_to_transform, + random_point_on_xy_plane)), + (pypp::call>( + "SphericalKerrSchildCoords", "cartesian_from_spherical_ks", + random_vector_to_transform, random_point_on_xy_plane))); +} + +template +void test_coord_transformation_along_z_axis(const DataType& used_for_size) { + MAKE_GENERATOR(generator); + std::uniform_real_distribution<> distribution(-13.0, 13.0); + + auto random_radial_vector_to_transform = + make_with_value>(used_for_size, 0.0); + get<0>(random_radial_vector_to_transform) = distribution(generator); + + const double z_coordinate = distribution(generator); + auto random_point_along_z_axis = + make_with_value>(used_for_size, 0.0); + get<2>(random_point_along_z_axis) = z_coordinate; + + gr::SphericalKerrSchildCoords sks_coords{1.532, 0.9375}; + + // New vector should point along z-axis, with v^z = sign(z) v^r + const auto transformed_vector = sks_coords.cartesian_from_spherical_ks( + random_radial_vector_to_transform, random_point_along_z_axis); + CHECK(get<0>(transformed_vector) == 0.0); + CHECK(get<1>(transformed_vector) == 0.0); + CHECK(get<2>(transformed_vector) == + (z_coordinate > 0.0 ? 1.0 : -1.0) * + get<0>(random_radial_vector_to_transform)); + + // Since we transform a uniform vector, the transformation should not modify + // the new vector upon flipping the sign of the z coordinate + get<2>(random_point_along_z_axis) *= -1.0; + // v^r flips sign in order to represent the same uniform vector at z < 0 + get<0>(random_radial_vector_to_transform) *= -1.0; + CHECK_ITERABLE_APPROX( + transformed_vector, + sks_coords.cartesian_from_spherical_ks(random_radial_vector_to_transform, + random_point_along_z_axis)); +} + +template +void test_theta_component(const DataType& used_for_size) { + MAKE_GENERATOR(generator); + std::uniform_real_distribution<> distribution(-13.0, 13.0); + + auto random_point_along_z_axis = + make_with_value>(used_for_size, 0.0); + get<2>(random_point_along_z_axis) = distribution(generator); + + // input vector has a nonvanishing theta component along the z-axis. + gr::SphericalKerrSchildCoords{3.123, 0.854}.cartesian_from_spherical_ks( + make_with_value>( + used_for_size, distribution(generator)), + random_point_along_z_axis); +} + +} // namespace + +SPECTRE_TEST_CASE( + "Unit.PointwiseFunctions.GeneralRelativity.SphericalKerrSchildCoords", + "[PointwiseFunctions][Unit]") { + pypp::SetupLocalPythonEnvironment local_python_env( + "PointwiseFunctions/GeneralRelativity/"); + + const double d(std::numeric_limits::signaling_NaN()); + const DataVector dv(5); + test_kerr_schild_coords(d); + test_kerr_schild_coords(dv); + test_coord_transformation_along_z_axis(d); + test_coord_transformation_along_z_axis(dv); + test_coord_transformation_on_xy_plane(d); + test_coord_transformation_on_xy_plane(dv); + + gr::SphericalKerrSchildCoords sph_kerr_schild_coords(0.8624, 0.151); + test_serialization(sph_kerr_schild_coords); + + gr::SphericalKerrSchildCoords sph_kerr_schild_coords_copy(0.8624, 0.151); + test_move_semantics(std::move(sph_kerr_schild_coords), + sph_kerr_schild_coords_copy); // NOLINT + +#ifdef SPECTRE_DEBUG + CHECK_THROWS_WITH( + (test_theta_component(std::numeric_limits::signaling_NaN())), + Catch::Matchers::ContainsSubstring( + "The input vector must have a vanishing theta component")); + CHECK_THROWS_WITH( + (test_theta_component(DataVector(5))), + Catch::Matchers::ContainsSubstring( + "The input vector must have a vanishing theta component")); + CHECK_THROWS_WITH( + (gr::SphericalKerrSchildCoords(-4.21, 0.999)), + Catch::Matchers::ContainsSubstring("The mass must be positive")); + CHECK_THROWS_WITH((gr::SphericalKerrSchildCoords(0.15, -1.3)), + Catch::Matchers::ContainsSubstring( + "The dimensionless spin must be in the range (-1, 1)")); + CHECK_THROWS_WITH((gr::SphericalKerrSchildCoords(1.532, 4.2)), + Catch::Matchers::ContainsSubstring( + "The dimensionless spin must be in the range (-1, 1)")); +#endif +}