From 15c9b2ce98933331cfce3acdedd2a0ed6aa11d62 Mon Sep 17 00:00:00 2001 From: Jooheon Yoo Date: Wed, 25 Oct 2023 11:32:07 -0700 Subject: [PATCH] Replace KerrSchild with SphericalKerrSchild as background spacetime for FM_disk --- src/DataStructures/CachedTempBuffer.hpp | 1 + src/DataStructures/TempBuffer.hpp | 2 + .../AnalyticData/GrMhd/MagnetizedFmDisk.cpp | 80 +++-- .../AnalyticData/GrMhd/MagnetizedFmDisk.hpp | 69 ++-- .../GeneralRelativity/SphericalKerrSchild.hpp | 55 ++- .../FishboneMoncriefDisk.cpp | 320 +++++++++--------- .../FishboneMoncriefDisk.hpp | 176 ++++------ .../AnalyticData/GrMhd/MagnetizedFmDisk.py | 31 +- .../GrMhd/Test_MagnetizedFmDisk.cpp | 234 ++++++------- .../GrMhd/Test_PolarMagnetizedFmDisk.cpp | 13 +- .../RelativisticEuler/FishboneMoncriefDisk.py | 114 +++++-- .../Test_FishboneMoncriefDisk.cpp | 18 +- 12 files changed, 593 insertions(+), 520 deletions(-) diff --git a/src/DataStructures/CachedTempBuffer.hpp b/src/DataStructures/CachedTempBuffer.hpp index 5574344e40a3..dcc437bca4b4 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 72d21a7df757..7b9912180938 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 2b8bb09b6b46..0f5d068771e0 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,53 @@ 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); + } + } + // 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 j = 0; j < 3; ++j) { + for (size_t i = 0; i < 3; ++i) { + magnetic_field_sks.get(j) += + inv_jacobians.get(j, i) * magnetic_field_ks.get(i); } } - return kerr_schild_coords_.cartesian_from_spherical_ks( - std::move(magnetic_field), x); + + 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 +284,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 819ed2a119d7..e97e068eb767 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp @@ -70,6 +70,21 @@ namespace grmhd::AnalyticData { * * is the norm of the magnetic field in the fluid frame, with \f$v_i\f$ being * the spatial velocity, and \f$W\f$ the Lorentz factor. + * + * \note When using Kerr-Schild coordinates, the horizon that is at + * constant \f$r\f$ is not spherical, but instead spheroidal. This could make + * application of boundary condition and computing various fluxes + * across the horizon more complicated than they need to be. + * Thus, similar to RelativisticEuler::Solutions::FishboneMoncriefDisk + * we use Spherical Kerr-Schild coordinates, + * see gr::Solutions::SphericalKerrSchild, in which constant \f$r\f$ + * is spherical. Because we compute variables in Kerr-Schild coordinates, + * there is a necessary extra step of transforming them back to + * Spherical Kerr-Schild coordinates. + * + * \warning Spherical Kerr-Schild coordinates and "spherical Kerr-Schild" + * coordinates are not same. + * */ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData, public MarkAsAnalyticData { @@ -153,38 +168,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 +190,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 +208,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 bb055182ffee..38faa043f9f8 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,46 @@ 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>; + + 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 7d0cdc555263..4cc5eb195c10 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 28bee599c04a..bdd9cfc3ff60 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 @@ -154,6 +151,17 @@ namespace Solutions { * \note Kozlowski et al. \cite Kozlowski1978aa denote * \f$l_* = u_\phi u^t\f$ in order to * distinguish this quantity from their own definition \f$l = - u_\phi/u_t\f$. + * + * \note When using Kerr-Schild coordinates, the horizon that is at + * constant \f$r\f$ is not spherical, but instead spheroidal. This could make + * application of boundary condition and computing various fluxes + * across the horizon more complicated than they need to be. + * Thus, we use Spherical Kerr-Schild coordinates, + * see gr::Solutions::SphericalKerrSchild, in which constant \f$r\f$ + * is spherical. Because we compute variables in Kerr-Schild coordinates, + * there is a necessary extra step of transforming them back to + * Spherical Kerr-Schild coordinates. + * */ class FishboneMoncriefDisk : public virtual evolution::initial_data::InitialData, @@ -161,7 +169,7 @@ class FishboneMoncriefDisk public AnalyticSolution<3>, public hydro::TemperatureInitialization { protected: - template + template struct IntermediateVariables; public: @@ -258,64 +266,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 +303,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 +413,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/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py index 2797c3e0e652..eb69ec600367 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,14 @@ def magnetic_field( polytropic_exponent, ) ) + 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 +234,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 disk parameters used + # in test_variables function of Test_MagnetizedFmDisk.cpp. + # Note: we are using lower (than default) + # b_field_normalization for a faster testing time. + normalization = 2.077412435947072 + 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 = inv_jac @ 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 f105d9e5eba6..f3e942d78f33 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,7 @@ void test_serialize() { template void test_variables(const DataType& used_for_size) { - const double bh_mass = 1.12; + const double bh_mass = 1.23; const double bh_dimless_spin = 0.97; const double inner_edge_radius = 6.2; const double max_pressure_radius = 11.6; @@ -119,11 +119,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, @@ -140,21 +141,20 @@ 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 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 +172,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 +197,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 f58968ac3022..dc4fa4b198f0 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 685d58a929ec..853681f83fd6 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py @@ -8,6 +8,45 @@ 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 from SKS to 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_ji = np.diag([rho / r, rho / r, 1.0]) + return P_ji + + +def jacobian_matrix(x, a): + # Jacbian matrix for coordinate transformation from SKS to 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_ik = (-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_ik, 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 + + def sigma(r_sqrd, sin_theta_sqrd, a): return r_sqrd + (1.0 - sin_theta_sqrd) * a**2 @@ -37,9 +76,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 +86,40 @@ 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) +def kerr_schild_spatial_null_vec(x, m, 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, + ] ) + return l_vec + + +def sph_kerr_schild_spatial_null_form(x, m, a): + jac = jacobian_matrix(x, a) + l_vec = kerr_schild_spatial_null_vec(x, m, a) + return jac @ l_vec + + +def sph_kerr_schild_spatial_null_vec(x, m, a): + inv_jac = inverse_jacobian_matrix(x, a) + l_vec = kerr_schild_spatial_null_vec(x, m, a) + return inv_jac.T @ l_vec def kerr_schild_lapse(x, m, a): @@ -68,20 +130,22 @@ def kerr_schild_lapse(x, m, a): ) -def kerr_schild_shift(x, m, a): +def sph_kerr_schild_shift(x, m, a): null_vector_0 = -1.0 return ( -2.0 * kerr_schild_h(x, m, a) * null_vector_0 * kerr_schild_lapse(x, m, a) ** 2 - ) * kerr_schild_spatial_null_form(x, m, a) + ) * sph_kerr_schild_spatial_null_vec(x, m, a) -def kerr_schild_spatial_metric(x, m, a): +def sph_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) + null_form = sph_kerr_schild_spatial_null_form(x, m, a) + T_ij = jacobian_matrix(x, a) + result = T_ij @ T_ij.T + prefactor * np.outer(null_form, null_form) + return result def angular_momentum(m, a, rmax): @@ -142,7 +206,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,21 +315,29 @@ 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 ) - + kerr_schild_shift(x, bh_mass, bh_spin_a) + ) + transport_velocity_sks = ( + inverse_jacobian_matrix(x, bh_spin_a) @ transport_velocity_ks + ) + result += ( + transport_velocity_sks + + sph_kerr_schild_shift(x, bh_mass, bh_spin_a) ) / kerr_schild_lapse(x, bh_mass, bh_spin_a) - return result + return np.array(result) def lorentz_factor( @@ -279,7 +351,7 @@ def lorentz_factor( polytropic_exponent, ): bh_spin_a = bh_mass * bh_dimless_a - spatial_metric = kerr_schild_spatial_metric(x, bh_mass, bh_spin_a) + spatial_metric = sph_kerr_schild_spatial_metric(x, bh_mass, bh_spin_a) velocity = spatial_velocity( x, t, diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp index a0749b237b8f..34c3f2643211 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.23; 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>{}));