diff --git a/src/Evolution/VariableFixing/FixToAtmosphere.cpp b/src/Evolution/VariableFixing/FixToAtmosphere.cpp index 5a50506c6b2f..3727074bdaca 100644 --- a/src/Evolution/VariableFixing/FixToAtmosphere.cpp +++ b/src/Evolution/VariableFixing/FixToAtmosphere.cpp @@ -4,24 +4,25 @@ #include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include +#include #include #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Options/ParseError.hpp" #include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/Serialization/PupStlCpp17.hpp" namespace VariableFixing { template FixToAtmosphere::FixToAtmosphere( const double density_of_atmosphere, const double density_cutoff, - const double transition_density_cutoff, const double max_velocity_magnitude, + const std::optional velocity_limiting, const Options::Context& context) : density_of_atmosphere_(density_of_atmosphere), density_cutoff_(density_cutoff), - transition_density_cutoff_(transition_density_cutoff), - max_velocity_magnitude_(max_velocity_magnitude) { + velocity_limiting_(velocity_limiting) { if (density_of_atmosphere_ > density_cutoff_) { PARSE_ERROR(context, "The cutoff density (" << density_cutoff_ @@ -29,18 +30,55 @@ FixToAtmosphere::FixToAtmosphere( "density value in the atmosphere (" << density_of_atmosphere_ << ')'); } - if (transition_density_cutoff_ < density_of_atmosphere_ or - transition_density_cutoff_ > 10.0 * density_of_atmosphere_) { - PARSE_ERROR(context, "The transition density must be in [" - << density_of_atmosphere_ << ", " - << 10 * density_of_atmosphere_ << "], but is " - << transition_density_cutoff_); - } - if (transition_density_cutoff_ <= density_cutoff_) { - PARSE_ERROR(context, "The transition density cutoff (" - << transition_density_cutoff_ - << ") must be bigger than the density cutoff (" - << density_cutoff_ << ")"); + + if (velocity_limiting_.has_value()) { + if (velocity_limiting_->atmosphere_max_velocity < 0.0) { + PARSE_ERROR(context, + "The AtmosphereMaxVelocity must be non-negative but is " + << velocity_limiting_->atmosphere_max_velocity); + } + if (velocity_limiting_->near_atmosphere_max_velocity < 0.0) { + PARSE_ERROR(context, + "The NearAtmosphereMaxVelocity must be non-negative but is " + << velocity_limiting_->near_atmosphere_max_velocity); + } + if (velocity_limiting_->atmosphere_max_velocity > + velocity_limiting_->near_atmosphere_max_velocity) { + PARSE_ERROR(context, + "The AtmosphereMaxVelocity (" + << velocity_limiting_->atmosphere_max_velocity + << ") must be smaller NearAtmosphereMaxVelocity (" + << velocity_limiting_->near_atmosphere_max_velocity + << ")."); + } + if (velocity_limiting_->atmosphere_density_cutoff < 0.0) { + PARSE_ERROR(context, + "The AtmosphereDensityCutoff must be non-negative but is " + << velocity_limiting_->atmosphere_density_cutoff); + } + if (velocity_limiting_->transition_density_bound < 0.0) { + PARSE_ERROR(context, + "The TransitionDensityBound must be non-negative but is " + << velocity_limiting_->transition_density_bound); + } + if (velocity_limiting_->atmosphere_density_cutoff < + density_of_atmosphere_) { + PARSE_ERROR( + context, + "The AtmosphereDensityCutoff (" + << velocity_limiting_->atmosphere_density_cutoff + << ") must be greater than or equal to the DensityOfAtmosphere (" + << density_of_atmosphere_ << ")."); + } + if (velocity_limiting_->transition_density_bound < + velocity_limiting_->atmosphere_density_cutoff) { + PARSE_ERROR(context, "The TransitionDensityBound (" + << velocity_limiting_->transition_density_bound + << ") must be greater than or equal to the " + "AtmosphereDensityCutoff (" + << velocity_limiting_->atmosphere_density_cutoff + << ")."); + } } } @@ -49,8 +87,7 @@ template void FixToAtmosphere::pup(PUP::er& p) { p | density_of_atmosphere_; p | density_cutoff_; - p | transition_density_cutoff_; - p | max_velocity_magnitude_; + p | velocity_limiting_; } template @@ -72,14 +109,10 @@ void FixToAtmosphere::operator()( set_density_to_atmosphere(rest_mass_density, specific_internal_energy, temperature, pressure, electron_fraction, equation_of_state, i); - for (size_t d = 0; d < Dim; ++d) { - spatial_velocity->get(d)[i] = 0.0; - } - get(*lorentz_factor)[i] = 1.0; - } else if (UNLIKELY(rest_mass_density->get()[i] < - transition_density_cutoff_)) { - set_to_magnetic_free_transition(spatial_velocity, lorentz_factor, - *rest_mass_density, spatial_metric, i); + } + if (velocity_limiting_.has_value()) { + apply_velocity_limit(spatial_velocity, lorentz_factor, *rest_mass_density, + spatial_metric, i); } // For 2D & 3D EoS, we also need to limit the temperature / energy @@ -179,13 +212,30 @@ void FixToAtmosphere::set_density_to_atmosphere( } template -void FixToAtmosphere::set_to_magnetic_free_transition( +void FixToAtmosphere::apply_velocity_limit( const gsl::not_null*> spatial_velocity, const gsl::not_null*> lorentz_factor, const Scalar& rest_mass_density, const tnsr::ii& spatial_metric, const size_t grid_index) const { + if (get(rest_mass_density)[grid_index] > + velocity_limiting_->transition_density_bound) { + return; + } + + if (get(rest_mass_density)[grid_index] < + velocity_limiting_->atmosphere_density_cutoff) { + for (size_t i = 0; i < Dim; ++i) { + spatial_velocity->get(i)[grid_index] = + velocity_limiting_->atmosphere_max_velocity; + } + if (LIKELY(velocity_limiting_->atmosphere_max_velocity == 0.0)) { + get(*lorentz_factor)[grid_index] = 1.0; + return; + } + } + double magnitude_of_velocity = 0.0; for (size_t j = 0; j < Dim; ++j) { magnitude_of_velocity += spatial_velocity->get(j)[grid_index] * @@ -197,15 +247,24 @@ void FixToAtmosphere::set_to_magnetic_free_transition( spatial_metric.get(j, k)[grid_index]; } } + if (UNLIKELY(get(rest_mass_density)[grid_index] < + velocity_limiting_->atmosphere_density_cutoff)) { + // Note: magnitude_of_velocity is squared still. + get(*lorentz_factor)[grid_index] = 1.0 / sqrt(1.0 - magnitude_of_velocity); + return; + } magnitude_of_velocity = sqrt(magnitude_of_velocity); - const double scale_factor = - (get(rest_mass_density)[grid_index] - density_cutoff_) / - (transition_density_cutoff_ - density_cutoff_); - if (const double max_mag_of_velocity = scale_factor * max_velocity_magnitude_; + const double scale_factor = (get(rest_mass_density)[grid_index] - + velocity_limiting_->atmosphere_density_cutoff) / + (velocity_limiting_->transition_density_bound - + velocity_limiting_->atmosphere_density_cutoff); + if (const double max_mag_of_velocity = + scale_factor * velocity_limiting_->near_atmosphere_max_velocity; magnitude_of_velocity > max_mag_of_velocity) { + const double one_over_max_mag_of_velocity = 1.0 / magnitude_of_velocity; for (size_t j = 0; j < Dim; ++j) { spatial_velocity->get(j)[grid_index] *= - max_mag_of_velocity / magnitude_of_velocity; + max_mag_of_velocity * one_over_max_mag_of_velocity; } get(*lorentz_factor)[grid_index] = 1.0 / sqrt(1.0 - max_mag_of_velocity * max_mag_of_velocity); @@ -217,8 +276,7 @@ bool operator==(const FixToAtmosphere& lhs, const FixToAtmosphere& rhs) { return lhs.density_of_atmosphere_ == rhs.density_of_atmosphere_ and lhs.density_cutoff_ == rhs.density_cutoff_ and - lhs.transition_density_cutoff_ == rhs.transition_density_cutoff_ and - lhs.max_velocity_magnitude_ == rhs.max_velocity_magnitude_; + lhs.velocity_limiting_ == rhs.velocity_limiting_; } template @@ -227,6 +285,29 @@ bool operator!=(const FixToAtmosphere& lhs, return not(lhs == rhs); } +template +void FixToAtmosphere::VelocityLimitingOptions::pup(PUP::er& p) { + p | atmosphere_max_velocity; + p | near_atmosphere_max_velocity; + p | atmosphere_density_cutoff; + p | transition_density_bound; +} + +template +bool FixToAtmosphere::VelocityLimitingOptions::operator==( + const VelocityLimitingOptions& rhs) const { + return atmosphere_max_velocity == rhs.atmosphere_max_velocity and + near_atmosphere_max_velocity == rhs.near_atmosphere_max_velocity and + atmosphere_density_cutoff == rhs.atmosphere_density_cutoff and + transition_density_bound == rhs.transition_density_bound; +} + +template +bool FixToAtmosphere::VelocityLimitingOptions::operator!=( + const VelocityLimitingOptions& rhs) const { + return not(*this == rhs); +} + #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) #define THERMO_DIM(data) BOOST_PP_TUPLE_ELEM(1, data) diff --git a/src/Evolution/VariableFixing/FixToAtmosphere.hpp b/src/Evolution/VariableFixing/FixToAtmosphere.hpp index a6878dcadcf0..e75581e7f0ea 100644 --- a/src/Evolution/VariableFixing/FixToAtmosphere.hpp +++ b/src/Evolution/VariableFixing/FixToAtmosphere.hpp @@ -5,8 +5,10 @@ #include #include +#include #include "DataStructures/Tensor/TypeAliases.hpp" +#include "Options/Auto.hpp" #include "Options/Context.hpp" #include "Options/String.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" @@ -35,23 +37,7 @@ namespace VariableFixing { * (DensityOfAtmosphere), and the pressure, and specific internal energy (for * one-dimensional equations of state) are adjusted to satisfy the equation of * state. For a two-dimensional equation of state, the specific internal energy - * is set to zero. In addition, the spatial velocity is set to zero, and the - * Lorentz factor is set to one. - * - * If the rest mass density is above \f$\rho_{\textrm{cutoff}}\f$ but below - * \f$\rho_{\textrm{transition}}\f$ (TransitionDensityCutoff) then the velocity - * is rescaled such that - * - * \f{align*}{ - * \sqrt{v^i v_i}\le \frac{(\rho-\rho_{\textrm{cutoff}})} - * {(\rho_{\textrm{transition}} - \rho_{\textrm{cutoff}})} v_{\max} - * \f} - * - * where \f$v_{\max}\f$ (MaxVelocityMagnitude) is the maximum allowed magnitude - * of the velocity. This prescription follows Appendix 2.d of - * \cite Muhlberger2014pja Note that we require - * \f$\rho_{\textrm{transition}}\in(\rho_{\textrm{cutoff}}, - * 10\rho_{\textrm{atm}}]\f$ + * is set to zero. */ template class FixToAtmosphere { @@ -70,34 +56,92 @@ class FixToAtmosphere { static constexpr Options::String help = { "Density to impose atmosphere at. Must be >= rho_atm"}; }; - /// \brief For densities between DensityOfAtmosphere and - /// TransitionDensityCutoff the velocity is transitioned away from atmosphere - /// to avoid abrupt cutoffs. - /// - /// This value must not be larger than `10 * DensityOfAtmosphere`. - struct TransitionDensityCutoff { - using type = double; - static type lower_bound() { return 0.0; } + + /*! + * \brief Limit the velocity in and near the atmosphere. + * + * Let $v_{\max}$ be the maximum magnitude of the + * velocity near the atmosphere, which we typically set to $10^{-4}$. + * We let $v_{\mathrm{atm}}$ be the maximum magnitude of the velocity + * in the atmosphere, which we typically set to $0$. We then define + * the maximum magnitude of the spatial velocity to be + * + * \f{align*}{ + * \tilde{v} + * &=\begin{cases} + * v_{\mathrm{atm}}, & \mathrm{if}\; \rho < \rho_{v^-} \ \ + * v_{\mathrm{atm}} + \left(v_{\max} - v_{\mathrm{atm}}\right) + * \frac{\rho - \rho_{v^-}}{\rho_{v^+} - \rho_{v^-}}, + * & \mathrm{if}\;\rho_{v^-} \le \rho < \rho_{v^+} + * \end{cases} + * \f} + * + * We then rescale the velocity by + * + * \f{align*}{ + * v^i\to v^i\frac{\tilde{v}}{\sqrt{v^i\gamma_{ij}v^j}}. + * \f} + */ + struct VelocityLimitingOptions { + struct AtmosphereMaxVelocity { + using type = double; + static constexpr Options::String help = { + "The maximum velocity magnitude IN the atmosphere. Typically set to " + "0."}; + }; + + struct NearAtmosphereMaxVelocity { + using type = double; + static constexpr Options::String help = { + "The maximum velocity magnitude NEAR the atmosphere. Typically set " + "to 1e-4."}; + }; + + struct AtmosphereDensityCutoff { + using type = double; + static constexpr Options::String help = { + "The rest mass density cutoff below which the velocity magnitude is " + "limited to AtmosphereMaxVelocity. Typically set to " + "(10 or 20)*DensityOfAtmosphere."}; + }; + + struct TransitionDensityBound { + using type = double; + static constexpr Options::String help = { + "The rest mass density above which no velocity limiting is done. " + "Between " + "this value and AtmosphereDensityCutoff a linear transition in the " + "maximum magnitude of the velocity between AtmosphereMaxVelocity and " + "NearAtmosphereMaxVelocity is done. Typically set to " + "10*AtmosphereDensityCutoff."}; + }; + using options = tmpl::list; static constexpr Options::String help = { - "For densities between DensityOfAtmosphere and TransitionDensityCutoff " - "the velocity is transitioned away from atmosphere to avoid abrupt " - "cutoffs.\n\n" - "This value must not be larger than 10 * DensityOfAtmosphere."}; + "Limit the velocity in and near the atmosphere."}; + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& p); + + bool operator==(const VelocityLimitingOptions& rhs) const; + bool operator!=(const VelocityLimitingOptions& rhs) const; + + double atmosphere_max_velocity{ + std::numeric_limits::signaling_NaN()}; + double near_atmosphere_max_velocity{ + std::numeric_limits::signaling_NaN()}; + double atmosphere_density_cutoff{ + std::numeric_limits::signaling_NaN()}; + double transition_density_bound{ + std::numeric_limits::signaling_NaN()}; }; - /// \brief The maximum magnitude of the velocity when the density is below - /// `TransitionDensityCutoff` - struct MaxVelocityMagnitude { - using type = double; - static type lower_bound() { return 0.0; } - static type upper_bound() { return 1.0; } - static constexpr Options::String help = { - "The maximum sqrt(v^i v^j gamma_{ij}) allowed when the density is " - "below TransitionDensityCutoff."}; + struct VelocityLimiting { + using type = Options::Auto; + static constexpr Options::String help = VelocityLimitingOptions::help; }; using options = - tmpl::list; + tmpl::list; static constexpr Options::String help = { "If the rest mass density is below DensityCutoff, it is set\n" "to DensityOfAtmosphere, and the pressure, and specific internal energy\n" @@ -108,8 +152,7 @@ class FixToAtmosphere { "factor is set to one.\n"}; FixToAtmosphere(double density_of_atmosphere, double density_cutoff, - double transition_density_cutoff, - double max_velocity_magnitude, + std::optional velocity_limiting, const Options::Context& context = {}); FixToAtmosphere() = default; @@ -160,7 +203,7 @@ class FixToAtmosphere { equation_of_state, size_t grid_index) const; - void set_to_magnetic_free_transition( + void apply_velocity_limit( gsl::not_null*> spatial_velocity, gsl::not_null*> lorentz_factor, @@ -175,9 +218,7 @@ class FixToAtmosphere { double density_of_atmosphere_{std::numeric_limits::signaling_NaN()}; double density_cutoff_{std::numeric_limits::signaling_NaN()}; - double transition_density_cutoff_{ - std::numeric_limits::signaling_NaN()}; - double max_velocity_magnitude_{std::numeric_limits::signaling_NaN()}; + std::optional velocity_limiting_{std::nullopt}; }; template diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml index effc521bc814..a8595c2bbf2e 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml @@ -4,6 +4,7 @@ Executable: EvolveGhValenciaDivCleanBns Testing: Check: parse + Timeout: 8 Priority: High --- @@ -89,8 +90,11 @@ VariableFixing: FixToAtmosphere: DensityOfAtmosphere: 1.0e-15 DensityCutoff: &AtmosphereDensityCutoff 1.1e-15 - TransitionDensityCutoff: 9.0e-15 - MaxVelocityMagnitude: 1.0e-4 + VelocityLimiting: + AtmosphereMaxVelocity: 0 + NearAtmosphereMaxVelocity: 1.0e-4 + AtmosphereDensityCutoff: 3.0e-15 + TransitionDensityBound: 3.0e-14 LimitLorentzFactor: Enable: false # Only needed when not using Kastaun for recovery MaxDensityCutoff: 1.0e-08 diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml index f7ff71a390de..4bd47a67e8e5 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdBondiMichel.yaml @@ -82,8 +82,11 @@ VariableFixing: FixToAtmosphere: DensityOfAtmosphere: 1.0e-12 DensityCutoff: 1.0e-12 - TransitionDensityCutoff: 1.0e-11 - MaxVelocityMagnitude: 1.0e-4 + VelocityLimiting: + AtmosphereMaxVelocity: 0 + NearAtmosphereMaxVelocity: 1.0e-4 + AtmosphereDensityCutoff: 3.0e-12 + TransitionDensityBound: 1.0e-11 LimitLorentzFactor: Enable: false # Only needed when not using Kastaun for recovery MaxDensityCutoff: 1.0e-12 diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml index 4d726575b5ee..a8fc3b6df7b2 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/GhMhdTovStar.yaml @@ -80,8 +80,11 @@ VariableFixing: FixToAtmosphere: DensityOfAtmosphere: 1.0e-15 DensityCutoff: &AtmosphereDensityCutoff 1.1e-15 - TransitionDensityCutoff: 9.0e-15 - MaxVelocityMagnitude: 1.0e-4 + VelocityLimiting: + AtmosphereMaxVelocity: 0 + NearAtmosphereMaxVelocity: 1.0e-4 + AtmosphereDensityCutoff: 3.0e-15 + TransitionDensityBound: 3.0e-14 LimitLorentzFactor: Enable: false # Only needed when not using Kastaun for recovery MaxDensityCutoff: 1.0e-12 diff --git a/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml b/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml index 23644081df90..3d8181796efc 100644 --- a/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml +++ b/tests/InputFiles/GrMhd/ValenciaDivClean/BlastWave.yaml @@ -4,6 +4,7 @@ Executable: EvolveValenciaDivClean Testing: Check: parse;execute + Timeout: 8 Priority: High --- @@ -113,8 +114,11 @@ VariableFixing: FixToAtmosphere: DensityOfAtmosphere: 1.0e-12 DensityCutoff: 1.0e-12 - TransitionDensityCutoff: 1.0e-11 - MaxVelocityMagnitude: 1.0e-4 + VelocityLimiting: + AtmosphereMaxVelocity: 0 + NearAtmosphereMaxVelocity: 1.0e-4 + AtmosphereDensityCutoff: 3.0e-12 + TransitionDensityBound: 1.0e-11 PrimitiveFromConservative: CutoffDForInversion: *CutoffD diff --git a/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml b/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml index 6be821c219c3..bcb2a8cdbc46 100644 --- a/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml +++ b/tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml @@ -131,8 +131,11 @@ VariableFixing: FixToAtmosphere: DensityOfAtmosphere: 1.0e-12 DensityCutoff: 1.0e-12 - TransitionDensityCutoff: 1.0e-11 - MaxVelocityMagnitude: 1.0e-4 + VelocityLimiting: + AtmosphereMaxVelocity: 0 + NearAtmosphereMaxVelocity: 1.0e-4 + AtmosphereDensityCutoff: 3.0e-12 + TransitionDensityBound: 1.0e-11 PrimitiveFromConservative: CutoffDForInversion: *CutoffD diff --git a/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp b/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp index 8efea9ae439f..f0df86dabf14 100644 --- a/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp +++ b/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp @@ -142,9 +142,11 @@ void test_variable_fixer( template void test_variable_fixer() { + using Vlo = + typename VariableFixing::FixToAtmosphere::VelocityLimitingOptions; // Test for representative 1-d equation of state - VariableFixing::FixToAtmosphere variable_fixer{1.e-12, 3.e-12, 1.e-11, - 1.e-4}; + const VariableFixing::FixToAtmosphere variable_fixer{ + 1.e-12, 3.e-12, Vlo{0.0, 1.e-4, 3.e-12, 1.e-11}}; EquationsOfState::PolytropicFluid polytrope{1.0, 2.0}; test_variable_fixer(variable_fixer, polytrope); test_serialization(variable_fixer); @@ -153,8 +155,11 @@ void test_variable_fixer() { TestHelpers::test_creation>( "DensityOfAtmosphere: 1.0e-12\n" "DensityCutoff: 3.0e-12\n" - "TransitionDensityCutoff: 1.0e-11\n" - "MaxVelocityMagnitude: 1.0e-4\n"); + "VelocityLimiting:\n" + " AtmosphereMaxVelocity: 0\n" + " NearAtmosphereMaxVelocity: 1.0e-4\n" + " AtmosphereDensityCutoff: 3.0e-12\n" + " TransitionDensityBound: 1.0e-11\n"); test_variable_fixer(fixer_from_options, polytrope); // Test for representative 2-d equation of state