Skip to content

Commit

Permalink
Merge branch 'development' into nse_second_order_refactor
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale authored Nov 29, 2023
2 parents 4d5df40 + f104eff commit e4e143d
Show file tree
Hide file tree
Showing 4 changed files with 99 additions and 37 deletions.
12 changes: 6 additions & 6 deletions Make.Microphysics_extern
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,15 +23,15 @@ endif

SCREEN_METHOD ?= screen5
ifeq ($(SCREEN_METHOD), null)
DEFINES += -DSCREEN_METHOD=0
DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_null
else ifeq ($(SCREEN_METHOD), screen5)
DEFINES += -DSCREEN_METHOD=1
DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_screen5
else ifeq ($(SCREEN_METHOD), chugunov2007)
DEFINES += -DSCREEN_METHOD=2
DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_chugunov2007
else ifeq ($(SCREEN_METHOD), chugunov2009)
DEFINES += -DSCREEN_METHOD=3
DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_chugunov2009
else ifeq ($(SCREEN_METHOD), chabrier1998)
DEFINES += -DSCREEN_METHOD=4
DEFINES += -DSCREEN_METHOD=SCREEN_METHOD_chabrier1998
else
$(error Invalid value for SCREEN_METHOD)
endif
Expand Down Expand Up @@ -155,4 +155,4 @@ endif
clean::
@if [ -L helm_table.dat ]; then rm -f helm_table.dat; fi
@if [ -L reaclib_rate_metadata.dat ]; then rm -f reaclib_rate_metadata.dat; fi
$(foreach t, $(wildcard *_betadecay.dat *_electroncapture.dat nse*.tbl), $(shell if [ -L $t ]; then rm -f $t; fi))
$(foreach t, $(wildcard *_betadecay.dat *_electroncapture.dat nse*.tbl), $(shell if [ -L $t ]; then rm -f $t; fi))
4 changes: 2 additions & 2 deletions nse_solver/nse_solver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ void apply_nse_exponent(T& nse_state) {
// if we use chabrier1998 screening
// Get the required terms to calculate coulomb correction term, u_c

#if SCREEN_METHOD == 4
#if SCREEN_METHOD == SCREEN_METHOD_chabrier1998

// Find n_e for original state;
const amrex::Real n_e = nse_state.rho * nse_state.y_e / C::m_u;
Expand All @@ -113,7 +113,7 @@ void apply_nse_exponent(T& nse_state) {
// term for calculating u_c

// if use chabrier1998 screening, calculate the coulomb correction term
#if SCREEN_METHOD == 4
#if SCREEN_METHOD == SCREEN_METHOD_chabrier1998
gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e;

// chemical potential for coulomb correction
Expand Down
49 changes: 30 additions & 19 deletions screening/screen.H
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,13 @@
#ifndef SCREEN_H
#define SCREEN_H

// these need to be defined before screen_data.H is included
#define SCREEN_METHOD_null 0
#define SCREEN_METHOD_screen5 1
#define SCREEN_METHOD_chugunov2007 2
#define SCREEN_METHOD_chugunov2009 3
#define SCREEN_METHOD_chabrier1998 4

#include <AMReX.H>
#include <AMReX_Algorithm.H>
#include <AMReX_Array.H>
Expand All @@ -16,15 +23,15 @@
using namespace amrex;
using namespace integrator_rp;

#if SCREEN_METHOD == 0
#if SCREEN_METHOD == SCREEN_METHOD_null
const std::string screen_name = "null";
#elif SCREEN_METHOD == 1
#elif SCREEN_METHOD == SCREEN_METHOD_screen5
const std::string screen_name = "screen5";
#elif SCREEN_METHOD == 2
#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007
const std::string screen_name = "chugunov2007";
#elif SCREEN_METHOD == 3
#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009
const std::string screen_name = "chugunov2009";
#elif SCREEN_METHOD == 4
#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998
const std::string screen_name = "chabrier1998";
#endif

Expand Down Expand Up @@ -169,6 +176,7 @@ fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Array
}
}

#if SCREEN_METHOD == SCREEN_METHOD_screen5
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 +424,7 @@ void actual_screen5 (const plasma_state_t& state,
}
}

#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007
template <int do_T_derivatives>
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void chugunov2007 (const plasma_state_t& state,
Expand Down Expand Up @@ -460,7 +469,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 +620,7 @@ void chugunov2007 (const plasma_state_t& state,
}
}

#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009
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 +708,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.z2_53;
Real Gamma_comp = Gamma_e * scn_fac.zs53;

Real Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde;

Expand Down Expand Up @@ -753,7 +763,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 +790,7 @@ void chugunov2009 (const plasma_state_t& state,
}
}

#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998
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 +841,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,25 +944,26 @@ void chabrier1998 (const plasma_state_t& state,
}

}
#endif

template <int do_T_derivatives>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void actual_screen(const plasma_state_t& state,
const scrn::screen_factors_t& scn_fac,
Real& scor, Real& scordt)
{
#if SCREEN_METHOD == 0
#if SCREEN_METHOD == SCREEN_METHOD_null
// null screening
amrex::ignore_unused(state, scn_fac);
scor = 1.0_rt;
scordt = 0.0_rt;
#elif SCREEN_METHOD == 1
#elif SCREEN_METHOD == SCREEN_METHOD_screen5
actual_screen5<do_T_derivatives>(state, scn_fac, scor, scordt);
#elif SCREEN_METHOD == 2
#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007
chugunov2007<do_T_derivatives>(state, scn_fac, scor, scordt);
#elif SCREEN_METHOD == 3
#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009
chugunov2009<do_T_derivatives>(state, scn_fac, scor, scordt);
#elif SCREEN_METHOD == 4
#elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998
chabrier1998<do_T_derivatives>(state, scn_fac, scor, scordt);
#endif
}
Expand Down
71 changes: 61 additions & 10 deletions screening/screen_data.H
View file
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 == SCREEN_METHOD_screen5
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 == SCREEN_METHOD_chugunov2007
amrex::Real ztilde = 0.0;
amrex::Real ztilde3 = 0.0;
#elif SCREEN_METHOD == SCREEN_METHOD_chugunov2009
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 == SCREEN_METHOD_chabrier1998
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 == SCREEN_METHOD_screen5
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 == SCREEN_METHOD_chugunov2007
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 == SCREEN_METHOD_chugunov2009
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 == SCREEN_METHOD_chabrier1998
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

0 comments on commit e4e143d

Please sign in to comment.