diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp index 8d7477fe0d25..ae6549b2e583 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp @@ -28,10 +28,10 @@ MagnetizedFmDisk::MagnetizedFmDisk( const double bh_mass, const double bh_dimless_spin, const double inner_edge_radius, const double max_pressure_radius, const double polytropic_constant, const double polytropic_exponent, - const double threshold_density, const double inverse_plasma_beta, - const size_t normalization_grid_res) + const double noise, const double threshold_density, + const double inverse_plasma_beta, const size_t normalization_grid_res) : fm_disk_(bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, - polytropic_constant, polytropic_exponent), + polytropic_constant, polytropic_exponent, noise), threshold_density_(threshold_density), inverse_plasma_beta_(inverse_plasma_beta), normalization_grid_res_(normalization_grid_res), diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp index 807dc3095651..9f50a358429c 100644 --- a/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp +++ b/src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp @@ -150,7 +150,7 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData, MagnetizedFmDisk( double bh_mass, double bh_dimless_spin, double inner_edge_radius, double max_pressure_radius, double polytropic_constant, - double polytropic_exponent, double threshold_density, + double polytropic_exponent, double noise, double threshold_density, double inverse_plasma_beta, size_t normalization_grid_res = BFieldNormGridRes::suggested_value()); diff --git a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp index b86b3ee06e78..2d96180de95e 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp +++ b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/EagerMath/DotProduct.hpp" @@ -31,13 +32,15 @@ FishboneMoncriefDisk::FishboneMoncriefDisk(const double bh_mass, const double inner_edge_radius, const double max_pressure_radius, const double polytropic_constant, - const double polytropic_exponent) + const double polytropic_exponent, + const double noise) : bh_mass_(bh_mass), bh_spin_a_(bh_mass * bh_dimless_spin), inner_edge_radius_(bh_mass * inner_edge_radius), max_pressure_radius_(bh_mass * max_pressure_radius), polytropic_constant_(polytropic_constant), polytropic_exponent_(polytropic_exponent), + noise_(noise), equation_of_state_{polytropic_constant_, polytropic_exponent_}, background_spacetime_{ bh_mass_, {{0.0, 0.0, bh_dimless_spin}}, {{0.0, 0.0, 0.0}}} { @@ -52,6 +55,13 @@ FishboneMoncriefDisk::FishboneMoncriefDisk(const double bh_mass, (square(bh_spin_a_) - 2.0 * a_sqrt_m * sqrt_rmax + rmax_squared) / (2.0 * a_sqrt_m * rmax_sqrt_rmax + (rmax - 3.0 * bh_mass_) * rmax_squared); + + // compute rho_max_ here: + const double potential_at_rmax = potential(rmax_squared, 1.0); + const double potential_at_rin = potential(square(inner_edge_radius_), 1.0); + const double h_at_rmax = exp(potential_at_rin - potential_at_rmax); + rho_max_ = get(equation_of_state_.rest_mass_density_from_enthalpy( + Scalar{h_at_rmax})); } std::unique_ptr @@ -68,6 +78,8 @@ void FishboneMoncriefDisk::pup(PUP::er& p) { p | polytropic_constant_; p | polytropic_exponent_; p | angular_momentum_; + p | rho_max_; + p | noise_; p | equation_of_state_; p | background_spacetime_; } @@ -140,6 +152,7 @@ FishboneMoncriefDisk::variables( 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) = + (1 / rho_max_) * get(equation_of_state_.rest_mass_density_from_enthalpy( Scalar{get_element(get(specific_enthalpy), s)})); }); @@ -183,14 +196,18 @@ FishboneMoncriefDisk::variables( const tnsr::I& x, tmpl::list> /*meta*/, gsl::not_null*> vars) const { + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution<> dis(-noise_, noise_); const auto rest_mass_density = get>( variables(x, tmpl::list>{}, vars)); auto pressure = make_with_value>(x, 0.0); - variables_impl(vars, [&pressure, &rest_mass_density, this]( + variables_impl(vars, [&pressure, &rest_mass_density, &gen, &dis, this]( const size_t s, const double /*potential_at_s*/) { get_element(get(pressure), s) = + (1 / rho_max_) * (1 + dis(gen)) * get(equation_of_state_.pressure_from_density( - Scalar{get_element(get(rest_mass_density), s)})); + Scalar{rho_max_ * get_element(get(rest_mass_density), s)})); }); return {std::move(pressure)}; } @@ -208,7 +225,7 @@ FishboneMoncriefDisk::variables( const size_t s, const double /*potential_at_s*/) { get_element(get(specific_internal_energy), s) = get(equation_of_state_.specific_internal_energy_from_density( - Scalar{get_element(get(rest_mass_density), s)})); + Scalar{rho_max_ * get_element(get(rest_mass_density), s)})); }); return {std::move(specific_internal_energy)}; } @@ -358,7 +375,8 @@ bool operator==(const FishboneMoncriefDisk& lhs, lhs.inner_edge_radius_ == rhs.inner_edge_radius_ and lhs.max_pressure_radius_ == rhs.max_pressure_radius_ and lhs.polytropic_constant_ == rhs.polytropic_constant_ and - lhs.polytropic_exponent_ == rhs.polytropic_exponent_; + lhs.polytropic_exponent_ == rhs.polytropic_exponent_ and + lhs.noise_ == rhs.noise_ and lhs.rho_max_ == rhs.rho_max_; } bool operator!=(const FishboneMoncriefDisk& lhs, diff --git a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp index 21bda856c652..6372e495dd5b 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp @@ -126,6 +126,9 @@ namespace Solutions { * p = K\rho^\gamma. * \f} * + * Following \cite Porth2016rfi, the rest mass density is normalized + * to be 1. at maximum. + * * Once all the variables are known in Boyer-Lindquist (or Kerr) coordinates, it * is straightforward to write them in Cartesian Kerr-Schild coordinates. The * coordinate transformation in gr::KerrSchildCoords helps read the Jacobian @@ -213,10 +216,20 @@ class FishboneMoncriefDisk "The polytropic exponent of the fluid."}; static type lower_bound() { return 1.; } }; + /// The magnitude of noise added to pressure/energy of the Disk + /// to drive MRI. + struct Noise { + using type = double; + static constexpr Options::String help = { + "The magnitude of the white noise perturbation added to " + "pressure to excite MRI in the disk."}; + static type lower_bound() { return 0.; } + static type upper_bound() { return 1.; } + }; using options = tmpl::list; + PolytropicConstant, PolytropicExponent, Noise>; static constexpr Options::String help = { "Fluid disk orbiting a Kerr black hole."}; @@ -230,7 +243,8 @@ class FishboneMoncriefDisk FishboneMoncriefDisk(double bh_mass, double bh_dimless_spin, double inner_edge_radius, double max_pressure_radius, - double polytropic_constant, double polytropic_exponent); + double polytropic_constant, double polytropic_exponent, + double noise); auto get_clone() const -> std::unique_ptr override; @@ -410,6 +424,8 @@ class FishboneMoncriefDisk double polytropic_constant_ = std::numeric_limits::signaling_NaN(); double polytropic_exponent_ = std::numeric_limits::signaling_NaN(); double angular_momentum_ = std::numeric_limits::signaling_NaN(); + double rho_max_ = std::numeric_limits::signaling_NaN(); + double noise_ = std::numeric_limits::signaling_NaN(); EquationsOfState::PolytropicFluid equation_of_state_{}; gr::Solutions::SphericalKerrSchild background_spacetime_{}; }; diff --git a/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml b/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml index 5fb19113fcbe..6be821c219c3 100644 --- a/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml +++ b/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml @@ -55,6 +55,7 @@ InitialData: &initial_data MaxPressureRadius: 12.0 PolytropicConstant: 0.001 PolytropicExponent: 1.3333333333333333333333 + Noise: 0.0 EquationOfState: FromInitialData diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py index eb69ec600367..90f22865d03b 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py @@ -14,6 +14,7 @@ def rest_mass_density( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -27,6 +28,7 @@ def rest_mass_density( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) @@ -38,6 +40,7 @@ def spatial_velocity( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -51,6 +54,7 @@ def spatial_velocity( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) @@ -62,6 +66,7 @@ def specific_internal_energy( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -75,6 +80,7 @@ def specific_internal_energy( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) @@ -86,6 +92,7 @@ def pressure( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -99,6 +106,7 @@ def pressure( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) @@ -111,6 +119,7 @@ def magnetic_potential( rmax, polytropic_constant, polytropic_exponent, + noise, threshold_rest_mass_density, ): l = fm_disk.angular_momentum(m, a, rmax) @@ -135,6 +144,7 @@ def magnetic_field( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -152,6 +162,7 @@ def magnetic_field( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) x_max = np.array([r_max, 0.0, 0.0]) threshold_rho = threshold_density * ( @@ -164,6 +175,7 @@ def magnetic_field( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) ) x_ks = [0.0, 0.0, 0.0] @@ -191,6 +203,7 @@ def magnetic_field( r_max, polytropic_constant, polytropic_exponent, + noise, threshold_rho, ) - magnetic_potential( @@ -202,6 +215,7 @@ def magnetic_field( r_max, polytropic_constant, polytropic_exponent, + noise, threshold_rho, ) ) @@ -220,6 +234,7 @@ def magnetic_field( r_max, polytropic_constant, polytropic_exponent, + noise, threshold_rho, ) - magnetic_potential( @@ -231,6 +246,7 @@ def magnetic_field( r_max, polytropic_constant, polytropic_exponent, + noise, threshold_rho, ) ) * prefactor @@ -242,7 +258,7 @@ def magnetic_field( # 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 + normalization = 8.49943510498341 result = normalization * ks_coords.cartesian_from_spherical_ks( result, x_ks, bh_mass, bh_dimless_spin ) @@ -259,6 +275,7 @@ def divergence_cleaning_field( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -273,6 +290,7 @@ def lorentz_factor( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -286,6 +304,7 @@ def lorentz_factor( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) @@ -297,6 +316,7 @@ def specific_enthalpy( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, threshold_density, plasma_beta, ): @@ -310,4 +330,5 @@ def specific_enthalpy( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp index 8e41be80d254..faa6a0defdde 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedFmDisk.cpp @@ -82,6 +82,7 @@ void test_create_from_options() { " MaxPressureRadius: 14.2\n" " PolytropicConstant: 0.065\n" " PolytropicExponent: 1.654\n" + " Noise: 0.0\n" " ThresholdDensity: 0.42\n" " InversePlasmaBeta: 85.0\n" " BFieldNormGridRes: 4") @@ -92,20 +93,20 @@ 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, 4)); + 1.3, 0.345, 6.123, 14.2, 0.065, 1.654, 0.0, 0.42, 85.0, 4)); } void test_move() { grmhd::AnalyticData::MagnetizedFmDisk disk(3.51, 0.87, 7.43, 15.3, 42.67, - 1.87, 0.13, 0.015, 4); - grmhd::AnalyticData::MagnetizedFmDisk disk_copy(3.51, 0.87, 7.43, 15.3, 42.67, - 1.87, 0.13, 0.015, 4); + 1.87, 0.0, 0.13, 0.015, 4); + const grmhd::AnalyticData::MagnetizedFmDisk disk_copy( + 3.51, 0.87, 7.43, 15.3, 42.67, 1.87, 0.0, 0.13, 0.015, 4); test_move_semantics(std::move(disk), disk_copy); // NOLINT } void test_serialize() { - grmhd::AnalyticData::MagnetizedFmDisk disk(3.51, 0.87, 7.43, 15.3, 42.67, - 1.87, 0.13, 0.015, 4); + const grmhd::AnalyticData::MagnetizedFmDisk disk( + 3.51, 0.87, 7.43, 15.3, 42.67, 1.87, 0.0, 0.13, 0.015, 4); test_serialization(disk); } @@ -117,17 +118,18 @@ void test_variables(const DataType& used_for_size) { const double max_pressure_radius = 11.6; const double polytropic_constant = 0.034; const double polytropic_exponent = 1.65; + const double noise = 0.0; 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, b_field_normalization); + const MagnetizedFmDiskProxy disk( + bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, + polytropic_constant, polytropic_exponent, noise, threshold_density, + 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, + polytropic_constant, polytropic_exponent, noise, threshold_density, inverse_plasma_beta); pypp::check_with_random_values<1>( @@ -162,9 +164,10 @@ void test_variables(const DataType& used_for_size) { // Check that when InversePlasmaBeta = 0, magnetic field vanishes and // we recover FishboneMoncriefDisk - MagnetizedFmDiskProxy another_disk( + const MagnetizedFmDiskProxy another_disk( bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, - polytropic_constant, polytropic_exponent, threshold_density, 0.0, 4); + polytropic_constant, polytropic_exponent, noise, threshold_density, 0.0, + 4); pypp::check_with_random_values<1>( &MagnetizedFmDiskProxy::hydro_variables, another_disk, @@ -201,36 +204,36 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.MagFmDisk", CHECK_THROWS_WITH( []() { grmhd::AnalyticData::MagnetizedFmDisk disk( - 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, -0.2, 0.023, 4); + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 0.0, -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); + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 0.0, 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); + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 0.0, 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); + grmhd::AnalyticData::MagnetizedFmDisk disk( + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 0.0, 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); + grmhd::AnalyticData::MagnetizedFmDisk disk( + 0.7, 0.61, 5.4, 9.182, 11.123, 1.44, 0.0, 0.2, 0.023, 4); }(), Catch::Matchers::ContainsSubstring("Max b squared is zero.")); #endif @@ -242,6 +245,7 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.MagFmDisk", "MaxPressureRadius: 7.6\n" "PolytropicConstant: 2.42\n" "PolytropicExponent: 1.33\n" + "Noise: 0.0\n" "ThresholdDensity: -0.01\n" "InversePlasmaBeta: 0.016\n" "BFieldNormGridRes: 4"), @@ -255,6 +259,7 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.MagFmDisk", "MaxPressureRadius: 8.2\n" "PolytropicConstant: 41.1\n" "PolytropicExponent: 1.8\n" + "Noise: 0.0\n" "ThresholdDensity: 4.1\n" "InversePlasmaBeta: 0.03\n" "BFieldNormGridRes: 4"), @@ -268,6 +273,7 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.MagFmDisk", "MaxPressureRadius: 7.8\n" "PolytropicConstant: 13.5\n" "PolytropicExponent: 1.54\n" + "Noise : 0.0\n" "ThresholdDensity: 0.22\n" "InversePlasmaBeta: -0.03\n" "BFieldNormGridRes: 4"), @@ -281,6 +287,7 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticData.GrMhd.MagFmDisk", "MaxPressureRadius: 7.8\n" "PolytropicConstant: 13.5\n" "PolytropicExponent: 1.54\n" + "Noise: 0.0\n" "ThresholdDensity: 0.22\n" "InversePlasmaBeta: 0.03\n" "BFieldNormGridRes: 2"), diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp index 0d2bbadabd3a..e63ba3d28eec 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_PolarMagnetizedFmDisk.cpp @@ -47,6 +47,7 @@ void test_create_from_options() { " MaxPressureRadius: 14.2\n" " PolytropicConstant: 0.065\n" " PolytropicExponent: 1.654\n" + " Noise: 0.0\n" " ThresholdDensity: 0.42\n" " InversePlasmaBeta: 85.0\n" " BFieldNormGridRes: 4\n" @@ -63,19 +64,19 @@ 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, 4), + grmhd::AnalyticData::MagnetizedFmDisk( + 1.3, 0.345, 6.123, 14.2, 0.065, 1.654, 0.0, 0.42, 85.0, 4), grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7, 0.0))); } 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, 4), + 1.654, 0.0, 0.42, 85.0, 4), grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7, 0.0)); grmhd::AnalyticData::PolarMagnetizedFmDisk disk_copy( grmhd::AnalyticData::MagnetizedFmDisk(1.3, 0.345, 6.123, 14.2, 0.065, - 1.654, 0.42, 85.0, 4), + 1.654, 0.0, 0.42, 85.0, 4), grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7, 0.0)); test_move_semantics(std::move(disk), disk_copy); } @@ -83,7 +84,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, 4), + 1.654, 0.0, 0.42, 85.0, 4), grmhd::AnalyticData::SphericalTorus(2.5, 3.6, 0.9, 0.7, 0.0)); test_serialization(disk); } @@ -100,11 +101,12 @@ struct Wrapper { template void test_variables(const DataType& used_for_size) { const double bh_mass = 1.12; - const double bh_dimless_spin = 0.97; + const double bh_dimless_spin = 0.9375; const double inner_edge_radius = 6.2; const double max_pressure_radius = 11.6; const double polytropic_constant = 0.034; const double polytropic_exponent = 1.65; + const double noise = 0.0; const double threshold_density = 0.14; const double inverse_plasma_beta = 0.023; const size_t b_field_normalization = 51; @@ -112,7 +114,7 @@ void test_variables(const DataType& used_for_size) { 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, + polytropic_constant, polytropic_exponent, noise, threshold_density, inverse_plasma_beta, b_field_normalization), grmhd::AnalyticData::SphericalTorus(3.0, 20.0, 1.0, 0.3, 0.0)); diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/GhRelativisticEuler/Test_InstantiateWrappedGr.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/GhRelativisticEuler/Test_InstantiateWrappedGr.cpp index b904ff329405..235ed8a8f284 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/GhRelativisticEuler/Test_InstantiateWrappedGr.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/GhRelativisticEuler/Test_InstantiateWrappedGr.cpp @@ -22,10 +22,10 @@ SPECTRE_TEST_CASE( const tnsr::I x_1d{DataVector{3.0, 4.0}}; check_wrapped_gr_solution_consistency( gh::Solutions::WrappedGr< - RelativisticEuler::Solutions::FishboneMoncriefDisk>{1.0, 0.23, 6.0, - 12.0, 0.001, 1.4}, + RelativisticEuler::Solutions::FishboneMoncriefDisk>{ + 1.0, 0.23, 6.0, 12.0, 0.001, 1.4, 0.0}, RelativisticEuler::Solutions::FishboneMoncriefDisk{1.0, 0.23, 6.0, 12.0, - 0.001, 1.4}, + 0.001, 1.4, 0.0}, x_3d, t); check_wrapped_gr_solution_consistency( diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py index 853681f83fd6..bbcd60758c51 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py @@ -4,6 +4,31 @@ import numpy as np +def rho_max( + bh_mass, + bh_dimless_a, + dimless_r_in, + dimless_r_max, + polytropic_constant, + polytropic_exponent, +): + r_in = bh_mass * dimless_r_in + r_max = bh_mass * dimless_r_max + 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) + Wmax = potential(l, r_max * r_max, 1.0, bh_mass, bh_spin_a) + h_max = np.exp(Win - Wmax) + return np.power( + ( + (polytropic_exponent - 1.0) + * (h_max - 1.0) + / (polytropic_constant * polytropic_exponent) + ), + 1.0 / (polytropic_exponent - 1.0), + ) + + def delta(r_sqrd, m, a): return r_sqrd - 2.0 * m * np.sqrt(r_sqrd) + a**2 @@ -201,6 +226,7 @@ def specific_enthalpy( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ): r_in = bh_mass * dimless_r_in bh_spin_a = bh_mass * bh_dimless_a @@ -226,8 +252,17 @@ def rest_mass_density( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ): - return np.power( + rho_max_ = rho_max( + bh_mass, + bh_dimless_a, + dimless_r_in, + dimless_r_max, + polytropic_constant, + polytropic_exponent, + ) + return (1 / rho_max_) * np.power( (polytropic_exponent - 1.0) * ( specific_enthalpy( @@ -239,6 +274,7 @@ def rest_mass_density( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) - 1.0 ) @@ -256,11 +292,21 @@ def specific_internal_energy( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ): + rho_max_ = rho_max( + bh_mass, + bh_dimless_a, + dimless_r_in, + dimless_r_max, + polytropic_constant, + polytropic_exponent, + ) return ( polytropic_constant * np.power( - rest_mass_density( + rho_max_ + * rest_mass_density( x, t, bh_mass, @@ -269,6 +315,7 @@ def specific_internal_energy( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ), polytropic_exponent - 1.0, ) @@ -285,20 +332,35 @@ def pressure( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ): - return polytropic_constant * np.power( - rest_mass_density( - x, - t, - bh_mass, - bh_dimless_a, - dimless_r_in, - dimless_r_max, - polytropic_constant, - polytropic_exponent, - ), + rho_max_ = rho_max( + bh_mass, + bh_dimless_a, + dimless_r_in, + dimless_r_max, + polytropic_constant, polytropic_exponent, ) + return ( + (1 / rho_max_) + * polytropic_constant + * np.power( + rho_max_ + * rest_mass_density( + x, + t, + bh_mass, + bh_dimless_a, + dimless_r_in, + dimless_r_max, + polytropic_constant, + polytropic_exponent, + noise, + ), + polytropic_exponent, + ) + ) def spatial_velocity( @@ -310,6 +372,7 @@ def spatial_velocity( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ): r_in = bh_mass * dimless_r_in bh_spin_a = bh_mass * bh_dimless_a @@ -349,6 +412,7 @@ def lorentz_factor( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ): bh_spin_a = bh_mass * bh_dimless_a spatial_metric = sph_kerr_schild_spatial_metric(x, bh_mass, bh_spin_a) @@ -361,6 +425,7 @@ def lorentz_factor( dimless_r_max, polytropic_constant, polytropic_exponent, + noise, ) return 1.0 / np.sqrt( 1.0 - np.einsum("i,j,ij", velocity, velocity, spatial_metric) diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp index 15cd6d7ac87f..332fc9b04195 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_FishboneMoncriefDisk.cpp @@ -78,22 +78,23 @@ void test_create_from_options() { "InnerEdgeRadius: 6.0\n" "MaxPressureRadius: 12.0\n" "PolytropicConstant: 0.001\n" - "PolytropicExponent: 1.4"); + "PolytropicExponent: 1.4\n" + "Noise: 0.0"); CHECK(disk == RelativisticEuler::Solutions::FishboneMoncriefDisk( - 1.0, 0.23, 6.0, 12.0, 0.001, 1.4)); + 1.0, 0.23, 6.0, 12.0, 0.001, 1.4, 0.0)); } void test_move() { RelativisticEuler::Solutions::FishboneMoncriefDisk disk(3.45, 0.23, 4.8, 8.6, - 0.02, 1.5); - RelativisticEuler::Solutions::FishboneMoncriefDisk disk_copy(3.45, 0.23, 4.8, - 8.6, 0.02, 1.5); + 0.02, 1.5, 0.0); + const RelativisticEuler::Solutions::FishboneMoncriefDisk disk_copy( + 3.45, 0.23, 4.8, 8.6, 0.02, 1.5, 0.0); test_move_semantics(std::move(disk), disk_copy); // NOLINT } void test_serialize() { RelativisticEuler::Solutions::FishboneMoncriefDisk disk(4.21, 0.65, 6.0, 12.0, - 0.001, 1.4); + 0.001, 1.4, 0.0); test_serialization(disk); } @@ -105,13 +106,14 @@ void test_variables(const DataType& used_for_size) { const double max_pressure_radius = 12.0; const double polytropic_constant = 0.001; const double polytropic_exponent = 4.0 / 3.0; + const double noise = 0.0; - FishboneMoncriefDiskProxy disk(bh_mass, bh_dimless_spin, inner_edge_radius, - max_pressure_radius, polytropic_constant, - polytropic_exponent); + const FishboneMoncriefDiskProxy disk( + bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, + polytropic_constant, polytropic_exponent, noise); const auto member_variables = std::make_tuple( bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius, - polytropic_constant, polytropic_exponent); + polytropic_constant, polytropic_exponent, noise); pypp::check_with_random_values<1>( &FishboneMoncriefDiskProxy::hydro_variables, disk, @@ -163,8 +165,8 @@ template void test_sin_theta_squared(const DataType& used_for_size) { // Numbers below reproduce the initial data the bug was spotted with, // along with the points where the FPEs were found. - FishboneMoncriefDiskProxy disk(1.0, 0.9375, 6.0, 12.0, 0.001, - 1.3333333333333333333333); + const FishboneMoncriefDiskProxy disk(1.0, 0.9375, 6.0, 12.0, 0.001, + 1.3333333333333333333333, 0.0); using variables_tags = FishboneMoncriefDiskProxy::grmhd_variables_tags; auto coords = make_with_value>(used_for_size, 0.0); @@ -192,7 +194,8 @@ void test_solution() { " InnerEdgeRadius: 6.0\n" " MaxPressureRadius: 12.0\n" " PolytropicConstant: 0.001\n" - " PolytropicExponent: 1.33333333333333333\n") + " PolytropicExponent: 1.33333333333333333\n" + " Noise: 0.0\n") ->get_clone(); const auto deserialized_option_solution = serialize_and_deserialize(option_solution); @@ -236,7 +239,8 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.RelEuler.FMDisk", "InnerEdgeRadius: 4.3\n" "MaxPressureRadius: 6.7\n" "PolytropicConstant: 0.12\n" - "PolytropicExponent: 1.5"), + "PolytropicExponent: 1.5\n" + "Noise : 0.0"), Catch::Matchers::ContainsSubstring( "Value -1.5 is below the lower bound of 0")); CHECK_THROWS_WITH(TestHelpers::test_creation< @@ -246,7 +250,8 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.RelEuler.FMDisk", "InnerEdgeRadius: 5.76\n" "MaxPressureRadius: 13.2\n" "PolytropicConstant: 0.002\n" - "PolytropicExponent: 1.34"), + "PolytropicExponent: 1.34\n" + "Noise: 0.0"), Catch::Matchers::ContainsSubstring( "Value -0.24 is below the lower bound of 0")); CHECK_THROWS_WITH(TestHelpers::test_creation< @@ -256,7 +261,8 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.RelEuler.FMDisk", "InnerEdgeRadius: 5.76\n" "MaxPressureRadius: 13.2\n" "PolytropicConstant: 0.002\n" - "PolytropicExponent: 1.34"), + "PolytropicExponent: 1.34\n" + "Noise: 0.0"), Catch::Matchers::ContainsSubstring( "Value 1.24 is above the upper bound of 1")); CHECK_THROWS_WITH(TestHelpers::test_creation< @@ -266,7 +272,8 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.RelEuler.FMDisk", "InnerEdgeRadius: 4.3\n" "MaxPressureRadius: 6.7\n" "PolytropicConstant: -0.12\n" - "PolytropicExponent: 1.5"), + "PolytropicExponent: 1.5\n" + "Noise: 0.0"), Catch::Matchers::ContainsSubstring( "Value -0.12 is below the lower bound of 0")); CHECK_THROWS_WITH(TestHelpers::test_creation< @@ -276,7 +283,30 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.RelEuler.FMDisk", "InnerEdgeRadius: 4.3\n" "MaxPressureRadius: 6.7\n" "PolytropicConstant: 0.123\n" - "PolytropicExponent: 0.25"), + "PolytropicExponent: 0.25\n" + "Noise: 0.0"), Catch::Matchers::ContainsSubstring( "Value 0.25 is below the lower bound of 1")); + CHECK_THROWS_WITH(TestHelpers::test_creation< + RelativisticEuler::Solutions::FishboneMoncriefDisk>( + "BhMass: 1.5\n" + "BhDimlessSpin: 0.3\n" + "InnerEdgeRadius: 4.3\n" + "MaxPressureRadius: 6.7\n" + "PolytropicConstant: 0.123\n" + "PolytropicExponent: 1.5\n" + "Noise: -1.25"), + Catch::Matchers::ContainsSubstring( + "Value -1.25 is below the lower bound of 0")); + CHECK_THROWS_WITH(TestHelpers::test_creation< + RelativisticEuler::Solutions::FishboneMoncriefDisk>( + "BhMass: 1.5\n" + "BhDimlessSpin: 0.3\n" + "InnerEdgeRadius: 4.3\n" + "MaxPressureRadius: 6.7\n" + "PolytropicConstant: 0.123\n" + "PolytropicExponent: 1.5\n" + "Noise: 1.5"), + Catch::Matchers::ContainsSubstring( + "Value 1.5 is above the upper bound of 1")); }