From 5a7415fe3c996b49fbd1c7a08c00407ebc45aba7 Mon Sep 17 00:00:00 2001 From: Piyush Sharda <34922596+psharda@users.noreply.github.com> Date: Fri, 21 Jul 2023 18:52:07 +0200 Subject: [PATCH 1/2] Use amu instead of 1/n_A for m_nucleon in gamma_law (#1240) this will cause small diffs for any problem that uses the gamma law EOS --- EOS/gamma_law/actual_eos.H | 4 +-- EOS/primordial_chem/actual_eos.H | 4 +-- .../eos_cell/ci-benchmarks/eos_gamma_law.out | 32 +++++++++---------- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/EOS/gamma_law/actual_eos.H b/EOS/gamma_law/actual_eos.H index d82bf12857..1be67f41da 100644 --- a/EOS/gamma_law/actual_eos.H +++ b/EOS/gamma_law/actual_eos.H @@ -52,8 +52,8 @@ void actual_eos (I input, T& state) { static_assert(std::is_same::value, "input must be an eos_input_t"); - // Get the mass of a nucleon from Avogadro's number. - const Real m_nucleon = 1.0 / C::n_A; + // Get the mass of a nucleon from m_u. + const Real m_nucleon = C::m_u; if constexpr (has_xn::value) { if (eos_assume_neutral) { diff --git a/EOS/primordial_chem/actual_eos.H b/EOS/primordial_chem/actual_eos.H index ec3a6dc6d5..650cb72636 100644 --- a/EOS/primordial_chem/actual_eos.H +++ b/EOS/primordial_chem/actual_eos.H @@ -139,8 +139,8 @@ void actual_eos (I input, T& state) { static_assert(std::is_same::value, "input must be either eos_input_rt or eos_input_re"); - const Real gasconstant = 8.3144725e7; - const Real protonmass = 1.672621637e-24; + const Real gasconstant = C::n_A * C::k_B; + const Real protonmass = C::m_p; // Special gamma factors Real sum_Abarinv = 0.0_rt; Real sum_gammasinv = 0.0_rt; diff --git a/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out b/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out index 076fc4f23c..bf346e3bde 100644 --- a/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out +++ b/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out @@ -1,24 +1,24 @@ -AMReX (22.10-75-g0a3ee1486cc3) initialized +AMReX (23.06-31-g9c256b12b89f) initialized calling the EOS on a single zone state... rho = 1000000 T = 1000000000 xn = 1 0 0 0 0 0 0 0 0 0 0 0 0 -p = 2.078615536e+22 -e = 3.117923305e+16 -h = 5.196538841e+16 -s = 315177040.1 -dpdT = 2.078615536e+13 -dpdr = 2.078615536e+16 -dedT = 31179233.05 +p = 2.078615354e+22 +e = 3.117923031e+16 +h = 5.196538385e+16 +s = 315177017.1 +dpdT = 2.078615354e+13 +dpdr = 2.078615354e+16 +dedT = 31179230.31 dedr = 0 -dhdT = 51965388.41 +dhdT = 51965383.85 dhdr = 0 -dsdT = 0.03117923305 -dsdr = -20.78615536 +dsdT = 0.03117923031 +dsdr = -20.78615354 dpde = 666666.6667 -dpdr_e = 2.078615536e+16 -cv = 31179233.05 -cp = 51965388.41 +dpdr_e = 2.078615354e+16 +cv = 31179230.31 +cp = 51965383.85 xne = 0 xnp = 0 eta = 0 @@ -28,10 +28,10 @@ mu = 4 mu_e = 2 y_e = 0.5 gam1 = 1.666666667 -cs = 186127892.2 +cs = 186127884.1 abar = 4 zbar = 2 Unused ParmParse Variables: [TOP]::unit_test.x4(nvals = 1) :: [0.0] -AMReX (22.10-75-g0a3ee1486cc3) finalized +AMReX (23.06-31-g9c256b12b89f) finalized From 9edc7c1ffd42d84fbe445053a881c76d133cef6f Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Fri, 21 Jul 2023 15:29:55 -0400 Subject: [PATCH 2/2] Renumber the screening methods so null is at index 0 (#1286) --- Make.Microphysics_extern | 8 ++++---- nse_solver/nse_solver.H | 4 ++-- screening/screen.H | 20 ++++++++++---------- sphinx_docs/source/screening.rst | 8 ++++---- 4 files changed, 20 insertions(+), 20 deletions(-) diff --git a/Make.Microphysics_extern b/Make.Microphysics_extern index 2ee4665125..984e461e4c 100644 --- a/Make.Microphysics_extern +++ b/Make.Microphysics_extern @@ -15,13 +15,13 @@ ifeq ($(USE_NONAKA_PLOT),TRUE) endif SCREEN_METHOD ?= screen5 -ifeq ($(SCREEN_METHOD), screen5) +ifeq ($(SCREEN_METHOD), null) DEFINES += -DSCREEN_METHOD=0 -else ifeq ($(SCREEN_METHOD), chugunov2007) +else ifeq ($(SCREEN_METHOD), screen5) DEFINES += -DSCREEN_METHOD=1 -else ifeq ($(SCREEN_METHOD), chugunov2009) +else ifeq ($(SCREEN_METHOD), chugunov2007) DEFINES += -DSCREEN_METHOD=2 -else ifeq ($(SCREEN_METHOD), null) +else ifeq ($(SCREEN_METHOD), chugunov2009) DEFINES += -DSCREEN_METHOD=3 else ifeq ($(SCREEN_METHOD), chabrier1998) DEFINES += -DSCREEN_METHOD=4 diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H index d5117ab67f..b8fadacbaf 100644 --- a/nse_solver/nse_solver.H +++ b/nse_solver/nse_solver.H @@ -89,7 +89,7 @@ void apply_nse_exponent(T& nse_state) { // if we use a non-null screening method. // Get the required terms to calculate coulomb correction term, u_c -#if SCREEN_METHOD != 3 +#if SCREEN_METHOD != 0 // three unit-less constants for calculating coulomb correction term // see Calder 2007, doi:10.1086/510709 for more detail @@ -119,7 +119,7 @@ void apply_nse_exponent(T& nse_state) { // term for calculating u_c // if use non-null screening, calculate the coulomb correction term -#if SCREEN_METHOD != 3 +#if SCREEN_METHOD != 0 gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e; // chemical potential for coulomb correction diff --git a/screening/screen.H b/screening/screen.H index 34599b27a8..868fa1060e 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -17,13 +17,13 @@ using namespace amrex; using namespace integrator_rp; #if SCREEN_METHOD == 0 -const std::string screen_name = "screen5"; +const std::string screen_name = "null"; #elif SCREEN_METHOD == 1 -const std::string screen_name = "chugunov2007"; +const std::string screen_name = "screen5"; #elif SCREEN_METHOD == 2 -const std::string screen_name = "chugunov2009"; +const std::string screen_name = "chugunov2007"; #elif SCREEN_METHOD == 3 -const std::string screen_name = "null"; +const std::string screen_name = "chugunov2009"; #elif SCREEN_METHOD == 4 const std::string screen_name = "chabrier1998"; #endif @@ -977,16 +977,16 @@ void actual_screen(const plasma_state_t& state, Real& scor, Real& scordt) { #if SCREEN_METHOD == 0 - actual_screen5(state, scn_fac, scor, scordt); -#elif SCREEN_METHOD == 1 - chugunov2007(state, scn_fac, scor, scordt); -#elif SCREEN_METHOD == 2 - chugunov2009(state, scn_fac, scor, scordt); -#elif SCREEN_METHOD == 3 // null screening amrex::ignore_unused(state, scn_fac); scor = 1.0_rt; scordt = 0.0_rt; +#elif SCREEN_METHOD == 1 + actual_screen5(state, scn_fac, scor, scordt); +#elif SCREEN_METHOD == 2 + chugunov2007(state, scn_fac, scor, scordt); +#elif SCREEN_METHOD == 3 + chugunov2009(state, scn_fac, scor, scordt); #elif SCREEN_METHOD == 4 chabrier1998(state, scn_fac, scor, scordt); #endif diff --git a/sphinx_docs/source/screening.rst b/sphinx_docs/source/screening.rst index fb5e837f1b..e54cd1a878 100644 --- a/sphinx_docs/source/screening.rst +++ b/sphinx_docs/source/screening.rst @@ -40,11 +40,11 @@ The options are: This includes the portion in the appendix that blends in the weak screening limit. +* ``chabrier1998``: + + This implements the screening of :cite:`Chabrier_1998` as well as the quantum corrections for strong screening according to screen5. This is suggested in the appendix of :cite:`Calder_2007`. This screening is compatible with NSE calculations unlike ``screen5``, ``chugunov2007``, and ``chugunov2009``. This screening should be valid in the weak screening regime, :math:`\Gamma < 0.1`, and strong screening regime, :math:`1 \lesssim \Gamma \lesssim 160`. + * ``null`` : This disables screening by always returning 1 for the screening enhancement factor. - -* ``chabrier1998``: - - This implements the screening of :cite:`Chabrier_1998` as well as the quantum corrections for strong screening according to screen5. This is suggested in the appendix of :cite:`Calder_2007`. This screening is compatible with NSE calculations unlike ``screen5``, ``chugunov2007``, and ``chugunov2009``. This screening should be valid in the weak screening regime, :math:`\Gamma < 0.1`, and strong screening regime, :math:`1 \lesssim \Gamma \lesssim 160`.