diff --git a/screening/screen.H b/screening/screen.H index 97a6d9f20..15691c28a 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -185,7 +185,7 @@ number_t debye_huckel (const plasma_state_t& state, template AMREX_GPU_HOST_DEVICE AMREX_INLINE number_t actual_screen5 (const plasma_state_t& state, - const scrn::screen_factors_t& scn_fac) + const scrn::screen_factors_t& scn_fac) { // this subroutine calculates screening factors and their derivatives // for nuclear reaction rates in the weak, intermediate and strong regimes. @@ -200,17 +200,18 @@ number_t actual_screen5 (const plasma_state_t& state, // fact = 2^(1/3) - const amrex::Real fact = 1.25992104989487e0_rt; - const amrex::Real gamefx = 0.3e0_rt; // lower gamma limit for intermediate screening - const amrex::Real gamefs = 0.8e0_rt; // upper gamma limit for intermediate screening - const amrex::Real h12_max = 300.e0_rt; + constexpr amrex::Real fact = 1.25992104989487e0_rt; + constexpr amrex::Real gamefx = 0.3e0_rt; // lower gamma limit for intermediate screening + constexpr amrex::Real gamefs = 0.8e0_rt; // upper gamma limit for intermediate screening + constexpr amrex::Real inv_dg = 1.0_rt / (gamefs - gamefx); + constexpr amrex::Real h12_max = 300.e0_rt; // Get the ion data based on the input index - amrex::Real z1 = scn_fac.z1; - amrex::Real z2 = scn_fac.z2; + const amrex::Real z1 = scn_fac.z1; + const amrex::Real z2 = scn_fac.z2; // calculate individual screening factors - amrex::Real bb = z1 * z2; + const amrex::Real bb = z1 * z2; number_t gamp = state.aa; // In Eq.4 in Itoh:1979, this term is 2*Z_1*Z_2/(Z_1^(1/3) + Z_2^(1/3)) @@ -252,7 +253,7 @@ number_t actual_screen5 (const plasma_state_t& state, // Full version of Eq. 19 in Graboske:1973 by considering weak regime // and Wallace:1982 Eq. A14. Here the degeneracy factor is assumed to be 1. - number_t h12w = bb * state.qlam0z; + const number_t h12w = bb * state.qlam0z; number_t h12 = h12w; @@ -260,34 +261,33 @@ number_t actual_screen5 (const plasma_state_t& state, if (gamef > gamefx) { - // gamma_ij^(1/4) - - number_t gamp14 = admath::pow(gamp, 0.25_rt); + // gamma_ij^(1/4) + const number_t gamp14 = admath::sqrt(admath::sqrt(gamp)); // Here we follow Eq. A9 in Wallace:1982 // See Eq. 25 Alastuey:1978, Eq. 16 and 17 in Jancovici:1977 for reference - number_t cc = 0.896434e0_rt * gamp * scn_fac.zhat + const number_t cc = 0.896434e0_rt * gamp * scn_fac.zhat - 3.44740e0_rt * gamp14 * scn_fac.zhat2 - 0.5551e0_rt * (admath::log(gamp) + scn_fac.lzav) - 2.996e0_rt; // (3gamma_ij/tau_ij)^3 - number_t a3 = alph12 * alph12 * alph12; + const number_t a3 = alph12 * alph12 * alph12; // Part of Eq. 28 in Alastuey:1978 - number_t rr = (5.0_rt/32.0_rt) - alph12*(0.014e0_rt + 0.0128e0_rt*alph12); + const number_t rr = (5.0_rt/32.0_rt) - alph12*(0.014e0_rt + 0.0128e0_rt*alph12); // Part of Eq. 28 in Alastuey:1978 - number_t ss = tau12*rr; + const number_t ss = tau12*rr; // Part of Eq. 31 in Alastuey:1978 - number_t tt = -0.0098e0_rt + 0.0048e0_rt*alph12; + const number_t tt = -0.0098e0_rt + 0.0048e0_rt*alph12; // Part of Eq. 31 in Alastuey:1978 - number_t uu = 0.0055e0_rt + alph12*tt; + const number_t uu = 0.0055e0_rt + alph12*tt; // Part of Eq. 31 in Alastuey:1978 - number_t vv = gamef * alph12 * uu; + const number_t vv = gamef * alph12 * uu; // Exponent of Eq. 32 in Alastuey:1978, which uses Eq.28 and Eq.31 // Strong screening factor @@ -295,34 +295,22 @@ number_t actual_screen5 (const plasma_state_t& state, // See conclusion and Eq. 34 in Alastuey:1978 // This is an extra factor to account for quantum effects - rr = 1.0_rt - 0.0562e0_rt*a3; - - number_t xlgfac; + // In extreme case, limit to 0.77, see conclusion in Alastuey:1978 - // In extreme case, rr is 0.77, see conclusion in Alastuey:1978 - if (rr >= 0.77e0_rt) { - xlgfac = rr; - } else { - xlgfac = 0.77e0_rt; - } + const number_t xlgfac = admath::max(1.0_rt - 0.0562e0_rt*a3, 0.77_rt); // Include the extra factor that accounts for quantum effects - h12 = admath::log(xlgfac) + h12; + h12 += admath::log(xlgfac); // If gamma_ij < upper limit of intermediate regime // then it is in the intermediate regime, else strong screening. if (gamef <= gamefs) { - amrex::Real dgamma = 1.0e0_rt/(gamefs - gamefx); - - rr = dgamma*(gamefs - gamef); - - ss = dgamma*(gamef - gamefx); - - vv = h12; + const number_t weight1 = inv_dg * (gamefs - gamef); + const number_t weight2 = inv_dg * (gamef - gamefx); // Then the screening factor is a combination // of the strong and weak screening factor. - h12 = h12w*rr + vv*ss; + h12 = h12w * weight1 + h12 * weight2; } // end of intermediate and strong screening diff --git a/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out b/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out index b6fe8d05f..7dc2ea262 100644 --- a/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out @@ -1,5 +1,5 @@ -Initializing AMReX (24.10-178-g5a2390e00837)... -AMReX (24.10-178-g5a2390e00837) initialized +Initializing AMReX (25.12-25-g10b67a795031)... +AMReX (25.12-25-g10b67a795031) initialized starting the single zone burn... Maximum Time (s): 0.01585 State Density (g/cm^3): 1000000000 @@ -13,21 +13,21 @@ RHS at t = 0 ash 0.01230280576 ------------------------------------ successful? 1 - - Hnuc = 5.277066077e+17 - - added e = 8.364149732e+15 - - final T = 1433686831 + - Hnuc = 5.277014299e+17 + - added e = 8.364067664e+15 + - final T = 1433682844 ------------------------------------ e initial = 1.253426044e+18 -e final = 1.261790194e+18 +e final = 1.261790112e+18 ------------------------------------ new mass fractions: -C12 0.9657916864 +C12 0.965792022 O16 1e-30 -ash 0.03420831361 +ash 0.03420797796 ------------------------------------ species creation rates: -omegadot(C12): -2.158253225 -omegadot(O16): 1.989224949e-43 -omegadot(ash): 2.158253225 +omegadot(C12): -2.158232048 +omegadot(O16): 7.735874803e-44 +omegadot(ash): 2.158232048 number of steps taken: 381 -AMReX (24.10-178-g5a2390e00837) finalized +AMReX (25.12-25-g10b67a795031) finalized diff --git a/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out b/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out index 530f28259..75357037c 100644 --- a/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out @@ -1,5 +1,5 @@ -Initializing AMReX (25.09-20-g6e2d7bbde65e)... -AMReX (25.09-20-g6e2d7bbde65e) initialized +Initializing AMReX (25.12-25-g10b67a795031)... +AMReX (25.12-25-g10b67a795031) initialized starting the single zone burn... reading in network electron-capture / beta-decay tables... Maximum Time (s): 0.01 @@ -30,37 +30,37 @@ RHS at t = 0 S32 1323519.174 ------------------------------------ successful? 1 - - Hnuc = 5.491366142e+19 - - added e = 5.491366142e+17 - - final T = 7126065275 + - Hnuc = 4.572458836e+19 + - added e = 4.572458836e+17 + - final T = 6740499652 ------------------------------------ e initial = 5.995956082e+17 -e final = 1.148732222e+18 +e final = 1.056841492e+18 ------------------------------------ new mass fractions: -H1 0.00900254878 -He4 6.899608988e-13 -O16 5.02310191e-07 -O20 0.009835982937 -F20 0.009834632311 -Ne20 4.012407296e-14 -Mg24 7.687349513e-06 -Al27 2.397057447e-17 -Si28 0.2644621097 -P31 1.618654609e-13 -S32 0.7068565366 +H1 0.008875172349 +He4 9.999999998e-31 +O16 6.258606287e-07 +O20 0.009862665409 +F20 0.009861311457 +Ne20 8.532947059e-14 +Mg24 1.247632983e-06 +Al27 1.79296072e-18 +Si28 0.3055707459 +P31 9.999999998e-31 +S32 0.6658282314 ------------------------------------ species creation rates: -omegadot(H1): -0.09974512199 +omegadot(H1): -0.1124827651 omegadot(He4): -3 -omegadot(O16): -49.99994977 -omegadot(O20): -0.01640170629 -omegadot(F20): -0.01653676888 +omegadot(O16): -49.99993741 +omegadot(O20): -0.01373345906 +omegadot(F20): -0.01386885428 omegadot(Ne20): -30 -omegadot(Mg24): -9.999231265 +omegadot(Mg24): -9.999875237 omegadot(Al27): -1 -omegadot(Si28): 25.44621097 +omegadot(Si28): 29.55707459 omegadot(P31): -1 -omegadot(S32): 69.68565366 -number of steps taken: 13271 -AMReX (25.09-20-g6e2d7bbde65e) finalized +omegadot(S32): 65.58282314 +number of steps taken: 16498 +AMReX (25.12-25-g10b67a795031) finalized