From 8b48c0281d9f18eac6a6ed167871bbb2bd433d37 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 13:09:04 -0400 Subject: [PATCH 1/3] move some Fermi integral coefficients to static arrays (#1371) no need to set these each time we enter the function --- neutrinos/sneut5.H | 176 ++++++++++++++++++++++----------------------- 1 file changed, 88 insertions(+), 88 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 9348f75f86..b5af627dcd 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -15,49 +15,48 @@ 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; - Array1D a1; - Array1D b1; - Array1D a2; - Array1D b2; + 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; - - 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; + 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}; + + const Array1D b1 = { + 1.771804140488e4_rt, + -2.014785161019e3_rt, + 9.130355392717e1_rt, + -1.670718177489e0_rt}; + + const 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}; + + 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}; if (f < 4.0e0_rt) { @@ -125,64 +124,65 @@ Real zfermim12(const Real x) // reference: antia apjs 84,101 1993 // declare work variables - int m1,k1,m2,k2; Real rn,den,xx; - Array1D a1,b1; - Array1D a2,b2; // return value Real zfermim12r; // load the coefficients of the expansion from Table 2 of Antia - m1 = 7; - k1 = 7; - 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; + constexpr int m1{7}; + constexpr int k1{7}; + 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) { From e052f6eaaef09bd58ff25540d8fcc814a9087eae Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 13:44:53 -0400 Subject: [PATCH 2/3] move common sneut factors into a struct (#1372) --- neutrinos/sneut5.H | 361 +++++++++++++++++++++++++-------------------- 1 file changed, 203 insertions(+), 158 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index b5af627dcd..7f676b8525 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -223,6 +223,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 @@ -249,22 +327,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, 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 @@ -301,7 +378,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; @@ -357,48 +433,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; - - // some composition variables - ye = zbar*abari; + auto sf = get_sneut_factors(den, temp, abar, zbar); - // 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; @@ -408,16 +453,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; @@ -438,22 +483,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; } @@ -467,24 +512,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; @@ -495,11 +540,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; @@ -528,22 +573,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); @@ -580,20 +625,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 @@ -673,12 +718,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; @@ -777,7 +822,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; @@ -818,7 +863,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); @@ -886,12 +931,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; } @@ -904,24 +949,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; @@ -936,19 +981,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); @@ -992,7 +1037,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) { @@ -1004,11 +1049,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; @@ -1055,13 +1100,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; @@ -1084,12 +1129,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; @@ -1108,11 +1153,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; @@ -1129,7 +1174,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); @@ -1204,11 +1249,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; @@ -1241,11 +1286,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; @@ -1260,11 +1305,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 @@ -1306,11 +1351,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); @@ -1353,50 +1398,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 6f60462ff9ea041d44e24ec9d09a813d76b9d282 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 13:48:16 -0400 Subject: [PATCH 3/3] SDC+NSE predictor-corrector update: make number of iters a runtime param (#1370) --- integration/_parameters | 3 +++ integration/nse_update_simplified_sdc.H | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/integration/_parameters b/integration/_parameters index 998ffa47a5..1579602917 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -77,3 +77,6 @@ sdc_burn_tol_factor real 1.d0 # for Strang, this simply means scaling e by the initial energy? scale_system integer 0 +# for the NSE update predictor-corrector, how many iterations +# do we take to get the new time NSE state +nse_iters integer 3 diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 533eb5171b..fe5fb20548 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -93,7 +93,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { rho_aux_new[n] = state.y[SFX+n] + dt * state.ydot_a[SFX+n] + dt * aux_source[n]; } - for (int iter = 0; iter < 3; iter++) { + for (int iter = 0; iter < integrator_rp::nse_iters; iter++) { // update (rho e)^{n+1} based on the new energy generation rate rhoe_new = state.y[SEINT] + dt * state.ydot_a[SEINT] + dt * rho_enucdot; @@ -258,7 +258,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { rhoX_new[n] = state.y[SFS+n] + dt * state.ydot_a[SFS+n] + dt * rhoX_source[n]; } - for (int iter = 0; iter < 3; iter++) { + for (int iter = 0; iter < integrator_rp::nse_iters; iter++) { // update (rho e)^{n+1} based on the new energy generation rate rhoe_new = state.y[SEINT] + dt * state.ydot_a[SEINT] + dt * rho_enucdot;