From 4720aa3ecb2e60cab3af03a5936fbe2d3f21bd4d Mon Sep 17 00:00:00 2001 From: Jooheon Yoo Date: Tue, 26 Sep 2023 10:42:19 -0700 Subject: [PATCH] Fix Kastaun for 1D EOS --- src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp | 4 +++- .../Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp index 5bdc7b1cce508..ebb0551c70321 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp @@ -249,7 +249,7 @@ Primitives FunctionOfMu::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::primitives(const double mu) const { if constexpr (ThermodynamicDim == 1) { p_hat = get(equation_of_state_.pressure_from_density(Scalar(rho_hat))); + epsilon_hat = get(equation_of_state_.specific_internal_energy_from_density( + Scalar(rho_hat))); } else if constexpr (ThermodynamicDim == 2) { p_hat = get(equation_of_state_.pressure_from_density_and_energy( Scalar(rho_hat), Scalar(epsilon_hat))); diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp index f49c3e41b88b9..163ca9ce34452 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp @@ -120,7 +120,7 @@ Primitives FunctionOfZ::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::primitives( if constexpr (ThermodynamicDim == 1) { p_hat = get(equation_of_state_.pressure_from_density(Scalar(rho_hat))); + epsilon_hat = get(equation_of_state_.specific_internal_energy_from_density( + Scalar(rho_hat))); } else if constexpr (ThermodynamicDim == 2) { p_hat = get(equation_of_state_.pressure_from_density_and_energy( Scalar(rho_hat), Scalar(epsilon_hat)));