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; }