From d2ba08128e8a44c393679e256870af11f85ad733 Mon Sep 17 00:00:00 2001 From: llibert94 Date: Tue, 1 Aug 2023 08:10:03 +0100 Subject: [PATCH] Bugs fixed in Initial Data and Primitive Recovery --- .../PerfectFluid/ConservativeRecovery.hpp | 4 +-- Examples/PerfectFluid/Fluxes.hpp | 2 +- Examples/PerfectFluid/InitialFluidData.hpp | 10 ++++-- Examples/PerfectFluid/PerfectFluid.impl.hpp | 32 +++++++++---------- Examples/PerfectFluid/PrimitiveRecovery.hpp | 4 +-- Examples/PerfectFluid/Sources.hpp | 2 +- Examples/PerfectFluid/WENODerivatives.hpp | 16 ++++++---- 7 files changed, 39 insertions(+), 31 deletions(-) diff --git a/Examples/PerfectFluid/ConservativeRecovery.hpp b/Examples/PerfectFluid/ConservativeRecovery.hpp index 441db72b2..59a604872 100644 --- a/Examples/PerfectFluid/ConservativeRecovery.hpp +++ b/Examples/PerfectFluid/ConservativeRecovery.hpp @@ -55,12 +55,12 @@ void CtoP(vars_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 index 10b739b7f..acb26d5c7 100644 --- a/Examples/PerfectFluid/Fluxes.hpp +++ b/Examples/PerfectFluid/Fluxes.hpp @@ -22,7 +22,7 @@ vars_t compute_flux(const vars_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 index f6b2d73ee..784dcaea3 100644 --- a/Examples/PerfectFluid/InitialFluidData.hpp +++ b/Examples/PerfectFluid/InitialFluidData.hpp @@ -51,7 +51,7 @@ class InitialFluidData // where am i? Coordinates coords(current_cell, m_dx, m_params.center); - MetricVars metric_vars; + const auto metric_vars = current_cell.template load_vars(); 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 index 9bdd1abba..9bce430c2 100644 --- a/Examples/PerfectFluid/PerfectFluid.impl.hpp +++ b/Examples/PerfectFluid/PerfectFluid.impl.hpp @@ -20,13 +20,13 @@ emtensor_t PerfectFluid::compute_emtensor( emtensor_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 PerfectFluid::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::add_matter_rhs( FOR(idir) { vars_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 flux_right_p = Fluxes::compute_num_flux(vars_right_p, idir, m_lambda); 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]; + // 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 flux_right_m = Fluxes::compute_num_flux(vars_right_m, idir, m_lambda); @@ -132,17 +132,17 @@ void PerfectFluid::add_matter_rhs( FOR(idir) { vars_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 flux_left_p = Fluxes::compute_num_flux(vars_left_p, idir, m_lambda); 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]; } + // 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 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 index 7dde2e5a4..b8da9e4ee 100644 --- 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 index 5171e91eb..ca3231cbc 100644 --- a/Examples/PerfectFluid/Sources.hpp +++ b/Examples/PerfectFluid/Sources.hpp @@ -62,4 +62,4 @@ vars_t compute_source(const vars_t &vars, } } // namespace Sources -#endif /* FLUXES_HPP_ */ +#endif /* SOURCES_HPP_ */ diff --git a/Examples/PerfectFluid/WENODerivatives.hpp b/Examples/PerfectFluid/WENODerivatives.hpp index aa852129e..ff958975c 100644 --- 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> 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( - current_cell.get_box_pointers().m_in_ptr[ivar], in_index, - strides[idir], dir_switch); - } - }); + FOR(idir) + { + var[idir] = get_Pface( + current_cell.get_box_pointers().m_in_ptr[ivar], + in_index, strides[idir], dir_switch); + } + }); return d1; }