Fix a bug in Equilibrium3D EoS

src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp

@@ -641,12 +641,47 @@ class EquationOfState<IsRelativistic, 3> : 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

src/PointwiseFunctions/Hydro/EquationsOfState/Equilibrium3D.cpp

@@ -119,14 +119,12 @@ Equilibrium3D<EquilEos>::sound_speed_squared_from_density_and_temperature_impl(
   // Cold part plus temperature dependent part, see the 3D EoS documentation for
   // expression.
   return Scalar<DataType>{
-      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<HybridEos<PolytropicFluid<true>>>;

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<true> cold_eos{100.0, 2.0};
-  const EoS::HybridEos<EoS::PolytropicFluid<true>> underlying_eos{cold_eos,
-                                                                  2.0};
+
+  const double ideal_adiabatic_index = 2.0;
+  const double cold_polytropic_index = 2.0;
+  const EoS::PolytropicFluid<true> cold_eos{100.0, cold_polytropic_index};
+  const EoS::IdealFluid<true> eos_ideal_fluid{ideal_adiabatic_index};
+
+  EoS::Equilibrium3D<EoS::IdealFluid<true>> eos_ideal_fluid_3d{eos_ideal_fluid};
+  const EoS::HybridEos<EoS::PolytropicFluid<true>> underlying_eos{
+      cold_eos, ideal_adiabatic_index};
   EoS::Equilibrium3D<EoS::HybridEos<EoS::PolytropicFluid<true>>> 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<DataVector> rest_mass_density{
+        DataVector{1e-10, 1e-6, 1e-4, 1e-3, 1e-2, 1.0}};
+    // electron fraction should be irrelevant
+    const Scalar<DataVector> electron_fraction{
+        DataVector{1e-10, 1e-6, 1e-4, 1e-3, 1e-2, 1.0}};
+    const Scalar<DataVector> 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<DataVector> pressure =
+        eos_ideal_fluid_3d.pressure_from_density_and_temperature(
+            rest_mass_density, temperature, electron_fraction);
+    const Scalar<DataVector> 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<double> 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<double> rest_mass_density{{1e-3}};
+    // temperature and electron fraction should be irrelevant
+    const Scalar<double> electron_fraction{{1e-1}};
+    const Scalar<double> temperature{{1e-1}};
+    const Scalar<double> pressure =
+        eos_ideal_fluid_3d.pressure_from_density_and_temperature(
+            rest_mass_density, temperature, electron_fraction);
+    const Scalar<double> 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<double>{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<double>{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<EoS::EquationOfState<true, 3>>();