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

organize screening, put z1 z2 related pow vars to screen_factors_t #1397

Merged
merged 8 commits into from
Nov 28, 2023
Merged

organize screening, put z1 z2 related pow vars to screen_factors_t #1397

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
7b3295c
put z1z2 related vars in screen_factors_t
zhichen3 Nov 27, 2023
7ee37b8
cleanup
zhichen3 Nov 27, 2023
41945ac
remove tab
zhichen3 Nov 27, 2023
9230bd8
fix var name
zhichen3 Nov 27, 2023
f6d5503
add header to screening functions
zhichen3 Nov 27, 2023
7c39914
Merge branch 'development' into reorganize_screening
zingale Nov 28, 2023
ef093e6
remove ztilde for screen5
zhichen3 Nov 28, 2023
486fffd
Merge branch 'development' into reorganize_screening
zhichen3 Nov 28, 2023
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: 13 additions & 9 deletions screening/screen.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,7 @@ fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Array
}
}

#if SCREEN_METHOD == 1
template <int do_T_derivatives>
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void actual_screen5 (const plasma_state_t& state,
Expand Down Expand Up @@ -416,6 +417,7 @@ void actual_screen5 (const plasma_state_t& state,
}
}

#elif SCREEN_METHOD == 2
template <int do_T_derivatives>
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void chugunov2007 (const plasma_state_t& state,
Expand Down Expand Up @@ -460,7 +462,7 @@ void chugunov2007 (const plasma_state_t& state,
// m_i -> m_u * abar
Real mu12 = scn_fac.a1 * scn_fac.a2 / (scn_fac.a1 + scn_fac.a2);
Real z_factor = scn_fac.z1 * scn_fac.z2;
Real n_i = state.n_e / gcem::pow(scn_fac.ztilde, 3);
Real n_i = state.n_e / scn_fac.ztilde3;
Real m_i = 2.0_rt * mu12 / C::n_A;

constexpr Real T_p_factor = C::hbar/C::k_B*C::q_e*gcem::sqrt(4.0_rt*GCEM_PI);
Expand Down Expand Up @@ -611,6 +613,7 @@ void chugunov2007 (const plasma_state_t& state,
}
}

#elif SCREEN_METHOD == 3
template <int do_T_derivatives>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void chugunov2009_f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT)
Expand Down Expand Up @@ -698,9 +701,9 @@ void chugunov2009 (const plasma_state_t& state,
Real dlog_Gamma_dT = -1.0_rt / state.temp;

// Coulomb coupling parameters for ions and compound nucleus, eqs. 7 & 9
Real Gamma_1 = Gamma_e * std::pow(scn_fac.z1, 5.0_rt/3.0_rt);
Real Gamma_2 = Gamma_e * std::pow(scn_fac.z2, 5.0_rt/3.0_rt);
Real Gamma_comp = Gamma_e * std::pow(zcomp, 5.0_rt/3.0_rt);
Real Gamma_1 = Gamma_e * scn_fac.z1_53;
Real Gamma_2 = Gamma_e * scn_fac.z1_53;
Real Gamma_comp = Gamma_e * scn_fac.zs53;

Real Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde;

Expand Down Expand Up @@ -753,7 +756,7 @@ void chugunov2009 (const plasma_state_t& state,

// weak screening correction term, eq. A3
Real corr_C = 3.0_rt*z1z2 * std::sqrt(state.z2bar/state.zbar) /
(std::pow(zcomp, 2.5_rt) - std::pow(scn_fac.z1, 2.5_rt) - std::pow(scn_fac.z2, 2.5_rt));
(scn_fac.zs52 - scn_fac.z1_52 - scn_fac.z2_52);

// corrected enhancement factor, eq. A4
Real Gamma_12_2 = Gamma_12 * Gamma_12;
Expand All @@ -780,6 +783,7 @@ void chugunov2009 (const plasma_state_t& state,
}
}

#elif SCREEN_METHOD == 4
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) {
Expand Down Expand Up @@ -830,13 +834,12 @@ void chabrier1998 (const plasma_state_t& state,
// Eq. 2 in Chabrier & Potekhin 1998

Real Gamma_e = state.gamma_e_fac / state.temp;
Real zcomp = scn_fac.z1 + scn_fac.z2;

// See Calder2007 appendix Eq. A9

Real Gamma1 = Gamma_e * std::pow(scn_fac.z1, 5.0_rt/3.0_rt);
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);
Real Gamma1 = Gamma_e * scn_fac.z1_53;
Real Gamma2 = Gamma_e * scn_fac.z2_53;
Real Gamma12 = Gamma_e * scn_fac.zs53;

Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{};

Expand Down Expand Up @@ -934,6 +937,7 @@ void chabrier1998 (const plasma_state_t& state,
}

}
#endif

