RePrimAnd Con2Prim

AsterX/interface.ccl

@@ -7,7 +7,7 @@ INHERITS: CarpetXRegrid ADMBaseX HydroBaseX TmunuBaseX AsterMasks
 USES INCLUDE HEADER: fixmath.hxx
 USES INCLUDE HEADER: loop_device.hxx
 USES INCLUDE HEADER: setup_eos.hxx
-USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DPalenzuela.hxx c2p_1DEntropy.hxx
+USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DRePrimAnd.hxx c2p_1DPalenzuela.hxx c2p_1DEntropy.hxx
 USES INCLUDE HEADER: reconstruct.hxx
 USES INCLUDE HEADER: aster_fd.hxx aster_interp.hxx aster_utils.hxx
 

AsterX/param.ccl

@@ -174,6 +174,8 @@ USES CCTK_REAL vw_lim
 USES BOOLEAN Ye_lenient
 USES CCTK_INT max_iter
 USES CCTK_REAL c2p_tol
+USES BOOLEAN soft_root_convergence
+USES CCTK_REAL soft_root_width_factor
 USES BOOLEAN use_press_atmo
 USES CCTK_REAL r_atmo
 USES CCTK_REAL n_rho_atmo

AsterX/src/con2prim.cxx

@@ -8,6 +8,7 @@
 #include "c2p.hxx"
 #include "c2p_1DEntropy.hxx"
 #include "c2p_1DPalenzuela.hxx"
+#include "c2p_1DRePrimAnd.hxx"
 #include "c2p_2DNoble.hxx"
 
 #include "aster_utils.hxx"
@@ -21,8 +22,8 @@ using namespace Con2PrimFactory;
 using namespace AsterUtils;
 
 enum class eos_3param { IdealGas, Hybrid, Tabulated };
-enum class c2p_first_t { None, Noble, Palenzuela, Entropy };
-enum class c2p_second_t { None, Noble, Palenzuela, Entropy };
+enum class c2p_first_t { None, Noble, RePrimAnd, Palenzuela, Entropy };
+enum class c2p_second_t { None, Noble, RePrimAnd, Palenzuela, Entropy };
 
 enum C2PFlag : CCTK_INT {
   C2P_INIT = 0,    // initial value
@@ -45,6 +46,8 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
 
   if (CCTK_EQUALS(c2p_prime, "Noble")) {
     c2p_fir = c2p_first_t::Noble;
+  } else if (CCTK_EQUALS(c2p_prime, "RePrimAnd")) {
+    c2p_fir = c2p_first_t::RePrimAnd;
   } else if (CCTK_EQUALS(c2p_prime, "Palenzuela")) {
     c2p_fir = c2p_first_t::Palenzuela;
   } else if (CCTK_EQUALS(c2p_prime, "Entropy")) {
@@ -57,6 +60,8 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
 
   if (CCTK_EQUALS(c2p_second, "Noble")) {
     c2p_sec = c2p_second_t::Noble;
+  } else if (CCTK_EQUALS(c2p_second, "RePrimAnd")) {
+    c2p_sec = c2p_second_t::RePrimAnd;
   } else if (CCTK_EQUALS(c2p_second, "Palenzuela")) {
     c2p_sec = c2p_second_t::Palenzuela;
   } else if (CCTK_EQUALS(c2p_second, "Entropy")) {
@@ -141,19 +146,29 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
     c2p_2DNoble c2p_Noble(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, vw_lim,
                           B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
                           inv_beta_max, Ye_lenient, use_z, use_temperature,
-                          use_press_atmo);
+                          use_press_atmo, soft_root_convergence,
+                          soft_root_width_factor);
+
+    // Construct RePrimAnd c2p object:
+    c2p_1DRePrimAnd c2p_RPA(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, vw_lim,
+                            B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
+                            inv_beta_max, Ye_lenient, use_z, use_temperature,
+                            use_press_atmo, soft_root_convergence,
+                            soft_root_width_factor);
 
     // Construct Palenzuela c2p object:
     c2p_1DPalenzuela c2p_Pal(eos_3p, atmo, max_iter, c2p_tol, alp_thresh,
                              vw_lim, B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
                              inv_beta_max, Ye_lenient, use_z, use_temperature,
-                             use_press_atmo);
+                             use_press_atmo, soft_root_convergence,
+                             soft_root_width_factor);
 
     // Construct Entropy c2p object:
     c2p_1DEntropy c2p_Ent(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, vw_lim,
                           B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
                           inv_beta_max, Ye_lenient, use_z, use_temperature,
-                          use_press_atmo);
+                          use_press_atmo, soft_root_convergence,
+                          soft_root_width_factor);
 
     // ----------
 
