From bd08c0b048bb7910c3846858c7b499f20ce380bc Mon Sep 17 00:00:00 2001 From: dinatraykova Date: Mon, 19 Feb 2024 03:55:14 +0100 Subject: [PATCH] Tidy + change EoS defn to P_over_rho --- Examples/Fluid_Kerr/ConservedQuantities.hpp | 6 +-- Examples/Fluid_Kerr/DefaultEoS.hpp | 8 +--- Examples/Fluid_Kerr/EoS.hpp | 8 +--- Examples/Fluid_Kerr/Fluxes.hpp | 15 +++--- Examples/Fluid_Kerr/PerfectFluid.impl.hpp | 53 ++++++++++----------- Examples/Fluid_Kerr/PrimitiveRecovery.hpp | 15 +++--- Examples/Fluid_Kerr/Sources.hpp | 21 ++++---- 7 files changed, 55 insertions(+), 71 deletions(-) diff --git a/Examples/Fluid_Kerr/ConservedQuantities.hpp b/Examples/Fluid_Kerr/ConservedQuantities.hpp index 47dfb2574..50df5dd8b 100644 --- a/Examples/Fluid_Kerr/ConservedQuantities.hpp +++ b/Examples/Fluid_Kerr/ConservedQuantities.hpp @@ -15,7 +15,7 @@ namespace ConservedQuantities { // Primitive to conservative variables template class vars_t> -void PtoC(const data_t P_of_rho, vars_t &vars) +void PtoC(const data_t P_over_rho, vars_t &vars) { data_t chi_regularised = simd_max(1e-6, vars.chi); Tensor<1, data_t> vi_D; @@ -29,10 +29,10 @@ void PtoC(const data_t P_of_rho, vars_t &vars) FOR(i) v2 += vars.vi[i] * vi_D[i]; data_t WW = 1. / (1. - v2); - data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + data_t hh = 1. + vars.eps + P_over_rho; vars.D = vars.rho * sqrt(WW); - vars.tau = vars.rho * hh * WW - P_of_rho - vars.D; + vars.tau = vars.rho * (hh * WW - P_over_rho) - vars.D; // S_j (note lower index) = - n^a T_ai FOR(i) { vars.Sj[i] = vars.rho * hh * WW * vi_D[i]; } diff --git a/Examples/Fluid_Kerr/DefaultEoS.hpp b/Examples/Fluid_Kerr/DefaultEoS.hpp index 950866f40..12505ffbe 100644 --- a/Examples/Fluid_Kerr/DefaultEoS.hpp +++ b/Examples/Fluid_Kerr/DefaultEoS.hpp @@ -17,14 +17,10 @@ class DefaultEoS //! Set the pressure of the perfect fluid here to zero template class vars_t> - void compute_eos(data_t &P_of_rho, data_t &dPdrho, - const vars_t &vars) const + void compute_eos(data_t &P_over_rho, const vars_t &vars) const { // The pressure value in function of rho - P_of_rho = 0.0; - - // The pressure gradient wrt rho0 - dPdrho = 0.0; + P_over_rho = 1. / 3.; } }; diff --git a/Examples/Fluid_Kerr/EoS.hpp b/Examples/Fluid_Kerr/EoS.hpp index f3d97957d..1769c44a5 100644 --- a/Examples/Fluid_Kerr/EoS.hpp +++ b/Examples/Fluid_Kerr/EoS.hpp @@ -25,14 +25,10 @@ class EoS //! Set the pressure of the perfect fluid here template class vars_t> - void compute_eos(data_t &P_of_rho, data_t &dPdrho, - const vars_t &vars) const + void compute_eos(data_t &P_over_rho, const vars_t &vars) const { // The pressure value in function of rho - P_of_rho = m_params.eos_w * vars.rho * (1. + vars.eps); - - // The pressure gradient wrt rho - dPdrho = m_params.eos_w * (1. + vars.eps); + P_over_rho = m_params.eos_w * (1. + vars.eps); } }; diff --git a/Examples/Fluid_Kerr/Fluxes.hpp b/Examples/Fluid_Kerr/Fluxes.hpp index c7df449d3..f38b80229 100644 --- a/Examples/Fluid_Kerr/Fluxes.hpp +++ b/Examples/Fluid_Kerr/Fluxes.hpp @@ -14,7 +14,7 @@ namespace Fluxes { template class vars_t> -vars_t compute_flux(const data_t P_of_rho, const vars_t &vars, +vars_t compute_flux(const data_t P_over_rho, const vars_t &vars, const int idir) { const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); @@ -38,32 +38,29 @@ vars_t compute_flux(const data_t P_of_rho, const vars_t &vars, data_t v2 = 0.; FOR(i) v2 += vars.vi[i] * vi_D[i]; - // data_t P = - // vars.rho * (1. + vars.eps) / 3.; // for now we assume a conformal fluid data_t WW = 1. / (1. - v2); - data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + data_t hh = 1. + vars.eps + P_over_rho; FOR(j) { out.Sj[j] = vars.lapse * vars.rho * hh * WW * vars.vi[idir] * vi_D[j] - vars.shift[idir] * vars.Sj[j]; FOR(k) - out.Sj[j] += vars.lapse * P_of_rho * h_UU[idir][k] * vars.h[j][k]; + out.Sj[j] += + vars.lapse * vars.rho * P_over_rho * h_UU[idir][k] * vars.h[j][k]; } out.tau = vars.lapse * (Sj_U[idir] - vars.D * vars.vi[idir]) - vars.shift[idir] * vars.tau; - // out.Jt = vars.nn * vars.vi[idir] * sqrt(WW); - return out; } template class vars_t> -vars_t compute_num_flux(const data_t P_of_rho, +vars_t compute_num_flux(const data_t P_over_rho, const vars_t &vars, const int idir, const double lambda, const int sign) { - vars_t out = compute_flux(P_of_rho, vars, idir); + vars_t out = compute_flux(P_over_rho, vars, idir); out.D += sign * lambda * vars.D; FOR(j) out.Sj[j] += sign * lambda * vars.Sj[j]; out.tau += sign * lambda * vars.tau; diff --git a/Examples/Fluid_Kerr/PerfectFluid.impl.hpp b/Examples/Fluid_Kerr/PerfectFluid.impl.hpp index cae192d12..4ea13ff4f 100644 --- a/Examples/Fluid_Kerr/PerfectFluid.impl.hpp +++ b/Examples/Fluid_Kerr/PerfectFluid.impl.hpp @@ -20,9 +20,8 @@ emtensor_t PerfectFluid::compute_emtensor( emtensor_t out; // Set up EoS - data_t P_of_rho = 0.; - data_t dPdrho = 0.; - my_eos.compute_eos(P_of_rho, dPdrho, vars); + data_t P_over_rho = 0.; + my_eos.compute_eos(P_over_rho, vars); // Useful quantities data_t chi_regularised = simd_max(1e-6, vars.chi); @@ -32,7 +31,7 @@ emtensor_t PerfectFluid::compute_emtensor( v2 += vars.h[i][j] * vars.vi[i] * vars.vi[j] / chi_regularised; } data_t WW = 1. / (1. - v2); - data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + data_t hh = 1. + vars.eps + P_over_rho; Tensor<1, data_t> vi_D; FOR(i) @@ -45,8 +44,9 @@ emtensor_t PerfectFluid::compute_emtensor( // S_ij = T_ij FOR(i, j) { - out.Sij[i][j] = vars.rho * hh * WW * vi_D[i] * vi_D[j] + - vars.h[i][j] * P_of_rho / chi_regularised; + out.Sij[i][j] = + vars.rho * (hh * WW * vi_D[i] * vi_D[j] + + vars.h[i][j] * P_over_rho / chi_regularised); } // S = Tr_S_ij @@ -56,7 +56,7 @@ emtensor_t PerfectFluid::compute_emtensor( FOR(i) { out.Si[i] = vars.rho * hh * WW * vi_D[i]; } // rho = n^a n^b T_ab - out.rho = vars.rho * hh * WW - P_of_rho; + out.rho = vars.rho * (hh * WW - P_over_rho); return out; } @@ -72,11 +72,10 @@ void PerfectFluid::add_matter_rhs( { using namespace TensorAlgebra; - data_t P_of_rho = 0.0; - data_t dPdrho = 0.0; - my_eos.compute_eos(P_of_rho, dPdrho, vars); + data_t P_over_rho = 0.0; + my_eos.compute_eos(P_over_rho, vars); - vars_t source = Sources::compute_source(P_of_rho, vars, d1); + vars_t source = Sources::compute_source(P_over_rho, vars, d1); // evolution equations for the fluid conservative variables data_t divshift = TensorAlgebra::compute_trace(d1.shift); data_t chi_regularised = simd_max(vars.chi, 1e-6); @@ -91,7 +90,7 @@ void PerfectFluid::add_matter_rhs( GR_SPACEDIM / 2. * advec_chi); FOR(i) { - vars_t flux = Fluxes::compute_flux(P_of_rho, vars, i); + vars_t flux = Fluxes::compute_flux(P_over_rho, vars, i); rhs.D -= GR_SPACEDIM / 2. * d1.chi[i] / chi_regularised * flux.D; FOR(j) rhs.Sj[j] -= @@ -105,19 +104,19 @@ void PerfectFluid::add_matter_rhs( vars_right_p.rho = rp.rho[idir]; vars_right_p.eps = rp.eps[idir]; FOR(j) { vars_right_p.vi[j] = rp.vi[j][idir]; } - my_eos.compute_eos(P_of_rho, dPdrho, vars_right_p); - ConservedQuantities::PtoC(P_of_rho, vars_right_p); + my_eos.compute_eos(P_over_rho, vars_right_p); + ConservedQuantities::PtoC(P_over_rho, vars_right_p); vars_t flux_right_p = Fluxes::compute_num_flux( - P_of_rho, vars_right_p, idir, m_lambda, -1); + P_over_rho, vars_right_p, idir, m_lambda, -1); vars_t vars_right_m = vars; vars_right_m.rho = rm.rho[idir]; vars_right_m.eps = rm.eps[idir]; FOR(j) vars_right_m.vi[j] = rm.vi[j][idir]; - my_eos.compute_eos(P_of_rho, dPdrho, vars_right_m); - ConservedQuantities::PtoC(P_of_rho, vars_right_m); - vars_t flux_right_m = - Fluxes::compute_num_flux(P_of_rho, vars_right_m, idir, m_lambda, 1); + my_eos.compute_eos(P_over_rho, vars_right_m); + ConservedQuantities::PtoC(P_over_rho, vars_right_m); + vars_t flux_right_m = Fluxes::compute_num_flux( + P_over_rho, vars_right_m, idir, m_lambda, 1); rhs.D += -1. / (2. * m_dx) * (flux_right_p.D + flux_right_m.D); FOR(j) @@ -134,19 +133,19 @@ void PerfectFluid::add_matter_rhs( vars_left_p.rho = lp.rho[idir]; vars_left_p.eps = lp.eps[idir]; FOR(j) { vars_left_p.vi[j] = lp.vi[j][idir]; } - my_eos.compute_eos(P_of_rho, dPdrho, vars_left_p); - ConservedQuantities::PtoC(P_of_rho, vars_left_p); - vars_t flux_left_p = - Fluxes::compute_num_flux(P_of_rho, vars_left_p, idir, m_lambda, -1); + my_eos.compute_eos(P_over_rho, vars_left_p); + ConservedQuantities::PtoC(P_over_rho, vars_left_p); + vars_t flux_left_p = Fluxes::compute_num_flux( + P_over_rho, vars_left_p, idir, m_lambda, -1); vars_t vars_left_m = vars; vars_left_m.rho = lm.rho[idir]; vars_left_m.eps = lm.eps[idir]; FOR(j) { vars_left_m.vi[j] = lm.vi[j][idir]; } - my_eos.compute_eos(P_of_rho, dPdrho, vars_left_m); - ConservedQuantities::PtoC(P_of_rho, vars_left_m); - vars_t flux_left_m = - Fluxes::compute_num_flux(P_of_rho, vars_left_m, idir, m_lambda, 1); + my_eos.compute_eos(P_over_rho, vars_left_m); + ConservedQuantities::PtoC(P_over_rho, vars_left_m); + vars_t flux_left_m = Fluxes::compute_num_flux( + P_over_rho, vars_left_m, idir, m_lambda, 1); rhs.D += 1. / (2. * m_dx) * (flux_left_p.D + flux_left_m.D); FOR(j) diff --git a/Examples/Fluid_Kerr/PrimitiveRecovery.hpp b/Examples/Fluid_Kerr/PrimitiveRecovery.hpp index 8f33d8c69..01110ee16 100644 --- a/Examples/Fluid_Kerr/PrimitiveRecovery.hpp +++ b/Examples/Fluid_Kerr/PrimitiveRecovery.hpp @@ -38,8 +38,7 @@ class PrimitiveRecovery data_t tolerance = 1e-8; data_t Wa2, Wa, xn, diff; - data_t P_of_rho = 0.; - data_t dPdrho = 0.; + data_t P_over_rho = 0.; data_t S2 = 0.; FOR(i, j) S2 += vars.chi * h_UU[i][j] * vars.Sj[i] * vars.Sj[j]; @@ -51,10 +50,10 @@ class PrimitiveRecovery vars.rho = vars.D / Wa; vars.eps = -1. + xa / Wa * (1. - Wa * Wa) + Wa * (1. + q); - // my_eos.compute_eos(P_of_rho, dPdrho, vars); - P_of_rho = vars.rho * (1. + vars.eps) / 3.; + // my_eos.compute_eos(P_over_rho, vars); + P_over_rho = (1. + vars.eps) / 3.; - xn = Wa * (1. + vars.eps + P_of_rho / vars.rho); + xn = Wa * (1. + vars.eps + P_over_rho); diff = abs(xn - xa); int i = 0; @@ -67,10 +66,10 @@ class PrimitiveRecovery vars.rho = vars.D / Wa; vars.eps = -1. + xa / Wa * (1. - Wa * Wa) + Wa * (1. + q); - // my_eos.compute_eos(P_of_rho, dPdrho, vars); - P_of_rho = vars.rho * (1. + vars.eps) / 3.; + // my_eos.compute_eos(P_over_rho, vars); + P_over_rho = (1. + vars.eps) / 3.; - xn = Wa * (1. + vars.eps + P_of_rho / vars.rho); + xn = Wa * (1. + vars.eps + P_over_rho); diff = abs(xn - xa); xa = xn; if (i >= 1000) diff --git a/Examples/Fluid_Kerr/Sources.hpp b/Examples/Fluid_Kerr/Sources.hpp index fa3478a07..17d7bc266 100644 --- a/Examples/Fluid_Kerr/Sources.hpp +++ b/Examples/Fluid_Kerr/Sources.hpp @@ -15,7 +15,8 @@ namespace Sources { // Primitive to conservative variables template class vars_t> -vars_t compute_source(const data_t P_of_rho, const vars_t &vars, +vars_t compute_source(const data_t P_over_rho, + const vars_t &vars, const vars_t> &d1) { data_t chi_regularised = simd_max(1e-6, vars.chi); @@ -25,10 +26,8 @@ vars_t compute_source(const data_t P_of_rho, const vars_t &vars, data_t v2 = 0.; FOR(i, j) v2 += vars.h[i][j] * vars.vi[j] * vars.vi[i] / chi_regularised; - // data_t P = - // vars.rho * (1. + vars.eps) / 3.; // for now we assume a conformal fluid data_t WW = 1. / (1. - v2); - data_t hh = 1. + vars.eps + P_of_rho / vars.rho; + data_t hh = 1. + vars.eps + P_over_rho; out.D = 0.; FOR(j) @@ -42,9 +41,9 @@ vars_t compute_source(const data_t P_of_rho, const vars_t &vars, out.Sj[j] += vars.lapse / 2. * (d1.h[i][k][j] - vars.h[i][k] * d1.chi[j] / chi_regularised) * - (vars.rho * hh * WW * vars.vi[i] * vars.vi[k] / - chi_regularised + - P_of_rho * h_UU[i][k]); + (vars.rho * (hh * WW * vars.vi[i] * vars.vi[k] / + chi_regularised + + P_over_rho * h_UU[i][k])); } } } @@ -52,14 +51,12 @@ vars_t compute_source(const data_t P_of_rho, const vars_t &vars, FOR(i, j) { out.tau += vars.lapse * (vars.A[i][j] + vars.h[i][j] / 3. * vars.K) * - (vars.rho * hh * WW * vars.vi[i] * vars.vi[j] / - chi_regularised + - P_of_rho * h_UU[i][j]) - + (vars.rho * + (hh * WW * vars.vi[i] * vars.vi[j] / chi_regularised + + P_over_rho * h_UU[i][j])) - vars.chi * h_UU[i][j] * vars.Sj[i] * d1.lapse[j]; } - // out.Jt = 0.; - return out; }