From 37a921e7b846afa8ce6d783ca5010ecc1e76a862 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 14:08:22 -0400 Subject: [PATCH 01/10] move some Fermi integral coefficients to static arrays no need to set these each time we enter the function --- neutrinos/sneut5.H | 192 +++++++++++++++++++++++---------------------- 1 file changed, 100 insertions(+), 92 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 9348f75f86..20eca2a1aa 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -6,6 +6,90 @@ using namespace amrex; +namespace ifermi12_coeffs { + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a1 = { + 1.999266880833e4_rt, + 5.702479099336e3_rt, + 6.610132843877e2_rt, + 3.818838129486e1_rt, + 1.0e0_rt}; + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b1 = { + 1.771804140488e4_rt, + -2.014785161019e3_rt, + 9.130355392717e1_rt, + -1.670718177489e0_rt}; + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a2 = { + -1.277060388085e-2_rt, + 7.187946804945e-2_rt, + -4.262314235106e-1_rt, + 4.997559426872e-1_rt, + -1.285579118012e0_rt, + -3.930805454272e-1_rt, + 1.0e0_rt}; + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b2 = { + -9.745794806288e-3_rt, + 5.485432756838e-2_rt, + -3.299466243260e-1_rt, + 4.077841975923e-1_rt, + -1.145531476975e0_rt, + -6.067091689181e-2_rt}; +}; + +namespace zfermim12_coeffs { + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a1 = { + 1.71446374704454e7_rt, + 3.88148302324068e7_rt, + 3.16743385304962e7_rt, + 1.14587609192151e7_rt, + 1.83696370756153e6_rt, + 1.14980998186874e5_rt, + 1.98276889924768e3_rt, + 1.0e0_rt}; + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b1 = { + 9.67282587452899e6_rt, + 2.87386436731785e7_rt, + 3.26070130734158e7_rt, + 1.77657027846367e7_rt, + 4.81648022267831e6_rt, + 6.13709569333207e5_rt, + 3.13595854332114e4_rt, + 4.35061725080755e2_rt}; + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a2 = { + -4.46620341924942e-15_rt, + -1.58654991146236e-12_rt, + -4.44467627042232e-10_rt, + -6.84738791621745e-8_rt, + -6.64932238528105e-6_rt, + -3.69976170193942e-4_rt, + -1.12295393687006e-2_rt, + -1.60926102124442e-1_rt, + -8.52408612877447e-1_rt, + -7.45519953763928e-1_rt, + 2.98435207466372e0_rt, + 1.0e0_rt}; + + HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b2 = { + -2.23310170962369e-15_rt, + -7.94193282071464e-13_rt, + -2.22564376956228e-10_rt, + -3.43299431079845e-8_rt, + -3.33919612678907e-6_rt, + -1.86432212187088e-4_rt, + -5.69764436880529e-3_rt, + -8.34904593067194e-2_rt, + -4.78770844009440e-1_rt, + -4.99759250374148e-1_rt, + 1.86795964993052e0_rt, + 4.16485970495288e-1_rt}; +}; + AMREX_GPU_HOST_DEVICE AMREX_INLINE Real ifermi12(const Real f) { @@ -17,10 +101,6 @@ Real ifermi12(const Real f) // declare work variables int m1,k1,m2,k2; Real an,rn,den,ff; - Array1D a1; - Array1D b1; - Array1D a2; - Array1D b2; // the return value Real ifermi12r; @@ -33,32 +113,6 @@ Real ifermi12(const Real f) m2 = 6; k2 = 5; - a1(1) = 1.999266880833e4_rt; - a1(2) = 5.702479099336e3_rt; - a1(3) = 6.610132843877e2_rt; - a1(4) = 3.818838129486e1_rt; - a1(5) = 1.0e0_rt; - - b1(1) = 1.771804140488e4_rt; - b1(2) = -2.014785161019e3_rt; - b1(3) = 9.130355392717e1_rt; - b1(4) = -1.670718177489e0_rt; - - a2(1) = -1.277060388085e-2_rt; - a2(2) = 7.187946804945e-2_rt; - a2(3) = -4.262314235106e-1_rt; - a2(4) = 4.997559426872e-1_rt; - a2(5) = -1.285579118012e0_rt; - a2(6) = -3.930805454272e-1_rt; - a2(7) = 1.0e0_rt; - - b2(1) = -9.745794806288e-3_rt; - b2(2) = 5.485432756838e-2_rt; - b2(3) = -3.299466243260e-1_rt; - b2(4) = 4.077841975923e-1_rt; - b2(5) = -1.145531476975e0_rt; - b2(6) = -6.067091689181e-2_rt; - if (f < 4.0e0_rt) { // build sum_{i=0, m1} a_i x**i. This is the numerator @@ -72,18 +126,18 @@ Real ifermi12(const Real f) // // in the starting rn here, the leading f is actually // a1(m1+1) * f, but that element is 1 - rn = f + a1(m1); + rn = f + ifermi12_coeffs::a1(m1); for (int i = m1 - 1; i >= 1; --i) { - rn = rn*f + a1(i); + rn = rn*f + ifermi12_coeffs::a1(i); } // now we do the denominator in Eq. 4. None of the coefficients // are 1, so we loop over all - den = b1(k1+1); + den = ifermi12_coeffs::b1(k1+1); for (int i = k1; i >= 1; --i) { - den = den*f + b1(i); + den = den*f + ifermi12_coeffs::b1(i); } // Eq. 6 of Antia @@ -97,16 +151,16 @@ Real ifermi12(const Real f) // this construction is the same as above, but using the // second set of coefficients - rn = ff + a2(m2); + rn = ff + ifermi12_coeffs::a2(m2); for (int i = m2 - 1; i >= 1; --i) { - rn = rn*ff + a2(i); + rn = rn*ff + ifermi12_coeffs::a2(i); } - den = b2(k2+1); + den = ifermi12_coeffs::b2(k2+1); for (int i = k2; i >= 1; --i) { - den = den*ff + b2(i); + den = den*ff + ifermi12_coeffs::b2(i); } ifermi12r = rn/(den*ff); @@ -127,8 +181,6 @@ Real zfermim12(const Real x) // declare work variables int m1,k1,m2,k2; Real rn,den,xx; - Array1D a1,b1; - Array1D a2,b2; // return value Real zfermim12r; @@ -140,63 +192,19 @@ Real zfermim12(const Real x) m2 = 11; k2 = 11; - a1(1) = 1.71446374704454e7_rt; - a1(2) = 3.88148302324068e7_rt; - a1(3) = 3.16743385304962e7_rt; - a1(4) = 1.14587609192151e7_rt; - a1(5) = 1.83696370756153e6_rt; - a1(6) = 1.14980998186874e5_rt; - a1(7) = 1.98276889924768e3_rt; - a1(8) = 1.0e0_rt; - - b1(1) = 9.67282587452899e6_rt; - b1(2) = 2.87386436731785e7_rt; - b1(3) = 3.26070130734158e7_rt; - b1(4) = 1.77657027846367e7_rt; - b1(5) = 4.81648022267831e6_rt; - b1(6) = 6.13709569333207e5_rt; - b1(7) = 3.13595854332114e4_rt; - b1(8) = 4.35061725080755e2_rt; - - a2(1) = -4.46620341924942e-15_rt; - a2(2) = -1.58654991146236e-12_rt; - a2(3) = -4.44467627042232e-10_rt; - a2(4) = -6.84738791621745e-8_rt; - a2(5) = -6.64932238528105e-6_rt; - a2(6) = -3.69976170193942e-4_rt; - a2(7) = -1.12295393687006e-2_rt; - a2(8) = -1.60926102124442e-1_rt; - a2(9) = -8.52408612877447e-1_rt; - a2(10) = -7.45519953763928e-1_rt; - a2(11) = 2.98435207466372e0_rt; - a2(12) = 1.0e0_rt; - - b2(1) = -2.23310170962369e-15_rt; - b2(2) = -7.94193282071464e-13_rt; - b2(3) = -2.22564376956228e-10_rt; - b2(4) = -3.43299431079845e-8_rt; - b2(5) = -3.33919612678907e-6_rt; - b2(6) = -1.86432212187088e-4_rt; - b2(7) = -5.69764436880529e-3_rt; - b2(8) = -8.34904593067194e-2_rt; - b2(9) = -4.78770844009440e-1_rt; - b2(10) = -4.99759250374148e-1_rt; - b2(11) = 1.86795964993052e0_rt; - b2(12) = 4.16485970495288e-1_rt; - if (x < 2.0e0_rt) { xx = std::exp(x); - rn = xx + a1(m1); + rn = xx + zfermim12_coeffs::a1(m1); for (int i = m1 - 1; i >= 1; --i) { - rn = rn*xx + a1(i); + rn = rn*xx + zfermim12_coeffs::a1(i); } - den = b1(k1+1); + den = zfermim12_coeffs::b1(k1+1); for (int i = k1; i >= 1; --i) { - den = den*xx + b1(i); + den = den*xx + zfermim12_coeffs::b1(i); } zfermim12r = xx * rn/den; @@ -204,16 +212,16 @@ Real zfermim12(const Real x) } else { xx = 1.0e0_rt/(x*x); - rn = xx + a2(m2); + rn = xx + zfermim12_coeffs::a2(m2); for (int i = m2 - 1; i >= 1; --i) { - rn = rn*xx + a2(i); + rn = rn*xx + zfermim12_coeffs::a2(i); } - den = b2(k2+1); + den = zfermim12_coeffs::b2(k2+1); for (int i = k2; i >= 1; --i) { - den = den*xx + b2(i); + den = den*xx + zfermim12_coeffs::b2(i); } zfermim12r = std::sqrt(x)*rn/den; From 97523064aebb0b8533e9e86bb0dd606c55dafb50 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 14:14:22 -0400 Subject: [PATCH 02/10] some more cleaning --- neutrinos/sneut5.H | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 20eca2a1aa..3225cc6822 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -99,19 +99,18 @@ Real ifermi12(const Real f) // maximum error is 4.19e-9_rt. reference: antia apjs 84,101 1993 // declare work variables - int m1,k1,m2,k2; - Real an,rn,den,ff; + Real rn,den,ff; // the return value Real ifermi12r; // load the coefficients of the expansion from Table 8 of Antia - an = 0.5e0_rt; - m1 = 4; - k1 = 3; - m2 = 6; - k2 = 5; + constexpr Real an{0.5e0_rt}; + constexpr int m1{4}; + constexpr int k1{3}; + constexpr int m2{6}; + constexpr int k2{5}; if (f < 4.0e0_rt) { @@ -179,7 +178,6 @@ Real zfermim12(const Real x) // reference: antia apjs 84,101 1993 // declare work variables - int m1,k1,m2,k2; Real rn,den,xx; // return value @@ -187,10 +185,10 @@ Real zfermim12(const Real x) // load the coefficients of the expansion from Table 2 of Antia - m1 = 7; - k1 = 7; - m2 = 11; - k2 = 11; + constexpr int m1{7}; + constexpr int k1{7}; + constexpr int m2{11}; + constexpr int k2{11}; if (x < 2.0e0_rt) { From 1d94c87dfc3620cb4aa63498c84c45943195c29e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 17:48:36 -0400 Subject: [PATCH 03/10] move common sneut factors into a struct --- neutrinos/sneut5.H | 361 +++++++++++++++++++++++++-------------------- 1 file changed, 203 insertions(+), 158 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 3225cc6822..0aa51d41ee 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -229,6 +229,84 @@ Real zfermim12(const Real x) return zfermim12r; } +struct sneutf_t { + + Real deni; + Real tempi; + Real abari; + Real zbari; + + Real ye; + + Real t9; + Real xl; + Real xldt; + Real xlp5; + Real xl2; + Real xl3; + Real xl4; + Real xl5; + Real xl6; + Real xl7; + Real xl8; + Real xl9; + Real xlmp5; + Real xlm1; + Real xlm2; + Real xlm3; + Real xlm4; + Real rm; + Real rmda; + Real rmdz; + Real rmi; + +}; + + +AMREX_GPU_HOST_DEVICE inline +sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { + + constexpr Real con1 = 1.0e0_rt/5.9302e0_rt; + + sneutf_t sf; + + // to avoid lots of divisions + sf.deni = 1.0e0_rt / den; + sf.tempi = 1.0e0_rt / temp; + sf.abari = 1.0e0_rt / abar; + sf.zbari = 1.0e0_rt / zbar; + + // some composition variables + sf.ye = zbar * sf.abari; + + // some frequent factors + sf.t9 = temp * 1.0e-9_rt; + sf.xl = sf.t9 * con1; + sf.xldt = 1.0e-9_rt * con1; + sf.xlp5 = std::sqrt(sf.xl); + sf.xl2 = sf.xl * sf.xl; + sf.xl3 = sf.xl2 * sf.xl; + sf.xl4 = sf.xl3 * sf.xl; + sf.xl5 = sf.xl4 * sf.xl; + sf.xl6 = sf.xl5 * sf.xl; + sf.xl7 = sf.xl6 * sf.xl; + sf.xl8 = sf.xl7 * sf.xl; + sf.xl9 = sf.xl8 * sf.xl; + sf.xlmp5 = 1.0e0_rt / sf.xlp5; + sf.xlm1 = 1.0e0_rt / sf.xl; + sf.xlm2 = sf.xlm1 * sf.xlm1; + sf.xlm3 = sf.xlm1 * sf.xlm2; + sf.xlm4 = sf.xlm1 * sf.xlm3; + + sf.rm = den * sf.ye; + sf.rmda = -sf.rm * sf.abari; + sf.rmdz = den * sf.abari; + sf.rmi = 1.0e0_rt / sf.rm; + + return sf; + +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -255,22 +333,21 @@ void sneut5(const Real temp, const Real den, dsnudz = derivative of snu with zbar */ - Real t9,xl,xldt,xlp5,xl2,xl3,xl4,xl5,xl6,xl7,xl8,xl9, - xlmp5,xlm1,xlm2,xlm3,xlm4,xlnt,cc,den6,tfermi, + Real xlnt,cc,den6,tfermi, a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06, c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d{0.0},f0, - f1{0.0},deni,tempi,abari,zbari,f2,f3,z,ye, + f1{0.0},f2,f3,z,ye, dum,dumdt,gum,gumdt,dumda,dumdz,gumda,gumdz; // pair production - Real rm,rmi,gl,gldt, + Real gl,gldt, zeta,zetadt,zeta2,zeta3, xnum,xnumdt,xden,xdendt, fpair,fpairdt,qpair,qpairdt, fpairda,fpairdz,qpairda,qpairdz, - rmda,rmdz,zetada,zetadz, + zetada,zetadz, xnumda,xnumdz,xdenda,xdendz; // plasma @@ -307,7 +384,6 @@ void sneut5(const Real temp, const Real den, constexpr Real fac3 = M_PI / 5.0e0_rt; constexpr Real oneth = 1.0e0_rt/3.0e0_rt; constexpr Real twoth = 2.0e0_rt/3.0e0_rt; - constexpr Real con1 = 1.0e0_rt/5.9302e0_rt; constexpr Real sixth = 1.0e0_rt/6.0e0_rt; constexpr Real iln10 = 4.342944819032518e-1_rt; @@ -363,48 +439,17 @@ void sneut5(const Real temp, const Real den, if (temp < 1.0e7_rt) return; - // to avoid lots of divisions - deni = 1.0e0_rt/den; - tempi = 1.0e0_rt/temp; - abari = 1.0e0_rt/abar; - zbari = 1.0e0_rt/zbar; + auto sf = get_sneut_factors(den, temp, abar, zbar); - // some composition variables - ye = zbar*abari; - - // some frequent factors - t9 = temp * 1.0e-9_rt; - xl = t9 * con1; - xldt = 1.0e-9_rt * con1; - xlp5 = std::sqrt(xl); - xl2 = xl*xl; - xl3 = xl2*xl; - xl4 = xl3*xl; - xl5 = xl4*xl; - xl6 = xl5*xl; - xl7 = xl6*xl; - xl8 = xl7*xl; - xl9 = xl8*xl; - xlmp5 = 1.0e0_rt/xlp5; - xlm1 = 1.0e0_rt/xl; - xlm2 = xlm1*xlm1; - xlm3 = xlm1*xlm2; - xlm4 = xlm1*xlm3; - - rm = den*ye; - rmda = -rm*abari; - rmdz = den*abari; - rmi = 1.0e0_rt/rm; - - a0 = rm * 1.0e-9_rt; + a0 = sf.rm * 1.0e-9_rt; a1 = std::pow(a0, oneth); - zeta = a1 * xlm1; + zeta = a1 * sf.xlm1; if constexpr (do_derivatives) { - a2 = oneth * a1*rmi * xlm1; - zetadt = -a1 * xlm2 * xldt; - zetada = a2 * rmda; - zetadz = a2 * rmdz; + a2 = oneth * a1*sf.rmi * sf.xlm1; + zetadt = -a1 * sf.xlm2 * sf.xldt; + zetada = a2 * sf.rmda; + zetadz = a2 * sf.rmdz; } zeta2 = zeta * zeta; @@ -414,16 +459,16 @@ void sneut5(const Real temp, const Real den, // for reactions like e+ + e- => nu_e + nubar_e // equation 2.8 - gl = 1.0e0_rt - 13.04e0_rt*xl2 +133.5e0_rt*xl4 +1534.0e0_rt*xl6 +918.6e0_rt*xl8; + gl = 1.0e0_rt - 13.04e0_rt*sf.xl2 +133.5e0_rt*sf.xl4 +1534.0e0_rt*sf.xl6 +918.6e0_rt*sf.xl8; if constexpr (do_derivatives) { - gldt = xldt*(-26.08e0_rt*xl +534.0e0_rt*xl3 +9204.0e0_rt*xl5 +7348.8e0_rt*xl7); + gldt = sf.xldt*(-26.08e0_rt*sf.xl +534.0e0_rt*sf.xl3 +9204.0e0_rt*sf.xl5 +7348.8e0_rt*sf.xl7); } // equation 2.7 a1 = 6.002e19_rt + 2.084e20_rt*zeta + 1.872e21_rt*zeta2; - if (t9 < 10.0_rt) { + if (sf.t9 < 10.0_rt) { b1 = std::exp(-5.5924e0_rt*zeta); if constexpr (do_derivatives) { b2 = -b1*5.5924e0_rt; @@ -444,22 +489,22 @@ void sneut5(const Real temp, const Real den, xnumdz = c*zetadz; } - if (t9 < 10.0_rt) { - a1 = 9.383e-1_rt*xlm1 - 4.141e-1_rt*xlm2 + 5.829e-2_rt*xlm3; + if (sf.t9 < 10.0_rt) { + a1 = 9.383e-1_rt*sf.xlm1 - 4.141e-1_rt*sf.xlm2 + 5.829e-2_rt*sf.xlm3; if constexpr (do_derivatives) { - a2 = -9.383e-1_rt*xlm2 + 2.0e0_rt*4.141e-1_rt*xlm3 - 3.0e0_rt*5.829e-2_rt*xlm4; + a2 = -9.383e-1_rt*sf.xlm2 + 2.0e0_rt*4.141e-1_rt*sf.xlm3 - 3.0e0_rt*5.829e-2_rt*sf.xlm4; } } else { - a1 = 1.2383e0_rt*xlm1 - 8.141e-1_rt*xlm2; + a1 = 1.2383e0_rt*sf.xlm1 - 8.141e-1_rt*sf.xlm2; if constexpr (do_derivatives) { - a2 = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3; + a2 = -1.2383e0_rt*sf.xlm2 + 2.0e0_rt*8.141e-1_rt*sf.xlm3; } } xden = zeta3 + a1; if constexpr (do_derivatives) { b1 = 3.0e0_rt*zeta2; - xdendt = b1*zetadt + a2*xldt; + xdendt = b1*zetadt + a2*sf.xldt; xdenda = b1*zetada; xdendz = b1*zetadz; } @@ -473,24 +518,24 @@ void sneut5(const Real temp, const Real den, } // equation 2.6 - a1 = 10.7480e0_rt*xl2 + 0.3967e0_rt*xlp5 + 1.005e0_rt; + a1 = 10.7480e0_rt*sf.xl2 + 0.3967e0_rt*sf.xlp5 + 1.005e0_rt; xnum = 1.0e0_rt/a1; if constexpr (do_derivatives) { - a2 = xldt*(2.0e0_rt*10.7480e0_rt*xl + 0.5e0_rt*0.3967e0_rt*xlmp5); + a2 = sf.xldt*(2.0e0_rt*10.7480e0_rt*sf.xl + 0.5e0_rt*0.3967e0_rt*sf.xlmp5); xnumdt = -xnum*xnum*a2; } - a1 = 7.692e7_rt*xl3 + 9.715e6_rt*xlp5; + a1 = 7.692e7_rt*sf.xl3 + 9.715e6_rt*sf.xlp5; c = 1.0e0_rt/a1; - b1 = 1.0e0_rt + rm*c; + b1 = 1.0e0_rt + sf.rm*c; xden = std::pow(b1, -0.3e0_rt); if constexpr (do_derivatives) { - a2 = xldt*(3.0e0_rt*7.692e7_rt*xl2 + 0.5e0_rt*9.715e6_rt*xlmp5); + a2 = sf.xldt*(3.0e0_rt*7.692e7_rt*sf.xl2 + 0.5e0_rt*9.715e6_rt*sf.xlmp5); d = -0.3e0_rt*xden/b1; - xdendt = -d*rm*c*c*a2; - xdenda = d*rmda*c; - xdendz = d*rmdz*c; + xdendt = -d*sf.rm*c*c*a2; + xdenda = d*sf.rmda*c; + xdendz = d*sf.rmdz*c; } qpair = xnum*xden; @@ -501,11 +546,11 @@ void sneut5(const Real temp, const Real den, } // equation 2.5 - a1 = std::exp(-2.0e0_rt*xlm1); + a1 = std::exp(-2.0e0_rt*sf.xlm1); spair = a1*fpair; if constexpr (do_derivatives) { - a2 = a1*2.0e0_rt*xlm2*xldt; + a2 = a1*2.0e0_rt*sf.xlm2*sf.xldt; spairdt = a2*fpair + a1*fpairdt; spairda = a1*fpairda; spairdz = a1*fpairdz; @@ -534,22 +579,22 @@ void sneut5(const Real temp, const Real den, // for collective reactions like gamma_plasmon => nu_e + nubar_e // equation 4.6 - a1 = 1.019e-6_rt*rm; + a1 = 1.019e-6_rt*sf.rm; a2 = std::pow(a1, twoth); b1 = std::sqrt(1.0e0_rt + a2); c00 = 1.0e0_rt/(temp*temp*b1); - gl2 = 1.1095e11_rt * rm * c00; + gl2 = 1.1095e11_rt * sf.rm * c00; if constexpr (do_derivatives) { a3 = twoth*a2/a1; b2 = 1.0e0_rt/b1; - gl2dt = -2.0e0_rt*gl2*tempi; - d = rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt; - gl2da = 1.1095e11_rt * (rmda*c00 - d*rmda); - gl2dz = 1.1095e11_rt * (rmdz*c00 - d*rmdz); + gl2dt = -2.0e0_rt*gl2*sf.tempi; + d = sf.rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt; + gl2da = 1.1095e11_rt * (sf.rmda*c00 - d*sf.rmda); + gl2dz = 1.1095e11_rt * (sf.rmdz*c00 - d*sf.rmdz); } gl = std::sqrt(gl2); @@ -586,20 +631,20 @@ void sneut5(const Real temp, const Real den, } // equation 4.9 and 4.10 - cc = std::log10(2.0e0_rt*rm); + cc = std::log10(2.0e0_rt*sf.rm); xlnt = std::log10(temp); xnum = sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt); xden = sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt); if constexpr (do_derivatives) { - xnumdt = -iln10*0.5e0_rt*tempi; - a2 = iln10*sixth*rmi; - xnumda = a2*rmda; - xnumdz = a2*rmdz; - - xdendt = iln10*0.5e0_rt*tempi; - xdenda = a2*rmda; - xdendz = a2*rmdz; + xnumdt = -iln10*0.5e0_rt*sf.tempi; + a2 = iln10*sixth*sf.rmi; + xnumda = a2*sf.rmda; + xnumdz = a2*sf.rmdz; + + xdendt = iln10*0.5e0_rt*sf.tempi; + xdenda = a2*sf.rmda; + xdendz = a2*sf.rmdz; } // equation 4.11 @@ -679,12 +724,12 @@ void sneut5(const Real temp, const Real den, splasdz = a2*splasdz + a3*gl2dz*a1; } - a2 = 0.93153e0_rt * 3.0e21_rt * xl9; + a2 = 0.93153e0_rt * 3.0e21_rt * sf.xl9; a1 = splas; splas = a2*a1; if constexpr (do_derivatives) { - a3 = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt; + a3 = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*sf.xl8*sf.xldt; splasdt = a2*splasdt + a3*a1; splasda = a2*splasda; splasdz = a2*splasdz; @@ -783,7 +828,7 @@ void sneut5(const Real temp, const Real den, // T > 1.e9 - tau = std::log10(t9); + tau = std::log10(sf.t9); cc = 1.5654e0_rt; c00 = 9.581e10_rt; c01 = 4.107e8_rt; @@ -824,7 +869,7 @@ void sneut5(const Real temp, const Real den, } - taudt = iln10*tempi; + taudt = iln10*sf.tempi; // equation 3.7, compute the expensive trig functions only one time cos1 = std::cos(fac1*tau); @@ -892,12 +937,12 @@ void sneut5(const Real temp, const Real den, xnumdz = dumdz*z - dum*z*cc*zetadz; } - xden = zeta3 + 6.290e-3_rt*xlm1 + 7.483e-3_rt*xlm2 + 3.061e-4_rt*xlm3; + xden = zeta3 + 6.290e-3_rt*sf.xlm1 + 7.483e-3_rt*sf.xlm2 + 3.061e-4_rt*sf.xlm3; dum = 3.0e0_rt*zeta2; if constexpr (do_derivatives) { - xdendt = dum*zetadt - xldt*(6.290e-3_rt*xlm2 - + 2.0e0_rt*7.483e-3_rt*xlm3 + 3.0e0_rt*3.061e-4_rt*xlm4); + xdendt = dum*zetadt - sf.xldt*(6.290e-3_rt*sf.xlm2 + + 2.0e0_rt*7.483e-3_rt*sf.xlm3 + 3.0e0_rt*3.061e-4_rt*sf.xlm4); xdenda = dum*zetada; xdendz = dum*zetadz; } @@ -910,24 +955,24 @@ void sneut5(const Real temp, const Real den, } // equation 3.3 - a0 = 1.0e0_rt + 2.045e0_rt * xl; + a0 = 1.0e0_rt + 2.045e0_rt * sf.xl; xnum = 0.666e0_rt*std::pow(a0, -2.066e0_rt); if constexpr (do_derivatives) { - xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*xldt; + xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*sf.xldt; } - dum = 1.875e8_rt*xl + 1.653e8_rt*xl2 + 8.499e8_rt*xl3 - 1.604e8_rt*xl4; + dum = 1.875e8_rt*sf.xl + 1.653e8_rt*sf.xl2 + 8.499e8_rt*sf.xl3 - 1.604e8_rt*sf.xl4; if constexpr (do_derivatives) { - dumdt = xldt*(1.875e8_rt + 2.0e0_rt*1.653e8_rt*xl + 3.0e0_rt*8.499e8_rt*xl2 - - 4.0e0_rt*1.604e8_rt*xl3); + dumdt = sf.xldt*(1.875e8_rt + 2.0e0_rt*1.653e8_rt*sf.xl + 3.0e0_rt*8.499e8_rt*sf.xl2 + - 4.0e0_rt*1.604e8_rt*sf.xl3); } z = 1.0e0_rt/dum; - xden = 1.0e0_rt + rm*z; + xden = 1.0e0_rt + sf.rm*z; if constexpr (do_derivatives) { - xdendt = -rm*z*z*dumdt; - xdenda = rmda*z; - xdendz = rmdz*z; + xdendt = -sf.rm*z*z*dumdt; + xdenda = sf.rmda*z; + xdendz = sf.rmdz*z; } z = 1.0e0_rt/xden; @@ -942,19 +987,19 @@ void sneut5(const Real temp, const Real den, } // equation 3.2 - sphot = xl5 * fphot; + sphot = sf.xl5 * fphot; if constexpr (do_derivatives) { - sphotdt = 5.0e0_rt*xl4*xldt*fphot + xl5*fphotdt; - sphotda = xl5*fphotda; - sphotdz = xl5*fphotdz; + sphotdt = 5.0e0_rt*sf.xl4*sf.xldt*fphot + sf.xl5*fphotdt; + sphotda = sf.xl5*fphotda; + sphotdz = sf.xl5*fphotdz; } a1 = sphot; - sphot = rm*a1; + sphot = sf.rm*a1; if constexpr (do_derivatives) { - sphotdt = rm*sphotdt; - sphotda = rm*sphotda + rmda*a1; - sphotdz = rm*sphotdz + rmdz*a1; + sphotdt = sf.rm*sphotdt; + sphotda = sf.rm*sphotda + sf.rmda*a1; + sphotdz = sf.rm*sphotdz + sf.rmdz*a1; } a1 = tfac4*(1.0e0_rt - tfac3 * qphot); @@ -998,7 +1043,7 @@ void sneut5(const Real temp, const Real den, t8m6 = t8m5*t8m1; - tfermi = 5.9302e9_rt*(std::sqrt(1.0e0_rt+1.018e0_rt*std::pow(den6*ye, twoth))-1.0e0_rt); + tfermi = 5.9302e9_rt*(std::sqrt(1.0e0_rt+1.018e0_rt*std::pow(den6*sf.ye, twoth))-1.0e0_rt); // "weak" degenerate electrons only if (temp > 0.3e0_rt * tfermi) { @@ -1010,11 +1055,11 @@ void sneut5(const Real temp, const Real den, } z = 1.0e0_rt/dum; - eta = rm*z; + eta = sf.rm*z; if constexpr (do_derivatives) { - etadt = -rm*z*z*dumdt; - etada = rmda*z; - etadz = rmdz*z; + etadt = -sf.rm*z*z*dumdt; + etada = sf.rmda*z; + etadz = sf.rmdz*z; } etam1 = 1.0e0_rt/eta; @@ -1061,13 +1106,13 @@ void sneut5(const Real temp, const Real den, a0 = 230.0e0_rt + 6.7e5_rt*t8m2 + 7.66e9_rt*t8m5; f0 = (-2.0e0_rt*6.7e5_rt*t8m3 - 5.0e0_rt*7.66e9_rt*t8m6)*1.0e-8_rt; - z = 1.0e0_rt + rm*1.0e-9_rt; + z = 1.0e0_rt + sf.rm*1.0e-9_rt; dum = a0*z; if constexpr (do_derivatives) { dumdt = f0*z; z = a0*1.0e-9_rt; - dumda = z*rmda; - dumdz = z*rmdz; + dumda = z*sf.rmda; + dumdz = z*sf.rmdz; } xnum = 1.0e0_rt/dum; @@ -1090,12 +1135,12 @@ void sneut5(const Real temp, const Real den, } z = std::pow(den, 0.656e0_rt); - dum = c00*rmi + c01 + c02*z; + dum = c00*sf.rmi + c01 + c02*z; if constexpr (do_derivatives) { - dumdt = dd00*rmi + dd01 + dd02*z; - z = -c00*rmi*rmi; - dumda = z*rmda; - dumdz = z*rmdz; + dumdt = dd00*sf.rmi + dd01 + dd02*z; + z = -c00*sf.rmi*sf.rmi; + dumda = z*sf.rmda; + dumdz = z*sf.rmdz; } xden = 1.0e0_rt/dum; @@ -1114,11 +1159,11 @@ void sneut5(const Real temp, const Real den, } // equation 5.1 - dum = 0.5738e0_rt*zbar*ye*t86*den; + dum = 0.5738e0_rt*zbar*sf.ye*t86*den; if constexpr (do_derivatives) { - dumdt = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt; - dumda = -dum*abari; - dumdz = 0.5738e0_rt*2.0e0_rt*ye*t86*den; + dumdt = 0.5738e0_rt*zbar*sf.ye*6.0e0_rt*t85*den*1.0e-8_rt; + dumda = -dum*sf.abari; + dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; } z = tfac4*fbrem - tfac5*gbrem; @@ -1135,7 +1180,7 @@ void sneut5(const Real temp, const Real den, } else { u = fac3 * (std::log10(den) - 3.0e0_rt); - //a0 = iln10*fac3*deni; + //a0 = iln10*fac3*sf.deni; // compute the expensive trig functions of equation 5.21 only once cos1 = std::cos(u); @@ -1210,11 +1255,11 @@ void sneut5(const Real temp, const Real den, // - 0.00031e0_rt*sin4*4.0e0_rt - 0.00018e0_rt*cos4*4.0e0_rt // - 0.00069e0_rt*sin5*5.0e0_rt); - dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*abari, oneth); + dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*sf.abari, oneth); if constexpr (do_derivatives) { - dumdt = -dum*tempi; - dumda = -oneth*dum*abari; - dumdz = 2.0e0_rt*dum*zbari; + dumdt = -dum*sf.tempi; + dumda = -oneth*dum*sf.abari; + dumdz = 2.0e0_rt*dum*sf.zbari; } gm1 = 1.0e0_rt/dum; @@ -1247,11 +1292,11 @@ void sneut5(const Real temp, const Real den, } // equation 5.17 - dum = 0.5738e0_rt*zbar*ye*t86*den; + dum = 0.5738e0_rt*zbar*sf.ye*t86*den; if constexpr (do_derivatives) { - dumdt = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt; - dumda = -dum*abari; - dumdz = 0.5738e0_rt*2.0e0_rt*ye*t86*den; + dumdt = 0.5738e0_rt*zbar*sf.ye*6.0e0_rt*t85*den*1.0e-8_rt; + dumda = -dum*sf.abari; + dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; } z = tfac4*fliq - tfac5*gliq; @@ -1266,11 +1311,11 @@ void sneut5(const Real temp, const Real den, // recombination neutrino section // for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e // equation 6.11 solved for nu - xnum = 1.10520e8_rt * den * ye /(temp*std::sqrt(temp)); + xnum = 1.10520e8_rt * den * sf.ye /(temp*std::sqrt(temp)); if constexpr (do_derivatives) { - xnumdt = -1.50e0_rt*xnum*tempi; - xnumda = -xnum*abari; - xnumdz = xnum*zbari; + xnumdt = -1.50e0_rt*xnum*sf.tempi; + xnumda = -xnum*sf.abari; + xnumdz = xnum*sf.zbari; } // the chemical potential @@ -1312,11 +1357,11 @@ void sneut5(const Real temp, const Real den, // equation 6.7, 6.13 and 6.14 if (nu >= -20.0_rt && nu <= 10.0_rt) { - zeta = 1.579e5_rt*zbar*zbar*tempi; + zeta = 1.579e5_rt*zbar*zbar*sf.tempi; if constexpr (do_derivatives) { - zetadt = -zeta*tempi; + zetadt = -zeta*sf.tempi; zetada = 0.0e0_rt; - zetadz = 2.0e0_rt*zeta*zbari; + zetadz = 2.0e0_rt*zeta*sf.zbari; } c00 = 1.0e0_rt/(1.0e0_rt + f1*nu + f2*nu2 + f3*nu3); @@ -1359,50 +1404,50 @@ void sneut5(const Real temp, const Real den, a1 = 1.0e0_rt/dum; a2 = 1.0e0_rt/bigj; - sreco = tfac6 * 2.649e-18_rt * ye * std::pow(zbar, 13.0_rt) * den * bigj*a1; + sreco = tfac6 * 2.649e-18_rt * sf.ye * std::pow(zbar, 13.0_rt) * den * bigj*a1; if constexpr (do_derivatives) { srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1); - srecoda = sreco*(-1.0e0_rt*abari + bigjda*a2 - z*(zetada+nuda)*a1); - srecodz = sreco*(14.0e0_rt*zbari + bigjdz*a2 - z*(zetadz+nudz)*a1); + srecoda = sreco*(-1.0e0_rt*sf.abari + bigjda*a2 - z*(zetada+nuda)*a1); + srecodz = sreco*(14.0e0_rt*sf.zbari + bigjdz*a2 - z*(zetadz+nudz)*a1); } } // convert from erg/cm^3/s to erg/g/s // comment these out to duplicate the itoh et al plots - spair = spair*deni; + spair = spair*sf.deni; if constexpr (do_derivatives) { - spairdt = spairdt*deni; - spairda = spairda*deni; - spairdz = spairdz*deni; + spairdt = spairdt*sf.deni; + spairda = spairda*sf.deni; + spairdz = spairdz*sf.deni; } - splas = splas*deni; + splas = splas*sf.deni; if constexpr (do_derivatives) { - splasdt = splasdt*deni; - splasda = splasda*deni; - splasdz = splasdz*deni; + splasdt = splasdt*sf.deni; + splasda = splasda*sf.deni; + splasdz = splasdz*sf.deni; } - sphot = sphot*deni; + sphot = sphot*sf.deni; if constexpr (do_derivatives) { - sphotdt = sphotdt*deni; - sphotda = sphotda*deni; - sphotdz = sphotdz*deni; + sphotdt = sphotdt*sf.deni; + sphotda = sphotda*sf.deni; + sphotdz = sphotdz*sf.deni; } - sbrem = sbrem*deni; + sbrem = sbrem*sf.deni; if constexpr (do_derivatives) { - sbremdt = sbremdt*deni; - sbremda = sbremda*deni; - sbremdz = sbremdz*deni; + sbremdt = sbremdt*sf.deni; + sbremda = sbremda*sf.deni; + sbremdz = sbremdz*sf.deni; } - sreco = sreco*deni; + sreco = sreco*sf.deni; if constexpr (do_derivatives) { - srecodt = srecodt*deni; - srecoda = srecoda*deni; - srecodz = srecodz*deni; + srecodt = srecodt*sf.deni; + srecoda = srecoda*sf.deni; + srecodz = srecodz*sf.deni; } // the total neutrino loss rate From 5902af642b746c6c2dde15623d39073797cd0341 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 17:51:05 -0400 Subject: [PATCH 04/10] fix --- neutrinos/sneut5.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 0aa51d41ee..07f8ca7c5d 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -338,7 +338,7 @@ void sneut5(const Real temp, const Real den, c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d{0.0},f0, - f1{0.0},f2,f3,z,ye, + f1{0.0},f2,f3,z, dum,dumdt,gum,gumdt,dumda,dumdz,gumda,gumdz; // pair production From 0a15fc3a2639bae1b225b96fb1d6f4e66c6d12ea Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 18:26:32 -0400 Subject: [PATCH 05/10] move recombination into its own function --- neutrinos/sneut5.H | 398 ++++++++++++++++++++++++--------------------- 1 file changed, 216 insertions(+), 182 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 07f8ca7c5d..726f1269c6 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -6,6 +6,39 @@ using namespace amrex; +namespace nu_constants { + + // numerical constants + constexpr Real fac1 = 5.0e0_rt * M_PI / 3.0e0_rt; + constexpr Real fac2 = 10.0e0_rt * M_PI; + constexpr Real fac3 = M_PI / 5.0e0_rt; + constexpr Real oneth = 1.0e0_rt/3.0e0_rt; + constexpr Real twoth = 2.0e0_rt/3.0e0_rt; + constexpr Real sixth = 1.0e0_rt/6.0e0_rt; + constexpr Real iln10 = 4.342944819032518e-1_rt; + + // theta is sin**2(theta_weinberg) = 0.2319 plus/minus 0.00005 (1996) + // xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998) + // change theta and xnufam if need be, and the changes will automatically + // propagate through the routine. cv and ca are the vector and axial currents. + + constexpr Real theta = 0.2319e0_rt; + constexpr Real xnufam = 3.0e0_rt; + constexpr Real cv = 0.5e0_rt + 2.0e0_rt * theta; + constexpr Real cvp = 1.0e0_rt - cv; + constexpr Real ca = 0.5e0_rt; + constexpr Real cap = 1.0e0_rt - ca; + constexpr Real tfac1 = cv*cv + ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp+cap*cap); + constexpr Real tfac2 = cv*cv - ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp - cap*cap); + constexpr Real tfac3 = tfac2/tfac1; + constexpr Real tfac4 = 0.5e0_rt * tfac1; + constexpr Real tfac5 = 0.5e0_rt * tfac2; + constexpr Real tfac6 = cv*cv + 1.5e0_rt*ca*ca + (xnufam - 1.0e0_rt)*(cvp*cvp + 1.5e0_rt*cap*cap); + + +} + + namespace ifermi12_coeffs { HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a1 = { @@ -231,6 +264,11 @@ Real zfermim12(const Real x) struct sneutf_t { + Real den; + Real temp; + Real abar; + Real zbar; + Real deni; Real tempi; Real abari; @@ -270,6 +308,11 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { sneutf_t sf; + sf.den = den; + sf.temp = temp; + sf.abar = abar; + sf.zbar = zbar; + // to avoid lots of divisions sf.deni = 1.0e0_rt / den; sf.tempi = 1.0e0_rt / temp; @@ -308,6 +351,130 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_recomb(const sneutf_t& sf, + Real& sreco, Real& srecodt, Real& srecoda, Real& srecodz) { + + + // recombination neutrino section + // for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e + // equation 6.11 solved for nu + + Real xnum = 1.10520e8_rt * sf.den * sf.ye / (sf.temp * std::sqrt(sf.temp)); + Real xnumdt{0.0}; + Real xnumda{0.0}; + Real xnumdz{0.0}; + + if constexpr (do_derivatives) { + xnumdt = -1.50e0_rt*xnum*sf.tempi; + xnumda = -xnum*sf.abari; + xnumdz = xnum*sf.zbari; + } + + // the chemical potential + Real nu = ifermi12(xnum); + + // a0 is d(nu)/d(xnum) + Real a0 = 1.0e0_rt/(0.5e0_rt*zfermim12(nu)); + Real nudt, nuda, nudz; + if constexpr (do_derivatives) { + nudt = a0*xnumdt; + nuda = a0*xnumda; + nudz = a0*xnumdz; + } + Real nu2 = nu * nu; + Real nu3 = nu2 * nu; + + Real a1, a2, a3, b, c, d, f1, f2, f3; + + // table 12 + if (nu >= -20.0_rt && nu < 0.0_rt) { + a1 = 1.51e-2_rt; + a2 = 2.42e-1_rt; + a3 = 1.21e0_rt; + b = 3.71e-2_rt; + c = 9.06e-1_rt; + d = 9.28e-1_rt; + f1 = 0.0e0_rt; + f2 = 0.0e0_rt; + f3 = 0.0e0_rt; + } else if (nu >= 0.0_rt && nu <= 10.0_rt) { + a1 = 1.23e-2_rt; + a2 = 2.66e-1_rt; + a3 = 1.30e0_rt; + b = 1.17e-1_rt; + c = 8.97e-1_rt; + d = 1.77e-1_rt; + f1 = -1.20e-2_rt; + f2 = 2.29e-2_rt; + f3 = -1.04e-3_rt; + } + + // equation 6.7, 6.13 and 6.14 + if (nu >= -20.0_rt && nu <= 10.0_rt) { + + Real zeta = 1.579e5_rt * sf.zbar * sf.zbar * sf.tempi; + Real zetadt, zetada, zetadz; + if constexpr (do_derivatives) { + zetadt = -zeta * sf.tempi; + zetada = 0.0e0_rt; + zetadz = 2.0e0_rt * zeta * sf.zbari; + } + + Real c00 = 1.0e0_rt / (1.0e0_rt + f1 * nu + f2 * nu2 + f3 * nu3); + Real c01 = f1 + f2 * 2.0e0_rt * nu + f3 * 3.0e0_rt * nu2; + Real dum = zeta * c00; + Real dumdt, dumda, dumdz; + if constexpr (do_derivatives) { + dumdt = zetadt * c00 + zeta * c01 * nudt; + dumda = zeta * c01 * nuda; + dumdz = zetadz * c00 + zeta * c01 * nudz; + } + + Real z = 1.0e0_rt / dum; + Real dd00 = std::pow(dum, -2.25_rt); + Real dd01 = std::pow(dum, -4.55_rt); + c00 = a1 * z + a2 * dd00 + a3 * dd01; + c01 = -(a1 * z + 2.25_rt * a2 * dd00 + 4.55_rt * a3 * dd01) * z; + + z = std::exp(c * nu); + dd00 = b * z * (1.0e0_rt + d * dum); + Real gum = 1.0e0_rt + dd00; + Real gumdt, gumda, gumdz; + if constexpr (do_derivatives) { + gumdt = dd00 * c * nudt + b * z * d * dumdt; + gumda = dd00 * c * nuda + b * z * d * dumda; + gumdz = dd00 * c * nudz + b * z * d * dumdz; + } + + z = std::exp(nu); + a1 = 1.0e0_rt / gum; + + Real bigj = c00 * z * a1; + Real bigjdt, bigjda, bigjdz; + if constexpr (do_derivatives) { + bigjdt = c01 * dumdt * z * a1 + c00 * z * nudt * a1 - c00 * z * a1 * a1 * gumdt; + bigjda = c01 * dumda * z * a1 + c00 * z * nuda * a1 - c00 * z * a1 * a1 * gumda; + bigjdz = c01 * dumdz * z * a1 + c00 * z * nudz * a1 - c00 * z * a1 * a1 * gumdz; + } + + // equation 6.5 + z = std::exp(zeta + nu); + dum = 1.0e0_rt + z; + a1 = 1.0e0_rt/dum; + a2 = 1.0e0_rt/bigj; + + sreco = nu_constants::tfac6 * 2.649e-18_rt * sf.ye * std::pow(sf.zbar, 13.0_rt) * sf.den * bigj * a1; + if constexpr (do_derivatives) { + srecodt = sreco * (bigjdt * a2 - z * (zetadt + nudt) * a1); + srecoda = sreco * (-1.0e0_rt * sf.abari + bigjda * a2 - z * (zetada + nuda) * a1); + srecodz = sreco * (14.0e0_rt * sf.zbari + bigjdz * a2 - z * (zetadz + nudz) * a1); + } + } + +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void sneut5(const Real temp, const Real den, @@ -374,37 +541,6 @@ void sneut5(const Real temp, const Real den, fbremda,fbremdz,gbremda,gbremdz, fliqda,fliqdz,gliqda,gliqdz; - // recomb - Real nu,nudt,nu2,nu3,bigj,bigjdt,nuda,nudz,bigjda,bigjdz; - - - // numerical constants - constexpr Real fac1 = 5.0e0_rt * M_PI / 3.0e0_rt; - constexpr Real fac2 = 10.0e0_rt * M_PI; - constexpr Real fac3 = M_PI / 5.0e0_rt; - constexpr Real oneth = 1.0e0_rt/3.0e0_rt; - constexpr Real twoth = 2.0e0_rt/3.0e0_rt; - constexpr Real sixth = 1.0e0_rt/6.0e0_rt; - constexpr Real iln10 = 4.342944819032518e-1_rt; - - // theta is sin**2(theta_weinberg) = 0.2319 plus/minus 0.00005 (1996) - // xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998) - // change theta and xnufam if need be, and the changes will automatically - // propagate through the routine. cv and ca are the vector and axial currents. - - constexpr Real theta = 0.2319e0_rt; - constexpr Real xnufam = 3.0e0_rt; - constexpr Real cv = 0.5e0_rt + 2.0e0_rt * theta; - constexpr Real cvp = 1.0e0_rt - cv; - constexpr Real ca = 0.5e0_rt; - constexpr Real cap = 1.0e0_rt - ca; - constexpr Real tfac1 = cv*cv + ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp+cap*cap); - constexpr Real tfac2 = cv*cv - ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp - cap*cap); - constexpr Real tfac3 = tfac2/tfac1; - constexpr Real tfac4 = 0.5e0_rt * tfac1; - constexpr Real tfac5 = 0.5e0_rt * tfac2; - constexpr Real tfac6 = cv*cv + 1.5e0_rt*ca*ca + (xnufam - 1.0e0_rt)*(cvp*cvp + 1.5e0_rt*cap*cap); - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; @@ -442,11 +578,11 @@ void sneut5(const Real temp, const Real den, auto sf = get_sneut_factors(den, temp, abar, zbar); a0 = sf.rm * 1.0e-9_rt; - a1 = std::pow(a0, oneth); + a1 = std::pow(a0, nu_constants::oneth); zeta = a1 * sf.xlm1; if constexpr (do_derivatives) { - a2 = oneth * a1*sf.rmi * sf.xlm1; + a2 = nu_constants::oneth * a1*sf.rmi * sf.xlm1; zetadt = -a1 * sf.xlm2 * sf.xldt; zetada = a2 * sf.rmda; zetadz = a2 * sf.rmdz; @@ -564,12 +700,12 @@ void sneut5(const Real temp, const Real den, spairdz = gl*spairdz; } - a1 = tfac4*(1.0e0_rt + tfac3 * qpair); + a1 = nu_constants::tfac4*(1.0e0_rt + nu_constants::tfac3 * qpair); a3 = spair; spair = a1*a3; if constexpr (do_derivatives) { - a2 = tfac4*tfac3; + a2 = nu_constants::tfac4 * nu_constants::tfac3; spairdt = a1*spairdt + a2*qpairdt*a3; spairda = a1*spairda + a2*qpairda*a3; spairdz = a1*spairdz + a2*qpairdz*a3; @@ -580,7 +716,7 @@ void sneut5(const Real temp, const Real den, // equation 4.6 a1 = 1.019e-6_rt*sf.rm; - a2 = std::pow(a1, twoth); + a2 = std::pow(a1, nu_constants::twoth); b1 = std::sqrt(1.0e0_rt + a2); @@ -589,7 +725,7 @@ void sneut5(const Real temp, const Real den, gl2 = 1.1095e11_rt * sf.rm * c00; if constexpr (do_derivatives) { - a3 = twoth*a2/a1; + a3 = nu_constants::twoth*a2/a1; b2 = 1.0e0_rt/b1; gl2dt = -2.0e0_rt*gl2*sf.tempi; d = sf.rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt; @@ -634,15 +770,15 @@ void sneut5(const Real temp, const Real den, cc = std::log10(2.0e0_rt*sf.rm); xlnt = std::log10(temp); - xnum = sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt); - xden = sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt); + xnum = nu_constants::sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt); + xden = nu_constants::sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt); if constexpr (do_derivatives) { - xnumdt = -iln10*0.5e0_rt*sf.tempi; - a2 = iln10*sixth*sf.rmi; + xnumdt = -nu_constants::iln10*0.5e0_rt*sf.tempi; + a2 = nu_constants::iln10*nu_constants::sixth*sf.rmi; xnumda = a2*sf.rmda; xnumdz = a2*sf.rmdz; - xdendt = iln10*0.5e0_rt*sf.tempi; + xdendt = nu_constants::iln10*0.5e0_rt*sf.tempi; xdenda = a2*sf.rmda; xdendz = a2*sf.rmdz; } @@ -869,22 +1005,22 @@ void sneut5(const Real temp, const Real den, } - taudt = iln10*sf.tempi; + taudt = nu_constants::iln10*sf.tempi; // equation 3.7, compute the expensive trig functions only one time - cos1 = std::cos(fac1*tau); - cos2 = std::cos(fac1*2.0e0_rt*tau); - cos3 = std::cos(fac1*3.0e0_rt*tau); - cos4 = std::cos(fac1*4.0e0_rt*tau); - cos5 = std::cos(fac1*5.0e0_rt*tau); - last = std::cos(fac2*tau); - - sin1 = std::sin(fac1*tau); - sin2 = std::sin(fac1*2.0e0_rt*tau); - sin3 = std::sin(fac1*3.0e0_rt*tau); - sin4 = std::sin(fac1*4.0e0_rt*tau); - sin5 = std::sin(fac1*5.0e0_rt*tau); - xast = std::sin(fac2*tau); + cos1 = std::cos(nu_constants::fac1*tau); + cos2 = std::cos(nu_constants::fac1*2.0e0_rt*tau); + cos3 = std::cos(nu_constants::fac1*3.0e0_rt*tau); + cos4 = std::cos(nu_constants::fac1*4.0e0_rt*tau); + cos5 = std::cos(nu_constants::fac1*5.0e0_rt*tau); + last = std::cos(nu_constants::fac2*tau); + + sin1 = std::sin(nu_constants::fac1*tau); + sin2 = std::sin(nu_constants::fac1*2.0e0_rt*tau); + sin3 = std::sin(nu_constants::fac1*3.0e0_rt*tau); + sin4 = std::sin(nu_constants::fac1*4.0e0_rt*tau); + sin5 = std::sin(nu_constants::fac1*5.0e0_rt*tau); + xast = std::sin(nu_constants::fac2*tau); a0 = 0.5e0_rt*c00 + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2 @@ -902,21 +1038,21 @@ void sneut5(const Real temp, const Real den, + c25*cos5 + dd25*sin5 + 0.5e0_rt*c26*last; if constexpr (do_derivatives) { - f0 = taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0e0_rt + f0 = taudt*nu_constants::fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0e0_rt + dd02*cos2*2.0e0_rt - c03*sin3*3.0e0_rt + dd03*cos3*3.0e0_rt - c04*sin4*4.0e0_rt + dd04*cos4*4.0e0_rt - c05*sin5*5.0e0_rt + dd05*cos5*5.0e0_rt) - - 0.5e0_rt*c06*xast*fac2*taudt; + - 0.5e0_rt*c06*xast*nu_constants::fac2*taudt; - f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0e0_rt + f1 = taudt*nu_constants::fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0e0_rt + dd12*cos2*2.0e0_rt - c13*sin3*3.0e0_rt + dd13*cos3*3.0e0_rt - c14*sin4*4.0e0_rt + dd14*cos4*4.0e0_rt - c15*sin5*5.0e0_rt - + dd15*cos5*5.0e0_rt) - 0.5e0_rt*c16*xast*fac2*taudt; + + dd15*cos5*5.0e0_rt) - 0.5e0_rt*c16*xast*nu_constants::fac2*taudt; - f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0e0_rt + f2 = taudt*nu_constants::fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0e0_rt + dd22*cos2*2.0e0_rt - c23*sin3*3.0e0_rt + dd23*cos3*3.0e0_rt - c24*sin4*4.0e0_rt + dd24*cos4*4.0e0_rt - c25*sin5*5.0e0_rt - + dd25*cos5*5.0e0_rt) - 0.5e0_rt*c26*xast*fac2*taudt; + + dd25*cos5*5.0e0_rt) - 0.5e0_rt*c26*xast*nu_constants::fac2*taudt; } @@ -1002,12 +1138,12 @@ void sneut5(const Real temp, const Real den, sphotdz = sf.rm*sphotdz + sf.rmdz*a1; } - a1 = tfac4*(1.0e0_rt - tfac3 * qphot); + a1 = nu_constants::tfac4*(1.0e0_rt - nu_constants::tfac3 * qphot); a3 = sphot; sphot = a1*a3; if constexpr (do_derivatives) { - a2 = -tfac4*tfac3; + a2 = -nu_constants::tfac4*nu_constants::tfac3; sphotdt = a1*sphotdt + a2*qphotdt*a3; sphotda = a1*sphotda + a2*qphotda*a3; sphotdz = a1*sphotdz + a2*qphotdz*a3; @@ -1043,7 +1179,7 @@ void sneut5(const Real temp, const Real den, t8m6 = t8m5*t8m1; - tfermi = 5.9302e9_rt*(std::sqrt(1.0e0_rt+1.018e0_rt*std::pow(den6*sf.ye, twoth))-1.0e0_rt); + tfermi = 5.9302e9_rt*(std::sqrt(1.0e0_rt+1.018e0_rt*std::pow(den6*sf.ye, nu_constants::twoth))-1.0e0_rt); // "weak" degenerate electrons only if (temp > 0.3e0_rt * tfermi) { @@ -1166,12 +1302,12 @@ void sneut5(const Real temp, const Real den, dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; } - z = tfac4*fbrem - tfac5*gbrem; + z = nu_constants::tfac4*fbrem - nu_constants::tfac5*gbrem; sbrem = dum * z; if constexpr (do_derivatives) { - sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt); - sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda); - sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz); + sbremdt = dumdt*z + dum*(nu_constants::tfac4*fbremdt - nu_constants::tfac5*gbremdt); + sbremda = dumda*z + dum*(nu_constants::tfac4*fbremda - nu_constants::tfac5*gbremda); + sbremdz = dumdz*z + dum*(nu_constants::tfac4*fbremdz - nu_constants::tfac5*gbremdz); } // liquid metal with c12 parameters (not too different for other elements) @@ -1179,7 +1315,7 @@ void sneut5(const Real temp, const Real den, } else { - u = fac3 * (std::log10(den) - 3.0e0_rt); + u = nu_constants::fac3 * (std::log10(den) - 3.0e0_rt); //a0 = iln10*fac3*sf.deni; // compute the expensive trig functions of equation 5.21 only once @@ -1255,26 +1391,26 @@ void sneut5(const Real temp, const Real den, // - 0.00031e0_rt*sin4*4.0e0_rt - 0.00018e0_rt*cos4*4.0e0_rt // - 0.00069e0_rt*sin5*5.0e0_rt); - dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*sf.abari, oneth); + dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*sf.abari, nu_constants::oneth); if constexpr (do_derivatives) { dumdt = -dum*sf.tempi; - dumda = -oneth*dum*sf.abari; + dumda = -nu_constants::oneth*dum*sf.abari; dumdz = 2.0e0_rt*dum*sf.zbari; } gm1 = 1.0e0_rt/dum; gm2 = gm1*gm1; - gm13 = std::pow(gm1, oneth); + gm13 = std::pow(gm1, nu_constants::oneth); gm23 = gm13 * gm13; gm43 = gm13*gm1; gm53 = gm23*gm1; // equation 5.25 and 5.26 v = -0.05483e0_rt - 0.01946e0_rt*gm13 + 1.86310e0_rt*gm23 - 0.78873e0_rt*gm1; - a0 = oneth*0.01946e0_rt*gm43 - twoth*1.86310e0_rt*gm53 + 0.78873e0_rt*gm2; + a0 = nu_constants::oneth*0.01946e0_rt*gm43 - nu_constants::twoth*1.86310e0_rt*gm53 + 0.78873e0_rt*gm2; w = -0.06711e0_rt + 0.06859e0_rt*gm13 + 1.74360e0_rt*gm23 - 0.74498e0_rt*gm1; - a1 = -oneth*0.06859e0_rt*gm43 - twoth*1.74360e0_rt*gm53 + 0.74498e0_rt*gm2; + a1 = -nu_constants::oneth*0.06859e0_rt*gm43 - nu_constants::twoth*1.74360e0_rt*gm53 + 0.74498e0_rt*gm2; // equation 5.19 and 5.20 fliq = v*fb + (1.0e0_rt - v)*ft; @@ -1299,118 +1435,16 @@ void sneut5(const Real temp, const Real den, dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; } - z = tfac4*fliq - tfac5*gliq; + z = nu_constants::tfac4*fliq - nu_constants::tfac5*gliq; sbrem = dum * z; if constexpr (do_derivatives) { - sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt); - sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda); - sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz); + sbremdt = dumdt*z + dum*(nu_constants::tfac4*fliqdt - nu_constants::tfac5*gliqdt); + sbremda = dumda*z + dum*(nu_constants::tfac4*fliqda - nu_constants::tfac5*gliqda); + sbremdz = dumdz*z + dum*(nu_constants::tfac4*fliqdz - nu_constants::tfac5*gliqdz); } } - // recombination neutrino section - // for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e - // equation 6.11 solved for nu - xnum = 1.10520e8_rt * den * sf.ye /(temp*std::sqrt(temp)); - if constexpr (do_derivatives) { - xnumdt = -1.50e0_rt*xnum*sf.tempi; - xnumda = -xnum*sf.abari; - xnumdz = xnum*sf.zbari; - } - - // the chemical potential - nu = ifermi12(xnum); - - // a0 is d(nu)/d(xnum) - a0 = 1.0e0_rt/(0.5e0_rt*zfermim12(nu)); - if constexpr (do_derivatives) { - nudt = a0*xnumdt; - nuda = a0*xnumda; - nudz = a0*xnumdz; - } - nu2 = nu * nu; - nu3 = nu2 * nu; - - // table 12 - if (nu >= -20.0_rt && nu < 0.0_rt) { - a1 = 1.51e-2_rt; - a2 = 2.42e-1_rt; - a3 = 1.21e0_rt; - b = 3.71e-2_rt; - c = 9.06e-1_rt; - d = 9.28e-1_rt; - f1 = 0.0e0_rt; - f2 = 0.0e0_rt; - f3 = 0.0e0_rt; - } else if (nu >= 0.0_rt && nu <= 10.0_rt) { - a1 = 1.23e-2_rt; - a2 = 2.66e-1_rt; - a3 = 1.30e0_rt; - b = 1.17e-1_rt; - c = 8.97e-1_rt; - d = 1.77e-1_rt; - f1 = -1.20e-2_rt; - f2 = 2.29e-2_rt; - f3 = -1.04e-3_rt; - } - - // equation 6.7, 6.13 and 6.14 - if (nu >= -20.0_rt && nu <= 10.0_rt) { - - zeta = 1.579e5_rt*zbar*zbar*sf.tempi; - if constexpr (do_derivatives) { - zetadt = -zeta*sf.tempi; - zetada = 0.0e0_rt; - zetadz = 2.0e0_rt*zeta*sf.zbari; - } - - c00 = 1.0e0_rt/(1.0e0_rt + f1*nu + f2*nu2 + f3*nu3); - c01 = f1 + f2*2.0e0_rt*nu + f3*3.0e0_rt*nu2; - dum = zeta*c00; - if constexpr (do_derivatives) { - dumdt = zetadt*c00 + zeta*c01*nudt; - dumda = zeta*c01*nuda; - dumdz = zetadz*c00 + zeta*c01*nudz; - } - - z = 1.0e0_rt/dum; - dd00 = std::pow(dum, -2.25_rt); - dd01 = std::pow(dum, -4.55_rt); - c00 = a1*z + a2*dd00 + a3*dd01; - c01 = -(a1*z + 2.25_rt*a2*dd00 + 4.55_rt*a3*dd01)*z; - - z = std::exp(c*nu); - dd00 = b*z*(1.0e0_rt + d*dum); - gum = 1.0e0_rt + dd00; - if constexpr (do_derivatives) { - gumdt = dd00*c*nudt + b*z*d*dumdt; - gumda = dd00*c*nuda + b*z*d*dumda; - gumdz = dd00*c*nudz + b*z*d*dumdz; - } - - z = std::exp(nu); - a1 = 1.0e0_rt/gum; - - bigj = c00 * z * a1; - if constexpr (do_derivatives) { - bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt; - bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda; - bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz; - } - - // equation 6.5 - z = std::exp(zeta + nu); - dum = 1.0e0_rt + z; - a1 = 1.0e0_rt/dum; - a2 = 1.0e0_rt/bigj; - - sreco = tfac6 * 2.649e-18_rt * sf.ye * std::pow(zbar, 13.0_rt) * den * bigj*a1; - if constexpr (do_derivatives) { - srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1); - srecoda = sreco*(-1.0e0_rt*sf.abari + bigjda*a2 - z*(zetada+nuda)*a1); - srecodz = sreco*(14.0e0_rt*sf.zbari + bigjdz*a2 - z*(zetadz+nudz)*a1); - } - } + nu_recomb(sf, sreco, srecodt, srecoda, srecodz); // convert from erg/cm^3/s to erg/g/s // comment these out to duplicate the itoh et al plots From 5ae7e6f29f6440002bb0d72636364ea974ce7d3b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 19:03:52 -0400 Subject: [PATCH 06/10] fix undefined --- neutrinos/sneut5.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 726f1269c6..c154b11404 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -504,9 +504,9 @@ void sneut5(const Real temp, const Real den, a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06, c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, - dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d{0.0},f0, - f1{0.0},f2,f3,z, - dum,dumdt,gum,gumdt,dumda,dumdz,gumda,gumdz; + dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,c,d{0.0},f0, + f1{0.0},f2,z, + dum,dumdt,gum,dumda,dumdz; // pair production Real gl,gldt, From 260fea4566c9f2d559ca1271adc8a067b4590efd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 19:42:59 -0400 Subject: [PATCH 07/10] move sneut bremsstrahlung into its own function --- neutrinos/sneut5.H | 602 +++++++++++++++++++++++---------------------- 1 file changed, 304 insertions(+), 298 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index c154b11404..d5ac326e93 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -35,7 +35,6 @@ namespace nu_constants { constexpr Real tfac5 = 0.5e0_rt * tfac2; constexpr Real tfac6 = cv*cv + 1.5e0_rt*ca*ca + (xnufam - 1.0e0_rt)*(cvp*cvp + 1.5e0_rt*cap*cap); - } @@ -350,6 +349,309 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_brem(const sneutf_t& sf, + Real& sbrem, Real& sbremdt, Real& sbremda, Real& sbremdz) { + + // bremsstrahlung neutrino section + // for reactions like e- + (z,a) => e- + (z,a) + nu + nubar + // n + n => n + n + nu + nubar + // n + p => n + p + nu + nubar + // equation 4.3 + + Real den6 = sf.den * 1.0e-6_rt; + Real t8 = sf.temp * 1.0e-8_rt; + Real t812 = std::sqrt(t8); + Real t832 = t8 * t812; + Real t82 = t8*t8; + Real t83 = t82*t8; + Real t85 = t82*t83; + Real t86 = t85*t8; + Real t8m1 = 1.0e0_rt/t8; + Real t8m2 = t8m1*t8m1; + Real t8m3 = t8m2*t8m1; + Real t8m5 = t8m3*t8m2; + Real t8m6 = t8m5*t8m1; + + Real tfermi = 5.9302e9_rt * (std::sqrt(1.0e0_rt + 1.018e0_rt * std::pow(den6 * sf.ye, nu_constants::twoth)) - 1.0e0_rt); + + // "weak" degenerate electrons only + if (sf.temp > 0.3e0_rt * tfermi) { + + // equation 5.3 + Real dum = 7.05e6_rt * t832 + 5.12e4_rt * t83; + Real dumdt; + if constexpr (do_derivatives) { + dumdt = (1.5e0_rt * 7.05e6_rt * t812 + 3.0e0_rt * 5.12e4_rt * t82) * 1.0e-8_rt; + } + + Real z = 1.0e0_rt / dum; + Real eta = sf.rm * z; + Real etadt, etada, etadz; + if constexpr (do_derivatives) { + etadt = -sf.rm*z*z*dumdt; + etada = sf.rmda*z; + etadz = sf.rmdz*z; + } + + Real etam1 = 1.0e0_rt/eta; + Real etam2 = etam1 * etam1; + Real etam3 = etam2 * etam1; + + // equation 5.2 + Real a0 = 23.5e0_rt + 6.83e4_rt * t8m2 + 7.81e8_rt * t8m5; + Real f0 = (-2.0e0_rt * 6.83e4_rt * t8m3 - 5.0e0_rt * 7.81e8_rt * t8m6) * 1.0e-8_rt; + Real xnum = 1.0e0_rt / a0; + + dum = 1.0e0_rt + 1.47e0_rt * etam1 + 3.29e-2_rt * etam2; + Real dumda, dumdz; + if constexpr (do_derivatives) { + z = -1.47e0_rt * etam2 - 2.0e0_rt * 3.29e-2_rt * etam3; + dumdt = z*etadt; + dumda = z*etada; + dumdz = z*etadz; + } + + Real c00 = 1.26e0_rt * (1.0e0_rt + etam1); + Real c01, c02, c03, c04; + if constexpr (do_derivatives) { + z = -1.26e0_rt*etam2; + c01 = z*etadt; + c03 = z*etada; + c04 = z*etadz; + } + + z = 1.0e0_rt/dum; + Real xden = c00 * z; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + xdendt = (c01 - xden * dumdt) * z; + xdenda = (c03 - xden * dumda) * z; + xdendz = (c04 - xden * dumdz) * z; + } + + Real fbrem = xnum + xden; + Real fbremdt, fbremda, fbremdz; + if constexpr (do_derivatives) { + fbremdt = -xnum*xnum*f0 + xdendt; + fbremda = xdenda; + fbremdz = xdendz; + } + + // equation 5.9 + a0 = 230.0e0_rt + 6.7e5_rt * t8m2 + 7.66e9_rt * t8m5; + f0 = (-2.0e0_rt * 6.7e5_rt * t8m3 - 5.0e0_rt * 7.66e9_rt * t8m6) * 1.0e-8_rt; + + z = 1.0e0_rt + sf.rm * 1.0e-9_rt; + dum = a0 * z; + if constexpr (do_derivatives) { + dumdt = f0 * z; + z = a0 * 1.0e-9_rt; + dumda = z * sf.rmda; + dumdz = z * sf.rmdz; + } + + xnum = 1.0e0_rt / dum; + Real xnumdt, xnumda, xnumdz; + if constexpr (do_derivatives) { + z = -xnum * xnum; + xnumdt = z * dumdt; + xnumda = z * dumda; + xnumdz = z * dumdz; + } + + c00 = 7.75e5_rt * t832 + 247.0e0_rt * std::pow(t8, 3.85e0_rt); + c01 = 4.07e0_rt + 0.0240e0_rt * std::pow(t8, 1.4e0_rt); + c02 = 4.59e-5_rt * std::pow(t8, -0.110e0_rt); + + Real dd00, dd01, dd02; + if constexpr (do_derivatives) { + dd00 = (1.5e0_rt * 7.75e5_rt * t812 + 3.85e0_rt * 247.0e0_rt * + std::pow(t8, 2.85e0_rt)) * 1.0e-8_rt; + + dd01 = 1.4e0_rt * 0.0240e0_rt * std::pow(t8, 0.4e0_rt) * 1.0e-8_rt; + dd02 = -0.11e0_rt * 4.59e-5_rt * std::pow(t8, -1.11e0_rt) * 1.0e-8_rt; + } + + z = std::pow(sf.den, 0.656e0_rt); + dum = c00 * sf.rmi + c01 + c02 * z; + if constexpr (do_derivatives) { + dumdt = dd00 * sf.rmi + dd01 + dd02 * z; + z = -c00 * sf.rmi * sf.rmi; + dumda = z * sf.rmda; + dumdz = z * sf.rmdz; + } + + xden = 1.0e0_rt / dum; + if constexpr (do_derivatives) { + z = -xden * xden; + xdendt = z * dumdt; + xdenda = z * dumda; + xdendz = z * dumdz; + } + + Real gbrem = xnum + xden; + Real gbremdt, gbremda, gbremdz; + if constexpr (do_derivatives) { + gbremdt = xnumdt + xdendt; + gbremda = xnumda + xdenda; + gbremdz = xnumdz + xdendz; + } + + // equation 5.1 + dum = 0.5738e0_rt * sf.zbar * sf.ye * t86 * sf.den; + if constexpr (do_derivatives) { + dumdt = 0.5738e0_rt * sf.zbar * sf.ye * 6.0e0_rt * t85 * sf.den * 1.0e-8_rt; + dumda = -dum * sf.abari; + dumdz = 0.5738e0_rt * 2.0e0_rt * sf.ye * t86 * sf.den; + } + + z = nu_constants::tfac4 * fbrem - nu_constants::tfac5 * gbrem; + sbrem = dum * z; + if constexpr (do_derivatives) { + sbremdt = dumdt * z + dum * (nu_constants::tfac4 * fbremdt - nu_constants::tfac5 * gbremdt); + sbremda = dumda * z + dum * (nu_constants::tfac4 * fbremda - nu_constants::tfac5 * gbremda); + sbremdz = dumdz * z + dum * (nu_constants::tfac4 * fbremdz - nu_constants::tfac5 * gbremdz); + } + + // liquid metal with c12 parameters (not too different for other elements) + // equation 5.18 and 5.16 + + } else { + + Real u = nu_constants::fac3 * (std::log10(sf.den) - 3.0e0_rt); + //a0 = iln10*fac3*sf.deni; + + // compute the expensive trig functions of equation 5.21 only once + Real cos1 = std::cos(u); + Real cos2 = std::cos(2.0e0_rt*u); + Real cos3 = std::cos(3.0e0_rt*u); + Real cos4 = std::cos(4.0e0_rt*u); + Real cos5 = std::cos(5.0e0_rt*u); + + Real sin1 = std::sin(u); + Real sin2 = std::sin(2.0e0_rt*u); + Real sin3 = std::sin(3.0e0_rt*u); + Real sin4 = std::sin(4.0e0_rt*u); + //sin5 = std::sin(5.0e0_rt*u); + + // equation 5.21 + Real fb = 0.5e0_rt * 0.17946e0_rt + 0.00945e0_rt * u + 0.34529e0_rt + - 0.05821e0_rt * cos1 - 0.04969e0_rt * sin1 + - 0.01089e0_rt * cos2 - 0.01584e0_rt * sin2 + - 0.01147e0_rt * cos3 - 0.00504e0_rt * sin3 + - 0.00656e0_rt * cos4 - 0.00281e0_rt * sin4 + - 0.00519e0_rt * cos5; + + // c00 = a0*(0.00945e0_rt + // + 0.05821e0_rt*sin1 - 0.04969e0_rt*cos1 + // + 0.01089e0_rt*sin2*2.0e0_rt - 0.01584e0_rt*cos2*2.0e0_rt + // + 0.01147e0_rt*sin3*3.0e0_rt - 0.00504e0_rt*cos3*3.0e0_rt + // + 0.00656e0_rt*sin4*4.0e0_rt - 0.00281e0_rt*cos4*4.0e0_rt + // + 0.00519e0_rt*sin5*5.0e0_rt); + + // equation 5.22 + Real ft = 0.5e0_rt * 0.06781e0_rt - 0.02342e0_rt * u + 0.24819e0_rt + - 0.00944e0_rt * cos1 - 0.02213e0_rt * sin1 + - 0.01289e0_rt * cos2 - 0.01136e0_rt * sin2 + - 0.00589e0_rt * cos3 - 0.00467e0_rt * sin3 + - 0.00404e0_rt * cos4 - 0.00131e0_rt * sin4 + - 0.00330e0_rt * cos5; + + // c01 = a0*(-0.02342e0_rt + // + 0.00944e0_rt*sin1 - 0.02213e0_rt*cos1 + // + 0.01289e0_rt*sin2*2.0e0_rt - 0.01136e0_rt*cos2*2.0e0_rt + // + 0.00589e0_rt*sin3*3.0e0_rt - 0.00467e0_rt*cos3*3.0e0_rt + // + 0.00404e0_rt*sin4*4.0e0_rt - 0.00131e0_rt*cos4*4.0e0_rt + // + 0.00330e0_rt*sin5*5.0e0_rt); + + // equation 5.23 + Real gb = 0.5e0_rt * 0.00766e0_rt - 0.01259e0_rt * u + 0.07917e0_rt + - 0.00710e0_rt * cos1 + 0.02300e0_rt * sin1 + - 0.00028e0_rt * cos2 - 0.01078e0_rt * sin2 + + 0.00232e0_rt * cos3 + 0.00118e0_rt * sin3 + + 0.00044e0_rt * cos4 - 0.00089e0_rt * sin4 + + 0.00158e0_rt * cos5; + + // c02 = a0*(-0.01259e0_rt + // + 0.00710e0_rt*sin1 + 0.02300e0_rt*cos1 + // + 0.00028e0_rt*sin2*2.0e0_rt - 0.01078e0_rt*cos2*2.0e0_rt + // - 0.00232e0_rt*sin3*3.0e0_rt + 0.00118e0_rt*cos3*3.0e0_rt + // - 0.00044e0_rt*sin4*4.0e0_rt - 0.00089e0_rt*cos4*4.0e0_rt + // - 0.00158e0_rt*sin5*5.0e0_rt); + + // equation 5.24 + Real gt = -0.5e0_rt * 0.00769e0_rt - 0.00829e0_rt * u + 0.05211e0_rt + + 0.00356e0_rt * cos1 + 0.01052e0_rt * sin1 + - 0.00184e0_rt * cos2 - 0.00354e0_rt * sin2 + + 0.00146e0_rt * cos3 - 0.00014e0_rt * sin3 + + 0.00031e0_rt * cos4 - 0.00018e0_rt * sin4 + + 0.00069e0_rt * cos5; + + // c03 = a0*(-0.00829e0_rt + // - 0.00356e0_rt*sin1 + 0.01052e0_rt*cos1 + // + 0.00184e0_rt*sin2*2.0e0_rt - 0.00354e0_rt*cos2*2.0e0_rt + // - 0.00146e0_rt*sin3*3.0e0_rt - 0.00014e0_rt*cos3*3.0e0_rt + // - 0.00031e0_rt*sin4*4.0e0_rt - 0.00018e0_rt*cos4*4.0e0_rt + // - 0.00069e0_rt*sin5*5.0e0_rt); + + Real dum = 2.275e-1_rt * sf.zbar * sf.zbar * t8m1 * std::pow(den6*sf.abari, nu_constants::oneth); + Real dumdt, dumda, dumdz; + if constexpr (do_derivatives) { + dumdt = -dum*sf.tempi; + dumda = -nu_constants::oneth*dum*sf.abari; + dumdz = 2.0e0_rt*dum*sf.zbari; + } + + Real gm1 = 1.0e0_rt / dum; + Real gm2 = gm1 * gm1; + Real gm13 = std::pow(gm1, nu_constants::oneth); + Real gm23 = gm13 * gm13; + Real gm43 = gm13 * gm1; + Real gm53 = gm23 * gm1; + + // equation 5.25 and 5.26 + Real v = -0.05483e0_rt - 0.01946e0_rt * gm13 + 1.86310e0_rt * gm23 - 0.78873e0_rt * gm1; + Real a0 = nu_constants::oneth * 0.01946e0_rt * gm43 - nu_constants::twoth * 1.86310e0_rt * gm53 + 0.78873e0_rt * gm2; + + Real w = -0.06711e0_rt + 0.06859e0_rt * gm13 + 1.74360e0_rt * gm23 - 0.74498e0_rt * gm1; + Real a1 = -nu_constants::oneth*0.06859e0_rt * gm43 - nu_constants::twoth * 1.74360e0_rt * gm53 + 0.74498e0_rt * gm2; + + // equation 5.19 and 5.20 + Real fliq = v * fb + (1.0e0_rt - v) * ft; + Real fliqdt, fliqda, fliqdz; + if constexpr (do_derivatives) { + fliqdt = a0 * dumdt * (fb - ft); + fliqda = a0 * dumda * (fb - ft); + fliqdz = a0 * dumdz * (fb - ft); + } + + Real gliq = w * gb + (1.0e0_rt - w) * gt; + Real gliqdt, gliqda, gliqdz; + if constexpr (do_derivatives) { + gliqdt = a1 * dumdt*(gb - gt); + gliqda = a1 * dumda*(gb - gt); + gliqdz = a1 * dumdz*(gb - gt); + } + + // equation 5.17 + dum = 0.5738e0_rt * sf.zbar * sf.ye * t86 * sf.den; + if constexpr (do_derivatives) { + dumdt = 0.5738e0_rt * sf.zbar * sf.ye * 6.0e0_rt * t85 * sf.den * 1.0e-8_rt; + dumda = -dum * sf.abari; + dumdz = 0.5738e0_rt * 2.0e0_rt * sf.ye * t86 * sf.den; + } + + Real z = nu_constants::tfac4*fliq - nu_constants::tfac5*gliq; + sbrem = dum * z; + if constexpr (do_derivatives) { + sbremdt = dumdt*z + dum*(nu_constants::tfac4*fliqdt - nu_constants::tfac5*gliqdt); + sbremda = dumda*z + dum*(nu_constants::tfac4*fliqda - nu_constants::tfac5*gliqda); + sbremdz = dumdz*z + dum*(nu_constants::tfac4*fliqdz - nu_constants::tfac5*gliqdz); + } + } +} template AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -529,18 +831,6 @@ void sneut5(const Real temp, const Real den, fphot,fphotdt,qphot,qphotdt, fphotda,fphotdz,qphotda,qphotdz; - // brem - Real t8,t812,t832,t82,t83,t85,t86,t8m1,t8m2,t8m3,t8m5, - t8m6, - eta,etadt,etam1,etam2,etam3, - fbrem,fbremdt, - gbrem,gbremdt, - u,gm1,gm2,gm13,gm23,gm43,gm53,v,w,fb,gt,gb, - fliq,fliqdt,gliq,gliqdt, - etada,etadz, - fbremda,fbremdz,gbremda,gbremdz, - fliqda,fliqdz,gliqda,gliqdz; - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; @@ -1158,291 +1448,7 @@ void sneut5(const Real temp, const Real den, } } - // bremsstrahlung neutrino section - // for reactions like e- + (z,a) => e- + (z,a) + nu + nubar - // n + n => n + n + nu + nubar - // n + p => n + p + nu + nubar - // equation 4.3 - - den6 = den * 1.0e-6_rt; - t8 = temp * 1.0e-8_rt; - t812 = std::sqrt(t8); - t832 = t8 * t812; - t82 = t8*t8; - t83 = t82*t8; - t85 = t82*t83; - t86 = t85*t8; - t8m1 = 1.0e0_rt/t8; - t8m2 = t8m1*t8m1; - t8m3 = t8m2*t8m1; - t8m5 = t8m3*t8m2; - t8m6 = t8m5*t8m1; - - - tfermi = 5.9302e9_rt*(std::sqrt(1.0e0_rt+1.018e0_rt*std::pow(den6*sf.ye, nu_constants::twoth))-1.0e0_rt); - - // "weak" degenerate electrons only - if (temp > 0.3e0_rt * tfermi) { - - // equation 5.3 - dum = 7.05e6_rt * t832 + 5.12e4_rt * t83; - if constexpr (do_derivatives) { - dumdt = (1.5e0_rt*7.05e6_rt*t812 + 3.0e0_rt*5.12e4_rt*t82)*1.0e-8_rt; - } - - z = 1.0e0_rt/dum; - eta = sf.rm*z; - if constexpr (do_derivatives) { - etadt = -sf.rm*z*z*dumdt; - etada = sf.rmda*z; - etadz = sf.rmdz*z; - } - - etam1 = 1.0e0_rt/eta; - etam2 = etam1 * etam1; - etam3 = etam2 * etam1; - - // equation 5.2 - a0 = 23.5e0_rt + 6.83e4_rt*t8m2 + 7.81e8_rt*t8m5; - f0 = (-2.0e0_rt*6.83e4_rt*t8m3 - 5.0e0_rt*7.81e8_rt*t8m6)*1.0e-8_rt; - xnum = 1.0e0_rt/a0; - - dum = 1.0e0_rt + 1.47e0_rt*etam1 + 3.29e-2_rt*etam2; - if constexpr (do_derivatives) { - z = -1.47e0_rt*etam2 - 2.0e0_rt*3.29e-2_rt*etam3; - dumdt = z*etadt; - dumda = z*etada; - dumdz = z*etadz; - } - - c00 = 1.26e0_rt * (1.0e0_rt+etam1); - if constexpr (do_derivatives) { - z = -1.26e0_rt*etam2; - c01 = z*etadt; - c03 = z*etada; - c04 = z*etadz; - } - - z = 1.0e0_rt/dum; - xden = c00*z; - if constexpr (do_derivatives) { - xdendt = (c01 - xden*dumdt)*z; - xdenda = (c03 - xden*dumda)*z; - xdendz = (c04 - xden*dumdz)*z; - } - - fbrem = xnum + xden; - if constexpr (do_derivatives) { - fbremdt = -xnum*xnum*f0 + xdendt; - fbremda = xdenda; - fbremdz = xdendz; - } - - // equation 5.9 - a0 = 230.0e0_rt + 6.7e5_rt*t8m2 + 7.66e9_rt*t8m5; - f0 = (-2.0e0_rt*6.7e5_rt*t8m3 - 5.0e0_rt*7.66e9_rt*t8m6)*1.0e-8_rt; - - z = 1.0e0_rt + sf.rm*1.0e-9_rt; - dum = a0*z; - if constexpr (do_derivatives) { - dumdt = f0*z; - z = a0*1.0e-9_rt; - dumda = z*sf.rmda; - dumdz = z*sf.rmdz; - } - - xnum = 1.0e0_rt/dum; - if constexpr (do_derivatives) { - z = -xnum*xnum; - xnumdt = z*dumdt; - xnumda = z*dumda; - xnumdz = z*dumdz; - } - - c00 = 7.75e5_rt*t832 + 247.0e0_rt * std::pow(t8, 3.85e0_rt); - c01 = 4.07e0_rt + 0.0240e0_rt * std::pow(t8, 1.4e0_rt); - c02 = 4.59e-5_rt * std::pow(t8, -0.110e0_rt); - - if constexpr (do_derivatives) { - dd00 = (1.5e0_rt*7.75e5_rt*t812 + 3.85e0_rt*247.0e0_rt*std::pow(t8, 2.85e0_rt))*1.0e-8_rt; - - dd01 = 1.4e0_rt*0.0240e0_rt * std::pow(t8, 0.4e0_rt)*1.0e-8_rt; - dd02 = -0.11e0_rt*4.59e-5_rt * std::pow(t8, -1.11e0_rt)*1.0e-8_rt; - } - - z = std::pow(den, 0.656e0_rt); - dum = c00*sf.rmi + c01 + c02*z; - if constexpr (do_derivatives) { - dumdt = dd00*sf.rmi + dd01 + dd02*z; - z = -c00*sf.rmi*sf.rmi; - dumda = z*sf.rmda; - dumdz = z*sf.rmdz; - } - - xden = 1.0e0_rt/dum; - if constexpr (do_derivatives) { - z = -xden*xden; - xdendt = z*dumdt; - xdenda = z*dumda; - xdendz = z*dumdz; - } - - gbrem = xnum + xden; - if constexpr (do_derivatives) { - gbremdt = xnumdt + xdendt; - gbremda = xnumda + xdenda; - gbremdz = xnumdz + xdendz; - } - - // equation 5.1 - dum = 0.5738e0_rt*zbar*sf.ye*t86*den; - if constexpr (do_derivatives) { - dumdt = 0.5738e0_rt*zbar*sf.ye*6.0e0_rt*t85*den*1.0e-8_rt; - dumda = -dum*sf.abari; - dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; - } - - z = nu_constants::tfac4*fbrem - nu_constants::tfac5*gbrem; - sbrem = dum * z; - if constexpr (do_derivatives) { - sbremdt = dumdt*z + dum*(nu_constants::tfac4*fbremdt - nu_constants::tfac5*gbremdt); - sbremda = dumda*z + dum*(nu_constants::tfac4*fbremda - nu_constants::tfac5*gbremda); - sbremdz = dumdz*z + dum*(nu_constants::tfac4*fbremdz - nu_constants::tfac5*gbremdz); - } - - // liquid metal with c12 parameters (not too different for other elements) - // equation 5.18 and 5.16 - - } else { - - u = nu_constants::fac3 * (std::log10(den) - 3.0e0_rt); - //a0 = iln10*fac3*sf.deni; - - // compute the expensive trig functions of equation 5.21 only once - cos1 = std::cos(u); - cos2 = std::cos(2.0e0_rt*u); - cos3 = std::cos(3.0e0_rt*u); - cos4 = std::cos(4.0e0_rt*u); - cos5 = std::cos(5.0e0_rt*u); - - sin1 = std::sin(u); - sin2 = std::sin(2.0e0_rt*u); - sin3 = std::sin(3.0e0_rt*u); - sin4 = std::sin(4.0e0_rt*u); - //sin5 = std::sin(5.0e0_rt*u); - - // equation 5.21 - fb = 0.5e0_rt * 0.17946e0_rt + 0.00945e0_rt*u + 0.34529e0_rt - - 0.05821e0_rt*cos1 - 0.04969e0_rt*sin1 - - 0.01089e0_rt*cos2 - 0.01584e0_rt*sin2 - - 0.01147e0_rt*cos3 - 0.00504e0_rt*sin3 - - 0.00656e0_rt*cos4 - 0.00281e0_rt*sin4 - - 0.00519e0_rt*cos5; - - // c00 = a0*(0.00945e0_rt - // + 0.05821e0_rt*sin1 - 0.04969e0_rt*cos1 - // + 0.01089e0_rt*sin2*2.0e0_rt - 0.01584e0_rt*cos2*2.0e0_rt - // + 0.01147e0_rt*sin3*3.0e0_rt - 0.00504e0_rt*cos3*3.0e0_rt - // + 0.00656e0_rt*sin4*4.0e0_rt - 0.00281e0_rt*cos4*4.0e0_rt - // + 0.00519e0_rt*sin5*5.0e0_rt); - - // equation 5.22 - ft = 0.5e0_rt * 0.06781e0_rt - 0.02342e0_rt*u + 0.24819e0_rt - - 0.00944e0_rt*cos1 - 0.02213e0_rt*sin1 - - 0.01289e0_rt*cos2 - 0.01136e0_rt*sin2 - - 0.00589e0_rt*cos3 - 0.00467e0_rt*sin3 - - 0.00404e0_rt*cos4 - 0.00131e0_rt*sin4 - - 0.00330e0_rt*cos5; - - // c01 = a0*(-0.02342e0_rt - // + 0.00944e0_rt*sin1 - 0.02213e0_rt*cos1 - // + 0.01289e0_rt*sin2*2.0e0_rt - 0.01136e0_rt*cos2*2.0e0_rt - // + 0.00589e0_rt*sin3*3.0e0_rt - 0.00467e0_rt*cos3*3.0e0_rt - // + 0.00404e0_rt*sin4*4.0e0_rt - 0.00131e0_rt*cos4*4.0e0_rt - // + 0.00330e0_rt*sin5*5.0e0_rt); - - // equation 5.23 - gb = 0.5e0_rt * 0.00766e0_rt - 0.01259e0_rt*u + 0.07917e0_rt - - 0.00710e0_rt*cos1 + 0.02300e0_rt*sin1 - - 0.00028e0_rt*cos2 - 0.01078e0_rt*sin2 - + 0.00232e0_rt*cos3 + 0.00118e0_rt*sin3 - + 0.00044e0_rt*cos4 - 0.00089e0_rt*sin4 - + 0.00158e0_rt*cos5; - - // c02 = a0*(-0.01259e0_rt - // + 0.00710e0_rt*sin1 + 0.02300e0_rt*cos1 - // + 0.00028e0_rt*sin2*2.0e0_rt - 0.01078e0_rt*cos2*2.0e0_rt - // - 0.00232e0_rt*sin3*3.0e0_rt + 0.00118e0_rt*cos3*3.0e0_rt - // - 0.00044e0_rt*sin4*4.0e0_rt - 0.00089e0_rt*cos4*4.0e0_rt - // - 0.00158e0_rt*sin5*5.0e0_rt); - - // equation 5.24 - gt = -0.5e0_rt * 0.00769e0_rt - 0.00829e0_rt*u + 0.05211e0_rt - + 0.00356e0_rt*cos1 + 0.01052e0_rt*sin1 - - 0.00184e0_rt*cos2 - 0.00354e0_rt*sin2 - + 0.00146e0_rt*cos3 - 0.00014e0_rt*sin3 - + 0.00031e0_rt*cos4 - 0.00018e0_rt*sin4 - + 0.00069e0_rt*cos5; - - // c03 = a0*(-0.00829e0_rt - // - 0.00356e0_rt*sin1 + 0.01052e0_rt*cos1 - // + 0.00184e0_rt*sin2*2.0e0_rt - 0.00354e0_rt*cos2*2.0e0_rt - // - 0.00146e0_rt*sin3*3.0e0_rt - 0.00014e0_rt*cos3*3.0e0_rt - // - 0.00031e0_rt*sin4*4.0e0_rt - 0.00018e0_rt*cos4*4.0e0_rt - // - 0.00069e0_rt*sin5*5.0e0_rt); - - dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*sf.abari, nu_constants::oneth); - if constexpr (do_derivatives) { - dumdt = -dum*sf.tempi; - dumda = -nu_constants::oneth*dum*sf.abari; - dumdz = 2.0e0_rt*dum*sf.zbari; - } - - gm1 = 1.0e0_rt/dum; - gm2 = gm1*gm1; - gm13 = std::pow(gm1, nu_constants::oneth); - gm23 = gm13 * gm13; - gm43 = gm13*gm1; - gm53 = gm23*gm1; - - // equation 5.25 and 5.26 - v = -0.05483e0_rt - 0.01946e0_rt*gm13 + 1.86310e0_rt*gm23 - 0.78873e0_rt*gm1; - a0 = nu_constants::oneth*0.01946e0_rt*gm43 - nu_constants::twoth*1.86310e0_rt*gm53 + 0.78873e0_rt*gm2; - - w = -0.06711e0_rt + 0.06859e0_rt*gm13 + 1.74360e0_rt*gm23 - 0.74498e0_rt*gm1; - a1 = -nu_constants::oneth*0.06859e0_rt*gm43 - nu_constants::twoth*1.74360e0_rt*gm53 + 0.74498e0_rt*gm2; - - // equation 5.19 and 5.20 - fliq = v*fb + (1.0e0_rt - v)*ft; - if constexpr (do_derivatives) { - fliqdt = a0*dumdt*(fb - ft); - fliqda = a0*dumda*(fb - ft); - fliqdz = a0*dumdz*(fb - ft); - } - - gliq = w*gb + (1.0e0_rt - w)*gt; - if constexpr (do_derivatives) { - gliqdt = a1*dumdt*(gb - gt); - gliqda = a1*dumda*(gb - gt); - gliqdz = a1*dumdz*(gb - gt); - } - - // equation 5.17 - dum = 0.5738e0_rt*zbar*sf.ye*t86*den; - if constexpr (do_derivatives) { - dumdt = 0.5738e0_rt*zbar*sf.ye*6.0e0_rt*t85*den*1.0e-8_rt; - dumda = -dum*sf.abari; - dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; - } - - z = nu_constants::tfac4*fliq - nu_constants::tfac5*gliq; - sbrem = dum * z; - if constexpr (do_derivatives) { - sbremdt = dumdt*z + dum*(nu_constants::tfac4*fliqdt - nu_constants::tfac5*gliqdt); - sbremda = dumda*z + dum*(nu_constants::tfac4*fliqda - nu_constants::tfac5*gliqda); - sbremdz = dumdz*z + dum*(nu_constants::tfac4*fliqdz - nu_constants::tfac5*gliqdz); - } - } + nu_brem(sf, sbrem, sbremdt, sbremda, sbremdz); nu_recomb(sf, sreco, srecodt, srecoda, srecodz); From 5a705426198fb9a90edd555935136513c731ecfb Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 19:45:41 -0400 Subject: [PATCH 08/10] fix comp --- neutrinos/sneut5.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index d5ac326e93..b48cc1a491 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -802,10 +802,10 @@ void sneut5(const Real temp, const Real den, dsnudz = derivative of snu with zbar */ - Real xlnt,cc,den6,tfermi, + Real xlnt,cc, a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06, c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, - c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, + c25,c26,dd01,dd02,dd03,dd04,dd05,dd11,dd12, dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,c,d{0.0},f0, f1{0.0},f2,z, dum,dumdt,gum,dumda,dumdz; From c8890c9e7de85de8dc56aecac26ad8adc1347cd6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 12:22:16 -0400 Subject: [PATCH 09/10] move back --- neutrinos/sneut5.H | 186 ++++++++++++++++++++++----------------------- 1 file changed, 90 insertions(+), 96 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 3225cc6822..2280ca6a1f 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -4,24 +4,44 @@ #include #include -using namespace amrex; -namespace ifermi12_coeffs { - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a1 = { +AMREX_GPU_HOST_DEVICE AMREX_INLINE +Real ifermi12(const Real f) +{ + + // this routine applies a rational function expansion to get the inverse + // fermi-dirac integral of order 1/2 when it is equal to f. + // maximum error is 4.19e-9_rt. reference: antia apjs 84,101 1993 + + // declare work variables + Real rn,den,ff; + + // the return value + Real ifermi12r; + + // load the coefficients of the expansion from Table 8 of Antia + + constexpr Real an{0.5e0_rt}; + constexpr int m1{4}; + constexpr int k1{3}; + constexpr int m2{6}; + constexpr int k2{5}; + + const Array1D a1 = { 1.999266880833e4_rt, 5.702479099336e3_rt, 6.610132843877e2_rt, 3.818838129486e1_rt, 1.0e0_rt}; - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b1 = { + const Array1D b1 = { 1.771804140488e4_rt, -2.014785161019e3_rt, 9.130355392717e1_rt, -1.670718177489e0_rt}; - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a2 = { + const Array1D a2 = { -1.277060388085e-2_rt, 7.187946804945e-2_rt, -4.262314235106e-1_rt, @@ -30,87 +50,13 @@ namespace ifermi12_coeffs { -3.930805454272e-1_rt, 1.0e0_rt}; - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b2 = { + const Array1D b2 = { -9.745794806288e-3_rt, 5.485432756838e-2_rt, -3.299466243260e-1_rt, 4.077841975923e-1_rt, -1.145531476975e0_rt, -6.067091689181e-2_rt}; -}; - -namespace zfermim12_coeffs { - - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a1 = { - 1.71446374704454e7_rt, - 3.88148302324068e7_rt, - 3.16743385304962e7_rt, - 1.14587609192151e7_rt, - 1.83696370756153e6_rt, - 1.14980998186874e5_rt, - 1.98276889924768e3_rt, - 1.0e0_rt}; - - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b1 = { - 9.67282587452899e6_rt, - 2.87386436731785e7_rt, - 3.26070130734158e7_rt, - 1.77657027846367e7_rt, - 4.81648022267831e6_rt, - 6.13709569333207e5_rt, - 3.13595854332114e4_rt, - 4.35061725080755e2_rt}; - - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a2 = { - -4.46620341924942e-15_rt, - -1.58654991146236e-12_rt, - -4.44467627042232e-10_rt, - -6.84738791621745e-8_rt, - -6.64932238528105e-6_rt, - -3.69976170193942e-4_rt, - -1.12295393687006e-2_rt, - -1.60926102124442e-1_rt, - -8.52408612877447e-1_rt, - -7.45519953763928e-1_rt, - 2.98435207466372e0_rt, - 1.0e0_rt}; - - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D b2 = { - -2.23310170962369e-15_rt, - -7.94193282071464e-13_rt, - -2.22564376956228e-10_rt, - -3.43299431079845e-8_rt, - -3.33919612678907e-6_rt, - -1.86432212187088e-4_rt, - -5.69764436880529e-3_rt, - -8.34904593067194e-2_rt, - -4.78770844009440e-1_rt, - -4.99759250374148e-1_rt, - 1.86795964993052e0_rt, - 4.16485970495288e-1_rt}; -}; - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -Real ifermi12(const Real f) -{ - - // this routine applies a rational function expansion to get the inverse - // fermi-dirac integral of order 1/2 when it is equal to f. - // maximum error is 4.19e-9_rt. reference: antia apjs 84,101 1993 - - // declare work variables - Real rn,den,ff; - - // the return value - Real ifermi12r; - - // load the coefficients of the expansion from Table 8 of Antia - - constexpr Real an{0.5e0_rt}; - constexpr int m1{4}; - constexpr int k1{3}; - constexpr int m2{6}; - constexpr int k2{5}; if (f < 4.0e0_rt) { @@ -125,18 +71,18 @@ Real ifermi12(const Real f) // // in the starting rn here, the leading f is actually // a1(m1+1) * f, but that element is 1 - rn = f + ifermi12_coeffs::a1(m1); + rn = f + a1(m1); for (int i = m1 - 1; i >= 1; --i) { - rn = rn*f + ifermi12_coeffs::a1(i); + rn = rn*f + a1(i); } // now we do the denominator in Eq. 4. None of the coefficients // are 1, so we loop over all - den = ifermi12_coeffs::b1(k1+1); + den = b1(k1+1); for (int i = k1; i >= 1; --i) { - den = den*f + ifermi12_coeffs::b1(i); + den = den*f + b1(i); } // Eq. 6 of Antia @@ -150,16 +96,16 @@ Real ifermi12(const Real f) // this construction is the same as above, but using the // second set of coefficients - rn = ff + ifermi12_coeffs::a2(m2); + rn = ff + a2(m2); for (int i = m2 - 1; i >= 1; --i) { - rn = rn*ff + ifermi12_coeffs::a2(i); + rn = rn*ff + a2(i); } - den = ifermi12_coeffs::b2(k2+1); + den = b2(k2+1); for (int i = k2; i >= 1; --i) { - den = den*ff + ifermi12_coeffs::b2(i); + den = den*ff + b2(i); } ifermi12r = rn/(den*ff); @@ -190,19 +136,67 @@ Real zfermim12(const Real x) constexpr int m2{11}; constexpr int k2{11}; + const Array1D a1 = { + 1.71446374704454e7_rt, + 3.88148302324068e7_rt, + 3.16743385304962e7_rt, + 1.14587609192151e7_rt, + 1.83696370756153e6_rt, + 1.14980998186874e5_rt, + 1.98276889924768e3_rt, + 1.0e0_rt}; + + const Array1D b1 = { + 9.67282587452899e6_rt, + 2.87386436731785e7_rt, + 3.26070130734158e7_rt, + 1.77657027846367e7_rt, + 4.81648022267831e6_rt, + 6.13709569333207e5_rt, + 3.13595854332114e4_rt, + 4.35061725080755e2_rt}; + + const Array1D a2 = { + -4.46620341924942e-15_rt, + -1.58654991146236e-12_rt, + -4.44467627042232e-10_rt, + -6.84738791621745e-8_rt, + -6.64932238528105e-6_rt, + -3.69976170193942e-4_rt, + -1.12295393687006e-2_rt, + -1.60926102124442e-1_rt, + -8.52408612877447e-1_rt, + -7.45519953763928e-1_rt, + 2.98435207466372e0_rt, + 1.0e0_rt}; + + const Array1D b2 = { + -2.23310170962369e-15_rt, + -7.94193282071464e-13_rt, + -2.22564376956228e-10_rt, + -3.43299431079845e-8_rt, + -3.33919612678907e-6_rt, + -1.86432212187088e-4_rt, + -5.69764436880529e-3_rt, + -8.34904593067194e-2_rt, + -4.78770844009440e-1_rt, + -4.99759250374148e-1_rt, + 1.86795964993052e0_rt, + 4.16485970495288e-1_rt}; + if (x < 2.0e0_rt) { xx = std::exp(x); - rn = xx + zfermim12_coeffs::a1(m1); + rn = xx + a1(m1); for (int i = m1 - 1; i >= 1; --i) { - rn = rn*xx + zfermim12_coeffs::a1(i); + rn = rn*xx + a1(i); } - den = zfermim12_coeffs::b1(k1+1); + den = b1(k1+1); for (int i = k1; i >= 1; --i) { - den = den*xx + zfermim12_coeffs::b1(i); + den = den*xx + b1(i); } zfermim12r = xx * rn/den; @@ -210,16 +204,16 @@ Real zfermim12(const Real x) } else { xx = 1.0e0_rt/(x*x); - rn = xx + zfermim12_coeffs::a2(m2); + rn = xx + a2(m2); for (int i = m2 - 1; i >= 1; --i) { - rn = rn*xx + zfermim12_coeffs::a2(i); + rn = rn*xx + a2(i); } - den = zfermim12_coeffs::b2(k2+1); + den = b2(k2+1); for (int i = k2; i >= 1; --i) { - den = den*xx + zfermim12_coeffs::b2(i); + den = den*xx + b2(i); } zfermim12r = std::sqrt(x)*rn/den; From 5edd1092c8a3918d6a88fc524f1ebed55c3e61f3 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 12:23:10 -0400 Subject: [PATCH 10/10] put back namespace --- neutrinos/sneut5.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 2280ca6a1f..b5af627dcd 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -4,7 +4,7 @@ #include #include - +using namespace amrex; AMREX_GPU_HOST_DEVICE AMREX_INLINE Real ifermi12(const Real f)