diff --git a/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp b/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp index 6a7810ed545e..22106ff9eeb8 100644 --- a/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp +++ b/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp @@ -641,12 +641,47 @@ class EquationOfState : public PUP::able { * Computes adiabatic sound speed squared * \f[ * c_s^2 \equiv \frac{\partial p}{\partial e} |_{s, Y_e} = - * \frac{\rho}{h}\frac{\partial p}{\rho} |_{e, Y_e} + + * \frac{\rho}{h}\frac{\partial p}{\partial \rho} |_{e, Y_e} + * \frac{\partial p}{\partial e}|_{\rho, Y_e} * \f]. * With \f$p, e\f$ the pressure and energy density respectively, * \f$s\f$ the entropy density, \f$Y_e\f$ the electron fraction - * and \f$\rho\f$ the rest-mass density. + * \f$\rho\f$ the rest-mass density, and \f$h\f$ the enthalpy density. + * Note that \f$e\f$ is the total energy density and not the internal energy, + * therefore + * \f[ + * \frac{\partial p}{\partial \rho} |_{e, Y_e} \neq \chi \equiv \frac{\partial + * p}{\partial \rho} |_{\epsilon, Y_e} + * \f] + * as defined in the 2-d EoS above. By definition + * \f$ e = (1+\epsilon) \rho \f$ so holding \f$e\f$ constant + * \f[ + * 0 = \frac{d e}{d \rho} = \frac{\partial e}{\partial \rho} + + * \frac{\partial e}{\partial \epsilon} \frac{\partial \epsilon}{\rho}. + * \f] + * (where we have suppressed \f$ Y_e\f$ dependence) + * So \f$ \partial \epsilon / \partial \rho |_{e} = (1 + \epsilon)/\rho \f$ + * and we can expand + * \f[ + * \frac{\partial p}{\partial \rho} |_{e, Y_e} = \frac{\partial e}{\partial + * \rho}_{\epsilon, Y_e} + \frac{(1 + \epsilon)}{\rho} \frac{\partial + * e}{\partial \epsilon}|_{\rho, Y_e} + * \f] + * Finally, we can rewrite the entire sound speed using only the rest-mass + * density, specific internal energy, and electron fraction as variables, + * by using \f$ \frac{\partial e}{\partial \epsilon}|_{\rho, Y_e} = 1 \f$ + * \f[ + * c_s^2 = + * \frac{\rho}{h}\frac{\partial p}{\partial \rho} |_{\epsilon, Y_e} + + * \frac{1}{\rho} \frac{\partial p}{\partial \epsilon}|_{\rho, Y_e} \left( + * 1 - \frac{(1 + \epsilon)\rho}{h}\right) + * \f] + * Which reduces to our preferred form + * \f[ + * c_s^2 = + * \frac{\rho}{h}(\chi + \kappa) + * \f] + * * Computed as a function of temperature, rest-mass density and electron * fraction. Note that this will break thermodynamic consistency if the * pressure and internal energy interpolated separately. The precise impact of diff --git a/src/PointwiseFunctions/Hydro/EquationsOfState/Equilibrium3D.cpp b/src/PointwiseFunctions/Hydro/EquationsOfState/Equilibrium3D.cpp index 8d8706c9e65f..b36b3253a810 100644 --- a/src/PointwiseFunctions/Hydro/EquationsOfState/Equilibrium3D.cpp +++ b/src/PointwiseFunctions/Hydro/EquationsOfState/Equilibrium3D.cpp @@ -119,14 +119,12 @@ Equilibrium3D::sound_speed_squared_from_density_and_temperature_impl( // Cold part plus temperature dependent part, see the 3D EoS documentation for // expression. return Scalar{ - get(rest_mass_density) * - get(underlying_eos_.chi_from_density_and_energy( - rest_mass_density, specific_internal_energy)) / - enthalpy_density + - get(rest_mass_density) / pressure * - get(underlying_eos_ - .kappa_times_p_over_rho_squared_from_density_and_energy( - rest_mass_density, specific_internal_energy))}; + get(rest_mass_density) / enthalpy_density * + (get(underlying_eos_.chi_from_density_and_energy( + rest_mass_density, specific_internal_energy)) + + get(underlying_eos_ + .kappa_times_p_over_rho_squared_from_density_and_energy( + rest_mass_density, specific_internal_energy)))}; } template class Equilibrium3D>>; diff --git a/tests/Unit/PointwiseFunctions/Hydro/EquationsOfState/Test_Equilibrium3D.cpp b/tests/Unit/PointwiseFunctions/Hydro/EquationsOfState/Test_Equilibrium3D.cpp index 728afe7ebda0..59c064a64b3b 100644 --- a/tests/Unit/PointwiseFunctions/Hydro/EquationsOfState/Test_Equilibrium3D.cpp +++ b/tests/Unit/PointwiseFunctions/Hydro/EquationsOfState/Test_Equilibrium3D.cpp @@ -16,15 +16,22 @@ #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/Equilibrium3D.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp" #include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" SPECTRE_TEST_CASE("Unit.PointwiseFunctions.EquationsOfState.Equilibrium3D", "[Unit][EquationsOfState]") { namespace EoS = EquationsOfState; - const EoS::PolytropicFluid cold_eos{100.0, 2.0}; - const EoS::HybridEos> underlying_eos{cold_eos, - 2.0}; + + const double ideal_adiabatic_index = 2.0; + const double cold_polytropic_index = 2.0; + const EoS::PolytropicFluid cold_eos{100.0, cold_polytropic_index}; + const EoS::IdealFluid eos_ideal_fluid{ideal_adiabatic_index}; + + EoS::Equilibrium3D> eos_ideal_fluid_3d{eos_ideal_fluid}; + const EoS::HybridEos> underlying_eos{ + cold_eos, ideal_adiabatic_index}; EoS::Equilibrium3D>> eos{ underlying_eos}; CHECK(eos.rest_mass_density_lower_bound() == @@ -83,22 +90,32 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.EquationsOfState.Equilibrium3D", eos.temperature_from_density_and_energy( rest_mass_density, specific_internal_energy, electron_fraction), temperature); + } + { + const Scalar rest_mass_density{ + DataVector{1e-10, 1e-6, 1e-4, 1e-3, 1e-2, 1.0}}; + // electron fraction should be irrelevant + const Scalar electron_fraction{ + DataVector{1e-10, 1e-6, 1e-4, 1e-3, 1e-2, 1.0}}; + const Scalar temperature{ + DataVector{1e-10, 1e-6, 1e-4, 1e-3, 1e-2, 1.0}}; + // For the ideal fluid the sound speed can be computed analytically + const Scalar pressure = + eos_ideal_fluid_3d.pressure_from_density_and_temperature( + rest_mass_density, temperature, electron_fraction); + const Scalar specific_internal_energy = + eos_ideal_fluid_3d + .specific_internal_energy_from_density_and_temperature( + rest_mass_density, temperature, electron_fraction); const DataVector enthalpy_density = get(rest_mass_density) * (1.0 + get(specific_internal_energy)) + get(pressure); + CHECK_ITERABLE_APPROX( - get(eos.sound_speed_squared_from_density_and_temperature( + get(eos_ideal_fluid_3d.sound_speed_squared_from_density_and_temperature( rest_mass_density, temperature, electron_fraction)), - get(rest_mass_density) * - get(underlying_eos.chi_from_density_and_energy( - rest_mass_density, specific_internal_energy)) / - enthalpy_density + - get(rest_mass_density) / get(pressure) * - get(underlying_eos - .kappa_times_p_over_rho_squared_from_density_and_energy( - rest_mass_density, specific_internal_energy))); + (ideal_adiabatic_index)*get(pressure) / enthalpy_density); } - { // double functions const Scalar rest_mass_density{{1e-3}}; @@ -121,42 +138,39 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.EquationsOfState.Equilibrium3D", rest_mass_density, specific_internal_energy, electron_fraction) == underlying_eos.temperature_from_density_and_energy( rest_mass_density, specific_internal_energy)); + } + { + // Using the ideal fluid for sound speed computations + const Scalar rest_mass_density{{1e-3}}; + // temperature and electron fraction should be irrelevant + const Scalar electron_fraction{{1e-1}}; + const Scalar temperature{{1e-1}}; + const Scalar pressure = + eos_ideal_fluid_3d.pressure_from_density_and_temperature( + rest_mass_density, temperature, electron_fraction); + const Scalar specific_internal_energy = + eos_ideal_fluid_3d + .specific_internal_energy_from_density_and_temperature( + rest_mass_density, temperature, electron_fraction); + const double enthalpy_density = get(rest_mass_density) * (1.0 + get(specific_internal_energy)) + get(pressure); - CHECK( - get(eos.sound_speed_squared_from_density_and_temperature( - rest_mass_density, temperature, electron_fraction)) == - approx( - get(rest_mass_density) * - get(underlying_eos.chi_from_density_and_energy( - rest_mass_density, specific_internal_energy)) / - enthalpy_density + - get(rest_mass_density) / get(pressure) * - get(underlying_eos - .kappa_times_p_over_rho_squared_from_density_and_energy( - rest_mass_density, specific_internal_energy)))); + // IdealFluid value is known analytically + CHECK_ITERABLE_APPROX( + get(eos_ideal_fluid_3d.sound_speed_squared_from_density_and_temperature( + rest_mass_density, temperature, electron_fraction)), + (ideal_adiabatic_index)*get(pressure) / enthalpy_density); // Check that the sound speed calculation works at zero temperature CHECK_ITERABLE_APPROX( get(eos.sound_speed_squared_from_density_and_temperature( - rest_mass_density, Scalar{0.0}, electron_fraction)) - - get(rest_mass_density) / - get(cold_eos.pressure_from_density(rest_mass_density)) * - get(underlying_eos - .kappa_times_p_over_rho_squared_from_density_and_energy( - rest_mass_density, - cold_eos.specific_internal_energy_from_density( - rest_mass_density))), - get(rest_mass_density) * - get(underlying_eos.chi_from_density_and_energy( - rest_mass_density, - cold_eos.specific_internal_energy_from_density( - rest_mass_density))) / - (get(rest_mass_density) + - get(cold_eos.pressure_from_density(rest_mass_density)) + - get(cold_eos.specific_internal_energy_from_density( - rest_mass_density)) * - get(rest_mass_density))); + rest_mass_density, Scalar{0.0}, electron_fraction)), + get(rest_mass_density) / + (get(cold_eos.pressure_from_density(rest_mass_density)) + + get(rest_mass_density) * + (1.0 + get(cold_eos.specific_internal_energy_from_density( + rest_mass_density)))) * + get(cold_eos.chi_from_density(rest_mass_density))); } { register_derived_classes_with_charm>();