factor out some terms to have less sqrt calls

AMReX-Astro · Jun 20, 2024 · f361dc3 · f361dc3
1 parent 30469dd
commit f361dc3
Showing 1 changed file with 9 additions and 6 deletions.
diff --git a/screening/screen.H b/screening/screen.H
@@ -802,15 +802,18 @@ void chabrier1998_helmholtz_F(const amrex::Real gamma, const amrex::Real dgamma_
     // Helmholtz free energy, See Chabrier & Potekhin 1998 Eq. 28
 
     // Fitted parameters, see Chabrier & Potekhin 1998 Sec.IV
+
     constexpr amrex::Real A_1 = -0.9052_rt;
     constexpr amrex::Real A_2 = 0.6322_rt;
-    constexpr amrex::Real A_3 = -0.5_rt * gcem::sqrt(3.0_rt) - A_1 / gcem::sqrt(A_2);
+    constexpr amrex::Real sqrt_A2 = gcem::sqrt(A_2);
+    constexpr amrex::Real A_3 = -0.5_rt * gcem::sqrt(3.0_rt) - A_1 / sqrt_A2;
 
     // Precompute some expressions that are reused in the derivative
+
     const amrex::Real sqrt_gamma = std::sqrt(gamma);
     const amrex::Real sqrt_1_gamma_A2 = std::sqrt(1.0_rt + gamma/A_2);
-    const amrex::Real sqrt_gamma_A2_gamma = std::sqrt(gamma * (A_2 + gamma));
-    const amrex::Real sqrt_gamma_A2 = std::sqrt(gamma/A_2);
+    const amrex::Real sqrt_gamma_A2_gamma = sqrt_gamma * sqrt_A2 * sqrt_1_gamma_A2;
+    const amrex::Real sqrt_gamma_A2 = sqrt_gamma / sqrt_A2;
 
     f = A_1 * (sqrt_gamma_A2_gamma -
                A_2 * std::log(sqrt_gamma_A2 + sqrt_1_gamma_A2)) +
@@ -873,8 +876,8 @@ void chabrier1998 (const plasma_state_t& state,
     // Now we add quantum correction terms discussed in Alastuey 1978.
     // Notice in Alastuey 1978, they have a different classical term,
     // which is implemented in the strong screening limit of our screen5 routine.
-    amrex::Real quantum_corr_1 = 0.0_rt;
 
+    amrex::Real quantum_corr_1 = 0.0_rt;
     amrex::Real quantum_corr_2 = 0.0_rt;
 
     [[maybe_unused]] amrex::Real quantum_corr_1_dT = 0.0_rt;
@@ -883,8 +886,8 @@ void chabrier1998 (const plasma_state_t& state,
     if (screening_rp::enable_chabrier1998_quantum_corr) {
         // See Wallace1982, Eq. A13
 
-        amrex::Real Gamma_eff = std::cbrt(2.0_rt) * scn_fac.z1 * scn_fac.z2 *
-                         scn_fac.zs13inv * Gamma_e;
+        amrex::Real Gamma_eff = 1.25992104989_rt * scn_fac.z1 * scn_fac.z2 *
+                                scn_fac.zs13inv * Gamma_e;
         [[maybe_unused]] amrex::Real Gamma_eff_dT;
 
         if constexpr (do_T_derivatives) {