From c8890c9e7de85de8dc56aecac26ad8adc1347cd6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 12:22:16 -0400 Subject: [PATCH 1/2] 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 2/2] 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)