Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Move Helmholtz free energy from chabrier1998 into a helper function #1334

Merged
merged 1 commit into from
Sep 13, 2023
Merged

Move Helmholtz free energy from chabrier1998 into a helper function #1334

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 8 additions & 14 deletions nse_solver/nse_solver.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include <eos_composition.H>
#include <microphysics_sort.H>
#include <hybrj.H>
#include <screen.H>
#include <cctype>
#include <algorithm>

Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
86 changes: 40 additions & 46 deletions screening/screen.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down