Fix Kastaun for 1D EOS

jyoo1042 · Sep 26, 2023 · 4720aa3 · 4720aa3
1 parent fb1f216
commit 4720aa3
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 2 deletions.
diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp
@@ -249,7 +249,7 @@ Primitives FunctionOfMu<ThermodynamicDim>::primitives(const double mu) const {
       q_ - 0.5 * b_squared_ -
       0.5 * square(mu * x) * (r_squared_ * b_squared_ - r_dot_b_squared_);
   // Equation (42) with bounds from Equation (6)
-  const double epsilon_hat = std::clamp(
+  double epsilon_hat = std::clamp(
       w_hat * (q_bar - mu * r_bar_squared) +
           v_hat_squared * square(w_hat) / (1.0 + w_hat),
       equation_of_state_.specific_internal_energy_lower_bound(rho_hat),
@@ -259,6 +259,8 @@ Primitives FunctionOfMu<ThermodynamicDim>::primitives(const double mu) const {
   if constexpr (ThermodynamicDim == 1) {
     p_hat =
         get(equation_of_state_.pressure_from_density(Scalar<double>(rho_hat)));
+    epsilon_hat = get(equation_of_state_.specific_internal_energy_from_density(
+        Scalar<double>(rho_hat)));
   } else if constexpr (ThermodynamicDim == 2) {
     p_hat = get(equation_of_state_.pressure_from_density_and_energy(
         Scalar<double>(rho_hat), Scalar<double>(epsilon_hat)));

diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp
@@ -120,7 +120,7 @@ Primitives FunctionOfZ<ThermodynamicDim, EnforcePhysicality>::primitives(
                  equation_of_state_.rest_mass_density_upper_bound());
 
   // Equation (C14) and (C16)
-  const double epsilon_hat = std::clamp(
+  double epsilon_hat = std::clamp(
       w_hat * q_ - z * (r_ - z / (1. + w_hat)),
       equation_of_state_.specific_internal_energy_lower_bound(rho_hat),
       equation_of_state_.specific_internal_energy_upper_bound(rho_hat));
@@ -129,6 +129,8 @@ Primitives FunctionOfZ<ThermodynamicDim, EnforcePhysicality>::primitives(
   if constexpr (ThermodynamicDim == 1) {
     p_hat =
         get(equation_of_state_.pressure_from_density(Scalar<double>(rho_hat)));
+    epsilon_hat = get(equation_of_state_.specific_internal_energy_from_density(
+        Scalar<double>(rho_hat)));
   } else if constexpr (ThermodynamicDim == 2) {
     p_hat = get(equation_of_state_.pressure_from_density_and_energy(
         Scalar<double>(rho_hat), Scalar<double>(epsilon_hat)));