From 635c36c6ff74638fd5275b356eb719fa6e359cf1 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Tue, 28 Nov 2023 17:10:08 -0500 Subject: [PATCH 1/3] organize screening, put z1 z2 related pow vars to screen_factors_t (#1397) --- screening/screen.H | 22 +++++++------ screening/screen_data.H | 71 +++++++++++++++++++++++++++++++++++------ 2 files changed, 74 insertions(+), 19 deletions(-) diff --git a/screening/screen.H b/screening/screen.H index f07fd9f4dd..f702a533cf 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -169,6 +169,7 @@ fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Array } } +#if SCREEN_METHOD == 1 template AMREX_GPU_HOST_DEVICE AMREX_INLINE void actual_screen5 (const plasma_state_t& state, @@ -416,6 +417,7 @@ void actual_screen5 (const plasma_state_t& state, } } +#elif SCREEN_METHOD == 2 template AMREX_GPU_HOST_DEVICE AMREX_INLINE void chugunov2007 (const plasma_state_t& state, @@ -460,7 +462,7 @@ void chugunov2007 (const plasma_state_t& state, // m_i -> m_u * abar Real mu12 = scn_fac.a1 * scn_fac.a2 / (scn_fac.a1 + scn_fac.a2); Real z_factor = scn_fac.z1 * scn_fac.z2; - Real n_i = state.n_e / gcem::pow(scn_fac.ztilde, 3); + Real n_i = state.n_e / scn_fac.ztilde3; Real m_i = 2.0_rt * mu12 / C::n_A; constexpr Real T_p_factor = C::hbar/C::k_B*C::q_e*gcem::sqrt(4.0_rt*GCEM_PI); @@ -611,6 +613,7 @@ void chugunov2007 (const plasma_state_t& state, } } +#elif SCREEN_METHOD == 3 template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void chugunov2009_f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT) @@ -698,9 +701,9 @@ void chugunov2009 (const plasma_state_t& state, Real dlog_Gamma_dT = -1.0_rt / state.temp; // Coulomb coupling parameters for ions and compound nucleus, eqs. 7 & 9 - Real Gamma_1 = Gamma_e * std::pow(scn_fac.z1, 5.0_rt/3.0_rt); - Real Gamma_2 = Gamma_e * std::pow(scn_fac.z2, 5.0_rt/3.0_rt); - Real Gamma_comp = Gamma_e * std::pow(zcomp, 5.0_rt/3.0_rt); + Real Gamma_1 = Gamma_e * scn_fac.z1_53; + Real Gamma_2 = Gamma_e * scn_fac.z1_53; + Real Gamma_comp = Gamma_e * scn_fac.zs53; Real Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde; @@ -753,7 +756,7 @@ void chugunov2009 (const plasma_state_t& state, // weak screening correction term, eq. A3 Real corr_C = 3.0_rt*z1z2 * std::sqrt(state.z2bar/state.zbar) / - (std::pow(zcomp, 2.5_rt) - std::pow(scn_fac.z1, 2.5_rt) - std::pow(scn_fac.z2, 2.5_rt)); + (scn_fac.zs52 - scn_fac.z1_52 - scn_fac.z2_52); // corrected enhancement factor, eq. A4 Real Gamma_12_2 = Gamma_12 * Gamma_12; @@ -780,6 +783,7 @@ void chugunov2009 (const plasma_state_t& state, } } +#elif SCREEN_METHOD == 4 template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void chabrier1998_helmholtz_F(const Real gamma, const Real dgamma_dT, Real& f, Real& df_dT) { @@ -830,13 +834,12 @@ void chabrier1998 (const plasma_state_t& state, // Eq. 2 in Chabrier & Potekhin 1998 Real Gamma_e = state.gamma_e_fac / state.temp; - Real zcomp = scn_fac.z1 + scn_fac.z2; // See Calder2007 appendix Eq. A9 - Real Gamma1 = Gamma_e * std::pow(scn_fac.z1, 5.0_rt/3.0_rt); - Real Gamma2 = Gamma_e * std::pow(scn_fac.z2, 5.0_rt/3.0_rt); - Real Gamma12 = Gamma_e * std::pow(zcomp, 5.0_rt/3.0_rt); + Real Gamma1 = Gamma_e * scn_fac.z1_53; + Real Gamma2 = Gamma_e * scn_fac.z2_53; + Real Gamma12 = Gamma_e * scn_fac.zs53; Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{}; @@ -934,6 +937,7 @@ void chabrier1998 (const plasma_state_t& state, } } +#endif template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/screening/screen_data.H b/screening/screen_data.H index 8bf2d255e9..77021fc581 100644 --- a/screening/screen_data.H +++ b/screening/screen_data.H @@ -17,20 +17,50 @@ namespace scrn { amrex::Real a1 = -1; amrex::Real a2 = -1; + // z1_53 = (z1)**(5./3.) + // z2_53 = (z2)**(5./3.) + // zs52 = (z1+z2)**(5./2.) + // z1_52 = (z1)**(5./2.) + // z2_52 = (z2)**(5./2.) // zs13 = (z1+z2)**(1./3.) + // zs53 = (z1+z2)**(5./3.) // zhat = combination of z1 and z2 raised to the 5/3 power // zhat2 = combination of z1 and z2 raised to the 5/12 power // lzav = log of effective charge // aznut = combination of a1,z1,a2,z2 raised to 1/3 power // ztilde = effective ion radius factor for a MCP + // ztilde3 = ztilde**3 +#if SCREEN_METHOD == 1 + amrex::Real zs53 = 0.0; + amrex::Real z1_53 = 0.0; + amrex::Real z2_53 = 0.0; amrex::Real zs13 = 0.0; amrex::Real zs13inv = 0.0; amrex::Real zhat = 0.0; amrex::Real zhat2 = 0.0; amrex::Real lzav = 0.0; amrex::Real aznut = 0.0; +#elif SCREEN_METHOD == 2 amrex::Real ztilde = 0.0; + amrex::Real ztilde3 = 0.0; +#elif SCREEN_METHOD == 3 + amrex::Real zs52 = 0.0; + amrex::Real z1_52 = 0.0; + amrex::Real z2_52 = 0.0; + amrex::Real zs53 = 0.0; + amrex::Real z1_53 = 0.0; + amrex::Real z2_53 = 0.0; + amrex::Real aznut = 0.0; + amrex::Real ztilde = 0.0; +#elif SCREEN_METHOD == 4 + amrex::Real zs53 = 0.0; + amrex::Real z1_53 = 0.0; + amrex::Real z2_53 = 0.0; + amrex::Real zs13 = 0.0; + amrex::Real zs13inv = 0.0; + amrex::Real aznut = 0.0; +#endif [[nodiscard]] bool validate_nuclei(const amrex::Real z1_pass, const amrex::Real a1_pass, @@ -57,25 +87,46 @@ namespace scrn { scn_fac.z2 = z2; scn_fac.a2 = a2; +#if SCREEN_METHOD == 1 + scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt); + scn_fac.z1_53 = gcem::pow(z1, 5.0_rt / 3.0_rt); + scn_fac.z2_53 = gcem::pow(z2, 5.0_rt / 3.0_rt); scn_fac.zs13 = gcem::pow(z1 + z2, 1.0_rt / 3.0_rt); - scn_fac.zs13inv = 1.0_rt / scn_fac.zs13; - - scn_fac.zhat = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt) - - gcem::pow(z1, 5.0_rt/3.0_rt) - - gcem::pow(z2, 5.0_rt/3.0_rt); - + scn_fac.zhat = scn_fac.zs53 - scn_fac.z1_53 - scn_fac.z2_53; scn_fac.zhat2 = gcem::pow(z1 + z2, 5.0_rt / 12.0_rt) - - gcem::pow(z1, 5.0_rt/12.0_rt) - - gcem::pow(z2, 5.0_rt/12.0_rt); - + gcem::pow(z1, 5.0_rt / 12.0_rt) - + gcem::pow(z2, 5.0_rt / 12.0_rt); scn_fac.lzav = (5.0_rt / 3.0_rt) * gcem::log(z1 * z2 / (z1 + z2)); - scn_fac.aznut = gcem::pow(z1 * z1 * z2 * z2 * a1 * a2 / (a1 + a2), 1.0_rt / 3.0_rt); +#elif SCREEN_METHOD == 2 scn_fac.ztilde = 0.5_rt * (gcem::pow(z1, 1.0_rt / 3.0_rt) + gcem::pow(z2, 1.0_rt / 3.0_rt)); + scn_fac.ztilde3 = gcem::pow(scn_fac.ztilde, 3); + +#elif SCREEN_METHOD == 3 + scn_fac.zs52 = gcem::pow(z1 + z2, 5.0_rt / 2.0_rt); + scn_fac.z1_52 = gcem::pow(z1, 5.0_rt / 2.0_rt); + scn_fac.z2_52 = gcem::pow(z2, 5.0_rt / 2.0_rt); + scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt); + scn_fac.z1_53 = gcem::pow(z1, 5.0_rt / 3.0_rt); + scn_fac.z2_53 = gcem::pow(z2, 5.0_rt / 3.0_rt); + scn_fac.aznut = gcem::pow(z1 * z1 * z2 * z2 * a1 * a2 / (a1 + a2), + 1.0_rt / 3.0_rt); + scn_fac.ztilde = 0.5_rt * (gcem::pow(z1, 1.0_rt / 3.0_rt) + + gcem::pow(z2, 1.0_rt / 3.0_rt)); + +#elif SCREEN_METHOD == 4 + scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt); + scn_fac.z1_53 = gcem::pow(z1, 5.0_rt/3.0_rt); + scn_fac.z2_53 = gcem::pow(z2, 5.0_rt/3.0_rt); + scn_fac.zs13 = gcem::pow(z1 + z2, 1.0_rt / 3.0_rt); + scn_fac.zs13inv = 1.0_rt / scn_fac.zs13; + scn_fac.aznut = gcem::pow(z1 * z1 * z2 * z2 * a1 * a2 / (a1 + a2), + 1.0_rt / 3.0_rt); +#endif return scn_fac; } From f03189021b9d0411f9c9479050ac1e5794cd6282 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 28 Nov 2023 19:13:09 -0500 Subject: [PATCH 2/3] fix a factor in chug2009 screening (#1402) --- screening/screen.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/screening/screen.H b/screening/screen.H index f702a533cf..81b1cf7701 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -702,7 +702,7 @@ void chugunov2009 (const plasma_state_t& state, // Coulomb coupling parameters for ions and compound nucleus, eqs. 7 & 9 Real Gamma_1 = Gamma_e * scn_fac.z1_53; - Real Gamma_2 = Gamma_e * scn_fac.z1_53; + Real Gamma_2 = Gamma_e * scn_fac.z2_53; Real Gamma_comp = Gamma_e * scn_fac.zs53; Real Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde; From f104eff31026ebaef3a4c2056c3ca1b49e3f2178 Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Tue, 28 Nov 2023 16:33:34 -0800 Subject: [PATCH 3/3] Add constants for the allowed values of SCREEN_METHOD (#1401) An enum would be much nicer here, but the preprocessor doesn't understand those. --- Make.Microphysics_extern | 12 ++++++------ nse_solver/nse_solver.H | 4 ++-- screening/screen.H | 35 +++++++++++++++++++++-------------- screening/screen_data.H | 16 ++++++++-------- 4 files changed, 37 insertions(+), 30 deletions(-) diff --git a/Make.Microphysics_extern b/Make.Microphysics_extern index 063e029934..28526b39d2 100644 --- a/Make.Microphysics_extern +++ b/Make.Microphysics_extern @@ -23,15 +23,15 @@ endif SCREEN_METHOD ?= screen5 ifeq ($(SCREEN_METHOD), null) - DEFINES += -DSCREEN_METHOD=0 + DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_null else ifeq ($(SCREEN_METHOD), screen5) - DEFINES += -DSCREEN_METHOD=1 + DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_screen5 else ifeq ($(SCREEN_METHOD), chugunov2007) - DEFINES += -DSCREEN_METHOD=2 + DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_chugunov2007 else ifeq ($(SCREEN_METHOD), chugunov2009) - DEFINES += -DSCREEN_METHOD=3 + DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_chugunov2009 else ifeq ($(SCREEN_METHOD), chabrier1998) - DEFINES += -DSCREEN_METHOD=4 + DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_chabrier1998 else $(error Invalid value for SCREEN_METHOD) endif @@ -155,4 +155,4 @@ endif clean:: @if [ -L helm_table.dat ]; then rm -f helm_table.dat; fi @if [ -L reaclib_rate_metadata.dat ]; then rm -f reaclib_rate_metadata.dat; fi - $(foreach t, $(wildcard *_betadecay.dat *_electroncapture.dat nse*.tbl), $(shell if [ -L $t ]; then rm -f $t; fi)) \ No newline at end of file + $(foreach t, $(wildcard *_betadecay.dat *_electroncapture.dat nse*.tbl), $(shell if [ -L $t ]; then rm -f $t; fi)) diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H index b357f2c16b..3b6578f40f 100644 --- a/nse_solver/nse_solver.H +++ b/nse_solver/nse_solver.H @@ -92,7 +92,7 @@ void apply_nse_exponent(T& nse_state) { // if we use chabrier1998 screening // Get the required terms to calculate coulomb correction term, u_c -#if SCREEN_METHOD == 4 +#if SCREEN_METHOD == SCREEN_METHOD_chabrier1998 // Find n_e for original state; const amrex::Real n_e = nse_state.rho * nse_state.y_e / C::m_u; @@ -113,7 +113,7 @@ void apply_nse_exponent(T& nse_state) { // term for calculating u_c // if use chabrier1998 screening, calculate the coulomb correction term -#if SCREEN_METHOD == 4 +#if SCREEN_METHOD == SCREEN_METHOD_chabrier1998 gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e; // chemical potential for coulomb correction diff --git a/screening/screen.H b/screening/screen.H index 81b1cf7701..ffb0b750e9 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -1,6 +1,13 @@ #ifndef SCREEN_H #define SCREEN_H +// these need to be defined before screen_data.H is included +#define SCREEN_METHOD_null 0 +#define SCREEN_METHOD_screen5 1 +#define SCREEN_METHOD_chugunov2007 2 +#define SCREEN_METHOD_chugunov2009 3 +#define SCREEN_METHOD_chabrier1998 4 + #include #include #include @@ -16,15 +23,15 @@ using namespace amrex; using namespace integrator_rp; -#if SCREEN_METHOD == 0 +#if SCREEN_METHOD == SCREEN_METHOD_null const std::string screen_name = "null"; -#elif SCREEN_METHOD == 1 +#elif SCREEN_METHOD == SCREEN_METHOD_screen5 const std::string screen_name = "screen5"; -#elif SCREEN_METHOD == 2 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007 const std::string screen_name = "chugunov2007"; -#elif SCREEN_METHOD == 3 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009 const std::string screen_name = "chugunov2009"; -#elif SCREEN_METHOD == 4 +#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998 const std::string screen_name = "chabrier1998"; #endif @@ -169,7 +176,7 @@ fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Array } } -#if SCREEN_METHOD == 1 +#if SCREEN_METHOD == SCREEN_METHOD_screen5 template AMREX_GPU_HOST_DEVICE AMREX_INLINE void actual_screen5 (const plasma_state_t& state, @@ -417,7 +424,7 @@ void actual_screen5 (const plasma_state_t& state, } } -#elif SCREEN_METHOD == 2 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007 template AMREX_GPU_HOST_DEVICE AMREX_INLINE void chugunov2007 (const plasma_state_t& state, @@ -613,7 +620,7 @@ void chugunov2007 (const plasma_state_t& state, } } -#elif SCREEN_METHOD == 3 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009 template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void chugunov2009_f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT) @@ -783,7 +790,7 @@ void chugunov2009 (const plasma_state_t& state, } } -#elif SCREEN_METHOD == 4 +#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998 template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void chabrier1998_helmholtz_F(const Real gamma, const Real dgamma_dT, Real& f, Real& df_dT) { @@ -945,18 +952,18 @@ void actual_screen(const plasma_state_t& state, const scrn::screen_factors_t& scn_fac, Real& scor, Real& scordt) { -#if SCREEN_METHOD == 0 +#if SCREEN_METHOD == SCREEN_METHOD_null // null screening amrex::ignore_unused(state, scn_fac); scor = 1.0_rt; scordt = 0.0_rt; -#elif SCREEN_METHOD == 1 +#elif SCREEN_METHOD == SCREEN_METHOD_screen5 actual_screen5(state, scn_fac, scor, scordt); -#elif SCREEN_METHOD == 2 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007 chugunov2007(state, scn_fac, scor, scordt); -#elif SCREEN_METHOD == 3 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009 chugunov2009(state, scn_fac, scor, scordt); -#elif SCREEN_METHOD == 4 +#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998 chabrier1998(state, scn_fac, scor, scordt); #endif } diff --git a/screening/screen_data.H b/screening/screen_data.H index 77021fc581..f217fec44b 100644 --- a/screening/screen_data.H +++ b/screening/screen_data.H @@ -31,7 +31,7 @@ namespace scrn { // ztilde = effective ion radius factor for a MCP // ztilde3 = ztilde**3 -#if SCREEN_METHOD == 1 +#if SCREEN_METHOD == SCREEN_METHOD_screen5 amrex::Real zs53 = 0.0; amrex::Real z1_53 = 0.0; amrex::Real z2_53 = 0.0; @@ -41,10 +41,10 @@ namespace scrn { amrex::Real zhat2 = 0.0; amrex::Real lzav = 0.0; amrex::Real aznut = 0.0; -#elif SCREEN_METHOD == 2 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007 amrex::Real ztilde = 0.0; amrex::Real ztilde3 = 0.0; -#elif SCREEN_METHOD == 3 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009 amrex::Real zs52 = 0.0; amrex::Real z1_52 = 0.0; amrex::Real z2_52 = 0.0; @@ -53,7 +53,7 @@ namespace scrn { amrex::Real z2_53 = 0.0; amrex::Real aznut = 0.0; amrex::Real ztilde = 0.0; -#elif SCREEN_METHOD == 4 +#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998 amrex::Real zs53 = 0.0; amrex::Real z1_53 = 0.0; amrex::Real z2_53 = 0.0; @@ -87,7 +87,7 @@ namespace scrn { scn_fac.z2 = z2; scn_fac.a2 = a2; -#if SCREEN_METHOD == 1 +#if SCREEN_METHOD == SCREEN_METHOD_screen5 scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt); scn_fac.z1_53 = gcem::pow(z1, 5.0_rt / 3.0_rt); scn_fac.z2_53 = gcem::pow(z2, 5.0_rt / 3.0_rt); @@ -101,12 +101,12 @@ namespace scrn { scn_fac.aznut = gcem::pow(z1 * z1 * z2 * z2 * a1 * a2 / (a1 + a2), 1.0_rt / 3.0_rt); -#elif SCREEN_METHOD == 2 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007 scn_fac.ztilde = 0.5_rt * (gcem::pow(z1, 1.0_rt / 3.0_rt) + gcem::pow(z2, 1.0_rt / 3.0_rt)); scn_fac.ztilde3 = gcem::pow(scn_fac.ztilde, 3); -#elif SCREEN_METHOD == 3 +#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009 scn_fac.zs52 = gcem::pow(z1 + z2, 5.0_rt / 2.0_rt); scn_fac.z1_52 = gcem::pow(z1, 5.0_rt / 2.0_rt); scn_fac.z2_52 = gcem::pow(z2, 5.0_rt / 2.0_rt); @@ -118,7 +118,7 @@ namespace scrn { scn_fac.ztilde = 0.5_rt * (gcem::pow(z1, 1.0_rt / 3.0_rt) + gcem::pow(z2, 1.0_rt / 3.0_rt)); -#elif SCREEN_METHOD == 4 +#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998 scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt); scn_fac.z1_53 = gcem::pow(z1, 5.0_rt/3.0_rt); scn_fac.z2_53 = gcem::pow(z2, 5.0_rt/3.0_rt);