template <int do_T_derivatives>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand Down
71 changes: 61 additions & 10 deletions screening/screen_data.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,20 +17,50 @@ namespace scrn {
amrex::Real a1 = -1;
amrex::Real a2 = -1;

// z1_53 = (z1)**(5./3.)
// z2_53 = (z2)**(5./3.)
// zs52 = (z1+z2)**(5./2.)
// z1_52 = (z1)**(5./2.)
// z2_52 = (z2)**(5./2.)
// zs13 = (z1+z2)**(1./3.)
// zs53 = (z1+z2)**(5./3.)
// zhat = combination of z1 and z2 raised to the 5/3 power
// zhat2 = combination of z1 and z2 raised to the 5/12 power
// lzav = log of effective charge
// aznut = combination of a1,z1,a2,z2 raised to 1/3 power
// ztilde = effective ion radius factor for a MCP
// ztilde3 = ztilde**3

#if SCREEN_METHOD == 1
amrex::Real zs53 = 0.0;
amrex::Real z1_53 = 0.0;
amrex::Real z2_53 = 0.0;
amrex::Real zs13 = 0.0;
amrex::Real zs13inv = 0.0;
amrex::Real zhat = 0.0;
amrex::Real zhat2 = 0.0;
amrex::Real lzav = 0.0;
amrex::Real aznut = 0.0;
#elif SCREEN_METHOD == 2
amrex::Real ztilde = 0.0;
amrex::Real ztilde3 = 0.0;
#elif SCREEN_METHOD == 3
amrex::Real zs52 = 0.0;
amrex::Real z1_52 = 0.0;
amrex::Real z2_52 = 0.0;
amrex::Real zs53 = 0.0;
amrex::Real z1_53 = 0.0;
amrex::Real z2_53 = 0.0;
amrex::Real aznut = 0.0;
amrex::Real ztilde = 0.0;
#elif SCREEN_METHOD == 4
amrex::Real zs53 = 0.0;
amrex::Real z1_53 = 0.0;
amrex::Real z2_53 = 0.0;
amrex::Real zs13 = 0.0;
amrex::Real zs13inv = 0.0;
amrex::Real aznut = 0.0;
#endif

[[nodiscard]]
bool validate_nuclei(const amrex::Real z1_pass, const amrex::Real a1_pass,
Expand All @@ -57,25 +87,46 @@ namespace scrn {
scn_fac.z2 = z2;
scn_fac.a2 = a2;

#if SCREEN_METHOD == 1
scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt);
scn_fac.z1_53 = gcem::pow(z1, 5.0_rt / 3.0_rt);
scn_fac.z2_53 = gcem::pow(z2, 5.0_rt / 3.0_rt);
scn_fac.zs13 = gcem::pow(z1 + z2, 1.0_rt / 3.0_rt);

scn_fac.zs13inv = 1.0_rt / scn_fac.zs13;

scn_fac.zhat = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt) -
gcem::pow(z1, 5.0_rt/3.0_rt) -
gcem::pow(z2, 5.0_rt/3.0_rt);

scn_fac.zhat = scn_fac.zs53 - scn_fac.z1_53 - scn_fac.z2_53;
scn_fac.zhat2 = gcem::pow(z1 + z2, 5.0_rt / 12.0_rt) -
gcem::pow(z1, 5.0_rt/12.0_rt) -
gcem::pow(z2, 5.0_rt/12.0_rt);

gcem::pow(z1, 5.0_rt / 12.0_rt) -
gcem::pow(z2, 5.0_rt / 12.0_rt);
scn_fac.lzav = (5.0_rt / 3.0_rt) * gcem::log(z1 * z2 / (z1 + z2));

scn_fac.aznut = gcem::pow(z1 * z1 * z2 * z2 * a1 * a2 / (a1 + a2),
1.0_rt / 3.0_rt);

#elif SCREEN_METHOD == 2
scn_fac.ztilde = 0.5_rt * (gcem::pow(z1, 1.0_rt / 3.0_rt) +
gcem::pow(z2, 1.0_rt / 3.0_rt));
scn_fac.ztilde3 = gcem::pow(scn_fac.ztilde, 3);

#elif SCREEN_METHOD == 3
scn_fac.zs52 = gcem::pow(z1 + z2, 5.0_rt / 2.0_rt);
scn_fac.z1_52 = gcem::pow(z1, 5.0_rt / 2.0_rt);
scn_fac.z2_52 = gcem::pow(z2, 5.0_rt / 2.0_rt);
scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt);
scn_fac.z1_53 = gcem::pow(z1, 5.0_rt / 3.0_rt);
scn_fac.z2_53 = gcem::pow(z2, 5.0_rt / 3.0_rt);
scn_fac.aznut = gcem::pow(z1 * z1 * z2 * z2 * a1 * a2 / (a1 + a2),
1.0_rt / 3.0_rt);
scn_fac.ztilde = 0.5_rt * (gcem::pow(z1, 1.0_rt / 3.0_rt) +
gcem::pow(z2, 1.0_rt / 3.0_rt));

#elif SCREEN_METHOD == 4
scn_fac.zs53 = gcem::pow(z1 + z2, 5.0_rt / 3.0_rt);
scn_fac.z1_53 = gcem::pow(z1, 5.0_rt/3.0_rt);
scn_fac.z2_53 = gcem::pow(z2, 5.0_rt/3.0_rt);
scn_fac.zs13 = gcem::pow(z1 + z2, 1.0_rt / 3.0_rt);
scn_fac.zs13inv = 1.0_rt / scn_fac.zs13;
scn_fac.aznut = gcem::pow(z1 * z1 * z2 * z2 * a1 * a2 / (a1 + a2),
1.0_rt / 3.0_rt);
#endif

return scn_fac;
}
Expand Down