diff --git a/screening/_parameters b/screening/_parameters new file mode 100644 index 0000000000..3e6d41cc63 --- /dev/null +++ b/screening/_parameters @@ -0,0 +1,3 @@ +@namespace: screening + +enable_chabrier1998_quantum_corr integer 0 diff --git a/screening/screen.H b/screening/screen.H index 0cb5edb28a..7f79f6fb14 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -22,6 +22,7 @@ using namespace amrex; using namespace integrator_rp; +using namespace screening_rp; #if SCREEN_METHOD == SCREEN_METHOD_null const std::string screen_name = "null"; @@ -37,40 +38,40 @@ const std::string screen_name = "chabrier1998"; struct plasma_state_t { - Real qlam0z; - Real qlam0zdt; - //Real qlam0zdd; + Real qlam0z; + Real qlam0zdt; + //Real qlam0zdd; - Real taufac; - Real taufacdt; + Real taufac; + Real taufacdt; - Real aa; - Real daadt; - //Real daadd; + Real aa; + Real daadt; + //Real daadd; - Real temp; - Real zbar; - Real z2bar; - Real n_e; - Real gamma_e_fac; + Real temp; + Real zbar; + Real z2bar; + Real n_e; + Real gamma_e_fac; }; inline std::ostream& operator<< (std::ostream& o, plasma_state_t const& pstate) { - o << "qlam0z = " << pstate.qlam0z << std::endl; - o << "qlam0zdt = " << pstate.qlam0zdt << std::endl; - o << "taufac = " << pstate.taufac << std::endl; - o << "taufacdt = " << pstate.taufacdt << std::endl; - o << "aa = " << pstate.aa << std::endl; - o << "daadt = " << pstate.daadt << std::endl; - o << "temp = " << pstate.temp << std::endl; - o << "zbar = " << pstate.zbar << std::endl; - o << "z2bar = " << pstate.z2bar << std::endl; - o << "n_e = " << pstate.n_e << std::endl; - o << "gamma_e_fac = " << pstate.gamma_e_fac << std::endl; - - return o; + o << "qlam0z = " << pstate.qlam0z << std::endl; + o << "qlam0zdt = " << pstate.qlam0zdt << std::endl; + o << "taufac = " << pstate.taufac << std::endl; + o << "taufacdt = " << pstate.taufacdt << std::endl; + o << "aa = " << pstate.aa << std::endl; + o << "daadt = " << pstate.daadt << std::endl; + o << "temp = " << pstate.temp << std::endl; + o << "zbar = " << pstate.zbar << std::endl; + o << "z2bar = " << pstate.z2bar << std::endl; + o << "n_e = " << pstate.n_e << std::endl; + o << "gamma_e_fac = " << pstate.gamma_e_fac << std::endl; + + return o; } AMREX_FORCE_INLINE @@ -88,85 +89,87 @@ screening_finalize() { template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Array1D const& y) { - - Real sum = 0.0_rt; - for (int n = 1; n <= NumSpec; n++) { - sum += y(n); - } - Real abar = 1.0_rt / sum; - Real ytot = sum; - - sum = 0.0_rt; - Real sum2 = 0.0_rt; - for (int n = 1; n <= NumSpec; n++) { - sum += zion[n-1]*y(n); - sum2 += zion[n-1]*zion[n-1]*y(n); - } - - // Part of Eq.6 in Itoh:1979 - // 4.248719e3 = (27*pi^2*e^4*m_u/(2*k_B*hbar^2))^(1/3) - // the extra (1/3) to make tau -> tau/3 - const Real co2 = (1.0_rt/3.0_rt) * 4.248719e3_rt; - - Real zbar = sum * abar; - Real z2bar = sum2 * abar; - - // ntot - Real rr = dens * ytot; - Real tempi = 1.0_rt / temp; - [[maybe_unused]] Real dtempi; - if constexpr (do_T_derivatives) { - dtempi = -tempi * tempi; - } - - // Part of Eq. 19 in Graboske:1973 - // pp = sqrt( \tilde{z}*(rho/u_I/T) ) - Real pp = std::sqrt(rr*tempi*(z2bar + zbar)); - [[maybe_unused]] Real dppdt; - if constexpr (do_T_derivatives) { - Real qq = 0.5_rt/pp *(z2bar + zbar); - dppdt = qq*rr*dtempi; - } - - // Part version of Eq. 19 in Graboske:1973 - state.qlam0z = 1.88e8_rt * tempi * pp; - if constexpr (do_T_derivatives) { - state.qlam0zdt = 1.88e8_rt * (dtempi*pp + tempi*dppdt); - } - - // Part of Eq.6 in Itoh:1979 - state.taufac = co2 * std::cbrt(tempi); - if constexpr (do_T_derivatives) { - state.taufacdt = -(1.0_rt/3.0_rt) * state.taufac * tempi; - } - - Real xni = std::cbrt(rr * zbar); - - // Part of Eq.4 in Itoh:1979 - // 2.27493e5 = e^2 / ( (3*m_u/(4pi))^(1/3) *k_B ) - state.aa = 2.27493e5_rt * tempi * xni; - if constexpr (do_T_derivatives) { - state.daadt = 2.27493e5_rt * dtempi * xni; - } - - state.temp = temp; - state.zbar = zbar; - state.z2bar = z2bar; - - // Electron number density - // zbar * ntot works out to sum(z[i] * n[i]), after cancelling terms - state.n_e = zbar * rr * C::n_A; - - // precomputed part of Gamma_e, from Chugunov 2009 eq. 6 - constexpr Real gamma_e_constants = - C::q_e*C::q_e/C::k_B * gcem::pow(4.0_rt/3.0_rt*M_PI, 1.0_rt/3.0_rt); - state.gamma_e_fac = gamma_e_constants * std::cbrt(state.n_e); +fill_plasma_state(plasma_state_t& state, const Real temp, + const Real dens, Array1D const& y) { + + Real sum = 0.0_rt; + for (int n = 1; n <= NumSpec; n++) { + sum += y(n); + } + Real abar = 1.0_rt / sum; + Real ytot = sum; + + sum = 0.0_rt; + Real sum2 = 0.0_rt; + for (int n = 1; n <= NumSpec; n++) { + sum += zion[n-1]*y(n); + sum2 += zion[n-1]*zion[n-1]*y(n); + } + + // Part of Eq.6 in Itoh:1979 + // 4.248719e3 = (27*pi^2*e^4*m_u/(2*k_B*hbar^2))^(1/3) + // the extra (1/3) to make tau -> tau/3 + const Real co2 = (1.0_rt/3.0_rt) * 4.248719e3_rt; + + Real zbar = sum * abar; + Real z2bar = sum2 * abar; + + // ntot + Real rr = dens * ytot; + Real tempi = 1.0_rt / temp; + [[maybe_unused]] Real dtempi; + if constexpr (do_T_derivatives) { + dtempi = -tempi * tempi; + } + + // Part of Eq. 19 in Graboske:1973 + // pp = sqrt( \tilde{z}*(rho/u_I/T) ) + Real pp = std::sqrt(rr*tempi*(z2bar + zbar)); + [[maybe_unused]] Real dppdt; + if constexpr (do_T_derivatives) { + Real qq = 0.5_rt/pp *(z2bar + zbar); + dppdt = qq*rr*dtempi; + } + + // Part version of Eq. 19 in Graboske:1973 + state.qlam0z = 1.88e8_rt * tempi * pp; + if constexpr (do_T_derivatives) { + state.qlam0zdt = 1.88e8_rt * (dtempi*pp + tempi*dppdt); + } + + // Part of Eq.6 in Itoh:1979 + state.taufac = co2 * std::cbrt(tempi); + if constexpr (do_T_derivatives) { + state.taufacdt = -(1.0_rt/3.0_rt) * state.taufac * tempi; + } + + Real xni = std::cbrt(rr * zbar); + + // Part of Eq.4 in Itoh:1979 + // 2.27493e5 = e^2 / ( (3*m_u/(4pi))^(1/3) *k_B ) + state.aa = 2.27493e5_rt * tempi * xni; + if constexpr (do_T_derivatives) { + state.daadt = 2.27493e5_rt * dtempi * xni; + } + + state.temp = temp; + state.zbar = zbar; + state.z2bar = z2bar; + + // Electron number density + // zbar * ntot works out to sum(z[i] * n[i]), after cancelling terms + state.n_e = zbar * rr * C::n_A; + + // precomputed part of Gamma_e, from Chugunov 2009 eq. 6 + constexpr Real gamma_e_constants = + C::q_e*C::q_e/C::k_B * gcem::pow(4.0_rt/3.0_rt*M_PI, 1.0_rt/3.0_rt); + state.gamma_e_fac = gamma_e_constants * std::cbrt(state.n_e); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Array1D const& y) { +fill_plasma_state(plasma_state_t& state, const Real temp, + const Real dens, Array1D const& y) { if (jacobian == 1) { constexpr int do_T_derivatives = 1; fill_plasma_state(state, temp, dens, y); @@ -715,9 +718,11 @@ void chugunov2009 (const plasma_state_t& state, Real Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde; // Coulomb barrier penetrability, eq. 10 + constexpr Real tau_factor = gcem::pow( 27.0_rt/2.0_rt * amrex::Math::powi<2>(M_PI*C::q_e*C::q_e/C::hbar) / (C::n_A*C::k_B), 1.0_rt/3.0_rt); + Real tau_12 = tau_factor * scn_fac.aznut / std::cbrt(state.temp); Real dlog_tau_12_dT; if constexpr (do_T_derivatives) { @@ -742,7 +747,8 @@ void chugunov2009 (const plasma_state_t& state, Real dlog_dT = 0.0_rt; if constexpr (do_T_derivatives) { Real dc3_dT = -1.8_rt / Gamma_12 * dlog_Gamma_dT; - Real dpoly_dT = (c1 + zeta*(2.0_rt*c2 + 3.0_rt*c3*zeta))*dzeta_dT + dc3_dT*zeta*zeta*zeta; + Real dpoly_dT = (c1 + zeta*(2.0_rt*c2 + 3.0_rt*c3*zeta))*dzeta_dT + + dc3_dT*zeta*zeta*zeta; Real dlog_t_12_dT = dpoly_dT / (3.0_rt * poly); dlog_dT = dlog_Gamma_dT - dlog_t_12_dT; } @@ -793,32 +799,33 @@ void chugunov2009 (const plasma_state_t& state, #elif SCREEN_METHOD == SCREEN_METHOD_chabrier1998 template 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; - } +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 @@ -827,122 +834,131 @@ void chabrier1998 (const plasma_state_t& state, const scrn::screen_factors_t& scn_fac, Real& scor, Real& scordt) { - // Calculates screening factors based on Chabrier & Potekhin 1998, - // Calder2007 and partly screen5 routine mentioned in Alastuey 1978. + // Calculates screening factors based on Chabrier & Potekhin 1998, + // Calder2007 and partly screen5 routine mentioned in Alastuey 1978. - // This screening is valid for weak screening: Gamma < 0.1 - // and strong screening: 1 <= Gamma <= 160 - // Reference: - // Chabrier & Potekhin 1998, PhRvE, 58, 4941 - // Calder 2007, doi:10.1086/510709 - // Wallace & Woosley 1982 - // Alastuey 1978 + // This screening is valid for weak screening: Gamma < 0.1 + // and strong screening: 1 <= Gamma <= 160 + // Reference: + // Chabrier & Potekhin 1998, PhRvE, 58, 4941 + // Calder 2007, doi:10.1086/510709 + // Wallace & Woosley 1982 + // Alastuey 1978 - // Eq. 2 in Chabrier & Potekhin 1998 + // Eq. 2 in Chabrier & Potekhin 1998 - Real Gamma_e = state.gamma_e_fac / state.temp; + Real Gamma_e = state.gamma_e_fac / state.temp; - // See Calder2007 appendix Eq. A9 + // See Calder2007 appendix Eq. A9 - Real Gamma1 = Gamma_e * scn_fac.z1_53; - Real Gamma2 = Gamma_e * scn_fac.z2_53; - Real Gamma12 = Gamma_e * scn_fac.zs53; + 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{}; + Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{}; - if constexpr (do_T_derivatives) { - Gamma1dT = -Gamma1 / state.temp; - Gamma2dT = -Gamma2 / state.temp; - Gamma12dT = -Gamma12 / state.temp; - } + if constexpr (do_T_derivatives) { + Gamma1dT = -Gamma1 / state.temp; + Gamma2dT = -Gamma2 / state.temp; + Gamma12dT = -Gamma12 / state.temp; + } - // Helmholtz free energy + // Helmholtz free energy - Real f1, f2, f12; - Real f1dT, f2dT, f12dT; + Real f1, f2, f12; + Real f1dT, f2dT, f12dT; - chabrier1998_helmholtz_F(Gamma1, Gamma1dT, f1, f1dT); - chabrier1998_helmholtz_F(Gamma2, Gamma2dT, f2, f2dT); - chabrier1998_helmholtz_F(Gamma12, Gamma12dT, f12, f12dT); + chabrier1998_helmholtz_F(Gamma1, Gamma1dT, f1, f1dT); + chabrier1998_helmholtz_F(Gamma2, Gamma2dT, f2, f2dT); + chabrier1998_helmholtz_F(Gamma12, Gamma12dT, f12, f12dT); - // 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. + // 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. + Real quantum_corr_1 = 0.0_rt; - // See Wallace1982, Eq. A13 + Real quantum_corr_2 = 0.0_rt; - Real Gamma_eff = std::cbrt(2.0_rt) * scn_fac.z1 * scn_fac.z2 * scn_fac.zs13inv * Gamma_e; - [[maybe_unused]] Real Gamma_eff_dT; + [[maybe_unused]] Real quantum_corr_1_dT = 0.0_rt; + [[maybe_unused]] Real quantum_corr_2_dT = 0.0_rt; - if constexpr (do_T_derivatives) { - Gamma_eff_dT = -Gamma_eff / state.temp; - } + if (enable_chabrier1998_quantum_corr) { + // See Wallace1982, Eq. A13 - // TAU/3, see Wallace1982, Eq. A2 + Real Gamma_eff = std::cbrt(2.0_rt) * scn_fac.z1 * scn_fac.z2 * + scn_fac.zs13inv * Gamma_e; + [[maybe_unused]] Real Gamma_eff_dT; - Real tau12 = state.taufac * scn_fac.aznut; + if constexpr (do_T_derivatives) { + Gamma_eff_dT = -Gamma_eff / state.temp; + } - [[maybe_unused]] Real tau12dT; - if constexpr (do_T_derivatives) { - tau12dT = state.taufacdt * scn_fac.aznut; - } + // TAU/3, see Wallace1982, Eq. A2 - // see Calder 2007 Eq. A8 + Real tau12 = state.taufac * scn_fac.aznut; - Real b_fac = Gamma_eff / tau12; + [[maybe_unused]] Real tau12dT; + if constexpr (do_T_derivatives) { + tau12dT = state.taufacdt * scn_fac.aznut; + } - [[maybe_unused]] Real b_fac_dT; - if constexpr (do_T_derivatives) { - b_fac_dT = (Gamma_eff_dT - b_fac * tau12dT) / tau12; - } + // see Calder 2007 Eq. A8 - // Quantum correction terms (same as screen5) - //see Calder 2007 Eq.A8 and Alastuey1978, Eq. 24 and 31 + Real b_fac = Gamma_eff / tau12; - Real quantum_corr_1 = -tau12 * (5.0_rt/32.0_rt * amrex::Math::powi<3>(b_fac) - 0.014_rt * amrex::Math::powi<4>(b_fac) - - 0.128_rt * amrex::Math::powi<5>(b_fac)); + [[maybe_unused]] Real b_fac_dT; + if constexpr (do_T_derivatives) { + b_fac_dT = (Gamma_eff_dT - b_fac * tau12dT) / tau12; + } - Real quantum_corr_2 = -Gamma_eff * (0.0055_rt * amrex::Math::powi<4>(b_fac) - 0.0098_rt * amrex::Math::powi<5>(b_fac) - + 0.0048_rt * amrex::Math::powi<6>(b_fac)); + // Quantum correction terms (same as screen5) + //see Calder 2007 Eq.A8 and Alastuey1978, Eq. 24 and 31 - [[maybe_unused]] Real quantum_corr_1_dT; - [[maybe_unused]] Real quantum_corr_2_dT; - if constexpr (do_T_derivatives) { - quantum_corr_1_dT = tau12dT / tau12 * quantum_corr_1 - tau12 * b_fac_dT - * (15.0_rt/32.0_rt * amrex::Math::powi<2>(b_fac) - 0.014_rt * 4.0_rt * amrex::Math::powi<3>(b_fac) - - 0.128_rt * 5.0_rt * amrex::Math::powi<4>(b_fac)); + quantum_corr_1 = -tau12 * (5.0_rt/32.0_rt * amrex::Math::powi<3>(b_fac) - + 0.014_rt * amrex::Math::powi<4>(b_fac) - + 0.128_rt * amrex::Math::powi<5>(b_fac)); - quantum_corr_2_dT = Gamma_eff_dT / Gamma_eff * quantum_corr_2 - Gamma_eff * b_fac_dT - * (0.0055_rt * 4.0_rt * amrex::Math::powi<3>(b_fac) - 0.0098_rt * 5.0_rt - * amrex::Math::powi<4>(b_fac) + 0.0048_rt * 6.0_rt * amrex::Math::powi<5>(b_fac)); - } + quantum_corr_2 = -Gamma_eff * (0.0055_rt * amrex::Math::powi<4>(b_fac) - + 0.0098_rt * amrex::Math::powi<5>(b_fac) + + 0.0048_rt * amrex::Math::powi<6>(b_fac)); - // See Calder2007 Appendix Eq. A8. - // f1 + f2 - f12 gives the classical terms - // The difference between this and strong screening of screen5 - // is that we replaced the classical term which is f1 + f2 - f12 - // using results from Chabrier&Potekhin1998. + if constexpr (do_T_derivatives) { + quantum_corr_1_dT = tau12dT / tau12 * quantum_corr_1 - tau12 * + b_fac_dT * (15.0_rt/32.0_rt * amrex::Math::powi<2>(b_fac) - + 0.014_rt * 4.0_rt * amrex::Math::powi<3>(b_fac) - + 0.128_rt * 5.0_rt * amrex::Math::powi<4>(b_fac)); + + quantum_corr_2_dT = Gamma_eff_dT / Gamma_eff * quantum_corr_2 - Gamma_eff * + b_fac_dT * (0.0055_rt * 4.0_rt * amrex::Math::powi<3>(b_fac) - + 0.0098_rt * 5.0_rt * amrex::Math::powi<4>(b_fac) + + 0.0048_rt * 6.0_rt * amrex::Math::powi<5>(b_fac)); + } + } + // See Calder2007 Appendix Eq. A8. + // f1 + f2 - f12 gives the classical terms + // The difference between this and strong screening of screen5 + // is that we replaced the classical term which is f1 + f2 - f12 + // using results from Chabrier&Potekhin1998. - Real h12 = f1 + f2 - f12 + quantum_corr_1 + quantum_corr_2; + Real h12 = f1 + f2 - f12 + quantum_corr_1 + quantum_corr_2; - [[maybe_unused]] Real dh12dT; - if constexpr (do_T_derivatives) { - dh12dT = f1dT + f2dT - f12dT + quantum_corr_1_dT + quantum_corr_2_dT; - } + [[maybe_unused]] Real dh12dT; + if constexpr (do_T_derivatives) { + dh12dT = f1dT + f2dT - f12dT + quantum_corr_1_dT + quantum_corr_2_dT; + } - Real h12_max = 300.0_rt; - h12 = amrex::min(h12_max, h12); + Real h12_max = 300.0_rt; + h12 = amrex::min(h12_max, h12); - scor = std::exp(h12); + scor = std::exp(h12); - if constexpr (do_T_derivatives) { - if (h12 == h12_max) { - scordt = 0.0_rt; - } else { - scordt = scor * dh12dT; + if constexpr (do_T_derivatives) { + if (h12 == h12_max) { + scordt = 0.0_rt; + } else { + scordt = scor * dh12dT; + } } - } - } #endif diff --git a/screening/screen_data.H b/screening/screen_data.H index 46d1322ab6..571d730860 100644 --- a/screening/screen_data.H +++ b/screening/screen_data.H @@ -14,66 +14,66 @@ namespace scrn { public: - amrex::Real z1 = -1; - amrex::Real z2 = -1; - 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 + amrex::Real z1 = -1; + amrex::Real z2 = -1; + 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; + 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; + 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; + 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; + 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, const amrex::Real z2_pass, const amrex::Real a2_pass) const { - // a simple function for unit testing / debug runs to - // ensure that we are accessing the proper screening info + // a simple function for unit testing / debug runs to + // ensure that we are accessing the proper screening info - return (z1_pass == z1) && - (z2_pass == z2) && - (a1_pass == a1) && - (a2_pass == a2); + return (z1_pass == z1) && + (z2_pass == z2) && + (a1_pass == a1) && + (a2_pass == a2); } }; diff --git a/sphinx_docs/source/screening.rst b/sphinx_docs/source/screening.rst index e54cd1a878..1d99a925c7 100644 --- a/sphinx_docs/source/screening.rst +++ b/sphinx_docs/source/screening.rst @@ -48,3 +48,7 @@ The options are: This disables screening by always returning 1 for the screening enhancement factor. + +Runtime Options +---------------- +* ``screening.enable_chabrier1998_quantum_corr = 1`` in the input file enables an additional quantum correction term added to the screening factor when ``SCREEN_METHOD=chabrier1998``. This is disabled by default since ``chabrier1998`` is often used along with ``USE_NSE_NET=TRUE``, and the NSE solver doesn't include quantum corrections.