Move Helmholtz free energy from chabrier1998 into a helper function (#…

…1334) This reduces code duplication and lets us reuse it in the NSE solver. Also pull some common subexpressions into variables so we only calculate them once.
AMReX-Astro · Sep 13, 2023 · 5e38c70 · 5e38c70
1 parent 94b016d
commit 5e38c70
Show file tree

Hide file tree

Showing 2 changed files with 48 additions and 60 deletions.
diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H
@@ -16,6 +16,7 @@
 #include <eos_composition.H>
 #include <microphysics_sort.H>
 #include <hybrj.H>
+#include <screen.H>
 #include <cctype>
 #include <algorithm>
 
@@ -93,15 +94,6 @@ void apply_nse_exponent(T& nse_state) {
 
 #if SCREEN_METHOD != 0
 
-  // three unit-less constants for calculating coulomb correction term
-  // see Calder 2007, doi:10.1086/510709 for more detail
-
-  const amrex::Real A1 = -0.9052_rt;
-  const amrex::Real A2 = 0.6322_rt;
-
-  // A3 = -0.5 * sqrt(3) - A1/sqrt(A2)
-  const amrex::Real A3 = 0.272433346667;
-
   // Find n_e for original state;
   const amrex::Real n_e = nse_state.rho * nse_state.y_e / C::m_u;
   const amrex::Real Gamma_e = C::q_e * C::q_e * std::cbrt(4.0_rt * M_PI * n_e / 3.0_rt)
@@ -125,12 +117,14 @@ void apply_nse_exponent(T& nse_state) {
       gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e;
 
       // chemical potential for coulomb correction
+      // see appendix of Calder 2007, doi:10.1086/510709 for more detail
+
+      // reuse existing implementation from screening routine
+      Real f, df;
+      constexpr int do_T_derivatives = 0;
+      chabrier1998_helmholtz_F<do_T_derivatives>(gamma, 0.0_rt, f, df);
 
-      u_c = C::k_B * T_in / C::Legacy::MeV2erg *
-          (A1 * (std::sqrt(gamma * (A2 + gamma)) -
-                 A2 * std::log(std::sqrt(gamma / A2) +
-                               std::sqrt(1.0_rt + gamma / A2))) +
-           2.0_rt * A3 * (std::sqrt(gamma) - std::atan(std::sqrt(gamma))));
+      u_c = C::k_B * T_in / C::Legacy::MeV2erg * f;
 #endif
 
       // find nse mass frac

diff --git a/screening/screen.H b/screening/screen.H
@@ -613,7 +613,7 @@ void chugunov2007 (const plasma_state_t& state,
 
 template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT)
+void chugunov2009_f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT)
 {
     // Calculates the free energy per ion in a OCP, from Chugunov and DeWitt 2009
     // equation 24.
@@ -742,9 +742,9 @@ void chugunov2009 (const plasma_state_t& state,
     // values similar to those from Chugunov 2007.
     Real term1, term2, term3;
     Real dterm1_dT = 0.0_rt, dterm2_dT = 0.0_rt, dterm3_dT = 0.0_rt;
-    f0<do_T_derivatives>(Gamma_1 / t_12, dlog_dT, term1, dterm1_dT);
-    f0<do_T_derivatives>(Gamma_2 / t_12, dlog_dT, term2, dterm2_dT);
-    f0<do_T_derivatives>(Gamma_comp / t_12, dlog_dT, term3, dterm3_dT);
+    chugunov2009_f0<do_T_derivatives>(Gamma_1 / t_12, dlog_dT, term1, dterm1_dT);
+    chugunov2009_f0<do_T_derivatives>(Gamma_2 / t_12, dlog_dT, term2, dterm2_dT);
+    chugunov2009_f0<do_T_derivatives>(Gamma_comp / t_12, dlog_dT, term3, dterm3_dT);
     Real h_fit = term1 + term2 - term3;
     Real dh_fit_dT;
     if constexpr (do_T_derivatives) {
@@ -780,6 +780,35 @@ void chugunov2009 (const plasma_state_t& state,
     }
 }
 
+template <int do_T_derivatives>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void chabrier1998_helmholtz_F(const Real gamma, const Real dgamma_dT, Real& f, Real& df_dT) {
+  // Helmholtz free energy, See Chabrier & Potekhin 1998 Eq. 28
+
+  // Fitted parameters, see Chabrier & Potekhin 1998 Sec.IV
+  constexpr Real A_1 = -0.9052_rt;
+  constexpr Real A_2 = 0.6322_rt;
+  constexpr Real A_3 = -0.5_rt * gcem::sqrt(3.0_rt) - A_1 / gcem::sqrt(A_2);
+
+  // Precompute some expressions that are reused in the derivative
+  const Real sqrt_gamma = std::sqrt(gamma);
+  const Real sqrt_1_gamma_A2 = std::sqrt(1.0_rt + gamma/A_2);
+  const Real sqrt_gamma_A2_gamma = std::sqrt(gamma * (A_2 + gamma));
+  const Real sqrt_gamma_A2 = std::sqrt(gamma/A_2);
+
+  f = A_1 * (sqrt_gamma_A2_gamma -
+             A_2 * std::log(sqrt_gamma_A2 + sqrt_1_gamma_A2))
+      + 2.0_rt * A_3 * (sqrt_gamma - std::atan(sqrt_gamma));
+
+  if constexpr (do_T_derivatives) {
+    df_dT = A_1 * ((A_2 + 2.0_rt * gamma) / (2.0_rt * sqrt_gamma_A2_gamma)
+                        - 0.5_rt / (sqrt_gamma_A2 + sqrt_1_gamma_A2)
+                        * (1.0_rt / sqrt_gamma_A2 + 1.0_rt / sqrt_1_gamma_A2))
+            + A_3 / sqrt_gamma * (1.0_rt - 1.0_rt / (1.0_rt + gamma));
+
+    df_dT *= dgamma_dT;
+  }
+}
 
 template <int do_T_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
@@ -809,57 +838,22 @@ void chabrier1998 (const plasma_state_t& state,
   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);
 
-  [[maybe_unused]] Real Gamma1dT, Gamma2dT, Gamma12dT;
+  Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{};
 
   if constexpr (do_T_derivatives) {
     Gamma1dT = -Gamma1 / state.temp;
     Gamma2dT = -Gamma2 / state.temp;
     Gamma12dT = -Gamma12 / state.temp;
   }
-  // Fitted parameters, see Chabrier & Potekhin 1998 Sec.IV
-
-  const Real A_1 = -0.9052_rt;
-  const Real A_2 = 0.6322_rt;
-  const Real A_3 = -0.5_rt * std::sqrt(3.0_rt) - A_1 / std::sqrt(A_2);
-
-  // Helmholtz free energy, See Chabrier & Potekhin 1998 Eq. 28
-
-  Real f1 = A_1 * (std::sqrt(Gamma1 * (A_2 + Gamma1)) -
-		   A_2 * std::log(std::sqrt(Gamma1/A_2) + std::sqrt(1.0_rt + Gamma1/A_2)))
-    + 2.0_rt * A_3 * (std::sqrt(Gamma1) - std::atan(std::sqrt(Gamma1)));
-
-  Real f2 = A_1 * (std::sqrt(Gamma2 * (A_2 + Gamma2)) -
-		   A_2 * std::log(std::sqrt(Gamma2/A_2) + std::sqrt(1.0_rt + Gamma2/A_2)))
-    + 2.0_rt * A_3 * (std::sqrt(Gamma2) - std::atan(std::sqrt(Gamma2)));
-
-  Real f12 = A_1 * (std::sqrt(Gamma12 * (A_2 + Gamma12)) -
-		   A_2 * std::log(std::sqrt(Gamma12/A_2) + std::sqrt(1.0_rt + Gamma12/A_2)))
-    + 2.0_rt * A_3 * (std::sqrt(Gamma12) - std::atan(std::sqrt(Gamma12)));
 
-  [[maybe_unused]] Real f1dT, f2dT, f12dT;
-
-  if constexpr (do_T_derivatives) {
-    f1dT = A_1 * ((A_2 + 2.0_rt * Gamma1) / (2.0_rt * std::sqrt(Gamma1 * (A_2 + Gamma1)))
-		       - 0.5_rt / (std::sqrt(Gamma1 / A_2) + std::sqrt(1.0_rt + Gamma1 / A_2))
-		       * (1.0_rt / std::sqrt(Gamma1 / A_2) + 1.0_rt / std::sqrt(1.0_rt + Gamma1 / A_2)))
-      + A_3 / std::sqrt(Gamma1) * (1.0_rt - 1.0_rt / (1.0_rt + Gamma1));
-
-    f1dT *= Gamma1dT;
+  // Helmholtz free energy
 
-    f2dT = A_1 * ((A_2 + 2.0_rt * Gamma2) / (2.0_rt * std::sqrt(Gamma2 * (A_2 + Gamma2)))
-		       - 0.5_rt / (std::sqrt(Gamma2 / A_2) + std::sqrt(1.0_rt + Gamma2 / A_2))
-		       * (1.0_rt / std::sqrt(Gamma2 / A_2) + 1.0_rt / std::sqrt(1.0_rt + Gamma2 / A_2)))
-      + A_3 / std::sqrt(Gamma2) * (1.0_rt - 1.0_rt / (1.0_rt + Gamma2));
+  Real f1, f2, f12;
+  Real f1dT, f2dT, f12dT;
 
-    f2dT *= Gamma2dT;
-
-    f12dT = A_1 * ((A_2 + 2.0_rt * Gamma12) / (2.0_rt * std::sqrt(Gamma12 * (A_2 + Gamma12)))
-		       - 0.5_rt / (std::sqrt(Gamma12 / A_2) + std::sqrt(1.0_rt + Gamma12 / A_2))
-		       * (1.0_rt / std::sqrt(Gamma12 / A_2) + 1.0_rt / std::sqrt(1.0_rt + Gamma12 / A_2)))
-      + A_3 / std::sqrt(Gamma12) * (1.0_rt - 1.0_rt / (1.0_rt + Gamma12));
-
-    f12dT *= Gamma12dT;
-  }
+  chabrier1998_helmholtz_F<do_T_derivatives>(Gamma1, Gamma1dT, f1, f1dT);
+  chabrier1998_helmholtz_F<do_T_derivatives>(Gamma2, Gamma2dT, f2, f2dT);
+  chabrier1998_helmholtz_F<do_T_derivatives>(Gamma12, Gamma12dT, f12, f12dT);
 
   // Now we add quantum correction terms discussed in Alastuey 1978.
   // Notice in Alastuey 1978, they have a different classical term,