From e78170f67de078c527e09cfaf64cc3363eec98ac Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Wed, 12 Jul 2023 16:56:18 +0000 Subject: [PATCH 1/3] ReconX: Keep rho positive also in limit of eppm --- AsterX/src/fluxes.cxx | 9 +++++++-- ReconX/src/eppm.hxx | 14 +++++++++----- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/AsterX/src/fluxes.cxx b/AsterX/src/fluxes.cxx index d42f1236..9277c80e 100644 --- a/AsterX/src/fluxes.cxx +++ b/AsterX/src/fluxes.cxx @@ -164,7 +164,7 @@ 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, false)}; + const vec rho_rc{reconstruct_pt(rho, p, true, true)}; const vec, 3> vels_rc([&](int i) ARITH_INLINE { return vec{reconstruct_pt(gf_vels(i), p, false, false)}; }); @@ -359,7 +359,8 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { isnan(fluxmomxs(dir)(p.I)) || isnan(fluxmomys(dir)(p.I)) || isnan(fluxmomzs(dir)(p.I)) || isnan(fluxtaus(dir)(p.I)) || isnan(fluxBxs(dir)(p.I)) || isnan(fluxBys(dir)(p.I)) || - isnan(fluxBzs(dir)(p.I))) { + isnan(fluxBzs(dir)(p.I)) || rho_rc(0) < 0.0 || rho_rc(1) < 0.0 || + press_rc(0) < 0.0 || press_rc(1) < 0.0) { printf("cctk_iteration = %i, dir = %i, ijk = %i, %i, %i, " "x, y, z = %16.8e, %16.8e, %16.8e.\n", cctk_iteration, dir, p.i, p.j, p.k, p.x, p.y, p.z); @@ -399,6 +400,10 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { printf(" bsq_rc = %16.8e, %16.8e \n", bsq_rc(0), bsq_rc(1)); printf(" press_rc = %16.8e, %16.8e \n", press_rc(0), press_rc(1)); printf(" eps_rc = %16.8e, %16.8e \n", eps_rc(0), eps_rc(1)); + printf(" rho = %16.8e, %16.8e, %16.8e, %16.8e, %16.8e, %16.8e;\n", + rho(p.I - p.DI[dir] * 3), rho(p.I - p.DI[dir] * 2), + rho(p.I - p.DI[dir]), rho(p.I), rho(p.I + p.DI[dir]), + rho(p.I + p.DI[dir] * 2)); printf(" press = %16.8e, %16.8e, %16.8e, %16.8e, %16.8e, %16.8e;\n", press(p.I - p.DI[dir] * 3), press(p.I - p.DI[dir] * 2), press(p.I - p.DI[dir]), press(p.I), press(p.I + p.DI[dir]), diff --git a/ReconX/src/eppm.hxx b/ReconX/src/eppm.hxx index fb2a3526..0897520c 100644 --- a/ReconX/src/eppm.hxx +++ b/ReconX/src/eppm.hxx @@ -38,7 +38,8 @@ 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 bool &gf_is_press) { + const CCTK_INT interface, const CCTK_REAL C, + const bool &keep_var_positive) { const CCTK_INT Im = interface; const CCTK_INT I = interface + 1; const CCTK_INT Ip = interface + 2; @@ -54,7 +55,8 @@ limit(const array &gf, const CCTK_REAL gf_rc, MIN3(C * fabs(D2aL), C * fabs(D2aR), fabs(D2a)) * 1.0 / 3.0; const CCTK_REAL gf_avg = 0.5 * (gf[I] + gf[Ip]); - if (D2a * D2aR < 0 || D2a * D2aL < 0 || (gf_is_press && D2aLim > gf_avg)) + if (D2a * D2aR < 0 || D2a * D2aL < 0 || + (keep_var_positive && D2aLim > gf_avg)) return gf_avg; else return gf_avg - D2aLim; @@ -111,7 +113,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 bool &gf_is_press, + const array, 5> &cells, const bool &keep_var_positive, const GF3D2 &gf_press, const GF3D2 &gf_vel_dir, const reconstruct_params_t &reconstruct_params) { @@ -139,8 +141,10 @@ eppm(const GF3D2 &gf_var, CCTK_REAL rc_plus = approx_at_cell_interface(gf_stencil, iplus); /* limit */ - 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); + rc_minus = + limit(gf_stencil, rc_minus, iminus, enhanced_ppm_C2, keep_var_positive); + rc_plus = + limit(gf_stencil, rc_plus, iplus, enhanced_ppm_C2, keep_var_positive); /* monotonize */ monotonize(gf_stencil, rc_minus, rc_plus, enhanced_ppm_C2); From 12c961dbb784ace473b1269dda2d205a1410a9f8 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Wed, 12 Jul 2023 18:54:02 +0000 Subject: [PATCH 2/3] AsterX: Use 1D Palenzuela if 2D Noble spit negative press --- AsterX/src/con2prim.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx index 88e17958..7595313f 100644 --- a/AsterX/src/con2prim.cxx +++ b/AsterX/src/con2prim.cxx @@ -17,7 +17,6 @@ #include #include "utils.hxx" -#include namespace AsterX { using namespace std; @@ -95,7 +94,7 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, CCTK_INT c2p_succeeded_Pal = 0; // false for now c2p_Noble.solve(eos_th, pv, pv_seeds, cv, glo, c2p_succeeded_Noble); - if (!c2p_succeeded_Noble) { + if (!c2p_succeeded_Noble || pv.press < 0.0) { c2p_Pal.solve(eos_th, pv, pv_seeds, cv, glo, c2p_succeeded_Pal); } From 50b5485130f885ebc2af448eebebe2f22472a51b Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Wed, 12 Jul 2023 15:31:32 -0500 Subject: [PATCH 3/3] AsterX: Set to atmo if press is less then p_atm in con2prim --- AsterX/src/con2prim.cxx | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx index 7595313f..28c14bf0 100644 --- a/AsterX/src/con2prim.cxx +++ b/AsterX/src/con2prim.cxx @@ -94,7 +94,7 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, CCTK_INT c2p_succeeded_Pal = 0; // false for now c2p_Noble.solve(eos_th, pv, pv_seeds, cv, glo, c2p_succeeded_Noble); - if (!c2p_succeeded_Noble || pv.press < 0.0) { + if (!c2p_succeeded_Noble) { c2p_Pal.solve(eos_th, pv, pv_seeds, cv, glo, c2p_succeeded_Pal); } @@ -170,6 +170,16 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, } } */ + + if(pv.press < p_atm) { + // set to atmo + cv.dBvec(0) = dBx(p.I); + cv.dBvec(1) = dBy(p.I); + cv.dBvec(2) = dBz(p.I); + pv.Bvec = cv.dBvec / sqrt_detg; + atmo.set(pv, cv, glo); + } + // dummy vars CCTK_REAL Ex, Ey, Ez;