Merge pull request #24 from lwJi/tuning-eppm

Tuning eppm

AsterX/src/con2prim.cxx

@@ -17,7 +17,6 @@
 #include <eos_idealgas.hxx>
 
 #include "utils.hxx"
-#include <boost/math/tools/roots.hpp>
 
 namespace AsterX {
 using namespace std;
@@ -171,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;
 

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<CCTK_REAL, 2> rho_rc{reconstruct_pt(rho, p, true, false)};
+    const vec<CCTK_REAL, 2> rho_rc{reconstruct_pt(rho, p, true, true)};
     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, 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]),

ReconX/src/eppm.hxx

@@ -38,7 +38,8 @@ 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 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<const CCTK_REAL, 5> &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<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 bool &gf_is_press,
+     const array<const vect<int, dim>, 5> &cells, const bool &keep_var_positive,
      const GF3D2<const CCTK_REAL> &gf_press,
      const GF3D2<const CCTK_REAL> &gf_vel_dir,
      const reconstruct_params_t &reconstruct_params) {
@@ -139,8 +141,10 @@ eppm(const GF3D2<const CCTK_REAL> &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);