diff --git a/AsterX/src/fluxes.cxx b/AsterX/src/fluxes.cxx index e3e823cb..d5c8096c 100644 --- a/AsterX/src/fluxes.cxx +++ b/AsterX/src/fluxes.cxx @@ -123,9 +123,10 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { const auto reconstruct_pt = [=] CCTK_DEVICE(const GF3D2 &var, const PointDesc &p, - const bool &gf_is_rho) CCTK_ATTRIBUTE_ALWAYS_INLINE { - return reconstruct(var, p, reconstruction, dir, gf_is_rho, press, - gf_vels(dir), reconstruct_params); + const bool &gf_is_rho, + const bool &gf_is_press) CCTK_ATTRIBUTE_ALWAYS_INLINE { + return reconstruct(var, p, reconstruction, dir, gf_is_rho, gf_is_press, + press, gf_vels(dir), reconstruct_params); }; const auto calcflux = [=] CCTK_DEVICE(vec, 2> lam, vec var, @@ -160,13 +161,14 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { /* Reconstruct primitives from the cells on left (indice 0) and right * (indice 1) side of this face rc = reconstructed variables or * computed from reconstructed variables */ - const vec rho_rc{reconstruct_pt(rho, p, true)}; + const vec rho_rc{reconstruct_pt(rho, p, true, false)}; const vec, 3> vels_rc([&](int i) ARITH_INLINE { - return vec{reconstruct_pt(gf_vels(i), p, false)}; + return vec{reconstruct_pt(gf_vels(i), p, false, false)}; }); - const vec press_rc{reconstruct_pt(press, p, false)}; + const vec press_rc{reconstruct_pt(press, p, false, true)}; const vec, 3> Bs_rc([&](int i) ARITH_INLINE { - return vec{reconstruct_pt(gf_Bvecs(i), p, false)}; + return vec{ + reconstruct_pt(gf_Bvecs(i), p, false, false)}; }); /* Interpolate metric components from vertices to faces */ diff --git a/ReconX/src/eppm.hxx b/ReconX/src/eppm.hxx index 21fbb176..fb2a3526 100644 --- a/ReconX/src/eppm.hxx +++ b/ReconX/src/eppm.hxx @@ -25,8 +25,6 @@ using namespace std; using namespace Arith; using namespace Loop; -enum InterfaceType { MINUS = 0, PLUS = 1 }; - inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST CCTK_REAL approx_at_cell_interface(const array &gf, const CCTK_INT interface) { @@ -40,7 +38,7 @@ approx_at_cell_interface(const array &gf, inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST CCTK_REAL limit(const array &gf, const CCTK_REAL gf_rc, - const CCTK_INT interface, const CCTK_REAL C) { + const CCTK_INT interface, const CCTK_REAL C, const bool &gf_is_press) { const CCTK_INT Im = interface; const CCTK_INT I = interface + 1; const CCTK_INT Ip = interface + 2; @@ -55,10 +53,11 @@ limit(const array &gf, const CCTK_REAL gf_rc, const CCTK_REAL D2aLim = copysign(1.0, D2a) * MIN3(C * fabs(D2aL), C * fabs(D2aR), fabs(D2a)) * 1.0 / 3.0; - if (D2a * D2aR >= 0 && D2a * D2aL >= 0) - return 0.5 * (gf[I] + gf[Ip]) - D2aLim; + const CCTK_REAL gf_avg = 0.5 * (gf[I] + gf[Ip]); + if (D2a * D2aR < 0 || D2a * D2aL < 0 || (gf_is_press && D2aLim > gf_avg)) + return gf_avg; else - return 0.5 * (gf[I] + gf[Ip]); + return gf_avg - D2aLim; } } @@ -112,7 +111,7 @@ monotonize(const array &gf, CCTK_REAL &rc_minus, & Colella (2011) */ inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array eppm(const GF3D2 &gf_var, - const array, 5> &cells, + const array, 5> &cells, const bool &gf_is_press, const GF3D2 &gf_press, const GF3D2 &gf_vel_dir, const reconstruct_params_t &reconstruct_params) { @@ -122,7 +121,7 @@ eppm(const GF3D2 &gf_var, const auto &I = cells[2]; const auto &Ip = cells[3]; const auto &Ipp = cells[4]; - + const array gf_stencil{gf_var(Imm), gf_var(Im), gf_var(I), gf_var(Ip), gf_var(Ipp)}; @@ -133,13 +132,15 @@ eppm(const GF3D2 &gf_var, const CCTK_REAL &ppm_omega2 = reconstruct_params.ppm_omega2; const CCTK_REAL &enhanced_ppm_C2 = reconstruct_params.enhanced_ppm_C2; + const int iminus = 0, iplus = 1; + /* approx at cell interface */ - CCTK_REAL rc_minus = approx_at_cell_interface(gf_stencil, MINUS); - CCTK_REAL rc_plus = approx_at_cell_interface(gf_stencil, PLUS); + CCTK_REAL rc_minus = approx_at_cell_interface(gf_stencil, iminus); + CCTK_REAL rc_plus = approx_at_cell_interface(gf_stencil, iplus); /* limit */ - rc_minus = limit(gf_stencil, rc_minus, MINUS, enhanced_ppm_C2); - rc_plus = limit(gf_stencil, rc_plus, PLUS, enhanced_ppm_C2); + rc_minus = limit(gf_stencil, rc_minus, iminus, enhanced_ppm_C2, gf_is_press); + rc_plus = limit(gf_stencil, rc_plus, iplus, enhanced_ppm_C2, gf_is_press); /* monotonize */ monotonize(gf_stencil, rc_minus, rc_plus, enhanced_ppm_C2); diff --git a/ReconX/src/reconstruct.hxx b/ReconX/src/reconstruct.hxx index ea42d7b9..937140c2 100644 --- a/ReconX/src/reconstruct.hxx +++ b/ReconX/src/reconstruct.hxx @@ -34,7 +34,8 @@ enum class reconstruction_t { inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_HOST CCTK_DEVICE array reconstruct(const GF3D2 &gf_var, const PointDesc &p, const reconstruction_t &reconstruction, const int &dir, - const bool &gf_is_rho, const GF3D2 &gf_press, + const bool &gf_is_rho, const bool &gf_is_press, + const GF3D2 &gf_press, const GF3D2 &gf_vel_dir, const reconstruct_params_t &reconstruct_params) { // Neighbouring "plus" and "minus" cell indices @@ -99,9 +100,11 @@ reconstruct(const GF3D2 &gf_var, const PointDesc &p, const array, 5> cells_Ip = {Imm, Im, Ip, Ipp, Ippp}; const array rc_Im = - eppm(gf_var, cells_Im, gf_press, gf_vel_dir, reconstruct_params); + eppm(gf_var, cells_Im, gf_is_press, gf_press, gf_vel_dir, + reconstruct_params); const array rc_Ip = - eppm(gf_var, cells_Ip, gf_press, gf_vel_dir, reconstruct_params); + eppm(gf_var, cells_Ip, gf_is_press, gf_press, gf_vel_dir, + reconstruct_params); return array{rc_Im[1], rc_Ip[0]}; }