Bugs fixed in Initial Data and Primitive Recovery

dinatraykova · Aug 1, 2023 · d2ba081 · d2ba081
1 parent c74e725
commit d2ba081
Show file tree

Hide file tree

Showing 7 changed files with 39 additions and 31 deletions.
diff --git a/Examples/PerfectFluid/ConservativeRecovery.hpp b/Examples/PerfectFluid/ConservativeRecovery.hpp
@@ -55,12 +55,12 @@ void CtoP(vars_t<data_t> &vars)
 
     // eps
     vars.eps =
-        sqrt2 / (sqrt(2.) * vars.D) - 1.; // Is it dangerous to divide by D?
+        sqrt2 / (2. * vars.D) - 1.; // Is it dangerous to divide by D?
     // rho
     vars.rho = (sqrt1 - E) / (1. + vars.eps);
     // vi_D
     Tensor<1, data_t> vi_D;
-    FOR(i) vi_D[i] = 3. / 2. * vars.Sj[i] * (sqrt1 - E) / (sqrt2 * sqrt2);
+    FOR(i) vi_D[i] = 3. * vars.Sj[i] * (sqrt1 - E) / (sqrt2 * sqrt2);
     // vi
     FOR(i)
     {

diff --git a/Examples/PerfectFluid/Fluxes.hpp b/Examples/PerfectFluid/Fluxes.hpp
@@ -22,7 +22,7 @@ vars_t<data_t> compute_flux(const vars_t<data_t> &vars, const int idir)
     Tensor<1, data_t> Sj_U;
     FOR(i)
     {
-        Sj_U[i] = 0;
+        Sj_U[i] = 0.;
         FOR(j) Sj_U[i] += vars.chi * h_UU[i][j] * vars.Sj[j];
     }
 

diff --git a/Examples/PerfectFluid/InitialFluidData.hpp b/Examples/PerfectFluid/InitialFluidData.hpp
@@ -51,7 +51,7 @@ class InitialFluidData
         // where am i?
         Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
 
-        MetricVars<data_t> metric_vars;
+        const auto metric_vars = current_cell.template load_vars<MetricVars>();
 
         data_t x = coords.x;
         double y = coords.y;
@@ -82,7 +82,13 @@ class InitialFluidData
 
         data_t D = rho * sqrt(WW);
         data_t tau = rho * hh * WW - P - D;
-        FOR(i) Sj[i] = rho * hh * WW * vi[i];
+        FOR(i)
+        {
+            Sj[i] = 0.;
+            FOR(j)
+            Sj[i] +=
+                rho * hh * WW * metric_vars.h[i][j] * vi[j] / chi_regularised;
+        }
 
         // store the vars
         current_cell.store_vars(rho, c_rho);

diff --git a/Examples/PerfectFluid/PerfectFluid.impl.hpp b/Examples/PerfectFluid/PerfectFluid.impl.hpp
@@ -20,13 +20,13 @@ emtensor_t<data_t> PerfectFluid<eos_t>::compute_emtensor(
     emtensor_t<data_t> out;
 
     // Set up EoS
-    data_t P_of_rho = 0.0;
-    data_t dPdrho = 0.0;
+    data_t P_of_rho = 0.;
+    data_t dPdrho = 0.;
     my_eos.compute_eos(P_of_rho, dPdrho, vars);
 
     // Useful quantities
     data_t chi_regularised = simd_max(1e-6, vars.chi);
-    data_t v2 = 0.0;
+    data_t v2 = 0.;
     FOR(i, j)
     {
         v2 += vars.h[i][j] * vars.vi[i] * vars.vi[j] / chi_regularised;
@@ -37,7 +37,7 @@ emtensor_t<data_t> PerfectFluid<eos_t>::compute_emtensor(
     Tensor<1, data_t> vi_D;
     FOR(i)
     {
-        vi_D[i] = 0;
+        vi_D[i] = 0.;
         FOR(j) { vi_D[i] += vars.h[i][j] * vars.vi[j] / chi_regularised; }
     }
 
@@ -105,17 +105,17 @@ void PerfectFluid<eos_t>::add_matter_rhs(
     FOR(idir)
     {
         vars_t<data_t> vars_right_p = vars;
-        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]; }
+        // 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]; }
         ConservativeRecovery::PtoC(vars_right_p);
         vars_t<data_t> flux_right_p =
             Fluxes::compute_num_flux(vars_right_p, idir, m_lambda);
 
         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];
+        // 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];
         ConservativeRecovery::PtoC(vars_right_m);
         vars_t<data_t> flux_right_m =
             Fluxes::compute_num_flux(vars_right_m, idir, m_lambda);
@@ -132,17 +132,17 @@ void PerfectFluid<eos_t>::add_matter_rhs(
     FOR(idir)
     {
         vars_t<data_t> vars_left_p = vars;
-        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]; }
+        // 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]; }
         ConservativeRecovery::PtoC(vars_left_p);
         vars_t<data_t> flux_left_p =
             Fluxes::compute_num_flux(vars_left_p, idir, m_lambda);
 
         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]; }
+        // 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]; }
         ConservativeRecovery::PtoC(vars_left_m);
         vars_t<data_t> flux_left_m =
             Fluxes::compute_num_flux(vars_left_m, idir, m_lambda);

diff --git a/Examples/PerfectFluid/PrimitiveRecovery.hpp b/Examples/PerfectFluid/PrimitiveRecovery.hpp
@@ -43,12 +43,12 @@ class PrimitiveRecovery
 
         // eps
         vars.eps =
-            sqrt2 / (sqrt(2.) * vars.D) - 1.; // Is it dangerous to divide by D?
+            sqrt2 / (2. * vars.D) - 1.; // Is it dangerous to divide by D?
         // rho
         vars.rho = (sqrt1 - E) / (1. + vars.eps);
         // vi_D
         Tensor<1, data_t> vi_D;
-        FOR(i) vi_D[i] = 3. / 2. * vars.Sj[i] * (sqrt1 - E) / (sqrt2 * sqrt2);
+        FOR(i) vi_D[i] = 3. * vars.Sj[i] * (sqrt1 - E) / (sqrt2 * sqrt2);
         // vi
         FOR(i)
         {

diff --git a/Examples/PerfectFluid/Sources.hpp b/Examples/PerfectFluid/Sources.hpp
@@ -62,4 +62,4 @@ vars_t<data_t> compute_source(const vars_t<data_t> &vars,
 }
 
 } // namespace Sources
-#endif /* FLUXES_HPP_ */
+#endif /* SOURCES_HPP_ */
diff --git a/Examples/PerfectFluid/WENODerivatives.hpp b/Examples/PerfectFluid/WENODerivatives.hpp
@@ -135,14 +135,16 @@ class WENODerivatives
         const auto in_index = current_cell.get_in_index();
         const auto strides = current_cell.get_box_pointers().m_in_stride;
         vars_t<Tensor<1, data_t>> d1;
-        d1.enum_mapping([&](const int &ivar, Tensor<1, data_t> &var) {
-            FOR(idir)
+        d1.enum_mapping(
+            [&](const int &ivar, Tensor<1, data_t> &var)
             {
-                var[idir] = get_Pface<data_t>(
-                    current_cell.get_box_pointers().m_in_ptr[ivar], in_index,
-                    strides[idir], dir_switch);
-            }
-        });
+                FOR(idir)
+                {
+                    var[idir] = get_Pface<data_t>(
+                        current_cell.get_box_pointers().m_in_ptr[ivar],
+                        in_index, strides[idir], dir_switch);
+                }
+            });
         return d1;
     }