@@ -270,6 +285,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
                         rep_first);
         break;
       }
+      case c2p_first_t::RePrimAnd: {
+        c2p_RPA.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_first);
+        break;
+      }
       case c2p_first_t::Palenzuela: {
         c2p_Pal.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_first);
         break;
@@ -302,6 +321,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
                           rep_second);
           break;
         }
+        case c2p_second_t::RePrimAnd: {
+          c2p_RPA.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_second);
+          break;
+        }
         case c2p_second_t::Palenzuela: {
           c2p_Pal.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_second);
           break;
@@ -487,6 +510,18 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTS_AsterX_Con2Prim;
   DECLARE_CCTK_PARAMETERS;
 
+  if (CCTK_EQUALS(evolution_eos, "Hybrid") && thermal_eos_atmo) {
+    CCTK_ERROR("Hybrid EOS does not implement *_from_rho_temp_ye; set "
+               "Con2PrimFactory::thermal_eos_atmo = no.");
+  }
+  if (CCTK_EQUALS(evolution_eos, "Tabulated3d") && !use_temperature) {
+    CCTK_ERROR("Tabulated3d requires Con2PrimFactory::use_temperature = yes.");
+  }
+  if (CCTK_EQUALS(evolution_eos, "Tabulated3d") && use_press_atmo) {
+    CCTK_ERROR("Tabulated3d does not support eps_from_rho_press_ye; set "
+               "Con2PrimFactory::use_press_atmo = no.");
+  }
+
   // defining EOS objects
   eos_3param eos_3p_type;
 
@@ -571,21 +606,50 @@ extern "C" void AsterX_Con2Prim_Interpolate_Failed(CCTK_ARGUMENTS) {
           const vec<CCTK_REAL, 6> vely_nbs = get_neighbors(vely, p);
           const vec<CCTK_REAL, 6> velz_nbs = get_neighbors(velz, p);
           const vec<CCTK_REAL, 6> eps_nbs = get_neighbors(eps, p);
+          const vec<CCTK_REAL, 6> press_nbs = get_neighbors(press, p);
           const vec<CCTK_REAL, 6> saved_rho_nbs = get_neighbors(saved_rho, p);
           const vec<CCTK_REAL, 6> saved_velx_nbs = get_neighbors(saved_velx, p);
           const vec<CCTK_REAL, 6> saved_vely_nbs = get_neighbors(saved_vely, p);
           const vec<CCTK_REAL, 6> saved_velz_nbs = get_neighbors(saved_velz, p);
           const vec<CCTK_REAL, 6> saved_eps_nbs = get_neighbors(saved_eps, p);
 
-          CCTK_REAL sum_nbs =
-              sum<6>([&](int i) ARITH_INLINE { return flag_nbs(i); });
-          assert(sum_nbs > 0);
-          rho(p.I) = calc_avg_neighbors(flag_nbs, rho_nbs, saved_rho_nbs);
-          velx(p.I) = calc_avg_neighbors(flag_nbs, velx_nbs, saved_velx_nbs);
-          vely(p.I) = calc_avg_neighbors(flag_nbs, vely_nbs, saved_vely_nbs);
-          velz(p.I) = calc_avg_neighbors(flag_nbs, velz_nbs, saved_velz_nbs);
-          eps(p.I) = calc_avg_neighbors(flag_nbs, eps_nbs, saved_eps_nbs);
-          press(p.I) = (gl_gamma - 1) * eps(p.I) * rho(p.I);
+          const auto good_nb = [&](int i) ARITH_INLINE -> CCTK_REAL {
+            const CCTK_INT flag_i = CCTK_INT(flag_nbs(i));
+            return ((flag_i != C2P_FAIL) && (flag_i != C2P_INIT))
+                       ? CCTK_REAL(1)
+                       : CCTK_REAL(0);
+          };
+
+          const CCTK_REAL sum_nbs =
+              sum<6>([&](int i) ARITH_INLINE { return good_nb(i); });
+          if (sum_nbs <= CCTK_REAL(0)) {
+            return;
+          }
+
+          rho(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                       return good_nb(i) * rho_nbs(i);
+                     }) /
+                     sum_nbs;
+          velx(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                        return good_nb(i) * velx_nbs(i);
+                      }) /
+                      sum_nbs;
+          vely(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                        return good_nb(i) * vely_nbs(i);
+                      }) /
+                      sum_nbs;
+          velz(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                        return good_nb(i) * velz_nbs(i);
+                      }) /
+                      sum_nbs;
+          eps(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                       return good_nb(i) * eps_nbs(i);
+                     }) /
+                     sum_nbs;
+          press(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                         return good_nb(i) * press_nbs(i);
+                       }) /
+                       sum_nbs;
 
           /* reset flag */
           con2prim_flag(p.I) = C2P_AVG;

