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; diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 2f69143dda..2734855984 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -38,22 +38,42 @@ namespace nu_constants { } -namespace ifermi12_coeffs { +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}; - HIP_CONSTEXPR static AMREX_GPU_MANAGED Array1D a1 = { + 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, @@ -62,87 +82,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) { @@ -157,18 +103,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 @@ -182,16 +128,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); @@ -222,19 +168,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; @@ -242,16 +236,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;