Skip to content

Commit

Permalink
Merge pull request #18 from lwJi/fix-eppm-rec-press
Browse files Browse the repository at this point in the history 
ReconX: Fix limit in eppm to keep press positive
  • Loading branch information
@lwJi
lwJi authored Jul 3, 2023
2 parents f6e5cfe + 72242cc commit b677386
Show file tree
Hide file tree
Showing 3 changed files with 28 additions and 22 deletions.
16 changes: 9 additions & 7 deletions AsterX/src/fluxes.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,9 +123,10 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {

const auto reconstruct_pt =
[=] CCTK_DEVICE(const GF3D2<const CCTK_REAL> &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<vec<CCTK_REAL, 4>, 2> lam, vec<CCTK_REAL, 2> var,
Expand Down Expand Up @@ -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<CCTK_REAL, 2> rho_rc{reconstruct_pt(rho, p, true)};
const vec<CCTK_REAL, 2> rho_rc{reconstruct_pt(rho, p, true, false)};
const vec<vec<CCTK_REAL, 2>, 3> vels_rc([&](int i) ARITH_INLINE {
return vec<CCTK_REAL, 2>{reconstruct_pt(gf_vels(i), p, false)};
return vec<CCTK_REAL, 2>{reconstruct_pt(gf_vels(i), p, false, false)};
});
const vec<CCTK_REAL, 2> press_rc{reconstruct_pt(press, p, false)};
const vec<CCTK_REAL, 2> press_rc{reconstruct_pt(press, p, false, true)};
const vec<vec<CCTK_REAL, 2>, 3> Bs_rc([&](int i) ARITH_INLINE {
return vec<CCTK_REAL, 2>{reconstruct_pt(gf_Bvecs(i), p, false)};
return vec<CCTK_REAL, 2>{
reconstruct_pt(gf_Bvecs(i), p, false, false)};
});

/* Interpolate metric components from vertices to faces */
Expand Down
25 changes: 13 additions & 12 deletions ReconX/src/eppm.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<const CCTK_REAL, 5> &gf,
const CCTK_INT interface) {
Expand All @@ -40,7 +38,7 @@ approx_at_cell_interface(const array<const CCTK_REAL, 5> &gf,

inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST CCTK_REAL
limit(const array<const CCTK_REAL, 5> &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;
Expand All @@ -55,10 +53,11 @@ limit(const array<const CCTK_REAL, 5> &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;
}
}

Expand Down Expand Up @@ -112,7 +111,7 @@ monotonize(const array<const CCTK_REAL, 5> &gf, CCTK_REAL &rc_minus,
& Colella (2011) */
inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array<CCTK_REAL, 2>
eppm(const GF3D2<const CCTK_REAL> &gf_var,
const array<const vect<int, dim>, 5> &cells,
const array<const vect<int, dim>, 5> &cells, const bool &gf_is_press,
const GF3D2<const CCTK_REAL> &gf_press,
const GF3D2<const CCTK_REAL> &gf_vel_dir,
const reconstruct_params_t &reconstruct_params) {
Expand All @@ -122,7 +121,7 @@ eppm(const GF3D2<const CCTK_REAL> &gf_var,
const auto &I = cells[2];
const auto &Ip = cells[3];
const auto &Ipp = cells[4];

const array<const CCTK_REAL, 5> gf_stencil{gf_var(Imm), gf_var(Im), gf_var(I),
gf_var(Ip), gf_var(Ipp)};

Expand All @@ -133,13 +132,15 @@ eppm(const GF3D2<const CCTK_REAL> &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);
Expand Down
9 changes: 6 additions & 3 deletions ReconX/src/reconstruct.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,8 @@ enum class reconstruction_t {
inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_HOST CCTK_DEVICE array<CCTK_REAL, 2>
reconstruct(const GF3D2<const CCTK_REAL> &gf_var, const PointDesc &p,
const reconstruction_t &reconstruction, const int &dir,
const bool &gf_is_rho, const GF3D2<const CCTK_REAL> &gf_press,
const bool &gf_is_rho, const bool &gf_is_press,
const GF3D2<const CCTK_REAL> &gf_press,
const GF3D2<const CCTK_REAL> &gf_vel_dir,
const reconstruct_params_t &reconstruct_params) {
// Neighbouring "plus" and "minus" cell indices
Expand Down Expand Up @@ -99,9 +100,11 @@ reconstruct(const GF3D2<const CCTK_REAL> &gf_var, const PointDesc &p,
const array<const vect<int, dim>, 5> cells_Ip = {Imm, Im, Ip, Ipp, Ippp};

const array<CCTK_REAL, 2> 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<CCTK_REAL, 2> 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<CCTK_REAL, 2>{rc_Im[1], rc_Ip[0]};
}
Expand Down

0 comments on commit b677386

Please sign in to comment.