diff --git a/Con2PrimFactory/interface.ccl b/Con2PrimFactory/interface.ccl index f6298a19..d5f4a793 100644 --- a/Con2PrimFactory/interface.ccl +++ b/Con2PrimFactory/interface.ccl @@ -6,6 +6,5 @@ INCLUDES HEADER: prims.hxx IN prims.hxx 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 -USES INCLUDE HEADER: eos.hxx -USES INCLUDE HEADER: eos_idealgas.hxx +USES INCLUDE HEADER: setup_eos.hxx USES INCLUDE HEADER: roots.hxx diff --git a/Con2PrimFactory/schedule.ccl b/Con2PrimFactory/schedule.ccl index 169aa277..ea1ffd15 100644 --- a/Con2PrimFactory/schedule.ccl +++ b/Con2PrimFactory/schedule.ccl @@ -1,7 +1,7 @@ #Schedule definitions for thorn Con2PrimFactory if (unit_test) { - SCHEDULE Con2PrimFactory_Test AT wragh { + SCHEDULE Con2PrimFactory_Test AT CCTK_WRAGH AFTER EOSX_Setup_EOS { LANG: C OPTIONS : meta } diff --git a/Con2PrimFactory/src/c2p.hxx b/Con2PrimFactory/src/c2p.hxx index b9644985..d8b999b3 100644 --- a/Con2PrimFactory/src/c2p.hxx +++ b/Con2PrimFactory/src/c2p.hxx @@ -13,8 +13,7 @@ c2p is effectively an interface to be used by different c2p implementations. #include #include #include "c2p_utils.hxx" -#include "eos.hxx" -#include "eos_idealgas.hxx" +#include "setup_eos.hxx" #include #include "prims.hxx" diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx index 44ba27b9..9d4e466d 100644 --- a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx +++ b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx @@ -13,7 +13,7 @@ public: /* Constructor */ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DPalenzuela( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len); CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL @@ -34,18 +34,18 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, CCTK_REAL Bsq, - CCTK_REAL BiSi, EOSType &eos_th, prim_vars &pv, + CCTK_REAL BiSi, EOSType &eos_3p, prim_vars &pv, cons_vars cv, const smat &gup, const smat &glo) const; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, - CCTK_REAL x, EOSType &eos_th, cons_vars &cv) const; + CCTK_REAL x, EOSType &eos_3p, cons_vars &cv) const; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv, + solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv, const smat &glo, c2p_report &rep) const; /* Destructor */ @@ -56,10 +56,10 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DPalenzuela::c2p_1DPalenzuela( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len) { - GammaIdealFluid = eos_th.gamma; + GammaIdealFluid = eos_3p.gamma; maxIterations = maxIter; tolerance = tol; rho_strict = rho_str; @@ -131,7 +131,7 @@ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void c2p_1DPalenzuela::xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, - EOSType &eos_th, prim_vars &pv, + EOSType &eos_3p, prim_vars &pv, cons_vars cv, const smat &gup, const smat &glo) const { const CCTK_REAL qPalenzuela = cv.tau / cv.dens; @@ -184,7 +184,7 @@ c2p_1DPalenzuela::xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, pv.Ye = cv.dYe / cv.dens; - pv.press = eos_th.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.Bvec = cv.dBvec; @@ -196,7 +196,7 @@ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL c2p_1DPalenzuela::funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, CCTK_REAL x, - EOSType &eos_th, cons_vars &cv) const { + EOSType &eos_3p, cons_vars &cv) const { // computes f(x) from x and q,r,s,t const CCTK_REAL qPalenzuela = cv.tau / cv.dens; const CCTK_REAL rPalenzuela = Ssq / pow(cv.dens, 2); @@ -223,14 +223,14 @@ c2p_1DPalenzuela::funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, // (iv) CCTK_REAL P_loc = - eos_th.press_from_valid_rho_eps_ye(rho_loc, eps_loc, Ye_loc); + eos_3p.press_from_valid_rho_eps_ye(rho_loc, eps_loc, Ye_loc); return (x - (1.0 + eps_loc + P_loc / rho_loc) * W_loc); } template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void -c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, +c2p_1DPalenzuela::solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv, const smat &glo, c2p_report &rep) const { @@ -311,7 +311,7 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, CCTK_REAL a = xPalenzuela_lowerBound; CCTK_REAL b = xPalenzuela_upperBound; auto fn = [&](auto x) { - return funcRoot_1DPalenzuela(Ssq, Bsq, BiSi, x, eos_th, cv); + return funcRoot_1DPalenzuela(Ssq, Bsq, BiSi, x, eos_3p, cv); }; auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters); @@ -334,7 +334,7 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, status = ROOTSTAT::NOT_CONVERGED; } - xPalenzuelaToPrim(xPalenzuela_Sol, Ssq, Bsq, BiSi, eos_th, pv, cv, gup, glo); + xPalenzuelaToPrim(xPalenzuela_Sol, Ssq, Bsq, BiSi, eos_3p, pv, cv, gup, glo); // set to atmo if computed rho is below floor density if (pv.rho < atmo.rho_cut) { @@ -344,7 +344,7 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, } // check the validity of the computed eps - auto rgeps = eos_th.range_eps_from_valid_rho_ye(pv.rho, pv.Ye); + auto rgeps = eos_3p.range_eps_from_valid_rho_ye(pv.rho, pv.Ye); if (pv.eps > rgeps.max) { //printf("(pv.eps > rgeps.max) is true, adjusting cons.. \n"); rep.adjust_cons = true; @@ -383,8 +383,8 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, } pv.vel *= v_lim / sol_v; pv.w_lor = w_lim; - pv.eps = std::min(std::max(eos_th.rgeps.min, pv.eps), eos_th.rgeps.max); - pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.eps = std::min(std::max(eos_3p.rgeps.min, pv.eps), eos_3p.rgeps.max); + pv.press = eos_3p.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); rep.adjust_cons = true; } diff --git a/Con2PrimFactory/src/c2p_2DNoble.hxx b/Con2PrimFactory/src/c2p_2DNoble.hxx index 87d0c538..39ab7301 100644 --- a/Con2PrimFactory/src/c2p_2DNoble.hxx +++ b/Con2PrimFactory/src/c2p_2DNoble.hxx @@ -15,7 +15,7 @@ public: /* Constructor */ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_2DNoble( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len); CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL get_Ssq_Exact(vec &mom, @@ -45,11 +45,11 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, CCTK_REAL BiSi, - EOSType &eos_th, prim_vars &pv, cons_vars cv, + EOSType &eos_3p, prim_vars &pv, cons_vars cv, const smat &gup) const; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv, + solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv, const smat &glo, c2p_report &rep) const; /* Destructor */ @@ -60,10 +60,10 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_2DNoble::c2p_2DNoble( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len) { - GammaIdealFluid = eos_th.gamma; + GammaIdealFluid = eos_3p.gamma; maxIterations = maxIter; tolerance = tol; rho_strict = rho_str; @@ -151,7 +151,7 @@ c2p_2DNoble::get_dPdVsq_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, - CCTK_REAL BiSi, EOSType &eos_th, prim_vars &pv, + CCTK_REAL BiSi, EOSType &eos_3p, prim_vars &pv, cons_vars cv, const smat &gup) const { CCTK_REAL W_Sol = 1.0 / sqrt(1.0 - vsq_Sol); @@ -178,7 +178,7 @@ c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, pv.Ye = cv.dYe / cv.dens; - pv.press = eos_th.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.Bvec = cv.dBvec; @@ -188,7 +188,7 @@ c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void -c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, +c2p_2DNoble::solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv, const smat &glo, c2p_report &rep) const { @@ -260,12 +260,12 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, pv_seeds.rho = cv.dens / pv_seeds.w_lor; CCTK_REAL eps_min = - eos_th.range_eps_from_valid_rho_ye(pv_seeds.rho, pv_seeds.Ye).min; + eos_3p.range_eps_from_valid_rho_ye(pv_seeds.rho, pv_seeds.Ye).min; CCTK_REAL eps_last = max({pv_seeds.eps, eps_min}); /* get pressure seed from updated pv_seeds.rho */ pv_seeds.press = - eos_th.press_from_valid_rho_eps_ye(pv_seeds.rho, eps_last, pv_seeds.Ye); + eos_3p.press_from_valid_rho_eps_ye(pv_seeds.rho, eps_last, pv_seeds.Ye); /* get Z seed */ CCTK_REAL Z_Seed = @@ -387,7 +387,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, CCTK_REAL vsq_Sol = x[1]; /* Write prims if C2P succeeded */ - WZ2Prim(Z_Sol, vsq_Sol, Bsq, BiSi, eos_th, pv, cv, gup); + WZ2Prim(Z_Sol, vsq_Sol, Bsq, BiSi, eos_3p, pv, cv, gup); // set to atmo if computed rho is below floor density if (pv.rho < atmo.rho_cut) { @@ -397,7 +397,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, } // check the validity of the computed eps - auto rgeps = eos_th.range_eps_from_valid_rho_ye(pv.rho, pv.Ye); + auto rgeps = eos_3p.range_eps_from_valid_rho_ye(pv.rho, pv.Ye); if (pv.eps > rgeps.max) { //printf("(pv.eps > rgeps.max) is true, adjusting cons.. \n"); rep.adjust_cons = true; @@ -434,8 +434,8 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, } pv.vel *= v_lim / sol_v; pv.w_lor = w_lim; - pv.eps = std::min(std::max(eos_th.rgeps.min, pv.eps), eos_th.rgeps.max); - pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.eps = std::min(std::max(eos_3p.rgeps.min, pv.eps), eos_3p.rgeps.max); + pv.press = eos_3p.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); rep.adjust_cons = true; } diff --git a/Con2PrimFactory/src/test.cxx b/Con2PrimFactory/src/test.cxx index f91fd845..1334212f 100644 --- a/Con2PrimFactory/src/test.cxx +++ b/Con2PrimFactory/src/test.cxx @@ -7,8 +7,8 @@ #include "c2p.hxx" #include "c2p_1DPalenzuela.hxx" #include "c2p_2DNoble.hxx" -#include -#include + +#include "setup_eos.hxx" namespace Con2PrimFactory { @@ -18,10 +18,6 @@ using namespace EOSX; extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; - // Initializing the IG EOS - eos::range rgeps(0, 100), rgrho(1e-13, 20), rgye(0.5, 0.5); - const eos_idealgas eos_th(2.0, 938.985, rgeps, rgrho, rgye); - // Setting up atmosphere atmosphere atmo(1e-10, 1e-8, 0.5, 1e-8, 0.001); @@ -30,8 +26,8 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { 1.0, 0.0, 1.0}; // xx, xy, xz, yy, yz, zz // Con2Prim objects - c2p_2DNoble c2p_Noble(eos_th, atmo, 100, 1e-8, 1e8, 1, 1, true); - c2p_1DPalenzuela c2p_Pal(eos_th, atmo, 100, 1e-8, 1e8, 1, 1, true); + c2p_2DNoble c2p_Noble(*eos_3p_ig, atmo, 100, 1e-8, 1e8, 1, 1, true); + c2p_1DPalenzuela c2p_Pal(*eos_3p_ig, atmo, 100, 1e-8, 1e8, 1, 1, true); // Construct error report object: c2p_report rep_Noble; @@ -49,7 +45,7 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { // Testing C2P Noble CCTK_VINFO("Testing C2P Noble..."); - c2p_Noble.solve(eos_th, pv, pv_seeds, cv, g, rep_Noble); + c2p_Noble.solve(*eos_3p_ig, pv, pv_seeds, cv, g, rep_Noble); printf("pv_seeds, pv: \n" "rho: %f, %f \n" @@ -90,7 +86,7 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { // Testing C2P Palenzuela CCTK_VINFO("Testing C2P Palenzuela..."); - c2p_Pal.solve(eos_th, pv, pv_seeds, cv, g, rep_Pal); + c2p_Pal.solve(*eos_3p_ig, pv, pv_seeds, cv, g, rep_Pal); printf("pv_seeds, pv: \n" "rho: %f, %f \n"