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

Factor out some sqrt terms in chabrier1998 screening. #1585

Merged
merged 4 commits into from
Jun 21, 2024
Merged

Factor out some sqrt terms in chabrier1998 screening. #1585

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f361dc3
factor out some terms to have less sqrt calls
zhichen3 Jun 20, 2024
d4c0184
update benchmark
zhichen3 Jun 20, 2024
a5fa160
update benchmnark again
zhichen3 Jun 20, 2024
c86a060
give CBRT_2 its own name
zhichen3 Jun 20, 2024
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
16 changes: 10 additions & 6 deletions screening/screen.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)) +
Expand Down Expand Up @@ -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;
Expand All @@ -883,8 +886,9 @@ 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;
constexpr amrex::Real CBRT_2 = gcem::pow(2.0_rt, 1.0_rt/3.0_rt);
amrex::Real Gamma_eff = CBRT_2 * scn_fac.z1 * scn_fac.z2 *
scn_fac.zs13inv * Gamma_e;
[[maybe_unused]] amrex::Real Gamma_eff_dT;

if constexpr (do_T_derivatives) {
Expand Down
Loading
Loading