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/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 4a1b5aaf19..81449010fb 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -808,32 +808,35 @@ 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; @@ -1039,13 +1042,15 @@ void sneut5(const Real temp, const Real den, } 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;