diff --git a/networks/CNO_extras/actual_rhs.H b/networks/CNO_extras/actual_rhs.H index 79e1367582..03501b4031 100644 --- a/networks/CNO_extras/actual_rhs.H +++ b/networks/CNO_extras/actual_rhs.H @@ -794,8 +794,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1197,7 +1197,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H index b7aa913121..0fb4bc3422 100644 --- a/networks/ECSN/actual_rhs.H +++ b/networks/ECSN/actual_rhs.H @@ -412,8 +412,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); // Append the energy equation (this is erg/g/s) @@ -629,7 +629,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz); diff --git a/networks/ase/actual_rhs.H b/networks/ase/actual_rhs.H index 26445601af..242edac898 100644 --- a/networks/ase/actual_rhs.H +++ b/networks/ase/actual_rhs.H @@ -1093,8 +1093,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1646,7 +1646,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H index 6a84156b9a..5422a92f35 100644 --- a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H +++ b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H @@ -212,8 +212,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -309,7 +309,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H index cc76ef7078..be1b24da39 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_rhs.H +++ b/networks/ignition_reaclib/URCA-simple/actual_rhs.H @@ -234,8 +234,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -343,7 +343,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/nova/actual_rhs.H b/networks/nova/actual_rhs.H index a9e79d1fcc..260a47e241 100644 --- a/networks/nova/actual_rhs.H +++ b/networks/nova/actual_rhs.H @@ -487,8 +487,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -746,7 +746,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/nova2/actual_rhs.H b/networks/nova2/actual_rhs.H index 49ea50923f..9cae820c4b 100644 --- a/networks/nova2/actual_rhs.H +++ b/networks/nova2/actual_rhs.H @@ -727,8 +727,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1058,7 +1058,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/partition_test/actual_rhs.H b/networks/partition_test/actual_rhs.H index 37d0cac6c0..cbff936ed8 100644 --- a/networks/partition_test/actual_rhs.H +++ b/networks/partition_test/actual_rhs.H @@ -202,8 +202,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -335,7 +335,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/rhs.H b/networks/rhs.H index 55a0060be2..e507d786aa 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -1348,9 +1348,10 @@ void rhs (burn_t& burn_state, Array1D& ydot) // Evaluate the neutrino cooling. #ifdef NEUTRINOS + constexpr int do_derivatives{0}; Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, - sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, + sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); #else Real sneut = 0.0; #endif @@ -1495,8 +1496,9 @@ void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) // Evaluate the neutrino cooling. #ifdef NEUTRINOS Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, - sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + constexpr int do_derivatives{1}; + sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, + sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); #else Real sneut = 0.0, dsneutdt = 0.0, dsneutdd = 0.0, dsnuda = 0.0, dsnudz = 0.0; amrex::ignore_unused(sneut, dsneutdd); diff --git a/networks/sn160/actual_rhs.H b/networks/sn160/actual_rhs.H index 8f81f200bb..647a980a70 100644 --- a/networks/sn160/actual_rhs.H +++ b/networks/sn160/actual_rhs.H @@ -10396,8 +10396,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -18028,7 +18028,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/subch_approx/actual_rhs.H b/networks/subch_approx/actual_rhs.H index a9b90e7e67..c4e168ec2c 100644 --- a/networks/subch_approx/actual_rhs.H +++ b/networks/subch_approx/actual_rhs.H @@ -1375,8 +1375,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1928,7 +1928,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/subch_full/actual_rhs.H b/networks/subch_full/actual_rhs.H index d082bd48de..c27ce0b3d7 100644 --- a/networks/subch_full/actual_rhs.H +++ b/networks/subch_full/actual_rhs.H @@ -1603,8 +1603,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -2390,7 +2390,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/subch_simple/actual_rhs.H b/networks/subch_simple/actual_rhs.H index 68e84fa45f..eae7f7bbbb 100644 --- a/networks/subch_simple/actual_rhs.H +++ b/networks/subch_simple/actual_rhs.H @@ -1227,8 +1227,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1723,7 +1723,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 91093437f9..9348f75f86 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -223,6 +223,8 @@ Real zfermim12(const Real x) return zfermim12r; } + +template AMREX_GPU_HOST_DEVICE AMREX_INLINE void sneut5(const Real temp, const Real den, const Real abar, const Real zbar, @@ -252,8 +254,8 @@ 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,f0, - f1,deni,tempi,abari,zbari,f2,f3,z,ye, + 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, dum,dumdt,gum,gumdt,dumda,dumdz,gumda,gumdz; // pair production @@ -387,15 +389,18 @@ void sneut5(const Real temp, const Real den, rmda = -rm*abari; rmdz = den*abari; rmi = 1.0e0_rt/rm; -; + a0 = rm * 1.0e-9_rt; a1 = std::pow(a0, oneth); zeta = a1 * xlm1; - zetadt = -a1 * xlm2 * xldt; - a2 = oneth * a1*rmi * xlm1; - zetada = a2 * rmda; - zetadz = a2 * rmdz; -; + + if constexpr (do_derivatives) { + a2 = oneth * a1*rmi * xlm1; + zetadt = -a1 * xlm2 * xldt; + zetada = a2 * rmda; + zetadz = a2 * rmdz; + } + zeta2 = zeta * zeta; zeta3 = zeta2 * zeta; @@ -404,95 +409,120 @@ void sneut5(const Real temp, const Real den, // equation 2.8 gl = 1.0e0_rt - 13.04e0_rt*xl2 +133.5e0_rt*xl4 +1534.0e0_rt*xl6 +918.6e0_rt*xl8; - gldt = xldt*(-26.08e0_rt*xl +534.0e0_rt*xl3 +9204.0e0_rt*xl5 +7348.8e0_rt*xl7); + if constexpr (do_derivatives) { + gldt = xldt*(-26.08e0_rt*xl +534.0e0_rt*xl3 +9204.0e0_rt*xl5 +7348.8e0_rt*xl7); + } // equation 2.7 a1 = 6.002e19_rt + 2.084e20_rt*zeta + 1.872e21_rt*zeta2; - a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta; if (t9 < 10.0_rt) { b1 = std::exp(-5.5924e0_rt*zeta); - b2 = -b1*5.5924e0_rt; + if constexpr (do_derivatives) { + b2 = -b1*5.5924e0_rt; + } } else { b1 = std::exp(-4.9924e0_rt*zeta); - b2 = -b1*4.9924e0_rt; + if constexpr (do_derivatives) { + b2 = -b1*4.9924e0_rt; + } } xnum = a1 * b1; - c = a2*b1 + a1*b2; - xnumdt = c*zetadt; - xnumda = c*zetada; - xnumdz = c*zetadz; + if constexpr (do_derivatives) { + a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta; + c = a2*b1 + a1*b2; + xnumdt = c*zetadt; + xnumda = c*zetada; + xnumdz = c*zetadz; + } if (t9 < 10.0_rt) { a1 = 9.383e-1_rt*xlm1 - 4.141e-1_rt*xlm2 + 5.829e-2_rt*xlm3; - a2 = -9.383e-1_rt*xlm2 + 2.0e0_rt*4.141e-1_rt*xlm3 - 3.0e0_rt*5.829e-2_rt*xlm4; + 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; + } } else { a1 = 1.2383e0_rt*xlm1 - 8.141e-1_rt*xlm2; - a2 = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3; + if constexpr (do_derivatives) { + a2 = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3; + } } - b1 = 3.0e0_rt*zeta2; - xden = zeta3 + a1; - xdendt = b1*zetadt + a2*xldt; - xdenda = b1*zetada; - xdendz = b1*zetadz; + if constexpr (do_derivatives) { + b1 = 3.0e0_rt*zeta2; + xdendt = b1*zetadt + a2*xldt; + xdenda = b1*zetada; + xdendz = b1*zetadz; + } a1 = 1.0e0_rt/xden; fpair = xnum*a1; - fpairdt = (xnumdt - fpair*xdendt)*a1; - fpairda = (xnumda - fpair*xdenda)*a1; - fpairdz = (xnumdz - fpair*xdendz)*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*xl2 + 0.3967e0_rt*xlp5 + 1.005e0_rt; - a2 = xldt*(2.0e0_rt*10.7480e0_rt*xl + 0.5e0_rt*0.3967e0_rt*xlmp5); xnum = 1.0e0_rt/a1; - xnumdt = -xnum*xnum*a2; + if constexpr (do_derivatives) { + a2 = xldt*(2.0e0_rt*10.7480e0_rt*xl + 0.5e0_rt*0.3967e0_rt*xlmp5); + xnumdt = -xnum*xnum*a2; + } a1 = 7.692e7_rt*xl3 + 9.715e6_rt*xlp5; - a2 = xldt*(3.0e0_rt*7.692e7_rt*xl2 + 0.5e0_rt*9.715e6_rt*xlmp5); - c = 1.0e0_rt/a1; b1 = 1.0e0_rt + rm*c; xden = std::pow(b1, -0.3e0_rt); - - d = -0.3e0_rt*xden/b1; - xdendt = -d*rm*c*c*a2; - xdenda = d*rmda*c; - xdendz = d*rmdz*c; + if constexpr (do_derivatives) { + a2 = xldt*(3.0e0_rt*7.692e7_rt*xl2 + 0.5e0_rt*9.715e6_rt*xlmp5); + d = -0.3e0_rt*xden/b1; + xdendt = -d*rm*c*c*a2; + xdenda = d*rmda*c; + xdendz = d*rmdz*c; + } qpair = xnum*xden; - qpairdt = xnumdt*xden + xnum*xdendt; - qpairda = xnum*xdenda; - qpairdz = xnum*xdendz; + if constexpr (do_derivatives) { + qpairdt = xnumdt*xden + xnum*xdendt; + qpairda = xnum*xdenda; + qpairdz = xnum*xdendz; + } // equation 2.5 a1 = std::exp(-2.0e0_rt*xlm1); - a2 = a1*2.0e0_rt*xlm2*xldt; spair = a1*fpair; - spairdt = a2*fpair + a1*fpairdt; - spairda = a1*fpairda; - spairdz = a1*fpairdz; + if constexpr (do_derivatives) { + a2 = a1*2.0e0_rt*xlm2*xldt; + spairdt = a2*fpair + a1*fpairdt; + spairda = a1*fpairda; + spairdz = a1*fpairdz; + } a1 = spair; spair = gl*a1; - spairdt = gl*spairdt + gldt*a1; - spairda = gl*spairda; - spairdz = gl*spairdz; + if constexpr (do_derivatives) { + spairdt = gl*spairdt + gldt*a1; + spairda = gl*spairda; + spairdz = gl*spairdz; + } a1 = tfac4*(1.0e0_rt + tfac3 * qpair); - a2 = tfac4*tfac3; a3 = spair; spair = a1*a3; - spairdt = a1*spairdt + a2*qpairdt*a3; - spairda = a1*spairda + a2*qpairda*a3; - spairdz = a1*spairdz + a2*qpairdz*a3; + if constexpr (do_derivatives) { + a2 = tfac4*tfac3; + spairdt = a1*spairdt + a2*qpairdt*a3; + spairda = a1*spairda + a2*qpairda*a3; + spairdz = a1*spairdz + a2*qpairdz*a3; + } // plasma neutrino section // for collective reactions like gamma_plasmon => nu_e + nubar_e @@ -500,19 +530,21 @@ void sneut5(const Real temp, const Real den, a1 = 1.019e-6_rt*rm; a2 = std::pow(a1, twoth); - a3 = twoth*a2/a1; b1 = std::sqrt(1.0e0_rt + a2); - b2 = 1.0e0_rt/b1; c00 = 1.0e0_rt/(temp*temp*b1); gl2 = 1.1095e11_rt * rm * c00; - 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); + 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); + } gl = std::sqrt(gl2); gl12 = std::sqrt(gl); @@ -523,40 +555,46 @@ void sneut5(const Real temp, const Real den, // equation 4.7 ft = 2.4e0_rt + 0.6e0_rt*gl12 + 0.51e0_rt*gl + 1.25e0_rt*gl32; gum = 1.0e0_rt/gl2; - 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; + 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; - a2 = 8.6e0_rt + 1.75e0_rt*1.35e0_rt*gl72*gum; b1 = 225.0e0_rt - 17.0e0_rt*gl + gl2; - b2 = -0.5e0_rt*17.0e0_rt*gl*gum + 1.0e0_rt; c = 1.0e0_rt/b1; fl = a1*c; - d = (a2 - fl*b2)*c; - fldt = d*gl2dt; - flda = d*gl2da; - fldz = d*gl2dz; + 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*rm); xlnt = std::log10(temp); xnum = sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt); - xnumdt = -iln10*0.5e0_rt*tempi; - a2 = iln10*sixth*rmi; - xnumda = a2*rmda; - xnumdz = a2*rmdz; - xden = sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt); - xdendt = iln10*0.5e0_rt*tempi; - xdenda = a2*rmda; - xdendz = a2*rmdz; + 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; + } // equation 4.11 if (std::abs(xnum) > 0.7e0_rt || xden < 0.0e0_rt) { @@ -569,70 +607,82 @@ void sneut5(const Real temp, const Real den, } else { a1 = 0.39e0_rt - 1.25e0_rt*xnum - 0.35e0_rt*std::sin(4.5e0_rt*xnum); - a2 = -1.25e0_rt - 4.5e0_rt*0.35e0_rt*std::cos(4.5e0_rt*xnum); b1 = 0.3e0_rt * std::exp(-std::pow(4.5e0_rt*xnum + 0.9e0_rt, 2.0_rt)); - b2 = -b1*2.0e0_rt*(4.5e0_rt*xnum + 0.9e0_rt)*4.5e0_rt; c = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt*xnum); - 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; + + 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); - 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); - fxy = 1.05e0_rt + (a1 - b1)*c00; - 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; - + 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; - splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt; - splasda = (ftda + flda)*fxy + (ft+fl)*fxyda; - splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz; + 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); - a3 = -0.5e0_rt*a2*gl*gum; a1 = splas; splas = a2*a1; - splasdt = a2*splasdt + a3*gl2dt*a1; - splasda = a2*splasda + a3*gl2da*a1; - splasdz = a2*splasdz + a3*gl2dz*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; - a3 = 3.0e0_rt*gl6*gum; a1 = splas; splas = a2*a1; - splasdt = a2*splasdt + a3*gl2dt*a1; - splasda = a2*splasda + a3*gl2da*a1; - splasdz = a2*splasdz + a3*gl2dz*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 * xl9; - a3 = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt; a1 = splas; splas = a2*a1; - splasdt = a2*splasdt + a3*a1; - splasda = a2*splasda; - splasdz = a2*splasdz; + if constexpr (do_derivatives) { + a3 = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt; + splasdt = a2*splasdt + a3*a1; + splasda = a2*splasda; + splasdz = a2*splasdz; + } // photoneutrino process section // for reactions like e- + gamma => e- + nu_e + nubar_e @@ -789,108 +839,136 @@ void sneut5(const Real temp, const Real den, + 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; - - f0 = taudt*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; - 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; - f1 = taudt*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; 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; - f2 = taudt*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; + if constexpr (do_derivatives) { + f0 = taudt*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; + + f1 = taudt*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; + + f2 = taudt*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; + } + // equation 3.4 dum = a0 + a1*zeta + a2*zeta2; - 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; + 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; - xnumdt = dumdt*z - dum*z*cc*zetadt; - xnumda = dumda*z - dum*z*cc*zetada; - xnumdz = dumdz*z - dum*z*cc*zetadz; + 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*xlm1 + 7.483e-3_rt*xlm2 + 3.061e-4_rt*xlm3; dum = 3.0e0_rt*zeta2; - 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); - xdenda = dum*zetada; - xdendz = dum*zetadz; - + 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); + xdenda = dum*zetada; + xdendz = dum*zetadz; + } dum = 1.0e0_rt/xden; fphot = xnum*dum; - fphotdt = (xnumdt - fphot*xdendt)*dum; - fphotda = (xnumda - fphot*xdenda)*dum; - fphotdz = (xnumdz - fphot*xdendz)*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 * xl; xnum = 0.666e0_rt*std::pow(a0, -2.066e0_rt); - xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*xldt; + if constexpr (do_derivatives) { + xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*xldt; + } dum = 1.875e8_rt*xl + 1.653e8_rt*xl2 + 8.499e8_rt*xl3 - 1.604e8_rt*xl4; - 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); + 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); + } z = 1.0e0_rt/dum; xden = 1.0e0_rt + rm*z; - xdendt = -rm*z*z*dumdt; - xdenda = rmda*z; - xdendz = rmdz*z; + if constexpr (do_derivatives) { + xdendt = -rm*z*z*dumdt; + xdenda = rmda*z; + xdendz = rmdz*z; + } z = 1.0e0_rt/xden; qphot = xnum*z; - qphotdt = (xnumdt - qphot*xdendt)*z; + if constexpr (do_derivatives) { + qphotdt = (xnumdt - qphot*xdendt)*z; + } dum = -qphot*z; - qphotda = dum*xdenda; - qphotdz = dum*xdendz; + if constexpr (do_derivatives) { + qphotda = dum*xdenda; + qphotdz = dum*xdendz; + } // equation 3.2 sphot = xl5 * fphot; - sphotdt = 5.0e0_rt*xl4*xldt*fphot + xl5*fphotdt; - sphotda = xl5*fphotda; - sphotdz = xl5*fphotdz; + if constexpr (do_derivatives) { + sphotdt = 5.0e0_rt*xl4*xldt*fphot + xl5*fphotdt; + sphotda = xl5*fphotda; + sphotdz = xl5*fphotdz; + } a1 = sphot; sphot = rm*a1; - sphotdt = rm*sphotdt; - sphotda = rm*sphotda + rmda*a1; - sphotdz = rm*sphotdz + rmdz*a1; + if constexpr (do_derivatives) { + sphotdt = rm*sphotdt; + sphotda = rm*sphotda + rmda*a1; + sphotdz = rm*sphotdz + rmdz*a1; + } a1 = tfac4*(1.0e0_rt - tfac3 * qphot); - a2 = -tfac4*tfac3; a3 = sphot; sphot = a1*a3; - sphotdt = a1*sphotdt + a2*qphotdt*a3; - sphotda = a1*sphotda + a2*qphotda*a3; - sphotdz = a1*sphotdz + a2*qphotdz*a3; + if constexpr (do_derivatives) { + a2 = -tfac4*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; - sphotdt = 0.0e0_rt; - sphotda = 0.0e0_rt; - sphotdz = 0.0e0_rt; + if constexpr (do_derivatives) { + sphotdt = 0.0e0_rt; + sphotda = 0.0e0_rt; + sphotdz = 0.0e0_rt; + } } // bremsstrahlung neutrino section @@ -921,13 +999,17 @@ void sneut5(const Real temp, const Real den, // equation 5.3 dum = 7.05e6_rt * t832 + 5.12e4_rt * t83; - dumdt = (1.5e0_rt*7.05e6_rt*t812 + 3.0e0_rt*5.12e4_rt*t82)*1.0e-8_rt; + 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 = rm*z; - etadt = -rm*z*z*dumdt; - etada = rmda*z; - etadz = rmdz*z; + if constexpr (do_derivatives) { + etadt = -rm*z*z*dumdt; + etada = rmda*z; + etadz = rmdz*z; + } etam1 = 1.0e0_rt/eta; etam2 = etam1 * etam1; @@ -939,27 +1021,35 @@ void sneut5(const Real temp, const Real den, xnum = 1.0e0_rt/a0; dum = 1.0e0_rt + 1.47e0_rt*etam1 + 3.29e-2_rt*etam2; - z = -1.47e0_rt*etam2 - 2.0e0_rt*3.29e-2_rt*etam3; - dumdt = z*etadt; - dumda = z*etada; - dumdz = z*etadz; + 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); - z = -1.26e0_rt*etam2; - c01 = z*etadt; - c03 = z*etada; - c04 = z*etadz; + 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; - xdendt = (c01 - xden*dumdt)*z; - xdenda = (c03 - xden*dumda)*z; - xdendz = (c04 - xden*dumdz)*z; + if constexpr (do_derivatives) { + xdendt = (c01 - xden*dumdt)*z; + xdenda = (c03 - xden*dumda)*z; + xdendz = (c04 - xden*dumdz)*z; + } fbrem = xnum + xden; - fbremdt = -xnum*xnum*f0 + xdendt; - fbremda = xdenda; - fbremdz = xdendz; + 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; @@ -967,55 +1057,71 @@ void sneut5(const Real temp, const Real den, z = 1.0e0_rt + rm*1.0e-9_rt; dum = a0*z; - dumdt = f0*z; - z = a0*1.0e-9_rt; - dumda = z*rmda; - dumdz = z*rmdz; + if constexpr (do_derivatives) { + dumdt = f0*z; + z = a0*1.0e-9_rt; + dumda = z*rmda; + dumdz = z*rmdz; + } xnum = 1.0e0_rt/dum; - z = -xnum*xnum; - xnumdt = z*dumdt; - xnumda = z*dumda; - xnumdz = z*dumdz; + 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); - dd00 = (1.5e0_rt*7.75e5_rt*t812 + 3.85e0_rt*247.0e0_rt*std::pow(t8, 2.85e0_rt))*1.0e-8_rt; - c01 = 4.07e0_rt + 0.0240e0_rt * std::pow(t8, 1.4e0_rt); - dd01 = 1.4e0_rt*0.0240e0_rt * std::pow(t8, 0.4e0_rt)*1.0e-8_rt; - c02 = 4.59e-5_rt * std::pow(t8, -0.110e0_rt); - dd02 = -0.11e0_rt*4.59e-5_rt * std::pow(t8, -1.11e0_rt)*1.0e-8_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*rmi + c01 + c02*z; - dumdt = dd00*rmi + dd01 + dd02*z; - z = -c00*rmi*rmi; - dumda = z*rmda; - dumdz = z*rmdz; + if constexpr (do_derivatives) { + dumdt = dd00*rmi + dd01 + dd02*z; + z = -c00*rmi*rmi; + dumda = z*rmda; + dumdz = z*rmdz; + } xden = 1.0e0_rt/dum; - z = -xden*xden; - xdendt = z*dumdt; - xdenda = z*dumda; - xdendz = z*dumdz; + if constexpr (do_derivatives) { + z = -xden*xden; + xdendt = z*dumdt; + xdenda = z*dumda; + xdendz = z*dumdz; + } gbrem = xnum + xden; - gbremdt = xnumdt + xdendt; - gbremda = xnumda + xdenda; - gbremdz = xnumdz + xdendz; + if constexpr (do_derivatives) { + gbremdt = xnumdt + xdendt; + gbremda = xnumda + xdenda; + gbremdz = xnumdz + xdendz; + } // equation 5.1 dum = 0.5738e0_rt*zbar*ye*t86*den; - 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; + 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; + } z = tfac4*fbrem - tfac5*gbrem; sbrem = dum * z; - sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt); - sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda); - sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz); + 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); + } // liquid metal with c12 parameters (not too different for other elements) // equation 5.18 and 5.16 @@ -1099,9 +1205,11 @@ void sneut5(const Real temp, const Real den, // - 0.00069e0_rt*sin5*5.0e0_rt); dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*abari, oneth); - dumdt = -dum*tempi; - dumda = -oneth*dum*abari; - dumdz = 2.0e0_rt*dum*zbari; + if constexpr (do_derivatives) { + dumdt = -dum*tempi; + dumda = -oneth*dum*abari; + dumdz = 2.0e0_rt*dum*zbari; + } gm1 = 1.0e0_rt/dum; gm2 = gm1*gm1; @@ -1119,46 +1227,56 @@ void sneut5(const Real temp, const Real den, // equation 5.19 and 5.20 fliq = v*fb + (1.0e0_rt - v)*ft; - fliqdt = a0*dumdt*(fb - ft); - fliqda = a0*dumda*(fb - ft); - fliqdz = a0*dumdz*(fb - 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; - gliqdt = a1*dumdt*(gb - gt); - gliqda = a1*dumda*(gb - gt); - gliqdz = a1*dumdz*(gb - 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*ye*t86*den; - 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; + 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; + } z = tfac4*fliq - tfac5*gliq; sbrem = dum * z; - sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt); - sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda); - sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz); - + 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); + } } // 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)); - xnumdt = -1.50e0_rt*xnum*tempi; - xnumda = -xnum*abari; - xnumdz = xnum*zbari; + if constexpr (do_derivatives) { + xnumdt = -1.50e0_rt*xnum*tempi; + xnumda = -xnum*abari; + xnumdz = xnum*zbari; + } // the chemical potential nu = ifermi12(xnum); // a0 is d(nu)/d(xnum) a0 = 1.0e0_rt/(0.5e0_rt*zfermim12(nu)); - nudt = a0*xnumdt; - nuda = a0*xnumda; - nudz = a0*xnumdz; - + if constexpr (do_derivatives) { + nudt = a0*xnumdt; + nuda = a0*xnumda; + nudz = a0*xnumdz; + } nu2 = nu * nu; nu3 = nu2 * nu; @@ -1189,16 +1307,20 @@ void sneut5(const Real temp, const Real den, if (nu >= -20.0_rt && nu <= 10.0_rt) { zeta = 1.579e5_rt*zbar*zbar*tempi; - zetadt = -zeta*tempi; - zetada = 0.0e0_rt; - zetadz = 2.0e0_rt*zeta*zbari; + if constexpr (do_derivatives) { + zetadt = -zeta*tempi; + zetada = 0.0e0_rt; + zetadz = 2.0e0_rt*zeta*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; - dumdt = zetadt*c00 + zeta*c01*nudt; - dumda = zeta*c01*nuda; - dumdz = zetadz*c00 + zeta*c01*nudz; + 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); @@ -1209,17 +1331,21 @@ void sneut5(const Real temp, const Real den, z = std::exp(c*nu); dd00 = b*z*(1.0e0_rt + d*dum); gum = 1.0e0_rt + dd00; - gumdt = dd00*c*nudt + b*z*d*dumdt; - gumda = dd00*c*nuda + b*z*d*dumda; - gumdz = dd00*c*nudz + b*z*d*dumdz; + 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; - 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; + 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); @@ -1228,45 +1354,58 @@ void sneut5(const Real temp, const Real den, a2 = 1.0e0_rt/bigj; sreco = tfac6 * 2.649e-18_rt * ye * std::pow(zbar, 13.0_rt) * den * bigj*a1; - 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); - + 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); + } } // convert from erg/cm^3/s to erg/g/s // comment these out to duplicate the itoh et al plots spair = spair*deni; - spairdt = spairdt*deni; - spairda = spairda*deni; - spairdz = spairdz*deni; + if constexpr (do_derivatives) { + spairdt = spairdt*deni; + spairda = spairda*deni; + spairdz = spairdz*deni; + } splas = splas*deni; - splasdt = splasdt*deni; - splasda = splasda*deni; - splasdz = splasdz*deni; + if constexpr (do_derivatives) { + splasdt = splasdt*deni; + splasda = splasda*deni; + splasdz = splasdz*deni; + } sphot = sphot*deni; - sphotdt = sphotdt*deni; - sphotda = sphotda*deni; - sphotdz = sphotdz*deni; + if constexpr (do_derivatives) { + sphotdt = sphotdt*deni; + sphotda = sphotda*deni; + sphotdz = sphotdz*deni; + } sbrem = sbrem*deni; - sbremdt = sbremdt*deni; - sbremda = sbremda*deni; - sbremdz = sbremdz*deni; + if constexpr (do_derivatives) { + sbremdt = sbremdt*deni; + sbremda = sbremda*deni; + sbremdz = sbremdz*deni; + } sreco = sreco*deni; - srecodt = srecodt*deni; - srecoda = srecoda*deni; - srecodz = srecodz*deni; + if constexpr (do_derivatives) { + srecodt = srecodt*deni; + srecoda = srecoda*deni; + srecodz = srecodz*deni; + } // the total neutrino loss rate snu = splas + spair + sphot + sbrem + sreco; - dsnudt = splasdt + spairdt + sphotdt + sbremdt + srecodt; - dsnuda = splasda + spairda + sphotda + sbremda + srecoda; - dsnudz = splasdz + spairdz + sphotdz + sbremdz + srecodz; + if constexpr (do_derivatives) { + dsnudt = splasdt + spairdt + sphotdt + sbremdt + srecodt; + dsnuda = splasda + spairda + sphotda + sbremda + srecoda; + dsnudz = splasdz + spairdz + sphotdz + sbremdz + srecodz; + } } diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H index 8ce878f196..1957dfe9e3 100644 --- a/nse_solver/nse_check.H +++ b/nse_solver/nse_check.H @@ -1070,7 +1070,8 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { composition(state); amrex::Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); #else amrex::Real sneut = 0.0_rt; #endif diff --git a/unit_test/test_neutrino_cooling/neutrino_util.cpp b/unit_test/test_neutrino_cooling/neutrino_util.cpp index 95a4799e96..efa6c70bb4 100644 --- a/unit_test/test_neutrino_cooling/neutrino_util.cpp +++ b/unit_test/test_neutrino_cooling/neutrino_util.cpp @@ -80,8 +80,10 @@ void neut_test_C(const Box& bx, Real dsnuda; Real dsnudz; - sneut5(temp_zone, dens_zone, abar, zbar, - snu, dsnudt, dsnudd, dsnuda, dsnudz); + constexpr int do_derivatives{1}; + + sneut5(temp_zone, dens_zone, abar, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); sp(i, j, k, vars.isneut) = snu; sp(i, j, k, vars.isneutdt) = dsnudt;