From 5dd1a910e5d662c94108cc29964be51071cbbfdb Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 24 Oct 2023 11:09:24 -0400 Subject: [PATCH 01/10] update Urca network with the p -> n e-capture rate from Langanke (#1359) --- .../URCA-simple/actual_network.H | 20 +-- .../URCA-simple/actual_network_data.cpp | 5 +- .../ignition_reaclib/URCA-simple/actual_rhs.H | 32 +++- .../URCA-simple/n-p_betadecay.dat | 148 ++++++++++++++++++ .../URCA-simple/p-n_electroncapture.dat | 148 ++++++++++++++++++ .../URCA-simple/reaclib_rates.H | 37 ----- .../URCA-simple/table_rates.H | 12 +- .../URCA-simple/table_rates_data.cpp | 26 +++ networks/ignition_reaclib/URCA-simple/urca.py | 22 +-- 9 files changed, 388 insertions(+), 62 deletions(-) create mode 100644 networks/ignition_reaclib/URCA-simple/n-p_betadecay.dat create mode 100644 networks/ignition_reaclib/URCA-simple/p-n_electroncapture.dat diff --git a/networks/ignition_reaclib/URCA-simple/actual_network.H b/networks/ignition_reaclib/URCA-simple/actual_network.H index e3c8577cab..f83ed5f7ff 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_network.H +++ b/networks/ignition_reaclib/URCA-simple/actual_network.H @@ -28,19 +28,20 @@ namespace Rates k_C12_C12_to_n_Mg23 = 2, k_C12_C12_to_p_Na23 = 3, k_He4_C12_to_O16 = 4, - k_n_to_p_weak_wc12 = 5, - k_Na23_to_Ne23 = 6, - k_Ne23_to_Na23 = 7, - NumRates = k_Ne23_to_Na23 + k_Na23_to_Ne23 = 5, + k_Ne23_to_Na23 = 6, + k_n_to_p = 7, + k_p_to_n = 8, + NumRates = k_p_to_n }; // number of reaclib rates - const int NrateReaclib = 5; + const int NrateReaclib = 4; // number of tabular rates - const int NrateTabular = 2; + const int NrateTabular = 4; // rate names -- note: the rates are 1-based, not zero-based, so we pad // this vector with rate_names[0] = "" so the indices line up with the @@ -52,9 +53,10 @@ namespace Rates "C12_C12_to_n_Mg23", // 2, "C12_C12_to_p_Na23", // 3, "He4_C12_to_O16", // 4, - "n_to_p_weak_wc12", // 5, - "Na23_to_Ne23", // 6, - "Ne23_to_Na23" // 7, + "Na23_to_Ne23", // 5, + "Ne23_to_Na23", // 6, + "n_to_p", // 7, + "p_to_n" // 8, }; } diff --git a/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp b/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp index 4a7e7368ae..39fcc09e09 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp +++ b/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp @@ -15,9 +15,10 @@ namespace NSE_INDEX -1, 3, 3, -1, 0, 8, -1, -1, 3, 3, -1, 1, 7, -1, -1, 2, 3, -1, -1, 4, -1, - -1, -1, 0, -1, -1, 1, -1, -1, -1, 7, -1, -1, 6, -1, - -1, -1, 6, -1, -1, 7, 6 + -1, -1, 6, -1, -1, 7, 5, + -1, -1, 0, -1, -1, 1, -1, + -1, -1, 1, -1, -1, 0, -1 }; } #endif diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H index 97f78124f3..a5768da770 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_rhs.H +++ b/networks/ignition_reaclib/URCA-simple/actual_rhs.H @@ -151,6 +151,22 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { } rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne23) * (edot_nu + edot_gamma); + tabular_evaluate(j_n_p_meta, j_n_p_rhoy, j_n_p_temp, j_n_p_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_n_to_p) = rate; + if constexpr (std::is_same::value) { + rate_eval.dscreened_rates_dT(k_n_to_p) = drate_dt; + } + rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (edot_nu + edot_gamma); + + tabular_evaluate(j_p_n_meta, j_p_n_rhoy, j_p_n_temp, j_p_n_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_p_to_n) = rate; + if constexpr (std::is_same::value) { + rate_eval.dscreened_rates_dT(k_p_to_n) = drate_dt; + } + rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma); + } @@ -163,12 +179,14 @@ void rhs_nuc(const burn_t& state, using namespace Rates; ydot_nuc(N) = - -screened_rates(k_n_to_p_weak_wc12)*Y(N) + + -screened_rates(k_n_to_p)*Y(N) + + screened_rates(k_p_to_n)*Y(H1) + 0.5*screened_rates(k_C12_C12_to_n_Mg23)*std::pow(Y(C12), 2)*state.rho; ydot_nuc(H1) = 0.5*screened_rates(k_C12_C12_to_p_Na23)*std::pow(Y(C12), 2)*state.rho + - screened_rates(k_n_to_p_weak_wc12)*Y(N); + screened_rates(k_n_to_p)*Y(N) + + -screened_rates(k_p_to_n)*Y(H1); ydot_nuc(He4) = 0.5*screened_rates(k_C12_C12_to_He4_Ne20)*std::pow(Y(C12), 2)*state.rho + @@ -254,15 +272,21 @@ void jac_nuc(const burn_t& state, Real scratch; - scratch = -screened_rates(k_n_to_p_weak_wc12); + scratch = -screened_rates(k_n_to_p); jac.set(N, N, scratch); + scratch = screened_rates(k_p_to_n); + jac.set(N, H1, scratch); + scratch = 1.0*screened_rates(k_C12_C12_to_n_Mg23)*Y(C12)*state.rho; jac.set(N, C12, scratch); - scratch = screened_rates(k_n_to_p_weak_wc12); + scratch = screened_rates(k_n_to_p); jac.set(H1, N, scratch); + scratch = -screened_rates(k_p_to_n); + jac.set(H1, H1, scratch); + scratch = 1.0*screened_rates(k_C12_C12_to_p_Na23)*Y(C12)*state.rho; jac.set(H1, C12, scratch); diff --git a/networks/ignition_reaclib/URCA-simple/n-p_betadecay.dat b/networks/ignition_reaclib/URCA-simple/n-p_betadecay.dat new file mode 100644 index 0000000000..31ddea3b8d --- /dev/null +++ b/networks/ignition_reaclib/URCA-simple/n-p_betadecay.dat @@ -0,0 +1,148 @@ +!n -> p, beta-decay +!Q=-1.2933 MeV +! +!Log(rhoY) Log(temp) mu dQ VS Log(beta-decay-rate) Log(nu-energy-loss) Log(gamma-energy) +!Log(g/cm^3) Log(K) erg erg erg Log(1/s) Log(erg/s) Log(erg/s) +1.000000 7.000000 -4.806530e-09 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +1.000000 8.000000 -9.292624e-08 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +1.000000 8.301030 -2.146917e-07 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +1.000000 8.602060 -4.902661e-07 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +1.000000 8.845098 -8.058948e-07 0.00 0.00 -2.958959 -9.072290e+00 -100.00 +1.000000 9.000000 -8.187123e-07 0.00 0.00 -2.957841 -9.067290e+00 -100.00 +1.000000 9.176091 -8.187123e-07 0.00 0.00 -2.941107 -8.997290e+00 -100.00 +1.000000 9.301030 -8.187123e-07 0.00 0.00 -2.874417 -8.765290e+00 -100.00 +1.000000 9.477121 -8.187123e-07 0.00 0.00 -2.540055 -8.086290e+00 -100.00 +1.000000 9.698970 -8.187123e-07 0.00 0.00 -1.720786 -7.027290e+00 -100.00 +1.000000 10.000000 -8.187123e-07 0.00 0.00 -0.414110 -5.490290e+00 -100.00 +1.000000 10.477121 -8.187123e-07 0.00 0.00 1.802004 -2.852290e+00 -100.00 +1.000000 11.000000 -8.187123e-07 0.00 0.00 4.347000 1.987104e-01 -100.00 +2.000000 7.000000 -1.602177e-09 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +2.000000 8.000000 -6.088271e-08 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +2.000000 8.301030 -1.522068e-07 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +2.000000 8.602060 -3.636941e-07 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +2.000000 8.845098 -7.145708e-07 0.00 0.00 -2.958984 -9.072290e+00 -100.00 +2.000000 9.000000 -8.107014e-07 0.00 0.00 -2.957898 -9.067290e+00 -100.00 +2.000000 9.176091 -8.187123e-07 0.00 0.00 -2.941192 -8.997290e+00 -100.00 +2.000000 9.301030 -8.187123e-07 0.00 0.00 -2.874417 -8.765290e+00 -100.00 +2.000000 9.477121 -8.187123e-07 0.00 0.00 -2.540055 -8.086290e+00 -100.00 +2.000000 9.698970 -8.187123e-07 0.00 0.00 -1.720786 -7.027290e+00 -100.00 +2.000000 10.000000 -8.187123e-07 0.00 0.00 -0.414110 -5.490290e+00 -100.00 +2.000000 10.477121 -8.187123e-07 0.00 0.00 1.802004 -2.852290e+00 -100.00 +2.000000 11.000000 -8.187123e-07 0.00 0.00 4.347000 1.987104e-01 -100.00 +3.000000 7.000000 3.204353e-09 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +3.000000 8.000000 -2.883918e-08 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +3.000000 8.301030 -8.811971e-08 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +3.000000 8.602060 -2.355200e-07 0.00 0.00 -2.959000 -9.072290e+00 -100.00 +3.000000 8.845098 -5.030835e-07 0.00 0.00 -2.958998 -9.072290e+00 -100.00 +3.000000 9.000000 -7.450121e-07 0.00 0.00 -2.958314 -9.069290e+00 -100.00 +3.000000 9.176091 -8.107014e-07 0.00 0.00 -2.941776 -9.000290e+00 -100.00 +3.000000 9.301030 -8.171101e-07 0.00 0.00 -2.874975 -8.767290e+00 -100.00 +3.000000 9.477121 -8.187123e-07 0.00 0.00 -2.540055 -8.087290e+00 -100.00 +3.000000 9.698970 -8.187123e-07 0.00 0.00 -1.721736 -7.027290e+00 -100.00 +3.000000 10.000000 -8.187123e-07 0.00 0.00 -0.414110 -5.490290e+00 -100.00 +3.000000 10.477121 -8.187123e-07 0.00 0.00 1.802004 -2.852290e+00 -100.00 +3.000000 11.000000 -8.187123e-07 0.00 0.00 4.347000 1.987104e-01 -100.00 +4.000000 7.000000 1.922612e-08 0.00 0.00 -2.961000 -9.076290e+00 -100.00 +4.000000 8.000000 8.010883e-09 0.00 0.00 -2.961000 -9.075290e+00 -100.00 +4.000000 8.301030 -2.082830e-08 0.00 0.00 -2.961000 -9.075290e+00 -100.00 +4.000000 8.602060 -1.073458e-07 0.00 0.00 -2.961000 -9.075290e+00 -100.00 +4.000000 8.845098 -2.787787e-07 0.00 0.00 -2.960000 -9.074290e+00 -100.00 +4.000000 9.000000 -4.838573e-07 0.00 0.00 -2.959896 -9.073290e+00 -100.00 +4.000000 9.176091 -7.434100e-07 0.00 0.00 -2.947722 -9.020290e+00 -100.00 +4.000000 9.301030 -7.962818e-07 0.00 0.00 -2.881381 -8.783290e+00 -100.00 +4.000000 9.477121 -8.123036e-07 0.00 0.00 -2.544233 -8.091290e+00 -100.00 +4.000000 9.698970 -8.171101e-07 0.00 0.00 -1.721736 -7.028290e+00 -100.00 +4.000000 10.000000 -8.187123e-07 0.00 0.00 -0.414110 -5.491290e+00 -100.00 +4.000000 10.477121 -8.187123e-07 0.00 0.00 1.802004 -2.852290e+00 -100.00 +4.000000 11.000000 -8.187123e-07 0.00 0.00 4.347000 1.987104e-01 -100.00 +5.000000 7.000000 8.491536e-08 0.00 0.00 -2.980000 -9.105290e+00 -100.00 +5.000000 8.000000 8.331318e-08 0.00 0.00 -2.979000 -9.104290e+00 -100.00 +5.000000 8.301030 7.530230e-08 0.00 0.00 -2.979000 -9.103290e+00 -100.00 +5.000000 8.602060 4.165659e-08 0.00 0.00 -2.978000 -9.100290e+00 -100.00 +5.000000 8.845098 -4.165659e-08 0.00 0.00 -2.975000 -9.096290e+00 -100.00 +5.000000 9.000000 -1.570133e-07 0.00 0.00 -2.972990 -9.092290e+00 -100.00 +5.000000 9.176091 -3.941355e-07 0.00 0.00 -2.967462 -9.076290e+00 -100.00 +5.000000 9.301030 -6.136337e-07 0.00 0.00 -2.925218 -8.906290e+00 -100.00 +5.000000 9.477121 -7.626361e-07 0.00 0.00 -2.577269 -8.139290e+00 -100.00 +5.000000 9.698970 -8.042927e-07 0.00 0.00 -1.729383 -7.036290e+00 -100.00 +5.000000 10.000000 -8.155079e-07 0.00 0.00 -0.415108 -5.491290e+00 -100.00 +5.000000 10.477121 -8.187123e-07 0.00 0.00 1.802004 -2.852290e+00 -100.00 +5.000000 11.000000 -8.187123e-07 0.00 0.00 4.347000 1.987104e-01 -100.00 +6.000000 7.000000 3.444680e-07 0.00 0.00 -3.140000 -9.345290e+00 -100.00 +6.000000 8.000000 3.428658e-07 0.00 0.00 -3.140000 -9.343290e+00 -100.00 +6.000000 8.301030 3.412636e-07 0.00 0.00 -3.138000 -9.340290e+00 -100.00 +6.000000 8.602060 3.316506e-07 0.00 0.00 -3.131000 -9.325290e+00 -100.00 +6.000000 8.845098 3.028114e-07 0.00 0.00 -3.116000 -9.295290e+00 -100.00 +6.000000 9.000000 2.579504e-07 0.00 0.00 -3.097999 -9.262290e+00 -100.00 +6.000000 9.176091 1.490024e-07 0.00 0.00 -3.069766 -9.215290e+00 -100.00 +6.000000 9.301030 1.602177e-09 0.00 0.00 -3.041782 -9.153290e+00 -100.00 +6.000000 9.477121 -3.396614e-07 0.00 0.00 -2.813017 -8.512290e+00 -100.00 +6.000000 9.698970 -6.729142e-07 0.00 0.00 -1.806664 -7.116290e+00 -100.00 +6.000000 10.000000 -7.866687e-07 0.00 0.00 -0.423096 -5.500290e+00 -100.00 +6.000000 10.477121 -8.155079e-07 0.00 0.00 1.801004 -2.852290e+00 -100.00 +6.000000 11.000000 -8.187123e-07 0.00 0.00 4.347000 1.987104e-01 -100.00 +7.000000 7.000000 1.140750e-06 0.00 0.00 -5.499000 -1.257129e+01 -100.00 +7.000000 8.000000 1.139148e-06 0.00 0.00 -5.442000 -1.246129e+01 -100.00 +7.000000 8.301030 1.139148e-06 0.00 0.00 -5.306000 -1.221829e+01 -100.00 +7.000000 8.602060 1.134341e-06 0.00 0.00 -4.989000 -1.170729e+01 -100.00 +7.000000 8.845098 1.123126e-06 0.00 0.00 -4.592000 -1.112229e+01 -100.00 +7.000000 9.000000 1.103900e-06 0.00 0.00 -4.297000 -1.071129e+01 -100.00 +7.000000 9.176091 1.060641e-06 0.00 0.00 -3.954978 -1.025629e+01 -100.00 +7.000000 9.301030 9.965539e-07 0.00 0.00 -3.726186 -9.962290e+00 -100.00 +7.000000 9.477121 8.171101e-07 0.00 0.00 -3.408689 -9.420290e+00 -100.00 +7.000000 9.698970 2.996070e-07 0.00 0.00 -2.366458 -7.714290e+00 -100.00 +7.000000 10.000000 -5.046856e-07 0.00 0.00 -0.508967 -5.587290e+00 -100.00 +7.000000 10.477121 -7.850666e-07 0.00 0.00 1.798004 -2.855290e+00 -100.00 +7.000000 11.000000 -8.155079e-07 0.00 0.00 4.347000 1.987104e-01 -100.00 +8.000000 7.000000 3.101814e-06 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +8.000000 8.000000 3.101814e-06 0.00 0.00 -65.570000 -7.295929e+01 -100.00 +8.000000 8.301030 3.100212e-06 0.00 0.00 -35.616000 -4.271029e+01 -100.00 +8.000000 8.602060 3.098610e-06 0.00 0.00 -20.209000 -2.701829e+01 -100.00 +8.000000 8.845098 3.093803e-06 0.00 0.00 -13.299000 -1.988929e+01 -100.00 +8.000000 9.000000 3.085792e-06 0.00 0.00 -10.394000 -1.685929e+01 -100.00 +8.000000 9.176091 3.064964e-06 0.00 0.00 -8.019984 -1.436829e+01 -100.00 +8.000000 9.301030 3.036125e-06 0.00 0.00 -6.772438 -1.305429e+01 -100.00 +8.000000 9.477121 2.954414e-06 0.00 0.00 -5.436532 -1.153729e+01 -100.00 +8.000000 9.698970 2.693259e-06 0.00 0.00 -3.791028 -9.205290e+00 -100.00 +8.000000 10.000000 1.573337e-06 0.00 0.00 -1.150873 -6.234290e+00 -100.00 +8.000000 10.477121 -4.838573e-07 0.00 0.00 1.768005 -2.886290e+00 -100.00 +8.000000 11.000000 -7.882709e-07 0.00 0.00 4.346000 1.977104e-01 -100.00 +9.000000 7.000000 7.480563e-06 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +9.000000 8.000000 7.480563e-06 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +9.000000 8.301030 7.480563e-06 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +9.000000 8.602060 7.480563e-06 0.00 0.00 -54.663000 -6.147129e+01 -100.00 +9.000000 8.845098 7.477358e-06 0.00 0.00 -33.000000 -3.959029e+01 -100.00 +9.000000 9.000000 7.474154e-06 0.00 0.00 -24.198000 -3.066429e+01 -100.00 +9.000000 9.176091 7.464541e-06 0.00 0.00 -17.245984 -2.359429e+01 -100.00 +9.000000 9.301030 7.450121e-06 0.00 0.00 -13.715440 -1.999829e+01 -100.00 +9.000000 9.477121 7.413271e-06 0.00 0.00 -10.109646 -1.621129e+01 -100.00 +9.000000 9.698970 7.291506e-06 0.00 0.00 -6.677964 -1.209729e+01 -100.00 +9.000000 10.000000 6.724335e-06 0.00 0.00 -2.765936 -7.852290e+00 -100.00 +9.000000 10.477121 2.340780e-06 0.00 0.00 1.478006 -3.179290e+00 -100.00 +9.000000 11.000000 -5.191052e-07 0.00 0.00 4.338000 1.897104e-01 -100.00 +10.000000 7.000000 1.699429e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +10.000000 8.000000 1.699429e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +10.000000 8.301030 1.699429e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +10.000000 8.602060 1.699429e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +10.000000 8.845098 1.699269e-05 0.00 0.00 -75.759000 -8.234929e+01 -100.00 +10.000000 9.000000 1.699108e-05 0.00 0.00 -54.136000 -6.060229e+01 -100.00 +10.000000 9.176091 1.698628e-05 0.00 0.00 -37.214984 -4.356329e+01 -100.00 +10.000000 9.301030 1.697987e-05 0.00 0.00 -28.703440 -3.498629e+01 -100.00 +10.000000 9.477121 1.696224e-05 0.00 0.00 -20.122703 -2.622429e+01 -100.00 +10.000000 9.698970 1.690617e-05 0.00 0.00 -12.726964 -1.814629e+01 -100.00 +10.000000 10.000000 1.664181e-05 0.00 0.00 -5.885901 -1.097129e+01 -100.00 +10.000000 10.477121 1.386203e-05 0.00 0.00 0.277010 -4.383290e+00 -100.00 +10.000000 11.000000 2.164541e-06 0.00 0.00 4.255000 1.057104e-01 -100.00 +11.000000 7.000000 3.752778e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +11.000000 8.000000 3.752778e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +11.000000 8.301030 3.752778e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +11.000000 8.602060 3.752778e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +11.000000 8.845098 3.752618e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +11.000000 9.000000 3.752618e-05 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +11.000000 9.176091 3.752298e-05 0.00 0.00 -80.281984 -8.663029e+01 -100.00 +11.000000 9.301030 3.752137e-05 0.00 0.00 -61.008440 -6.729129e+01 -100.00 +11.000000 9.477121 3.751176e-05 0.00 0.00 -41.669646 -4.777129e+01 -100.00 +11.000000 9.698970 3.748613e-05 0.00 0.00 -25.673964 -3.109329e+01 -100.00 +11.000000 10.000000 3.736436e-05 0.00 0.00 -12.403901 -1.748929e+01 -100.00 +11.000000 10.477121 3.605538e-05 0.00 0.00 -2.049990 -6.710290e+00 -100.00 +11.000000 11.000000 2.244810e-05 0.00 0.00 3.625000 -5.272896e-01 -100.00 diff --git a/networks/ignition_reaclib/URCA-simple/p-n_electroncapture.dat b/networks/ignition_reaclib/URCA-simple/p-n_electroncapture.dat new file mode 100644 index 0000000000..21f333657d --- /dev/null +++ b/networks/ignition_reaclib/URCA-simple/p-n_electroncapture.dat @@ -0,0 +1,148 @@ +!p -> n, e- capture +!Q=1.2933 MeV +! +!Log(rhoY) Log(temp) mu dQ Vs Log(e-cap-rate) Log(nu-energy-loss) Log(gamma-energy) +!Log(g/cm^3) Log(K) erg erg erg Log(1/s) Log(erg/s) Log(erg/s) +1.000000 7.000000 -4.806530e-09 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +1.000000 8.000000 -9.292624e-08 0.00 0.00 -49.750000 -5.712729e+01 -100.00 +1.000000 8.301030 -2.146917e-07 0.00 0.00 -29.580000 -3.665029e+01 -100.00 +1.000000 8.602060 -4.902661e-07 0.00 0.00 -19.262000 -2.602029e+01 -100.00 +1.000000 8.845098 -8.058948e-07 0.00 0.00 -14.019000 -2.051929e+01 -100.00 +1.000000 9.000000 -8.187123e-07 0.00 0.00 -10.766000 -1.709829e+01 -100.00 +1.000000 9.176091 -8.187123e-07 0.00 0.00 -7.990000 -1.412829e+01 -100.00 +1.000000 9.301030 -8.187123e-07 0.00 0.00 -6.458000 -1.245629e+01 -100.00 +1.000000 9.477121 -8.187123e-07 0.00 0.00 -4.715000 -1.051429e+01 -100.00 +1.000000 9.698970 -8.187123e-07 0.00 0.00 -2.968000 -8.516290e+00 -100.00 +1.000000 10.000000 -8.187123e-07 0.00 0.00 -1.035000 -6.245290e+00 -100.00 +1.000000 10.477121 -8.187123e-07 0.00 0.00 1.600000 -3.099290e+00 -100.00 +1.000000 11.000000 -8.187123e-07 0.00 0.00 4.293000 1.317104e-01 -100.00 +2.000000 7.000000 -1.602177e-09 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +2.000000 8.000000 -6.088271e-08 0.00 0.00 -48.749000 -5.612529e+01 -100.00 +2.000000 8.301030 -1.522068e-07 0.00 0.00 -28.580000 -3.564929e+01 -100.00 +2.000000 8.602060 -3.636941e-07 0.00 0.00 -18.262000 -2.502029e+01 -100.00 +2.000000 8.845098 -7.145708e-07 0.00 0.00 -13.606000 -2.010629e+01 -100.00 +2.000000 9.000000 -8.107014e-07 0.00 0.00 -10.744000 -1.707729e+01 -100.00 +2.000000 9.176091 -8.187123e-07 0.00 0.00 -7.989000 -1.412729e+01 -100.00 +2.000000 9.301030 -8.187123e-07 0.00 0.00 -6.458000 -1.245629e+01 -100.00 +2.000000 9.477121 -8.187123e-07 0.00 0.00 -4.715000 -1.051429e+01 -100.00 +2.000000 9.698970 -8.187123e-07 0.00 0.00 -2.968000 -8.516290e+00 -100.00 +2.000000 10.000000 -8.187123e-07 0.00 0.00 -1.035000 -6.245290e+00 -100.00 +2.000000 10.477121 -8.187123e-07 0.00 0.00 1.600000 -3.099290e+00 -100.00 +2.000000 11.000000 -8.187123e-07 0.00 0.00 4.293000 1.317104e-01 -100.00 +3.000000 7.000000 3.204353e-09 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +3.000000 8.000000 -2.883918e-08 0.00 0.00 -47.732000 -5.510929e+01 -100.00 +3.000000 8.301030 -8.811971e-08 0.00 0.00 -27.574000 -3.464429e+01 -100.00 +3.000000 8.602060 -2.355200e-07 0.00 0.00 -17.260000 -2.401829e+01 -100.00 +3.000000 8.845098 -5.030835e-07 0.00 0.00 -12.658000 -1.915829e+01 -100.00 +3.000000 9.000000 -7.450121e-07 0.00 0.00 -10.538000 -1.687029e+01 -100.00 +3.000000 9.176091 -8.107014e-07 0.00 0.00 -7.974000 -1.411229e+01 -100.00 +3.000000 9.301030 -8.171101e-07 0.00 0.00 -6.454000 -1.245329e+01 -100.00 +3.000000 9.477121 -8.187123e-07 0.00 0.00 -4.714000 -1.051429e+01 -100.00 +3.000000 9.698970 -8.187123e-07 0.00 0.00 -2.968000 -8.515290e+00 -100.00 +3.000000 10.000000 -8.187123e-07 0.00 0.00 -1.035000 -6.245290e+00 -100.00 +3.000000 10.477121 -8.187123e-07 0.00 0.00 1.600000 -3.099290e+00 -100.00 +3.000000 11.000000 -8.187123e-07 0.00 0.00 4.293000 1.317104e-01 -100.00 +4.000000 7.000000 1.922612e-08 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +4.000000 8.000000 8.010883e-09 0.00 0.00 -46.572000 -5.394829e+01 -100.00 +4.000000 8.301030 -2.082830e-08 0.00 0.00 -26.519000 -3.358929e+01 -100.00 +4.000000 8.602060 -1.073458e-07 0.00 0.00 -16.243000 -2.300029e+01 -100.00 +4.000000 8.845098 -2.787787e-07 0.00 0.00 -11.652000 -1.815229e+01 -100.00 +4.000000 9.000000 -4.838573e-07 0.00 0.00 -9.716000 -1.604929e+01 -100.00 +4.000000 9.176091 -7.434100e-07 0.00 0.00 -7.833000 -1.397129e+01 -100.00 +4.000000 9.301030 -7.962818e-07 0.00 0.00 -6.423000 -1.242129e+01 -100.00 +4.000000 9.477121 -8.123036e-07 0.00 0.00 -4.709000 -1.050829e+01 -100.00 +4.000000 9.698970 -8.171101e-07 0.00 0.00 -2.967000 -8.515290e+00 -100.00 +4.000000 10.000000 -8.187123e-07 0.00 0.00 -1.035000 -6.245290e+00 -100.00 +4.000000 10.477121 -8.187123e-07 0.00 0.00 1.600000 -3.099290e+00 -100.00 +4.000000 11.000000 -8.187123e-07 0.00 0.00 4.293000 1.317104e-01 -100.00 +5.000000 7.000000 8.491536e-08 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +5.000000 8.000000 8.331318e-08 0.00 0.00 -44.223000 -5.160029e+01 -100.00 +5.000000 8.301030 7.530230e-08 0.00 0.00 -25.017000 -3.208629e+01 -100.00 +5.000000 8.602060 4.165659e-08 0.00 0.00 -15.072000 -2.182929e+01 -100.00 +5.000000 8.845098 -4.165659e-08 0.00 0.00 -10.585000 -1.708529e+01 -100.00 +5.000000 9.000000 -1.570133e-07 0.00 0.00 -8.685000 -1.501729e+01 -100.00 +5.000000 9.176091 -3.941355e-07 0.00 0.00 -7.099000 -1.323729e+01 -100.00 +5.000000 9.301030 -6.136337e-07 0.00 0.00 -6.134000 -1.213329e+01 -100.00 +5.000000 9.477121 -7.626361e-07 0.00 0.00 -4.656000 -1.045529e+01 -100.00 +5.000000 9.698970 -8.042927e-07 0.00 0.00 -2.959000 -8.506290e+00 -100.00 +5.000000 10.000000 -8.155079e-07 0.00 0.00 -1.034000 -6.244290e+00 -100.00 +5.000000 10.477121 -8.187123e-07 0.00 0.00 1.600000 -3.099290e+00 -100.00 +5.000000 11.000000 -8.187123e-07 0.00 0.00 4.293000 1.317104e-01 -100.00 +6.000000 7.000000 3.444680e-07 0.00 0.00 -99.697970 -1.057943e+02 -100.00 +6.000000 8.000000 3.428658e-07 0.00 0.00 -36.029000 -4.340629e+01 -100.00 +6.000000 8.301030 3.412636e-07 0.00 0.00 -20.834000 -2.790429e+01 -100.00 +6.000000 8.602060 3.316506e-07 0.00 0.00 -12.798000 -1.955529e+01 -100.00 +6.000000 8.845098 3.028114e-07 0.00 0.00 -9.036000 -1.553629e+01 -100.00 +6.000000 9.000000 2.579504e-07 0.00 0.00 -7.382000 -1.371429e+01 -100.00 +6.000000 9.176091 1.490024e-07 0.00 0.00 -5.962000 -1.210029e+01 -100.00 +6.000000 9.301030 1.602177e-09 0.00 0.00 -5.168000 -1.116629e+01 -100.00 +6.000000 9.477121 -3.396614e-07 0.00 0.00 -4.214000 -1.001329e+01 -100.00 +6.000000 9.698970 -6.729142e-07 0.00 0.00 -2.876000 -8.424290e+00 -100.00 +6.000000 10.000000 -7.866687e-07 0.00 0.00 -1.025000 -6.235290e+00 -100.00 +6.000000 10.477121 -8.155079e-07 0.00 0.00 1.601000 -3.098290e+00 -100.00 +6.000000 11.000000 -8.187123e-07 0.00 0.00 4.293000 1.317104e-01 -100.00 +7.000000 7.000000 1.140750e-06 0.00 0.00 -46.082000 -5.446429e+01 -100.00 +7.000000 8.000000 1.139148e-06 0.00 0.00 -10.983000 -1.835929e+01 -100.00 +7.000000 8.301030 1.139148e-06 0.00 0.00 -8.290000 -1.535929e+01 -100.00 +7.000000 8.602060 1.134341e-06 0.00 0.00 -6.487000 -1.324229e+01 -100.00 +7.000000 8.845098 1.123126e-06 0.00 0.00 -5.363000 -1.185729e+01 -100.00 +7.000000 9.000000 1.103900e-06 0.00 0.00 -4.733000 -1.105829e+01 -100.00 +7.000000 9.176091 1.060641e-06 0.00 0.00 -4.066000 -1.019629e+01 -100.00 +7.000000 9.301030 9.965539e-07 0.00 0.00 -3.616000 -9.607290e+00 -100.00 +7.000000 9.477121 8.171101e-07 0.00 0.00 -3.010000 -8.803290e+00 -100.00 +7.000000 9.698970 2.996070e-07 0.00 0.00 -2.270000 -7.815290e+00 -100.00 +7.000000 10.000000 -5.046856e-07 0.00 0.00 -0.937000 -6.147290e+00 -100.00 +7.000000 10.477121 -7.850666e-07 0.00 0.00 1.604000 -3.095290e+00 -100.00 +7.000000 11.000000 -8.155079e-07 0.00 0.00 4.293000 1.317104e-01 -100.00 +8.000000 7.000000 3.101814e-06 0.00 0.00 -1.350000 -7.185290e+00 -100.00 +8.000000 8.000000 3.101814e-06 0.00 0.00 -1.349000 -7.185290e+00 -100.00 +8.000000 8.301030 3.100212e-06 0.00 0.00 -1.348000 -7.183290e+00 -100.00 +8.000000 8.602060 3.098610e-06 0.00 0.00 -1.345000 -7.176290e+00 -100.00 +8.000000 8.845098 3.093803e-06 0.00 0.00 -1.335000 -7.158290e+00 -100.00 +8.000000 9.000000 3.085792e-06 0.00 0.00 -1.320000 -7.129290e+00 -100.00 +8.000000 9.176091 3.064964e-06 0.00 0.00 -1.285000 -7.065290e+00 -100.00 +8.000000 9.301030 3.036125e-06 0.00 0.00 -1.240000 -6.983290e+00 -100.00 +8.000000 9.477121 2.954414e-06 0.00 0.00 -1.128000 -6.787290e+00 -100.00 +8.000000 9.698970 2.693259e-06 0.00 0.00 -0.875000 -6.369290e+00 -100.00 +8.000000 10.000000 1.573337e-06 0.00 0.00 -0.301000 -5.502290e+00 -100.00 +8.000000 10.477121 -4.838573e-07 0.00 0.00 1.635000 -3.064290e+00 -100.00 +8.000000 11.000000 -7.882709e-07 0.00 0.00 4.294000 1.327104e-01 -100.00 +9.000000 7.000000 7.480563e-06 0.00 0.00 0.831000 -4.464290e+00 -100.00 +9.000000 8.000000 7.480563e-06 0.00 0.00 0.831000 -4.464290e+00 -100.00 +9.000000 8.301030 7.480563e-06 0.00 0.00 0.832000 -4.464290e+00 -100.00 +9.000000 8.602060 7.480563e-06 0.00 0.00 0.832000 -4.463290e+00 -100.00 +9.000000 8.845098 7.477358e-06 0.00 0.00 0.833000 -4.461290e+00 -100.00 +9.000000 9.000000 7.474154e-06 0.00 0.00 0.835000 -4.458290e+00 -100.00 +9.000000 9.176091 7.464541e-06 0.00 0.00 0.838000 -4.451290e+00 -100.00 +9.000000 9.301030 7.450121e-06 0.00 0.00 0.844000 -4.442290e+00 -100.00 +9.000000 9.477121 7.413271e-06 0.00 0.00 0.859000 -4.414290e+00 -100.00 +9.000000 9.698970 7.291506e-06 0.00 0.00 0.904000 -4.333290e+00 -100.00 +9.000000 10.000000 6.724335e-06 0.00 0.00 1.074000 -4.035290e+00 -100.00 +9.000000 10.477121 2.340780e-06 0.00 0.00 1.922000 -2.772290e+00 -100.00 +9.000000 11.000000 -5.191052e-07 0.00 0.00 4.302000 1.407104e-01 -100.00 +10.000000 7.000000 1.699429e-05 0.00 0.00 2.676000 -2.211290e+00 -100.00 +10.000000 8.000000 1.699429e-05 0.00 0.00 2.676000 -2.211290e+00 -100.00 +10.000000 8.301030 1.699429e-05 0.00 0.00 2.676000 -2.211290e+00 -100.00 +10.000000 8.602060 1.699429e-05 0.00 0.00 2.676000 -2.211290e+00 -100.00 +10.000000 8.845098 1.699269e-05 0.00 0.00 2.676000 -2.211290e+00 -100.00 +10.000000 9.000000 1.699108e-05 0.00 0.00 2.676000 -2.210290e+00 -100.00 +10.000000 9.176091 1.698628e-05 0.00 0.00 2.677000 -2.209290e+00 -100.00 +10.000000 9.301030 1.697987e-05 0.00 0.00 2.678000 -2.208290e+00 -100.00 +10.000000 9.477121 1.696224e-05 0.00 0.00 2.680000 -2.203290e+00 -100.00 +10.000000 9.698970 1.690617e-05 0.00 0.00 2.688000 -2.188290e+00 -100.00 +10.000000 10.000000 1.664181e-05 0.00 0.00 2.725000 -2.122290e+00 -100.00 +10.000000 10.477121 1.386203e-05 0.00 0.00 3.006000 -1.638290e+00 -100.00 +10.000000 11.000000 2.164541e-06 0.00 0.00 4.385000 2.237104e-01 -100.00 +11.000000 7.000000 3.752778e-05 0.00 0.00 4.416000 -1.062896e-01 -100.00 +11.000000 8.000000 3.752778e-05 0.00 0.00 4.416000 -1.062896e-01 -100.00 +11.000000 8.301030 3.752778e-05 0.00 0.00 4.416000 -1.062896e-01 -100.00 +11.000000 8.602060 3.752778e-05 0.00 0.00 4.416000 -1.062896e-01 -100.00 +11.000000 8.845098 3.752618e-05 0.00 0.00 4.416000 -1.062896e-01 -100.00 +11.000000 9.000000 3.752618e-05 0.00 0.00 4.416000 -1.052896e-01 -100.00 +11.000000 9.176091 3.752298e-05 0.00 0.00 4.416000 -1.052896e-01 -100.00 +11.000000 9.301030 3.752137e-05 0.00 0.00 4.416000 -1.052896e-01 -100.00 +11.000000 9.477121 3.751176e-05 0.00 0.00 4.417000 -1.042896e-01 -100.00 +11.000000 9.698970 3.748613e-05 0.00 0.00 4.419000 -1.012896e-01 -100.00 +11.000000 10.000000 3.736436e-05 0.00 0.00 4.426000 -8.828961e-02 -100.00 +11.000000 10.477121 3.605538e-05 0.00 0.00 4.499000 4.271039e-02 -100.00 +11.000000 11.000000 2.244810e-05 0.00 0.00 4.989000 8.437104e-01 -100.00 diff --git a/networks/ignition_reaclib/URCA-simple/reaclib_rates.H b/networks/ignition_reaclib/URCA-simple/reaclib_rates.H index 1e72dfcf2f..3cc069c026 100644 --- a/networks/ignition_reaclib/URCA-simple/reaclib_rates.H +++ b/networks/ignition_reaclib/URCA-simple/reaclib_rates.H @@ -168,37 +168,6 @@ void rate_He4_C12_to_O16(const tf_t& tfactors, Real& rate, Real& drate_dT) { } -template -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void rate_n_to_p_weak_wc12(const tf_t& tfactors, Real& rate, Real& drate_dT) { - - // n --> p - - rate = 0.0; - drate_dT = 0.0; - - Real ln_set_rate{0.0}; - Real dln_set_rate_dT9{0.0}; - Real set_rate{0.0}; - - // wc12w - ln_set_rate = -6.78161; - amrex::ignore_unused(tfactors); - - if constexpr (do_T_derivatives) { - dln_set_rate_dT9 = 0.0; - } - - // avoid underflows by zeroing rates in [0.0, 1.e-100] - ln_set_rate = std::max(ln_set_rate, -230.0); - set_rate = std::exp(ln_set_rate); - rate += set_rate; - if constexpr (do_T_derivatives) { - drate_dT += set_rate * dln_set_rate_dT9 / 1.0e9; - } - -} - template @@ -234,12 +203,6 @@ fill_reaclib_rates(const tf_t& tfactors, T& rate_eval) rate_eval.dscreened_rates_dT(k_He4_C12_to_O16) = drate_dT; } - rate_n_to_p_weak_wc12(tfactors, rate, drate_dT); - rate_eval.screened_rates(k_n_to_p_weak_wc12) = rate; - if constexpr (std::is_same::value) { - rate_eval.dscreened_rates_dT(k_n_to_p_weak_wc12) = drate_dT; - - } } diff --git a/networks/ignition_reaclib/URCA-simple/table_rates.H b/networks/ignition_reaclib/URCA-simple/table_rates.H index 20da307673..43ca83dc25 100644 --- a/networks/ignition_reaclib/URCA-simple/table_rates.H +++ b/networks/ignition_reaclib/URCA-simple/table_rates.H @@ -20,7 +20,7 @@ void init_tabular(); // Log(g/cm^3) Log(K) erg erg erg Log(1/s) Log(erg/s) Log(erg/s) // -const int num_tables = 2; +const int num_tables = 4; enum TableVars { @@ -65,6 +65,16 @@ namespace rate_tables extern AMREX_GPU_MANAGED Array1D j_Ne23_Na23_rhoy; extern AMREX_GPU_MANAGED Array1D j_Ne23_Na23_temp; + extern AMREX_GPU_MANAGED table_t j_n_p_meta; + extern AMREX_GPU_MANAGED Array3D j_n_p_data; + extern AMREX_GPU_MANAGED Array1D j_n_p_rhoy; + extern AMREX_GPU_MANAGED Array1D j_n_p_temp; + + extern AMREX_GPU_MANAGED table_t j_p_n_meta; + extern AMREX_GPU_MANAGED Array3D j_p_n_data; + extern AMREX_GPU_MANAGED Array1D j_p_n_rhoy; + extern AMREX_GPU_MANAGED Array1D j_p_n_temp; + } template diff --git a/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp b/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp index 2e7b8130be..12ed4a635a 100644 --- a/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp +++ b/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp @@ -18,6 +18,16 @@ namespace rate_tables AMREX_GPU_MANAGED Array1D j_Ne23_Na23_rhoy; AMREX_GPU_MANAGED Array1D j_Ne23_Na23_temp; + AMREX_GPU_MANAGED table_t j_n_p_meta; + AMREX_GPU_MANAGED Array3D j_n_p_data; + AMREX_GPU_MANAGED Array1D j_n_p_rhoy; + AMREX_GPU_MANAGED Array1D j_n_p_temp; + + AMREX_GPU_MANAGED table_t j_p_n_meta; + AMREX_GPU_MANAGED Array3D j_p_n_data; + AMREX_GPU_MANAGED Array1D j_p_n_rhoy; + AMREX_GPU_MANAGED Array1D j_p_n_temp; + } @@ -45,5 +55,21 @@ void init_tabular() init_tab_info(j_Ne23_Na23_meta, "23ne-23na_betadecay.dat", j_Ne23_Na23_rhoy, j_Ne23_Na23_temp, j_Ne23_Na23_data); + j_n_p_meta.ntemp = 13; + j_n_p_meta.nrhoy = 11; + j_n_p_meta.nvars = 6; + j_n_p_meta.nheader = 5; + + init_tab_info(j_n_p_meta, "n-p_betadecay.dat", j_n_p_rhoy, j_n_p_temp, j_n_p_data); + + + j_p_n_meta.ntemp = 13; + j_p_n_meta.nrhoy = 11; + j_p_n_meta.nvars = 6; + j_p_n_meta.nheader = 5; + + init_tab_info(j_p_n_meta, "p-n_electroncapture.dat", j_p_n_rhoy, j_p_n_temp, j_p_n_data); + + } diff --git a/networks/ignition_reaclib/URCA-simple/urca.py b/networks/ignition_reaclib/URCA-simple/urca.py index 5db0666b34..4c8e1c8288 100644 --- a/networks/ignition_reaclib/URCA-simple/urca.py +++ b/networks/ignition_reaclib/URCA-simple/urca.py @@ -1,14 +1,18 @@ # C-burning with A=23 URCA rate module generator -from pynucastro.networks import AmrexAstroCxxNetwork +import pynucastro as pyna -files = ["c12-c12a-ne20-cf88", - "c12-c12n-mg23-cf88", - "c12-c12p-na23-cf88", - "c12-ag-o16-nac2", - "na23--ne23-toki", - "ne23--na23-toki", - "n--p-wc12"] +rl = pyna.ReacLibLibrary() +rl_rates = rl.get_rate_by_name(["c12(c12,a)ne20", + "c12(c12,n)mg23", + "c12(c12,p)na23", + "c12(a,g)o16"]) -urca_net = AmrexAstroCxxNetwork(files) +tl = pyna.TabularLibrary() +tl_rates = tl.get_rate_by_name(["na23(,)ne23", + "ne23(,)na23", + "n(,)p", + "p(,)n"]) + +urca_net = pyna.AmrexAstroCxxNetwork(rates=rl_rates+tl_rates) urca_net.write_network() From 8b48c0281d9f18eac6a6ed167871bbb2bd433d37 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 13:09:04 -0400 Subject: [PATCH 02/10] 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 03/10] 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 04/10] 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; From fc40d670ae0e09e303eb4b4558a634a4926fe742 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 18:45:17 -0400 Subject: [PATCH 05/10] move sneut5 recombination section into its own function (#1373) --- neutrinos/sneut5.H | 404 ++++++++++++++++++++++++--------------------- 1 file changed, 219 insertions(+), 185 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 7f676b8525..e6a8f7aedc 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -6,6 +6,39 @@ using namespace amrex; +namespace nu_constants { + + // numerical constants + constexpr Real fac1 = 5.0e0_rt * M_PI / 3.0e0_rt; + constexpr Real fac2 = 10.0e0_rt * M_PI; + 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 sixth = 1.0e0_rt/6.0e0_rt; + constexpr Real iln10 = 4.342944819032518e-1_rt; + + // theta is sin**2(theta_weinberg) = 0.2319 plus/minus 0.00005 (1996) + // xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998) + // change theta and xnufam if need be, and the changes will automatically + // propagate through the routine. cv and ca are the vector and axial currents. + + constexpr Real theta = 0.2319e0_rt; + constexpr Real xnufam = 3.0e0_rt; + constexpr Real cv = 0.5e0_rt + 2.0e0_rt * theta; + constexpr Real cvp = 1.0e0_rt - cv; + constexpr Real ca = 0.5e0_rt; + constexpr Real cap = 1.0e0_rt - ca; + constexpr Real tfac1 = cv*cv + ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp+cap*cap); + constexpr Real tfac2 = cv*cv - ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp - cap*cap); + constexpr Real tfac3 = tfac2/tfac1; + constexpr Real tfac4 = 0.5e0_rt * tfac1; + constexpr Real tfac5 = 0.5e0_rt * tfac2; + constexpr Real tfac6 = cv*cv + 1.5e0_rt*ca*ca + (xnufam - 1.0e0_rt)*(cvp*cvp + 1.5e0_rt*cap*cap); + + +} + + AMREX_GPU_HOST_DEVICE AMREX_INLINE Real ifermi12(const Real f) { @@ -225,6 +258,11 @@ Real zfermim12(const Real x) struct sneutf_t { + Real den; + Real temp; + Real abar; + Real zbar; + Real deni; Real tempi; Real abari; @@ -264,6 +302,11 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { sneutf_t sf; + sf.den = den; + sf.temp = temp; + sf.abar = abar; + sf.zbar = zbar; + // to avoid lots of divisions sf.deni = 1.0e0_rt / den; sf.tempi = 1.0e0_rt / temp; @@ -302,6 +345,130 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_recomb(const sneutf_t& sf, + Real& sreco, Real& srecodt, Real& srecoda, Real& srecodz) { + + + // recombination neutrino section + // for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e + // equation 6.11 solved for nu + + Real xnum = 1.10520e8_rt * sf.den * sf.ye / (sf.temp * std::sqrt(sf.temp)); + Real xnumdt{0.0}; + Real xnumda{0.0}; + Real xnumdz{0.0}; + + if constexpr (do_derivatives) { + xnumdt = -1.50e0_rt*xnum*sf.tempi; + xnumda = -xnum*sf.abari; + xnumdz = xnum*sf.zbari; + } + + // the chemical potential + Real nu = ifermi12(xnum); + + // a0 is d(nu)/d(xnum) + Real a0 = 1.0e0_rt/(0.5e0_rt*zfermim12(nu)); + Real nudt, nuda, nudz; + if constexpr (do_derivatives) { + nudt = a0*xnumdt; + nuda = a0*xnumda; + nudz = a0*xnumdz; + } + Real nu2 = nu * nu; + Real nu3 = nu2 * nu; + + Real a1, a2, a3, b, c, d, f1, f2, f3; + + // table 12 + if (nu >= -20.0_rt && nu < 0.0_rt) { + a1 = 1.51e-2_rt; + a2 = 2.42e-1_rt; + a3 = 1.21e0_rt; + b = 3.71e-2_rt; + c = 9.06e-1_rt; + d = 9.28e-1_rt; + f1 = 0.0e0_rt; + f2 = 0.0e0_rt; + f3 = 0.0e0_rt; + } else if (nu >= 0.0_rt && nu <= 10.0_rt) { + a1 = 1.23e-2_rt; + a2 = 2.66e-1_rt; + a3 = 1.30e0_rt; + b = 1.17e-1_rt; + c = 8.97e-1_rt; + d = 1.77e-1_rt; + f1 = -1.20e-2_rt; + f2 = 2.29e-2_rt; + f3 = -1.04e-3_rt; + } + + // equation 6.7, 6.13 and 6.14 + if (nu >= -20.0_rt && nu <= 10.0_rt) { + + Real zeta = 1.579e5_rt * sf.zbar * sf.zbar * sf.tempi; + Real zetadt, zetada, zetadz; + if constexpr (do_derivatives) { + zetadt = -zeta * sf.tempi; + zetada = 0.0e0_rt; + zetadz = 2.0e0_rt * zeta * sf.zbari; + } + + Real c00 = 1.0e0_rt / (1.0e0_rt + f1 * nu + f2 * nu2 + f3 * nu3); + Real c01 = f1 + f2 * 2.0e0_rt * nu + f3 * 3.0e0_rt * nu2; + Real dum = zeta * c00; + Real dumdt, dumda, dumdz; + if constexpr (do_derivatives) { + dumdt = zetadt * c00 + zeta * c01 * nudt; + dumda = zeta * c01 * nuda; + dumdz = zetadz * c00 + zeta * c01 * nudz; + } + + Real z = 1.0e0_rt / dum; + Real dd00 = std::pow(dum, -2.25_rt); + Real dd01 = std::pow(dum, -4.55_rt); + c00 = a1 * z + a2 * dd00 + a3 * dd01; + c01 = -(a1 * z + 2.25_rt * a2 * dd00 + 4.55_rt * a3 * dd01) * z; + + z = std::exp(c * nu); + dd00 = b * z * (1.0e0_rt + d * dum); + Real gum = 1.0e0_rt + dd00; + Real gumdt, gumda, gumdz; + 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; + + Real bigj = c00 * z * a1; + Real bigjdt, bigjda, bigjdz; + 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); + dum = 1.0e0_rt + z; + a1 = 1.0e0_rt/dum; + a2 = 1.0e0_rt/bigj; + + sreco = nu_constants::tfac6 * 2.649e-18_rt * sf.ye * std::pow(sf.zbar, 13.0_rt) * sf.den * bigj * a1; + if constexpr (do_derivatives) { + srecodt = sreco * (bigjdt * a2 - z * (zetadt + nudt) * 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); + } + } + +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void sneut5(const Real temp, const Real den, @@ -331,9 +498,9 @@ 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{0.0},f0, - f1{0.0},f2,f3,z, - dum,dumdt,gum,gumdt,dumda,dumdz,gumda,gumdz; + dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,c,d{0.0},f0, + f1{0.0},f2,z, + dum,dumdt,gum,dumda,dumdz; // pair production Real gl,gldt, @@ -368,37 +535,6 @@ void sneut5(const Real temp, const Real den, fbremda,fbremdz,gbremda,gbremdz, fliqda,fliqdz,gliqda,gliqdz; - // recomb - Real nu,nudt,nu2,nu3,bigj,bigjdt,nuda,nudz,bigjda,bigjdz; - - - // numerical constants - constexpr Real fac1 = 5.0e0_rt * M_PI / 3.0e0_rt; - constexpr Real fac2 = 10.0e0_rt * M_PI; - 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 sixth = 1.0e0_rt/6.0e0_rt; - constexpr Real iln10 = 4.342944819032518e-1_rt; - - // theta is sin**2(theta_weinberg) = 0.2319 plus/minus 0.00005 (1996) - // xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998) - // change theta and xnufam if need be, and the changes will automatically - // propagate through the routine. cv and ca are the vector and axial currents. - - constexpr Real theta = 0.2319e0_rt; - constexpr Real xnufam = 3.0e0_rt; - constexpr Real cv = 0.5e0_rt + 2.0e0_rt * theta; - constexpr Real cvp = 1.0e0_rt - cv; - constexpr Real ca = 0.5e0_rt; - constexpr Real cap = 1.0e0_rt - ca; - constexpr Real tfac1 = cv*cv + ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp+cap*cap); - constexpr Real tfac2 = cv*cv - ca*ca + (xnufam-1.0e0_rt) * (cvp*cvp - cap*cap); - constexpr Real tfac3 = tfac2/tfac1; - constexpr Real tfac4 = 0.5e0_rt * tfac1; - constexpr Real tfac5 = 0.5e0_rt * tfac2; - constexpr Real tfac6 = cv*cv + 1.5e0_rt*ca*ca + (xnufam - 1.0e0_rt)*(cvp*cvp + 1.5e0_rt*cap*cap); - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; @@ -436,11 +572,11 @@ void sneut5(const Real temp, const Real den, auto sf = get_sneut_factors(den, temp, abar, zbar); a0 = sf.rm * 1.0e-9_rt; - a1 = std::pow(a0, oneth); + a1 = std::pow(a0, nu_constants::oneth); zeta = a1 * sf.xlm1; if constexpr (do_derivatives) { - a2 = oneth * a1*sf.rmi * sf.xlm1; + a2 = nu_constants::oneth * a1*sf.rmi * sf.xlm1; zetadt = -a1 * sf.xlm2 * sf.xldt; zetada = a2 * sf.rmda; zetadz = a2 * sf.rmdz; @@ -558,12 +694,12 @@ void sneut5(const Real temp, const Real den, spairdz = gl*spairdz; } - a1 = tfac4*(1.0e0_rt + tfac3 * qpair); + a1 = nu_constants::tfac4*(1.0e0_rt + nu_constants::tfac3 * qpair); a3 = spair; spair = a1*a3; if constexpr (do_derivatives) { - a2 = tfac4*tfac3; + a2 = nu_constants::tfac4 * nu_constants::tfac3; spairdt = a1*spairdt + a2*qpairdt*a3; spairda = a1*spairda + a2*qpairda*a3; spairdz = a1*spairdz + a2*qpairdz*a3; @@ -574,7 +710,7 @@ void sneut5(const Real temp, const Real den, // equation 4.6 a1 = 1.019e-6_rt*sf.rm; - a2 = std::pow(a1, twoth); + a2 = std::pow(a1, nu_constants::twoth); b1 = std::sqrt(1.0e0_rt + a2); @@ -583,7 +719,7 @@ void sneut5(const Real temp, const Real den, gl2 = 1.1095e11_rt * sf.rm * c00; if constexpr (do_derivatives) { - a3 = twoth*a2/a1; + a3 = nu_constants::twoth*a2/a1; b2 = 1.0e0_rt/b1; gl2dt = -2.0e0_rt*gl2*sf.tempi; d = sf.rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt; @@ -628,15 +764,15 @@ void sneut5(const Real temp, const Real den, 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); + xnum = nu_constants::sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt); + xden = nu_constants::sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt); if constexpr (do_derivatives) { - xnumdt = -iln10*0.5e0_rt*sf.tempi; - a2 = iln10*sixth*sf.rmi; + xnumdt = -nu_constants::iln10*0.5e0_rt*sf.tempi; + a2 = nu_constants::iln10*nu_constants::sixth*sf.rmi; xnumda = a2*sf.rmda; xnumdz = a2*sf.rmdz; - xdendt = iln10*0.5e0_rt*sf.tempi; + xdendt = nu_constants::iln10*0.5e0_rt*sf.tempi; xdenda = a2*sf.rmda; xdendz = a2*sf.rmdz; } @@ -863,22 +999,22 @@ void sneut5(const Real temp, const Real den, } - taudt = iln10*sf.tempi; + taudt = nu_constants::iln10*sf.tempi; // equation 3.7, compute the expensive trig functions only one time - cos1 = std::cos(fac1*tau); - cos2 = std::cos(fac1*2.0e0_rt*tau); - cos3 = std::cos(fac1*3.0e0_rt*tau); - cos4 = std::cos(fac1*4.0e0_rt*tau); - cos5 = std::cos(fac1*5.0e0_rt*tau); - last = std::cos(fac2*tau); - - sin1 = std::sin(fac1*tau); - sin2 = std::sin(fac1*2.0e0_rt*tau); - sin3 = std::sin(fac1*3.0e0_rt*tau); - sin4 = std::sin(fac1*4.0e0_rt*tau); - sin5 = std::sin(fac1*5.0e0_rt*tau); - xast = std::sin(fac2*tau); + cos1 = std::cos(nu_constants::fac1*tau); + cos2 = std::cos(nu_constants::fac1*2.0e0_rt*tau); + cos3 = std::cos(nu_constants::fac1*3.0e0_rt*tau); + cos4 = std::cos(nu_constants::fac1*4.0e0_rt*tau); + cos5 = std::cos(nu_constants::fac1*5.0e0_rt*tau); + last = std::cos(nu_constants::fac2*tau); + + sin1 = std::sin(nu_constants::fac1*tau); + sin2 = std::sin(nu_constants::fac1*2.0e0_rt*tau); + sin3 = std::sin(nu_constants::fac1*3.0e0_rt*tau); + sin4 = std::sin(nu_constants::fac1*4.0e0_rt*tau); + sin5 = std::sin(nu_constants::fac1*5.0e0_rt*tau); + xast = std::sin(nu_constants::fac2*tau); a0 = 0.5e0_rt*c00 + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2 @@ -896,21 +1032,21 @@ void sneut5(const Real temp, const Real den, + c25*cos5 + dd25*sin5 + 0.5e0_rt*c26*last; if constexpr (do_derivatives) { - f0 = taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0e0_rt + f0 = taudt*nu_constants::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; + - 0.5e0_rt*c06*xast*nu_constants::fac2*taudt; - f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0e0_rt + f1 = taudt*nu_constants::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; + + dd15*cos5*5.0e0_rt) - 0.5e0_rt*c16*xast*nu_constants::fac2*taudt; - f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0e0_rt + f2 = taudt*nu_constants::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; + + dd25*cos5*5.0e0_rt) - 0.5e0_rt*c26*xast*nu_constants::fac2*taudt; } @@ -996,12 +1132,12 @@ void sneut5(const Real temp, const Real den, sphotdz = sf.rm*sphotdz + sf.rmdz*a1; } - a1 = tfac4*(1.0e0_rt - tfac3 * qphot); + a1 = nu_constants::tfac4*(1.0e0_rt - nu_constants::tfac3 * qphot); a3 = sphot; sphot = a1*a3; if constexpr (do_derivatives) { - a2 = -tfac4*tfac3; + a2 = -nu_constants::tfac4*nu_constants::tfac3; sphotdt = a1*sphotdt + a2*qphotdt*a3; sphotda = a1*sphotda + a2*qphotda*a3; sphotdz = a1*sphotdz + a2*qphotdz*a3; @@ -1037,7 +1173,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*sf.ye, twoth))-1.0e0_rt); + tfermi = 5.9302e9_rt*(std::sqrt(1.0e0_rt+1.018e0_rt*std::pow(den6*sf.ye, nu_constants::twoth))-1.0e0_rt); // "weak" degenerate electrons only if (temp > 0.3e0_rt * tfermi) { @@ -1160,12 +1296,12 @@ void sneut5(const Real temp, const Real den, dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; } - z = tfac4*fbrem - tfac5*gbrem; + z = nu_constants::tfac4*fbrem - nu_constants::tfac5*gbrem; sbrem = dum * z; 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); + sbremdt = dumdt*z + dum*(nu_constants::tfac4*fbremdt - nu_constants::tfac5*gbremdt); + sbremda = dumda*z + dum*(nu_constants::tfac4*fbremda - nu_constants::tfac5*gbremda); + sbremdz = dumdz*z + dum*(nu_constants::tfac4*fbremdz - nu_constants::tfac5*gbremdz); } // liquid metal with c12 parameters (not too different for other elements) @@ -1173,7 +1309,7 @@ void sneut5(const Real temp, const Real den, } else { - u = fac3 * (std::log10(den) - 3.0e0_rt); + u = nu_constants::fac3 * (std::log10(den) - 3.0e0_rt); //a0 = iln10*fac3*sf.deni; // compute the expensive trig functions of equation 5.21 only once @@ -1249,26 +1385,26 @@ 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*sf.abari, oneth); + dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*sf.abari, nu_constants::oneth); if constexpr (do_derivatives) { dumdt = -dum*sf.tempi; - dumda = -oneth*dum*sf.abari; + dumda = -nu_constants::oneth*dum*sf.abari; dumdz = 2.0e0_rt*dum*sf.zbari; } gm1 = 1.0e0_rt/dum; gm2 = gm1*gm1; - gm13 = std::pow(gm1, oneth); + gm13 = std::pow(gm1, nu_constants::oneth); gm23 = gm13 * gm13; gm43 = gm13*gm1; gm53 = gm23*gm1; // equation 5.25 and 5.26 v = -0.05483e0_rt - 0.01946e0_rt*gm13 + 1.86310e0_rt*gm23 - 0.78873e0_rt*gm1; - a0 = oneth*0.01946e0_rt*gm43 - twoth*1.86310e0_rt*gm53 + 0.78873e0_rt*gm2; + a0 = nu_constants::oneth*0.01946e0_rt*gm43 - nu_constants::twoth*1.86310e0_rt*gm53 + 0.78873e0_rt*gm2; w = -0.06711e0_rt + 0.06859e0_rt*gm13 + 1.74360e0_rt*gm23 - 0.74498e0_rt*gm1; - a1 = -oneth*0.06859e0_rt*gm43 - twoth*1.74360e0_rt*gm53 + 0.74498e0_rt*gm2; + a1 = -nu_constants::oneth*0.06859e0_rt*gm43 - nu_constants::twoth*1.74360e0_rt*gm53 + 0.74498e0_rt*gm2; // equation 5.19 and 5.20 fliq = v*fb + (1.0e0_rt - v)*ft; @@ -1293,118 +1429,16 @@ void sneut5(const Real temp, const Real den, dumdz = 0.5738e0_rt*2.0e0_rt*sf.ye*t86*den; } - z = tfac4*fliq - tfac5*gliq; + z = nu_constants::tfac4*fliq - nu_constants::tfac5*gliq; sbrem = dum * z; 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); + sbremdt = dumdt*z + dum*(nu_constants::tfac4*fliqdt - nu_constants::tfac5*gliqdt); + sbremda = dumda*z + dum*(nu_constants::tfac4*fliqda - nu_constants::tfac5*gliqda); + sbremdz = dumdz*z + dum*(nu_constants::tfac4*fliqdz - nu_constants::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 * sf.ye /(temp*std::sqrt(temp)); - if constexpr (do_derivatives) { - xnumdt = -1.50e0_rt*xnum*sf.tempi; - xnumda = -xnum*sf.abari; - xnumdz = xnum*sf.zbari; - } - - // the chemical potential - nu = ifermi12(xnum); - - // a0 is d(nu)/d(xnum) - a0 = 1.0e0_rt/(0.5e0_rt*zfermim12(nu)); - if constexpr (do_derivatives) { - nudt = a0*xnumdt; - nuda = a0*xnumda; - nudz = a0*xnumdz; - } - nu2 = nu * nu; - nu3 = nu2 * nu; - - // table 12 - if (nu >= -20.0_rt && nu < 0.0_rt) { - a1 = 1.51e-2_rt; - a2 = 2.42e-1_rt; - a3 = 1.21e0_rt; - b = 3.71e-2_rt; - c = 9.06e-1_rt; - d = 9.28e-1_rt; - f1 = 0.0e0_rt; - f2 = 0.0e0_rt; - f3 = 0.0e0_rt; - } else if (nu >= 0.0_rt && nu <= 10.0_rt) { - a1 = 1.23e-2_rt; - a2 = 2.66e-1_rt; - a3 = 1.30e0_rt; - b = 1.17e-1_rt; - c = 8.97e-1_rt; - d = 1.77e-1_rt; - f1 = -1.20e-2_rt; - f2 = 2.29e-2_rt; - f3 = -1.04e-3_rt; - } - - // equation 6.7, 6.13 and 6.14 - if (nu >= -20.0_rt && nu <= 10.0_rt) { - - zeta = 1.579e5_rt*zbar*zbar*sf.tempi; - if constexpr (do_derivatives) { - zetadt = -zeta*sf.tempi; - zetada = 0.0e0_rt; - zetadz = 2.0e0_rt*zeta*sf.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; - 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); - dd01 = std::pow(dum, -4.55_rt); - c00 = a1*z + a2*dd00 + a3*dd01; - c01 = -(a1*z + 2.25_rt*a2*dd00 + 4.55_rt*a3*dd01)*z; - - z = std::exp(c*nu); - dd00 = b*z*(1.0e0_rt + d*dum); - gum = 1.0e0_rt + dd00; - 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; - 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); - dum = 1.0e0_rt + z; - a1 = 1.0e0_rt/dum; - a2 = 1.0e0_rt/bigj; - - 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*sf.abari + bigjda*a2 - z*(zetada+nuda)*a1); - srecodz = sreco*(14.0e0_rt*sf.zbari + bigjdz*a2 - z*(zetadz+nudz)*a1); - } - } + nu_recomb(sf, sreco, srecodt, srecoda, srecodz); // convert from erg/cm^3/s to erg/g/s // comment these out to duplicate the itoh et al plots From df21ae5a4491625d1d6835f1b2c7cafb4e1330fc Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 28 Oct 2023 19:21:47 -0400 Subject: [PATCH 06/10] move sneut bremsstrahlung into its own function (#1374) --- neutrinos/sneut5.H | 606 +++++++++++++++++++++++---------------------- 1 file changed, 306 insertions(+), 300 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index e6a8f7aedc..5b48671608 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -35,7 +35,6 @@ namespace nu_constants { constexpr Real tfac5 = 0.5e0_rt * tfac2; constexpr Real tfac6 = cv*cv + 1.5e0_rt*ca*ca + (xnufam - 1.0e0_rt)*(cvp*cvp + 1.5e0_rt*cap*cap); - } @@ -344,6 +343,309 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_brem(const sneutf_t& sf, + Real& sbrem, Real& sbremdt, Real& sbremda, Real& sbremdz) { + + // bremsstrahlung neutrino section + // for reactions like e- + (z,a) => e- + (z,a) + nu + nubar + // n + n => n + n + nu + nubar + // n + p => n + p + nu + nubar + // equation 4.3 + + Real den6 = sf.den * 1.0e-6_rt; + Real t8 = sf.temp * 1.0e-8_rt; + Real t812 = std::sqrt(t8); + Real t832 = t8 * t812; + Real t82 = t8*t8; + Real t83 = t82*t8; + Real t85 = t82*t83; + Real t86 = t85*t8; + Real t8m1 = 1.0e0_rt/t8; + Real t8m2 = t8m1*t8m1; + Real t8m3 = t8m2*t8m1; + Real t8m5 = t8m3*t8m2; + Real t8m6 = t8m5*t8m1; + + Real tfermi = 5.9302e9_rt * (std::sqrt(1.0e0_rt + 1.018e0_rt * std::pow(den6 * sf.ye, nu_constants::twoth)) - 1.0e0_rt); + + // "weak" degenerate electrons only + if (sf.temp > 0.3e0_rt * tfermi) { + + // equation 5.3 + Real dum = 7.05e6_rt * t832 + 5.12e4_rt * t83; + Real dumdt; + if constexpr (do_derivatives) { + dumdt = (1.5e0_rt * 7.05e6_rt * t812 + 3.0e0_rt * 5.12e4_rt * t82) * 1.0e-8_rt; + } + + Real z = 1.0e0_rt / dum; + Real eta = sf.rm * z; + Real etadt, etada, etadz; + if constexpr (do_derivatives) { + etadt = -sf.rm*z*z*dumdt; + etada = sf.rmda*z; + etadz = sf.rmdz*z; + } + + Real etam1 = 1.0e0_rt/eta; + Real etam2 = etam1 * etam1; + Real etam3 = etam2 * etam1; + + // equation 5.2 + Real a0 = 23.5e0_rt + 6.83e4_rt * t8m2 + 7.81e8_rt * t8m5; + Real f0 = (-2.0e0_rt * 6.83e4_rt * t8m3 - 5.0e0_rt * 7.81e8_rt * t8m6) * 1.0e-8_rt; + Real xnum = 1.0e0_rt / a0; + + dum = 1.0e0_rt + 1.47e0_rt * etam1 + 3.29e-2_rt * etam2; + Real dumda, dumdz; + 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; + } + + Real c00 = 1.26e0_rt * (1.0e0_rt + etam1); + Real c01, c02, c03, c04; + if constexpr (do_derivatives) { + z = -1.26e0_rt*etam2; + c01 = z*etadt; + c03 = z*etada; + c04 = z*etadz; + } + + z = 1.0e0_rt/dum; + Real xden = c00 * z; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + xdendt = (c01 - xden * dumdt) * z; + xdenda = (c03 - xden * dumda) * z; + xdendz = (c04 - xden * dumdz) * z; + } + + Real fbrem = xnum + xden; + Real fbremdt, fbremda, fbremdz; + 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; + f0 = (-2.0e0_rt * 6.7e5_rt * t8m3 - 5.0e0_rt * 7.66e9_rt * t8m6) * 1.0e-8_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 * sf.rmda; + dumdz = z * sf.rmdz; + } + + xnum = 1.0e0_rt / dum; + Real xnumdt, xnumda, xnumdz; + 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); + c01 = 4.07e0_rt + 0.0240e0_rt * std::pow(t8, 1.4e0_rt); + c02 = 4.59e-5_rt * std::pow(t8, -0.110e0_rt); + + Real dd00, dd01, dd02; + 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(sf.den, 0.656e0_rt); + dum = c00 * sf.rmi + c01 + c02 * z; + if constexpr (do_derivatives) { + 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; + if constexpr (do_derivatives) { + z = -xden * xden; + xdendt = z * dumdt; + xdenda = z * dumda; + xdendz = z * dumdz; + } + + Real gbrem = xnum + xden; + Real gbremdt, gbremda, gbremdz; + if constexpr (do_derivatives) { + gbremdt = xnumdt + xdendt; + gbremda = xnumda + xdenda; + gbremdz = xnumdz + xdendz; + } + + // equation 5.1 + dum = 0.5738e0_rt * sf.zbar * sf.ye * t86 * sf.den; + if constexpr (do_derivatives) { + dumdt = 0.5738e0_rt * sf.zbar * sf.ye * 6.0e0_rt * t85 * sf.den * 1.0e-8_rt; + dumda = -dum * sf.abari; + dumdz = 0.5738e0_rt * 2.0e0_rt * sf.ye * t86 * sf.den; + } + + z = nu_constants::tfac4 * fbrem - nu_constants::tfac5 * gbrem; + sbrem = dum * z; + if constexpr (do_derivatives) { + sbremdt = dumdt * z + dum * (nu_constants::tfac4 * fbremdt - nu_constants::tfac5 * gbremdt); + sbremda = dumda * z + dum * (nu_constants::tfac4 * fbremda - nu_constants::tfac5 * gbremda); + sbremdz = dumdz * z + dum * (nu_constants::tfac4 * fbremdz - nu_constants::tfac5 * gbremdz); + } + + // liquid metal with c12 parameters (not too different for other elements) + // equation 5.18 and 5.16 + + } else { + + Real u = nu_constants::fac3 * (std::log10(sf.den) - 3.0e0_rt); + //a0 = iln10*fac3*sf.deni; + + // compute the expensive trig functions of equation 5.21 only once + Real cos1 = std::cos(u); + Real cos2 = std::cos(2.0e0_rt*u); + Real cos3 = std::cos(3.0e0_rt*u); + Real cos4 = std::cos(4.0e0_rt*u); + Real cos5 = std::cos(5.0e0_rt*u); + + Real sin1 = std::sin(u); + Real sin2 = std::sin(2.0e0_rt*u); + Real sin3 = std::sin(3.0e0_rt*u); + Real sin4 = std::sin(4.0e0_rt*u); + //sin5 = std::sin(5.0e0_rt*u); + + // equation 5.21 + Real fb = 0.5e0_rt * 0.17946e0_rt + 0.00945e0_rt * u + 0.34529e0_rt + - 0.05821e0_rt * cos1 - 0.04969e0_rt * sin1 + - 0.01089e0_rt * cos2 - 0.01584e0_rt * sin2 + - 0.01147e0_rt * cos3 - 0.00504e0_rt * sin3 + - 0.00656e0_rt * cos4 - 0.00281e0_rt * sin4 + - 0.00519e0_rt * cos5; + + // c00 = a0*(0.00945e0_rt + // + 0.05821e0_rt*sin1 - 0.04969e0_rt*cos1 + // + 0.01089e0_rt*sin2*2.0e0_rt - 0.01584e0_rt*cos2*2.0e0_rt + // + 0.01147e0_rt*sin3*3.0e0_rt - 0.00504e0_rt*cos3*3.0e0_rt + // + 0.00656e0_rt*sin4*4.0e0_rt - 0.00281e0_rt*cos4*4.0e0_rt + // + 0.00519e0_rt*sin5*5.0e0_rt); + + // equation 5.22 + Real ft = 0.5e0_rt * 0.06781e0_rt - 0.02342e0_rt * u + 0.24819e0_rt + - 0.00944e0_rt * cos1 - 0.02213e0_rt * sin1 + - 0.01289e0_rt * cos2 - 0.01136e0_rt * sin2 + - 0.00589e0_rt * cos3 - 0.00467e0_rt * sin3 + - 0.00404e0_rt * cos4 - 0.00131e0_rt * sin4 + - 0.00330e0_rt * cos5; + + // c01 = a0*(-0.02342e0_rt + // + 0.00944e0_rt*sin1 - 0.02213e0_rt*cos1 + // + 0.01289e0_rt*sin2*2.0e0_rt - 0.01136e0_rt*cos2*2.0e0_rt + // + 0.00589e0_rt*sin3*3.0e0_rt - 0.00467e0_rt*cos3*3.0e0_rt + // + 0.00404e0_rt*sin4*4.0e0_rt - 0.00131e0_rt*cos4*4.0e0_rt + // + 0.00330e0_rt*sin5*5.0e0_rt); + + // equation 5.23 + Real gb = 0.5e0_rt * 0.00766e0_rt - 0.01259e0_rt * u + 0.07917e0_rt + - 0.00710e0_rt * cos1 + 0.02300e0_rt * sin1 + - 0.00028e0_rt * cos2 - 0.01078e0_rt * sin2 + + 0.00232e0_rt * cos3 + 0.00118e0_rt * sin3 + + 0.00044e0_rt * cos4 - 0.00089e0_rt * sin4 + + 0.00158e0_rt * cos5; + + // c02 = a0*(-0.01259e0_rt + // + 0.00710e0_rt*sin1 + 0.02300e0_rt*cos1 + // + 0.00028e0_rt*sin2*2.0e0_rt - 0.01078e0_rt*cos2*2.0e0_rt + // - 0.00232e0_rt*sin3*3.0e0_rt + 0.00118e0_rt*cos3*3.0e0_rt + // - 0.00044e0_rt*sin4*4.0e0_rt - 0.00089e0_rt*cos4*4.0e0_rt + // - 0.00158e0_rt*sin5*5.0e0_rt); + + // equation 5.24 + Real gt = -0.5e0_rt * 0.00769e0_rt - 0.00829e0_rt * u + 0.05211e0_rt + + 0.00356e0_rt * cos1 + 0.01052e0_rt * sin1 + - 0.00184e0_rt * cos2 - 0.00354e0_rt * sin2 + + 0.00146e0_rt * cos3 - 0.00014e0_rt * sin3 + + 0.00031e0_rt * cos4 - 0.00018e0_rt * sin4 + + 0.00069e0_rt * cos5; + + // c03 = a0*(-0.00829e0_rt + // - 0.00356e0_rt*sin1 + 0.01052e0_rt*cos1 + // + 0.00184e0_rt*sin2*2.0e0_rt - 0.00354e0_rt*cos2*2.0e0_rt + // - 0.00146e0_rt*sin3*3.0e0_rt - 0.00014e0_rt*cos3*3.0e0_rt + // - 0.00031e0_rt*sin4*4.0e0_rt - 0.00018e0_rt*cos4*4.0e0_rt + // - 0.00069e0_rt*sin5*5.0e0_rt); + + Real dum = 2.275e-1_rt * sf.zbar * sf.zbar * t8m1 * std::pow(den6*sf.abari, nu_constants::oneth); + Real dumdt, dumda, dumdz; + if constexpr (do_derivatives) { + dumdt = -dum*sf.tempi; + dumda = -nu_constants::oneth*dum*sf.abari; + dumdz = 2.0e0_rt*dum*sf.zbari; + } + + Real gm1 = 1.0e0_rt / dum; + Real gm2 = gm1 * gm1; + Real gm13 = std::pow(gm1, nu_constants::oneth); + Real gm23 = gm13 * gm13; + Real gm43 = gm13 * gm1; + Real gm53 = gm23 * gm1; + + // equation 5.25 and 5.26 + Real v = -0.05483e0_rt - 0.01946e0_rt * gm13 + 1.86310e0_rt * gm23 - 0.78873e0_rt * gm1; + Real a0 = nu_constants::oneth * 0.01946e0_rt * gm43 - nu_constants::twoth * 1.86310e0_rt * gm53 + 0.78873e0_rt * gm2; + + Real w = -0.06711e0_rt + 0.06859e0_rt * gm13 + 1.74360e0_rt * gm23 - 0.74498e0_rt * gm1; + Real a1 = -nu_constants::oneth*0.06859e0_rt * gm43 - nu_constants::twoth * 1.74360e0_rt * gm53 + 0.74498e0_rt * gm2; + + // equation 5.19 and 5.20 + Real fliq = v * fb + (1.0e0_rt - v) * ft; + Real fliqdt, fliqda, fliqdz; + if constexpr (do_derivatives) { + fliqdt = a0 * dumdt * (fb - ft); + fliqda = a0 * dumda * (fb - ft); + fliqdz = a0 * dumdz * (fb - ft); + } + + Real gliq = w * gb + (1.0e0_rt - w) * gt; + Real gliqdt, gliqda, gliqdz; + 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 * sf.zbar * sf.ye * t86 * sf.den; + if constexpr (do_derivatives) { + dumdt = 0.5738e0_rt * sf.zbar * sf.ye * 6.0e0_rt * t85 * sf.den * 1.0e-8_rt; + dumda = -dum * sf.abari; + dumdz = 0.5738e0_rt * 2.0e0_rt * sf.ye * t86 * sf.den; + } + + Real z = nu_constants::tfac4*fliq - nu_constants::tfac5*gliq; + sbrem = dum * z; + if constexpr (do_derivatives) { + sbremdt = dumdt*z + dum*(nu_constants::tfac4*fliqdt - nu_constants::tfac5*gliqdt); + sbremda = dumda*z + dum*(nu_constants::tfac4*fliqda - nu_constants::tfac5*gliqda); + sbremdz = dumdz*z + dum*(nu_constants::tfac4*fliqdz - nu_constants::tfac5*gliqdz); + } + } +} template AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -494,10 +796,10 @@ void sneut5(const Real temp, const Real den, dsnudz = derivative of snu with zbar */ - Real xlnt,cc,den6,tfermi, + Real xlnt,cc, 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, + c25,c26,dd01,dd02,dd03,dd04,dd05,dd11,dd12, dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,c,d{0.0},f0, f1{0.0},f2,z, dum,dumdt,gum,dumda,dumdz; @@ -523,18 +825,6 @@ void sneut5(const Real temp, const Real den, fphot,fphotdt,qphot,qphotdt, fphotda,fphotdz,qphotda,qphotdz; - // brem - Real t8,t812,t832,t82,t83,t85,t86,t8m1,t8m2,t8m3,t8m5, - t8m6, - eta,etadt,etam1,etam2,etam3, - fbrem,fbremdt, - gbrem,gbremdt, - u,gm1,gm2,gm13,gm23,gm43,gm53,v,w,fb,gt,gb, - fliq,fliqdt,gliq,gliqdt, - etada,etadz, - fbremda,fbremdz,gbremda,gbremdz, - fliqda,fliqdz,gliqda,gliqdz; - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; @@ -1152,291 +1442,7 @@ void sneut5(const Real temp, const Real den, } } - // bremsstrahlung neutrino section - // for reactions like e- + (z,a) => e- + (z,a) + nu + nubar - // n + n => n + n + nu + nubar - // n + p => n + p + nu + nubar - // equation 4.3 - - den6 = den * 1.0e-6_rt; - t8 = temp * 1.0e-8_rt; - t812 = std::sqrt(t8); - t832 = t8 * t812; - t82 = t8*t8; - t83 = t82*t8; - t85 = t82*t83; - t86 = t85*t8; - t8m1 = 1.0e0_rt/t8; - t8m2 = t8m1*t8m1; - t8m3 = t8m2*t8m1; - t8m5 = t8m3*t8m2; - t8m6 = t8m5*t8m1; - - - tfermi = 5.9302e9_rt*(std::sqrt(1.0e0_rt+1.018e0_rt*std::pow(den6*sf.ye, nu_constants::twoth))-1.0e0_rt); - - // "weak" degenerate electrons only - if (temp > 0.3e0_rt * tfermi) { - - // equation 5.3 - dum = 7.05e6_rt * t832 + 5.12e4_rt * t83; - 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 = sf.rm*z; - if constexpr (do_derivatives) { - etadt = -sf.rm*z*z*dumdt; - etada = sf.rmda*z; - etadz = sf.rmdz*z; - } - - etam1 = 1.0e0_rt/eta; - etam2 = etam1 * etam1; - etam3 = etam2 * etam1; - - // equation 5.2 - a0 = 23.5e0_rt + 6.83e4_rt*t8m2 + 7.81e8_rt*t8m5; - f0 = (-2.0e0_rt*6.83e4_rt*t8m3 - 5.0e0_rt*7.81e8_rt*t8m6)*1.0e-8_rt; - xnum = 1.0e0_rt/a0; - - dum = 1.0e0_rt + 1.47e0_rt*etam1 + 3.29e-2_rt*etam2; - 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); - 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; - if constexpr (do_derivatives) { - xdendt = (c01 - xden*dumdt)*z; - xdenda = (c03 - xden*dumda)*z; - xdendz = (c04 - xden*dumdz)*z; - } - - fbrem = xnum + xden; - 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; - f0 = (-2.0e0_rt*6.7e5_rt*t8m3 - 5.0e0_rt*7.66e9_rt*t8m6)*1.0e-8_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*sf.rmda; - dumdz = z*sf.rmdz; - } - - xnum = 1.0e0_rt/dum; - 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); - c01 = 4.07e0_rt + 0.0240e0_rt * std::pow(t8, 1.4e0_rt); - c02 = 4.59e-5_rt * std::pow(t8, -0.110e0_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*sf.rmi + c01 + c02*z; - if constexpr (do_derivatives) { - 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; - if constexpr (do_derivatives) { - z = -xden*xden; - xdendt = z*dumdt; - xdenda = z*dumda; - xdendz = z*dumdz; - } - - gbrem = xnum + xden; - if constexpr (do_derivatives) { - gbremdt = xnumdt + xdendt; - gbremda = xnumda + xdenda; - gbremdz = xnumdz + xdendz; - } - - // equation 5.1 - dum = 0.5738e0_rt*zbar*sf.ye*t86*den; - if constexpr (do_derivatives) { - 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 = nu_constants::tfac4*fbrem - nu_constants::tfac5*gbrem; - sbrem = dum * z; - if constexpr (do_derivatives) { - sbremdt = dumdt*z + dum*(nu_constants::tfac4*fbremdt - nu_constants::tfac5*gbremdt); - sbremda = dumda*z + dum*(nu_constants::tfac4*fbremda - nu_constants::tfac5*gbremda); - sbremdz = dumdz*z + dum*(nu_constants::tfac4*fbremdz - nu_constants::tfac5*gbremdz); - } - - // liquid metal with c12 parameters (not too different for other elements) - // equation 5.18 and 5.16 - - } else { - - u = nu_constants::fac3 * (std::log10(den) - 3.0e0_rt); - //a0 = iln10*fac3*sf.deni; - - // compute the expensive trig functions of equation 5.21 only once - cos1 = std::cos(u); - cos2 = std::cos(2.0e0_rt*u); - cos3 = std::cos(3.0e0_rt*u); - cos4 = std::cos(4.0e0_rt*u); - cos5 = std::cos(5.0e0_rt*u); - - sin1 = std::sin(u); - sin2 = std::sin(2.0e0_rt*u); - sin3 = std::sin(3.0e0_rt*u); - sin4 = std::sin(4.0e0_rt*u); - //sin5 = std::sin(5.0e0_rt*u); - - // equation 5.21 - fb = 0.5e0_rt * 0.17946e0_rt + 0.00945e0_rt*u + 0.34529e0_rt - - 0.05821e0_rt*cos1 - 0.04969e0_rt*sin1 - - 0.01089e0_rt*cos2 - 0.01584e0_rt*sin2 - - 0.01147e0_rt*cos3 - 0.00504e0_rt*sin3 - - 0.00656e0_rt*cos4 - 0.00281e0_rt*sin4 - - 0.00519e0_rt*cos5; - - // c00 = a0*(0.00945e0_rt - // + 0.05821e0_rt*sin1 - 0.04969e0_rt*cos1 - // + 0.01089e0_rt*sin2*2.0e0_rt - 0.01584e0_rt*cos2*2.0e0_rt - // + 0.01147e0_rt*sin3*3.0e0_rt - 0.00504e0_rt*cos3*3.0e0_rt - // + 0.00656e0_rt*sin4*4.0e0_rt - 0.00281e0_rt*cos4*4.0e0_rt - // + 0.00519e0_rt*sin5*5.0e0_rt); - - // equation 5.22 - ft = 0.5e0_rt * 0.06781e0_rt - 0.02342e0_rt*u + 0.24819e0_rt - - 0.00944e0_rt*cos1 - 0.02213e0_rt*sin1 - - 0.01289e0_rt*cos2 - 0.01136e0_rt*sin2 - - 0.00589e0_rt*cos3 - 0.00467e0_rt*sin3 - - 0.00404e0_rt*cos4 - 0.00131e0_rt*sin4 - - 0.00330e0_rt*cos5; - - // c01 = a0*(-0.02342e0_rt - // + 0.00944e0_rt*sin1 - 0.02213e0_rt*cos1 - // + 0.01289e0_rt*sin2*2.0e0_rt - 0.01136e0_rt*cos2*2.0e0_rt - // + 0.00589e0_rt*sin3*3.0e0_rt - 0.00467e0_rt*cos3*3.0e0_rt - // + 0.00404e0_rt*sin4*4.0e0_rt - 0.00131e0_rt*cos4*4.0e0_rt - // + 0.00330e0_rt*sin5*5.0e0_rt); - - // equation 5.23 - gb = 0.5e0_rt * 0.00766e0_rt - 0.01259e0_rt*u + 0.07917e0_rt - - 0.00710e0_rt*cos1 + 0.02300e0_rt*sin1 - - 0.00028e0_rt*cos2 - 0.01078e0_rt*sin2 - + 0.00232e0_rt*cos3 + 0.00118e0_rt*sin3 - + 0.00044e0_rt*cos4 - 0.00089e0_rt*sin4 - + 0.00158e0_rt*cos5; - - // c02 = a0*(-0.01259e0_rt - // + 0.00710e0_rt*sin1 + 0.02300e0_rt*cos1 - // + 0.00028e0_rt*sin2*2.0e0_rt - 0.01078e0_rt*cos2*2.0e0_rt - // - 0.00232e0_rt*sin3*3.0e0_rt + 0.00118e0_rt*cos3*3.0e0_rt - // - 0.00044e0_rt*sin4*4.0e0_rt - 0.00089e0_rt*cos4*4.0e0_rt - // - 0.00158e0_rt*sin5*5.0e0_rt); - - // equation 5.24 - gt = -0.5e0_rt * 0.00769e0_rt - 0.00829e0_rt*u + 0.05211e0_rt - + 0.00356e0_rt*cos1 + 0.01052e0_rt*sin1 - - 0.00184e0_rt*cos2 - 0.00354e0_rt*sin2 - + 0.00146e0_rt*cos3 - 0.00014e0_rt*sin3 - + 0.00031e0_rt*cos4 - 0.00018e0_rt*sin4 - + 0.00069e0_rt*cos5; - - // c03 = a0*(-0.00829e0_rt - // - 0.00356e0_rt*sin1 + 0.01052e0_rt*cos1 - // + 0.00184e0_rt*sin2*2.0e0_rt - 0.00354e0_rt*cos2*2.0e0_rt - // - 0.00146e0_rt*sin3*3.0e0_rt - 0.00014e0_rt*cos3*3.0e0_rt - // - 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*sf.abari, nu_constants::oneth); - if constexpr (do_derivatives) { - dumdt = -dum*sf.tempi; - dumda = -nu_constants::oneth*dum*sf.abari; - dumdz = 2.0e0_rt*dum*sf.zbari; - } - - gm1 = 1.0e0_rt/dum; - gm2 = gm1*gm1; - gm13 = std::pow(gm1, nu_constants::oneth); - gm23 = gm13 * gm13; - gm43 = gm13*gm1; - gm53 = gm23*gm1; - - // equation 5.25 and 5.26 - v = -0.05483e0_rt - 0.01946e0_rt*gm13 + 1.86310e0_rt*gm23 - 0.78873e0_rt*gm1; - a0 = nu_constants::oneth*0.01946e0_rt*gm43 - nu_constants::twoth*1.86310e0_rt*gm53 + 0.78873e0_rt*gm2; - - w = -0.06711e0_rt + 0.06859e0_rt*gm13 + 1.74360e0_rt*gm23 - 0.74498e0_rt*gm1; - a1 = -nu_constants::oneth*0.06859e0_rt*gm43 - nu_constants::twoth*1.74360e0_rt*gm53 + 0.74498e0_rt*gm2; - - // equation 5.19 and 5.20 - fliq = v*fb + (1.0e0_rt - v)*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; - 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*sf.ye*t86*den; - if constexpr (do_derivatives) { - 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 = nu_constants::tfac4*fliq - nu_constants::tfac5*gliq; - sbrem = dum * z; - if constexpr (do_derivatives) { - sbremdt = dumdt*z + dum*(nu_constants::tfac4*fliqdt - nu_constants::tfac5*gliqdt); - sbremda = dumda*z + dum*(nu_constants::tfac4*fliqda - nu_constants::tfac5*gliqda); - sbremdz = dumdz*z + dum*(nu_constants::tfac4*fliqdz - nu_constants::tfac5*gliqdz); - } - } + nu_brem(sf, sbrem, sbremdt, sbremda, sbremdz); nu_recomb(sf, sreco, srecodt, srecoda, srecodz); From 9066232c3a6b3c1e94f1de2020d879a61cdd7cf7 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Sun, 29 Oct 2023 15:22:16 -0400 Subject: [PATCH 07/10] add weak rate neutrino loss for self-consistent nse check (#1369) --- nse_solver/nse_check.H | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H index 3ec6fd2017..fb7457ee63 100644 --- a/nse_solver/nse_check.H +++ b/nse_solver/nse_check.H @@ -1065,6 +1065,9 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { amrex::Real enuc; ener_gener_rate(ydot, enuc); + // include any weak rate neutrino losses + enuc += rate_eval.enuc_weak; + #ifdef NEUTRINOS // get abar and zbar composition(state); From 53d979817c87aeeb0c7007bc20de7fdadf2b7b63 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 16:25:50 -0400 Subject: [PATCH 08/10] move sneut photo section into its own function (#1375) --- neutrinos/sneut5.H | 683 +++++++++++++++++++++++---------------------- 1 file changed, 357 insertions(+), 326 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 5b48671608..38a66f29bd 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -291,9 +291,16 @@ struct sneutf_t { Real rmdz; Real rmi; + Real zeta; + Real zeta2; + Real zeta3; + Real zetadt; + Real zetada; + Real zetadz; }; +template AMREX_GPU_HOST_DEVICE inline sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { @@ -339,10 +346,342 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { sf.rmdz = den * sf.abari; sf.rmi = 1.0e0_rt / sf.rm; + Real a0 = sf.rm * 1.0e-9_rt; + Real a1 = std::pow(a0, nu_constants::oneth); + sf.zeta = a1 * sf.xlm1; + + if constexpr (do_derivatives) { + Real a2 = nu_constants::oneth * a1 * sf.rmi * sf.xlm1; + sf.zetadt = -a1 * sf.xlm2 * sf.xldt; + sf.zetada = a2 * sf.rmda; + sf.zetadz = a2 * sf.rmdz; + } + + sf.zeta2 = sf.zeta * sf.zeta; + sf.zeta3 = sf.zeta2 * sf.zeta; + return sf; } +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_photo(const sneutf_t& sf, + Real& sphot, Real& sphotdt, Real& sphotda, Real& sphotdz) { + + // photoneutrino process section + // for reactions like e- + gamma => e- + nu_e + nubar_e + // e+ + gamma => e+ + nu_e + nubar_e + // equation 3.8 for tau, equation 3.6 for cc, + // and table 2 written out for speed + + Real tau; + Real cc; + Real c00, c01, c02, c03, c04, c05, c06; + Real c10, c11, c12, c13, c14, c15, c16; + Real c20, c21, c22, c23, c24, c25, c26; + Real dd01, dd02, dd03, dd04, dd05; + Real dd11, dd12, dd13, dd14, dd15; + Real dd21, dd22, dd23, dd24, dd25; + + if (sf.temp < 1.0e8_rt) { + + // note: we already bailed above for T < 1.e7, so this is really 1.e7 <= T < 1.e8 + + tau = std::log10(sf.temp * 1.0e-7_rt); + cc = 0.5654e0_rt + tau; + c00 = 1.008e11_rt; + c01 = 0.0e0_rt; + c02 = 0.0e0_rt; + c03 = 0.0e0_rt; + c04 = 0.0e0_rt; + c05 = 0.0e0_rt; + c06 = 0.0e0_rt; + c10 = 8.156e10_rt; + c11 = 9.728e8_rt; + c12 = -3.806e9_rt; + c13 = -4.384e9_rt; + c14 = -5.774e9_rt; + c15 = -5.249e9_rt; + c16 = -5.153e9_rt; + c20 = 1.067e11_rt; + c21 = -9.782e9_rt; + c22 = -7.193e9_rt; + c23 = -6.936e9_rt; + c24 = -6.893e9_rt; + c25 = -7.041e9_rt; + c26 = -7.193e9_rt; + dd01 = 0.0e0_rt; + dd02 = 0.0e0_rt; + dd03 = 0.0e0_rt; + dd04 = 0.0e0_rt; + dd05 = 0.0e0_rt; + dd11 = -1.879e10_rt; + dd12 = -9.667e9_rt; + dd13 = -5.602e9_rt; + dd14 = -3.370e9_rt; + dd15 = -1.825e9_rt; + dd21 = -2.919e10_rt; + dd22 = -1.185e10_rt; + dd23 = -7.270e9_rt; + dd24 = -4.222e9_rt; + dd25 = -1.560e9_rt; + + } else if (sf.temp >= 1.0e8_rt && sf.temp < 1.0e9_rt) { + + tau = std::log10(sf.temp * 1.0e-8_rt); + cc = 1.5654e0_rt; + c00 = 9.889e10_rt; + c01 = -4.524e8_rt; + c02 = -6.088e6_rt; + c03 = 4.269e7_rt; + c04 = 5.172e7_rt; + c05 = 4.910e7_rt; + c06 = 4.388e7_rt; + c10 = 1.813e11_rt; + c11 = -7.556e9_rt; + c12 = -3.304e9_rt; + c13 = -1.031e9_rt; + c14 = -1.764e9_rt; + c15 = -1.851e9_rt; + c16 = -1.928e9_rt; + c20 = 9.750e10_rt; + c21 = 3.484e10_rt; + c22 = 5.199e9_rt; + c23 = -1.695e9_rt; + c24 = -2.865e9_rt; + c25 = -3.395e9_rt; + c26 = -3.418e9_rt; + dd01 = -1.135e8_rt; + dd02 = 1.256e8_rt; + dd03 = 5.149e7_rt; + dd04 = 3.436e7_rt; + dd05 = 1.005e7_rt; + dd11 = 1.652e9_rt; + dd12 = -3.119e9_rt; + dd13 = -1.839e9_rt; + dd14 = -1.458e9_rt; + dd15 = -8.956e8_rt; + dd21 = -1.548e10_rt; + dd22 = -9.338e9_rt; + dd23 = -5.899e9_rt; + dd24 = -3.035e9_rt; + dd25 = -1.598e9_rt; + + } else { + + // T > 1.e9 + + tau = std::log10(sf.t9); + cc = 1.5654e0_rt; + c00 = 9.581e10_rt; + c01 = 4.107e8_rt; + c02 = 2.305e8_rt; + c03 = 2.236e8_rt; + c04 = 1.580e8_rt; + c05 = 2.165e8_rt; + c06 = 1.721e8_rt; + c10 = 1.459e12_rt; + c11 = 1.314e11_rt; + c12 = -1.169e11_rt; + c13 = -1.765e11_rt; + c14 = -1.867e11_rt; + c15 = -1.983e11_rt; + c16 = -1.896e11_rt; + c20 = 2.424e11_rt; + c21 = -3.669e9_rt; + c22 = -8.691e9_rt; + c23 = -7.967e9_rt; + c24 = -7.932e9_rt; + c25 = -7.987e9_rt; + c26 = -8.333e9_rt; + dd01 = 4.724e8_rt; + dd02 = 2.976e8_rt; + dd03 = 2.242e8_rt; + dd04 = 7.937e7_rt; + dd05 = 4.859e7_rt; + dd11 = -7.094e11_rt; + dd12 = -3.697e11_rt; + dd13 = -2.189e11_rt; + dd14 = -1.273e11_rt; + dd15 = -5.705e10_rt; + dd21 = -2.254e10_rt; + dd22 = -1.551e10_rt; + dd23 = -7.793e9_rt; + dd24 = -4.489e9_rt; + dd25 = -2.185e9_rt; + + } + + Real taudt = nu_constants::iln10 * sf.tempi; + + // equation 3.7, compute the expensive trig functions only one time + + Real cos1 = std::cos(nu_constants::fac1 * tau); + Real cos2 = std::cos(nu_constants::fac1 * 2.0e0_rt * tau); + Real cos3 = std::cos(nu_constants::fac1 * 3.0e0_rt * tau); + Real cos4 = std::cos(nu_constants::fac1 * 4.0e0_rt * tau); + Real cos5 = std::cos(nu_constants::fac1 * 5.0e0_rt * tau); + Real last = std::cos(nu_constants::fac2 * tau); + + Real sin1 = std::sin(nu_constants::fac1*tau); + Real sin2 = std::sin(nu_constants::fac1*2.0e0_rt*tau); + Real sin3 = std::sin(nu_constants::fac1*3.0e0_rt*tau); + Real sin4 = std::sin(nu_constants::fac1*4.0e0_rt*tau); + Real sin5 = std::sin(nu_constants::fac1*5.0e0_rt*tau); + Real xast = std::sin(nu_constants::fac2*tau); + + Real a0 = 0.5e0_rt * c00 + + 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; + Real 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; + Real 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; + + Real f0, f1, f2; + if constexpr (do_derivatives) { + f0 = taudt * nu_constants::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 * nu_constants::fac2 * taudt; + + f1 = taudt * nu_constants::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*nu_constants::fac2*taudt; + + f2 = taudt * nu_constants::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 * nu_constants::fac2 * taudt; + } + + // equation 3.4 + Real dum = a0 + a1 * sf.zeta + a2 * sf.zeta2; + Real dumdt, dumda, dumdz; + if constexpr (do_derivatives) { + dumdt = f0 + f1 * sf.zeta + a1 * sf.zetadt + f2 * sf.zeta2 + 2.0e0_rt * a2 * sf.zeta * sf.zetadt; + dumda = a1 * sf.zetada + 2.0e0_rt * a2 * sf.zeta * sf.zetada; + dumdz = a1 * sf.zetadz + 2.0e0_rt * a2 * sf.zeta * sf.zetadz; + } + + Real z = std::exp(-cc * sf.zeta); + + Real xnum = dum * z; + Real xnumdt, xnumda, xnumdz; + if constexpr (do_derivatives) { + xnumdt = dumdt * z - dum * z * cc * sf.zetadt; + xnumda = dumda * z - dum * z * cc * sf.zetada; + xnumdz = dumdz * z - dum * z * cc * sf.zetadz; + } + + Real xden = sf.zeta3 + 6.290e-3_rt * sf.xlm1 + 7.483e-3_rt * sf.xlm2 + 3.061e-4_rt * sf.xlm3; + + dum = 3.0e0_rt * sf.zeta2; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + xdendt = dum * sf.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 * sf.zetada; + xdendz = dum * sf.zetadz; + } + dum = 1.0e0_rt / xden; + Real fphot = xnum * dum; + Real fphotdt, fphotda, fphotdz; + 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 * 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 * sf.xldt; + } + + 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 = 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 + sf.rm * z; + if constexpr (do_derivatives) { + xdendt = -sf.rm * z * z * dumdt; + xdenda = sf.rmda * z; + xdendz = sf.rmdz * z; + } + + z = 1.0e0_rt / xden; + Real qphot = xnum * z; + Real qphotdt; + if constexpr (do_derivatives) { + qphotdt = (xnumdt - qphot * xdendt) * z; + } + dum = -qphot * z; + Real qphotda, qphotdz; + if constexpr (do_derivatives) { + qphotda = dum * xdenda; + qphotdz = dum * xdendz; + } + + // equation 3.2 + sphot = sf.xl5 * fphot; + if constexpr (do_derivatives) { + sphotdt = 5.0e0_rt * sf.xl4 * sf.xldt * fphot + sf.xl5 * fphotdt; + sphotda = sf.xl5 * fphotda; + sphotdz = sf.xl5 * fphotdz; + } + + a1 = sphot; + sphot = sf.rm * a1; + if constexpr (do_derivatives) { + sphotdt = sf.rm * sphotdt; + sphotda = sf.rm * sphotda + sf.rmda * a1; + sphotdz = sf.rm * sphotdz + sf.rmdz * a1; + } + + a1 = nu_constants::tfac4 * (1.0e0_rt - nu_constants::tfac3 * qphot); + Real a3 = sphot; + sphot = a1 * a3; + if constexpr (do_derivatives) { + a2 = -nu_constants::tfac4 * nu_constants::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; + if constexpr (do_derivatives) { + sphotdt = 0.0e0_rt; + sphotda = 0.0e0_rt; + sphotdz = 0.0e0_rt; + } + } +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void nu_brem(const sneutf_t& sf, @@ -797,20 +1136,16 @@ void sneut5(const Real temp, const Real den, */ Real xlnt,cc, - 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,dd01,dd02,dd03,dd04,dd05,dd11,dd12, - dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,c,d{0.0},f0, - f1{0.0},f2,z, - dum,dumdt,gum,dumda,dumdz; + a1,a2,a3,b1,b2,c00,c01,c03,c04, + c,d{0.0}, + f1{0.0}, + dumdt,gum,dumda,dumdz; // pair production Real gl,gldt, - zeta,zetadt,zeta2,zeta3, xnum,xnumdt,xden,xdendt, fpair,fpairdt,qpair,qpairdt, fpairda,fpairdz,qpairda,qpairdz, - zetada,zetadz, xnumda,xnumdz,xdenda,xdendz; // plasma @@ -819,11 +1154,6 @@ void sneut5(const Real temp, const Real den, flda,fldz,ftda,ftdz, fxyda,fxydz,gl2da,gl2dz; - // photo - Real tau,taudt,cos1,cos2,cos3,cos4,cos5,sin1,sin2, - sin3,sin4,sin5,last,xast, - fphot,fphotdt,qphot,qphotdt, - fphotda,fphotdz,qphotda,qphotdz; // initialize Real spair{0.0e0_rt}; @@ -859,21 +1189,7 @@ void sneut5(const Real temp, const Real den, if (temp < 1.0e7_rt) return; - auto sf = get_sneut_factors(den, temp, abar, zbar); - - a0 = sf.rm * 1.0e-9_rt; - a1 = std::pow(a0, nu_constants::oneth); - zeta = a1 * sf.xlm1; - - if constexpr (do_derivatives) { - a2 = nu_constants::oneth * a1*sf.rmi * sf.xlm1; - zetadt = -a1 * sf.xlm2 * sf.xldt; - zetada = a2 * sf.rmda; - zetadz = a2 * sf.rmdz; - } - - zeta2 = zeta * zeta; - zeta3 = zeta2 * zeta; + auto sf = get_sneut_factors(den, temp, abar, zbar); // pair neutrino section // for reactions like e+ + e- => nu_e + nubar_e @@ -886,15 +1202,15 @@ void sneut5(const Real temp, const Real den, // equation 2.7 - a1 = 6.002e19_rt + 2.084e20_rt*zeta + 1.872e21_rt*zeta2; + a1 = 6.002e19_rt + 2.084e20_rt * sf.zeta + 1.872e21_rt * sf.zeta2; if (sf.t9 < 10.0_rt) { - b1 = std::exp(-5.5924e0_rt*zeta); + b1 = std::exp(-5.5924e0_rt * sf.zeta); if constexpr (do_derivatives) { b2 = -b1*5.5924e0_rt; } } else { - b1 = std::exp(-4.9924e0_rt*zeta); + b1 = std::exp(-4.9924e0_rt * sf.zeta); if constexpr (do_derivatives) { b2 = -b1*4.9924e0_rt; } @@ -902,11 +1218,11 @@ void sneut5(const Real temp, const Real den, xnum = a1 * b1; if constexpr (do_derivatives) { - a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta; + a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt * sf.zeta; c = a2*b1 + a1*b2; - xnumdt = c*zetadt; - xnumda = c*zetada; - xnumdz = c*zetadz; + xnumdt = c * sf.zetadt; + xnumda = c * sf.zetada; + xnumdz = c * sf.zetadz; } if (sf.t9 < 10.0_rt) { @@ -921,12 +1237,12 @@ void sneut5(const Real temp, const Real den, } } - xden = zeta3 + a1; + xden = sf.zeta3 + a1; if constexpr (do_derivatives) { - b1 = 3.0e0_rt*zeta2; - xdendt = b1*zetadt + a2*sf.xldt; - xdenda = b1*zetada; - xdendz = b1*zetadz; + b1 = 3.0e0_rt*sf.zeta2; + xdendt = b1*sf.zetadt + a2*sf.xldt; + xdenda = b1*sf.zetada; + xdendz = b1*sf.zetadz; } a1 = 1.0e0_rt/xden; @@ -1155,292 +1471,7 @@ void sneut5(const Real temp, const Real den, splasdz = a2*splasdz; } - // photoneutrino process section - // for reactions like e- + gamma => e- + nu_e + nubar_e - // e+ + gamma => e+ + nu_e + nubar_e - // equation 3.8 for tau, equation 3.6 for cc, - // and table 2 written out for speed - if (temp < 1.0e8_rt) { - - // note: we already bailed above for T < 1.e7, so this is really 1.e7 <= T < 1.e8 - - tau = std::log10(temp * 1.0e-7_rt); - cc = 0.5654e0_rt + tau; - c00 = 1.008e11_rt; - c01 = 0.0e0_rt; - c02 = 0.0e0_rt; - c03 = 0.0e0_rt; - c04 = 0.0e0_rt; - c05 = 0.0e0_rt; - c06 = 0.0e0_rt; - c10 = 8.156e10_rt; - c11 = 9.728e8_rt; - c12 = -3.806e9_rt; - c13 = -4.384e9_rt; - c14 = -5.774e9_rt; - c15 = -5.249e9_rt; - c16 = -5.153e9_rt; - c20 = 1.067e11_rt; - c21 = -9.782e9_rt; - c22 = -7.193e9_rt; - c23 = -6.936e9_rt; - c24 = -6.893e9_rt; - c25 = -7.041e9_rt; - c26 = -7.193e9_rt; - dd01 = 0.0e0_rt; - dd02 = 0.0e0_rt; - dd03 = 0.0e0_rt; - dd04 = 0.0e0_rt; - dd05 = 0.0e0_rt; - dd11 = -1.879e10_rt; - dd12 = -9.667e9_rt; - dd13 = -5.602e9_rt; - dd14 = -3.370e9_rt; - dd15 = -1.825e9_rt; - dd21 = -2.919e10_rt; - dd22 = -1.185e10_rt; - dd23 = -7.270e9_rt; - dd24 = -4.222e9_rt; - dd25 = -1.560e9_rt; - - } else if (temp >= 1.0e8_rt && temp < 1.0e9_rt) { - - tau = std::log10(temp * 1.0e-8_rt); - cc = 1.5654e0_rt; - c00 = 9.889e10_rt; - c01 = -4.524e8_rt; - c02 = -6.088e6_rt; - c03 = 4.269e7_rt; - c04 = 5.172e7_rt; - c05 = 4.910e7_rt; - c06 = 4.388e7_rt; - c10 = 1.813e11_rt; - c11 = -7.556e9_rt; - c12 = -3.304e9_rt; - c13 = -1.031e9_rt; - c14 = -1.764e9_rt; - c15 = -1.851e9_rt; - c16 = -1.928e9_rt; - c20 = 9.750e10_rt; - c21 = 3.484e10_rt; - c22 = 5.199e9_rt; - c23 = -1.695e9_rt; - c24 = -2.865e9_rt; - c25 = -3.395e9_rt; - c26 = -3.418e9_rt; - dd01 = -1.135e8_rt; - dd02 = 1.256e8_rt; - dd03 = 5.149e7_rt; - dd04 = 3.436e7_rt; - dd05 = 1.005e7_rt; - dd11 = 1.652e9_rt; - dd12 = -3.119e9_rt; - dd13 = -1.839e9_rt; - dd14 = -1.458e9_rt; - dd15 = -8.956e8_rt; - dd21 = -1.548e10_rt; - dd22 = -9.338e9_rt; - dd23 = -5.899e9_rt; - dd24 = -3.035e9_rt; - dd25 = -1.598e9_rt; - - } else { - - // T > 1.e9 - - tau = std::log10(sf.t9); - cc = 1.5654e0_rt; - c00 = 9.581e10_rt; - c01 = 4.107e8_rt; - c02 = 2.305e8_rt; - c03 = 2.236e8_rt; - c04 = 1.580e8_rt; - c05 = 2.165e8_rt; - c06 = 1.721e8_rt; - c10 = 1.459e12_rt; - c11 = 1.314e11_rt; - c12 = -1.169e11_rt; - c13 = -1.765e11_rt; - c14 = -1.867e11_rt; - c15 = -1.983e11_rt; - c16 = -1.896e11_rt; - c20 = 2.424e11_rt; - c21 = -3.669e9_rt; - c22 = -8.691e9_rt; - c23 = -7.967e9_rt; - c24 = -7.932e9_rt; - c25 = -7.987e9_rt; - c26 = -8.333e9_rt; - dd01 = 4.724e8_rt; - dd02 = 2.976e8_rt; - dd03 = 2.242e8_rt; - dd04 = 7.937e7_rt; - dd05 = 4.859e7_rt; - dd11 = -7.094e11_rt; - dd12 = -3.697e11_rt; - dd13 = -2.189e11_rt; - dd14 = -1.273e11_rt; - dd15 = -5.705e10_rt; - dd21 = -2.254e10_rt; - dd22 = -1.551e10_rt; - dd23 = -7.793e9_rt; - dd24 = -4.489e9_rt; - dd25 = -2.185e9_rt; - - } - - taudt = nu_constants::iln10*sf.tempi; - - // equation 3.7, compute the expensive trig functions only one time - cos1 = std::cos(nu_constants::fac1*tau); - cos2 = std::cos(nu_constants::fac1*2.0e0_rt*tau); - cos3 = std::cos(nu_constants::fac1*3.0e0_rt*tau); - cos4 = std::cos(nu_constants::fac1*4.0e0_rt*tau); - cos5 = std::cos(nu_constants::fac1*5.0e0_rt*tau); - last = std::cos(nu_constants::fac2*tau); - - sin1 = std::sin(nu_constants::fac1*tau); - sin2 = std::sin(nu_constants::fac1*2.0e0_rt*tau); - sin3 = std::sin(nu_constants::fac1*3.0e0_rt*tau); - sin4 = std::sin(nu_constants::fac1*4.0e0_rt*tau); - sin5 = std::sin(nu_constants::fac1*5.0e0_rt*tau); - xast = std::sin(nu_constants::fac2*tau); - - a0 = 0.5e0_rt*c00 - + 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; - 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; - - - 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; - - if constexpr (do_derivatives) { - f0 = taudt*nu_constants::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*nu_constants::fac2*taudt; - - f1 = taudt*nu_constants::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*nu_constants::fac2*taudt; - - f2 = taudt*nu_constants::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*nu_constants::fac2*taudt; - } - - - // equation 3.4 - dum = a0 + a1*zeta + a2*zeta2; - 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; - 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*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 - 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; - } - dum = 1.0e0_rt/xden; - fphot = xnum*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 * 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*sf.xldt; - } - - 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 = 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 + sf.rm*z; - if constexpr (do_derivatives) { - xdendt = -sf.rm*z*z*dumdt; - xdenda = sf.rmda*z; - xdendz = sf.rmdz*z; - } - - z = 1.0e0_rt/xden; - qphot = xnum*z; - if constexpr (do_derivatives) { - qphotdt = (xnumdt - qphot*xdendt)*z; - } - dum = -qphot*z; - if constexpr (do_derivatives) { - qphotda = dum*xdenda; - qphotdz = dum*xdendz; - } - - // equation 3.2 - sphot = sf.xl5 * fphot; - if constexpr (do_derivatives) { - sphotdt = 5.0e0_rt*sf.xl4*sf.xldt*fphot + sf.xl5*fphotdt; - sphotda = sf.xl5*fphotda; - sphotdz = sf.xl5*fphotdz; - } - - a1 = sphot; - sphot = sf.rm*a1; - if constexpr (do_derivatives) { - sphotdt = sf.rm*sphotdt; - sphotda = sf.rm*sphotda + sf.rmda*a1; - sphotdz = sf.rm*sphotdz + sf.rmdz*a1; - } - - a1 = nu_constants::tfac4*(1.0e0_rt - nu_constants::tfac3 * qphot); - - a3 = sphot; - sphot = a1*a3; - if constexpr (do_derivatives) { - a2 = -nu_constants::tfac4*nu_constants::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; - if constexpr (do_derivatives) { - sphotdt = 0.0e0_rt; - sphotda = 0.0e0_rt; - sphotdz = 0.0e0_rt; - } - } + nu_photo(sf, sphot, sphotdt, sphotda, sphotdz); nu_brem(sf, sbrem, sbremdt, sbremda, sbremdz); From 5f987dd130b3fed27c045f67dd6c9752bff8fc6e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 18:03:11 -0400 Subject: [PATCH 09/10] move plasma neutrino contribution into its own function (#1377) --- neutrinos/sneut5.H | 348 +++++++++++++++++++++++---------------------- 1 file changed, 178 insertions(+), 170 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 38a66f29bd..6a9fecab43 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -364,6 +364,180 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_plasma(const sneutf_t& sf, + Real& splas, Real& splasdt, Real& splasda, Real& splasdz) { + + // plasma neutrino section + // for collective reactions like gamma_plasmon => nu_e + nubar_e + // equation 4.6 + + Real a1 = 1.019e-6_rt * sf.rm; + Real a2 = std::pow(a1, nu_constants::twoth); + + Real b1 = std::sqrt(1.0e0_rt + a2); + + Real c00 = 1.0e0_rt / (sf.temp * sf.temp * b1); + + Real gl2 = 1.1095e11_rt * sf.rm * c00; + Real gl2dt, gl2da, gl2dz; + if constexpr (do_derivatives) { + Real a3 = nu_constants::twoth * a2 / a1; + Real b2 = 1.0e0_rt / b1; + gl2dt = -2.0e0_rt * gl2 * sf.tempi; + Real 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); + } + + Real gl = std::sqrt(gl2); + Real gl12 = std::sqrt(gl); + Real gl32 = gl * gl12; + Real gl72 = gl2 * gl32; + Real gl6 = gl2 * gl2 * gl2; + + // equation 4.7 + Real ft = 2.4e0_rt + 0.6e0_rt * gl12 + 0.51e0_rt * gl + 1.25e0_rt * gl32; + Real gum = 1.0e0_rt / gl2; + Real ftdt, ftda, ftdz; + 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; + + b1 = 225.0e0_rt - 17.0e0_rt * gl + gl2; + + Real c = 1.0e0_rt / b1; + Real fl = a1 * c; + Real fldt, flda, fldz; + if constexpr (do_derivatives) { + a2 = 8.6e0_rt + 1.75e0_rt * 1.35e0_rt * gl72 * gum; + Real b2 = -0.5e0_rt * 17.0e0_rt * gl * gum + 1.0e0_rt; + Real d = (a2 - fl * b2) * c; + fldt = d * gl2dt; + flda = d * gl2da; + fldz = d * gl2dz; + } + + // equation 4.9 and 4.10 + Real cc = std::log10(2.0e0_rt * sf.rm); + Real xlnt = std::log10(sf.temp); + + Real xnum = nu_constants::sixth * (17.5e0_rt + cc - 3.0e0_rt * xlnt); + Real xden = nu_constants::sixth * (-24.5e0_rt + cc + 3.0e0_rt * xlnt); + Real xnumdt, xnumda, xnumdz; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + xnumdt = -nu_constants::iln10 * 0.5e0_rt * sf.tempi; + a2 = nu_constants::iln10 * nu_constants::sixth * sf.rmi; + xnumda = a2 * sf.rmda; + xnumdz = a2 * sf.rmdz; + + xdendt = nu_constants::iln10 * 0.5e0_rt * sf.tempi; + xdenda = a2 * sf.rmda; + xdendz = a2 * sf.rmdz; + } + + // equation 4.11 + Real fxy, fxydt, fxyda, fxydz; + Real dumdt, dumda, dumdz; + if (std::abs(xnum) > 0.7e0_rt || xden < 0.0e0_rt) { + fxy = 1.0e0_rt; + fxydt = 0.0e0_rt; + fxydz = 0.0e0_rt; + fxyda = 0.0e0_rt; + + } else { + + a1 = 0.39e0_rt - 1.25e0_rt * xnum - 0.35e0_rt * std::sin(4.5e0_rt * xnum); + + b1 = 0.3e0_rt * std::exp(-std::pow(4.5e0_rt*xnum + 0.9e0_rt, 2.0_rt)); + + c = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt * xnum); + + Real b2; + 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; + } + } + + Real d = 0.57e0_rt - 0.25e0_rt * xnum; + Real a3 = c / d; + c00 = std::exp(-a3 * a3); + + fxy = 1.05e0_rt + (a1 - b1) * c00; + if constexpr (do_derivatives) { + Real f1 = -c00 * 2.0e0_rt * a3 / d; + Real c01 = f1 * (dumdt + a3 * 0.25e0_rt * xnumdt); + Real c03 = f1 * (dumda + a3 * 0.25e0_rt * xnumda); + Real 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; + 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); + + a1 = splas; + splas = a2*a1; + if constexpr (do_derivatives) { + Real 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; + + a1 = splas; + splas = a2 * a1; + if constexpr (do_derivatives) { + Real 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 * sf.xl9; + + a1 = splas; + splas = a2 * a1; + if constexpr (do_derivatives) { + Real 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; + } +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void nu_photo(const sneutf_t& sf, @@ -1135,11 +1309,9 @@ void sneut5(const Real temp, const Real den, dsnudz = derivative of snu with zbar */ - Real xlnt,cc, - a1,a2,a3,b1,b2,c00,c01,c03,c04, - c,d{0.0}, - f1{0.0}, - dumdt,gum,dumda,dumdz; + Real a1,a2,a3,b1,b2, + c,d{0.0}; + // pair production Real gl,gldt, @@ -1148,13 +1320,6 @@ void sneut5(const Real temp, const Real den, fpairda,fpairdz,qpairda,qpairdz, xnumda,xnumdz,xdenda,xdendz; - // plasma - Real gl2,gl2dt,gl12,gl32,gl72,gl6, - ft,ftdt,fl,fldt,fxy,fxydt, - flda,fldz,ftda,ftdz, - fxyda,fxydz,gl2da,gl2dz; - - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; @@ -1311,165 +1476,8 @@ void sneut5(const Real temp, const Real den, spairdz = a1*spairdz + a2*qpairdz*a3; } - // plasma neutrino section - // for collective reactions like gamma_plasmon => nu_e + nubar_e - // equation 4.6 - - a1 = 1.019e-6_rt*sf.rm; - a2 = std::pow(a1, nu_constants::twoth); - - b1 = std::sqrt(1.0e0_rt + a2); - - c00 = 1.0e0_rt/(temp*temp*b1); - - gl2 = 1.1095e11_rt * sf.rm * c00; - - if constexpr (do_derivatives) { - a3 = nu_constants::twoth*a2/a1; - b2 = 1.0e0_rt/b1; - 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); - gl12 = std::sqrt(gl); - gl32 = gl * gl12; - gl72 = gl2 * gl32; - gl6 = gl2 * gl2 * gl2; - - // equation 4.7 - ft = 2.4e0_rt + 0.6e0_rt*gl12 + 0.51e0_rt*gl + 1.25e0_rt*gl32; - gum = 1.0e0_rt/gl2; - 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; - - b1 = 225.0e0_rt - 17.0e0_rt*gl + gl2; - - c = 1.0e0_rt/b1; - fl = a1*c; - - 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*sf.rm); - xlnt = std::log10(temp); - - xnum = nu_constants::sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt); - xden = nu_constants::sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt); - if constexpr (do_derivatives) { - xnumdt = -nu_constants::iln10*0.5e0_rt*sf.tempi; - a2 = nu_constants::iln10*nu_constants::sixth*sf.rmi; - xnumda = a2*sf.rmda; - xnumdz = a2*sf.rmdz; - - xdendt = nu_constants::iln10*0.5e0_rt*sf.tempi; - xdenda = a2*sf.rmda; - xdendz = a2*sf.rmdz; - } - - // equation 4.11 - if (std::abs(xnum) > 0.7e0_rt || xden < 0.0e0_rt) { - fxy = 1.0e0_rt; - fxydt = 0.0e0_rt; - fxydz = 0.0e0_rt; - fxyda = 0.0e0_rt; - - } else { - - a1 = 0.39e0_rt - 1.25e0_rt*xnum - 0.35e0_rt*std::sin(4.5e0_rt*xnum); - - b1 = 0.3e0_rt * std::exp(-std::pow(4.5e0_rt*xnum + 0.9e0_rt, 2.0_rt)); - - c = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt*xnum); - - 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); - - fxy = 1.05e0_rt + (a1 - b1)*c00; - 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; - 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); - - a1 = splas; - splas = a2*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; - - a1 = splas; - splas = a2*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 * sf.xl9; - - a1 = splas; - splas = a2*a1; - if constexpr (do_derivatives) { - 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; - } + nu_plasma(sf, splas, splasdt, splasda, splasdz); nu_photo(sf, sphot, sphotdt, sphotda, sphotdz); From 04ba99af58522437f7d0f7ad33b551486653bebb Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 29 Oct 2023 19:20:59 -0400 Subject: [PATCH 10/10] move pair neutrinos into it's own function in sneut5 (#1378) --- neutrinos/sneut5.H | 267 +++++++++++++++++++++++---------------------- 1 file changed, 136 insertions(+), 131 deletions(-) diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 6a9fecab43..4b26c89469 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -365,6 +365,141 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) { } + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void nu_pair(const sneutf_t& sf, + Real& spair, Real& spairdt, Real& spairda, Real& spairdz) { + + // pair neutrino section + // for reactions like e+ + e- => nu_e + nubar_e + + // equation 2.8 + Real gl = 1.0e0_rt - 13.04e0_rt * sf.xl2 +133.5e0_rt * sf.xl4 +1534.0e0_rt * sf.xl6 + 918.6e0_rt * sf.xl8; + Real gldt; + if constexpr (do_derivatives) { + 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 + + Real a1 = 6.002e19_rt + 2.084e20_rt * sf.zeta + 1.872e21_rt * sf.zeta2; + Real b1, b2; + if (sf.t9 < 10.0_rt) { + b1 = std::exp(-5.5924e0_rt * sf.zeta); + if constexpr (do_derivatives) { + b2 = -b1 * 5.5924e0_rt; + } + } else { + b1 = std::exp(-4.9924e0_rt * sf.zeta); + if constexpr (do_derivatives) { + b2 = -b1 * 4.9924e0_rt; + } + } + + Real a2, c; + + Real xnum = a1 * b1; + Real xnumdt, xnumda, xnumdz; + if constexpr (do_derivatives) { + a2 = 2.084e20_rt + 2.0e0_rt * 1.872e21_rt * sf.zeta; + c = a2 * b1 + a1 * b2; + xnumdt = c * sf.zetadt; + xnumda = c * sf.zetada; + xnumdz = c * sf.zetadz; + } + + 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 * 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 * sf.xlm1 - 8.141e-1_rt * sf.xlm2; + if constexpr (do_derivatives) { + a2 = -1.2383e0_rt * sf.xlm2 + 2.0e0_rt * 8.141e-1_rt * sf.xlm3; + } + } + + Real xden = sf.zeta3 + a1; + Real xdendt, xdenda, xdendz; + if constexpr (do_derivatives) { + b1 = 3.0e0_rt * sf.zeta2; + xdendt = b1 * sf.zetadt + a2 * sf.xldt; + xdenda = b1 * sf.zetada; + xdendz = b1 * sf.zetadz; + } + + a1 = 1.0e0_rt / xden; + Real fpair = xnum * a1; + Real fpairdt, fpairda, fpairdz; + 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 * sf.xl2 + 0.3967e0_rt * sf.xlp5 + 1.005e0_rt; + xnum = 1.0e0_rt / a1; + if constexpr (do_derivatives) { + 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 * sf.xl3 + 9.715e6_rt * sf.xlp5; + c = 1.0e0_rt / a1; + b1 = 1.0e0_rt + sf.rm * c; + + xden = std::pow(b1, -0.3e0_rt); + if constexpr (do_derivatives) { + a2 = sf.xldt * (3.0e0_rt * 7.692e7_rt * sf.xl2 + 0.5e0_rt * 9.715e6_rt * sf.xlmp5); + Real d = -0.3e0_rt * xden / b1; + xdendt = -d * sf.rm * c * c * a2; + xdenda = d * sf.rmda * c; + xdendz = d * sf.rmdz * c; + } + + Real qpair = xnum * xden; + Real qpairdt, qpairda, qpairdz; + if constexpr (do_derivatives) { + qpairdt = xnumdt * xden + xnum * xdendt; + qpairda = xnum * xdenda; + qpairdz = xnum * xdendz; + } + + // equation 2.5 + a1 = std::exp(-2.0e0_rt * sf.xlm1); + + spair = a1 * fpair; + if constexpr (do_derivatives) { + a2 = a1 * 2.0e0_rt * sf.xlm2 * sf.xldt; + spairdt = a2 * fpair + a1 * fpairdt; + spairda = a1 * fpairda; + spairdz = a1 * fpairdz; + } + + a1 = spair; + spair = gl*a1; + if constexpr (do_derivatives) { + spairdt = gl * spairdt + gldt * a1; + spairda = gl * spairda; + spairdz = gl * spairdz; + } + + a1 = nu_constants::tfac4 * (1.0e0_rt + nu_constants::tfac3 * qpair); + + Real a3 = spair; + spair = a1 * a3; + if constexpr (do_derivatives) { + a2 = nu_constants::tfac4 * nu_constants::tfac3; + spairdt = a1 * spairdt + a2 * qpairdt * a3; + spairda = a1 * spairda + a2 * qpairda * a3; + spairdz = a1 * spairdz + a2 * qpairdz * a3; + } + +} + template AMREX_GPU_HOST_DEVICE AMREX_INLINE void nu_plasma(const sneutf_t& sf, @@ -1309,17 +1444,6 @@ void sneut5(const Real temp, const Real den, dsnudz = derivative of snu with zbar */ - Real a1,a2,a3,b1,b2, - c,d{0.0}; - - - // pair production - Real gl,gldt, - xnum,xnumdt,xden,xdendt, - fpair,fpairdt,qpair,qpairdt, - fpairda,fpairdz,qpairda,qpairdz, - xnumda,xnumdz,xdenda,xdendz; - // initialize Real spair{0.0e0_rt}; Real spairdt{0.0e0_rt}; @@ -1356,126 +1480,7 @@ void sneut5(const Real temp, const Real den, auto sf = get_sneut_factors(den, temp, abar, zbar); - // pair neutrino section - // for reactions like e+ + e- => nu_e + nubar_e - - // equation 2.8 - 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 = 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 * sf.zeta + 1.872e21_rt * sf.zeta2; - - if (sf.t9 < 10.0_rt) { - b1 = std::exp(-5.5924e0_rt * sf.zeta); - if constexpr (do_derivatives) { - b2 = -b1*5.5924e0_rt; - } - } else { - b1 = std::exp(-4.9924e0_rt * sf.zeta); - if constexpr (do_derivatives) { - b2 = -b1*4.9924e0_rt; - } - } - - xnum = a1 * b1; - if constexpr (do_derivatives) { - a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt * sf.zeta; - c = a2*b1 + a1*b2; - xnumdt = c * sf.zetadt; - xnumda = c * sf.zetada; - xnumdz = c * sf.zetadz; - } - - 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*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*sf.xlm1 - 8.141e-1_rt*sf.xlm2; - if constexpr (do_derivatives) { - a2 = -1.2383e0_rt*sf.xlm2 + 2.0e0_rt*8.141e-1_rt*sf.xlm3; - } - } - - xden = sf.zeta3 + a1; - if constexpr (do_derivatives) { - b1 = 3.0e0_rt*sf.zeta2; - xdendt = b1*sf.zetadt + a2*sf.xldt; - xdenda = b1*sf.zetada; - xdendz = b1*sf.zetadz; - } - - a1 = 1.0e0_rt/xden; - fpair = xnum*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*sf.xl2 + 0.3967e0_rt*sf.xlp5 + 1.005e0_rt; - xnum = 1.0e0_rt/a1; - if constexpr (do_derivatives) { - 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*sf.xl3 + 9.715e6_rt*sf.xlp5; - c = 1.0e0_rt/a1; - b1 = 1.0e0_rt + sf.rm*c; - - xden = std::pow(b1, -0.3e0_rt); - if constexpr (do_derivatives) { - 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*sf.rm*c*c*a2; - xdenda = d*sf.rmda*c; - xdendz = d*sf.rmdz*c; - } - - qpair = xnum*xden; - if constexpr (do_derivatives) { - qpairdt = xnumdt*xden + xnum*xdendt; - qpairda = xnum*xdenda; - qpairdz = xnum*xdendz; - } - - // equation 2.5 - a1 = std::exp(-2.0e0_rt*sf.xlm1); - - spair = a1*fpair; - if constexpr (do_derivatives) { - a2 = a1*2.0e0_rt*sf.xlm2*sf.xldt; - spairdt = a2*fpair + a1*fpairdt; - spairda = a1*fpairda; - spairdz = a1*fpairdz; - } - - a1 = spair; - spair = gl*a1; - if constexpr (do_derivatives) { - spairdt = gl*spairdt + gldt*a1; - spairda = gl*spairda; - spairdz = gl*spairdz; - } - - a1 = nu_constants::tfac4*(1.0e0_rt + nu_constants::tfac3 * qpair); - - a3 = spair; - spair = a1*a3; - if constexpr (do_derivatives) { - a2 = nu_constants::tfac4 * nu_constants::tfac3; - spairdt = a1*spairdt + a2*qpairdt*a3; - spairda = a1*spairda + a2*qpairda*a3; - spairdz = a1*spairdz + a2*qpairdz*a3; - } - + nu_pair(sf, spair, spairdt, spairda, spairdz); nu_plasma(sf, splas, splasdt, splasda, splasdz);