AsterX/src/estimate_error.cxx

@@ -21,7 +21,7 @@ enum class regrid_t {
 };
 regrid_t regridMethod;
 vector<int> regridVarsI;
-vector<array<int, Loop::dim>> indextypeRegridVars;
+vector<array<int, Loop::dim> > indextypeRegridVars;
 
 extern "C" void AsterX_EstimateError_Setup(CCTK_ARGUMENTS) {
   DECLARE_CCTK_PARAMETERS;

AsterX/src/fluxes.cxx

@@ -113,19 +113,20 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType *eos_3p, const rec_var_t rec_var,
   // Prebind the velocity slice for this direction once
   const auto gf_vel_dir_i = gf_vels(dir_i);
 
-  const auto reconstruct_pt =
-      [=] CCTK_DEVICE(const GF3D2<const CCTK_REAL> &var, const PointDesc &p,
-                      bool gf_is_rho, bool gf_is_press) {
-        return reconstruct<vec<CCTK_REAL, 2>>(var, p, reconstruction, dir_i,
-                                              gf_is_rho, gf_is_press, press,
-                                              gf_vel_dir_i, reconstruct_params);
-      };
+  const auto reconstruct_pt = [=] CCTK_DEVICE(const GF3D2<const CCTK_REAL> &var,
+                                              const PointDesc &p,
+                                              bool gf_is_rho,
+                                              bool gf_is_press) {
+    return reconstruct<vec<CCTK_REAL, 2> >(var, p, reconstruction, dir_i,
+                                           gf_is_rho, gf_is_press, press,
+                                           gf_vel_dir_i, reconstruct_params);
+  };
   const auto reconstruct_loworder =
       [=] CCTK_DEVICE(const GF3D2<const CCTK_REAL> &var, const PointDesc &p,
                       bool gf_is_rho, bool gf_is_press) {
-        return reconstruct<vec<CCTK_REAL, 2>>(var, p, reconstruction_LO, dir_i,
-                                              gf_is_rho, gf_is_press, press,
-                                              gf_vel_dir_i, reconstruct_params);
+        return reconstruct<vec<CCTK_REAL, 2> >(
+            var, p, reconstruction_LO, dir_i, gf_is_rho, gf_is_press, press,
+            gf_vel_dir_i, reconstruct_params);
       };
 
   const auto calcflux =

Con2PrimFactory/interface.ccl

@@ -8,6 +8,7 @@ INCLUDES HEADER: cons.hxx IN cons.hxx
 INCLUDES HEADER: c2p_2DNoble.hxx IN c2p_2DNoble.hxx 
 INCLUDES HEADER: c2p_1DPalenzuela.hxx IN c2p_1DPalenzuela.hxx 
 INCLUDES HEADER: c2p_1DEntropy.hxx IN c2p_1DEntropy.hxx 
+INCLUDES HEADER: c2p_1DRePrimAnd.hxx IN c2p_1DRePrimAnd.hxx
 USES INCLUDE HEADER: setup_eos.hxx 
 USES INCLUDE HEADER: roots.hxx
 USES INCLUDE HEADER: aster_utils.hxx

Con2PrimFactory/param.ccl

@@ -6,6 +6,7 @@ no
 restricted :
 
 KEYWORD c2p_prime "Name of the main con2prim scheme" STEERABLE=ALWAYS {
+  "RePrimAnd" :: "RePrimAnd"
   "Noble" ::"Noble"
   "Palenzuela" ::"Palenzuela"
   "Entropy" ::"Entropy"
@@ -14,6 +15,7 @@ KEYWORD c2p_prime "Name of the main con2prim scheme" STEERABLE=ALWAYS {
 "Noble"
 
 KEYWORD c2p_second "Name of the backup con2prim scheme" STEERABLE=ALWAYS {
+  "RePrimAnd" :: "RePrimAnd"
   "Noble" ::"Noble"
   "Palenzuela" ::"Palenzuela"
   "Entropy" ::"Entropy"
@@ -52,6 +54,7 @@ CCTK_REAL B_lim "Upper limit for the value of magnetization" STEERABLE=ALWAYS
   0:            :: "Must be positive"
 }  100.0
 
+# Safe default for C2P robustness
 CCTK_REAL vw_lim "Upper limit for the value of velocity * lorentz_factor " STEERABLE=ALWAYS
 {
   0:            :: "Must be positive"
@@ -71,6 +74,19 @@ CCTK_REAL c2p_tol "c2p torelance for root finding" STEERABLE=ALWAYS
   0:* :: "Larger than zero"
 } 1e-10
 
+BOOLEAN soft_root_convergence
+"Allow near-converged root solutions to be accepted as success (warning-only)."
+STEERABLE=ALWAYS
+{
+} "no"
+
+CCTK_REAL soft_root_width_factor
+"Relative bracket-width factor (x c2p_tol) used for near-converged root acceptance. Ignored unless soft_root_convergence=yes."
+STEERABLE=ALWAYS
+{
+  1.0:* :: "Must be >= 1"
+} 1.0
+
 CCTK_REAL tauFluid_atmo "Minimum conserved internal energy of fluid" STEERABLE=ALWAYS
 {
   [0.0:* :: "Larger than zero"

Con2PrimFactory/src/c2p.hxx

@@ -51,6 +51,8 @@ protected:
   bool use_zprim;
   bool use_temp;
   bool use_press_atmo;
+  bool soft_root_convergence{false};
+  CCTK_REAL soft_root_width_factor{1.0};
 
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
   get_Ssq_Exact(const vec<CCTK_REAL, 3> &mom,
@@ -447,9 +449,12 @@ c2p::cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv,
   // Compute Bsq
   vec<CCTK_REAL, 3> B_low = calc_contraction(glo, cv.dBvec);
   const CCTK_REAL BsqL = calc_contraction(B_low, cv.dBvec);
+  //const CCTK_REAL tauF_atmo =
+  //    std::max(cv.dens * atmo.eps_atmo, sqrt_detg * tauFluid_atmo);
   const CCTK_REAL tau_lim = 0.5 * BsqL / sqrt_detg;
 
   if (cv.tau <= tau_lim) {
+    //cv.tau = tau_lim + tauF_atmo;
     cv.tau = tau_lim + sqrt_detg * tauFluid_atmo;
   }
 

Con2PrimFactory/src/c2p_1DEntropy.hxx

@@ -17,7 +17,8 @@ public:
       CCTK_REAL tol, CCTK_REAL alp_thresh_in, CCTK_REAL vwlim, CCTK_REAL B_lim,
       CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in,
       CCTK_REAL sigma_max_in, CCTK_REAL inv_beta_max_in, bool ye_len,
-      bool use_z, bool use_temperature, bool use_pressure_atmo);
+      bool use_z, bool use_temperature, bool use_pressure_atmo,
+      bool soft_root_conv, CCTK_REAL soft_root_width_factor_in);
 
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
   get_Ssq_Exact(const vec<CCTK_REAL, 3> &mom,
@@ -66,7 +67,8 @@ CCTK_HOST CCTK_DEVICE
         CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in,
         CCTK_REAL vwlim_BH_in, CCTK_REAL sigma_max_in,
         CCTK_REAL inv_beta_max_in, bool ye_len, bool use_z,
-        bool use_temperature, bool use_pressure_atmo) {
+        bool use_temperature, bool use_pressure_atmo, bool soft_root_conv,
+        CCTK_REAL soft_root_width_factor_in) {
 
   // Base
   atmo = atm;
@@ -86,6 +88,8 @@ CCTK_HOST CCTK_DEVICE
   use_zprim = use_z;
   use_temp = use_temperature;
   use_press_atmo = use_pressure_atmo;
+  soft_root_convergence = soft_root_conv;
+  soft_root_width_factor = fmax(CCTK_REAL(1.0), soft_root_width_factor_in);
 
   // Derived
   GammaIdealFluid = eos_3p->gamma;
@@ -270,10 +274,11 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
                      const CCTK_REAL alp, const vec<CCTK_REAL, 3> &beta,
                      const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const {
 
-  // ROOTSTAT status = ROOTSTAT::SUCCESS;
+  ROOTSTAT status = ROOTSTAT::SUCCESS;
   rep.iters = 0;
   rep.adjust_cons = false;
   rep.set_atmo = false;
+  rep.soft_root_conv = false;
   rep.status = c2p_report::SUCCESS;
 
   /* Check validity of the 3-metric and compute its inverse */
@@ -361,6 +366,12 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
   // const CCTK_INT minbits = std::numeric_limits<CCTK_REAL>::digits - 4;
   const CCTK_INT maxiters = maxIterations;
 
+  const CCTK_REAL f_a0 = fn(a);
+  const CCTK_REAL f_b0 = fn(b);
+  if ((!isfinite(f_a0)) || (!isfinite(f_b0)) || (f_a0 * f_b0 > 0.0)) {
+    status = ROOTSTAT::NOT_BRACKETED;
+  }
+
   auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters);
 
   CCTK_REAL xEntropy_Sol = 0.5 * (result.first + result.second);
@@ -389,14 +400,32 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
     return;
   }
 
-  if (abs(result.first - result.second) >
-      tolerance_0 * min(abs(result.first), abs(result.second))) {
-
-    // set status to root not converged
-    rep.set_root_conv();
-    cv = cv_const;
-    // status = ROOTSTAT::NOT_CONVERGED;
-    return;
+  const CCTK_REAL root_width = abs(result.first - result.second);
+  const CCTK_REAL strict_width_tol =
+      tolerance_0 * min(abs(result.first), abs(result.second));
+  if (root_width > strict_width_tol) {
+    bool accept_soft = false;
+    if (soft_root_convergence) {
+      const CCTK_REAL scale =
+          fmax(CCTK_REAL(0.0), fmax(abs(result.first), abs(result.second)));
+      const CCTK_REAL soft_width_tol =
+          soft_root_width_factor * tolerance * scale;
+      accept_soft = std::isfinite(root_width) && std::isfinite(soft_width_tol) &&
+                    (root_width <= soft_width_tol);
+    }
+    if (!accept_soft) {
+      // set status to root not converged / not bracketed
+      if (status == ROOTSTAT::NOT_BRACKETED) {
+        rep.set_root_bracket();
+      } else {
+        rep.set_root_conv();
+        status = ROOTSTAT::NOT_CONVERGED;
+      }
+      cv = cv_const;
+      (void)status;
+      return;
+    }
+    rep.set_soft_root_conv();
   }
 
   // set to atmo if computed rho is below floor density