From c946b61e208c583fbb06c04a07e9a774cc072313 Mon Sep 17 00:00:00 2001 From: Nils Deppe Date: Tue, 10 Oct 2023 15:12:35 -0700 Subject: [PATCH 1/5] Add Enable option to LimitLorentzFactor --- .../VariableFixing/LimitLorentzFactor.cpp | 14 ++++-- .../VariableFixing/LimitLorentzFactor.hpp | 12 ++++- .../GhValenciaDivClean/BinaryNeutronStar.yaml | 1 + .../GhValenciaDivClean/GhMhdBondiMichel.yaml | 1 + .../GhValenciaDivClean/GhMhdTovStar.yaml | 1 + .../Test_LimitLorentzFactor.cpp | 50 ++++++++++++------- 6 files changed, 55 insertions(+), 24 deletions(-) diff --git a/src/Evolution/VariableFixing/LimitLorentzFactor.cpp b/src/Evolution/VariableFixing/LimitLorentzFactor.cpp index 7bb818a45349..3ec5e8596c3e 100644 --- a/src/Evolution/VariableFixing/LimitLorentzFactor.cpp +++ b/src/Evolution/VariableFixing/LimitLorentzFactor.cpp @@ -16,13 +16,16 @@ namespace VariableFixing { LimitLorentzFactor::LimitLorentzFactor(const double max_density_cutoff, - const double lorentz_factor_cap) + const double lorentz_factor_cap, + const bool enable) : max_density_cuttoff_(max_density_cutoff), - lorentz_factor_cap_(lorentz_factor_cap) {} + lorentz_factor_cap_(lorentz_factor_cap), + enable_(enable) {} void LimitLorentzFactor::pup(PUP::er& p) { p | max_density_cuttoff_; p | lorentz_factor_cap_; + p | enable_; } void LimitLorentzFactor::operator()( @@ -30,6 +33,10 @@ void LimitLorentzFactor::operator()( const gsl::not_null*> spatial_velocity, const Scalar& rest_mass_density) const { + if (not enable_) { + return; + } + constexpr size_t dim = 3; const size_t number_of_grid_points = get(rest_mass_density).size(); for (size_t s = 0; s < number_of_grid_points; ++s) { @@ -49,7 +56,8 @@ void LimitLorentzFactor::operator()( bool operator==(const LimitLorentzFactor& lhs, const LimitLorentzFactor& rhs) { return lhs.max_density_cuttoff_ == rhs.max_density_cuttoff_ and - lhs.lorentz_factor_cap_ == rhs.lorentz_factor_cap_; + lhs.lorentz_factor_cap_ == rhs.lorentz_factor_cap_ and + lhs.enable_ == rhs.enable_; } bool operator!=(const LimitLorentzFactor& lhs, const LimitLorentzFactor& rhs) { diff --git a/src/Evolution/VariableFixing/LimitLorentzFactor.hpp b/src/Evolution/VariableFixing/LimitLorentzFactor.hpp index 272885233f86..47ee95d9e0f5 100644 --- a/src/Evolution/VariableFixing/LimitLorentzFactor.hpp +++ b/src/Evolution/VariableFixing/LimitLorentzFactor.hpp @@ -51,15 +51,22 @@ class LimitLorentzFactor { static type lower_bound() { return 1.0; } static constexpr Options::String help = {"Largest Lorentz factor allowed."}; }; + /// Whether or not the limiting is enabled + struct Enable { + using type = bool; + static constexpr Options::String help = { + "If true then the limiting is applied."}; + }; - using options = tmpl::list; + using options = tmpl::list; static constexpr Options::String help = { "Limit the maximum Lorentz factor to LorentzFactorCap in regions where " "the\n" "density is below MaxDensityCutoff. The Lorentz factor is set to\n" "LorentzFactorCap and the spatial velocity is adjusted accordingly."}; - LimitLorentzFactor(double max_density_cutoff, double lorentz_factor_cap); + LimitLorentzFactor(double max_density_cutoff, double lorentz_factor_cap, + bool enable); LimitLorentzFactor() = default; LimitLorentzFactor(const LimitLorentzFactor& /*rhs*/) = default; @@ -87,6 +94,7 @@ class LimitLorentzFactor { double max_density_cuttoff_; double lorentz_factor_cap_; + bool enable_; }; bool operator!=(const LimitLorentzFactor& lhs, const LimitLorentzFactor& rhs); diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml index 2f1f4d03fdad..d17a3c5278c3 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml @@ -79,6 +79,7 @@ VariableFixing: TransitionDensityCutoff: 9.0e-15 MaxVelocityMagnitude: 1.0e-4 LimitLorentzFactor: + Enable: false # Only needed when not using Kastaun for recovery MaxDensityCutoff: 1.0e-08 LorentzFactorCap: 10.0 diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml index c2f8cb19ccfd..aba7d4c3ef96 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml @@ -80,6 +80,7 @@ VariableFixing: TransitionDensityCutoff: 1.0e-11 MaxVelocityMagnitude: 1.0e-4 LimitLorentzFactor: + Enable: false # Only needed when not using Kastaun for recovery MaxDensityCutoff: 1.0e-12 LorentzFactorCap: 1.0 diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml index 03f6c81b5ccc..cc35359835f1 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml @@ -73,6 +73,7 @@ VariableFixing: TransitionDensityCutoff: 1.0e-11 MaxVelocityMagnitude: 1.0e-4 LimitLorentzFactor: + Enable: false # Only needed when not using Kastaun for recovery MaxDensityCutoff: 1.0e-12 LorentzFactorCap: 1.0 diff --git a/tests/Unit/Evolution/VariableFixing/Test_LimitLorentzFactor.cpp b/tests/Unit/Evolution/VariableFixing/Test_LimitLorentzFactor.cpp index 6d2a6a3c0728..ecc96855e191 100644 --- a/tests/Unit/Evolution/VariableFixing/Test_LimitLorentzFactor.cpp +++ b/tests/Unit/Evolution/VariableFixing/Test_LimitLorentzFactor.cpp @@ -21,7 +21,8 @@ namespace VariableFixing { namespace { -void test_variable_fixer(const LimitLorentzFactor& variable_fixer) { +void test_variable_fixer(const LimitLorentzFactor& variable_fixer, + const bool enable) { const Scalar density{DataVector{1.0, 2.0, 1.0e-5, 1.0e-6}}; auto spatial_metric = make_with_value>(density, 1.0); @@ -39,12 +40,13 @@ void test_variable_fixer(const LimitLorentzFactor& variable_fixer) { auto lorentz_factor = hydro::lorentz_factor( dot_product(spatial_velocity, spatial_velocity, spatial_metric)); - const Scalar expected_lorentz_factor{ - DataVector{get(lorentz_factor)[0], get(lorentz_factor)[1], 50.0, - get(lorentz_factor)[3]}}; + const Scalar expected_lorentz_factor{DataVector{ + get(lorentz_factor)[0], get(lorentz_factor)[1], + enable ? 50.0 : get(lorentz_factor)[2], get(lorentz_factor)[3]}}; const double rescale_velocity_factor = - sqrt((1.0 - 1.0 / square(get(expected_lorentz_factor)[2])) / - (1.0 - 1.0 / square(get(lorentz_factor)[2]))); + enable ? sqrt((1.0 - 1.0 / square(get(expected_lorentz_factor)[2])) / + (1.0 - 1.0 / square(get(lorentz_factor)[2]))) + : 1.0; const tnsr::I expected_spatial_velocity{ {{DataVector{0.4999, 0.4999, 0.4999 * rescale_velocity_factor, 0.4999}, DataVector{1.3567, 1.3566, 1.3567 * rescale_velocity_factor, 1.3566}, @@ -66,25 +68,35 @@ void test_variable_fixer(const LimitLorentzFactor& variable_fixer) { SPECTRE_TEST_CASE("Unit.Evolution.VariableFixing.LimitLorentzFactor", "[VariableFixing][Unit]") { SECTION("operator== and operator!=") { - CHECK(LimitLorentzFactor{1.0, 8.0} == LimitLorentzFactor{1.0, 8.0}); - CHECK_FALSE(LimitLorentzFactor{2., 8.0} == LimitLorentzFactor{1.0, 8.0}); - CHECK_FALSE(LimitLorentzFactor{1.0, 9.0} == LimitLorentzFactor{1.0, 8.0}); - CHECK_FALSE(LimitLorentzFactor{1.0, 8.0} != LimitLorentzFactor{1.0, 8.0}); - CHECK(LimitLorentzFactor{2.0, 8.0} != LimitLorentzFactor{1.0, 8.0}); - CHECK(LimitLorentzFactor{1.0, 9.0} != LimitLorentzFactor{1.0, 8.0}); + CHECK(LimitLorentzFactor{1.0, 8.0, true} == + LimitLorentzFactor{1.0, 8.0, true}); + CHECK_FALSE(LimitLorentzFactor{2., 8.0, true} == + LimitLorentzFactor{1.0, 8.0, true}); + CHECK_FALSE(LimitLorentzFactor{1.0, 9.0, true} == + LimitLorentzFactor{1.0, 8.0, true}); + CHECK_FALSE(LimitLorentzFactor{1.0, 8.0, true} != + LimitLorentzFactor{1.0, 8.0, true}); + CHECK_FALSE(LimitLorentzFactor{1.0, 8.0, true} == + LimitLorentzFactor{1.0, 8.0, false}); + CHECK(LimitLorentzFactor{2.0, 8.0, true} != + LimitLorentzFactor{1.0, 8.0, true}); + CHECK(LimitLorentzFactor{1.0, 9.0, true} != + LimitLorentzFactor{1.0, 8.0, true}); } SECTION("variable fixing") { - LimitLorentzFactor variable_fixer{1.0e-4, 50.0}; - test_variable_fixer(variable_fixer); - test_serialization(variable_fixer); - test_variable_fixer(serialize_and_deserialize(variable_fixer)); - + for (const bool enable : {true, false}) { + CAPTURE(enable); + LimitLorentzFactor variable_fixer{1.0e-4, 50.0, enable}; + test_serialization(variable_fixer); + test_variable_fixer(serialize_and_deserialize(variable_fixer), enable); + } const auto fixer_from_options = TestHelpers::test_creation( "MaxDensityCutoff: 1.0e-4\n" - "LorentzFactorCap: 50.0\n"); - test_variable_fixer(fixer_from_options); + "LorentzFactorCap: 50.0\n" + "Enable: true\n"); + test_variable_fixer(fixer_from_options, true); } } } // namespace From 98e16d2a7446bd7f876805590ef0540d6ba65bac Mon Sep 17 00:00:00 2001 From: Nils Deppe Date: Tue, 10 Oct 2023 15:37:46 -0700 Subject: [PATCH 2/5] Add Enable option to FixConservatives --- .../ValenciaDivClean/FixConservatives.cpp | 12 ++++-- .../ValenciaDivClean/FixConservatives.hpp | 30 ++++++++----- .../GhValenciaDivClean/BinaryNeutronStar.yaml | 1 + .../GhValenciaDivClean/GhMhdBondiMichel.yaml | 1 + .../GhValenciaDivClean/GhMhdTovStar.yaml | 1 + .../GrMhd/ValenciaDivClean/BlastWave.yaml | 1 + .../FishboneMoncriefDisk.yaml | 1 + .../Test_FixConservativesAndComputePrims.cpp | 2 +- .../Test_FixConservativesAndComputePrims.cpp | 2 +- .../Test_FixConservatives.cpp | 42 ++++++++++++------- 10 files changed, 62 insertions(+), 31 deletions(-) diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp index 9f241f821a01..1d23f71360ce 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.cpp @@ -104,7 +104,7 @@ FixConservatives::FixConservatives( const double safety_factor_for_magnetic_field, const double safety_factor_for_momentum_density, const double safety_factor_for_momentum_density_cutoff_d, - const double safety_factor_for_momentum_density_slope, + const double safety_factor_for_momentum_density_slope, const bool enable, const Options::Context& context) : minimum_rest_mass_density_times_lorentz_factor_( minimum_rest_mass_density_times_lorentz_factor), @@ -119,7 +119,8 @@ FixConservatives::FixConservatives( safety_factor_for_momentum_density_cutoff_d_( safety_factor_for_momentum_density_cutoff_d), safety_factor_for_momentum_density_slope_( - safety_factor_for_momentum_density_slope) { + safety_factor_for_momentum_density_slope), + enable_(enable) { if (minimum_rest_mass_density_times_lorentz_factor_ > rest_mass_density_times_lorentz_factor_cutoff_) { PARSE_ERROR(context, @@ -161,6 +162,7 @@ void FixConservatives::pup(PUP::er& p) { p | one_minus_safety_factor_for_momentum_density_; p | safety_factor_for_momentum_density_cutoff_d_; p | safety_factor_for_momentum_density_slope_; + p | enable_; } // WARNING! @@ -182,6 +184,9 @@ bool FixConservatives::operator()( const tnsr::II& inv_spatial_metric, const Scalar& sqrt_det_spatial_metric) const { bool needed_fixing = false; + if (not enable_) { + return needed_fixing; + } const size_t size = get<0>(tilde_b).size(); Variables, ::Tags::TempScalar<1>, ::Tags::TempScalar<2>, ::Tags::TempScalar<3>, @@ -395,7 +400,8 @@ bool operator==(const FixConservatives& lhs, const FixConservatives& rhs) { lhs.safety_factor_for_momentum_density_cutoff_d_ == rhs.safety_factor_for_momentum_density_cutoff_d_ and lhs.safety_factor_for_momentum_density_slope_ == - rhs.safety_factor_for_momentum_density_slope_; + rhs.safety_factor_for_momentum_density_slope_ and + lhs.enable_ == rhs.enable_; } bool operator!=(const FixConservatives& lhs, const FixConservatives& rhs) { diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.hpp index 0949a9ccb154..8ac5259601e3 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.hpp @@ -154,21 +154,28 @@ class FixConservatives { "Slope of safety factor for momentum density bound below " "SafetyFactorForSCutoffD, express as a function of log10(rho*W)."}; }; + /// Whether or not the limiting is enabled + struct Enable { + using type = bool; + static constexpr Options::String help = { + "If true then the limiting is applied."}; + }; - using options = tmpl::list; + using options = + tmpl::list; static constexpr Options::String help = { "Variable fixing used in Foucart's thesis.\n"}; - FixConservatives(const double minimum_rest_mass_density_times_lorentz_factor, - const double rest_mass_density_times_lorentz_factor_cutoff, - const double minimum_electron_fraction, - const double electron_fraction_cutoff, - const double safety_factor_for_magnetic_field, - const double safety_factor_for_momentum_density, - const double safety_factor_for_momentum_density_cutoff_d, - const double safety_factor_for_momentum_density_slope, + FixConservatives(double minimum_rest_mass_density_times_lorentz_factor, + double rest_mass_density_times_lorentz_factor_cutoff, + double minimum_electron_fraction, + double electron_fraction_cutoff, + double safety_factor_for_magnetic_field, + double safety_factor_for_momentum_density, + double safety_factor_for_momentum_density_cutoff_d, + double safety_factor_for_momentum_density_slope, bool enable, const Options::Context& context = {}); FixConservatives() = default; @@ -222,6 +229,7 @@ class FixConservatives { std::numeric_limits::signaling_NaN()}; double safety_factor_for_momentum_density_slope_{ std::numeric_limits::signaling_NaN()}; + bool enable_{true}; }; bool operator!=(const FixConservatives& lhs, const FixConservatives& rhs); diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml index d17a3c5278c3..f6fbdc55abb2 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml @@ -65,6 +65,7 @@ DomainCreator: VariableFixing: FixConservatives: + Enable: false # Only needed when not using Kastaun for recovery CutoffD: &CutoffD 1.1e-15 MinimumValueOfD: &MinimumD 1.0e-15 CutoffYe: 0.0 diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml index aba7d4c3ef96..ae35b24bf926 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml @@ -66,6 +66,7 @@ DomainCreator: VariableFixing: FixConservatives: + Enable: false # Only needed when not using Kastaun for recovery CutoffD: &CutoffD 1.0e-12 MinimumValueOfD: &MinimumD 1.0e-12 CutoffYe: 0.0 diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml index cc35359835f1..8d10d4f88e5a 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml @@ -59,6 +59,7 @@ DomainCreator: VariableFixing: FixConservatives: + Enable: false # Only needed when not using Kastaun for recovery CutoffD: &CutoffD 1.0e-12 MinimumValueOfD: &MinimumD 1.0e-12 CutoffYe: 0.0 diff --git a/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml b/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml index d9349d74b860..e519f29ff04a 100644 --- a/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml +++ b/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml @@ -90,6 +90,7 @@ EvolutionSystem: VariableFixing: FixConservatives: + Enable: false # Only needed when not using Kastaun for recovery CutoffD: &CutoffD 1.0e-12 MinimumValueOfD: &MinimumD 1.0e-12 CutoffYe: 0.0 diff --git a/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml b/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml index 0e37af80b15d..1b63dded4293 100644 --- a/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml +++ b/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml @@ -99,6 +99,7 @@ EvolutionSystem: VariableFixing: FixConservatives: + Enable: false # Only needed when not using Kastaun for recovery CutoffD: &CutoffD 1.0e-12 MinimumValueOfD: &MinimumD 1.0e-12 CutoffYe: 0.0 diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp index 84ecc3bfbe03..c6c3079d1238 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp @@ -50,7 +50,7 @@ SPECTRE_TEST_CASE( sqrt(get(determinant(spatial_metric)))}; const grmhd::ValenciaDivClean::FixConservatives variable_fixer{ - 1.e-7, 1.0e-7, 1.0e-10, 1.0e-10, 0.0, 0.0, 1.e-7, 0.0}; + 1.e-7, 1.0e-7, 1.0e-10, 1.0e-10, 0.0, 0.0, 1.e-7, 0.0, true}; typename System::variables_tag::type cons_vars{num_pts, 0.0}; get(get(cons_vars))[0] = 2.e-12; get(get(cons_vars))[0] = 2.e-13; diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp index 3b9157e722d6..5d4de75dd23e 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_FixConservativesAndComputePrims.cpp @@ -34,7 +34,7 @@ SPECTRE_TEST_CASE( } const Scalar sqrt_det_spatial_metric{num_pts, 1.0}; const grmhd::ValenciaDivClean::FixConservatives variable_fixer{ - 1.e-7, 1.0e-7, 1.0e-10, 1.0e-10, 0.0, 0.0, 1.e-7, 0.0}; + 1.e-7, 1.0e-7, 1.0e-10, 1.0e-10, 0.0, 0.0, 1.e-7, 0.0, true}; typename System::variables_tag::type cons_vars{num_pts, 0.0}; get(get(cons_vars))[0] = 2.e-12; get(get(cons_vars))[0] = 2.e-13; diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Test_FixConservatives.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Test_FixConservatives.cpp index 7ca09d499bea..716102752df6 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Test_FixConservatives.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Test_FixConservatives.cpp @@ -22,7 +22,8 @@ namespace { void test_variable_fixer( - const grmhd::ValenciaDivClean::FixConservatives& variable_fixer) { + const grmhd::ValenciaDivClean::FixConservatives& variable_fixer, + const bool enable) { // Call variable fixer at five points // [0]: tilde_d is too small, should be raised to limit // [1]: tilde_ye is too small, should be raised to limit @@ -45,17 +46,25 @@ void test_variable_fixer( } auto expected_tilde_d = tilde_d; - get(expected_tilde_d)[0] = 1.e-12; + if (enable) { + get(expected_tilde_d)[0] = 1.e-12; + } auto expected_tilde_ye = tilde_ye; - get(expected_tilde_ye)[0] = 1.e-13; // since Y_e = 0.1 - get(expected_tilde_ye)[1] = 1.e-10; + if (enable) { + get(expected_tilde_ye)[0] = 1.e-13; // since Y_e = 0.1 + get(expected_tilde_ye)[1] = 1.e-10; + } auto expected_tilde_tau = tilde_tau; - get(expected_tilde_tau)[2] = 2.0; + if (enable) { + get(expected_tilde_tau)[2] = 2.0; + } auto expected_tilde_s = tilde_s; - expected_tilde_s.get(0)[3] = sqrt(27.0); + if (enable) { + expected_tilde_s.get(0)[3] = sqrt(27.0); + } auto spatial_metric = make_with_value>(tilde_d, 0.0); @@ -68,9 +77,9 @@ void test_variable_fixer( inv_spatial_metric.get(d, d) = get(sqrt_det_spatial_metric); } - CHECK(variable_fixer(&tilde_d, &tilde_ye, &tilde_tau, &tilde_s, tilde_b, - spatial_metric, inv_spatial_metric, - sqrt_det_spatial_metric)); + CHECK(enable == variable_fixer(&tilde_d, &tilde_ye, &tilde_tau, &tilde_s, + tilde_b, spatial_metric, inv_spatial_metric, + sqrt_det_spatial_metric)); CHECK_ITERABLE_APPROX(tilde_d, expected_tilde_d); CHECK_ITERABLE_APPROX(tilde_ye, expected_tilde_ye); @@ -81,10 +90,12 @@ void test_variable_fixer( SPECTRE_TEST_CASE("Unit.Evolution.GrMhd.ValenciaDivClean.FixConservatives", "[VariableFixing][Unit]") { - grmhd::ValenciaDivClean::FixConservatives variable_fixer{ - 1.e-12, 1.0e-11, 1.0e-10, 1.0e-9, 0.0, 0.0, 1.e-12, 0.0}; - test_variable_fixer(variable_fixer); - test_serialization(variable_fixer); + for (const bool enable : {true, false}) { + grmhd::ValenciaDivClean::FixConservatives variable_fixer{ + 1.e-12, 1.0e-11, 1.0e-10, 1.0e-9, 0.0, 0.0, 1.e-12, 0.0, enable}; + test_variable_fixer(serialize_and_deserialize(variable_fixer), enable); + test_serialization(variable_fixer); + } const auto fixer_from_options = TestHelpers::test_creation( @@ -95,6 +106,7 @@ SPECTRE_TEST_CASE("Unit.Evolution.GrMhd.ValenciaDivClean.FixConservatives", "SafetyFactorForB: 0.0\n" "SafetyFactorForS: 0.0\n" "SafetyFactorForSCutoffD: 1.0e-12\n" - "SafetyFactorForSSlope: 0.0\n"); - test_variable_fixer(fixer_from_options); + "SafetyFactorForSSlope: 0.0\n" + "Enable: true\n"); + test_variable_fixer(fixer_from_options, true); } From c9d5cfc2e1b5345783b485ca0dc088c47c066edb Mon Sep 17 00:00:00 2001 From: Nils Deppe Date: Tue, 10 Oct 2023 15:40:25 -0700 Subject: [PATCH 3/5] Switch to 2d EOS for BNS --- .../GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp | 2 +- .../GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp b/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp index 9797e6b85101..2f8d76e754fb 100644 --- a/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp +++ b/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp @@ -273,7 +273,7 @@ struct get_thermodynamic_dim; template struct get_thermodynamic_dim { // Controls the thermodynamic dim used for numeric initial data. - static constexpr size_t value = 1; + static constexpr size_t value = 2; }; template diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml index f6fbdc55abb2..7f2ac901a88d 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml @@ -25,9 +25,11 @@ PhaseChangeAndTriggers: # Values taken from the Xcts/HeadOnBns.yaml input file EquationOfState: &eos - PolytropicFluid: - PolytropicConstant: 123.6489 + HybridEos(PolytropicFluid): + PolytropicFluid: + PolytropicConstant: 123.6 PolytropicExponent: 2. + ThermalAdiabaticIndex: 2.0 DomainCreator: BinaryCompactObject: From 1ba73fa5a080f913b16e14221e6dd7f567687f4d Mon Sep 17 00:00:00 2001 From: Nils Deppe Date: Mon, 9 Oct 2023 10:59:53 -0700 Subject: [PATCH 4/5] Add WCNS5Z GHMHD --- .../FiniteDifference/CMakeLists.txt | 2 + .../FiniteDifference/Factory.hpp | 1 + .../FiniteDifference/Reconstructor.hpp | 6 +- .../FiniteDifference/Wcns5z.cpp | 389 ++++++++++++++++++ .../FiniteDifference/Wcns5z.hpp | 204 +++++++++ .../GrMhd/GhValenciaDivClean/CMakeLists.txt | 1 + .../FiniteDifference/Test_Wcns5z.cpp | 41 ++ 7 files changed, 642 insertions(+), 2 deletions(-) create mode 100644 src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp create mode 100644 src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp create mode 100644 tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_Wcns5z.cpp diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/CMakeLists.txt b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/CMakeLists.txt index aff715bcda13..17ca92862912 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/CMakeLists.txt +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/CMakeLists.txt @@ -11,6 +11,7 @@ spectre_target_sources( PositivityPreservingAdaptiveOrder.cpp Reconstructor.cpp RegisterDerivedWithCharm.cpp + Wcns5z.cpp ) spectre_target_headers( @@ -30,4 +31,5 @@ spectre_target_headers( Reconstructor.hpp RegisterDerivedWithCharm.hpp Tag.hpp + Wcns5z.hpp ) diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Factory.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Factory.hpp index 50600272754d..05909d847ebb 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Factory.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Factory.hpp @@ -6,3 +6,4 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp" diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp index 144eef9f5e90..aacc87d2cf1f 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp @@ -19,6 +19,7 @@ namespace grmhd::GhValenciaDivClean::fd { /// \cond class MonotonisedCentralPrim; class PositivityPreservingAdaptiveOrderPrim; +class Wcns5zPrim; /// \endcond /*! @@ -38,8 +39,9 @@ class Reconstructor : public PUP::able { WRAPPED_PUPable_abstract(Reconstructor); // NOLINT /// \endcond - using creatable_classes = tmpl::list; + using creatable_classes = + tmpl::list; virtual std::unique_ptr get_clone() const = 0; diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp new file mode 100644 index 000000000000..a11c3a487def --- /dev/null +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp @@ -0,0 +1,389 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/FixedHashMap.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Structure/MaxNumberOfNeighbors.hpp" +#include "Domain/Structure/Side.hpp" +#include "Evolution/DgSubcell/GhostData.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.tpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" +#include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" +#include "NumericalAlgorithms/FiniteDifference/Reconstruct.tpp" +#include "NumericalAlgorithms/FiniteDifference/Unlimited.hpp" +#include "NumericalAlgorithms/FiniteDifference/Wcns5z.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "PointwiseFunctions/GeneralRelativity/Lapse.hpp" +#include "PointwiseFunctions/GeneralRelativity/Shift.hpp" +#include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/Gsl.hpp" + +namespace grmhd::GhValenciaDivClean::fd { + +Wcns5zPrim::Wcns5zPrim(const size_t nonlinear_weight_exponent, + const double epsilon, + const ::fd::reconstruction::FallbackReconstructorType + fallback_reconstructor, + const size_t max_number_of_extrema) + : nonlinear_weight_exponent_(nonlinear_weight_exponent), + epsilon_(epsilon), + fallback_reconstructor_(fallback_reconstructor), + max_number_of_extrema_(max_number_of_extrema) { + std::tie(reconstruct_, reconstruct_lower_neighbor_, + reconstruct_upper_neighbor_) = + ::fd::reconstruction::wcns5z_function_pointers<3>( + nonlinear_weight_exponent_, fallback_reconstructor_); +} + +Wcns5zPrim::Wcns5zPrim(CkMigrateMessage* const msg) : Reconstructor(msg) {} + +std::unique_ptr Wcns5zPrim::get_clone() const { + return std::make_unique(*this); +} + +void Wcns5zPrim::pup(PUP::er& p) { + Reconstructor::pup(p); + p | nonlinear_weight_exponent_; + p | epsilon_; + p | fallback_reconstructor_; + p | max_number_of_extrema_; + if (p.isUnpacking()) { + std::tie(reconstruct_, reconstruct_lower_neighbor_, + reconstruct_upper_neighbor_) = + ::fd::reconstruction::wcns5z_function_pointers<3>( + nonlinear_weight_exponent_, fallback_reconstructor_); + } +} + +// NOLINTNEXTLINE +PUP::able::PUP_ID Wcns5zPrim::my_PUP_ID = 0; + +template +void Wcns5zPrim::reconstruct( + const gsl::not_null, dim>*> + vars_on_lower_face, + const gsl::not_null, dim>*> + vars_on_upper_face, + const Variables>& volume_prims, + const Variables< + typename GhValenciaDivClean::System::variables_tag::type::tags_list>& + volume_spacetime_and_cons_vars, + const EquationsOfState::EquationOfState& eos, + const Element& element, + const FixedHashMap, ElementId>, + evolution::dg::subcell::GhostData, + boost::hash, ElementId>>>& + ghost_data, + const Mesh& subcell_mesh) const { + using all_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags:: + primitive_grmhd_and_spacetime_reconstruction_tags; + + FixedHashMap, ElementId>, + Variables, + boost::hash, ElementId>>> + neighbor_variables_data{}; + ::fd::neighbor_data_as_variables(make_not_null(&neighbor_variables_data), + ghost_data, ghost_zone_size(), + subcell_mesh); + + reconstruct_prims_work>, + prims_to_reconstruct_tags>( + vars_on_lower_face, vars_on_upper_face, + [this](auto upper_face_vars_ptr, auto lower_face_vars_ptr, + const auto& volume_vars, const auto& ghost_cell_vars, + const auto& subcell_extents, const size_t number_of_variables) { + reconstruct_(upper_face_vars_ptr, lower_face_vars_ptr, volume_vars, + ghost_cell_vars, subcell_extents, number_of_variables, + epsilon_, max_number_of_extrema_); + }, + [](auto upper_face_vars_ptr, auto lower_face_vars_ptr, + const auto& volume_vars, const auto& ghost_cell_vars, + const auto& subcell_extents, const size_t number_of_variables) { + ::fd::reconstruction::unlimited<4>( + upper_face_vars_ptr, lower_face_vars_ptr, volume_vars, + ghost_cell_vars, subcell_extents, number_of_variables); + }, + [](const auto vars_on_face_ptr) { + const auto& spacetime_metric = + get>(*vars_on_face_ptr); + auto& spatial_metric = + get>(*vars_on_face_ptr); + gr::spatial_metric(make_not_null(&spatial_metric), spacetime_metric); + auto& inverse_spatial_metric = + get>( + *vars_on_face_ptr); + auto& sqrt_det_spatial_metric = + get>(*vars_on_face_ptr); + + determinant_and_inverse(make_not_null(&sqrt_det_spatial_metric), + make_not_null(&inverse_spatial_metric), + spatial_metric); + get(sqrt_det_spatial_metric) = sqrt(get(sqrt_det_spatial_metric)); + + auto& shift = get>(*vars_on_face_ptr); + gr::shift(make_not_null(&shift), spacetime_metric, + inverse_spatial_metric); + gr::lapse( + make_not_null(&get>(*vars_on_face_ptr)), + shift, spacetime_metric); + }, + volume_prims, volume_spacetime_and_cons_vars, eos, element, + neighbor_variables_data, subcell_mesh, ghost_zone_size(), true); +} + +template +void Wcns5zPrim::reconstruct_fd_neighbor( + const gsl::not_null*> vars_on_face, + const Variables>& subcell_volume_prims, + const Variables< + grmhd::GhValenciaDivClean::Tags::spacetime_reconstruction_tags>& + subcell_volume_spacetime_metric, + const EquationsOfState::EquationOfState& eos, + const Element<3>& element, + const FixedHashMap< + maximum_number_of_neighbors(3), std::pair, ElementId<3>>, + evolution::dg::subcell::GhostData, + boost::hash, ElementId<3>>>>& ghost_data, + const Mesh<3>& subcell_mesh, + const Direction<3>& direction_to_reconstruct) const { + using prim_tags_for_reconstruction = + grmhd::GhValenciaDivClean::Tags::primitive_grmhd_reconstruction_tags; + using all_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags:: + primitive_grmhd_and_spacetime_reconstruction_tags; + + reconstruct_fd_neighbor_work( + vars_on_face, + [this](const auto tensor_component_on_face_ptr, + const auto& tensor_component_volume, + const auto& tensor_component_neighbor, + const Index& subcell_extents, + const Index& ghost_data_extents, + const Direction& local_direction_to_reconstruct) { + reconstruct_lower_neighbor_( + tensor_component_on_face_ptr, tensor_component_volume, + tensor_component_neighbor, subcell_extents, ghost_data_extents, + local_direction_to_reconstruct, epsilon_, max_number_of_extrema_); + }, + [](const auto tensor_component_on_face_ptr, + const auto& tensor_component_volume, + const auto& tensor_component_neighbor, + const Index& subcell_extents, + const Index& ghost_data_extents, + const Direction& local_direction_to_reconstruct) { + ::fd::reconstruction::reconstruct_neighbor< + Side::Lower, + ::fd::reconstruction::detail::UnlimitedReconstructor<4>>( + tensor_component_on_face_ptr, tensor_component_volume, + tensor_component_neighbor, subcell_extents, ghost_data_extents, + local_direction_to_reconstruct); + }, + [this](const auto tensor_component_on_face_ptr, + const auto& tensor_component_volume, + const auto& tensor_component_neighbor, + const Index& subcell_extents, + const Index& ghost_data_extents, + const Direction& local_direction_to_reconstruct) { + reconstruct_upper_neighbor_( + tensor_component_on_face_ptr, tensor_component_volume, + tensor_component_neighbor, subcell_extents, ghost_data_extents, + local_direction_to_reconstruct, epsilon_, max_number_of_extrema_); + }, + [](const auto tensor_component_on_face_ptr, + const auto& tensor_component_volume, + const auto& tensor_component_neighbor, + const Index& subcell_extents, + const Index& ghost_data_extents, + const Direction& local_direction_to_reconstruct) { + ::fd::reconstruction::reconstruct_neighbor< + Side::Upper, + ::fd::reconstruction::detail::UnlimitedReconstructor<4>>( + tensor_component_on_face_ptr, tensor_component_volume, + tensor_component_neighbor, subcell_extents, ghost_data_extents, + local_direction_to_reconstruct); + }, + [](const auto vars_on_face_ptr) { + const auto& spacetime_metric = + get>(*vars_on_face_ptr); + auto& spatial_metric = + get>(*vars_on_face_ptr); + gr::spatial_metric(make_not_null(&spatial_metric), spacetime_metric); + auto& inverse_spatial_metric = + get>( + *vars_on_face_ptr); + auto& sqrt_det_spatial_metric = + get>(*vars_on_face_ptr); + + determinant_and_inverse(make_not_null(&sqrt_det_spatial_metric), + make_not_null(&inverse_spatial_metric), + spatial_metric); + get(sqrt_det_spatial_metric) = sqrt(get(sqrt_det_spatial_metric)); + + auto& shift = get>(*vars_on_face_ptr); + gr::shift(make_not_null(&shift), spacetime_metric, + inverse_spatial_metric); + gr::lapse( + make_not_null(&get>(*vars_on_face_ptr)), + shift, spacetime_metric); + }, + subcell_volume_prims, subcell_volume_spacetime_metric, eos, element, + ghost_data, subcell_mesh, direction_to_reconstruct, ghost_zone_size(), + true); +} + +bool operator==(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { + // Don't check function pointers since they are set from + // nonlinear_weight_exponent_ and fallback_reconstructor_ + return lhs.nonlinear_weight_exponent_ == rhs.nonlinear_weight_exponent_ and + lhs.epsilon_ == rhs.epsilon_ and + lhs.fallback_reconstructor_ == rhs.fallback_reconstructor_ and + lhs.max_number_of_extrema_ == rhs.max_number_of_extrema_; +} + +bool operator!=(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { + return not(lhs == rhs); +} + +#define THERMO_DIM(data) BOOST_PP_TUPLE_ELEM(0, data) +#define TAGS_LIST_FD(data) \ + tmpl::list, \ + gh::Tags::Pi, gh::Tags::Phi, \ + ValenciaDivClean::Tags::TildeD, ValenciaDivClean::Tags::TildeYe, \ + ValenciaDivClean::Tags::TildeTau, \ + ValenciaDivClean::Tags::TildeS, \ + ValenciaDivClean::Tags::TildeB, \ + ValenciaDivClean::Tags::TildePhi, \ + hydro::Tags::RestMassDensity, \ + hydro::Tags::ElectronFraction, \ + hydro::Tags::SpecificInternalEnergy, \ + hydro::Tags::SpatialVelocity, \ + hydro::Tags::MagneticField, \ + hydro::Tags::DivergenceCleaningField, \ + hydro::Tags::LorentzFactor, \ + hydro::Tags::Pressure, \ + hydro::Tags::SpecificEnthalpy, \ + hydro::Tags::Temperature, \ + hydro::Tags::LorentzFactorTimesSpatialVelocity, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + ::Tags::Flux, \ + tmpl::size_t<3>, Frame::Inertial>, \ + ::Tags::Flux, \ + tmpl::size_t<3>, Frame::Inertial>, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + gr::Tags::Lapse, gr::Tags::Shift, \ + gr::Tags::SpatialMetric, \ + gr::Tags::SqrtDetSpatialMetric, \ + gr::Tags::InverseSpatialMetric, \ + evolution::dg::Actions::detail::NormalVector<3>> + +#define TAGS_LIST_DG_FD_INTERFACE(data) \ + tmpl::list, \ + gh::Tags::Pi, gh::Tags::Phi, \ + ValenciaDivClean::Tags::TildeD, ValenciaDivClean::Tags::TildeYe, \ + ValenciaDivClean::Tags::TildeTau, \ + ValenciaDivClean::Tags::TildeS, \ + ValenciaDivClean::Tags::TildeB, \ + ValenciaDivClean::Tags::TildePhi, \ + hydro::Tags::RestMassDensity, \ + hydro::Tags::ElectronFraction, \ + hydro::Tags::SpecificInternalEnergy, \ + hydro::Tags::SpatialVelocity, \ + hydro::Tags::MagneticField, \ + hydro::Tags::DivergenceCleaningField, \ + hydro::Tags::LorentzFactor, \ + hydro::Tags::Pressure, \ + hydro::Tags::SpecificEnthalpy, \ + hydro::Tags::Temperature, \ + hydro::Tags::LorentzFactorTimesSpatialVelocity, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + ::Tags::Flux, \ + tmpl::size_t<3>, Frame::Inertial>, \ + ::Tags::Flux, \ + tmpl::size_t<3>, Frame::Inertial>, \ + ::Tags::Flux, \ + Frame::Inertial>, \ + gh::ConstraintDamping::Tags::ConstraintGamma1, \ + gh::ConstraintDamping::Tags::ConstraintGamma2, \ + gr::Tags::Lapse, gr::Tags::Shift, \ + gr::Tags::SpatialMetric, \ + gr::Tags::SqrtDetSpatialMetric, \ + gr::Tags::InverseSpatialMetric, \ + evolution::dg::Actions::detail::NormalVector<3>> + +#define INSTANTIATION(r, data) \ + template void Wcns5zPrim::reconstruct( \ + gsl::not_null, 3>*> \ + vars_on_lower_face, \ + gsl::not_null, 3>*> \ + vars_on_upper_face, \ + const Variables>& volume_prims, \ + const Variables& \ + volume_spacetime_and_cons_vars, \ + const EquationsOfState::EquationOfState& eos, \ + const Element<3>& element, \ + const FixedHashMap, ElementId<3>>, \ + evolution::dg::subcell::GhostData, \ + boost::hash, ElementId<3>>>>& \ + ghost_data, \ + const Mesh<3>& subcell_mesh) const; \ + template void Wcns5zPrim::reconstruct_fd_neighbor( \ + gsl::not_null*> vars_on_face, \ + const Variables>& subcell_volume_prims, \ + const Variables< \ + grmhd::GhValenciaDivClean::Tags::spacetime_reconstruction_tags>& \ + subcell_volume_spacetime_metric, \ + const EquationsOfState::EquationOfState& eos, \ + const Element<3>& element, \ + const FixedHashMap, ElementId<3>>, \ + evolution::dg::subcell::GhostData, \ + boost::hash, ElementId<3>>>>& \ + ghost_data, \ + const Mesh<3>& subcell_mesh, \ + const Direction<3>& direction_to_reconstruct) const; + +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2)) + +#undef INSTANTIATION +#undef TAGS_LIST +#undef THERMO_DIM + +} // namespace grmhd::GhValenciaDivClean::fd diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp new file mode 100644 index 000000000000..88a82549963c --- /dev/null +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp @@ -0,0 +1,204 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include + +#include "DataStructures/FixedHashMap.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Domain/Structure/MaxNumberOfNeighbors.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" +#include "Options/Auto.hpp" +#include "Options/Context.hpp" +#include "Options/String.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/TMPL.hpp" + +/// \cond +class DataVector; +template +class Direction; +template +class Element; +template +class ElementId; +namespace EquationsOfState { +template +class EquationOfState; +} // namespace EquationsOfState +template +class Mesh; +namespace gsl { +template +class not_null; +} // namespace gsl +namespace PUP { +class er; +} // namespace PUP +template +class Variables; +namespace evolution::dg::subcell { +class GhostData; +} // namespace evolution::dg::subcell +/// \endcond + +namespace grmhd::GhValenciaDivClean::fd { +/*! + * \brief Fifth order weighted nonlinear compact scheme reconstruction using the + * Z oscillation indicator. See ::fd::reconstruction::wcns5z() for details. + * + */ +class Wcns5zPrim : public Reconstructor { + private: + using prims_to_reconstruct_tags = + tmpl::list, + hydro::Tags::ElectronFraction, + hydro::Tags::Temperature, + hydro::Tags::LorentzFactorTimesSpatialVelocity, + hydro::Tags::MagneticField, + hydro::Tags::DivergenceCleaningField>; + + using FallbackReconstructorType = + ::fd::reconstruction::FallbackReconstructorType; + + public: + static constexpr size_t dim = 3; + + struct NonlinearWeightExponent { + using type = size_t; + static constexpr Options::String help = { + "The exponent q to which the oscillation indicator term is raised"}; + }; + struct Epsilon { + using type = double; + static constexpr Options::String help = { + "The parameter added to the oscillation indicators to avoid division " + "by zero"}; + }; + struct FallbackReconstructor { + using type = FallbackReconstructorType; + static constexpr Options::String help = { + "A reconstruction scheme to fallback to adaptively. Finite difference " + "will switch to this reconstruction scheme if there are more extrema " + "in a FD stencil than a specified number. See also the option " + "'MaxNumberOfExtrema' below. Adaptive fallback is disabled if 'None'."}; + }; + struct MaxNumberOfExtrema { + using type = size_t; + static constexpr Options::String help = { + "The maximum allowed number of extrema in FD stencil for using Wcns5z " + "reconstruction before switching to a low-order reconstruction. If " + "FallbackReconstructor=None, this option is ignored"}; + }; + + using options = tmpl::list; + + static constexpr Options::String help{ + "WCNS 5Z reconstruction scheme using primitive variables."}; + + Wcns5zPrim() = default; + Wcns5zPrim(Wcns5zPrim&&) = default; + Wcns5zPrim& operator=(Wcns5zPrim&&) = default; + Wcns5zPrim(const Wcns5zPrim&) = default; + Wcns5zPrim& operator=(const Wcns5zPrim&) = default; + ~Wcns5zPrim() override = default; + + Wcns5zPrim(size_t nonlinear_weight_exponent, double epsilon, + FallbackReconstructorType fallback_reconstructor, + size_t max_number_of_extrema); + + explicit Wcns5zPrim(CkMigrateMessage* msg); + + WRAPPED_PUPable_decl_base_template(Reconstructor, Wcns5zPrim); + + auto get_clone() const -> std::unique_ptr override; + + static constexpr bool use_adaptive_order = false; + + void pup(PUP::er& p) override; + + size_t ghost_zone_size() const override { return 3; } + +using reconstruction_argument_tags = + tmpl::list<::Tags::Variables>, + typename System::variables_tag, + hydro::Tags::EquationOfStateBase, domain::Tags::Element, + evolution::dg::subcell::Tags::GhostDataForReconstruction, + evolution::dg::subcell::Tags::Mesh>; + + template + void reconstruct( + gsl::not_null, dim>*> vars_on_lower_face, + gsl::not_null, dim>*> vars_on_upper_face, + const Variables>& volume_prims, + const Variables& + volume_spacetime_and_cons_vars, + const EquationsOfState::EquationOfState& eos, + const Element& element, + const FixedHashMap< + maximum_number_of_neighbors(dim), + std::pair, ElementId>, + evolution::dg::subcell::GhostData, + boost::hash, ElementId>>>& ghost_data, + const Mesh& subcell_mesh) const; + + template + void reconstruct_fd_neighbor( + gsl::not_null*> vars_on_face, + const Variables>& subcell_volume_prims, + const Variables< + grmhd::GhValenciaDivClean::Tags::spacetime_reconstruction_tags>& + subcell_volume_spacetime_metric, + const EquationsOfState::EquationOfState& eos, + const Element& element, + const FixedHashMap< + maximum_number_of_neighbors(dim), + std::pair, ElementId>, + evolution::dg::subcell::GhostData, + boost::hash, ElementId>>>& ghost_data, + const Mesh& subcell_mesh, + const Direction& direction_to_reconstruct) const; + + private: + // NOLINTNEXTLINE(readability-redundant-declaration) + friend bool operator==(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs); + friend bool operator!=(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs); + + size_t nonlinear_weight_exponent_ = 0; + double epsilon_ = std::numeric_limits::signaling_NaN(); + FallbackReconstructorType fallback_reconstructor_ = + FallbackReconstructorType::None; + size_t max_number_of_extrema_ = 0; + + void (*reconstruct_)(gsl::not_null, dim>*>, + gsl::not_null, dim>*>, + const gsl::span&, + const DirectionMap>&, + const Index&, size_t, double, size_t) = nullptr; + void (*reconstruct_lower_neighbor_)(gsl::not_null, + const DataVector&, const DataVector&, + const Index&, const Index&, + const Direction&, const double&, + const size_t&) = nullptr; + void (*reconstruct_upper_neighbor_)(gsl::not_null, + const DataVector&, const DataVector&, + const Index&, const Index&, + const Direction&, const double&, + const size_t&) = nullptr; +}; + +} // namespace grmhd::GhValenciaDivClean::fd diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt index 573a8b552c60..642e1232e666 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt @@ -15,6 +15,7 @@ set(LIBRARY_SOURCES FiniteDifference/Test_FilterOptions.cpp FiniteDifference/Test_Filters.cpp FiniteDifference/Test_MonotonisedCentral.cpp + FiniteDifference/Test_Wcns5z.cpp Subcell/Test_FixConservativesAndComputePrims.cpp Subcell/Test_NeighborPackagedData.cpp Subcell/Test_PrimitiveGhostData.cpp diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_Wcns5z.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_Wcns5z.cpp new file mode 100644 index 000000000000..932f73895522 --- /dev/null +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_Wcns5z.cpp @@ -0,0 +1,41 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Tag.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp" +#include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" + +SPECTRE_TEST_CASE("Unit.Evolution.Systems.GrMhd.GhValenciaDivClean.Fd.Wcns5z", + "[Unit][Evolution]") { + namespace helpers = TestHelpers::grmhd::GhValenciaDivClean::fd; + PUPable_reg(SINGLE_ARG(grmhd::GhValenciaDivClean::fd::Wcns5zPrim)); + const auto wcns5z_from_options_base = TestHelpers::test_factory_creation< + grmhd::GhValenciaDivClean::fd::Reconstructor, + grmhd::GhValenciaDivClean::fd::OptionTags::Reconstructor>( + "Wcns5zPrim:\n" + " NonlinearWeightExponent: 2\n" + " Epsilon: 1.e-42\n" + " FallbackReconstructor: None\n" + " MaxNumberOfExtrema: 0\n"); + const auto wcns5z_deserialized = + serialize_and_deserialize(wcns5z_from_options_base); + auto* const wcns5z_from_options = + dynamic_cast( + wcns5z_deserialized.get()); + REQUIRE(wcns5z_from_options != nullptr); + CHECK(*wcns5z_from_options == + grmhd::GhValenciaDivClean::fd::Wcns5zPrim{ + 2, 1.e-42, fd::reconstruction::FallbackReconstructorType::None, 0}); + test_move_semantics( + grmhd::GhValenciaDivClean::fd::Wcns5zPrim{ + 2, 1.e-42, fd::reconstruction::FallbackReconstructorType::None, 0}, + grmhd::GhValenciaDivClean::fd::Wcns5zPrim{ + 2, 1.e-42, fd::reconstruction::FallbackReconstructorType::None, 0}); + helpers::test_prim_reconstructor(7, *wcns5z_from_options); +} From 30cfebdc442886aa319664a5fb94bc7cfe9739b4 Mon Sep 17 00:00:00 2001 From: Nils Deppe Date: Tue, 10 Oct 2023 19:18:54 -0700 Subject: [PATCH 5/5] Add PPAO test for GHMHD --- .../PositivityPreservingAdaptiveOrder.cpp | 2 + .../GrMhd/GhValenciaDivClean/CMakeLists.txt | 1 + ...Test_PositivityPreservingAdaptiveOrder.cpp | 47 +++++++++++++++++++ .../FiniteDifference/PrimReconstructor.hpp | 42 +++++++++++++++-- 4 files changed, 88 insertions(+), 4 deletions(-) create mode 100644 tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_PositivityPreservingAdaptiveOrder.cpp diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp index 333fad36621c..e0d270a120cb 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp @@ -66,9 +66,11 @@ PositivityPreservingAdaptiveOrderPrim::PositivityPreservingAdaptiveOrderPrim( PARSE_ERROR(context, "None is not an allowed low-order reconstructor."); } if (alpha_7.has_value()) { + PARSE_ERROR(context, "Alpha7 hasn't been tested."); six_to_the_alpha_7_ = pow(6.0, alpha_7.value()); } if (alpha_9.has_value()) { + PARSE_ERROR(context, "Alpha9 hasn't been tested."); eight_to_the_alpha_9_ = pow(8.0, alpha_9.value()); } set_function_pointers(); diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt index 642e1232e666..3b4100e24f08 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/CMakeLists.txt @@ -15,6 +15,7 @@ set(LIBRARY_SOURCES FiniteDifference/Test_FilterOptions.cpp FiniteDifference/Test_Filters.cpp FiniteDifference/Test_MonotonisedCentral.cpp + FiniteDifference/Test_PositivityPreservingAdaptiveOrder.cpp FiniteDifference/Test_Wcns5z.cpp Subcell/Test_FixConservativesAndComputePrims.cpp Subcell/Test_NeighborPackagedData.cpp diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_PositivityPreservingAdaptiveOrder.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_PositivityPreservingAdaptiveOrder.cpp new file mode 100644 index 000000000000..bde2d310eab9 --- /dev/null +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_PositivityPreservingAdaptiveOrder.cpp @@ -0,0 +1,47 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Tag.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp" +#include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" + +SPECTRE_TEST_CASE("Unit.Evolution.Systems.GrMhd.GhValenciaDivClean.Fd.Ppao", + "[Unit][Evolution]") { + namespace helpers = TestHelpers::grmhd::GhValenciaDivClean::fd; + PUPable_reg(SINGLE_ARG( + grmhd::GhValenciaDivClean::fd::PositivityPreservingAdaptiveOrderPrim)); + const auto ppao_from_options_base = TestHelpers::test_factory_creation< + grmhd::GhValenciaDivClean::fd::Reconstructor, + grmhd::GhValenciaDivClean::fd::OptionTags::Reconstructor>( + "PositivityPreservingAdaptiveOrderPrim:\n" + " Alpha5: 3.7\n" + " Alpha7: None\n" + " Alpha9: None\n" + " LowOrderReconstructor: MonotonisedCentral\n"); + const auto ppao_deserialized = + serialize_and_deserialize(ppao_from_options_base); + auto* const ppao_from_options = + dynamic_cast( + ppao_deserialized.get()); + REQUIRE(ppao_from_options != nullptr); + CHECK(*ppao_from_options == + grmhd::GhValenciaDivClean::fd::PositivityPreservingAdaptiveOrderPrim{ + 3.7, std::nullopt, std::nullopt, + fd::reconstruction::FallbackReconstructorType::MonotonisedCentral}); + test_move_semantics( + grmhd::GhValenciaDivClean::fd::PositivityPreservingAdaptiveOrderPrim{ + 3.7, std::nullopt, std::nullopt, + fd::reconstruction::FallbackReconstructorType::MonotonisedCentral}, + grmhd::GhValenciaDivClean::fd::PositivityPreservingAdaptiveOrderPrim{ + 3.7, std::nullopt, std::nullopt, + fd::reconstruction::FallbackReconstructorType::MonotonisedCentral}); + helpers::test_prim_reconstructor(10, *ppao_from_options); +} diff --git a/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp b/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp index acce8b77c589..8e8956cd5a4f 100644 --- a/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp +++ b/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp @@ -298,10 +298,44 @@ void test_prim_reconstructor_impl( } // Now we have everything to call the reconstruction - dynamic_cast(reconstructor) - .reconstruct(make_not_null(&vars_on_lower_face), - make_not_null(&vars_on_upper_face), volume_prims, - volume_cons_vars, eos, element, ghost_data, subcell_mesh); + if constexpr (Reconstructor::use_adaptive_order) { + std::array, 3> reconstruction_order_storage{}; + std::optional, 3>> + reconstruction_order{}; + reconstruction_order.emplace(); + for (size_t i = 0; i < 3; ++i) { + auto order_extents = subcell_mesh.extents(); + order_extents[i] += 2; + gsl::at(reconstruction_order_storage, i).resize(order_extents.product()); + // Ensure we have reset the values to max so the min calls are fine. + std::fill_n(gsl::at(reconstruction_order_storage, i).begin(), + order_extents.product(), + std::numeric_limits::max()); + gsl::at(reconstruction_order.value(), i) = gsl::span{ + gsl::at(reconstruction_order_storage, i).data(), + gsl::at(reconstruction_order_storage, i).size()}; + } + + dynamic_cast(reconstructor) + .reconstruct(make_not_null(&vars_on_lower_face), + make_not_null(&vars_on_upper_face), + make_not_null(&reconstruction_order), volume_prims, + volume_cons_vars, eos, element, ghost_data, subcell_mesh); + for (size_t d = 0; d < 3; ++d) { + CAPTURE(d); + for (size_t i = 0; i < gsl::at(reconstruction_order_storage, d).size(); + ++i) { + CAPTURE(i); + CHECK(gsl::at(reconstruction_order_storage, d)[i] >= 1); + CHECK(gsl::at(reconstruction_order_storage, d)[i] <= 9); + } + } + } else { + dynamic_cast(reconstructor) + .reconstruct(make_not_null(&vars_on_lower_face), + make_not_null(&vars_on_upper_face), volume_prims, + volume_cons_vars, eos, element, ghost_data, subcell_mesh); + } for (size_t dim = 0; dim < 3; ++dim) { CAPTURE(dim);