Tidy + change EoS defn to P_over_rho

dinatraykova · Feb 19, 2024 · bd08c0b · bd08c0b
1 parent 56f5ecf
commit bd08c0b
Show file tree

Hide file tree

Showing 7 changed files with 55 additions and 71 deletions.
diff --git a/Examples/Fluid_Kerr/ConservedQuantities.hpp b/Examples/Fluid_Kerr/ConservedQuantities.hpp
@@ -15,7 +15,7 @@ namespace ConservedQuantities
 {
 // Primitive to conservative variables
 template <class data_t, template <typename> class vars_t>
-void PtoC(const data_t P_of_rho, vars_t<data_t> &vars)
+void PtoC(const data_t P_over_rho, vars_t<data_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<data_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
@@ -17,14 +17,10 @@ class DefaultEoS
 
     //! Set the pressure of the perfect fluid here to zero
     template <class data_t, template <typename> class vars_t>
-    void compute_eos(data_t &P_of_rho, data_t &dPdrho,
-                     const vars_t<data_t> &vars) const
+    void compute_eos(data_t &P_over_rho, const vars_t<data_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
@@ -25,14 +25,10 @@ class EoS
 
     //! Set the pressure of the perfect fluid here
     template <class data_t, template <typename> class vars_t>
-    void compute_eos(data_t &P_of_rho, data_t &dPdrho,
-                     const vars_t<data_t> &vars) const
+    void compute_eos(data_t &P_over_rho, const vars_t<data_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
@@ -14,7 +14,7 @@
 namespace Fluxes
 {
 template <class data_t, template <typename> class vars_t>
-vars_t<data_t> compute_flux(const data_t P_of_rho, const vars_t<data_t> &vars,
+vars_t<data_t> compute_flux(const data_t P_over_rho, const vars_t<data_t> &vars,
                             const int idir)
 {
     const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h);
@@ -38,32 +38,29 @@ vars_t<data_t> compute_flux(const data_t P_of_rho, const vars_t<data_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 data_t, template <typename> class vars_t>
-vars_t<data_t> compute_num_flux(const data_t P_of_rho,
+vars_t<data_t> compute_num_flux(const data_t P_over_rho,
                                 const vars_t<data_t> &vars, const int idir,
                                 const double lambda, const int sign)
 {
-    vars_t<data_t> out = compute_flux(P_of_rho, vars, idir);
+    vars_t<data_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
@@ -20,9 +20,8 @@ emtensor_t<data_t> PerfectFluid<eos_t>::compute_emtensor(
     emtensor_t<data_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<data_t> PerfectFluid<eos_t>::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<data_t> PerfectFluid<eos_t>::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<data_t> PerfectFluid<eos_t>::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<eos_t>::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<data_t> source = Sources::compute_source(P_of_rho, vars, d1);
+    vars_t<data_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<eos_t>::add_matter_rhs(
                                        GR_SPACEDIM / 2. * advec_chi);
     FOR(i)
     {
-        vars_t<data_t> flux = Fluxes::compute_flux(P_of_rho, vars, i);
+        vars_t<data_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<eos_t>::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<data_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<data_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<data_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<data_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<eos_t>::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<data_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<data_t> flux_left_p = Fluxes::compute_num_flux(
+            P_over_rho, vars_left_p, idir, m_lambda, -1);
 
         vars_t<data_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<data_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<data_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
@@ -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
@@ -15,7 +15,8 @@ namespace Sources
 {
 // Primitive to conservative variables
 template <class data_t, template <typename> class vars_t>
-vars_t<data_t> compute_source(const data_t P_of_rho, const vars_t<data_t> &vars,
+vars_t<data_t> compute_source(const data_t P_over_rho,
+                              const vars_t<data_t> &vars,
                               const vars_t<Tensor<1, data_t>> &d1)
 {
     data_t chi_regularised = simd_max(1e-6, vars.chi);
@@ -25,10 +26,8 @@ vars_t<data_t> compute_source(const data_t P_of_rho, const vars_t<data_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,24 +41,22 @@ vars_t<data_t> compute_source(const data_t P_of_rho, const vars_t<data_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]));
             }
         }
     }
     out.tau = 0.;
     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;
 }