From a480c6f6d8758a44a063e7a20367b44e2cc40f82 Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 15 Oct 2025 07:51:29 -0500 Subject: [PATCH 01/10] Con2PrimFactory: add code for RePrimAnd C2P AsterX: add code for RePrimAnd C2P --- AsterX/interface.ccl | 2 +- AsterX/src/con2prim.cxx | 23 +- Con2PrimFactory/interface.ccl | 1 + Con2PrimFactory/param.ccl | 2 + Con2PrimFactory/src/c2p_1DRePrimAnd.hxx | 268 ++++++++++++++++++ .../src/c2p_1DRePrimAnd_internals.hxx | 251 ++++++++++++++++ .../src/c2p_1DRePrimAnd_intervals.hxx | 25 ++ .../src/c2p_1DRePrimAnd_rootfinder.hxx | 67 +++++ 8 files changed, 636 insertions(+), 3 deletions(-) create mode 100644 Con2PrimFactory/src/c2p_1DRePrimAnd.hxx create mode 100644 Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx create mode 100644 Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx create mode 100644 Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx diff --git a/AsterX/interface.ccl b/AsterX/interface.ccl index 3a287bed1..7bfc14043 100644 --- a/AsterX/interface.ccl +++ b/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 diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx index 73ebd3326..c734bd765 100644 --- a/AsterX/src/con2prim.cxx +++ b/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")) { @@ -146,6 +151,12 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p, vwlim_BH, Ye_lenient, use_z, use_temperature, use_press_atmo); + // Construct RePrimAnd c2p object: + c2p_1DRePrimAnd c2p_RPA(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, + cons_error_limit, vw_lim, B_lim, rho_BH, eps_BH, + vwlim_BH, Ye_lenient, use_z, use_temperature, + use_press_atmo); + // Construct Palenzuela c2p object: c2p_1DPalenzuela c2p_Pal(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, cons_error_limit, vw_lim, B_lim, rho_BH, eps_BH, @@ -268,6 +279,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p, c2p_Noble.solve(eos_3p, pv, pv_seeds, cv, glo, rep_first); break; } + case c2p_first_t::RePrimAnd: { + c2p_RPA.solve(eos_3p, pv, cv, glo, rep_first); + break; + } case c2p_first_t::Palenzuela: { c2p_Pal.solve(eos_3p, pv, cv, glo, rep_first); break; @@ -299,6 +314,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p, c2p_Noble.solve(eos_3p, pv, pv_seeds, cv, glo, rep_second); break; } + case c2p_second_t::RePrimAnd: { + c2p_RPA.solve(eos_3p, pv, cv, glo, rep_second); + break; + } case c2p_second_t::Palenzuela: { c2p_Pal.solve(eos_3p, pv, cv, glo, rep_second); break; diff --git a/Con2PrimFactory/interface.ccl b/Con2PrimFactory/interface.ccl index 1ae03a775..35e41db47 100644 --- a/Con2PrimFactory/interface.ccl +++ b/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 diff --git a/Con2PrimFactory/param.ccl b/Con2PrimFactory/param.ccl index 86a936058..bb3ff5f8c 100644 --- a/Con2PrimFactory/param.ccl +++ b/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" diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx new file mode 100644 index 000000000..f837cdef5 --- /dev/null +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx @@ -0,0 +1,268 @@ +#ifndef C2P_1DREPRIMAND_HXX +#define C2P_1DREPRIMAND_HXX + +#include "c2p.hxx" +#include "c2p_report.hxx" +#include "prims.hxx" +#include "cons.hxx" +#include "roots.hxx" + +#include "c2p_1DRePrimAnd_intervals.hxx" +#include "c2p_1DRePrimAnd_rootfinder.hxx" +#include "c2p_1DRePrimAnd_internals.hxx" + +namespace Con2PrimFactory { + +class c2p_1DRePrimAnd : public c2p { +public: + CCTK_REAL GammaIdealFluid; + + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline + c2p_1DRePrimAnd(const EOSType *eos_3p, const atmosphere &atm_in, + CCTK_INT maxIter, CCTK_REAL tol, CCTK_REAL alp_thresh_in, + CCTK_REAL consError, CCTK_REAL vwlim, CCTK_REAL B_lim, + CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, + CCTK_REAL vwlim_BH_in, bool ye_len, bool use_z, + bool use_temperature, bool use_pressure_atmo) { + + atmo = atm_in; + maxIterations = maxIter; + tolerance = tol; + alp_thresh = alp_thresh_in; + cons_error = consError; + vw_lim = vwlim; + w_lim = sqrt(1.0 + vw_lim * vw_lim); + v_lim = vw_lim / w_lim; + Bsq_lim = B_lim * B_lim; + rho_BH = rho_BH_in; + eps_BH = eps_BH_in; + vwlim_BH = vwlim_BH_in; + ye_lenient = ye_len; + use_zprim = use_z; + use_temp = use_temperature; + use_press_atmo = use_pressure_atmo; + + GammaIdealFluid = eos_3p->gamma; + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + set_to_nan(prim_vars &pv, cons_vars &cv) const { + pv.set_to_nan(); + cv.set_to_nan(); + } + + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + solve(const EOSType* eos_3p, + prim_vars& pv, cons_vars& cv, + const smat& glo, + c2p_report& rep) const + { + rep.iters = 0; + rep.adjust_cons = false; + rep.set_atmo = false; + rep.status = c2p_report::SUCCESS; + + const CCTK_REAL spatial_detg = calc_det(glo); + const CCTK_REAL sqrt_detg = sqrt(spatial_detg); + const bool minor1{glo(X,X) > 0.0}; + const bool minor2{glo(X,X)*glo(Y,Y) - glo(X,Y)*glo(X,Y) > 0.0}; + const bool minor3{spatial_detg > 0.0}; + if (!(minor1 && minor2 && minor3)) { + rep.set_invalid_detg(sqrt_detg); + set_to_nan(pv, cv); + return; + } + + const smat gup = calc_inv(glo, spatial_detg); + + const cons_vars cv_const = cv; + + cv.dens /= sqrt_detg; + cv.tau /= sqrt_detg; + cv.mom /= sqrt_detg; + cv.dBvec /= sqrt_detg; + cv.DYe /= sqrt_detg; + cv.DEnt /= sqrt_detg; + + const CCTK_REAL Ssq = calc_contraction(calc_contraction(gup, cv.mom), cv.mom); + const CCTK_REAL Bsq = calc_contraction(calc_contraction(glo, cv.dBvec), cv.dBvec); + const CCTK_REAL BiSi = cv.dBvec(X)*cv.mom(X) + cv.dBvec(Y)*cv.mom(Y) + cv.dBvec(Z)*cv.mom(Z); + + if ((!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) || + (!isfinite(BiSi)) || (!isfinite(cv.DYe)) || (!isfinite(cv.DEnt))) { + rep.set_nans_in_cons(cv.dens, Ssq, Bsq, BiSi, cv.DYe); + set_to_nan(pv, cv); + return; + } + + if (Bsq > Bsq_lim) { + rep.set_B_limit(Bsq); + set_to_nan(pv, cv); + return; + } + + // Ye needs to be a non-const lvalue for EOS calls + CCTK_REAL Ye_raw = cv.DYe / cv.dens; + CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max); + + typename RePrimAnd::froot::cache cache{}; + RePrimAnd::froot f(eos_3p, Ye, + cv.dens, + cv.tau / cv.dens, + Ssq / (cv.dens*cv.dens), + (BiSi*BiSi) / (cv.dens*cv.dens*cv.dens), + Bsq / cv.dens, + cache); + + interval mu_br = f.initial_bracket(); + + interval rgrho_iv{eos_3p->rgrho.min, eos_3p->rgrho.max}; + RePrimAnd::rarecase rc(mu_br, rgrho_iv, f); + mu_br = rc.bracket; + + const CCTK_REAL log2 = std::log(2.0); + const CCTK_INT minbits = int(std::abs(std::log(tolerance)) / log2); + const CCTK_REAL tolerance_0 = std::ldexp(double(1.0), -minbits); + const CCTK_INT maxiters = maxIterations; + + auto fn = [&](CCTK_REAL mu){ return f(mu); }; + if (fn(mu_br.min()) * fn(mu_br.max()) > 0) { + const CCTK_REAL qtot = cv.tau / cv.dens; + const CCTK_REAL sPal = Bsq / cv.dens; + // widen upper bound like Palenzuela’s fallback + CCTK_REAL new_hi = CCTK_REAL(3.0) + CCTK_REAL(3.0) * qtot - CCTK_REAL(1.5) * sPal; + // interval might present mutators; if not, rebuild: + mu_br = interval{mu_br.min(), new_hi}; + } + + auto result = Algo::brent(fn, mu_br.min(), mu_br.max(), minbits, maxiters, rep.iters); + + const CCTK_REAL a_root = result.first; + const CCTK_REAL b_root = result.second; + const CCTK_REAL fa = fn(a_root); + const CCTK_REAL fb = fn(b_root); + const CCTK_REAL mu = + (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) ? b_root : + (std::abs(fa) < std::abs(fb)) ? a_root : + CCTK_REAL(0.5) * (a_root + b_root); + + // RePrimAnd recovery + const CCTK_REAL rhomin = eos_3p->rgrho.min; + CCTK_REAL epsmin = eos_3p->rgeps.min; + CCTK_REAL Ye_tmp = Ye; // EOS may want non-const Ye + const CCTK_REAL pmin = eos_3p->press_from_valid_rho_eps_ye(rhomin, epsmin, Ye_tmp); + const CCTK_REAL h0 = CCTK_REAL(1) + epsmin + pmin / rhomin; + + const CCTK_REAL bsqr = Bsq / cv.dens; + const CCTK_REAL x = CCTK_REAL(1) / (CCTK_REAL(1) + mu * bsqr); + const CCTK_REAL rsqr = Ssq / (cv.dens * cv.dens); + const CCTK_REAL rbsq = (BiSi * BiSi) / (cv.dens * cv.dens * cv.dens); + const CCTK_REAL rfsq = x * (rsqr * x + mu * (x + CCTK_REAL(1)) * rbsq); + const CCTK_REAL W = std::sqrt(h0*h0 + rfsq) / h0; + + pv.rho = cv.dens / (W + 1e-300); + pv.Ye = Ye; + + CCTK_REAL eps_guess = cv.tau / cv.dens; + eps_guess = fmin(fmax(eos_3p->rgeps.min, eps_guess), eos_3p->rgeps.max); + pv.eps = eps_guess; + + CCTK_REAL Ye_tmp2 = pv.Ye; + pv.press = eos_3p->press_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp2); + CCTK_REAL Ye_tmp3 = pv.Ye; + pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp3); + CCTK_REAL Ye_tmp4 = pv.Ye; + pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp4); + + if (use_zprim) { + const vec mom_up = calc_contraction(gup, cv.mom); + const CCTK_REAL mom2 = calc_contraction(mom_up, calc_contraction(glo, mom_up)); + if (mom2 > 0) { + const CCTK_REAL vsq_target = (W*W - 1.0) / (W*W); + const CCTK_REAL vm = std::sqrt(fmax(CCTK_REAL(0), vsq_target)); + pv.vel = mom_up * (vm / (std::sqrt(mom2) + 1e-300)); + } else { + pv.vel = vec{0,0,0}; + } + pv.w_lor = W; + } else { + pv.vel = vec{0,0,0}; + pv.w_lor = W; + } + + vec v_low = calc_contraction(glo, pv.vel); + CCTK_REAL vsq = calc_contraction(v_low, pv.vel); + if (vsq >= v_lim*v_lim) { + const CCTK_REAL scale = v_lim / (std::sqrt(vsq) + 1e-300); + pv.vel *= scale; + vsq = v_lim*v_lim; + } + pv.w_lor = CCTK_REAL(1) / std::sqrt(std::max(CCTK_REAL(1e-16), CCTK_REAL(1) - vsq)); + + pv.Bvec = cv.dBvec; + const vec Elow = calc_cross_product(pv.Bvec, pv.vel); + pv.E = calc_contraction(gup, Elow); + + if (std::abs(result.first - result.second) > + tolerance_0 * std::min(std::abs(result.first), std::abs(result.second))) { + + cons_vars cv_check; + cv_check.from_prim(pv, glo); + + cv_check.dens /= sqrt_detg; + cv_check.tau /= sqrt_detg; + cv_check.mom /= sqrt_detg; + cv_check.dBvec /= sqrt_detg; + cv_check.DYe /= sqrt_detg; + cv_check.DEnt /= sqrt_detg; + + const CCTK_REAL small = 1e-50; + const CCTK_REAL max_error = + sqrt(max({ pow((cv_check.dens - cv.dens )/(cv.dens + small), 2.0), + pow((cv_check.mom(0) - cv.mom(0))/(cv.mom(0) + small), 2.0), + pow((cv_check.mom(1) - cv.mom(1))/(cv.mom(1) + small), 2.0), + pow((cv_check.mom(2) - cv.mom(2))/(cv.mom(2) + small), 2.0), + pow((cv_check.tau - cv.tau )/(cv.tau + small), 2.0) })); + + if (max_error > cons_error) { + rep.set_root_conv(); + cv = cv_const; + return; + } + } + + if (pv.rho < atmo.rho_cut) { + rep.set_atmo_set(); + atmo.set(pv, cv, glo); + return; + } + + c2p::prims_floors_and_ceilings(eos_3p, pv, cv, glo, rep); + + if (rep.adjust_cons) { + cv.from_prim(pv, glo); + cv.dBvec = cv_const.dBvec; + } else { + cv = cv_const; + cv.DEnt = cv.dens * pv.entropy; + } + } + + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + operator()(const EOSType* eos_3p, + prim_vars& pv, cons_vars& cv, + const smat& /*gup_unused*/, + const smat& glo, + c2p_report& rep) const + { + solve(eos_3p, pv, cv, glo, rep); + } +}; + +} // namespace Con2PrimFactory + +#endif + diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx new file mode 100644 index 000000000..0dafae954 --- /dev/null +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx @@ -0,0 +1,251 @@ +#ifndef C2P_1DREPRIMAND_INTERNALS_HXX +#define C2P_1DREPRIMAND_INTERNALS_HXX + +#include +#include +#include "prims.hxx" +#include "cons.hxx" +#include "c2p_1DRePrimAnd_intervals.hxx" +#include "c2p_1DRePrimAnd_rootfinder.hxx" + +namespace Con2PrimFactory { +namespace RePrimAnd { + +// --------------------------- f_upper --------------------------------- + +class f_upper { +public: + using value_t = CCTK_REAL; + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + f_upper(CCTK_REAL h0_, CCTK_REAL rsqr_, CCTK_REAL rbsqr_, CCTK_REAL bsqr_) + : h0(h0_), h0sqr(h0_*h0_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_) {} + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + CCTK_REAL x_from_mu(CCTK_REAL mu) const { return CCTK_REAL(1) / (CCTK_REAL(1) + mu * bsqr); } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + CCTK_REAL rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + return x * (rsqr * x + mu * (x + CCTK_REAL(1)) * rbsqr); + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + CCTK_REAL new_h0w_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + return std::sqrt(h0sqr + rfsqr_from_mu_x(mu, x)); + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + auto operator()(CCTK_REAL mu) const -> std::pair { + const CCTK_REAL x = x_from_mu(mu); + const CCTK_REAL xsq = x * x; + const CCTK_REAL hw = new_h0w_from_mu_x(mu, x); + const CCTK_REAL b = x * (xsq * rsqr + mu * (CCTK_REAL(1) + x + xsq) * rbsqr); + const CCTK_REAL f = mu * hw - CCTK_REAL(1); + const CCTK_REAL df = (h0sqr + b) / hw; + return { f, df }; + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + interval initial_bracket() const { + CCTK_REAL mu_min = CCTK_REAL(1) / std::sqrt(h0sqr + rsqr); + const CCTK_REAL mu0 = CCTK_REAL(1) / h0; + const CCTK_REAL rfs_min = rfsqr_from_mu_x(mu0, x_from_mu(mu0)); + CCTK_REAL mu_max = CCTK_REAL(1) / std::sqrt(h0sqr + rfs_min); + const CCTK_REAL margin = 10 * std::numeric_limits::epsilon(); + mu_max *= CCTK_REAL(1) + margin; + mu_min *= CCTK_REAL(1) - margin; + if (mu_max <= mu_min) { mu_min = CCTK_REAL(0); mu_max = mu0 * (CCTK_REAL(1) + margin); } + return { mu_min, mu_max }; + } + + CCTK_REAL h0, h0sqr, rsqr, rbsqr, bsqr; +}; + +// --------------------------- froot ----------------------------------- + +template +class froot { +public: + using value_t = CCTK_REAL; + + struct cache { + CCTK_REAL ye{0}, lmu{0}, x{0}; + CCTK_REAL rho{0}, rho_raw{0}; + CCTK_REAL eps{0}, eps_raw{0}; + CCTK_REAL press{0}; + CCTK_REAL vsqr{0}; + CCTK_REAL w{1}; + unsigned int calls{0}; + }; + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + froot(const EOSType* eos_, CCTK_REAL valid_ye, + CCTK_REAL d_, CCTK_REAL qtot_, CCTK_REAL rsqr_, CCTK_REAL rbsqr_, CCTK_REAL bsqr_, + cache &last_) + : eos(eos_), d(d_), qtot(qtot_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_), + brosqr(rsqr_ * bsqr_ - rbsqr_), last(last_) + { + last.ye = valid_ye; + last.calls = 0; + + // Build h0 from EOS min bounds (no minimal_h() in EOSX) + const CCTK_REAL rhomin = eos->rgrho.min; + CCTK_REAL epsmin = eos->rgeps.min; + CCTK_REAL ye_tmp = valid_ye; // EOS expects non-const lvalue ref for Ye + const CCTK_REAL pmin = eos->press_from_valid_rho_eps_ye(rhomin, epsmin, ye_tmp); + h0 = CCTK_REAL(1) + epsmin + pmin / rhomin; + + const CCTK_REAL zsqrinf = rsqr / (h0*h0); + const CCTK_REAL wsqrinf = CCTK_REAL(1) + zsqrinf; + winf = std::sqrt(wsqrinf); + vsqrinf = zsqrinf / wsqrinf; + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + CCTK_REAL operator()(CCTK_REAL mu) { + cache &c = last; + c.lmu = mu; + c.x = x_from_mu(mu); + + const CCTK_REAL rfsqr = rfsqr_from_mu_x(mu, c.x); + const CCTK_REAL qf = qf_from_mu_x(mu, c.x); + + c.vsqr = rfsqr * mu * mu; + if (c.vsqr >= vsqrinf) { c.vsqr = vsqrinf; c.w = winf; } + else { c.w = CCTK_REAL(1) / std::sqrt(CCTK_REAL(1) - c.vsqr); } + + c.rho_raw = d / c.w; + c.rho = fmin(fmax(eos->rgrho.min, c.rho_raw), eos->rgrho.max); + + c.eps_raw = get_eps_raw(mu, qf, rfsqr, c.w); + c.eps = fmin(fmax(eos->rgeps.min, c.eps_raw), eos->rgeps.max); + + // EOS again may take Ye by non-const ref; pass a local + CCTK_REAL ye_tmp = c.ye; + c.press = eos->press_from_valid_rho_eps_ye(c.rho, c.eps, ye_tmp); + ++c.calls; + + const CCTK_REAL a = c.press / (c.rho * (CCTK_REAL(1) + c.eps)); + const CCTK_REAL h = (CCTK_REAL(1) + c.eps) * (CCTK_REAL(1) + a); + + const CCTK_REAL hbw_raw = (CCTK_REAL(1) + a) * (CCTK_REAL(1) + qf - mu * rfsqr); + const CCTK_REAL hbw = fmax(hbw_raw, h / c.w); + const CCTK_REAL newmu = CCTK_REAL(1) / (hbw + rfsqr * mu); + + return mu - newmu; + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + interval initial_bracket() const { + CCTK_REAL mu_max = CCTK_REAL(1) / h0; // fallback + { + const int ndigits2 = 36; + RePrimAnd::f_upper g(h0, rsqr, rbsqr, bsqr); + ROOTSTAT status{}; + mu_max = findroot_using_deriv(g, status, ndigits2, ndigits2 + 4); + } + const CCTK_REAL margin = std::ldexp(CCTK_REAL(1), 3 - 36); + mu_max *= (CCTK_REAL(1) + margin); + const CCTK_REAL mu_min = CCTK_REAL(0); + return { mu_min, mu_max }; + } + + CCTK_REAL h0{1}, d{0}, qtot{0}, rsqr{0}, rbsqr{0}, bsqr{0}, brosqr{0}; + CCTK_REAL winf{1}, vsqrinf{0}; + +private: + const EOSType* eos; + cache &last; + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL x_from_mu(CCTK_REAL mu) const { + return CCTK_REAL(1) / (CCTK_REAL(1) + mu * bsqr); + } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + return x * (rsqr * x + mu * (x + CCTK_REAL(1)) * rbsqr); + } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL qf_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + const CCTK_REAL mux = mu * x; + return qtot - (bsqr + mux * mux * brosqr) * CCTK_REAL(0.5); + } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + static CCTK_REAL get_eps_raw(CCTK_REAL mu, CCTK_REAL qf, CCTK_REAL rfsqr, CCTK_REAL w) { + return w * (qf - mu * rfsqr * (CCTK_REAL(1) - mu * w / (CCTK_REAL(1) + w))); + } +}; + +// -------------------- f_rare and rarecase unchanged ------------------- + +template +class f_rare { + const CCTK_REAL v2targ; + const froot &f; +public: + using value_t = CCTK_REAL; + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + f_rare(CCTK_REAL wtarg, const froot &f_) : v2targ(CCTK_REAL(1) - CCTK_REAL(1)/(wtarg*wtarg)), f(f_) {} + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + auto operator()(CCTK_REAL mu) const -> std::pair { + const CCTK_REAL x = CCTK_REAL(1) / (CCTK_REAL(1) + mu * f.bsqr); + const CCTK_REAL xsq = x * x; + const CCTK_REAL rfsq = x * (f.rsqr * x + mu * (x + CCTK_REAL(1)) * f.rbsqr); + const CCTK_REAL vsq = mu * mu * rfsq; + const CCTK_REAL y = vsq - v2targ; + const CCTK_REAL dy = CCTK_REAL(2) * mu * x * (xsq * f.rsqr + mu * (xsq + x + CCTK_REAL(1)) * f.rbsqr); + return { y, dy }; + } +}; + +template +class rarecase { +public: + interval bracket; + bool rho_too_big{false}, rho_big{false}, rho_too_small{false}, rho_small{false}; + + CCTK_HOST CCTK_DEVICE + rarecase(const interval ibracket, + const interval rgrho, + const froot &f) + { + CCTK_REAL muc0 = ibracket.min(); + CCTK_REAL muc1 = ibracket.max(); + const int ndigits = 30; + + if (f.d > rgrho.max()) { + const CCTK_REAL wc = f.d / rgrho.max(); + if (wc > f.winf) { rho_too_big = true; } + else { + f_rare g(wc, f); + if (g(muc1).first <= CCTK_REAL(0)) { rho_too_big = true; } + else if (g(muc0).first < CCTK_REAL(0)) { + ROOTSTAT status{}; + CCTK_REAL mucc = findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); + if (status == ROOTSTAT::SUCCESS) { muc0 = fmax(muc0, mucc); rho_big = true; } else { rho_too_big = true; } + } + } + } + + if (f.d < rgrho.min()) { + const CCTK_REAL wc = f.d / rgrho.min(); + if (wc < CCTK_REAL(1)) { rho_too_small = true; } + else { + f_rare g(wc, f); + if (g(muc0).first >= CCTK_REAL(0)) { rho_too_small = true; } + else if (g(muc1).first > CCTK_REAL(0)) { + ROOTSTAT status{}; + CCTK_REAL mucc = findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); + if (status == ROOTSTAT::SUCCESS) { muc1 = fmin(muc1, mucc); rho_small = true; } else { rho_too_small = true; } + } + } + } + + bracket = interval{ muc0, muc1 }; + } +}; + +} // namespace RePrimAnd +} // namespace Con2PrimFactory + +#endif + diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx new file mode 100644 index 000000000..71f3b0b7c --- /dev/null +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx @@ -0,0 +1,25 @@ +#ifndef C2P_1DREPRIMAND_INTERVALS_HXX +#define C2P_1DREPRIMAND_INTERVALS_HXX + +namespace Con2PrimFactory { + +template +class interval { + T lo_, hi_; +public: + using value_t = T; + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval() : lo_(T(0)), hi_(T(0)) {} + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval(T lo, T hi) : lo_(lo), hi_(hi) {} + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T min() const { return lo_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T max() const { return hi_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T& min() { return lo_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T& max() { return hi_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T width() const { return hi_ - lo_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T mid() const { return (lo_ + hi_) * T(0.5); } +}; + +} // namespace Con2PrimFactory +#endif + diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx new file mode 100644 index 000000000..923d70b47 --- /dev/null +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx @@ -0,0 +1,67 @@ +#ifndef C2P_1DREPRIMAND_ROOTFINDER_HXX +#define C2P_1DREPRIMAND_ROOTFINDER_HXX + +#include +#include "roots.hxx" // internal Brent (Algo::brent) +#include "c2p_1DRePrimAnd_intervals.hxx" +#include "c2p_utils.hxx" + +namespace Con2PrimFactory { + +// Derivative-assisted adapter (we still just use Brent on the residual) +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +T findroot_using_deriv(F& f, interval bracket, ROOTSTAT& errs, + int digits, unsigned int max_calls = 64) +{ + const int minbits = digits; + int iters = 0; + auto fn = [&](T x){ + // Support both f(x)->value and f(x)->pair + auto val = f(x); + if constexpr (std::is_same>::value) return val.first; + else return val; + }; + const T fa = fn(bracket.min()), fb = fn(bracket.max()); + if (fa * fb > T(0)) { errs = ROOTSTAT::NOT_BRACKETED; return std::numeric_limits::quiet_NaN(); } + auto ab = Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); + errs = (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + const T a_root = ab.first, b_root = ab.second; + const T fa2 = fn(a_root), fb2 = fn(b_root); + if (fb2 == T(0) || std::abs(fb2) < std::abs(fa2)) return b_root; + if (std::abs(fa2) < std::abs(fb2)) return a_root; + return (a_root + b_root) * T(0.5); +} + +// Derivative-free Brent adapter +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +interval findroot_no_deriv(F& f, interval bracket, T /*acc*/, + unsigned int max_calls, ROOTSTAT& errs) +{ + // Default fine tolerance like Palenzuela does via minbits + const int minbits = std::numeric_limits::digits - 4; + int iters = 0; + auto fn = [&](T x){ return f(x); }; + const T fa = fn(bracket.min()), fb = fn(bracket.max()); + if (fa * fb > T(0)) { errs = ROOTSTAT::NOT_BRACKETED; return { bracket.min(), bracket.max() }; } + auto ab = Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); + errs = (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + return { ab.first, ab.second }; +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +T findroot_using_deriv(F& f, ROOTSTAT& errs, int digits, unsigned int max_calls = 64) { + return findroot_using_deriv(f, f.initial_bracket(), errs, digits, max_calls); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +interval findroot_no_deriv(F& f, T tol, unsigned int max_calls, ROOTSTAT& errs) { + return findroot_no_deriv(f, f.initial_bracket(), tol, max_calls, errs); +} + +} // namespace Con2PrimFactory +#endif + From 78891dab344a60f0833c2b08acadd93666f98631 Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Thu, 29 Jan 2026 20:01:11 -0600 Subject: [PATCH 02/10] Con2PrimFactory: Bug fixes in RPA, clang-format --- Con2PrimFactory/src/c2p.hxx | 96 +++++++++---------- Con2PrimFactory/src/c2p_1DEntropy.hxx | 5 +- Con2PrimFactory/src/c2p_1DPalenzuela.hxx | 29 +++--- Con2PrimFactory/src/c2p_1DRePrimAnd.hxx | 79 ++++++++------- .../src/c2p_1DRePrimAnd_internals.hxx | 29 ++++-- .../src/c2p_1DRePrimAnd_intervals.hxx | 32 +++++-- .../src/c2p_1DRePrimAnd_rootfinder.hxx | 34 +++++-- Con2PrimFactory/src/c2p_2DNoble.hxx | 6 +- Con2PrimFactory/src/c2p_report.hxx | 6 +- 9 files changed, 181 insertions(+), 135 deletions(-) diff --git a/Con2PrimFactory/src/c2p.hxx b/Con2PrimFactory/src/c2p.hxx index cf3c827d1..7676fa6cc 100644 --- a/Con2PrimFactory/src/c2p.hxx +++ b/Con2PrimFactory/src/c2p.hxx @@ -79,9 +79,9 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv, + cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv, const smat &glo, - const CCTK_REAL &tauFluid_atm) const; + const CCTK_REAL &tauFluid_atm) const; }; template @@ -134,7 +134,8 @@ c2p::prims_floors_and_ceilings(const EOSType *eos_3p, prim_vars &pv, // keeps pressure, changes eps recomp_eps_press_entropy = false; pv.eps = eos_3p->eps_from_valid_rho_press_ye(pv.rho, pv.press, pv.Ye); - pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.temperature = + eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); } } @@ -151,7 +152,8 @@ c2p::prims_floors_and_ceilings(const EOSType *eos_3p, prim_vars &pv, // keeps pressure, changes eps recomp_eps_press_entropy = false; pv.eps = eos_3p->eps_from_valid_rho_press_ye(pv.rho, pv.press, pv.Ye); - pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.temperature = + eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); } } @@ -160,13 +162,12 @@ c2p::prims_floors_and_ceilings(const EOSType *eos_3p, prim_vars &pv, // Ceiling for temperature // Keeps rho the same and changes press // ---------- - + if (pv.temperature > eos_3p->rgtemp.max) { pv.temperature = eos_3p->rgtemp.max; recomp_eps_press_entropy = true; rep.adjust_cons = true; - } // ---------- @@ -184,36 +185,34 @@ c2p::prims_floors_and_ceilings(const EOSType *eos_3p, prim_vars &pv, pv.press = atmo.press_atmo; pv.eps = eos_3p->eps_from_valid_rho_press_ye(pv.rho, pv.press, pv.Ye); - pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.temperature = + eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); recomp_eps_press_entropy = false; rep.adjust_cons = true; - } - + } else { // ---------- // Temperature floor // Keeps rho the same and changes press // ---------- - + if (pv.temperature < atmo.temp_atmo) { - + pv.temperature = atmo.temp_atmo; recomp_eps_press_entropy = true; rep.adjust_cons = true; - } - } if (recomp_eps_press_entropy) { pv.eps = eos_3p->eps_from_valid_rho_temp_ye(pv.rho, pv.temperature, pv.Ye); - pv.press = eos_3p->press_from_valid_rho_temp_ye(pv.rho, pv.temperature, pv.Ye); + pv.press = + eos_3p->press_from_valid_rho_temp_ye(pv.rho, pv.temperature, pv.Ye); pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); } - } template @@ -253,26 +252,27 @@ c2p::bh_interior(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, }; if (recomp_flag) { - - pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - pv.press = eos_3p->press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + + pv.temperature = + eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.press = eos_3p->press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - + cv.from_prim(pv, glo); }; - + } else { - + pv.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of initial // NS or disk pv.eps = eps_BH; pv.Ye = atmo.ye_atmo; - + pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - pv.press = eos_3p->press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.press = eos_3p->press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - // Set velocity such that new conserved momentum has same + // Set velocity such that new conserved momentum has same // direction as before // Inverse metric @@ -280,7 +280,8 @@ c2p::bh_interior(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, const smat gup = calc_inv(glo, spatial_detg); // Compute Z = rho * h * W * W - const CCTK_REAL Z_loc = ( pv.rho * ( 1.0 + pv.eps ) + pv.press ) * wlim_BH * wlim_BH; + const CCTK_REAL Z_loc = + (pv.rho * (1.0 + pv.eps) + pv.press) * wlim_BH * wlim_BH; // Get Bsq const vec B_low = calc_contraction(glo, pv.Bvec); @@ -288,7 +289,7 @@ c2p::bh_interior(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, // Norm of conserved momentum, undensitize here vec mom_low = cv.mom / sqrt(spatial_detg); - vec mom_up = calc_contraction(gup, mom_low); + vec mom_up = calc_contraction(gup, mom_low); const CCTK_REAL Ssq_old = calc_contraction(mom_low, mom_up); const CCTK_REAL S_old = sqrt(Ssq_old) + 1e-50; @@ -299,40 +300,37 @@ c2p::bh_interior(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, const CCTK_REAL BiEsi = BiSi_old / S_old; // Compute magnitude of new conserved momentum - const CCTK_REAL Ssq_new = ( (Z_loc + Bsq)*(Z_loc + Bsq)*vlim_BH*vlim_BH ) / - ( 1.0 + BiEsi * BiEsi * ( 2.0 * Z_loc + Bsq ) / ( Z_loc * Z_loc ) ); + const CCTK_REAL Ssq_new = + ((Z_loc + Bsq) * (Z_loc + Bsq) * vlim_BH * vlim_BH) / + (1.0 + BiEsi * BiEsi * (2.0 * Z_loc + Bsq) / (Z_loc * Z_loc)); const CCTK_REAL S_new = sqrt(Ssq_new); // Rescale momenta mom_low *= S_new / S_old; - mom_up *= S_new / S_old; + mom_up *= S_new / S_old; - // Finally, compute velocity + // Finally, compute velocity // This is (24) from https://arxiv.org/pdf/1712.07538 - pv.vel(X) = mom_up(X) / - (Z_loc + Bsq); + pv.vel(X) = mom_up(X) / (Z_loc + Bsq); pv.vel(X) += BiEsi * S_new * pv.Bvec(X) / (Z_loc * (Z_loc + Bsq)); - pv.vel(Y) = mom_up(Y) / - (Z_loc + Bsq); + pv.vel(Y) = mom_up(Y) / (Z_loc + Bsq); pv.vel(Y) += BiEsi * S_new * pv.Bvec(Y) / (Z_loc * (Z_loc + Bsq)); - pv.vel(Z) = mom_up(Z) / - (Z_loc + Bsq); + pv.vel(Z) = mom_up(Z) / (Z_loc + Bsq); pv.vel(Z) += BiEsi * S_new * pv.Bvec(Z) / (Z_loc * (Z_loc + Bsq)); pv.w_lor = wlim_BH; cv.from_prim(pv, glo); - }; }; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void -c2p::cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv, +c2p::cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv, const smat &glo, - const CCTK_REAL &tauFluid_atmo) const { + const CCTK_REAL &tauFluid_atmo) const { // Limit conservative variables // Note that conservatives are densitized @@ -346,29 +344,29 @@ c2p::cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv, // Based on Appendix A of https://arxiv.org/pdf/1112.0568 // Compute Bsq - vec B_low = calc_contraction(glo, cv.dBvec); + vec 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; + 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; } - // Dominant energy condition + // Dominant energy condition // (A5) from https://arxiv.org/pdf/1505.01607 - vec mom_up = calc_contraction(gup, cv.mom); + vec mom_up = calc_contraction(gup, cv.mom); const CCTK_REAL mom2L = calc_contraction(cv.mom, mom_up); - const CCTK_REAL slim = cv.dens + cv.tau; - const CCTK_REAL slim2 = slim*slim; + const CCTK_REAL slim = cv.dens + cv.tau; + const CCTK_REAL slim2 = slim * slim; if (mom2L > slim2) { - // (A51) from https://arxiv.org/pdf/1112.0568 - cv.mom = cv.mom * sqrt(slim2/mom2L); + // (A51) from https://arxiv.org/pdf/1112.0568 + cv.mom = cv.mom * sqrt(slim2 / mom2L); } - }; } // namespace Con2PrimFactory diff --git a/Con2PrimFactory/src/c2p_1DEntropy.hxx b/Con2PrimFactory/src/c2p_1DEntropy.hxx index ce2e7cdc0..2f034e54d 100644 --- a/Con2PrimFactory/src/c2p_1DEntropy.hxx +++ b/Con2PrimFactory/src/c2p_1DEntropy.hxx @@ -446,8 +446,9 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, // Computing b^mu b_mu const CCTK_REAL bs2 = (Bsq + bst * bst) / (pv.w_lor * pv.w_lor); // Recompute tau - cv.tau = sqrt_detg * (pv.w_lor * pv.w_lor * (pv.rho * (1.0 + pv.eps) + pv.press + bs2) - - (pv.press + 0.5 * bs2) - bst * bst) - + cv.tau = sqrt_detg * (pv.w_lor * pv.w_lor * + (pv.rho * (1.0 + pv.eps) + pv.press + bs2) - + (pv.press + 0.5 * bs2) - bst * bst) - cv.dens; } } diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx index 800f5c60e..f7f098aa2 100644 --- a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx +++ b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx @@ -173,11 +173,11 @@ c2p_1DPalenzuela::xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, sPalenzuela / (2.0 * W_sol * W_sol)); // TODO: Using this check here can lead to corrections of negative eps - // which could be accepted in certain cases. Thus, these cases will - // not be marked as failure after the solving for the root. Move - // the check to the tabulated EOS framework later. + // which could be accepted in certain cases. Thus, these cases will + // not be marked as failure after the solving for the root. Move + // the check to the tabulated EOS framework later. if (use_temp) { - pv.eps = std::max(pv.eps, atmo.eps_atmo); // check on lower bound + pv.eps = std::max(pv.eps, atmo.eps_atmo); // check on lower bound pv.eps = std::min(pv.eps, eos_3p->rgeps.max); // check on upper bound } @@ -288,11 +288,11 @@ c2p_1DPalenzuela::funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, sPalenzuela / (2 * W_loc * W_loc)); // TODO: Using this check here can lead to corrections of negative eps - // which could be accepted in certain cases. Thus, these cases will - // not be marked as failure after the solving for the root. Move - // the check to the tabulated EOS framework later. + // which could be accepted in certain cases. Thus, these cases will + // not be marked as failure after the solving for the root. Move + // the check to the tabulated EOS framework later. if (use_temp) { - eps_loc = std::max(eps_loc, atmo.eps_atmo); // check on lower bound + eps_loc = std::max(eps_loc, atmo.eps_atmo); // check on lower bound eps_loc = std::min(eps_loc, eos_3p->rgeps.max); // check on upper bound } @@ -408,13 +408,14 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, // Dominant energy check if (fn(a) * fn(b) > 0) { - //printf("for fn(a)*fn(b)>0, fa, fb: %26.16e, %26.16e \n\n", fn(a), fn(b)); + // printf("for fn(a)*fn(b)>0, fa, fb: %26.16e, %26.16e \n\n", fn(a), fn(b)); b = 3.0 + 3.0 * qPalenzuela - 1.5 * sPalenzuela; - //if (fn(a) * fn(b) < 0) { - // printf("for fn(a)*fn(b)<0, fa, fb: %26.16e, %26.16e \n\n", fn(a), fn(b)); - // printf( - // "Palenzuela C2P: dominant energy condition has been violated!\n\n"); - //} + // if (fn(a) * fn(b) < 0) { + // printf("for fn(a)*fn(b)<0, fa, fb: %26.16e, %26.16e \n\n", fn(a), + // fn(b)); printf( + // "Palenzuela C2P: dominant energy condition has been + // violated!\n\n"); + // } } auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters); diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx index f837cdef5..2de45b2d5 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx @@ -1,6 +1,9 @@ #ifndef C2P_1DREPRIMAND_HXX #define C2P_1DREPRIMAND_HXX +#include +#include + #include "c2p.hxx" #include "c2p_report.hxx" #include "prims.hxx" @@ -131,9 +134,7 @@ public: if (fn(mu_br.min()) * fn(mu_br.max()) > 0) { const CCTK_REAL qtot = cv.tau / cv.dens; const CCTK_REAL sPal = Bsq / cv.dens; - // widen upper bound like Palenzuela’s fallback CCTK_REAL new_hi = CCTK_REAL(3.0) + CCTK_REAL(3.0) * qtot - CCTK_REAL(1.5) * sPal; - // interval might present mutators; if not, rebuild: mu_br = interval{mu_br.min(), new_hi}; } @@ -148,48 +149,54 @@ public: (std::abs(fa) < std::abs(fb)) ? a_root : CCTK_REAL(0.5) * (a_root + b_root); - // RePrimAnd recovery - const CCTK_REAL rhomin = eos_3p->rgrho.min; - CCTK_REAL epsmin = eos_3p->rgeps.min; - CCTK_REAL Ye_tmp = Ye; // EOS may want non-const Ye - const CCTK_REAL pmin = eos_3p->press_from_valid_rho_eps_ye(rhomin, epsmin, Ye_tmp); - const CCTK_REAL h0 = CCTK_REAL(1) + epsmin + pmin / rhomin; - - const CCTK_REAL bsqr = Bsq / cv.dens; - const CCTK_REAL x = CCTK_REAL(1) / (CCTK_REAL(1) + mu * bsqr); - const CCTK_REAL rsqr = Ssq / (cv.dens * cv.dens); - const CCTK_REAL rbsq = (BiSi * BiSi) / (cv.dens * cv.dens * cv.dens); - const CCTK_REAL rfsq = x * (rsqr * x + mu * (x + CCTK_REAL(1)) * rbsq); - const CCTK_REAL W = std::sqrt(h0*h0 + rfsq) / h0; - - pv.rho = cv.dens / (W + 1e-300); - pv.Ye = Ye; - - CCTK_REAL eps_guess = cv.tau / cv.dens; - eps_guess = fmin(fmax(eos_3p->rgeps.min, eps_guess), eos_3p->rgeps.max); - pv.eps = eps_guess; - - CCTK_REAL Ye_tmp2 = pv.Ye; - pv.press = eos_3p->press_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp2); + // ------------------------------------------------------------------ + // IMPORTANT: ensure cache corresponds to the final chosen mu + // ------------------------------------------------------------------ + (void)f(mu); + + // ------------------------------------------------------------------ + // Use the EOS-consistent RePrimAnd cached primitives + // ------------------------------------------------------------------ + pv.rho = cache.rho; + pv.Ye = Ye; + pv.eps = cache.eps; + pv.press = cache.press; + pv.w_lor = cache.w; + CCTK_REAL Ye_tmp3 = pv.Ye; pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp3); CCTK_REAL Ye_tmp4 = pv.Ye; pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp4); + // ------------------------------------------------------------------ + // Velocity: RePrimAnd magnetic-aware reconstruction + // v^i = mu * x * ( r^i + (r.b) * mu * b^i ) + // where r^i = S^i / D, b^i = B^i / sqrt(D), (r.b) = (B^i S_i)/(D^(3/2)). + // ------------------------------------------------------------------ if (use_zprim) { const vec mom_up = calc_contraction(gup, cv.mom); - const CCTK_REAL mom2 = calc_contraction(mom_up, calc_contraction(glo, mom_up)); - if (mom2 > 0) { - const CCTK_REAL vsq_target = (W*W - 1.0) / (W*W); - const CCTK_REAL vm = std::sqrt(fmax(CCTK_REAL(0), vsq_target)); - pv.vel = mom_up * (vm / (std::sqrt(mom2) + 1e-300)); - } else { - pv.vel = vec{0,0,0}; - } - pv.w_lor = W; + + const CCTK_REAL D = cv.dens; + const CCTK_REAL sqD = std::sqrt(D + CCTK_REAL(0)); + + const vec r_u = mom_up * (CCTK_REAL(1) / (D + CCTK_REAL(1e-300))); + const vec b_u = cv.dBvec * (CCTK_REAL(1) / (sqD + CCTK_REAL(1e-300))); + const CCTK_REAL rb = BiSi / ((D + CCTK_REAL(1e-300)) * (sqD + CCTK_REAL(1e-300))); + + pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u); + pv.w_lor = cache.w; } else { - pv.vel = vec{0,0,0}; - pv.w_lor = W; + const vec mom_up = calc_contraction(gup, cv.mom); + + const CCTK_REAL D = cv.dens; + const CCTK_REAL sqD = std::sqrt(D + CCTK_REAL(0)); + + const vec r_u = mom_up * (CCTK_REAL(1) / (D + CCTK_REAL(1e-300))); + const vec b_u = cv.dBvec * (CCTK_REAL(1) / (sqD + CCTK_REAL(1e-300))); + const CCTK_REAL rb = BiSi / ((D + CCTK_REAL(1e-300)) * (sqD + CCTK_REAL(1e-300))); + + pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u); + pv.w_lor = cache.w; } vec v_low = calc_contraction(glo, pv.vel); diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx index 0dafae954..61fc35625 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx @@ -3,6 +3,7 @@ #include #include + #include "prims.hxx" #include "cons.hxx" #include "c2p_1DRePrimAnd_intervals.hxx" @@ -88,7 +89,7 @@ public: last.ye = valid_ye; last.calls = 0; - // Build h0 from EOS min bounds (no minimal_h() in EOSX) + // Build h0 from EOS min bounds const CCTK_REAL rhomin = eos->rgrho.min; CCTK_REAL epsmin = eos->rgeps.min; CCTK_REAL ye_tmp = valid_ye; // EOS expects non-const lvalue ref for Ye @@ -120,7 +121,6 @@ public: c.eps_raw = get_eps_raw(mu, qf, rfsqr, c.w); c.eps = fmin(fmax(eos->rgeps.min, c.eps_raw), eos->rgeps.max); - // EOS again may take Ye by non-const ref; pass a local CCTK_REAL ye_tmp = c.ye; c.press = eos->press_from_valid_rho_eps_ye(c.rho, c.eps, ye_tmp); ++c.calls; @@ -143,6 +143,7 @@ public: RePrimAnd::f_upper g(h0, rsqr, rbsqr, bsqr); ROOTSTAT status{}; mu_max = findroot_using_deriv(g, status, ndigits2, ndigits2 + 4); + if (status != ROOTSTAT::SUCCESS || !std::isfinite(mu_max)) mu_max = CCTK_REAL(1) / h0; } const CCTK_REAL margin = std::ldexp(CCTK_REAL(1), 3 - 36); mu_max *= (CCTK_REAL(1) + margin); @@ -150,6 +151,7 @@ public: return { mu_min, mu_max }; } + // Public parameters accessed by rarecase/f_rare CCTK_REAL h0{1}, d{0}, qtot{0}, rsqr{0}, rbsqr{0}, bsqr{0}, brosqr{0}; CCTK_REAL winf{1}, vsqrinf{0}; @@ -157,23 +159,29 @@ private: const EOSType* eos; cache &last; - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL x_from_mu(CCTK_REAL mu) const { + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + CCTK_REAL x_from_mu(CCTK_REAL mu) const { return CCTK_REAL(1) / (CCTK_REAL(1) + mu * bsqr); } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + CCTK_REAL rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { return x * (rsqr * x + mu * (x + CCTK_REAL(1)) * rbsqr); } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL qf_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + CCTK_REAL qf_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { const CCTK_REAL mux = mu * x; return qtot - (bsqr + mux * mux * brosqr) * CCTK_REAL(0.5); } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE static CCTK_REAL get_eps_raw(CCTK_REAL mu, CCTK_REAL qf, CCTK_REAL rfsqr, CCTK_REAL w) { return w * (qf - mu * rfsqr * (CCTK_REAL(1) - mu * w / (CCTK_REAL(1) + w))); } }; -// -------------------- f_rare and rarecase unchanged ------------------- +// -------------------- f_rare and rarecase ------------------- template class f_rare { @@ -183,7 +191,8 @@ public: using value_t = CCTK_REAL; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - f_rare(CCTK_REAL wtarg, const froot &f_) : v2targ(CCTK_REAL(1) - CCTK_REAL(1)/(wtarg*wtarg)), f(f_) {} + f_rare(CCTK_REAL wtarg, const froot &f_) + : v2targ(CCTK_REAL(1) - CCTK_REAL(1)/(wtarg*wtarg)), f(f_) {} CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE auto operator()(CCTK_REAL mu) const -> std::pair { @@ -221,7 +230,8 @@ public: else if (g(muc0).first < CCTK_REAL(0)) { ROOTSTAT status{}; CCTK_REAL mucc = findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); - if (status == ROOTSTAT::SUCCESS) { muc0 = fmax(muc0, mucc); rho_big = true; } else { rho_too_big = true; } + if (status == ROOTSTAT::SUCCESS) { muc0 = fmax(muc0, mucc); rho_big = true; } + else { rho_too_big = true; } } } } @@ -235,7 +245,8 @@ public: else if (g(muc1).first > CCTK_REAL(0)) { ROOTSTAT status{}; CCTK_REAL mucc = findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); - if (status == ROOTSTAT::SUCCESS) { muc1 = fmin(muc1, mucc); rho_small = true; } else { rho_too_small = true; } + if (status == ROOTSTAT::SUCCESS) { muc1 = fmin(muc1, mucc); rho_small = true; } + else { rho_too_small = true; } } } } diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx index 71f3b0b7c..7706bb58f 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx @@ -9,15 +9,29 @@ class interval { public: using value_t = T; - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval() : lo_(T(0)), hi_(T(0)) {} - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval(T lo, T hi) : lo_(lo), hi_(hi) {} - - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T min() const { return lo_; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T max() const { return hi_; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T& min() { return lo_; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T& max() { return hi_; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T width() const { return hi_ - lo_; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T mid() const { return (lo_ + hi_) * T(0.5); } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + interval() : lo_(T(0)), hi_(T(0)) {} + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + interval(T lo, T hi) : lo_(lo), hi_(hi) {} + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + T min() const { return lo_; } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + T max() const { return hi_; } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + T& min() { return lo_; } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + T& max() { return hi_; } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + T width() const { return hi_ - lo_; } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE + T mid() const { return (lo_ + hi_) * T(0.5); } }; } // namespace Con2PrimFactory diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx index 923d70b47..8ddc2374e 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx @@ -2,7 +2,11 @@ #define C2P_1DREPRIMAND_ROOTFINDER_HXX #include -#include "roots.hxx" // internal Brent (Algo::brent) +#include +#include +#include + +#include "roots.hxx" // Algo::brent #include "c2p_1DRePrimAnd_intervals.hxx" #include "c2p_utils.hxx" @@ -10,24 +14,31 @@ namespace Con2PrimFactory { // Derivative-assisted adapter (we still just use Brent on the residual) template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T findroot_using_deriv(F& f, interval bracket, ROOTSTAT& errs, int digits, unsigned int max_calls = 64) { const int minbits = digits; int iters = 0; - auto fn = [&](T x){ - // Support both f(x)->value and f(x)->pair + + auto fn = [&](T x) -> T { auto val = f(x); if constexpr (std::is_same>::value) return val.first; else return val; }; + const T fa = fn(bracket.min()), fb = fn(bracket.max()); - if (fa * fb > T(0)) { errs = ROOTSTAT::NOT_BRACKETED; return std::numeric_limits::quiet_NaN(); } + if (fa * fb > T(0)) { + errs = ROOTSTAT::NOT_BRACKETED; + return std::numeric_limits::quiet_NaN(); + } + auto ab = Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); errs = (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + const T a_root = ab.first, b_root = ab.second; const T fa2 = fn(a_root), fb2 = fn(b_root); + if (fb2 == T(0) || std::abs(fb2) < std::abs(fa2)) return b_root; if (std::abs(fa2) < std::abs(fb2)) return a_root; return (a_root + b_root) * T(0.5); @@ -35,29 +46,32 @@ T findroot_using_deriv(F& f, interval bracket, ROOTSTAT& errs, // Derivative-free Brent adapter template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval findroot_no_deriv(F& f, interval bracket, T /*acc*/, unsigned int max_calls, ROOTSTAT& errs) { - // Default fine tolerance like Palenzuela does via minbits const int minbits = std::numeric_limits::digits - 4; int iters = 0; - auto fn = [&](T x){ return f(x); }; + + auto fn = [&](T x) -> T { return f(x); }; + const T fa = fn(bracket.min()), fb = fn(bracket.max()); if (fa * fb > T(0)) { errs = ROOTSTAT::NOT_BRACKETED; return { bracket.min(), bracket.max() }; } + auto ab = Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); errs = (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + return { ab.first, ab.second }; } template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T findroot_using_deriv(F& f, ROOTSTAT& errs, int digits, unsigned int max_calls = 64) { return findroot_using_deriv(f, f.initial_bracket(), errs, digits, max_calls); } template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval findroot_no_deriv(F& f, T tol, unsigned int max_calls, ROOTSTAT& errs) { return findroot_no_deriv(f, f.initial_bracket(), tol, max_calls, errs); } diff --git a/Con2PrimFactory/src/c2p_2DNoble.hxx b/Con2PrimFactory/src/c2p_2DNoble.hxx index 28e5b4047..944802983 100644 --- a/Con2PrimFactory/src/c2p_2DNoble.hxx +++ b/Con2PrimFactory/src/c2p_2DNoble.hxx @@ -402,7 +402,7 @@ c2p_2DNoble::solve(const EOSType *eos_3p, prim_vars &pv, prim_vars &pv_seeds, /* Start Recovery with 2D NR Solver */ constexpr CCTK_INT n = 2; constexpr CCTK_REAL dv = (1. - 1.e-10); - //constexpr CCTK_REAL dw = 1. / (1. - dv); + // constexpr CCTK_REAL dw = 1. / (1. - dv); CCTK_REAL dx[n]; CCTK_REAL fjac[n][n]; @@ -424,7 +424,7 @@ c2p_2DNoble::solve(const EOSType *eos_3p, prim_vars &pv, prim_vars &pv_seeds, } if (x[0] <= 0.0) { - x[0] = fabs(x[0])+Zmin; + x[0] = fabs(x[0]) + Zmin; } else { if (x[0] > 1e20) { x[0] = x_old[0]; @@ -490,7 +490,7 @@ c2p_2DNoble::solve(const EOSType *eos_3p, prim_vars &pv, prim_vars &pv_seeds, } if (x[0] <= 0.0) { - x[0] = fabs(x[0])+Zmin; + x[0] = fabs(x[0]) + Zmin; } else { if (x[0] > 1e20) { x[0] = x_old[0]; diff --git a/Con2PrimFactory/src/c2p_report.hxx b/Con2PrimFactory/src/c2p_report.hxx index f612015b5..dfc2905ed 100644 --- a/Con2PrimFactory/src/c2p_report.hxx +++ b/Con2PrimFactory/src/c2p_report.hxx @@ -233,8 +233,7 @@ public: printf("Density out of range, dens = %16.8e, rho = %16.8e \n", dens, rho); break; case RANGE_EPS: - printf("Specific energy was out of range! eps = %16.8e \n", - eps); + printf("Specific energy was out of range! eps = %16.8e \n", eps); break; case SPEED_LIMIT: printf("Speed limit exceeded, vx, vy, vz = %16.8e, %16.8e, %16.8e \n", @@ -259,7 +258,8 @@ public: printf("Preparatory root finding failed (faulty bracketing) \n"); break; case ERR_CODE_NOT_SET: - printf("Error code not set! Report was initialized but never changed. \n"); + printf( + "Error code not set! Report was initialized but never changed. \n"); break; default: assert(false); From 4ef6a959eebf54bbf2022d516644551b124ff9ac Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 4 Feb 2026 09:58:59 -0600 Subject: [PATCH 03/10] Con2PrimFactory: clean up --- Con2PrimFactory/src/c2p_1DRePrimAnd.hxx | 159 +++++------- .../src/c2p_1DRePrimAnd_internals.hxx | 242 ++++++++++-------- .../src/c2p_1DRePrimAnd_intervals.hxx | 39 +-- .../src/c2p_1DRePrimAnd_rootfinder.hxx | 65 +++-- 4 files changed, 255 insertions(+), 250 deletions(-) diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx index 2de45b2d5..e9b656013 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx @@ -21,13 +21,12 @@ public: CCTK_REAL GammaIdealFluid; template - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline - c2p_1DRePrimAnd(const EOSType *eos_3p, const atmosphere &atm_in, - CCTK_INT maxIter, CCTK_REAL tol, CCTK_REAL alp_thresh_in, - CCTK_REAL consError, CCTK_REAL vwlim, CCTK_REAL B_lim, - CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, - CCTK_REAL vwlim_BH_in, bool ye_len, bool use_z, - bool use_temperature, bool use_pressure_atmo) { + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DRePrimAnd( + const EOSType *eos_3p, const atmosphere &atm_in, CCTK_INT maxIter, + CCTK_REAL tol, CCTK_REAL alp_thresh_in, CCTK_REAL consError, + CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, + CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, bool ye_len, bool use_z, + bool use_temperature, bool use_pressure_atmo) { atmo = atm_in; maxIterations = maxIter; @@ -57,11 +56,8 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - solve(const EOSType* eos_3p, - prim_vars& pv, cons_vars& cv, - const smat& glo, - c2p_report& rep) const - { + solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, + const smat &glo, c2p_report &rep) const { rep.iters = 0; rep.adjust_cons = false; rep.set_atmo = false; @@ -69,8 +65,8 @@ public: const CCTK_REAL spatial_detg = calc_det(glo); const CCTK_REAL sqrt_detg = sqrt(spatial_detg); - const bool minor1{glo(X,X) > 0.0}; - const bool minor2{glo(X,X)*glo(Y,Y) - glo(X,Y)*glo(X,Y) > 0.0}; + const bool minor1{glo(X, X) > 0.0}; + const bool minor2{glo(X, X) * glo(Y, Y) - glo(X, Y) * glo(X, Y) > 0.0}; const bool minor3{spatial_detg > 0.0}; if (!(minor1 && minor2 && minor3)) { rep.set_invalid_detg(sqrt_detg); @@ -78,20 +74,23 @@ public: return; } - const smat gup = calc_inv(glo, spatial_detg); + const smat gup = calc_inv(glo, spatial_detg); const cons_vars cv_const = cv; - cv.dens /= sqrt_detg; - cv.tau /= sqrt_detg; - cv.mom /= sqrt_detg; + cv.dens /= sqrt_detg; + cv.tau /= sqrt_detg; + cv.mom /= sqrt_detg; cv.dBvec /= sqrt_detg; - cv.DYe /= sqrt_detg; - cv.DEnt /= sqrt_detg; + cv.DYe /= sqrt_detg; + cv.DEnt /= sqrt_detg; - const CCTK_REAL Ssq = calc_contraction(calc_contraction(gup, cv.mom), cv.mom); - const CCTK_REAL Bsq = calc_contraction(calc_contraction(glo, cv.dBvec), cv.dBvec); - const CCTK_REAL BiSi = cv.dBvec(X)*cv.mom(X) + cv.dBvec(Y)*cv.mom(Y) + cv.dBvec(Z)*cv.mom(Z); + const CCTK_REAL Ssq = + calc_contraction(calc_contraction(gup, cv.mom), cv.mom); + const CCTK_REAL Bsq = + calc_contraction(calc_contraction(glo, cv.dBvec), cv.dBvec); + const CCTK_REAL BiSi = cv.dBvec(X) * cv.mom(X) + cv.dBvec(Y) * cv.mom(Y) + + cv.dBvec(Z) * cv.mom(Z); if ((!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) || (!isfinite(BiSi)) || (!isfinite(cv.DYe)) || (!isfinite(cv.DEnt))) { @@ -108,16 +107,12 @@ public: // Ye needs to be a non-const lvalue for EOS calls CCTK_REAL Ye_raw = cv.DYe / cv.dens; - CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max); + CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max); typename RePrimAnd::froot::cache cache{}; - RePrimAnd::froot f(eos_3p, Ye, - cv.dens, - cv.tau / cv.dens, - Ssq / (cv.dens*cv.dens), - (BiSi*BiSi) / (cv.dens*cv.dens*cv.dens), - Bsq / cv.dens, - cache); + RePrimAnd::froot f( + eos_3p, Ye, cv.dens, cv.tau / cv.dens, Ssq / (cv.dens * cv.dens), + (BiSi * BiSi) / (cv.dens * cv.dens * cv.dens), Bsq / cv.dens, cache); interval mu_br = f.initial_bracket(); @@ -126,28 +121,29 @@ public: mu_br = rc.bracket; const CCTK_REAL log2 = std::log(2.0); - const CCTK_INT minbits = int(std::abs(std::log(tolerance)) / log2); + const CCTK_INT minbits = int(std::abs(std::log(tolerance)) / log2); const CCTK_REAL tolerance_0 = std::ldexp(double(1.0), -minbits); - const CCTK_INT maxiters = maxIterations; + const CCTK_INT maxiters = maxIterations; - auto fn = [&](CCTK_REAL mu){ return f(mu); }; + auto fn = [&](CCTK_REAL mu) { return f(mu); }; if (fn(mu_br.min()) * fn(mu_br.max()) > 0) { const CCTK_REAL qtot = cv.tau / cv.dens; const CCTK_REAL sPal = Bsq / cv.dens; - CCTK_REAL new_hi = CCTK_REAL(3.0) + CCTK_REAL(3.0) * qtot - CCTK_REAL(1.5) * sPal; + CCTK_REAL new_hi = 3.0 + 3.0 * qtot - 1.5 * sPal; mu_br = interval{mu_br.min(), new_hi}; } - auto result = Algo::brent(fn, mu_br.min(), mu_br.max(), minbits, maxiters, rep.iters); + auto result = + Algo::brent(fn, mu_br.min(), mu_br.max(), minbits, maxiters, rep.iters); const CCTK_REAL a_root = result.first; const CCTK_REAL b_root = result.second; const CCTK_REAL fa = fn(a_root); const CCTK_REAL fb = fn(b_root); const CCTK_REAL mu = - (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) ? b_root : - (std::abs(fa) < std::abs(fb)) ? a_root : - CCTK_REAL(0.5) * (a_root + b_root); + (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) ? b_root + : (std::abs(fa) < std::abs(fb)) ? a_root + : CCTK_REAL(0.5) * (a_root + b_root); // ------------------------------------------------------------------ // IMPORTANT: ensure cache corresponds to the final chosen mu @@ -157,81 +153,66 @@ public: // ------------------------------------------------------------------ // Use the EOS-consistent RePrimAnd cached primitives // ------------------------------------------------------------------ - pv.rho = cache.rho; - pv.Ye = Ye; - pv.eps = cache.eps; + pv.rho = cache.rho; + pv.Ye = cache.ye; + pv.eps = cache.eps; pv.press = cache.press; pv.w_lor = cache.w; - CCTK_REAL Ye_tmp3 = pv.Ye; - pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp3); - CCTK_REAL Ye_tmp4 = pv.Ye; - pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, Ye_tmp4); + pv.temperature = eos_3p->temp_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_3p->kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); // ------------------------------------------------------------------ // Velocity: RePrimAnd magnetic-aware reconstruction // v^i = mu * x * ( r^i + (r.b) * mu * b^i ) // where r^i = S^i / D, b^i = B^i / sqrt(D), (r.b) = (B^i S_i)/(D^(3/2)). // ------------------------------------------------------------------ - if (use_zprim) { - const vec mom_up = calc_contraction(gup, cv.mom); - const CCTK_REAL D = cv.dens; - const CCTK_REAL sqD = std::sqrt(D + CCTK_REAL(0)); + const vec mom_up = calc_contraction(gup, cv.mom); - const vec r_u = mom_up * (CCTK_REAL(1) / (D + CCTK_REAL(1e-300))); - const vec b_u = cv.dBvec * (CCTK_REAL(1) / (sqD + CCTK_REAL(1e-300))); - const CCTK_REAL rb = BiSi / ((D + CCTK_REAL(1e-300)) * (sqD + CCTK_REAL(1e-300))); + const CCTK_REAL D = cv.dens; + const CCTK_REAL sqD = std::sqrt(D); - pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u); - pv.w_lor = cache.w; - } else { - const vec mom_up = calc_contraction(gup, cv.mom); - - const CCTK_REAL D = cv.dens; - const CCTK_REAL sqD = std::sqrt(D + CCTK_REAL(0)); + const vec r_u = mom_up * (1.0 / (D + 1e-300)); + const vec b_u = cv.dBvec * (1.0 / (sqD + 1e-300)); + const CCTK_REAL rb = BiSi / (D * (sqD + 1e-300)); - const vec r_u = mom_up * (CCTK_REAL(1) / (D + CCTK_REAL(1e-300))); - const vec b_u = cv.dBvec * (CCTK_REAL(1) / (sqD + CCTK_REAL(1e-300))); - const CCTK_REAL rb = BiSi / ((D + CCTK_REAL(1e-300)) * (sqD + CCTK_REAL(1e-300))); + pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u); - pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u); - pv.w_lor = cache.w; - } - - vec v_low = calc_contraction(glo, pv.vel); + vec v_low = calc_contraction(glo, pv.vel); CCTK_REAL vsq = calc_contraction(v_low, pv.vel); - if (vsq >= v_lim*v_lim) { + if (vsq >= v_lim * v_lim) { const CCTK_REAL scale = v_lim / (std::sqrt(vsq) + 1e-300); pv.vel *= scale; - vsq = v_lim*v_lim; + vsq = v_lim * v_lim; } - pv.w_lor = CCTK_REAL(1) / std::sqrt(std::max(CCTK_REAL(1e-16), CCTK_REAL(1) - vsq)); + pv.w_lor = CCTK_REAL(1) / std::sqrt(std::max(1e-32, 1.0 - vsq)); pv.Bvec = cv.dBvec; - const vec Elow = calc_cross_product(pv.Bvec, pv.vel); + const vec Elow = calc_cross_product(pv.Bvec, pv.vel); pv.E = calc_contraction(gup, Elow); if (std::abs(result.first - result.second) > - tolerance_0 * std::min(std::abs(result.first), std::abs(result.second))) { + tolerance_0 * + std::min(std::abs(result.first), std::abs(result.second))) { cons_vars cv_check; cv_check.from_prim(pv, glo); - cv_check.dens /= sqrt_detg; - cv_check.tau /= sqrt_detg; - cv_check.mom /= sqrt_detg; + cv_check.dens /= sqrt_detg; + cv_check.tau /= sqrt_detg; + cv_check.mom /= sqrt_detg; cv_check.dBvec /= sqrt_detg; - cv_check.DYe /= sqrt_detg; - cv_check.DEnt /= sqrt_detg; + cv_check.DYe /= sqrt_detg; + cv_check.DEnt /= sqrt_detg; const CCTK_REAL small = 1e-50; - const CCTK_REAL max_error = - sqrt(max({ pow((cv_check.dens - cv.dens )/(cv.dens + small), 2.0), - pow((cv_check.mom(0) - cv.mom(0))/(cv.mom(0) + small), 2.0), - pow((cv_check.mom(1) - cv.mom(1))/(cv.mom(1) + small), 2.0), - pow((cv_check.mom(2) - cv.mom(2))/(cv.mom(2) + small), 2.0), - pow((cv_check.tau - cv.tau )/(cv.tau + small), 2.0) })); + const CCTK_REAL max_error = sqrt( + max({pow((cv_check.dens - cv.dens) / (cv.dens + small), 2.0), + pow((cv_check.mom(0) - cv.mom(0)) / (cv.mom(0) + small), 2.0), + pow((cv_check.mom(1) - cv.mom(1)) / (cv.mom(1) + small), 2.0), + pow((cv_check.mom(2) - cv.mom(2)) / (cv.mom(2) + small), 2.0), + pow((cv_check.tau - cv.tau) / (cv.tau + small), 2.0)})); if (max_error > cons_error) { rep.set_root_conv(); @@ -259,12 +240,9 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - operator()(const EOSType* eos_3p, - prim_vars& pv, cons_vars& cv, - const smat& /*gup_unused*/, - const smat& glo, - c2p_report& rep) const - { + operator()(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, + const smat & /*gup_unused*/, + const smat &glo, c2p_report &rep) const { solve(eos_3p, pv, cv, glo, rep); } }; @@ -272,4 +250,3 @@ public: } // namespace Con2PrimFactory #endif - diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx index 61fc35625..4bd4d5dcf 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx @@ -20,43 +20,48 @@ public: CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE f_upper(CCTK_REAL h0_, CCTK_REAL rsqr_, CCTK_REAL rbsqr_, CCTK_REAL bsqr_) - : h0(h0_), h0sqr(h0_*h0_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_) {} + : h0(h0_), h0sqr(h0_ * h0_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_) {} - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - CCTK_REAL x_from_mu(CCTK_REAL mu) const { return CCTK_REAL(1) / (CCTK_REAL(1) + mu * bsqr); } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL + x_from_mu(CCTK_REAL mu) const { + return 1.0 / (1.0 + mu * bsqr); + } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - CCTK_REAL rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { - return x * (rsqr * x + mu * (x + CCTK_REAL(1)) * rbsqr); + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL + rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + return x * (rsqr * x + mu * (x + 1.0) * rbsqr); } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - CCTK_REAL new_h0w_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL + new_h0w_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { return std::sqrt(h0sqr + rfsqr_from_mu_x(mu, x)); } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - auto operator()(CCTK_REAL mu) const -> std::pair { - const CCTK_REAL x = x_from_mu(mu); + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE auto + operator()(CCTK_REAL mu) const -> std::pair { + const CCTK_REAL x = x_from_mu(mu); const CCTK_REAL xsq = x * x; - const CCTK_REAL hw = new_h0w_from_mu_x(mu, x); - const CCTK_REAL b = x * (xsq * rsqr + mu * (CCTK_REAL(1) + x + xsq) * rbsqr); - const CCTK_REAL f = mu * hw - CCTK_REAL(1); - const CCTK_REAL df = (h0sqr + b) / hw; - return { f, df }; + const CCTK_REAL hw = new_h0w_from_mu_x(mu, x); + const CCTK_REAL b = x * (xsq * rsqr + mu * (1.0 + x + xsq) * rbsqr); + const CCTK_REAL f = mu * hw - 1.0; + const CCTK_REAL df = (h0sqr + b) / hw; + return {f, df}; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - interval initial_bracket() const { - CCTK_REAL mu_min = CCTK_REAL(1) / std::sqrt(h0sqr + rsqr); - const CCTK_REAL mu0 = CCTK_REAL(1) / h0; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval + initial_bracket() const { + CCTK_REAL mu_min = 1.0 / std::sqrt(h0sqr + rsqr); + const CCTK_REAL mu0 = 1.0 / h0; const CCTK_REAL rfs_min = rfsqr_from_mu_x(mu0, x_from_mu(mu0)); - CCTK_REAL mu_max = CCTK_REAL(1) / std::sqrt(h0sqr + rfs_min); + CCTK_REAL mu_max = 1.0 / std::sqrt(h0sqr + rfs_min); const CCTK_REAL margin = 10 * std::numeric_limits::epsilon(); - mu_max *= CCTK_REAL(1) + margin; - mu_min *= CCTK_REAL(1) - margin; - if (mu_max <= mu_min) { mu_min = CCTK_REAL(0); mu_max = mu0 * (CCTK_REAL(1) + margin); } - return { mu_min, mu_max }; + mu_max *= 1.0 + margin; + mu_min *= 1.0 - margin; + if (mu_max <= mu_min) { + mu_min = 0.0; + mu_max = mu0 * (1.0 + margin); + } + return {mu_min, mu_max}; } CCTK_REAL h0, h0sqr, rsqr, rbsqr, bsqr; @@ -64,8 +69,7 @@ public: // --------------------------- froot ----------------------------------- -template -class froot { +template class froot { public: using value_t = CCTK_REAL; @@ -80,75 +84,77 @@ public: }; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - froot(const EOSType* eos_, CCTK_REAL valid_ye, - CCTK_REAL d_, CCTK_REAL qtot_, CCTK_REAL rsqr_, CCTK_REAL rbsqr_, CCTK_REAL bsqr_, - cache &last_) - : eos(eos_), d(d_), qtot(qtot_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_), - brosqr(rsqr_ * bsqr_ - rbsqr_), last(last_) - { - last.ye = valid_ye; + froot(const EOSType *eos_, CCTK_REAL valid_ye, CCTK_REAL d_, CCTK_REAL qtot_, + CCTK_REAL rsqr_, CCTK_REAL rbsqr_, CCTK_REAL bsqr_, cache &last_) + : eos(eos_), d(d_), qtot(qtot_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_), + brosqr(rsqr_ * bsqr_ - rbsqr_), last(last_) { + last.ye = valid_ye; last.calls = 0; // Build h0 from EOS min bounds const CCTK_REAL rhomin = eos->rgrho.min; CCTK_REAL epsmin = eos->rgeps.min; - CCTK_REAL ye_tmp = valid_ye; // EOS expects non-const lvalue ref for Ye - const CCTK_REAL pmin = eos->press_from_valid_rho_eps_ye(rhomin, epsmin, ye_tmp); - h0 = CCTK_REAL(1) + epsmin + pmin / rhomin; + const CCTK_REAL pmin = + eos->press_from_valid_rho_eps_ye(rhomin, epsmin, valid_ye); + h0 = 1.0 + epsmin + pmin / rhomin; - const CCTK_REAL zsqrinf = rsqr / (h0*h0); - const CCTK_REAL wsqrinf = CCTK_REAL(1) + zsqrinf; - winf = std::sqrt(wsqrinf); + const CCTK_REAL zsqrinf = rsqr / (h0 * h0); + const CCTK_REAL wsqrinf = 1.0 + zsqrinf; + winf = std::sqrt(wsqrinf); vsqrinf = zsqrinf / wsqrinf; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - CCTK_REAL operator()(CCTK_REAL mu) { + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL + operator()(CCTK_REAL mu) { cache &c = last; - c.lmu = mu; - c.x = x_from_mu(mu); + c.lmu = mu; + c.x = x_from_mu(mu); const CCTK_REAL rfsqr = rfsqr_from_mu_x(mu, c.x); - const CCTK_REAL qf = qf_from_mu_x(mu, c.x); + const CCTK_REAL qf = qf_from_mu_x(mu, c.x); c.vsqr = rfsqr * mu * mu; - if (c.vsqr >= vsqrinf) { c.vsqr = vsqrinf; c.w = winf; } - else { c.w = CCTK_REAL(1) / std::sqrt(CCTK_REAL(1) - c.vsqr); } + if (c.vsqr >= vsqrinf) { + c.vsqr = vsqrinf; + c.w = winf; + } else { + c.w = 1.0 / std::sqrt(1.0 - c.vsqr); + } c.rho_raw = d / c.w; - c.rho = fmin(fmax(eos->rgrho.min, c.rho_raw), eos->rgrho.max); + c.rho = fmin(fmax(eos->rgrho.min, c.rho_raw), eos->rgrho.max); c.eps_raw = get_eps_raw(mu, qf, rfsqr, c.w); - c.eps = fmin(fmax(eos->rgeps.min, c.eps_raw), eos->rgeps.max); + c.eps = fmin(fmax(eos->rgeps.min, c.eps_raw), eos->rgeps.max); - CCTK_REAL ye_tmp = c.ye; - c.press = eos->press_from_valid_rho_eps_ye(c.rho, c.eps, ye_tmp); + c.press = eos->press_from_valid_rho_eps_ye(c.rho, c.eps, c.ye); ++c.calls; - const CCTK_REAL a = c.press / (c.rho * (CCTK_REAL(1) + c.eps)); - const CCTK_REAL h = (CCTK_REAL(1) + c.eps) * (CCTK_REAL(1) + a); + const CCTK_REAL a = c.press / (c.rho * (1.0 + c.eps)); + const CCTK_REAL h = (1.0 + c.eps) * (1.0 + a); - const CCTK_REAL hbw_raw = (CCTK_REAL(1) + a) * (CCTK_REAL(1) + qf - mu * rfsqr); - const CCTK_REAL hbw = fmax(hbw_raw, h / c.w); - const CCTK_REAL newmu = CCTK_REAL(1) / (hbw + rfsqr * mu); + const CCTK_REAL hbw_raw = (1.0 + a) * (1.0 + qf - mu * rfsqr); + const CCTK_REAL hbw = fmax(hbw_raw, h / c.w); + const CCTK_REAL newmu = 1.0 / (hbw + rfsqr * mu); return mu - newmu; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - interval initial_bracket() const { - CCTK_REAL mu_max = CCTK_REAL(1) / h0; // fallback + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval + initial_bracket() const { + CCTK_REAL mu_max = 1.0 / h0; // fallback { const int ndigits2 = 36; RePrimAnd::f_upper g(h0, rsqr, rbsqr, bsqr); ROOTSTAT status{}; mu_max = findroot_using_deriv(g, status, ndigits2, ndigits2 + 4); - if (status != ROOTSTAT::SUCCESS || !std::isfinite(mu_max)) mu_max = CCTK_REAL(1) / h0; + if (status != ROOTSTAT::SUCCESS || !std::isfinite(mu_max)) + mu_max = 1.0 / h0; } - const CCTK_REAL margin = std::ldexp(CCTK_REAL(1), 3 - 36); - mu_max *= (CCTK_REAL(1) + margin); - const CCTK_REAL mu_min = CCTK_REAL(0); - return { mu_min, mu_max }; + const CCTK_REAL margin = std::ldexp(1.0, 3 - 36); + mu_max *= (1.0 + margin); + const CCTK_REAL mu_min = 0.0; + return {mu_min, mu_max}; } // Public parameters accessed by rarecase/f_rare @@ -156,102 +162,115 @@ public: CCTK_REAL winf{1}, vsqrinf{0}; private: - const EOSType* eos; + const EOSType *eos; cache &last; - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - CCTK_REAL x_from_mu(CCTK_REAL mu) const { - return CCTK_REAL(1) / (CCTK_REAL(1) + mu * bsqr); + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL + x_from_mu(CCTK_REAL mu) const { + return 1.0 / (1.0 + mu * bsqr); } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - CCTK_REAL rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { - return x * (rsqr * x + mu * (x + CCTK_REAL(1)) * rbsqr); + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL + rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + return x * (rsqr * x + mu * (x + 1.0) * rbsqr); } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - CCTK_REAL qf_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL + qf_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const { const CCTK_REAL mux = mu * x; return qtot - (bsqr + mux * mux * brosqr) * CCTK_REAL(0.5); } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - static CCTK_REAL get_eps_raw(CCTK_REAL mu, CCTK_REAL qf, CCTK_REAL rfsqr, CCTK_REAL w) { - return w * (qf - mu * rfsqr * (CCTK_REAL(1) - mu * w / (CCTK_REAL(1) + w))); + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE static CCTK_REAL + get_eps_raw(CCTK_REAL mu, CCTK_REAL qf, CCTK_REAL rfsqr, CCTK_REAL w) { + return w * (qf - mu * rfsqr * (1.0 - mu * w / (1.0 + w))); } }; // -------------------- f_rare and rarecase ------------------- -template -class f_rare { +template class f_rare { const CCTK_REAL v2targ; const froot &f; + public: using value_t = CCTK_REAL; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE f_rare(CCTK_REAL wtarg, const froot &f_) - : v2targ(CCTK_REAL(1) - CCTK_REAL(1)/(wtarg*wtarg)), f(f_) {} + : v2targ(1.0 - 1.0 / (wtarg * wtarg)), f(f_) {} - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - auto operator()(CCTK_REAL mu) const -> std::pair { - const CCTK_REAL x = CCTK_REAL(1) / (CCTK_REAL(1) + mu * f.bsqr); - const CCTK_REAL xsq = x * x; - const CCTK_REAL rfsq = x * (f.rsqr * x + mu * (x + CCTK_REAL(1)) * f.rbsqr); - const CCTK_REAL vsq = mu * mu * rfsq; - const CCTK_REAL y = vsq - v2targ; - const CCTK_REAL dy = CCTK_REAL(2) * mu * x * (xsq * f.rsqr + mu * (xsq + x + CCTK_REAL(1)) * f.rbsqr); - return { y, dy }; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE auto + operator()(CCTK_REAL mu) const -> std::pair { + const CCTK_REAL x = 1.0 / (1.0 + mu * f.bsqr); + const CCTK_REAL xsq = x * x; + const CCTK_REAL rfsq = x * (f.rsqr * x + mu * (x + 1.0) * f.rbsqr); + const CCTK_REAL vsq = mu * mu * rfsq; + const CCTK_REAL y = vsq - v2targ; + const CCTK_REAL dy = + CCTK_REAL(2) * mu * x * (xsq * f.rsqr + mu * (xsq + x + 1.0) * f.rbsqr); + return {y, dy}; } }; -template -class rarecase { +template class rarecase { public: interval bracket; - bool rho_too_big{false}, rho_big{false}, rho_too_small{false}, rho_small{false}; + bool rho_too_big{false}, rho_big{false}, rho_too_small{false}, + rho_small{false}; - CCTK_HOST CCTK_DEVICE - rarecase(const interval ibracket, - const interval rgrho, - const froot &f) - { + CCTK_HOST CCTK_DEVICE rarecase(const interval ibracket, + const interval rgrho, + const froot &f) { CCTK_REAL muc0 = ibracket.min(); CCTK_REAL muc1 = ibracket.max(); const int ndigits = 30; if (f.d > rgrho.max()) { const CCTK_REAL wc = f.d / rgrho.max(); - if (wc > f.winf) { rho_too_big = true; } - else { + if (wc > f.winf) { + rho_too_big = true; + } else { f_rare g(wc, f); - if (g(muc1).first <= CCTK_REAL(0)) { rho_too_big = true; } - else if (g(muc0).first < CCTK_REAL(0)) { + if (g(muc1).first <= 0.0) { + rho_too_big = true; + } else if (g(muc0).first < 0.0) { ROOTSTAT status{}; - CCTK_REAL mucc = findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); - if (status == ROOTSTAT::SUCCESS) { muc0 = fmax(muc0, mucc); rho_big = true; } - else { rho_too_big = true; } + CCTK_REAL mucc = + findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); + if (status == ROOTSTAT::SUCCESS) { + muc0 = fmax(muc0, mucc); + rho_big = true; + } else { + rho_too_big = true; + } } } } if (f.d < rgrho.min()) { const CCTK_REAL wc = f.d / rgrho.min(); - if (wc < CCTK_REAL(1)) { rho_too_small = true; } - else { + if (wc < 1.0) { + rho_too_small = true; + } else { f_rare g(wc, f); - if (g(muc0).first >= CCTK_REAL(0)) { rho_too_small = true; } - else if (g(muc1).first > CCTK_REAL(0)) { + if (g(muc0).first >= 0.0) { + rho_too_small = true; + } else if (g(muc1).first > 0.0) { ROOTSTAT status{}; - CCTK_REAL mucc = findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); - if (status == ROOTSTAT::SUCCESS) { muc1 = fmin(muc1, mucc); rho_small = true; } - else { rho_too_small = true; } + CCTK_REAL mucc = + findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2); + if (status == ROOTSTAT::SUCCESS) { + muc1 = fmin(muc1, mucc); + rho_small = true; + } else { + rho_too_small = true; + } } } } - bracket = interval{ muc0, muc1 }; + bracket = interval{muc0, muc1}; } }; @@ -259,4 +278,3 @@ public: } // namespace Con2PrimFactory #endif - diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx index 7706bb58f..ee516a59b 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx @@ -3,37 +3,38 @@ namespace Con2PrimFactory { -template -class interval { +template class interval { T lo_, hi_; + public: using value_t = T; - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - interval() : lo_(T(0)), hi_(T(0)) {} + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval() + : lo_(T(0)), hi_(T(0)) {} - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - interval(T lo, T hi) : lo_(lo), hi_(hi) {} + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval(T lo, T hi) + : lo_(lo), hi_(hi) {} - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - T min() const { return lo_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T min() const { + return lo_; + } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - T max() const { return hi_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T max() const { + return hi_; + } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - T& min() { return lo_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T &min() { return lo_; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - T& max() { return hi_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T &max() { return hi_; } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - T width() const { return hi_ - lo_; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T width() const { + return hi_ - lo_; + } - CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE - T mid() const { return (lo_ + hi_) * T(0.5); } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T mid() const { + return (lo_ + hi_) * T(0.5); + } }; } // namespace Con2PrimFactory #endif - diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx index 8ddc2374e..1737141e3 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx @@ -6,25 +6,26 @@ #include #include -#include "roots.hxx" // Algo::brent +#include "roots.hxx" // Algo::brent #include "c2p_1DRePrimAnd_intervals.hxx" #include "c2p_utils.hxx" namespace Con2PrimFactory { // Derivative-assisted adapter (we still just use Brent on the residual) -template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline -T findroot_using_deriv(F& f, interval bracket, ROOTSTAT& errs, - int digits, unsigned int max_calls = 64) -{ +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T +findroot_using_deriv(F &f, interval bracket, ROOTSTAT &errs, int digits, + unsigned int max_calls = 256) { const int minbits = digits; int iters = 0; auto fn = [&](T x) -> T { auto val = f(x); - if constexpr (std::is_same>::value) return val.first; - else return val; + if constexpr (std::is_same >::value) + return val.first; + else + return val; }; const T fa = fn(bracket.min()), fb = fn(bracket.max()); @@ -33,49 +34,57 @@ T findroot_using_deriv(F& f, interval bracket, ROOTSTAT& errs, return std::numeric_limits::quiet_NaN(); } - auto ab = Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); - errs = (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + auto ab = + Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); + errs = + (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; const T a_root = ab.first, b_root = ab.second; const T fa2 = fn(a_root), fb2 = fn(b_root); - if (fb2 == T(0) || std::abs(fb2) < std::abs(fa2)) return b_root; - if (std::abs(fa2) < std::abs(fb2)) return a_root; + if (fb2 == T(0) || std::abs(fb2) < std::abs(fa2)) + return b_root; + if (std::abs(fa2) < std::abs(fb2)) + return a_root; return (a_root + b_root) * T(0.5); } // Derivative-free Brent adapter -template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline -interval findroot_no_deriv(F& f, interval bracket, T /*acc*/, - unsigned int max_calls, ROOTSTAT& errs) -{ +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval +findroot_no_deriv(F &f, interval bracket, T /*acc*/, unsigned int max_calls, + ROOTSTAT &errs) { const int minbits = std::numeric_limits::digits - 4; int iters = 0; auto fn = [&](T x) -> T { return f(x); }; const T fa = fn(bracket.min()), fb = fn(bracket.max()); - if (fa * fb > T(0)) { errs = ROOTSTAT::NOT_BRACKETED; return { bracket.min(), bracket.max() }; } + if (fa * fb > T(0)) { + errs = ROOTSTAT::NOT_BRACKETED; + return {bracket.min(), bracket.max()}; + } - auto ab = Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); - errs = (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + auto ab = + Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); + errs = + (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; - return { ab.first, ab.second }; + return {ab.first, ab.second}; } -template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline -T findroot_using_deriv(F& f, ROOTSTAT& errs, int digits, unsigned int max_calls = 64) { +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T +findroot_using_deriv(F &f, ROOTSTAT &errs, int digits, + unsigned int max_calls = 256) { return findroot_using_deriv(f, f.initial_bracket(), errs, digits, max_calls); } -template -CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline -interval findroot_no_deriv(F& f, T tol, unsigned int max_calls, ROOTSTAT& errs) { +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval +findroot_no_deriv(F &f, T tol, unsigned int max_calls, ROOTSTAT &errs) { return findroot_no_deriv(f, f.initial_bracket(), tol, max_calls, errs); } } // namespace Con2PrimFactory #endif - From 152641fdcf764d42c47e6a3b4da4c87d400c999b Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 4 Mar 2026 01:48:44 -0600 Subject: [PATCH 04/10] AsterX: add soft root convergence parameters to param.ccl --- AsterX/param.ccl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/AsterX/param.ccl b/AsterX/param.ccl index 30fb5f8ae..11dfb461f 100644 --- a/AsterX/param.ccl +++ b/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 From e84ac674383010d0fc25cc85d01a34c5f6fdf1aa Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 4 Mar 2026 01:51:47 -0600 Subject: [PATCH 05/10] AsterX: add soft root convergence parameters in c2p object definitions --- AsterX/src/con2prim.cxx | 91 ++++++++++++++++++++++++++++++----------- 1 file changed, 68 insertions(+), 23 deletions(-) diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx index aaf83e716..886e55be0 100644 --- a/AsterX/src/con2prim.cxx +++ b/AsterX/src/con2prim.cxx @@ -47,7 +47,7 @@ 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; + 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")) { @@ -61,7 +61,7 @@ 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; + 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")) { @@ -146,25 +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, - cons_error_limit, vw_lim, B_lim, rho_BH, eps_BH, - vwlim_BH, Ye_lenient, use_z, use_temperature, - use_press_atmo); + 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); // ---------- @@ -282,9 +286,9 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p, break; } case c2p_first_t::RePrimAnd: { - c2p_RPA.solve(eos_3p, pv, cv, glo, rep_first); + 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; @@ -317,10 +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, glo, rep_second); + 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; @@ -506,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; @@ -590,21 +606,50 @@ extern "C" void AsterX_Con2Prim_Interpolate_Failed(CCTK_ARGUMENTS) { const vec vely_nbs = get_neighbors(vely, p); const vec velz_nbs = get_neighbors(velz, p); const vec eps_nbs = get_neighbors(eps, p); + const vec press_nbs = get_neighbors(press, p); const vec saved_rho_nbs = get_neighbors(saved_rho, p); const vec saved_velx_nbs = get_neighbors(saved_velx, p); const vec saved_vely_nbs = get_neighbors(saved_vely, p); const vec saved_velz_nbs = get_neighbors(saved_velz, p); const vec 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; From 5cc7103161884769ee0ea54eda5fc9c1c0d6d1fc Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 4 Mar 2026 01:59:25 -0600 Subject: [PATCH 06/10] AsterX: clang-format --- AsterX/src/estimate_error.cxx | 2 +- AsterX/src/fluxes.cxx | 21 +++++++++++---------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/AsterX/src/estimate_error.cxx b/AsterX/src/estimate_error.cxx index d8e3d49f3..968097e77 100644 --- a/AsterX/src/estimate_error.cxx +++ b/AsterX/src/estimate_error.cxx @@ -21,7 +21,7 @@ enum class regrid_t { }; regrid_t regridMethod; vector regridVarsI; -vector> indextypeRegridVars; +vector > indextypeRegridVars; extern "C" void AsterX_EstimateError_Setup(CCTK_ARGUMENTS) { DECLARE_CCTK_PARAMETERS; diff --git a/AsterX/src/fluxes.cxx b/AsterX/src/fluxes.cxx index 113b703a3..b337bd4ab 100644 --- a/AsterX/src/fluxes.cxx +++ b/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 &var, const PointDesc &p, - bool gf_is_rho, bool gf_is_press) { - return reconstruct>(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 &var, + const PointDesc &p, + bool gf_is_rho, + bool gf_is_press) { + return reconstruct >(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 &var, const PointDesc &p, bool gf_is_rho, bool gf_is_press) { - return reconstruct>(var, p, reconstruction_LO, dir_i, - gf_is_rho, gf_is_press, press, - gf_vel_dir_i, reconstruct_params); + return reconstruct >( + var, p, reconstruction_LO, dir_i, gf_is_rho, gf_is_press, press, + gf_vel_dir_i, reconstruct_params); }; const auto calcflux = From e39b988fcf82030a6abf21b055026db66184023a Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 4 Mar 2026 02:00:39 -0600 Subject: [PATCH 07/10] Con2PrimFactory: add params for soft root convergence --- Con2PrimFactory/param.ccl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/Con2PrimFactory/param.ccl b/Con2PrimFactory/param.ccl index 1e52cba5e..cf35d95ab 100644 --- a/Con2PrimFactory/param.ccl +++ b/Con2PrimFactory/param.ccl @@ -54,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" @@ -73,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" From 5370de687a8f91e9074c51fedb701e717996fa1e Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 4 Mar 2026 02:01:40 -0600 Subject: [PATCH 08/10] Con2PrimFactory: add soft root convergence mechanism, update rootstat diagnostic, add vwlim condition in Palenzuela C2P --- Con2PrimFactory/src/c2p_1DEntropy.hxx | 51 ++++++++--- Con2PrimFactory/src/c2p_1DPalenzuela.hxx | 104 +++++++++++++++-------- Con2PrimFactory/src/c2p_2DNoble.hxx | 31 +++++-- Con2PrimFactory/src/c2p_report.hxx | 12 +++ 4 files changed, 143 insertions(+), 55 deletions(-) diff --git a/Con2PrimFactory/src/c2p_1DEntropy.hxx b/Con2PrimFactory/src/c2p_1DEntropy.hxx index 1a2b46b3a..adbd5ac1c 100644 --- a/Con2PrimFactory/src/c2p_1DEntropy.hxx +++ b/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 &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 &beta, const smat &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::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(1.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 diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx index 66d64f100..05039584c 100644 --- a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx +++ b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx @@ -18,7 +18,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 &mom, @@ -68,7 +69,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; @@ -88,6 +90,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; @@ -311,10 +315,11 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, const CCTK_REAL alp, const vec &beta, const smat &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 */ @@ -420,34 +425,30 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, //} } - auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters); - - // hybrid: prefer endpoint with smaller |f|, else midpoint - CCTK_REAL a_root = result.first; - CCTK_REAL b_root = result.second; - CCTK_REAL fa = fn(a_root); - CCTK_REAL fb = fn(b_root); - - CCTK_REAL xPalenzuela_Sol; - if (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) { - // exact root or smaller residual at b - xPalenzuela_Sol = b_root; - } else if (std::abs(fa) < std::abs(fb)) { - // smaller residual at a - xPalenzuela_Sol = a_root; - } else { - // fall back to midpoint - xPalenzuela_Sol = CCTK_REAL(0.5) * (a_root + b_root); + 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; } - // for now, we do not use the above hybrid scheme - xPalenzuela_Sol = 0.5 * (result.first + result.second); + auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters); - // if (abs(fn(result.first)) < abs(fn(result.second))) { - // xPalenzuela_Sol = result.first; - //} else { - // xPalenzuela_Sol = result.second; - //} + // Legacy endpoint-preference selector kept for reference; below we use the + // midpoint rule for xPalenzuela_Sol. + // CCTK_REAL a_root = result.first; + // CCTK_REAL b_root = result.second; + // CCTK_REAL fa = fn(a_root); + // CCTK_REAL fb = fn(b_root); + // CCTK_REAL xPalenzuela_Sol; + // if (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) { + // xPalenzuela_Sol = b_root; + // } else if (std::abs(fa) < std::abs(fb)) { + // xPalenzuela_Sol = a_root; + // } else { + // xPalenzuela_Sol = CCTK_REAL(0.5) * (a_root + b_root); + // } + + CCTK_REAL xPalenzuela_Sol = CCTK_REAL(0.5) * (result.first + result.second); xPalenzuelaToPrim(xPalenzuela_Sol, Ssq, Bsq, BiSi, eos_3p, pv, cv, gup, glo); @@ -483,14 +484,47 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, // TODO: have an explicit check on max_iters, e.g.: // if (rep.iters >= maxiters || abs(fn(xPalenzuela_Sol)) > tolerance) { - if (abs(result.first - result.second) > - tolerance_0 * min(abs(result.first), abs(result.second))) { + 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(1.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 status to root not converged - rep.set_root_conv(); - cv = cv_const; - // status = ROOTSTAT::NOT_CONVERGED; - return; + const vec v_low_pre = calc_contraction(glo, pv.vel); + const CCTK_REAL vsq_pre = calc_contraction(v_low_pre, pv.vel); + const CCTK_REAL sol_v = std::sqrt(std::max(CCTK_REAL(0), vsq_pre)); + if (sol_v > v_lim) { + pv.rho = cv.dens / w_lim; + pv.vel *= v_lim / (sol_v + CCTK_REAL(1e-300)); + pv.w_lor = w_lim; + const auto rgeps_lim = eos_3p->range_eps_from_rho_ye(pv.rho, pv.Ye); + pv.eps = fmin(fmax(rgeps_lim.min, pv.eps), rgeps_lim.max); + pv.press = eos_3p->press_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + rep.adjust_cons = true; } // set to atmo if computed rho is below floor density diff --git a/Con2PrimFactory/src/c2p_2DNoble.hxx b/Con2PrimFactory/src/c2p_2DNoble.hxx index 7de979257..d3397dc2e 100644 --- a/Con2PrimFactory/src/c2p_2DNoble.hxx +++ b/Con2PrimFactory/src/c2p_2DNoble.hxx @@ -20,7 +20,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 &mom, @@ -78,7 +79,8 @@ CCTK_HOST 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; @@ -98,6 +100,8 @@ CCTK_HOST 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; @@ -286,6 +290,7 @@ c2p_2DNoble::solve(const EOSType *eos_3p, prim_vars &pv, prim_vars &pv_seeds, 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 */ @@ -522,18 +527,26 @@ c2p_2DNoble::solve(const EOSType *eos_3p, prim_vars &pv, prim_vars &pv_seeds, //} if (fabs(errx) > tolerance) { - // set status to root not converged - rep.set_root_conv(); - // status = ROOTSTAT::NOT_CONVERGED; - cv = cv_const; - return; + bool accept_soft = false; + if (soft_root_convergence && std::isfinite(errx)) { + const CCTK_REAL soft_tol = soft_root_width_factor * tolerance; + accept_soft = (fabs(errx) <= soft_tol); + } + if (!accept_soft) { + // set status to root not converged + rep.set_root_conv(); + // status = ROOTSTAT::NOT_CONVERGED; + cv = cv_const; + return; + } + rep.set_soft_root_conv(); } // Check for bad untrapped divergences if ((!isfinite(f)) || (!isfinite(df))) { rep.set_root_bracket(); - // TODO: currently, divergences in f and df are marked as failed bracketing. - // Change this. + cv = cv_const; + return; } /* Calculate primitives from Z and W */ diff --git a/Con2PrimFactory/src/c2p_report.hxx b/Con2PrimFactory/src/c2p_report.hxx index dfc2905ed..6c6824b5d 100644 --- a/Con2PrimFactory/src/c2p_report.hxx +++ b/Con2PrimFactory/src/c2p_report.hxx @@ -65,6 +65,9 @@ public: /// Whether the artificial atmosphere was enforced. bool set_atmo{false}; + /// Root solve accepted via near-convergence policy (warning-only). + bool soft_root_conv{false}; + /// Number of calls to the EOS needed for the root finding. CCTK_INT iters; @@ -205,6 +208,12 @@ public: adjust_cons = true; } + CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + set_soft_root_conv() { + status = SUCCESS; + soft_root_conv = true; + } + CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void debug_message() const { switch (status) { @@ -216,6 +225,9 @@ public: if (adjust_cons) { printf("Conserved variables have been changed. \n"); } + if (soft_root_conv) { + printf("Root solve accepted via near-convergence policy. \n"); + } break; case INVALID_DETG: printf("Invalid metric determinant: detg = %16.8e \n", detg); From bd6df87fabe6b21e4b0ff3dc9dc9372fcce56fdd Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Wed, 4 Mar 2026 02:06:12 -0600 Subject: [PATCH 09/10] Con2PrimFactory: add modifications to RPA; add toms748 solver --- Con2PrimFactory/src/c2p.hxx | 5 + Con2PrimFactory/src/c2p_1DRePrimAnd.hxx | 216 +++++++++----- .../src/c2p_1DRePrimAnd_internals.hxx | 63 ++-- .../src/c2p_1DRePrimAnd_rootfinder.hxx | 206 +++++++++++-- Con2PrimFactory/src/c2p_toms748.hxx | 271 ++++++++++++++++++ Con2PrimFactory/src/test.cxx | 50 +++- 6 files changed, 696 insertions(+), 115 deletions(-) create mode 100644 Con2PrimFactory/src/c2p_toms748.hxx diff --git a/Con2PrimFactory/src/c2p.hxx b/Con2PrimFactory/src/c2p.hxx index c14e9e7d0..e92a40816 100644 --- a/Con2PrimFactory/src/c2p.hxx +++ b/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 &mom, @@ -447,9 +449,12 @@ c2p::cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv, // Compute Bsq vec 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; } diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx index 21f1b6056..50bbeebda 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx @@ -23,16 +23,17 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DRePrimAnd( const EOSType *eos_3p, const atmosphere &atm_in, CCTK_INT maxIter, - CCTK_REAL tol, CCTK_REAL alp_thresh_in, CCTK_REAL consError, - CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, - CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, bool ye_len, bool use_z, - bool use_temperature, bool use_pressure_atmo) { + 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 soft_root_conv, + CCTK_REAL soft_root_width_factor_in) { atmo = atm_in; maxIterations = maxIter; tolerance = tol; alp_thresh = alp_thresh_in; - cons_error = consError; vw_lim = vwlim; w_lim = sqrt(1.0 + vw_lim * vw_lim); v_lim = vw_lim / w_lim; @@ -44,6 +45,10 @@ public: 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); + sigma_max = sigma_max_in; + inv_beta_max = inv_beta_max_in; GammaIdealFluid = eos_3p->gamma; } @@ -57,10 +62,12 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, + const CCTK_REAL alp, const vec &beta, const smat &glo, c2p_report &rep) const { rep.iters = 0; rep.adjust_cons = false; rep.set_atmo = false; + rep.soft_root_conv = false; rep.status = c2p_report::SUCCESS; const CCTK_REAL spatial_detg = calc_det(glo); @@ -85,6 +92,20 @@ public: cv.DYe /= sqrt_detg; cv.DEnt /= sqrt_detg; + if ((!isfinite(cv.dens)) || (cv.dens <= CCTK_REAL(0))) { + rep.set_nans_in_cons(cv.dens, CCTK_REAL(0), CCTK_REAL(0), CCTK_REAL(0), + cv.DYe); + set_to_nan(pv, cv); + return; + } + + if (cv.dens <= atmo.rho_cut) { + rep.set_atmo_set(); + pv.Bvec = cv.dBvec; + atmo.set(pv, cv, glo); + return; + } + const CCTK_REAL Ssq = calc_contraction(calc_contraction(gup, cv.mom), cv.mom); const CCTK_REAL Bsq = @@ -107,48 +128,108 @@ public: // Ye needs to be a non-const lvalue for EOS calls CCTK_REAL Ye_raw = cv.DYe / cv.dens; - CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max); + const CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max); + const bool ye_clipped = + (Ye_raw < eos_3p->rgye.min) || (Ye_raw > eos_3p->rgye.max); typename RePrimAnd::froot::cache cache{}; RePrimAnd::froot f( eos_3p, Ye, cv.dens, cv.tau / cv.dens, Ssq / (cv.dens * cv.dens), (BiSi * BiSi) / (cv.dens * cv.dens * cv.dens), Bsq / cv.dens, cache); - interval mu_br = f.initial_bracket(); + ROOTSTAT prep_stat{}; + interval mu_br = f.initial_bracket(prep_stat); + if (prep_stat != ROOTSTAT::SUCCESS) { + if (prep_stat == ROOTSTAT::NOT_BRACKETED) + rep.set_prep_root_bracket(); + else + rep.set_prep_root_conv(); + cv = cv_const; + return; + } interval rgrho_iv{eos_3p->rgrho.min, eos_3p->rgrho.max}; RePrimAnd::rarecase rc(mu_br, rgrho_iv, f); + if (rc.rho_too_big) { + rep.set_range_rho(cv.dens, cv.dens); + cv = cv_const; + return; + } + if (rc.rho_too_small) { + rep.set_atmo_set(); + pv.Bvec = cv.dBvec; + atmo.set(pv, cv, glo); + return; + } mu_br = rc.bracket; - const CCTK_REAL log2 = std::log(2.0); - const CCTK_INT minbits = int(std::abs(std::log(tolerance)) / log2); - const CCTK_REAL tolerance_0 = std::ldexp(double(1.0), -minbits); const CCTK_INT maxiters = maxIterations; - auto fn = [&](CCTK_REAL mu) { return f(mu); }; - if (fn(mu_br.min()) * fn(mu_br.max()) > 0) { - const CCTK_REAL qtot = cv.tau / cv.dens; - const CCTK_REAL sPal = Bsq / cv.dens; - CCTK_REAL new_hi = 3.0 + 3.0 * qtot - 1.5 * sPal; - mu_br = interval{mu_br.min(), new_hi}; + ROOTSTAT root_stat{}; + interval result = + findroot_no_deriv(f, mu_br, tolerance, (unsigned int)maxiters, root_stat); + rep.iters = (CCTK_INT)cache.calls; + + if (root_stat == ROOTSTAT::NOT_CONVERGED) { + bool accept_soft = false; + if (soft_root_convergence) { + const CCTK_REAL a = fmin(result.min(), result.max()); + const CCTK_REAL b = fmax(result.min(), result.max()); + if (std::isfinite(a) && std::isfinite(b) && b >= a) { + const CCTK_REAL width = b - a; + const CCTK_REAL scale = + fmax(CCTK_REAL(1.0), fmax(std::abs(a), std::abs(b))); + const CCTK_REAL width_tol = + soft_root_width_factor * tolerance * scale; + const CCTK_REAL mu_ref = + fmax(CCTK_REAL(0.5) * (a + b), CCTK_REAL(1.0e-300)); + const bool stopif_ok = f.stopif(mu_ref, width, tolerance); + accept_soft = (width <= width_tol) || stopif_ok; + } + } + if (!accept_soft) { + rep.set_root_conv(); + cv = cv_const; + return; + } + rep.set_soft_root_conv(); + root_stat = ROOTSTAT::SUCCESS; } - auto result = - Algo::brent(fn, mu_br.min(), mu_br.max(), minbits, maxiters, rep.iters); - - const CCTK_REAL a_root = result.first; - const CCTK_REAL b_root = result.second; - const CCTK_REAL fa = fn(a_root); - const CCTK_REAL fb = fn(b_root); - const CCTK_REAL mu = - (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) ? b_root - : (std::abs(fa) < std::abs(fb)) ? a_root - : CCTK_REAL(0.5) * (a_root + b_root); + if (root_stat != ROOTSTAT::SUCCESS) { + if (rc.rho_big) { + rep.set_range_rho(cv.dens, cv.dens); + cv = cv_const; + return; + } + if (rc.rho_small) { + rep.set_atmo_set(); + pv.Bvec = cv.dBvec; + atmo.set(pv, cv, glo); + return; + } + rep.set_root_bracket(); + cv = cv_const; + return; + } - // ------------------------------------------------------------------ - // IMPORTANT: ensure cache corresponds to the final chosen mu - // ------------------------------------------------------------------ - (void)f(mu); + CCTK_REAL mu = cache.lmu; + const CCTK_REAL mu_lo = fmin(result.min(), result.max()); + const CCTK_REAL mu_hi = fmax(result.min(), result.max()); + + if (!std::isfinite(mu) || mu < mu_lo || mu > mu_hi) { + auto fn = [&](CCTK_REAL x) { return f(x); }; + const CCTK_REAL a_root = result.min(); + const CCTK_REAL b_root = result.max(); + const CCTK_REAL fa = fn(a_root); + const CCTK_REAL fb = fn(b_root); + mu = (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) ? b_root + : (std::abs(fa) < std::abs(fb)) ? a_root + : CCTK_REAL(0.5) * (a_root + b_root); + + // Ensure cache corresponds to the final chosen mu. + (void)f(mu); + } // ------------------------------------------------------------------ // Use the EOS-consistent RePrimAnd cached primitives @@ -159,6 +240,18 @@ public: pv.press = cache.press; pv.w_lor = cache.w; + // If root-cache thermodynamics were clipped to EOS bounds, force + // conservative recomputation for consistency with adjusted primitives. + if (ye_clipped) { + rep.adjust_cons = true; + } + { + const auto rgeps = eos_3p->range_eps_from_rho_ye(pv.rho, pv.Ye); + if (cache.eps_raw < rgeps.min || cache.eps_raw > rgeps.max) { + rep.adjust_cons = true; + } + } + pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); @@ -179,55 +272,32 @@ public: pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u); - vec v_low = calc_contraction(glo, pv.vel); - CCTK_REAL vsq = calc_contraction(v_low, pv.vel); - if (vsq >= v_lim * v_lim) { - const CCTK_REAL scale = v_lim / (std::sqrt(vsq) + 1e-300); - pv.vel *= scale; - vsq = v_lim * v_lim; + const CCTK_REAL sol_v = std::sqrt(std::max(CCTK_REAL(0), cache.vsqr)); + if (sol_v > v_lim) { + pv.rho = cv.dens / w_lim; + pv.vel *= v_lim / (sol_v + CCTK_REAL(1e-300)); + pv.w_lor = w_lim; + const auto rgeps_lim = eos_3p->range_eps_from_rho_ye(pv.rho, pv.Ye); + pv.eps = fmin(fmax(rgeps_lim.min, pv.eps), rgeps_lim.max); + pv.press = eos_3p->press_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + rep.adjust_cons = true; + } else { + pv.w_lor = cache.w; } - pv.w_lor = CCTK_REAL(1) / std::sqrt(std::max(1e-32, 1.0 - vsq)); pv.Bvec = cv.dBvec; const vec Elow = calc_cross_product(pv.Bvec, pv.vel); pv.E = calc_contraction(gup, Elow); - if (std::abs(result.first - result.second) > - tolerance_0 * - std::min(std::abs(result.first), std::abs(result.second))) { - - cons_vars cv_check; - cv_check.from_prim(pv, glo); - - cv_check.dens /= sqrt_detg; - cv_check.tau /= sqrt_detg; - cv_check.mom /= sqrt_detg; - cv_check.dBvec /= sqrt_detg; - cv_check.DYe /= sqrt_detg; - cv_check.DEnt /= sqrt_detg; - - const CCTK_REAL small = 1e-50; - const CCTK_REAL max_error = sqrt( - max({pow((cv_check.dens - cv.dens) / (cv.dens + small), 2.0), - pow((cv_check.mom(0) - cv.mom(0)) / (cv.mom(0) + small), 2.0), - pow((cv_check.mom(1) - cv.mom(1)) / (cv.mom(1) + small), 2.0), - pow((cv_check.mom(2) - cv.mom(2)) / (cv.mom(2) + small), 2.0), - pow((cv_check.tau - cv.tau) / (cv.tau + small), 2.0)})); - - if (max_error > cons_error) { - rep.set_root_conv(); - cv = cv_const; - return; - } - } - if (pv.rho < atmo.rho_cut) { rep.set_atmo_set(); atmo.set(pv, cv, glo); return; } - c2p::prims_floors_and_ceilings(eos_3p, pv, cv, glo, rep); + c2p::prims_floors_and_ceilings(eos_3p, pv, cv, alp, beta, glo, rep); if (rep.adjust_cons) { cv.from_prim(pv, glo); @@ -238,6 +308,16 @@ public: } } + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, + const smat &glo, c2p_report &rep) const { + // Backward-compatible overload for legacy call sites. + const CCTK_REAL alp = 1.0; + const vec beta{0.0, 0.0, 0.0}; + solve(eos_3p, pv, cv, alp, beta, glo, rep); + } + template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void operator()(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx index d51702d9d..a00e0f9db 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx @@ -77,6 +77,7 @@ public: CCTK_REAL ye{0}, lmu{0}, x{0}; CCTK_REAL rho{0}, rho_raw{0}; CCTK_REAL eps{0}, eps_raw{0}; + CCTK_REAL etot{0}, etot_raw{0}; CCTK_REAL press{0}; CCTK_REAL vsqr{0}; CCTK_REAL w{1}; @@ -93,10 +94,12 @@ public: // Build h0 from EOS min bounds const CCTK_REAL rhomin = eos->rgrho.min; - CCTK_REAL epsmin = eos->rgeps.min; + const auto rgeps0 = eos->range_eps_from_rho_ye(rhomin, valid_ye); + CCTK_REAL epsmin = rgeps0.min; const CCTK_REAL pmin = eos->press_from_rho_eps_ye(rhomin, epsmin, valid_ye); h0 = 1.0 + epsmin + pmin / rhomin; + h0 = fmax(h0, CCTK_REAL(1.0) + std::numeric_limits::epsilon()); const CCTK_REAL zsqrinf = rsqr / (h0 * h0); const CCTK_REAL wsqrinf = 1.0 + zsqrinf; @@ -125,15 +128,22 @@ public: c.rho = fmin(fmax(eos->rgrho.min, c.rho_raw), eos->rgrho.max); c.eps_raw = get_eps_raw(mu, qf, rfsqr, c.w); - c.eps = fmin(fmax(eos->rgeps.min, c.eps_raw), eos->rgeps.max); + const auto rgeps = eos->range_eps_from_rho_ye(c.rho, c.ye); + c.eps = fmin(fmax(rgeps.min, c.eps_raw), rgeps.max); + + c.etot_raw = c.rho_raw * (1.0 + c.eps_raw); + c.etot = c.rho * (1.0 + c.eps); c.press = eos->press_from_rho_eps_ye(c.rho, c.eps, c.ye); ++c.calls; - const CCTK_REAL a = c.press / (c.rho * (1.0 + c.eps)); - const CCTK_REAL h = (1.0 + c.eps) * (1.0 + a); + const CCTK_REAL inv_etot = + CCTK_REAL(1.0) / + fmax(c.etot, std::numeric_limits::min()); + const CCTK_REAL one_plus_a = (c.etot + c.press) * inv_etot; + const CCTK_REAL h = (c.etot + c.press) / c.rho; - const CCTK_REAL hbw_raw = (1.0 + a) * (1.0 + qf - mu * rfsqr); + const CCTK_REAL hbw_raw = one_plus_a * (1.0 + qf - mu * rfsqr); const CCTK_REAL hbw = fmax(hbw_raw, h / c.w); const CCTK_REAL newmu = 1.0 / (hbw + rfsqr * mu); @@ -141,20 +151,39 @@ public: } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval - initial_bracket() const { - CCTK_REAL mu_max = 1.0 / h0; // fallback - { + initial_bracket(ROOTSTAT &status) const { + status = ROOTSTAT::SUCCESS; + CCTK_REAL mu_max = CCTK_REAL(1) / h0; + + if (rsqr >= h0 * h0) { const int ndigits2 = 36; + const CCTK_REAL margin = std::ldexp(CCTK_REAL(1), 3 - ndigits2); RePrimAnd::f_upper g(h0, rsqr, rbsqr, bsqr); - ROOTSTAT status{}; - mu_max = findroot_using_deriv(g, status, ndigits2, ndigits2 + 4); - if (status != ROOTSTAT::SUCCESS || !std::isfinite(mu_max)) - mu_max = 1.0 / h0; + const CCTK_REAL mu_try = + findroot_using_deriv(g, status, ndigits2, ndigits2 + 4); + + if (status != ROOTSTAT::SUCCESS || !std::isfinite(mu_try)) { + if (status == ROOTSTAT::SUCCESS) + status = ROOTSTAT::NOT_CONVERGED; + return {CCTK_REAL(0), CCTK_REAL(1) / h0}; + } + + mu_max = mu_try * (CCTK_REAL(1) + margin); } - const CCTK_REAL margin = std::ldexp(1.0, 3 - 36); - mu_max *= (1.0 + margin); - const CCTK_REAL mu_min = 0.0; - return {mu_min, mu_max}; + + return {CCTK_REAL(0), mu_max}; + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval + initial_bracket() const { + ROOTSTAT status{}; + return initial_bracket(status); + } + + // RePrimAnd stop condition used by the TOMS748-style bracket solver. + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE bool + stopif(CCTK_REAL mu, CCTK_REAL dmu, CCTK_REAL acc) const { + return std::abs(dmu) * last.w * last.w < mu * acc; } // Public parameters accessed by rarecase/f_rare @@ -248,7 +277,7 @@ public: } } - if (f.d < rgrho.min()) { + if (f.d < f.winf * rgrho.min()) { const CCTK_REAL wc = f.d / rgrho.min(); if (wc < 1.0) { rho_too_small = true; diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx index 1737141e3..017286415 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx @@ -1,6 +1,7 @@ #ifndef C2P_1DREPRIMAND_ROOTFINDER_HXX #define C2P_1DREPRIMAND_ROOTFINDER_HXX +#include #include #include #include @@ -8,39 +9,164 @@ #include "roots.hxx" // Algo::brent #include "c2p_1DRePrimAnd_intervals.hxx" +#include "c2p_toms748.hxx" #include "c2p_utils.hxx" namespace Con2PrimFactory { -// Derivative-assisted adapter (we still just use Brent on the residual) +namespace RePrimAnd_rootdetail { + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool same_sign(T a, + T b) { + return (a > T(0) && b > T(0)) || (a < T(0) && b < T(0)); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool +in_open_interval(T x, T a, T b) { + return (x > a) && (x < b); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T midpoint(T a, + T b) { + return (a + b) * T(0.5); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T secant(T a, T fa, + T b, T fb) { + const T denom = (fb - fa); + if (denom == T(0)) + return std::numeric_limits::quiet_NaN(); + return (a * fb - b * fa) / denom; +} + +} // namespace RePrimAnd_rootdetail + +// Derivative-assisted Newton with bracketing, Brent fallback on pathology. template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T findroot_using_deriv(F &f, interval bracket, ROOTSTAT &errs, int digits, unsigned int max_calls = 256) { - const int minbits = digits; - int iters = 0; + using namespace RePrimAnd_rootdetail; - auto fn = [&](T x) -> T { - auto val = f(x); - if constexpr (std::is_same >::value) - return val.first; - else - return val; + if (max_calls < 4u) { + errs = ROOTSTAT::NOT_CONVERGED; + return std::numeric_limits::quiet_NaN(); + } + + T a = bracket.min(); + T b = bracket.max(); + if (b < a) { + const T tmp = a; + a = b; + b = tmp; + } + + auto eval_pair = [&](T x) -> std::pair { + const auto val = f(x); + return {T(val.first), T(val.second)}; }; - const T fa = fn(bracket.min()), fb = fn(bracket.max()); - if (fa * fb > T(0)) { + auto va = eval_pair(a); + auto vb = eval_pair(b); + T fa = va.first; + T fb = vb.first; + + // Match RePrimAnd endpoint precedence: right endpoint first. + if (fb == T(0)) { + errs = ROOTSTAT::SUCCESS; + return b; + } + if (fa == T(0)) { + errs = ROOTSTAT::SUCCESS; + return a; + } + if (!std::isfinite(fa) || !std::isfinite(fb) || same_sign(fa, fb)) { + errs = ROOTSTAT::NOT_BRACKETED; + return std::numeric_limits::quiet_NaN(); + } + + const int minbits = + (digits > 0) ? digits : (std::numeric_limits::digits - 4); + const T tol = std::ldexp(T(1), -minbits); + const unsigned int newton_calls = max_calls - 2u; + + // Same initial guess as RePrimAnd reference. + T x = secant(a, fa, b, fb); + if (!std::isfinite(x) || !in_open_interval(x, a, b)) + x = midpoint(a, b); + + int stalled = 0; + for (unsigned int calls = 0; calls < newton_calls; ++calls) { + const auto vx = eval_pair(x); + const T fx = vx.first; + const T dfx = vx.second; + + if (!std::isfinite(fx) || !std::isfinite(dfx)) + break; + if (fx == T(0)) { + errs = ROOTSTAT::SUCCESS; + return x; + } + + if (same_sign(fa, fx)) { + a = x; + fa = fx; + } else { + b = x; + fb = fx; + } + + const T width = b - a; + const T scale = std::max(T(1), std::max(std::abs(a), std::abs(b))); + if (std::abs(width) <= tol * scale) { + errs = ROOTSTAT::SUCCESS; + return x; + } + + T x_new = std::numeric_limits::quiet_NaN(); + const T dfx_floor = std::sqrt(std::numeric_limits::epsilon()); + if (std::abs(dfx) > dfx_floor) + x_new = x - fx / dfx; + if (!std::isfinite(x_new) || !in_open_interval(x_new, a, b)) + x_new = secant(a, fa, b, fb); + if (!std::isfinite(x_new) || !in_open_interval(x_new, a, b)) + x_new = midpoint(a, b); + + const T step_scale = std::max(T(1), std::abs(x)); + if (std::abs(x_new - x) <= std::numeric_limits::epsilon() * step_scale) + ++stalled; + else + stalled = 0; + x = x_new; + + if (stalled >= 2) + break; + } + + // Fallback to robust bracketing in pathological derivative/stall cases. + auto fn = [&](T x0) -> T { return eval_pair(x0).first; }; + if (!std::isfinite(a) || !std::isfinite(b) || !std::isfinite(fa) || + !std::isfinite(fb) || same_sign(fa, fb)) { errs = ROOTSTAT::NOT_BRACKETED; return std::numeric_limits::quiet_NaN(); } - auto ab = - Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); - errs = - (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + int iters = 0; + auto ab = Algo::brent(fn, a, b, minbits, int(max_calls), iters); + errs = (iters >= int(max_calls)) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + if (errs != ROOTSTAT::SUCCESS) + return std::numeric_limits::quiet_NaN(); const T a_root = ab.first, b_root = ab.second; const T fa2 = fn(a_root), fb2 = fn(b_root); + if (!std::isfinite(fa2) || !std::isfinite(fb2)) { + errs = ROOTSTAT::NOT_CONVERGED; + return std::numeric_limits::quiet_NaN(); + } if (fb2 == T(0) || std::abs(fb2) < std::abs(fa2)) return b_root; @@ -49,26 +175,52 @@ findroot_using_deriv(F &f, interval bracket, ROOTSTAT &errs, int digits, return (a_root + b_root) * T(0.5); } -// Derivative-free Brent adapter +// Derivative-free TOMS748 adapter, Brent fallback on pathological stall. template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval -findroot_no_deriv(F &f, interval bracket, T /*acc*/, unsigned int max_calls, +findroot_no_deriv(F &f, interval bracket, T acc, unsigned int max_calls, ROOTSTAT &errs) { - const int minbits = std::numeric_limits::digits - 4; - int iters = 0; + if (max_calls < 10u) { + errs = ROOTSTAT::NOT_CONVERGED; + T a = bracket.min(); + T b = bracket.max(); + if (b < a) { + const T tmp = a; + a = b; + b = tmp; + } + return {a, b}; + } - auto fn = [&](T x) -> T { return f(x); }; + interval res = toms748_solve(f, bracket, acc, max_calls, errs); + if (errs == ROOTSTAT::SUCCESS || errs == ROOTSTAT::NOT_BRACKETED) + return res; - const T fa = fn(bracket.min()), fb = fn(bracket.max()); - if (fa * fb > T(0)) { + T a = res.min(); + T b = res.max(); + if (b < a) { + const T tmp = a; + a = b; + b = tmp; + } + if (!std::isfinite(a) || !std::isfinite(b)) + return res; + + auto fn = [&](T x) -> T { return f(x); }; + const T fa = fn(a), fb = fn(b); + if (fb == T(0)) + return {b, b}; + if (fa == T(0)) + return {a, a}; + if (!std::isfinite(fa) || !std::isfinite(fb) || fa * fb > T(0)) { errs = ROOTSTAT::NOT_BRACKETED; - return {bracket.min(), bracket.max()}; + return {std::numeric_limits::lowest(), std::numeric_limits::max()}; } - auto ab = - Algo::brent(fn, bracket.min(), bracket.max(), minbits, max_calls, iters); - errs = - (iters >= (int)max_calls) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; + const int minbits = std::numeric_limits::digits - 4; + int iters = 0; + auto ab = Algo::brent(fn, a, b, minbits, int(max_calls), iters); + errs = (iters >= int(max_calls)) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS; return {ab.first, ab.second}; } diff --git a/Con2PrimFactory/src/c2p_toms748.hxx b/Con2PrimFactory/src/c2p_toms748.hxx new file mode 100644 index 000000000..114b4d6f8 --- /dev/null +++ b/Con2PrimFactory/src/c2p_toms748.hxx @@ -0,0 +1,271 @@ +#ifndef C2P_TOMS748_HXX +#define C2P_TOMS748_HXX + +#include +#include +#include +#include +#include + +#include "c2p_1DRePrimAnd_intervals.hxx" +#include "c2p_utils.hxx" + +namespace Con2PrimFactory { + +namespace toms748_detail { + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool same_sign(T a, + T b) { + return (a > T(0) && b > T(0)) || (a < T(0) && b < T(0)); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool +in_open_interval(T x, T a, T b) { + return (x > a) && (x < b); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T midpoint(T a, + T b) { + return (a + b) * T(0.5); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T +secant_step(T a, T fa, T b, T fb) { + const T denom = fb - fa; + if (denom == T(0)) + return std::numeric_limits::quiet_NaN(); + return (a * fb - b * fa) / denom; +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T +iqi_step(T a, T fa, T b, T fb, T c, T fc) { + const T dab = (fa - fb); + const T dac = (fa - fc); + const T dba = (fb - fa); + const T dbc = (fb - fc); + const T dca = (fc - fa); + const T dcb = (fc - fb); + + if (dab == T(0) || dac == T(0) || dba == T(0) || dbc == T(0) || + dca == T(0) || dcb == T(0)) + return std::numeric_limits::quiet_NaN(); + + return a * fb * fc / (dab * dac) + b * fa * fc / (dba * dbc) + + c * fa * fb / (dca * dcb); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval +invalid_interval() { + return {std::numeric_limits::lowest(), std::numeric_limits::max()}; +} + +template struct has_stopif { +private: + template + static auto test(int) + -> decltype(std::declval().stopif(std::declval(), + std::declval(), + std::declval()), + std::true_type{}); + template static std::false_type test(...); + +public: + static constexpr bool value = decltype(test(0))::value; +}; + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool +call_stopif_impl(F &f, T left, T right, T tol, std::true_type) { + // Preserve RePrimAnd contract: stopif(left, right-left, tol) + return f.stopif(left, right - left, tol); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool +call_stopif_impl(F &, T left, T right, T tol, std::false_type) { + const T width = std::abs(right - left); + const T scale = std::max(T(1), std::max(std::abs(left), std::abs(right))); + return width <= tol * scale; +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool +call_stopif(F &f, T left, T right, T tol) { + using has_stopif_t = std::integral_constant::value>; + return call_stopif_impl(f, left, right, tol, has_stopif_t{}); +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +update_bracket(T x, T fx, T &a, T &fa, T &b, T &fb) { + if (same_sign(fa, fx)) { + a = x; + fa = fx; + } else { + b = x; + fb = fx; + } +} + +} // namespace toms748_detail + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval +toms748_solve(F &f, interval bracket, T tol, unsigned int max_calls, + ROOTSTAT &errs) { + using namespace toms748_detail; + + T a = bracket.min(); + T b = bracket.max(); + if (b < a) { + const T tmp = a; + a = b; + b = tmp; + } + + if (!std::isfinite(a) || !std::isfinite(b)) { + errs = ROOTSTAT::NOT_BRACKETED; + return invalid_interval(); + } + + T fa = f(a); + T fb = f(b); + unsigned int calls = 2u; + + // Match RePrimAnd endpoint precedence: right endpoint first. + if (fb == T(0)) { + errs = ROOTSTAT::SUCCESS; + return {b, b}; + } + if (fa == T(0)) { + errs = ROOTSTAT::SUCCESS; + return {a, a}; + } + + if (!std::isfinite(fa) || !std::isfinite(fb) || same_sign(fa, fb)) { + errs = ROOTSTAT::NOT_BRACKETED; + return invalid_interval(); + } + + if (max_calls <= 2u) { + errs = ROOTSTAT::NOT_CONVERGED; + return {a, b}; + } + + if (call_stopif(f, a, b, tol)) { + errs = ROOTSTAT::SUCCESS; + return {a, b}; + } + + T prev_width = b - a; + + while (calls < max_calls) { + const T a0 = a, b0 = b; + const T fa0 = fa, fb0 = fb; + const T width0 = b - a; + + T x1 = secant_step(a, fa, b, fb); + if (!std::isfinite(x1) || !in_open_interval(x1, a, b)) + x1 = midpoint(a, b); + + T fx1 = f(x1); + ++calls; + + if (!std::isfinite(fx1)) { + if (calls >= max_calls) + break; + x1 = midpoint(a, b); + fx1 = f(x1); + ++calls; + if (!std::isfinite(fx1)) { + errs = ROOTSTAT::NOT_CONVERGED; + return {a, b}; + } + } + + if (fx1 == T(0)) { + errs = ROOTSTAT::SUCCESS; + return {x1, x1}; + } + + update_bracket(x1, fx1, a, fa, b, fb); + if (call_stopif(f, a, b, tol)) { + errs = ROOTSTAT::SUCCESS; + return {a, b}; + } + if (calls >= max_calls) + break; + + T x2 = iqi_step(a0, fa0, b0, fb0, x1, fx1); + if (!std::isfinite(x2) || !in_open_interval(x2, a, b)) + x2 = midpoint(a, b); + + const T width1 = b - a; + if ((x2 - a) < T(0.125) * width1 || (b - x2) < T(0.125) * width1) + x2 = midpoint(a, b); + + T fx2 = f(x2); + ++calls; + + if (!std::isfinite(fx2)) { + if (calls >= max_calls) + break; + x2 = midpoint(a, b); + fx2 = f(x2); + ++calls; + if (!std::isfinite(fx2)) { + errs = ROOTSTAT::NOT_CONVERGED; + return {a, b}; + } + } + + if (fx2 == T(0)) { + errs = ROOTSTAT::SUCCESS; + return {x2, x2}; + } + + update_bracket(x2, fx2, a, fa, b, fb); + if (call_stopif(f, a, b, tol)) { + errs = ROOTSTAT::SUCCESS; + return {a, b}; + } + + const T new_width = b - a; + if (calls < max_calls && + (new_width > T(0.5) * width0 || new_width > T(0.9) * prev_width)) { + const T xm = midpoint(a, b); + const T fxm = f(xm); + ++calls; + + if (!std::isfinite(fxm)) { + errs = ROOTSTAT::NOT_CONVERGED; + return {a, b}; + } + if (fxm == T(0)) { + errs = ROOTSTAT::SUCCESS; + return {xm, xm}; + } + + update_bracket(xm, fxm, a, fa, b, fb); + if (call_stopif(f, a, b, tol)) { + errs = ROOTSTAT::SUCCESS; + return {a, b}; + } + } + + prev_width = b - a; + } + + errs = ROOTSTAT::NOT_CONVERGED; + return {a, b}; +} + +} // namespace Con2PrimFactory + +#endif diff --git a/Con2PrimFactory/src/test.cxx b/Con2PrimFactory/src/test.cxx index b8f87764a..59cb61c67 100644 --- a/Con2PrimFactory/src/test.cxx +++ b/Con2PrimFactory/src/test.cxx @@ -6,6 +6,7 @@ #include "c2p.hxx" #include "c2p_1DEntropy.hxx" #include "c2p_1DPalenzuela.hxx" +#include "c2p_1DRePrimAnd.hxx" #include "c2p_2DNoble.hxx" #include "c2p_utils.hxx" @@ -66,17 +67,21 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { // Ye_lenient, use_z, use_temperature) c2p_2DNoble c2p_Noble(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH, eps_BH, vwlim_BH, sigma_max, inv_beta_max, true, false, - false, false); + false, false, false, 1.0); c2p_1DPalenzuela c2p_Pal(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH, eps_BH, vwlim_BH, sigma_max, inv_beta_max, true, - false, false, false); + false, false, false, false, 1.0); + c2p_1DRePrimAnd c2p_RPA(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH, + eps_BH, vwlim_BH, sigma_max, inv_beta_max, true, + false, false, false, false, 1.0); c2p_1DEntropy c2p_Ent(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH, eps_BH, vwlim_BH, sigma_max, inv_beta_max, true, false, - false, false); + false, false, false, 1.0); // Construct error report object: c2p_report rep_Noble; c2p_report rep_Pal; + c2p_report rep_RPA; c2p_report rep_Ent; // Set primitive seeds @@ -211,6 +216,45 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { */ rep_Pal.debug_message(); + // Testing C2P RePrimAnd + CCTK_VINFO("Testing C2P RePrimAnd..."); + c2p_RPA.solve(eos_3p_ig, pv, cv_all, alp, beta, g, rep_RPA); + + printf("pv_seeds, pv: \n" + "rho: %f, %f \n" + "eps: %f, %f \n" + "Ye: %f, %f \n" + "press: %f, %f \n" + "temperature: %f, %f \n" + "entropy: %f, %f \n" + "velx: %f, %f \n" + "vely: %f, %f \n" + "velz: %f, %f \n" + "Bx: %f, %f \n" + "By: %f, %f \n" + "Bz: %f, %f \n", + pv_seeds.rho, pv.rho, pv_seeds.eps, pv.eps, pv_seeds.Ye, pv.Ye, + pv_seeds.press, pv.press, pv_seeds.temperature, pv.temperature, + pv_seeds.entropy, pv.entropy, pv_seeds.vel(0), pv.vel(0), + pv_seeds.vel(1), pv.vel(1), pv_seeds.vel(2), pv.vel(2), + pv_seeds.Bvec(0), pv.Bvec(0), pv_seeds.Bvec(1), pv.Bvec(1), + pv_seeds.Bvec(2), pv.Bvec(2)); + printf("cv: \n" + "dens: %f \n" + "tau: %f \n" + "momx: %f \n" + "momy: %f \n" + "momz: %f \n" + "DYe: %f \n" + "dBx: %f \n" + "dBy: %f \n" + "dBz: %f \n" + "DEnt: %f \n", + cv_all.dens, cv_all.tau, cv_all.mom(0), cv_all.mom(1), cv_all.mom(2), + cv_all.DYe, cv_all.dBvec(0), cv_all.dBvec(1), cv_all.dBvec(2), + cv_all.DEnt); + rep_RPA.debug_message(); + // Testing C2P Entropy CCTK_VINFO("Testing C2P Entropy..."); c2p_Ent.solve(eos_3p_ig, pv, cv_all, alp, beta, g, rep_Ent); From 6271517a5bafb04363a830707e11444a68bd9981 Mon Sep 17 00:00:00 2001 From: jaykalinani Date: Mon, 6 Apr 2026 18:55:10 -0500 Subject: [PATCH 10/10] Con2PrimFactory: Clean-up based on Michi's PR comments --- Con2PrimFactory/src/c2p_1DEntropy.hxx | 2 +- Con2PrimFactory/src/c2p_1DPalenzuela.hxx | 13 ++++++++++++- Con2PrimFactory/src/c2p_1DRePrimAnd.hxx | 16 +--------------- 3 files changed, 14 insertions(+), 17 deletions(-) diff --git a/Con2PrimFactory/src/c2p_1DEntropy.hxx b/Con2PrimFactory/src/c2p_1DEntropy.hxx index adbd5ac1c..558d15dd9 100644 --- a/Con2PrimFactory/src/c2p_1DEntropy.hxx +++ b/Con2PrimFactory/src/c2p_1DEntropy.hxx @@ -407,7 +407,7 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, bool accept_soft = false; if (soft_root_convergence) { const CCTK_REAL scale = - fmax(CCTK_REAL(1.0), fmax(abs(result.first), abs(result.second))); + 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) && diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx index 05039584c..f1cd37bdd 100644 --- a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx +++ b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx @@ -384,6 +384,13 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, return; } + // Clamp Ye in the local working conservatives so the root solve and prim + // reconstruction both use EOS-valid Ye. + CCTK_REAL Ye_raw = cv.DYe / cv.dens; + const CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max); + const bool ye_clipped = (Ye_raw < eos_3p->rgye.min) || (Ye_raw > eos_3p->rgye.max); + cv.DYe = cv.dens * Ye; + // Find x, this is the recovery process // const CCTK_INT minbits = std::numeric_limits::digits - 4; // const CCTK_INT maxiters = maxIterations; @@ -468,6 +475,10 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, return; } + if (ye_clipped) { + rep.adjust_cons = true; + } + // General comment: // One could think of expressing the following condition // in a way that is "safe" against NaNs and infs. First, this only makes @@ -491,7 +502,7 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv, bool accept_soft = false; if (soft_root_convergence) { const CCTK_REAL scale = - fmax(CCTK_REAL(1.0), fmax(abs(result.first), abs(result.second))); + 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) && diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx index 50bbeebda..d791b7a37 100644 --- a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx +++ b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx @@ -92,20 +92,6 @@ public: cv.DYe /= sqrt_detg; cv.DEnt /= sqrt_detg; - if ((!isfinite(cv.dens)) || (cv.dens <= CCTK_REAL(0))) { - rep.set_nans_in_cons(cv.dens, CCTK_REAL(0), CCTK_REAL(0), CCTK_REAL(0), - cv.DYe); - set_to_nan(pv, cv); - return; - } - - if (cv.dens <= atmo.rho_cut) { - rep.set_atmo_set(); - pv.Bvec = cv.dBvec; - atmo.set(pv, cv, glo); - return; - } - const CCTK_REAL Ssq = calc_contraction(calc_contraction(gup, cv.mom), cv.mom); const CCTK_REAL Bsq = @@ -178,7 +164,7 @@ public: if (std::isfinite(a) && std::isfinite(b) && b >= a) { const CCTK_REAL width = b - a; const CCTK_REAL scale = - fmax(CCTK_REAL(1.0), fmax(std::abs(a), std::abs(b))); + fmax(CCTK_REAL(0.0), fmax(std::abs(a), std::abs(b))); const CCTK_REAL width_tol = soft_root_width_factor * tolerance * scale; const CCTK_REAL mu_ref =