some cleaning of screen5

screening/screen.H

@@ -185,7 +185,7 @@ number_t debye_huckel (const plasma_state_t<number_t>& state,
 template <typename number_t>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 number_t actual_screen5 (const plasma_state_t<number_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<number_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,77 +253,64 @@ number_t actual_screen5 (const plasma_state_t<number_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;
 
     // intermediate and strong sceening regime
 
     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
         h12 = cc - a3 * (ss + vv);
 
         // 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

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

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