From 37a921e7b846afa8ce6d783ca5010ecc1e76a862 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 14:08:22 -0400 Subject: [PATCH 01/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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 f719ff8bed9fdc7588f01f71850b8fa6b9ebd870 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 20:32:52 -0400 Subject: [PATCH 09/28] no do photo --- neutrinos/sneut5.H | 681 ++++++++++++++++++++++++--------------------- 1 file changed, 357 insertions(+), 324 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index b48cc1a491..d6aa409b78 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -297,9 +297,16 @@ struct sneutf_t { Real rmdz; Real rmi; + Real zeta; + Real zeta2; + Real zeta3; + Real zetadt; + Real zetada; + Real zetadz; }; +template AMREX_GPU_HOST_DEVICE inline sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { @@ -345,10 +352,342 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { sf.rmdz = den * sf.abari; sf.rmi = 1.0e0_rt / sf.rm; + Real a0 = sf.rm * 1.0e-9_rt; + Real a1 = std::pow(a0, nu_constants::oneth); + sf.zeta = a1 * sf.xlm1; + + if constexpr (do_derivatives) { + Real a2 = nu_constants::oneth * a1 * sf.rmi * sf.xlm1; + sf.zetadt = -a1 * sf.xlm2 * sf.xldt; + sf.zetada = a2 * sf.rmda; + sf.zetadz = a2 * sf.rmdz; + } + + sf.zeta2 = sf.zeta * sf.zeta; + sf.zeta3 = sf.zeta2 * sf.zeta; + return sf; } +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_photo(const sneutf_t& sf, + Real& sphot, Real& sphotdt, Real& sphotda, Real& sphotdz) { + + // photoneutrino process section + // for reactions like e- + gamma => e- + nu_e + nubar_e + // e+ + gamma => e+ + nu_e + nubar_e + // equation 3.8 for tau, equation 3.6 for cc, + // and table 2 written out for speed + + Real tau; + Real cc; + Real c00, c01, c02, c03, c04, c05, c06; + Real c10, c11, c12, c13, c14, c15, c16; + Real c20, c21, c22, c23, c24, c25, c26; + Real dd01, dd02, dd03, dd04, dd05; + Real dd11, dd12, dd13, dd14, dd15; + Real dd21, dd22, dd23, dd24, dd25; + + if (sf.temp < 1.0e8_rt) { + + // note: we already bailed above for T < 1.e7, so this is really 1.e7 <= T < 1.e8 + + tau = std::log10(sf.temp * 1.0e-7_rt); + cc = 0.5654e0_rt + tau; + c00 = 1.008e11_rt; + c01 = 0.0e0_rt; + c02 = 0.0e0_rt; + c03 = 0.0e0_rt; + c04 = 0.0e0_rt; + c05 = 0.0e0_rt; + c06 = 0.0e0_rt; + c10 = 8.156e10_rt; + c11 = 9.728e8_rt; + c12 = -3.806e9_rt; + c13 = -4.384e9_rt; + c14 = -5.774e9_rt; + c15 = -5.249e9_rt; + c16 = -5.153e9_rt; + c20 = 1.067e11_rt; + c21 = -9.782e9_rt; + c22 = -7.193e9_rt; + c23 = -6.936e9_rt; + c24 = -6.893e9_rt; + c25 = -7.041e9_rt; + c26 = -7.193e9_rt; + dd01 = 0.0e0_rt; + dd02 = 0.0e0_rt; + dd03 = 0.0e0_rt; + dd04 = 0.0e0_rt; + dd05 = 0.0e0_rt; + dd11 = -1.879e10_rt; + dd12 = -9.667e9_rt; + dd13 = -5.602e9_rt; + dd14 = -3.370e9_rt; + dd15 = -1.825e9_rt; + dd21 = -2.919e10_rt; + dd22 = -1.185e10_rt; + dd23 = -7.270e9_rt; + dd24 = -4.222e9_rt; + dd25 = -1.560e9_rt; + + } else if (sf.temp >= 1.0e8_rt && sf.temp < 1.0e9_rt) { + + tau = std::log10(sf.temp * 1.0e-8_rt); + cc = 1.5654e0_rt; + c00 = 9.889e10_rt; + c01 = -4.524e8_rt; + c02 = -6.088e6_rt; + c03 = 4.269e7_rt; + c04 = 5.172e7_rt; + c05 = 4.910e7_rt; + c06 = 4.388e7_rt; + c10 = 1.813e11_rt; + c11 = -7.556e9_rt; + c12 = -3.304e9_rt; + c13 = -1.031e9_rt; + c14 = -1.764e9_rt; + c15 = -1.851e9_rt; + c16 = -1.928e9_rt; + c20 = 9.750e10_rt; + c21 = 3.484e10_rt; + c22 = 5.199e9_rt; + c23 = -1.695e9_rt; + c24 = -2.865e9_rt; + c25 = -3.395e9_rt; + c26 = -3.418e9_rt; + dd01 = -1.135e8_rt; + dd02 = 1.256e8_rt; + dd03 = 5.149e7_rt; + dd04 = 3.436e7_rt; + dd05 = 1.005e7_rt; + dd11 = 1.652e9_rt; + dd12 = -3.119e9_rt; + dd13 = -1.839e9_rt; + dd14 = -1.458e9_rt; + dd15 = -8.956e8_rt; + dd21 = -1.548e10_rt; + dd22 = -9.338e9_rt; + dd23 = -5.899e9_rt; + dd24 = -3.035e9_rt; + dd25 = -1.598e9_rt; + + } else { + + // T > 1.e9 + + tau = std::log10(sf.t9); + cc = 1.5654e0_rt; + c00 = 9.581e10_rt; + c01 = 4.107e8_rt; + c02 = 2.305e8_rt; + c03 = 2.236e8_rt; + c04 = 1.580e8_rt; + c05 = 2.165e8_rt; + c06 = 1.721e8_rt; + c10 = 1.459e12_rt; + c11 = 1.314e11_rt; + c12 = -1.169e11_rt; + c13 = -1.765e11_rt; + c14 = -1.867e11_rt; + c15 = -1.983e11_rt; + c16 = -1.896e11_rt; + c20 = 2.424e11_rt; + c21 = -3.669e9_rt; + c22 = -8.691e9_rt; + c23 = -7.967e9_rt; + c24 = -7.932e9_rt; + c25 = -7.987e9_rt; + c26 = -8.333e9_rt; + dd01 = 4.724e8_rt; + dd02 = 2.976e8_rt; + dd03 = 2.242e8_rt; + dd04 = 7.937e7_rt; + dd05 = 4.859e7_rt; + dd11 = -7.094e11_rt; + dd12 = -3.697e11_rt; + dd13 = -2.189e11_rt; + dd14 = -1.273e11_rt; + dd15 = -5.705e10_rt; + dd21 = -2.254e10_rt; + dd22 = -1.551e10_rt; + dd23 = -7.793e9_rt; + dd24 = -4.489e9_rt; + dd25 = -2.185e9_rt; + + } + + Real taudt = nu_constants::iln10 * sf.tempi; + + // equation 3.7, compute the expensive trig functions only one time + + Real cos1 = std::cos(nu_constants::fac1 * tau); + Real cos2 = std::cos(nu_constants::fac1 * 2.0e0_rt * tau); + Real cos3 = std::cos(nu_constants::fac1 * 3.0e0_rt * tau); + Real cos4 = std::cos(nu_constants::fac1 * 4.0e0_rt * tau); + Real cos5 = std::cos(nu_constants::fac1 * 5.0e0_rt * tau); + Real last = std::cos(nu_constants::fac2 * tau); + + Real sin1 = std::sin(nu_constants::fac1*tau); + Real sin2 = std::sin(nu_constants::fac1*2.0e0_rt*tau); + Real sin3 = std::sin(nu_constants::fac1*3.0e0_rt*tau); + Real sin4 = std::sin(nu_constants::fac1*4.0e0_rt*tau); + Real sin5 = std::sin(nu_constants::fac1*5.0e0_rt*tau); + Real xast = std::sin(nu_constants::fac2*tau); + + Real a0 = 0.5e0_rt * c00 + + c01 * cos1 + dd01 * sin1 + c02 * cos2 + dd02 * sin2 + + c03 * cos3 + dd03 * sin3 + c04 * cos4 + dd04 * sin4 + + c05 * cos5 + dd05 * sin5 + 0.5e0_rt * c06 * last; + Real a1 = 0.5e0_rt * c10 + + c11 * cos1 + dd11 * sin1 + c12 * cos2 + dd12 * sin2 + + c13 * cos3 + dd13 * sin3 + c14 * cos4 + dd14 * sin4 + + c15 * cos5 + dd15 * sin5 + 0.5e0_rt * c16 * last; + Real a2 = 0.5e0_rt * c20 + + c21 * cos1 + dd21 * sin1 + c22 * cos2 + dd22 * sin2 + + c23 * cos3 + dd23 * sin3 + c24 * cos4 + dd24 * sin4 + + c25 * cos5 + dd25 * sin5 + 0.5e0_rt * c26 * last; + + Real f0, f1, f2; + if constexpr (do_derivatives) { + 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 * nu_constants::fac2 * taudt; + + 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*nu_constants::fac2*taudt; + + 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 * nu_constants::fac2 * taudt; + } + + // equation 3.4 + Real dum = a0 + a1 * sf.zeta + a2 * sf.zeta2; + Real dumdt, dumda, dumdz; + if constexpr (do_derivatives) { + dumdt = f0 + f1 * sf.zeta + a1 * sf.zetadt + f2 * sf.zeta2 + 2.0e0_rt * a2 * sf.zeta * sf.zetadt; + dumda = a1 * sf.zetada + 2.0e0_rt * a2 * sf.zeta * sf.zetada; + dumdz = a1 * sf.zetadz + 2.0e0_rt * a2 * sf.zeta * sf.zetadz; + } + + Real z = std::exp(-cc * sf.zeta); + + Real xnum = dum * z; + Real xnumdt, xnumda, xnumdz; + if constexpr (do_derivatives) { + xnumdt = dumdt * z - dum * z * cc * sf.zetadt; + xnumda = dumda * z - dum * z * cc * sf.zetada; + xnumdz = dumdz * z - dum * z * cc * sf.zetadz; + } + + Real xden = sf.zeta3 + 6.290e-3_rt * sf.xlm1 + 7.483e-3_rt * sf.xlm2 + 3.061e-4_rt * sf.xlm3; + + dum = 3.0e0_rt * sf.zeta2; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + xdendt = dum * sf.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 * sf.zetada; + xdendz = dum * sf.zetadz; + } + dum = 1.0e0_rt / xden; + Real fphot = xnum * dum; + Real fphotdt, fphotda, fphotdz; + if constexpr (do_derivatives) { + fphotdt = (xnumdt - fphot * xdendt) * dum; + fphotda = (xnumda - fphot * xdenda) * dum; + fphotdz = (xnumdz - fphot * xdendz) * dum; + } + + // equation 3.3 + 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 * sf.xldt; + } + + 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 = 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 + sf.rm * z; + if constexpr (do_derivatives) { + xdendt = -sf.rm * z * z * dumdt; + xdenda = sf.rmda * z; + xdendz = sf.rmdz * z; + } + + z = 1.0e0_rt / xden; + Real qphot = xnum * z; + Real qphotdt; + if constexpr (do_derivatives) { + qphotdt = (xnumdt - qphot * xdendt) * z; + } + dum = -qphot * z; + Real qphotda, qphotdz; + if constexpr (do_derivatives) { + qphotda = dum * xdenda; + qphotdz = dum * xdendz; + } + + // equation 3.2 + sphot = sf.xl5 * fphot; + if constexpr (do_derivatives) { + sphotdt = 5.0e0_rt * sf.xl4 * sf.xldt * fphot + sf.xl5 * fphotdt; + sphotda = sf.xl5 * fphotda; + sphotdz = sf.xl5 * fphotdz; + } + + a1 = sphot; + sphot = sf.rm * a1; + if constexpr (do_derivatives) { + sphotdt = sf.rm * sphotdt; + sphotda = sf.rm * sphotda + sf.rmda * a1; + sphotdz = sf.rm * sphotdz + sf.rmdz * a1; + } + + a1 = nu_constants::tfac4 * (1.0e0_rt - nu_constants::tfac3 * qphot); + Real a3 = sphot; + sphot = a1 * a3; + if constexpr (do_derivatives) { + 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; + } + + if (sphot <= 0.0_rt) { + sphot = 0.0e0_rt; + if constexpr (do_derivatives) { + sphotdt = 0.0e0_rt; + sphotda = 0.0e0_rt; + sphotdz = 0.0e0_rt; + } + } +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void nu_brem(const sneutf_t& sf, @@ -762,16 +1101,16 @@ void nu_recomb(const sneutf_t& sf, } // equation 6.5 - z = std::exp(zeta + nu); + z = std::exp(sf.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); + srecodt = sreco * (bigjdt * a2 - z * (sf.zetadt + nudt) * a1); + srecoda = sreco * (-1.0e0_rt * sf.abari + bigjda * a2 - z * (sf.zetada + nuda) * a1); + srecodz = sreco * (14.0e0_rt * sf.zbari + bigjdz * a2 - z * (sf.zetadz + nudz) * a1); } } @@ -812,11 +1151,9 @@ void sneut5(const Real temp, const Real den, // pair production Real gl,gldt, - zeta,zetadt,zeta2,zeta3, xnum,xnumdt,xden,xdendt, fpair,fpairdt,qpair,qpairdt, fpairda,fpairdz,qpairda,qpairdz, - zetada,zetadz, xnumda,xnumdz,xdenda,xdendz; // plasma @@ -825,11 +1162,6 @@ void sneut5(const Real temp, const Real den, flda,fldz,ftda,ftdz, fxyda,fxydz,gl2da,gl2dz; - // photo - Real tau,taudt,cos1,cos2,cos3,cos4,cos5,sin1,sin2, - sin3,sin4,sin5,last,xast, - fphot,fphotdt,qphot,qphotdt, - fphotda,fphotdz,qphotda,qphotdz; // initialize Real spair{0.0e0_rt}; @@ -865,21 +1197,7 @@ void sneut5(const Real temp, const Real den, if (temp < 1.0e7_rt) return; - auto sf = get_sneut_factors(den, temp, abar, zbar); - - a0 = sf.rm * 1.0e-9_rt; - a1 = std::pow(a0, nu_constants::oneth); - zeta = a1 * sf.xlm1; - - if constexpr (do_derivatives) { - a2 = nu_constants::oneth * a1*sf.rmi * sf.xlm1; - zetadt = -a1 * sf.xlm2 * sf.xldt; - zetada = a2 * sf.rmda; - zetadz = a2 * sf.rmdz; - } - - zeta2 = zeta * zeta; - zeta3 = zeta2 * zeta; + auto sf = get_sneut_factors(den, temp, abar, zbar); // pair neutrino section // for reactions like e+ + e- => nu_e + nubar_e @@ -892,15 +1210,15 @@ void sneut5(const Real temp, const Real den, // equation 2.7 - a1 = 6.002e19_rt + 2.084e20_rt*zeta + 1.872e21_rt*zeta2; + a1 = 6.002e19_rt + 2.084e20_rt * sf.zeta + 1.872e21_rt * sf.zeta2; if (sf.t9 < 10.0_rt) { - b1 = std::exp(-5.5924e0_rt*zeta); + b1 = std::exp(-5.5924e0_rt * sf.zeta); if constexpr (do_derivatives) { b2 = -b1*5.5924e0_rt; } } else { - b1 = std::exp(-4.9924e0_rt*zeta); + b1 = std::exp(-4.9924e0_rt * sf.zeta); if constexpr (do_derivatives) { b2 = -b1*4.9924e0_rt; } @@ -908,11 +1226,11 @@ void sneut5(const Real temp, const Real den, xnum = a1 * b1; if constexpr (do_derivatives) { - a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta; + a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt * sf.zeta; c = a2*b1 + a1*b2; - xnumdt = c*zetadt; - xnumda = c*zetada; - xnumdz = c*zetadz; + xnumdt = c * sf.zetadt; + xnumda = c * sf.zetada; + xnumdz = c * sf.zetadz; } if (sf.t9 < 10.0_rt) { @@ -927,12 +1245,12 @@ void sneut5(const Real temp, const Real den, } } - xden = zeta3 + a1; + xden = sf.zeta3 + a1; if constexpr (do_derivatives) { - b1 = 3.0e0_rt*zeta2; - xdendt = b1*zetadt + a2*sf.xldt; - xdenda = b1*zetada; - xdendz = b1*zetadz; + b1 = 3.0e0_rt*sf.zeta2; + xdendt = b1*sf.zetadt + a2*sf.xldt; + xdenda = b1*sf.zetada; + xdendz = b1*sf.zetadz; } a1 = 1.0e0_rt/xden; @@ -1161,292 +1479,7 @@ void sneut5(const Real temp, const Real den, splasdz = a2*splasdz; } - // photoneutrino process section - // for reactions like e- + gamma => e- + nu_e + nubar_e - // e+ + gamma => e+ + nu_e + nubar_e - // equation 3.8 for tau, equation 3.6 for cc, - // and table 2 written out for speed - if (temp < 1.0e8_rt) { - - // note: we already bailed above for T < 1.e7, so this is really 1.e7 <= T < 1.e8 - - tau = std::log10(temp * 1.0e-7_rt); - cc = 0.5654e0_rt + tau; - c00 = 1.008e11_rt; - c01 = 0.0e0_rt; - c02 = 0.0e0_rt; - c03 = 0.0e0_rt; - c04 = 0.0e0_rt; - c05 = 0.0e0_rt; - c06 = 0.0e0_rt; - c10 = 8.156e10_rt; - c11 = 9.728e8_rt; - c12 = -3.806e9_rt; - c13 = -4.384e9_rt; - c14 = -5.774e9_rt; - c15 = -5.249e9_rt; - c16 = -5.153e9_rt; - c20 = 1.067e11_rt; - c21 = -9.782e9_rt; - c22 = -7.193e9_rt; - c23 = -6.936e9_rt; - c24 = -6.893e9_rt; - c25 = -7.041e9_rt; - c26 = -7.193e9_rt; - dd01 = 0.0e0_rt; - dd02 = 0.0e0_rt; - dd03 = 0.0e0_rt; - dd04 = 0.0e0_rt; - dd05 = 0.0e0_rt; - dd11 = -1.879e10_rt; - dd12 = -9.667e9_rt; - dd13 = -5.602e9_rt; - dd14 = -3.370e9_rt; - dd15 = -1.825e9_rt; - dd21 = -2.919e10_rt; - dd22 = -1.185e10_rt; - dd23 = -7.270e9_rt; - dd24 = -4.222e9_rt; - dd25 = -1.560e9_rt; - - } else if (temp >= 1.0e8_rt && temp < 1.0e9_rt) { - - tau = std::log10(temp * 1.0e-8_rt); - cc = 1.5654e0_rt; - c00 = 9.889e10_rt; - c01 = -4.524e8_rt; - c02 = -6.088e6_rt; - c03 = 4.269e7_rt; - c04 = 5.172e7_rt; - c05 = 4.910e7_rt; - c06 = 4.388e7_rt; - c10 = 1.813e11_rt; - c11 = -7.556e9_rt; - c12 = -3.304e9_rt; - c13 = -1.031e9_rt; - c14 = -1.764e9_rt; - c15 = -1.851e9_rt; - c16 = -1.928e9_rt; - c20 = 9.750e10_rt; - c21 = 3.484e10_rt; - c22 = 5.199e9_rt; - c23 = -1.695e9_rt; - c24 = -2.865e9_rt; - c25 = -3.395e9_rt; - c26 = -3.418e9_rt; - dd01 = -1.135e8_rt; - dd02 = 1.256e8_rt; - dd03 = 5.149e7_rt; - dd04 = 3.436e7_rt; - dd05 = 1.005e7_rt; - dd11 = 1.652e9_rt; - dd12 = -3.119e9_rt; - dd13 = -1.839e9_rt; - dd14 = -1.458e9_rt; - dd15 = -8.956e8_rt; - dd21 = -1.548e10_rt; - dd22 = -9.338e9_rt; - dd23 = -5.899e9_rt; - dd24 = -3.035e9_rt; - dd25 = -1.598e9_rt; - - } else { - - // T > 1.e9 - - tau = std::log10(sf.t9); - cc = 1.5654e0_rt; - c00 = 9.581e10_rt; - c01 = 4.107e8_rt; - c02 = 2.305e8_rt; - c03 = 2.236e8_rt; - c04 = 1.580e8_rt; - c05 = 2.165e8_rt; - c06 = 1.721e8_rt; - c10 = 1.459e12_rt; - c11 = 1.314e11_rt; - c12 = -1.169e11_rt; - c13 = -1.765e11_rt; - c14 = -1.867e11_rt; - c15 = -1.983e11_rt; - c16 = -1.896e11_rt; - c20 = 2.424e11_rt; - c21 = -3.669e9_rt; - c22 = -8.691e9_rt; - c23 = -7.967e9_rt; - c24 = -7.932e9_rt; - c25 = -7.987e9_rt; - c26 = -8.333e9_rt; - dd01 = 4.724e8_rt; - dd02 = 2.976e8_rt; - dd03 = 2.242e8_rt; - dd04 = 7.937e7_rt; - dd05 = 4.859e7_rt; - dd11 = -7.094e11_rt; - dd12 = -3.697e11_rt; - dd13 = -2.189e11_rt; - dd14 = -1.273e11_rt; - dd15 = -5.705e10_rt; - dd21 = -2.254e10_rt; - dd22 = -1.551e10_rt; - dd23 = -7.793e9_rt; - dd24 = -4.489e9_rt; - dd25 = -2.185e9_rt; - - } - - taudt = nu_constants::iln10*sf.tempi; - - // equation 3.7, compute the expensive trig functions only one time - 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 - + c03*cos3 + dd03*sin3 + c04*cos4 + dd04*sin4 - + c05*cos5 + dd05*sin5 + 0.5e0_rt*c06*last; - a1 = 0.5e0_rt*c10 - + c11*cos1 + dd11*sin1 + c12*cos2 + dd12*sin2 - + c13*cos3 + dd13*sin3 + c14*cos4 + dd14*sin4 - + c15*cos5 + dd15*sin5 + 0.5e0_rt*c16*last; - - - a2 = 0.5e0_rt*c20 - + c21*cos1 + dd21*sin1 + c22*cos2 + dd22*sin2 - + c23*cos3 + dd23*sin3 + c24*cos4 + dd24*sin4 - + c25*cos5 + dd25*sin5 + 0.5e0_rt*c26*last; - - if constexpr (do_derivatives) { - 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*nu_constants::fac2*taudt; - - 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*nu_constants::fac2*taudt; - - 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*nu_constants::fac2*taudt; - } - - - // equation 3.4 - dum = a0 + a1*zeta + a2*zeta2; - if constexpr (do_derivatives) { - dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0e0_rt*a2*zeta*zetadt; - dumda = a1*zetada + 2.0e0_rt*a2*zeta*zetada; - dumdz = a1*zetadz + 2.0e0_rt*a2*zeta*zetadz; - } - - z = std::exp(-cc*zeta); - - xnum = dum*z; - if constexpr (do_derivatives) { - xnumdt = dumdt*z - dum*z*cc*zetadt; - xnumda = dumda*z - dum*z*cc*zetada; - xnumdz = dumdz*z - dum*z*cc*zetadz; - } - - 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 - 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; - } - dum = 1.0e0_rt/xden; - fphot = xnum*dum; - if constexpr (do_derivatives) { - fphotdt = (xnumdt - fphot*xdendt)*dum; - fphotda = (xnumda - fphot*xdenda)*dum; - fphotdz = (xnumdz - fphot*xdendz)*dum; - } - - // equation 3.3 - 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*sf.xldt; - } - - 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 = 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 + sf.rm*z; - if constexpr (do_derivatives) { - xdendt = -sf.rm*z*z*dumdt; - xdenda = sf.rmda*z; - xdendz = sf.rmdz*z; - } - - z = 1.0e0_rt/xden; - qphot = xnum*z; - if constexpr (do_derivatives) { - qphotdt = (xnumdt - qphot*xdendt)*z; - } - dum = -qphot*z; - if constexpr (do_derivatives) { - qphotda = dum*xdenda; - qphotdz = dum*xdendz; - } - - // equation 3.2 - sphot = sf.xl5 * fphot; - if constexpr (do_derivatives) { - sphotdt = 5.0e0_rt*sf.xl4*sf.xldt*fphot + sf.xl5*fphotdt; - sphotda = sf.xl5*fphotda; - sphotdz = sf.xl5*fphotdz; - } - - a1 = sphot; - sphot = sf.rm*a1; - if constexpr (do_derivatives) { - sphotdt = sf.rm*sphotdt; - sphotda = sf.rm*sphotda + sf.rmda*a1; - sphotdz = sf.rm*sphotdz + sf.rmdz*a1; - } - - a1 = nu_constants::tfac4*(1.0e0_rt - nu_constants::tfac3 * qphot); - - a3 = sphot; - sphot = a1*a3; - if constexpr (do_derivatives) { - 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; - } - - if (sphot <= 0.0_rt) { - sphot = 0.0e0_rt; - if constexpr (do_derivatives) { - sphotdt = 0.0e0_rt; - sphotda = 0.0e0_rt; - sphotdz = 0.0e0_rt; - } - } + nu_photo(sf, sphot, sphotdt, sphotda, sphotdz); nu_brem(sf, sbrem, sbremdt, sbremda, sbremdz); From aa0a60874e03b083505b44d1a6d93042abcf9676 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 Oct 2023 20:37:04 -0400 Subject: [PATCH 10/28] fix compt --- neutrinos/sneut5.H | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index d6aa409b78..6c6fdf1e66 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -1142,12 +1142,10 @@ void sneut5(const Real temp, const Real den, */ 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,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; + a1,a2,a3,b1,b2,c00,c01,c03,c04, + c,d{0.0}, + f1{0.0}, + dumdt,gum,dumda,dumdz; // pair production Real gl,gldt, From c8890c9e7de85de8dc56aecac26ad8adc1347cd6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 12:22:16 -0400 Subject: [PATCH 11/28] 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 12/28] 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) From 75fe182ab4b2ecc109844772c88fed3efe4e3072 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 15:11:53 -0400 Subject: [PATCH 13/28] fix --- neutrinos/sneut5.H | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index d6780dd7ed..38a66f29bd 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -1095,16 +1095,16 @@ void nu_recomb(const sneutf_t& sf, } // equation 6.5 - z = std::exp(sf.zeta + nu); + 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 * (sf.zetadt + nudt) * a1); - srecoda = sreco * (-1.0e0_rt * sf.abari + bigjda * a2 - z * (sf.zetada + nuda) * a1); - srecodz = sreco * (14.0e0_rt * sf.zbari + bigjdz * a2 - z * (sf.zetadz + nudz) * a1); + 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); } } From 6dd58008a573fd8b5a427d83e77310f07ce5feb0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 16:38:10 -0400 Subject: [PATCH 14/28] move plasma neutrino contribution into its own function --- neutrinos/sneut5.H | 333 ++++++++++++++++++++++++--------------------- 1 file changed, 175 insertions(+), 158 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 38a66f29bd..0717160668 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -364,6 +364,180 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_plasma(const sneutf_t& sf, + Real& splas, Real& splasdt, Real& splasda, Real& splasdz) { + + // plasma neutrino section + // for collective reactions like gamma_plasmon => nu_e + nubar_e + // equation 4.6 + + Real a1 = 1.019e-6_rt * sf.rm; + Real a2 = std::pow(a1, nu_constants::twoth); + + Real b1 = std::sqrt(1.0e0_rt + a2); + + Real c00 = 1.0e0_rt / (sf.temp * sf.temp * b1); + + Real gl2 = 1.1095e11_rt * sf.rm * c00; + Real gl2dt, gl2da, gl2dz; + if constexpr (do_derivatives) { + Real a3 = nu_constants::twoth * a2 / a1; + Real b2 = 1.0e0_rt / b1; + gl2dt = -2.0e0_rt * gl2 * sf.tempi; + Real 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); + } + + Real gl = std::sqrt(gl2); + Real gl12 = std::sqrt(gl); + Real gl32 = gl * gl12; + Real gl72 = gl2 * gl32; + Real gl6 = gl2 * gl2 * gl2; + + // equation 4.7 + Real ft = 2.4e0_rt + 0.6e0_rt * gl12 + 0.51e0_rt * gl + 1.25e0_rt * gl32; + Real gum = 1.0e0_rt / gl2; + Real ftdt, ftda, ftdz; + if constexpr (do_derivatives) { + a1 = (0.25e0_rt * 0.6e0_rt * gl12 + + 0.5e0_rt * 0.51e0_rt * gl + + 0.75e0_rt * 1.25e0_rt * gl32) * gum; + ftdt = a1 * gl2dt; + ftda = a1 * gl2da; + ftdz = a1 * gl2dz; + } + + // equation 4.8 + a1 = 8.6e0_rt * gl2 + 1.35e0_rt * gl72; + + b1 = 225.0e0_rt - 17.0e0_rt * gl + gl2; + + Real c = 1.0e0_rt / b1; + Real fl = a1 * c; + Real fldt, flda, fldz; + if constexpr (do_derivatives) { + a2 = 8.6e0_rt + 1.75e0_rt * 1.35e0_rt * gl72 * gum; + Real b2 = -0.5e0_rt * 17.0e0_rt * gl * gum + 1.0e0_rt; + Real d = (a2 - fl * b2) * c; + fldt = d * gl2dt; + flda = d * gl2da; + fldz = d * gl2dz; + } + + // equation 4.9 and 4.10 + Real cc = std::log10(2.0e0_rt * sf.rm); + Real xlnt = std::log10(sf.temp); + + Real xnum = nu_constants::sixth * (17.5e0_rt + cc - 3.0e0_rt * xlnt); + Real xden = nu_constants::sixth * (-24.5e0_rt + cc + 3.0e0_rt * xlnt); + Real xnumdt, xnumda, xnumdz; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + 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 = nu_constants::iln10 * 0.5e0_rt * sf.tempi; + xdenda = a2 * sf.rmda; + xdendz = a2 * sf.rmdz; + } + + // equation 4.11 + Real fxy, fxydt, fxyda, fxydz; + Real dumdt, dumda, dumdz; + if (std::abs(xnum) > 0.7e0_rt || xden < 0.0e0_rt) { + fxy = 1.0e0_rt; + fxydt = 0.0e0_rt; + fxydz = 0.0e0_rt; + fxyda = 0.0e0_rt; + + } else { + + a1 = 0.39e0_rt - 1.25e0_rt * xnum - 0.35e0_rt * std::sin(4.5e0_rt * xnum); + + b1 = 0.3e0_rt * std::exp(-std::pow(4.5e0_rt*xnum + 0.9e0_rt, 2.0_rt)); + + c = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt * xnum); + + Real b2; + if constexpr (do_derivatives) { + a2 = -1.25e0_rt - 4.5e0_rt * 0.35e0_rt * std::cos(4.5e0_rt * xnum); + b2 = -b1 * 2.0e0_rt * (4.5e0_rt * xnum + 0.9e0_rt) * 4.5e0_rt; + if (c == 0.0_rt) { + dumdt = 0.0e0_rt; + dumda = 0.0e0_rt; + dumdz = 0.0e0_rt; + } else { + dumdt = xdendt + 1.25e0_rt * xnumdt; + dumda = xdenda + 1.25e0_rt * xnumda; + dumdz = xdendz + 1.25e0_rt * xnumdz; + } + } + + Real d = 0.57e0_rt - 0.25e0_rt * xnum; + Real a3 = c / d; + c00 = std::exp(-a3 * a3); + + fxy = 1.05e0_rt + (a1 - b1) * c00; + if constexpr (do_derivatives) { + Real f1 = -c00 * 2.0e0_rt * a3 / d; + Real c01 = f1 * (dumdt + a3 * 0.25e0_rt * xnumdt); + Real c03 = f1 * (dumda + a3 * 0.25e0_rt * xnumda); + Real c04 = f1 * (dumdz + a3 * 0.25e0_rt * xnumdz); + + fxydt = (a2 * xnumdt - b2 * xnumdt) * c00 + (a1 - b1) * c01; + fxyda = (a2 * xnumda - b2 * xnumda) * c00 + (a1 - b1) * c03; + fxydz = (a2 * xnumdz - b2 * xnumdz) * c00 + (a1 - b1) * c04; + } + } + + // equation 4.1 and 4.5 + splas = (ft + fl) * fxy; + if constexpr (do_derivatives) { + splasdt = (ftdt + fldt) * fxy + (ft + fl) * fxydt; + splasda = (ftda + flda) * fxy + (ft + fl) * fxyda; + splasdz = (ftdz + fldz) * fxy + (ft + fl) * fxydz; + } + + a2 = std::exp(-gl); + + a1 = splas; + splas = a2*a1; + if constexpr (do_derivatives) { + Real a3 = -0.5e0_rt * a2 * gl * gum; + splasdt = a2 * splasdt + a3 * gl2dt * a1; + splasda = a2 * splasda + a3 * gl2da * a1; + splasdz = a2 * splasdz + a3 * gl2dz * a1; + } + + a2 = gl6; + + a1 = splas; + splas = a2 * a1; + if constexpr (do_derivatives) { + Real a3 = 3.0e0_rt * gl6 * gum; + splasdt = a2 * splasdt + a3 * gl2dt * a1; + splasda = a2 * splasda + a3 * gl2da * a1; + splasdz = a2 * splasdz + a3 * gl2dz * a1; + } + + a2 = 0.93153e0_rt * 3.0e21_rt * sf.xl9; + + a1 = splas; + splas = a2 * a1; + if constexpr (do_derivatives) { + Real 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; + } +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void nu_photo(const sneutf_t& sf, @@ -1311,165 +1485,8 @@ void sneut5(const Real temp, const Real den, spairdz = a1*spairdz + a2*qpairdz*a3; } - // plasma neutrino section - // for collective reactions like gamma_plasmon => nu_e + nubar_e - // equation 4.6 - - a1 = 1.019e-6_rt*sf.rm; - a2 = std::pow(a1, nu_constants::twoth); - - b1 = std::sqrt(1.0e0_rt + a2); - - c00 = 1.0e0_rt/(temp*temp*b1); - - gl2 = 1.1095e11_rt * sf.rm * c00; - - if constexpr (do_derivatives) { - 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; - gl2da = 1.1095e11_rt * (sf.rmda*c00 - d*sf.rmda); - gl2dz = 1.1095e11_rt * (sf.rmdz*c00 - d*sf.rmdz); - } - - gl = std::sqrt(gl2); - gl12 = std::sqrt(gl); - gl32 = gl * gl12; - gl72 = gl2 * gl32; - gl6 = gl2 * gl2 * gl2; - - // equation 4.7 - ft = 2.4e0_rt + 0.6e0_rt*gl12 + 0.51e0_rt*gl + 1.25e0_rt*gl32; - gum = 1.0e0_rt/gl2; - if constexpr (do_derivatives) { - a1 =(0.25e0_rt*0.6e0_rt*gl12 +0.5e0_rt*0.51e0_rt*gl +0.75e0_rt*1.25e0_rt*gl32)*gum; - ftdt = a1*gl2dt; - ftda = a1*gl2da; - ftdz = a1*gl2dz; - } - - // equation 4.8 - a1 = 8.6e0_rt*gl2 + 1.35e0_rt*gl72; - - b1 = 225.0e0_rt - 17.0e0_rt*gl + gl2; - - c = 1.0e0_rt/b1; - fl = a1*c; - - if constexpr (do_derivatives) { - a2 = 8.6e0_rt + 1.75e0_rt*1.35e0_rt*gl72*gum; - b2 = -0.5e0_rt*17.0e0_rt*gl*gum + 1.0e0_rt; - d = (a2 - fl*b2)*c; - fldt = d*gl2dt; - flda = d*gl2da; - fldz = d*gl2dz; - } - - // equation 4.9 and 4.10 - cc = std::log10(2.0e0_rt*sf.rm); - xlnt = std::log10(temp); - - 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 = -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 = nu_constants::iln10*0.5e0_rt*sf.tempi; - xdenda = a2*sf.rmda; - xdendz = a2*sf.rmdz; - } - - // equation 4.11 - if (std::abs(xnum) > 0.7e0_rt || xden < 0.0e0_rt) { - - fxy = 1.0e0_rt; - fxydt = 0.0e0_rt; - fxydz = 0.0e0_rt; - fxyda = 0.0e0_rt; - } else { - - a1 = 0.39e0_rt - 1.25e0_rt*xnum - 0.35e0_rt*std::sin(4.5e0_rt*xnum); - - b1 = 0.3e0_rt * std::exp(-std::pow(4.5e0_rt*xnum + 0.9e0_rt, 2.0_rt)); - - c = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt*xnum); - - if constexpr (do_derivatives) { - a2 = -1.25e0_rt - 4.5e0_rt*0.35e0_rt*std::cos(4.5e0_rt*xnum); - b2 = -b1*2.0e0_rt*(4.5e0_rt*xnum + 0.9e0_rt)*4.5e0_rt; - if (c == 0.0_rt) { - dumdt = 0.0e0_rt; - dumda = 0.0e0_rt; - dumdz = 0.0e0_rt; - } else { - dumdt = xdendt + 1.25e0_rt*xnumdt; - dumda = xdenda + 1.25e0_rt*xnumda; - dumdz = xdendz + 1.25e0_rt*xnumdz; - } - } - - d = 0.57e0_rt - 0.25e0_rt*xnum; - a3 = c/d; - c00 = std::exp(-a3*a3); - - fxy = 1.05e0_rt + (a1 - b1)*c00; - if constexpr (do_derivatives) { - f1 = -c00*2.0e0_rt*a3/d; - c01 = f1*(dumdt + a3*0.25e0_rt*xnumdt); - c03 = f1*(dumda + a3*0.25e0_rt*xnumda); - c04 = f1*(dumdz + a3*0.25e0_rt*xnumdz); - - fxydt = (a2*xnumdt - b2*xnumdt)*c00 + (a1-b1)*c01; - fxyda = (a2*xnumda - b2*xnumda)*c00 + (a1-b1)*c03; - fxydz = (a2*xnumdz - b2*xnumdz)*c00 + (a1-b1)*c04; - } - } - - // equation 4.1 and 4.5 - splas = (ft + fl) * fxy; - if constexpr (do_derivatives) { - splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt; - splasda = (ftda + flda)*fxy + (ft+fl)*fxyda; - splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz; - } - - a2 = std::exp(-gl); - - a1 = splas; - splas = a2*a1; - if constexpr (do_derivatives) { - a3 = -0.5e0_rt*a2*gl*gum; - splasdt = a2*splasdt + a3*gl2dt*a1; - splasda = a2*splasda + a3*gl2da*a1; - splasdz = a2*splasdz + a3*gl2dz*a1; - } - - a2 = gl6; - - a1 = splas; - splas = a2*a1; - if constexpr (do_derivatives) { - a3 = 3.0e0_rt*gl6*gum; - splasdt = a2*splasdt + a3*gl2dt*a1; - splasda = a2*splasda + a3*gl2da*a1; - splasdz = a2*splasdz + a3*gl2dz*a1; - } - - 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*sf.xl8*sf.xldt; - splasdt = a2*splasdt + a3*a1; - splasda = a2*splasda; - splasdz = a2*splasdz; - } + nu_plasma(sf, splas, splasdt, splasda, splasdz); nu_photo(sf, sphot, sphotdt, sphotda, sphotdz); From b067b49371e74fd7796f3289f6a09b5d820886a4 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 17:29:16 -0400 Subject: [PATCH 15/28] fix comp --- neutrinos/sneut5.H | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 0717160668..6a9fecab43 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -1309,11 +1309,9 @@ void sneut5(const Real temp, const Real den, dsnudz = derivative of snu with zbar */ - Real xlnt,cc, - a1,a2,a3,b1,b2,c00,c01,c03,c04, - c,d{0.0}, - f1{0.0}, - dumdt,gum,dumda,dumdz; + Real a1,a2,a3,b1,b2, + c,d{0.0}; + // pair production Real gl,gldt, @@ -1322,13 +1320,6 @@ void sneut5(const Real temp, const Real den, fpairda,fpairdz,qpairda,qpairdz, xnumda,xnumdz,xdenda,xdendz; - // plasma - Real gl2,gl2dt,gl12,gl32,gl72,gl6, - ft,ftdt,fl,fldt,fxy,fxydt, - flda,fldz,ftda,ftdz, - fxyda,fxydz,gl2da,gl2dz; - - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; From 2c92d416ab5c0c1de9ae9353cf5b506407add476 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 17:47:54 -0400 Subject: [PATCH 16/28] finish up with spair --- neutrinos/sneut5.H | 267 +++++++++++++++++++++++---------------------- 1 file changed, 136 insertions(+), 131 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 6a9fecab43..4b26c89469 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -365,6 +365,141 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_pair(const sneutf_t& sf, + Real& spair, Real& spairdt, Real& spairda, Real& spairdz) { + + // pair neutrino section + // for reactions like e+ + e- => nu_e + nubar_e + + // equation 2.8 + Real gl = 1.0e0_rt - 13.04e0_rt * sf.xl2 +133.5e0_rt * sf.xl4 +1534.0e0_rt * sf.xl6 + 918.6e0_rt * sf.xl8; + Real gldt; + if constexpr (do_derivatives) { + 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 + + Real a1 = 6.002e19_rt + 2.084e20_rt * sf.zeta + 1.872e21_rt * sf.zeta2; + Real b1, b2; + if (sf.t9 < 10.0_rt) { + b1 = std::exp(-5.5924e0_rt * sf.zeta); + if constexpr (do_derivatives) { + b2 = -b1 * 5.5924e0_rt; + } + } else { + b1 = std::exp(-4.9924e0_rt * sf.zeta); + if constexpr (do_derivatives) { + b2 = -b1 * 4.9924e0_rt; + } + } + + Real a2, c; + + Real xnum = a1 * b1; + Real xnumdt, xnumda, xnumdz; + if constexpr (do_derivatives) { + a2 = 2.084e20_rt + 2.0e0_rt * 1.872e21_rt * sf.zeta; + c = a2 * b1 + a1 * b2; + xnumdt = c * sf.zetadt; + xnumda = c * sf.zetada; + xnumdz = c * sf.zetadz; + } + + 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 * 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 * sf.xlm1 - 8.141e-1_rt * sf.xlm2; + if constexpr (do_derivatives) { + a2 = -1.2383e0_rt * sf.xlm2 + 2.0e0_rt * 8.141e-1_rt * sf.xlm3; + } + } + + Real xden = sf.zeta3 + a1; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + b1 = 3.0e0_rt * sf.zeta2; + xdendt = b1 * sf.zetadt + a2 * sf.xldt; + xdenda = b1 * sf.zetada; + xdendz = b1 * sf.zetadz; + } + + a1 = 1.0e0_rt / xden; + Real fpair = xnum * a1; + Real fpairdt, fpairda, fpairdz; + if constexpr (do_derivatives) { + fpairdt = (xnumdt - fpair * xdendt) * a1; + fpairda = (xnumda - fpair * xdenda) * a1; + fpairdz = (xnumdz - fpair * xdendz) * a1; + } + + // equation 2.6 + a1 = 10.7480e0_rt * sf.xl2 + 0.3967e0_rt * sf.xlp5 + 1.005e0_rt; + xnum = 1.0e0_rt / a1; + if constexpr (do_derivatives) { + 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 * sf.xl3 + 9.715e6_rt * sf.xlp5; + c = 1.0e0_rt / a1; + b1 = 1.0e0_rt + sf.rm * c; + + xden = std::pow(b1, -0.3e0_rt); + if constexpr (do_derivatives) { + a2 = sf.xldt * (3.0e0_rt * 7.692e7_rt * sf.xl2 + 0.5e0_rt * 9.715e6_rt * sf.xlmp5); + Real d = -0.3e0_rt * xden / b1; + xdendt = -d * sf.rm * c * c * a2; + xdenda = d * sf.rmda * c; + xdendz = d * sf.rmdz * c; + } + + Real qpair = xnum * xden; + Real qpairdt, qpairda, qpairdz; + if constexpr (do_derivatives) { + qpairdt = xnumdt * xden + xnum * xdendt; + qpairda = xnum * xdenda; + qpairdz = xnum * xdendz; + } + + // equation 2.5 + a1 = std::exp(-2.0e0_rt * sf.xlm1); + + spair = a1 * fpair; + if constexpr (do_derivatives) { + a2 = a1 * 2.0e0_rt * sf.xlm2 * sf.xldt; + spairdt = a2 * fpair + a1 * fpairdt; + spairda = a1 * fpairda; + spairdz = a1 * fpairdz; + } + + a1 = spair; + spair = gl*a1; + if constexpr (do_derivatives) { + spairdt = gl * spairdt + gldt * a1; + spairda = gl * spairda; + spairdz = gl * spairdz; + } + + a1 = nu_constants::tfac4 * (1.0e0_rt + nu_constants::tfac3 * qpair); + + Real a3 = spair; + spair = a1 * a3; + if constexpr (do_derivatives) { + 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; + } + +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void nu_plasma(const sneutf_t& sf, @@ -1309,17 +1444,6 @@ void sneut5(const Real temp, const Real den, dsnudz = derivative of snu with zbar */ - Real a1,a2,a3,b1,b2, - c,d{0.0}; - - - // pair production - Real gl,gldt, - xnum,xnumdt,xden,xdendt, - fpair,fpairdt,qpair,qpairdt, - fpairda,fpairdz,qpairda,qpairdz, - xnumda,xnumdz,xdenda,xdendz; - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; @@ -1356,126 +1480,7 @@ void sneut5(const Real temp, const Real den, auto sf = get_sneut_factors(den, temp, abar, zbar); - // pair neutrino section - // for reactions like e+ + e- => nu_e + nubar_e - - // equation 2.8 - 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 = 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 * sf.zeta + 1.872e21_rt * sf.zeta2; - - if (sf.t9 < 10.0_rt) { - b1 = std::exp(-5.5924e0_rt * sf.zeta); - if constexpr (do_derivatives) { - b2 = -b1*5.5924e0_rt; - } - } else { - b1 = std::exp(-4.9924e0_rt * sf.zeta); - if constexpr (do_derivatives) { - b2 = -b1*4.9924e0_rt; - } - } - - xnum = a1 * b1; - if constexpr (do_derivatives) { - a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt * sf.zeta; - c = a2*b1 + a1*b2; - xnumdt = c * sf.zetadt; - xnumda = c * sf.zetada; - xnumdz = c * sf.zetadz; - } - - 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*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*sf.xlm1 - 8.141e-1_rt*sf.xlm2; - if constexpr (do_derivatives) { - a2 = -1.2383e0_rt*sf.xlm2 + 2.0e0_rt*8.141e-1_rt*sf.xlm3; - } - } - - xden = sf.zeta3 + a1; - if constexpr (do_derivatives) { - b1 = 3.0e0_rt*sf.zeta2; - xdendt = b1*sf.zetadt + a2*sf.xldt; - xdenda = b1*sf.zetada; - xdendz = b1*sf.zetadz; - } - - a1 = 1.0e0_rt/xden; - fpair = xnum*a1; - if constexpr (do_derivatives) { - fpairdt = (xnumdt - fpair*xdendt)*a1; - fpairda = (xnumda - fpair*xdenda)*a1; - fpairdz = (xnumdz - fpair*xdendz)*a1; - } - - // equation 2.6 - a1 = 10.7480e0_rt*sf.xl2 + 0.3967e0_rt*sf.xlp5 + 1.005e0_rt; - xnum = 1.0e0_rt/a1; - if constexpr (do_derivatives) { - 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*sf.xl3 + 9.715e6_rt*sf.xlp5; - c = 1.0e0_rt/a1; - b1 = 1.0e0_rt + sf.rm*c; - - xden = std::pow(b1, -0.3e0_rt); - if constexpr (do_derivatives) { - 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*sf.rm*c*c*a2; - xdenda = d*sf.rmda*c; - xdendz = d*sf.rmdz*c; - } - - qpair = xnum*xden; - if constexpr (do_derivatives) { - qpairdt = xnumdt*xden + xnum*xdendt; - qpairda = xnum*xdenda; - qpairdz = xnum*xdendz; - } - - // equation 2.5 - a1 = std::exp(-2.0e0_rt*sf.xlm1); - - spair = a1*fpair; - if constexpr (do_derivatives) { - a2 = a1*2.0e0_rt*sf.xlm2*sf.xldt; - spairdt = a2*fpair + a1*fpairdt; - spairda = a1*fpairda; - spairdz = a1*fpairdz; - } - - a1 = spair; - spair = gl*a1; - if constexpr (do_derivatives) { - spairdt = gl*spairdt + gldt*a1; - spairda = gl*spairda; - spairdz = gl*spairdz; - } - - a1 = nu_constants::tfac4*(1.0e0_rt + nu_constants::tfac3 * qpair); - - a3 = spair; - spair = a1*a3; - if constexpr (do_derivatives) { - 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; - } - + nu_pair(sf, spair, spairdt, spairda, spairdz); nu_plasma(sf, splas, splasdt, splasda, splasdz); From b95912da29587537a48d2c55d5ac35225a1c757d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 19:08:23 -0400 Subject: [PATCH 17/28] use multiple angle formula for sin/cos in sneut --- neutrinos/sneut5.H | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 4b26c89469..b45c5167eb 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -824,21 +824,28 @@ void nu_photo(const sneutf_t& sf, Real taudt = nu_constants::iln10 * sf.tempi; - // equation 3.7, compute the expensive trig functions only one time - - Real cos1 = std::cos(nu_constants::fac1 * tau); - Real cos2 = std::cos(nu_constants::fac1 * 2.0e0_rt * tau); - Real cos3 = std::cos(nu_constants::fac1 * 3.0e0_rt * tau); - Real cos4 = std::cos(nu_constants::fac1 * 4.0e0_rt * tau); - Real cos5 = std::cos(nu_constants::fac1 * 5.0e0_rt * tau); - Real last = std::cos(nu_constants::fac2 * tau); - - Real sin1 = std::sin(nu_constants::fac1*tau); - Real sin2 = std::sin(nu_constants::fac1*2.0e0_rt*tau); - Real sin3 = std::sin(nu_constants::fac1*3.0e0_rt*tau); - Real sin4 = std::sin(nu_constants::fac1*4.0e0_rt*tau); - Real sin5 = std::sin(nu_constants::fac1*5.0e0_rt*tau); - Real xast = std::sin(nu_constants::fac2*tau); + // equation 3.7 + + const Real [cos1, sin1] = amrex::sincos(nu_constants::fac1 * tau); + + // double, triple, etc. angle formulas + // sin/cos (2 fac1 tau) + const Real sin2 = 2.0_rt * sin1 * cos1; + const Real cos2 = 2.0_rt * cos1 * cos1 - 1.0_rt; + + // sin/cos (3 fac1 tau) + const Real sin3 = 3.0_rt * sin1 - 4.0_rt * sin1 * sin1 * sin1; + const Real cos3 = 4.0_rt * cos1 * cos1 - cos1 - 3.0_rt * cos1; + + // sin/cos (4 fac1 tau) -- use double angle on sin2/cos2 + const Real sin4 = 2.0_rt * sin2 * cos2; + const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; + + // sin/cos (5 fac1 tau) + const Real sin5 = 5.0_rt * sin1 - 20.0_rt * sin1 * sin1 * sin1 + 16.0_rt * sin1 * sin1 * sin1 * sin1 * sin; + const Real cos5 = 16.0_rt * cos1 * cos1 * cos1 * cos1 * cos1 - 20.0_rt * cos1 * cos1 * cos1 + 5.0_rt * cos1; + + const Real [xast, last] = amrex::sincos(fac2 * tau); Real a0 = 0.5e0_rt * c00 + c01 * cos1 + dd01 * sin1 + c02 * cos2 + dd02 * sin2 From 26f79d3aa58aafcc9fba67b08c2d72ac034b9a28 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 19:20:14 -0400 Subject: [PATCH 18/28] fix compilation --- neutrinos/sneut5.H | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index b45c5167eb..2c25818d3d 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -3,6 +3,7 @@ #include #include +#include using namespace amrex; @@ -826,7 +827,7 @@ void nu_photo(const sneutf_t& sf, // equation 3.7 - const Real [cos1, sin1] = amrex::sincos(nu_constants::fac1 * tau); + const auto [cos1, sin1] = amrex::Math::sincos(nu_constants::fac1 * tau); // double, triple, etc. angle formulas // sin/cos (2 fac1 tau) @@ -842,10 +843,10 @@ void nu_photo(const sneutf_t& sf, const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; // sin/cos (5 fac1 tau) - const Real sin5 = 5.0_rt * sin1 - 20.0_rt * sin1 * sin1 * sin1 + 16.0_rt * sin1 * sin1 * sin1 * sin1 * sin; + const Real sin5 = 5.0_rt * sin1 - 20.0_rt * sin1 * sin1 * sin1 + 16.0_rt * sin1 * sin1 * sin1 * sin1 * sin1; const Real cos5 = 16.0_rt * cos1 * cos1 * cos1 * cos1 * cos1 - 20.0_rt * cos1 * cos1 * cos1 + 5.0_rt * cos1; - const Real [xast, last] = amrex::sincos(fac2 * tau); + const auto [xast, last] = amrex::Math::sincos(nu_constants::fac2 * tau); Real a0 = 0.5e0_rt * c00 + c01 * cos1 + dd01 * sin1 + c02 * cos2 + dd02 * sin2 From 0dff5470490bccaaec0ff55e2097f627c09b5224 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 19:36:53 -0400 Subject: [PATCH 19/28] some simplifying --- neutrinos/sneut5.H | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 2c25818d3d..fdf11e160c 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -835,16 +835,16 @@ void nu_photo(const sneutf_t& sf, const Real cos2 = 2.0_rt * cos1 * cos1 - 1.0_rt; // sin/cos (3 fac1 tau) - const Real sin3 = 3.0_rt * sin1 - 4.0_rt * sin1 * sin1 * sin1; - const Real cos3 = 4.0_rt * cos1 * cos1 - cos1 - 3.0_rt * cos1; + const Real sin3 = sin1 * (3.0_rt - 4.0_rt * sin1 * sin1); + const Real cos3 = cos1 * (4.0_rt * cos1 * cos1 - 3.0_rt); // sin/cos (4 fac1 tau) -- use double angle on sin2/cos2 const Real sin4 = 2.0_rt * sin2 * cos2; const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; // sin/cos (5 fac1 tau) - const Real sin5 = 5.0_rt * sin1 - 20.0_rt * sin1 * sin1 * sin1 + 16.0_rt * sin1 * sin1 * sin1 * sin1 * sin1; - const Real cos5 = 16.0_rt * cos1 * cos1 * cos1 * cos1 * cos1 - 20.0_rt * cos1 * cos1 * cos1 + 5.0_rt * cos1; + const Real sin5 = 5.0_rt * sin1 - sin1 * sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1); + const Real cos5 = cos1 * cos1 * cos1 * (16.0_rt * cos1 * cos1 - 20.0_rt) + 5.0_rt * cos1; const auto [xast, last] = amrex::Math::sincos(nu_constants::fac2 * tau); From 0f7b90fbbad4f45452f1908ef1dc7b755d178ba9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 19:43:18 -0400 Subject: [PATCH 20/28] a little more factoring --- neutrinos/sneut5.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index fdf11e160c..ec8281bf97 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -843,8 +843,8 @@ void nu_photo(const sneutf_t& sf, const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; // sin/cos (5 fac1 tau) - const Real sin5 = 5.0_rt * sin1 - sin1 * sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1); - const Real cos5 = cos1 * cos1 * cos1 * (16.0_rt * cos1 * cos1 - 20.0_rt) + 5.0_rt * cos1; + const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1); + const Real cos5 = cos1 * (cos1 * cos1 * (16.0_rt * cos1 * cos1 - 20.0_rt) + 5.0_rt); const auto [xast, last] = amrex::Math::sincos(nu_constants::fac2 * tau); From d524e5ccfa3b645db35fd2eb03a3b65cfc19595c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 19:45:55 -0400 Subject: [PATCH 21/28] fix compilation --- neutrinos/sneut5.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index ec8281bf97..1ae4147e1a 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -843,7 +843,7 @@ void nu_photo(const sneutf_t& sf, const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; // sin/cos (5 fac1 tau) - const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1); + const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1)); const Real cos5 = cos1 * (cos1 * cos1 * (16.0_rt * cos1 * cos1 - 20.0_rt) + 5.0_rt); const auto [xast, last] = amrex::Math::sincos(nu_constants::fac2 * tau); From 93e344ccd5e5c31bfd67b18ebe3ae109c2f7a1e5 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 20:30:56 -0400 Subject: [PATCH 22/28] Update neutrinos/sneut5.H Co-authored-by: Eric T. Johnson --- neutrinos/sneut5.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 1ae4147e1a..87c9711d3e 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -843,7 +843,7 @@ void nu_photo(const sneutf_t& sf, const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; // sin/cos (5 fac1 tau) - const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1)); + const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt - 16.0_rt * sin1 * sin1)); const Real cos5 = cos1 * (cos1 * cos1 * (16.0_rt * cos1 * cos1 - 20.0_rt) + 5.0_rt); const auto [xast, last] = amrex::Math::sincos(nu_constants::fac2 * tau); From 8f46a783aedb81c9521afd4997715fe775673bb2 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 30 Oct 2023 08:32:28 -0400 Subject: [PATCH 23/28] finish sincos --- neutrinos/sneut5.H | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 1ae4147e1a..340fc29979 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -1174,17 +1174,25 @@ void nu_brem(const sneutf_t& sf, //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); + + const auto [cos1, sin1] = amrex::Math::sincos(u); + + // double, triple, etc. angle formulas + // sin/cos (2 fac1 tau) + const Real sin2 = 2.0_rt * sin1 * cos1; + const Real cos2 = 2.0_rt * cos1 * cos1 - 1.0_rt; + + // sin/cos (3 fac1 tau) + const Real sin3 = sin1 * (3.0_rt - 4.0_rt * sin1 * sin1); + const Real cos3 = cos1 * (4.0_rt * cos1 * cos1 - 3.0_rt); + + // sin/cos (4 fac1 tau) -- use double angle on sin2/cos2 + const Real sin4 = 2.0_rt * sin2 * cos2; + const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; + + // sin/cos (5 fac1 tau) + //const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1)); + const Real cos5 = cos1 * (cos1 * cos1 * (16.0_rt * cos1 * cos1 - 20.0_rt) + 5.0_rt); // equation 5.21 Real fb = 0.5e0_rt * 0.17946e0_rt + 0.00945e0_rt * u + 0.34529e0_rt From 0e0c3af8e1b8029de800018348b5aec5052964d6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 30 Oct 2023 08:33:19 -0400 Subject: [PATCH 24/28] finish --- neutrinos/sneut5.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index c96876e381..9770050deb 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -1191,7 +1191,7 @@ void nu_brem(const sneutf_t& sf, const Real cos4 = 2.0_rt * cos2 * cos2 - 1.0_rt; // sin/cos (5 fac1 tau) - //const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt + 16.0_rt * sin1 * sin1)); + //const Real sin5 = sin1 * (5.0_rt - sin1 * sin1 * (20.0_rt - 16.0_rt * sin1 * sin1)); const Real cos5 = cos1 * (cos1 * cos1 * (16.0_rt * cos1 * cos1 - 20.0_rt) + 5.0_rt); // equation 5.21 From d7a51e23304884ae06bd074ee6f74cd724edca3b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 30 Oct 2023 08:53:14 -0400 Subject: [PATCH 25/28] update test_rhs --- unit_test/test_rhs/ci-benchmarks/aprox13.out | 12 ++++++------ unit_test/test_rhs/ci-benchmarks/aprox19.out | 12 ++++++------ unit_test/test_rhs/ci-benchmarks/aprox21.out | 14 +++++++------- unit_test/test_rhs/ci-benchmarks/ecsn.out | 14 +++++++------- unit_test/test_rhs/ci-benchmarks/iso7.out | 4 ++-- 5 files changed, 28 insertions(+), 28 deletions(-) diff --git a/unit_test/test_rhs/ci-benchmarks/aprox13.out b/unit_test/test_rhs/ci-benchmarks/aprox13.out index 2a3cb6c119..56c4ec226c 100644 --- a/unit_test/test_rhs/ci-benchmarks/aprox13.out +++ b/unit_test/test_rhs/ci-benchmarks/aprox13.out @@ -43,7 +43,7 @@ J_chromium-48_helium-4 6.8935846419e-69 5498823.3125 J_iron-52_helium-4 5.1159900994e-72 9496238.1457 J_nickel-56_helium-4 1.2369463711e-79 6465258.5107 - J_E_helium-4 -37.96418269 5.6963666435e+27 + J_E_helium-4 -37.961886363 5.6963666435e+27 J_helium-4_carbon-12 -549068.13437 16350655002 J_carbon-12_carbon-12 -32736863803 -4.3536076636e-25 J_oxygen-16_carbon-12 -215702563.02 53050483.792 @@ -57,7 +57,7 @@ J_chromium-48_carbon-12 0 0 J_iron-52_carbon-12 0 0 J_nickel-56_carbon-12 0 0 - J_E_carbon-12 -7.4476136442e+22 7.3277248697e+28 + J_E_carbon-12 -7.4476136443e+22 7.3277248697e+28 J_helium-4_oxygen-16 -196341285.62 116442916.35 J_carbon-12_oxygen-16 -295329776.6 144805.77929 J_oxygen-16_oxygen-16 -328244197.95 -5.8277893387e-34 @@ -71,7 +71,7 @@ J_chromium-48_oxygen-16 0 0 J_iron-52_oxygen-16 0 0 J_nickel-56_oxygen-16 0 0 - J_E_oxygen-16 -9.7641101938e+23 3.501361507e+27 + J_E_oxygen-16 -9.7641101939e+23 3.501361507e+27 J_helium-4_neon-20 -10923354794 114705723.28 J_carbon-12_neon-20 0 0 J_oxygen-16_neon-20 0 144408911.44 @@ -99,7 +99,7 @@ J_chromium-48_magnesium-24 0 0 J_iron-52_magnesium-24 0 0 J_nickel-56_magnesium-24 0 0 - J_E_magnesium-24 -65676.173201 1.652673424e+29 + J_E_magnesium-24 -65734.443865 1.652673424e+29 J_helium-4_silicon-28 -2556508974.5 20661.78792 J_carbon-12_silicon-28 0 0 J_oxygen-16_silicon-28 0 0 @@ -113,7 +113,7 @@ J_chromium-48_silicon-28 0 0 J_iron-52_silicon-28 0 0 J_nickel-56_silicon-28 0 0 - J_E_silicon-28 -2.904944194e+23 1.7139561264e+28 + J_E_silicon-28 -2.9049441941e+23 1.7139561264e+28 J_helium-4_sulfur-32 -4251265603.3 8881309.5672 J_carbon-12_sulfur-32 0 0 J_oxygen-16_sulfur-32 0 0 @@ -211,7 +211,7 @@ J_chromium-48_nickel-56 0 0 J_iron-52_nickel-56 0 381082.60535 J_nickel-56_nickel-56 -381082.60535 0 - J_E_nickel-56 -2.9411320544e+24 8.0628253938e+14 + J_E_nickel-56 -2.9411320544e+24 8.0550862785e+14 J_helium-4_E -1.375359511e-10 9.1794090251e-09 J_carbon-12_E -1.8412378137e-08 6.8308749195e-13 J_oxygen-16_E -4.9194133409e-11 4.971951643e-12 diff --git a/unit_test/test_rhs/ci-benchmarks/aprox19.out b/unit_test/test_rhs/ci-benchmarks/aprox19.out index 201258eb2d..91c2d9702a 100644 --- a/unit_test/test_rhs/ci-benchmarks/aprox19.out +++ b/unit_test/test_rhs/ci-benchmarks/aprox19.out @@ -201,7 +201,7 @@ J_nickel-56_magnesium-24 0 0 J_neutron_magnesium-24 0 0 J_proton_magnesium-24 0 0 - J_E_magnesium-24 -65836.025628 1.6435693578e+29 + J_E_magnesium-24 -65859.317238 1.6435693578e+29 J_hydrogen-1_silicon-28 0 0 J_helium-3_silicon-28 0 0 J_helium-4_silicon-28 -2539600300.1 20658.714693 @@ -221,7 +221,7 @@ J_nickel-56_silicon-28 0 0 J_neutron_silicon-28 0 0 J_proton_silicon-28 0 0 - J_E_silicon-28 -2.9044184049e+23 1.7026200525e+28 + J_E_silicon-28 -2.904418405e+23 1.7026200525e+28 J_hydrogen-1_sulfur-32 0 0 J_helium-3_sulfur-32 0 0 J_helium-4_sulfur-32 -4220786865.5 8879079.2827 @@ -341,7 +341,7 @@ J_nickel-56_iron-52 4.3001564738e-78 552709472.84 J_neutron_iron-52 -6.4991725027e+12 0 J_proton_iron-52 0 431335119.97 - J_E_iron-52 -2.3084042815e+11 7.5442123176e+31 + J_E_iron-52 -2.2996831358e+11 7.5442123176e+31 J_hydrogen-1_iron-54 0 0 J_helium-3_iron-54 0 0 J_helium-4_iron-54 0 14675163.343 @@ -361,7 +361,7 @@ J_nickel-56_iron-54 0 14319319597 J_neutron_iron-54 0 604.52588062 J_proton_iron-54 -28667989521 0 - J_E_iron-54 -6.4881996912e+21 1.6907715174e+29 + J_E_iron-54 -6.488199697e+21 1.6907715174e+29 J_hydrogen-1_nickel-56 0 0 J_helium-3_nickel-56 0 0 J_helium-4_nickel-56 0 334649.39629 @@ -381,7 +381,7 @@ J_nickel-56_nickel-56 -247479689.06 0 J_neutron_nickel-56 0 0 J_proton_nickel-56 0 494955656.77 - J_E_nickel-56 -2.9211034564e+27 1.0288702015e+17 + J_E_nickel-56 -2.9211034564e+27 1.0288701986e+17 J_hydrogen-1_neutron 0 0 J_helium-3_neutron 0 0 J_helium-4_neutron 0 1.5520970787e+11 @@ -421,7 +421,7 @@ J_nickel-56_proton 0 530345170.26 J_neutron_proton -3.1041941575e+11 0.01005978836 J_proton_proton -3.1047599221e+11 -2.3827182398e-87 - J_E_proton -7.593501065e+15 4.237794952e+30 + J_E_proton -7.593501066e+15 4.237794952e+30 J_hydrogen-1_E 2.0563926691e-19 5.7718078776e-06 J_helium-3_E -5.7716864292e-06 -2.0648367761e-19 J_helium-4_E 1.026556653e-19 2.885937043e-06 diff --git a/unit_test/test_rhs/ci-benchmarks/aprox21.out b/unit_test/test_rhs/ci-benchmarks/aprox21.out index bc513f9e62..b8c6ee7383 100644 --- a/unit_test/test_rhs/ci-benchmarks/aprox21.out +++ b/unit_test/test_rhs/ci-benchmarks/aprox21.out @@ -221,7 +221,7 @@ J_nickel-56_magnesium-24 0 0 J_neutron_magnesium-24 0 0 J_proton_magnesium-24 0 0 - J_E_magnesium-24 -71759.536816 1.6357279567e+29 + J_E_magnesium-24 -71764.239466 1.6357279567e+29 J_hydrogen-1_silicon-28 0 0 J_helium-3_silicon-28 0 0 J_helium-4_silicon-28 -2525055183.6 20658.704958 @@ -243,7 +243,7 @@ J_nickel-56_silicon-28 0 0 J_neutron_silicon-28 0 0 J_proton_silicon-28 0 0 - J_E_silicon-28 -2.9044167997e+23 1.6928685775e+28 + J_E_silicon-28 -2.9044167998e+23 1.6928685775e+28 J_hydrogen-1_sulfur-32 0 0 J_helium-3_sulfur-32 0 0 J_helium-4_sulfur-32 -4194578244.6 8879072.2185 @@ -375,7 +375,7 @@ J_nickel-56_chromium-56 0 0 J_neutron_chromium-56 0 0 J_proton_chromium-56 0 0 - J_E_chromium-56 -1.8093693678e+15 2.256453157e+17 + J_E_chromium-56 -1.8292382081e+15 2.2563985305e+17 J_hydrogen-1_iron-52 0 0 J_helium-3_iron-52 0 0 J_helium-4_iron-52 -743179296.84 493082.01674 @@ -397,7 +397,7 @@ J_nickel-56_iron-52 4.5459700269e-78 539742599.12 J_neutron_iron-52 -5.7770422246e+12 0 J_proton_iron-52 0 463978849.94 - J_E_iron-52 -2.4561809471e+11 6.7059665045e+31 + J_E_iron-52 -2.4468743487e+11 6.7059665045e+31 J_hydrogen-1_iron-54 0 0 J_helium-3_iron-54 0 0 J_helium-4_iron-54 -416714108.26 89938.022334 @@ -419,7 +419,7 @@ J_nickel-56_iron-54 0 12628008866 J_neutron_iron-54 -5.8420788754e+13 0 J_proton_iron-54 -25242800365 49264281.385 - J_E_iron-54 -2.5384087204e+11 5.7722247891e+32 + J_E_iron-54 -2.529102122e+11 5.7722247891e+32 J_hydrogen-1_iron-56 0 0 J_helium-3_iron-56 0 0 J_helium-4_iron-56 0 1.0853220404e+11 @@ -463,7 +463,7 @@ J_nickel-56_nickel-56 -247303719.77 0 J_neutron_nickel-56 0 0 J_proton_nickel-56 0 494604136 - J_E_nickel-56 -2.9190272674e+27 3.8630535107e+17 + J_E_nickel-56 -2.9190272674e+27 3.863053733e+17 J_hydrogen-1_neutron 0 0 J_helium-3_neutron 0 0 J_helium-4_neutron 0 1.3482425469e+11 @@ -507,7 +507,7 @@ J_nickel-56_proton 0 467704032.09 J_neutron_proton -2.6964850938e+11 0.034588911424 J_proton_proton -2.6983923086e+11 -2.7428164302e-87 - J_E_proton -2.6108991385e+16 3.6817330507e+30 + J_E_proton -2.6108991386e+16 3.6817330507e+30 J_hydrogen-1_E -1.0006504976e-07 4.2872891104e-06 J_helium-3_E -4.2866803458e-06 -7.5769167747e-20 J_helium-4_E 3.7850094111e-20 2.1441607162e-06 diff --git a/unit_test/test_rhs/ci-benchmarks/ecsn.out b/unit_test/test_rhs/ci-benchmarks/ecsn.out index 9489ba5c4d..29c2eedc18 100644 --- a/unit_test/test_rhs/ci-benchmarks/ecsn.out +++ b/unit_test/test_rhs/ci-benchmarks/ecsn.out @@ -37,7 +37,7 @@ J_silicon-28_hydrogen-1 1.9049142873e-28 1.4457599197e+16 J_phosphorus-31_hydrogen-1 -1.4363730396e+16 -2.3890907984e-32 J_sulfur-32_hydrogen-1 2.3875895211e-32 3.8373494834e+15 - J_E_hydrogen-1 7.2668109011e-07 5.3233598076e+34 + J_E_hydrogen-1 7.3588532259e-07 5.3233598076e+34 J_hydrogen-1_helium-4 6.3597516852e-96 2.3050360212e+15 J_helium-4_helium-4 -2.3251430676e+15 -8.7617347065e-57 J_oxygen-16_helium-4 -5.9750641231e+12 -8.7617347064e-57 @@ -61,7 +61,7 @@ J_silicon-28_oxygen-16 5.3544432665e-93 1.4407576856e+26 J_phosphorus-31_oxygen-16 5.3544432665e-93 1.4407576856e+26 J_sulfur-32_oxygen-16 0 0 - J_E_oxygen-16 -54524.316756 2.4009435853e+45 + J_E_oxygen-16 -54556.178444 2.4009435853e+45 J_hydrogen-1_oxygen-20 0 0 J_helium-4_oxygen-20 0 0 J_oxygen-16_oxygen-20 0 0 @@ -73,7 +73,7 @@ J_silicon-28_oxygen-20 0 0 J_phosphorus-31_oxygen-20 0 0 J_sulfur-32_oxygen-20 0 0 - J_E_oxygen-20 -4.2481653689e+13 1.7104309459e+17 + J_E_oxygen-20 -4.2470934331e+13 1.7103809712e+17 J_hydrogen-1_fluorine-20 0 0 J_helium-4_fluorine-20 0 0 J_oxygen-16_fluorine-20 0 0 @@ -85,7 +85,7 @@ J_silicon-28_fluorine-20 0 0 J_phosphorus-31_fluorine-20 0 0 J_sulfur-32_fluorine-20 0 0 - J_E_fluorine-20 -1.5940355629e+21 4.7992907772e+17 + J_E_fluorine-20 -1.5940355628e+21 4.7992345557e+17 J_hydrogen-1_neon-20 0 0 J_helium-4_neon-20 -3.9073008309e+13 87637435.641 J_oxygen-16_neon-20 3.8849945775e-100 88180925.52 @@ -121,7 +121,7 @@ J_silicon-28_aluminum-27 5.1432681704e-27 3.7588755703e+16 J_phosphorus-31_aluminum-27 6.4778539534e-78 6.4064066131e+13 J_sulfur-32_aluminum-27 0 0 - J_E_aluminum-27 -2.3076325355e-05 4.5339945455e+35 + J_E_aluminum-27 -2.3473767915e-05 4.5339945455e+35 J_hydrogen-1_silicon-28 1.3633533854e-95 4.2286361267e+15 J_helium-4_silicon-28 -4.2420142918e+15 -1.9556484036e-82 J_oxygen-16_silicon-28 0 0 @@ -145,7 +145,7 @@ J_silicon-28_phosphorus-31 4.6539596654e-34 4.4390283525e+17 J_phosphorus-31_phosphorus-31 -4.4527564228e+17 -7.406181475e-31 J_sulfur-32_phosphorus-31 7.4015275153e-31 1.1895783399e+17 - J_E_phosphorus-31 -0.00018025525389 1.1296769894e+36 + J_E_phosphorus-31 -0.00021865701493 1.1296769894e+36 J_hydrogen-1_sulfur-32 0 0 J_helium-4_sulfur-32 0 0 J_oxygen-16_sulfur-32 0 0 @@ -157,7 +157,7 @@ J_silicon-28_sulfur-32 0 0 J_phosphorus-31_sulfur-32 0 0 J_sulfur-32_sulfur-32 0 0 - J_E_sulfur-32 -1.2066886642e+15 7.8006383744e+14 + J_E_sulfur-32 -1.2200141114e+15 7.734440401e+14 J_hydrogen-1_E -0.12580793235 1.2089803642e+12 J_helium-4_E -0.0078065301852 1.4905283338e+12 J_oxygen-16_E -5.3990173959e+12 1.2739583515e-06 diff --git a/unit_test/test_rhs/ci-benchmarks/iso7.out b/unit_test/test_rhs/ci-benchmarks/iso7.out index f77f18f11e..86f9a2a563 100644 --- a/unit_test/test_rhs/ci-benchmarks/iso7.out +++ b/unit_test/test_rhs/ci-benchmarks/iso7.out @@ -25,7 +25,7 @@ J_magnesium-24_helium-4 -4.6598486093e-09 14984957882 J_silicon-28_helium-4 -1.0350391329e+35 7196378096.4 J_nickel-56_helium-4 0 1.0350391329e+35 - J_E_helium-4 -87.837947379 4.932675415e+54 + J_E_helium-4 -87.77240322 4.932675415e+54 J_helium-4_carbon-12 -46065373.422 1.8369106708e+13 J_carbon-12_carbon-12 -3.6773326651e+13 -7.0444246309e-43 J_oxygen-16_carbon-12 -7.4641202755e+11 25858392627 @@ -73,7 +73,7 @@ J_magnesium-24_nickel-56 0 0 J_silicon-28_nickel-56 0 10000000000 J_nickel-56_nickel-56 -10000000000 0 - J_E_nickel-56 -4.7656897776e+29 3.4603369325e+17 + J_E_nickel-56 -4.7656897776e+29 3.4732921378e+17 J_helium-4_E -5.2043289985e-09 1.0072596721e+20 J_carbon-12_E -8.1026490435e-07 6.1135331864e-10 J_oxygen-16_E -1.0316135993e-08 2.5415462665e-08 From c46453dcfeef9be95cf9832c3e347663f17ea8db Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 30 Oct 2023 09:01:57 -0400 Subject: [PATCH 26/28] more fixes --- .../ci-benchmarks/aprox13_QSS_unit_test.out | 60 +++++------ .../ci-benchmarks/aprox13_RKC_unit_test.out | 40 +++---- .../ci-benchmarks/ecsn_unit_test.out | 50 ++++----- .../subch_approx_BE_unit_test.out | 94 ++++++++-------- .../ci-benchmarks/subch_approx_unit_test.out | 100 +++++++++--------- .../ci-benchmarks/test_sneut5.out | 4 +- 6 files changed, 174 insertions(+), 174 deletions(-) diff --git a/unit_test/burn_cell/ci-benchmarks/aprox13_QSS_unit_test.out b/unit_test/burn_cell/ci-benchmarks/aprox13_QSS_unit_test.out index ac6df607e6..955440d3aa 100644 --- a/unit_test/burn_cell/ci-benchmarks/aprox13_QSS_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/aprox13_QSS_unit_test.out @@ -1,4 +1,4 @@ -AMReX (22.11-1-gc4a4811c373d) initialized +AMReX (23.10-13-gecaa14456e3f) initialized starting the single zone burn... Maximum Time (s): 0.01 State Density (g/cm^3): 1000000 @@ -32,41 +32,41 @@ RHS at t = 0 Ni56 4.07828597e-29 ------------------------------------ successful? 1 - - Hnuc = 3.448624864e+19 - - added e = 3.448624864e+17 - - final T = 3173085674 + - Hnuc = 3.448624858e+19 + - added e = 3.448624858e+17 + - final T = 3173085673 ------------------------------------ e initial = 1.389440599e+18 -e final = 1.734303086e+18 +e final = 1.734303085e+18 ------------------------------------ new mass fractions: -He4 0.7718339061 -C12 0.0001227479966 -O16 8.071653239e-06 -Ne20 3.560425108e-07 -Mg24 6.162605681e-07 -Si28 5.774651995e-06 -S32 5.291893286e-06 -Ar36 8.531081471e-06 -Ca40 0.0001668187948 -Ti44 0.0004124607789 -Cr48 0.000710075351 -Fe52 0.003115481694 +He4 0.771833906 +C12 0.0001227479967 +O16 8.071653242e-06 +Ne20 3.56042511e-07 +Mg24 6.162605686e-07 +Si28 5.774651999e-06 +S32 5.291893291e-06 +Ar36 8.531081479e-06 +Ca40 0.000166818795 +Ti44 0.0004124607795 +Cr48 0.0007100753521 +Fe52 0.003115481699 Ni56 0.2236098678 ------------------------------------ species creation rates: -omegadot(He4): -22.81660939 -omegadot(C12): 0.01227479966 -omegadot(O16): 0.0008071653239 -omegadot(Ne20): 3.560425108e-05 -omegadot(Mg24): 6.162605681e-05 -omegadot(Si28): 0.0005774651995 -omegadot(S32): 0.0005291893286 -omegadot(Ar36): 0.0008531081471 -omegadot(Ca40): 0.01668187948 -omegadot(Ti44): 0.04124607789 -omegadot(Cr48): 0.0710075351 -omegadot(Fe52): 0.3115481694 +omegadot(He4): -22.8166094 +omegadot(C12): 0.01227479967 +omegadot(O16): 0.0008071653242 +omegadot(Ne20): 3.56042511e-05 +omegadot(Mg24): 6.162605686e-05 +omegadot(Si28): 0.0005774651999 +omegadot(S32): 0.0005291893291 +omegadot(Ar36): 0.0008531081479 +omegadot(Ca40): 0.0166818795 +omegadot(Ti44): 0.04124607795 +omegadot(Cr48): 0.07100753521 +omegadot(Fe52): 0.3115481699 omegadot(Ni56): 22.36098678 number of steps taken: 3695 -AMReX (22.11-1-gc4a4811c373d) finalized +AMReX (23.10-13-gecaa14456e3f) finalized diff --git a/unit_test/burn_cell/ci-benchmarks/aprox13_RKC_unit_test.out b/unit_test/burn_cell/ci-benchmarks/aprox13_RKC_unit_test.out index 3d539cc6e8..4a26dea63c 100644 --- a/unit_test/burn_cell/ci-benchmarks/aprox13_RKC_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/aprox13_RKC_unit_test.out @@ -1,4 +1,4 @@ -AMReX (23.06-6-g1e1c433b401c) initialized +AMReX (23.10-13-gecaa14456e3f) initialized starting the single zone burn... Maximum Time (s): 0.01 State Density (g/cm^3): 1000000 @@ -32,8 +32,8 @@ RHS at t = 0 Ni56 3.938787868e-40 ------------------------------------ successful? 1 - - Hnuc = 8.558390239e+18 - - added e = 8.558390239e+16 + - Hnuc = 8.558390237e+18 + - added e = 8.558390237e+16 - final T = 1516425860 ------------------------------------ e initial = 1.284393683e+17 @@ -42,31 +42,31 @@ e final = 2.140232707e+17 new mass fractions: He4 0.8760723967 C12 0.1064099566 -O16 0.0001403204362 +O16 0.0001403204361 Ne20 8.129701233e-05 Mg24 0.0002972369574 Si28 0.01113151728 -S32 0.005297014797 -Ar36 0.0005641483161 -Ca40 6.109847027e-06 -Ti44 2.004306621e-09 -Cr48 3.272104798e-14 -Fe52 7.12614193e-20 -Ni56 1.224820021e-26 +S32 0.005297014786 +Ar36 0.0005641483117 +Ca40 6.109846921e-06 +Ti44 2.004306558e-09 +Cr48 3.272104624e-14 +Fe52 7.126141319e-20 +Ni56 1.224819862e-26 ------------------------------------ species creation rates: omegadot(He4): -12.39276033 omegadot(C12): 10.64099566 -omegadot(O16): 0.01403204362 +omegadot(O16): 0.01403204361 omegadot(Ne20): 0.008129701233 omegadot(Mg24): 0.02972369574 omegadot(Si28): 1.113151728 -omegadot(S32): 0.5297014797 -omegadot(Ar36): 0.05641483161 -omegadot(Ca40): 0.0006109847027 -omegadot(Ti44): 2.004306621e-07 -omegadot(Cr48): 3.272104798e-12 -omegadot(Fe52): 7.12614193e-18 -omegadot(Ni56): 1.224720021e-24 +omegadot(S32): 0.5297014786 +omegadot(Ar36): 0.05641483117 +omegadot(Ca40): 0.0006109846921 +omegadot(Ti44): 2.004306558e-07 +omegadot(Cr48): 3.272104624e-12 +omegadot(Fe52): 7.126141319e-18 +omegadot(Ni56): 1.224719862e-24 number of steps taken: 255 -AMReX (23.06-6-g1e1c433b401c) finalized +AMReX (23.10-13-gecaa14456e3f) finalized diff --git a/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out b/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out index b50ff63d29..ea545c1a30 100644 --- a/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out @@ -1,4 +1,4 @@ -AMReX (23.05-4-ga393d7ff7e32) initialized +AMReX (23.10-13-gecaa14456e3f) initialized starting the single zone burn... reading in network electron-capture / beta-decay tables... Maximum Time (s): 0.01 @@ -29,37 +29,37 @@ RHS at t = 0 S32 1323519.174 ------------------------------------ successful? 1 - - Hnuc = 4.4691188e+19 - - added e = 4.4691188e+17 - - final T = 6711517341 + - Hnuc = 4.555404693e+19 + - added e = 4.555404693e+17 + - final T = 6735967493 ------------------------------------ e initial = 5.995956082e+17 -e final = 1.046507488e+18 +e final = 1.055136077e+18 ------------------------------------ new mass fractions: -H1 0.007036728874 -He4 4.569221711e-13 -O16 6.11671909e-07 -O20 0.009363850335 -F20 0.009362566168 -Ne20 8.595667728e-14 -Mg24 8.576669062e-06 -Al27 9.999999996e-31 -Si28 0.2286832772 -P31 9.999999996e-31 -S32 0.745544389 +H1 0.008595815549 +He4 4.152700815e-13 +O16 6.210300397e-07 +O20 0.009741776606 +F20 0.009740438851 +Ne20 9.999999868e-31 +Mg24 4.491783724e-09 +Al27 7.648139704e-21 +Si28 0.2708503134 +P31 8.806134354e-14 +S32 0.7010710301 ------------------------------------ species creation rates: -omegadot(H1): -0.2963271126 +omegadot(H1): -0.1404184451 omegadot(He4): -3 -omegadot(O16): -49.99993883 -omegadot(O20): -0.0636149665 -omegadot(F20): -0.06374338321 +omegadot(O16): -49.9999379 +omegadot(O20): -0.02582233936 +omegadot(F20): -0.02595611491 omegadot(Ne20): -30 -omegadot(Mg24): -9.999142333 +omegadot(Mg24): -9.999999551 omegadot(Al27): -1 -omegadot(Si28): 21.86832772 +omegadot(Si28): 26.08503134 omegadot(P31): -1 -omegadot(S32): 73.5544389 -number of steps taken: 14714 -AMReX (23.05-4-ga393d7ff7e32) finalized +omegadot(S32): 69.10710301 +number of steps taken: 22812 +AMReX (23.10-13-gecaa14456e3f) finalized diff --git a/unit_test/burn_cell/ci-benchmarks/subch_approx_BE_unit_test.out b/unit_test/burn_cell/ci-benchmarks/subch_approx_BE_unit_test.out index df93b6b562..e26539ee56 100644 --- a/unit_test/burn_cell/ci-benchmarks/subch_approx_BE_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/subch_approx_BE_unit_test.out @@ -1,4 +1,4 @@ -AMReX (23.05-4-ga393d7ff7e32) initialized +AMReX (23.10-13-gecaa14456e3f) initialized starting the single zone burn... reading in network electron-capture / beta-decay tables... Maximum Time (s): 1e-05 @@ -51,59 +51,59 @@ RHS at t = 0 Ni56 2.574225196e-29 ------------------------------------ successful? 1 - - Hnuc = 7.242769111e+22 - - added e = 7.242769111e+17 - - final T = 3332749634 + - Hnuc = 7.242779628e+22 + - added e = 7.242779628e+17 + - final T = 3332750042 ------------------------------------ e initial = 1.396711859e+18 -e final = 2.12098877e+18 +e final = 2.120989822e+18 ------------------------------------ new mass fractions: H1 0.1076931361 -He4 0.08267201555 -C12 1.78524719e-05 -N13 3.393590165e-08 -N14 3.304045126e-07 -O16 0.1086086904 -F18 1.418025583e-07 -Ne20 0.007001757674 -Ne21 5.733676716e-06 +He4 0.08267143674 +C12 1.78525873e-05 +N13 3.39360867e-08 +N14 3.30410921e-07 +O16 0.1086082971 +F18 1.418040871e-07 +Ne20 0.007001727083 +Ne21 5.733694226e-06 Na22 0.1571361579 -Na23 4.819437609e-07 -Mg24 0.004453897399 -Al27 6.192300093e-06 -Si28 0.0614675563 -P31 9.538998065e-06 -S32 0.1152074078 -Ar36 0.09902731773 -Ca40 0.2507582723 -Ti44 0.005256238334 -Cr48 0.0006643598898 -Fe52 1.284019674e-05 -Ni56 4.684206046e-08 +Na23 4.819387697e-07 +Mg24 0.004453876825 +Al27 6.192232453e-06 +Si28 0.06146737281 +P31 9.538910626e-06 +S32 0.1152076517 +Ar36 0.09902790304 +Ca40 0.2507587192 +Ti44 0.005256188598 +Cr48 0.0006643410328 +Fe52 1.283949753e-05 +Ni56 4.683798207e-08 ------------------------------------ species creation rates: -omegadot(H1): 769.3136122 -omegadot(He4): -41732.79844 -omegadot(C12): -9998.214753 +omegadot(H1): 769.3136117 +omegadot(He4): -41732.85633 +omegadot(C12): -9998.214741 omegadot(N13): -9999.996606 -omegadot(N14): -9999.96696 -omegadot(O16): 860.8690421 -omegadot(F18): 0.01418025583 -omegadot(Ne20): 700.1757674 -omegadot(Ne21): 0.5733676716 +omegadot(N14): -9999.966959 +omegadot(O16): 860.8297149 +omegadot(F18): 0.01418040871 +omegadot(Ne20): 700.1727083 +omegadot(Ne21): 0.5733694226 omegadot(Na22): 15713.61579 -omegadot(Na23): 0.04819437609 -omegadot(Mg24): 445.3897399 -omegadot(Al27): 0.6192300093 -omegadot(Si28): 6146.75563 -omegadot(P31): 0.9538998065 -omegadot(S32): 11520.74078 -omegadot(Ar36): 9902.731773 -omegadot(Ca40): 25075.82723 -omegadot(Ti44): 525.6238334 -omegadot(Cr48): 66.43598898 -omegadot(Fe52): 1.284019674 -omegadot(Ni56): 0.004684206046 -number of steps taken: 9600 -AMReX (23.05-4-ga393d7ff7e32) finalized +omegadot(Na23): 0.04819387697 +omegadot(Mg24): 445.3876825 +omegadot(Al27): 0.6192232453 +omegadot(Si28): 6146.737281 +omegadot(P31): 0.9538910626 +omegadot(S32): 11520.76517 +omegadot(Ar36): 9902.790304 +omegadot(Ca40): 25075.87192 +omegadot(Ti44): 525.6188598 +omegadot(Cr48): 66.43410328 +omegadot(Fe52): 1.283949753 +omegadot(Ni56): 0.004683798207 +number of steps taken: 9526 +AMReX (23.10-13-gecaa14456e3f) finalized diff --git a/unit_test/burn_cell/ci-benchmarks/subch_approx_unit_test.out b/unit_test/burn_cell/ci-benchmarks/subch_approx_unit_test.out index ed779a64dc..252b1ebd4c 100644 --- a/unit_test/burn_cell/ci-benchmarks/subch_approx_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/subch_approx_unit_test.out @@ -1,4 +1,4 @@ -AMReX (23.05-4-ga393d7ff7e32) initialized +AMReX (23.10-13-gecaa14456e3f) initialized starting the single zone burn... reading in network electron-capture / beta-decay tables... Maximum Time (s): 1e-05 @@ -51,59 +51,59 @@ RHS at t = 0 Ni56 2.574225196e-29 ------------------------------------ successful? 1 - - Hnuc = 7.242504668e+22 - - added e = 7.242504668e+17 - - final T = 3332738253 + - Hnuc = 7.242795607e+22 + - added e = 7.242795607e+17 + - final T = 3332749137 ------------------------------------ e initial = 1.396711859e+18 -e final = 2.120962326e+18 +e final = 2.12099142e+18 ------------------------------------ new mass fractions: -H1 0.1076931392 -He4 0.08271604755 -C12 1.784659301e-05 -N13 3.392571632e-08 -N14 3.299221207e-07 -O16 0.108645045 -F18 1.417109763e-07 -Ne20 0.007004350626 -Ne21 5.73319595e-06 -Na22 0.1571361633 -Na23 4.823680511e-07 -Mg24 0.004455549374 -Al27 6.197826902e-06 -Si28 0.06147651501 -P31 9.545310692e-06 -S32 0.1151668114 -Ar36 0.09896332458 -Ca40 0.2507603092 -Ti44 0.005262916908 -Cr48 0.0006665524196 -Fe52 1.291730016e-05 -Ni56 4.728298219e-08 +H1 0.1076931417 +He4 0.08271029045 +C12 1.784869828e-05 +N13 3.39288153e-08 +N14 3.299674554e-07 +O16 0.1086428974 +F18 1.41735561e-07 +Ne20 0.007004116589 +Ne21 5.733659437e-06 +Na22 0.1571361627 +Na23 4.823302158e-07 +Mg24 0.00445536025 +Al27 6.197224432e-06 +Si28 0.06147398898 +P31 9.544436584e-06 +S32 0.1151631603 +Ar36 0.0989622079 +Ca40 0.2507752683 +Ti44 0.005263479859 +Cr48 0.0006666468945 +Fe52 1.291933247e-05 +Ni56 4.7290363e-08 ------------------------------------ species creation rates: -omegadot(H1): 769.3139202 -omegadot(He4): -41728.39524 -omegadot(C12): -9998.215341 +omegadot(H1): 769.3141707 +omegadot(He4): -41728.97096 +omegadot(C12): -9998.21513 omegadot(N13): -9999.996607 -omegadot(N14): -9999.967008 -omegadot(O16): 864.5045001 -omegadot(F18): 0.01417109763 -omegadot(Ne20): 700.4350626 -omegadot(Ne21): 0.573319595 -omegadot(Na22): 15713.61633 -omegadot(Na23): 0.04823680511 -omegadot(Mg24): 445.5549374 -omegadot(Al27): 0.6197826902 -omegadot(Si28): 6147.651501 -omegadot(P31): 0.9545310692 -omegadot(S32): 11516.68114 -omegadot(Ar36): 9896.332458 -omegadot(Ca40): 25076.03092 -omegadot(Ti44): 526.2916908 -omegadot(Cr48): 66.65524196 -omegadot(Fe52): 1.291730016 -omegadot(Ni56): 0.004728298219 -number of steps taken: 594 -AMReX (23.05-4-ga393d7ff7e32) finalized +omegadot(N14): -9999.967003 +omegadot(O16): 864.289744 +omegadot(F18): 0.0141735561 +omegadot(Ne20): 700.4116589 +omegadot(Ne21): 0.5733659437 +omegadot(Na22): 15713.61627 +omegadot(Na23): 0.04823302158 +omegadot(Mg24): 445.536025 +omegadot(Al27): 0.6197224432 +omegadot(Si28): 6147.398898 +omegadot(P31): 0.9544436584 +omegadot(S32): 11516.31603 +omegadot(Ar36): 9896.22079 +omegadot(Ca40): 25077.52683 +omegadot(Ti44): 526.3479859 +omegadot(Cr48): 66.66468945 +omegadot(Fe52): 1.291933247 +omegadot(Ni56): 0.0047290363 +number of steps taken: 544 +AMReX (23.10-13-gecaa14456e3f) finalized diff --git a/unit_test/test_neutrino_cooling/ci-benchmarks/test_sneut5.out b/unit_test/test_neutrino_cooling/ci-benchmarks/test_sneut5.out index c314665d07..43e9d72478 100644 --- a/unit_test/test_neutrino_cooling/ci-benchmarks/test_sneut5.out +++ b/unit_test/test_neutrino_cooling/ci-benchmarks/test_sneut5.out @@ -26,6 +26,6 @@ X_proton 0 0.026315789474 sneut 0 3.3206297956e+23 dsneutdt 0 2.9114557633e+14 - dsneutda -2.0028760264e+16 1.6034657416e+20 - dsneutdz -1.832532276e+20 2.2890025857e+16 + dsneutda -1.9892659343e+16 1.6034657636e+20 + dsneutdz -1.8325323011e+20 2.2734481947e+16 From 48bfe6bcd961ca88db56d926d88e1fab5e10753e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 30 Oct 2023 09:03:17 -0400 Subject: [PATCH 27/28] another fix --- .../ci-benchmarks/aprox19_nse_unit_test.out | 88 +++++++++---------- 1 file changed, 44 insertions(+), 44 deletions(-) diff --git a/unit_test/burn_cell/ci-benchmarks/aprox19_nse_unit_test.out b/unit_test/burn_cell/ci-benchmarks/aprox19_nse_unit_test.out index f945932970..44a86cb701 100644 --- a/unit_test/burn_cell/ci-benchmarks/aprox19_nse_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/aprox19_nse_unit_test.out @@ -1,4 +1,4 @@ -AMReX (23.01-12-ga05384c0b7c4) initialized +AMReX (23.10-13-gecaa14456e3f) initialized starting the single zone burn... reading the NSE table (C++) ... Maximum Time (s): 0.01 @@ -43,57 +43,57 @@ RHS at t = 0 Ni56 -4.649371884e-24 n 6.323463259e-06 p 6.323463259e-06 -burn entered NSE during integration (after 280 steps), zone = (0, 0, 0) +burn entered NSE during integration (after 407 steps), zone = (0, 0, 0) recovering burn failure in NSE, zone = (0, 0, 0) ------------------------------------ successful? 1 - - Hnuc = 2.552804555e+21 - - added e = 2.552804555e+19 - - final T = 1.060964458e+10 + - Hnuc = 2.577064219e+21 + - added e = 2.577064219e+19 + - final T = 1.063350508e+10 ------------------------------------ e initial = 1.516256085e+18 -e final = 2.704430164e+19 +e final = 2.728689827e+19 ------------------------------------ new mass fractions: -H1 0.0001057720428 -He3 1.320599831e-10 -He4 0.06412077048 -C12 1.016297715e-11 -N14 9.999991592e-31 -O16 9.999991592e-31 -Ne20 3.890552725e-17 -Mg24 9.999991592e-31 -Si28 1.194568144e-21 -S32 1.956720348e-25 -Ar36 4.189728368e-29 -Ca40 9.999991592e-31 -Ti44 3.531480561e-30 -Cr48 9.999991592e-31 -Fe52 9.999991592e-31 -Fe54 9.999991592e-31 -Ni56 9.999991592e-31 -n 0.4665061227 -p 0.4692673346 +H1 6.505608868e-05 +He3 1.038927825e-10 +He4 0.06022070651 +C12 9.99999158e-31 +N14 1.731352613e-13 +O16 2.450459695e-15 +Ne20 9.99999158e-31 +Mg24 6.790701604e-20 +Si28 3.428830286e-22 +S32 9.99999158e-31 +Ar36 9.99999158e-31 +Ca40 2.030520937e-29 +Ti44 9.99999158e-31 +Cr48 5.120465763e-30 +Fe52 9.99999158e-31 +Fe54 3.385938672e-28 +Ni56 9.99999158e-31 +n 0.4684052388 +p 0.4713089985 ------------------------------------ species creation rates: -omegadot(H1): -9.989422796 -omegadot(He3): -2.499999987 -omegadot(He4): -73.58792295 -omegadot(C12): -2.499999999 +omegadot(H1): -9.993494391 +omegadot(He3): -2.49999999 +omegadot(He4): -73.97792935 +omegadot(C12): -2.5 omegadot(N14): -2.5 omegadot(O16): -2.5 -omegadot(Ne20): 3.890552725e-15 -omegadot(Mg24): -8.408373504e-35 -omegadot(Si28): 1.194568143e-19 -omegadot(S32): 1.956710348e-23 -omegadot(Ar36): 4.089728368e-27 -omegadot(Ca40): -8.408373504e-35 -omegadot(Ti44): 2.531480561e-28 -omegadot(Cr48): -8.408373504e-35 -omegadot(Fe52): -8.408373504e-35 -omegadot(Fe54): -8.408373504e-35 -omegadot(Ni56): -8.408373504e-35 -omegadot(n): 46.65061227 -omegadot(p): 46.92673346 -number of steps taken: 27384 -AMReX (23.01-12-ga05384c0b7c4) finalized +omegadot(Ne20): -8.420134623e-35 +omegadot(Mg24): 6.790701604e-18 +omegadot(Si28): 3.428830276e-20 +omegadot(S32): -8.420134623e-35 +omegadot(Ar36): -8.420134623e-35 +omegadot(Ca40): 1.930520937e-27 +omegadot(Ti44): -8.420134623e-35 +omegadot(Cr48): 4.120465763e-28 +omegadot(Fe52): -8.420134623e-35 +omegadot(Fe54): 3.375938672e-26 +omegadot(Ni56): -8.420134623e-35 +omegadot(n): 46.84052388 +omegadot(p): 47.13089985 +number of steps taken: 22341 +AMReX (23.10-13-gecaa14456e3f) finalized From 368c00c910388c7073219dd85f9e6c1cfb572f4c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 30 Oct 2023 09:07:48 -0400 Subject: [PATCH 28/28] more benchmark updates --- .../aprox19_BE_state_over_time.txt | 22 +++++++++---------- .../ci-benchmarks/aprox19_state_over_time.txt | 12 +++++----- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_BE_state_over_time.txt b/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_BE_state_over_time.txt index 69d51e0bea..c3efe2b44e 100644 --- a/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_BE_state_over_time.txt +++ b/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_BE_state_over_time.txt @@ -1,13 +1,13 @@ # Time Density Temperature H1 He3 He4 C12 N14 O16 Ne20 Mg24 Si28 S32 Ar36 Ca40 Ti44 Cr48 Fe52 Fe54 Ni56 n p 0 1e+06 3e+09 0.7 0.025 0.2 0.025 0.025 0.025 0 0 0 0 0 0 0 0 0 0 0 0 0 - 1e-10 1e+06 3.02856e+09 0.706582 0.00189497 0.215397 0.0182827 0.0328219 0.0249997 2.17185e-05 4.00707e-09 1.44235e-11 1.21931e-15 1.85797e-20 1.27593e-25 1.03595e-30 1e-30 1e-30 1e-30 1e-30 2.32624e-30 2.32624e-30 -3.16228e-10 1e+06 3.03476e+09 0.705509 0.000628581 0.216222 0.00931244 0.0432434 0.024999 8.512e-05 4.85269e-08 4.61301e-11 5.41599e-15 2.49264e-19 5.35054e-24 5.79192e-30 1e-30 1e-30 1e-30 1e-30 5.98461e-30 5.98461e-30 - 1e-09 1e+06 3.03969e+09 0.704286 0.000200506 0.21643 0.00111892 0.0526236 0.0249968 0.000342966 6.23243e-07 5.38214e-10 5.62735e-14 5.68043e-18 3.2007e-22 8.23321e-28 1.00066e-30 1e-30 4.24516e-28 1e-30 1.95328e-29 1.95328e-29 -3.16228e-09 1e+06 3.04074e+09 0.704146 6.42104e-05 0.216256 2.27397e-06 0.0533151 0.0249897 0.00121946 7.46561e-06 1.94541e-08 4.61129e-12 1.05971e-15 1.37192e-19 8.46684e-25 2.71189e-30 1.00001e-30 4.24516e-28 1e-30 6.51863e-29 6.51863e-29 - 1e-08 1e+06 3.04176e+09 0.70416 2.09463e-05 0.215457 2.67096e-11 0.0514278 0.0249673 0.00388704 7.92663e-05 6.95194e-07 5.53831e-10 4.29729e-13 1.88314e-16 3.94438e-21 2.71834e-26 1.10323e-30 4.24516e-28 1e-30 2.12447e-28 2.12447e-28 -3.16228e-08 1e+06 3.04471e+09 0.704165 6.50597e-06 0.212971 2.32879e-11 0.0459107 0.0248981 0.0112961 0.000733137 2.01775e-05 4.919e-08 1.14434e-10 1.47598e-13 8.93693e-18 1.75295e-22 1.86294e-27 4.24521e-28 1e-30 6.99839e-28 6.9984e-28 - 1e-07 1e+06 3.05275e+09 0.704166 2.0132e-06 0.206109 2.10593e-11 0.032257 0.0246919 0.0266086 0.00565053 0.000510727 3.96209e-06 2.88476e-08 1.15228e-10 2.13986e-14 1.28346e-18 4.13423e-23 7.01099e-28 1.05157e-30 2.45097e-27 2.46121e-27 -3.16228e-07 1e+06 3.07032e+09 0.704166 6.23619e-07 0.191035 1.67132e-11 0.0110708 0.024113 0.0345257 0.0260972 0.00875392 0.000232005 5.60651e-06 7.35601e-08 4.44762e-11 8.74683e-15 9.19961e-19 2.00534e-23 3.69734e-27 9.94337e-27 7.52646e-25 - 1e-06 1e+06 3.09565e+09 0.704167 1.92199e-07 0.169122 1.16179e-11 0.000490962 0.022467 0.00934994 0.0309577 0.0569389 0.00595735 0.000524875 2.45272e-05 5.14426e-08 3.52397e-11 1.27949e-14 9.58907e-19 1.73607e-22 4.89224e-26 3.55151e-20 -3.16228e-06 1e+06 3.11504e+09 0.704167 5.937e-08 0.150126 8.16958e-12 7.84999e-08 0.0180619 0.000893683 0.00277257 0.0734901 0.0352149 0.0128016 0.00245239 1.95824e-05 4.98184e-08 6.56227e-11 1.7601e-14 3.14095e-18 3.3984e-25 6.51888e-16 - 1e-05 1e+06 3.14013e+09 0.704167 1.8749e-08 0.120661 4.27519e-12 4.8354e-09 0.00988207 0.000461506 0.000861871 0.0268741 0.0368383 0.0508397 0.0478301 0.00156971 1.54457e-05 7.53675e-08 7.2984e-11 1.28002e-14 2.10934e-21 2.70311e-12 + 1e-10 1e+06 3.02856e+09 0.706582 0.00189497 0.215397 0.0182827 0.0328219 0.0249997 2.17185e-05 4.00707e-09 1.44235e-11 1.21931e-15 1.85797e-20 1.27593e-25 1.02919e-30 1.35535e-26 1.0011e-30 1e-30 1e-30 2.32624e-30 2.32624e-30 +3.16228e-10 1e+06 3.03476e+09 0.705509 0.000628581 0.216222 0.00931244 0.0432434 0.024999 8.512e-05 4.85269e-08 4.61301e-11 5.41599e-15 2.49264e-19 5.35054e-24 5.78516e-30 1.35535e-26 1.00859e-30 1e-30 1e-30 5.98461e-30 5.98461e-30 + 1e-09 1e+06 3.03969e+09 0.704286 0.000200506 0.21643 0.00111892 0.0526236 0.0249968 0.000342966 6.23243e-07 5.38214e-10 5.62735e-14 5.68043e-18 3.2007e-22 8.23314e-28 1.35535e-26 1.03309e-30 1e-30 1e-30 1.95328e-29 1.95328e-29 +3.16228e-09 1e+06 3.04074e+09 0.704146 6.42104e-05 0.216256 2.27397e-06 0.0533151 0.0249897 0.00121946 7.46561e-06 1.94541e-08 4.61129e-12 1.05971e-15 1.37192e-19 8.46684e-25 1.35551e-26 1.1116e-30 1e-30 1e-30 6.51863e-29 6.51863e-29 + 1e-08 1e+06 3.04176e+09 0.70416 2.09463e-05 0.215457 2.67096e-11 0.0514278 0.0249673 0.00388704 7.92663e-05 6.95194e-07 5.53831e-10 4.29729e-13 1.88314e-16 3.94438e-21 4.07356e-26 1.46363e-30 1.00001e-30 1e-30 2.12447e-28 2.12447e-28 +3.16228e-08 1e+06 3.04471e+09 0.704165 6.5068e-06 0.212971 2.32879e-11 0.0459108 0.0248981 0.011296 0.000733175 2.01812e-05 4.92086e-08 1.14507e-10 1.47741e-13 8.94913e-18 1.75628e-22 1.86847e-27 1.0042e-30 1e-30 6.99842e-28 6.99842e-28 + 1e-07 1e+06 3.05275e+09 0.704166 2.01328e-06 0.206109 2.10593e-11 0.0322571 0.0246919 0.0266085 0.00565055 0.000510734 3.96223e-06 2.88493e-08 1.15238e-10 2.14014e-14 1.28369e-18 4.13518e-23 2.77662e-28 1.05158e-30 2.45097e-27 2.46122e-27 +3.16228e-07 1e+06 3.07032e+09 0.704166 6.23626e-07 0.191035 1.67132e-11 0.0110708 0.024113 0.0345257 0.0260972 0.00875393 0.000232006 5.60654e-06 7.35607e-08 4.44768e-11 8.74697e-15 9.19981e-19 2.00535e-23 3.69744e-27 9.94337e-27 7.52667e-25 + 1e-06 1e+06 3.09565e+09 0.704167 1.92198e-07 0.169122 1.16179e-11 0.000490925 0.022467 0.00934976 0.0309578 0.0569392 0.00595731 0.00052486 2.45257e-05 5.14375e-08 3.52346e-11 1.27923e-14 9.58651e-19 1.73561e-22 4.89223e-26 3.55056e-20 +3.16228e-06 1e+06 3.11504e+09 0.704167 5.93699e-08 0.150126 8.16958e-12 7.84947e-08 0.0180619 0.000893682 0.00277256 0.0734901 0.035215 0.0128016 0.00245239 1.95822e-05 4.98179e-08 6.56216e-11 1.76006e-14 3.14089e-18 3.39838e-25 6.51874e-16 + 1e-05 1e+06 3.14013e+09 0.704167 1.8749e-08 0.120661 4.27519e-12 4.8354e-09 0.00988207 0.000461506 0.000861871 0.0268741 0.0368383 0.0508397 0.0478301 0.00156971 1.54457e-05 7.53674e-08 7.29839e-11 1.28002e-14 2.10934e-21 2.70311e-12 diff --git a/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_state_over_time.txt b/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_state_over_time.txt index d7c744215a..39d42da8cc 100644 --- a/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_state_over_time.txt +++ b/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_state_over_time.txt @@ -4,10 +4,10 @@ 3.16228e-10 1e+06 3.03478e+09 0.705507 0.000621921 0.216227 0.0092854 0.043275 0.024999 8.50466e-05 4.81558e-08 4.59635e-11 5.40157e-15 2.46986e-19 5.24295e-24 5.63193e-30 1e-30 1e-30 1e-30 1e-30 5.98338e-30 5.98338e-30 1e-09 1e+06 3.03971e+09 0.704282 0.000198228 0.216432 0.00109041 0.0526569 0.0249968 0.000342898 6.19859e-07 5.31673e-10 5.54672e-14 5.58299e-18 3.12982e-22 7.99114e-28 1.00063e-30 1e-30 1e-30 1e-30 1.95324e-29 1.95324e-29 3.16228e-09 1e+06 3.04074e+09 0.704146 6.28602e-05 0.216257 1.73328e-06 0.0533156 0.0249897 0.00121971 7.35522e-06 1.86066e-08 4.22894e-12 9.22256e-16 1.12401e-19 6.48983e-25 2.22129e-30 1.00001e-30 1e-30 9.99999e-31 6.51843e-29 6.51843e-29 - 1e-08 1e+06 3.04176e+09 0.70416 2.04828e-05 0.215457 0 0.0514256 0.0249673 0.00389241 7.69853e-05 6.38368e-07 4.73546e-10 3.40222e-13 1.3847e-16 2.72235e-21 1.7909e-26 1.06638e-30 1.00001e-30 9.99996e-31 2.12345e-28 2.12345e-28 -3.16228e-08 1e+06 3.04471e+09 0.704164 6.34237e-06 0.212969 2.41781e-11 0.0459028 0.0248981 0.0113145 0.000725596 1.98421e-05 4.89394e-08 1.18275e-10 1.63701e-13 1.10302e-17 2.50276e-22 3.20374e-27 1.00897e-30 9.99991e-31 6.99475e-28 6.99475e-28 - 1e-07 1e+06 3.05276e+09 0.704166 1.97272e-06 0.206101 2.10567e-11 0.0322273 0.0246918 0.0266706 0.00563202 0.00050569 3.92045e-06 2.8844e-08 1.18231e-10 2.29592e-14 1.4715e-18 5.18729e-23 3.90644e-28 1.07265e-30 2.4494e-27 2.46383e-27 -3.16228e-07 1e+06 3.07034e+09 0.704166 6.16355e-07 0.191015 1.67081e-11 0.0110228 0.0241129 0.0345787 0.0261422 0.00872611 0.000229534 5.49795e-06 7.15038e-08 4.2921e-11 8.40418e-15 8.83784e-19 1.93663e-23 3.5714e-27 9.93669e-27 7.27207e-25 + 1e-08 1e+06 3.04176e+09 0.70416 2.04835e-05 0.215457 0 0.0514256 0.0249673 0.00389241 7.69853e-05 6.38372e-07 4.73546e-10 3.40208e-13 1.3845e-16 2.72143e-21 1.78973e-26 1.0663e-30 1.00001e-30 9.99996e-31 2.12345e-28 2.12345e-28 +3.16228e-08 1e+06 3.04471e+09 0.704164 6.34257e-06 0.212969 2.41718e-11 0.0459028 0.0248981 0.0113145 0.000725593 1.98427e-05 4.8946e-08 1.18312e-10 1.63798e-13 1.10409e-17 2.50638e-22 3.21027e-27 1.009e-30 9.99991e-31 6.99482e-28 6.99482e-28 + 1e-07 1e+06 3.05276e+09 0.704166 1.97274e-06 0.206101 2.10567e-11 0.0322273 0.0246918 0.0266706 0.00563202 0.000505691 3.92047e-06 2.88443e-08 1.18234e-10 2.29599e-14 1.47157e-18 5.18764e-23 3.90679e-28 1.07266e-30 2.4494e-27 2.46384e-27 +3.16228e-07 1e+06 3.07034e+09 0.704166 6.16357e-07 0.191015 1.67081e-11 0.0110228 0.0241129 0.0345787 0.0261422 0.00872611 0.000229534 5.49795e-06 7.15038e-08 4.2921e-11 8.4042e-15 8.83786e-19 1.93664e-23 3.57141e-27 9.93669e-27 7.2721e-25 1e-06 1e+06 3.09568e+09 0.704166 1.90977e-07 0.169097 1.16129e-11 0.000479845 0.0224668 0.00930122 0.0309807 0.0570159 0.00594684 0.000520989 2.41528e-05 5.01772e-08 3.40108e-11 1.22121e-14 9.05195e-19 1.63882e-22 4.89062e-26 3.35258e-20 -3.16228e-06 1e+06 3.11505e+09 0.704166 5.90985e-08 0.150118 8.16827e-12 5.39796e-08 0.0180596 0.000891114 0.00273936 0.0735111 0.0352666 0.0127917 0.00243636 1.92951e-05 4.85823e-08 6.32194e-11 1.6726e-14 2.98487e-18 3.33563e-25 6.19482e-16 - 1e-05 1e+06 3.14015e+09 0.704166 1.8682e-08 0.120635 4.27245e-12 4.83365e-09 0.00987352 0.000461095 0.00086113 0.0268146 0.036848 0.0509385 0.0478247 0.00156184 1.52568e-05 7.37586e-08 7.06632e-11 1.23934e-14 2.01895e-21 2.61715e-12 +3.16228e-06 1e+06 3.11505e+09 0.704166 5.90985e-08 0.150118 8.16827e-12 5.39787e-08 0.0180596 0.000891114 0.00273936 0.0735111 0.0352666 0.0127917 0.00243636 1.92951e-05 4.85823e-08 6.32193e-11 1.6726e-14 2.98487e-18 3.33563e-25 6.19482e-16 + 1e-05 1e+06 3.14015e+09 0.704166 1.8682e-08 0.120635 4.27245e-12 4.83365e-09 0.00987352 0.000461095 0.00086113 0.0268146 0.036848 0.0509385 0.0478247 0.00156184 1.52568e-05 7.37585e-08 7.06631e-11 1.23934e-14 2.01895e-21 2.61715e-12