diff --git a/.github/workflows/test_neutrinos.yml b/.github/workflows/test_neutrinos.yml new file mode 100644 index 0000000000..53506e44ac --- /dev/null +++ b/.github/workflows/test_neutrinos.yml @@ -0,0 +1,49 @@ +name: test_neutrinos + +on: [pull_request] +jobs: + test_neutrinos: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + with: + fetch-depth: 0 + + - name: Get AMReX + run: | + mkdir external + cd external + git clone https://github.com/AMReX-Codes/amrex.git + cd amrex + git checkout development + echo 'AMREX_HOME=$(GITHUB_WORKSPACE)/external/amrex' >> $GITHUB_ENV + echo $AMREX_HOME + if [[ -n "${AMREX_HOME}" ]]; then exit 1; fi + cd ../.. + + - name: Install dependencies + run: | + sudo apt-get update -y -qq + sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0 + + - name: Build the fextrema tool + run: | + cd external/amrex/Tools/Plotfile + make programs=fextrema -j 2 + + - name: Compile + run: | + cd unit_test/test_neutrino_cooling + make realclean + make -j 4 + + - name: Run test_neutrino_cooling + run: | + cd unit_test/test_neutrino_cooling + ./main3d.gnu.ex inputs + ../../external/amrex/Tools/Plotfile/fextrema.gnu.ex test_sneut5 > test.out + + - name: Compare to stored output + run: | + cd unit_test/test_neutrino_cooling + diff test.out ci-benchmarks/test_sneut5.out diff --git a/CHANGES.md b/CHANGES.md index 0c7dfe56c8..82f4cea0b0 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -3,6 +3,11 @@ * The file NETWORK_PROPERTIES has been removed from each network, as the legacy screening method is no longer used. (#1310) + * The rprox network was updated and the Jacobian was fixed (#1300) + + * The primordial_chem EOS now can take density and pressure as + inputs (#1302) + # 23.07 * The preprocessor variable EXTRA_THERMO has been removed. @@ -17,7 +22,7 @@ for estimating the spectral radius of our ODE system (#1222) * removed SDC_EVOLVE_ENTHALPY -- this was not being used (#1204) - + # 23.06 * Added a new Runge-Kutta-Chebyshev integrator (#1191) diff --git a/CMakeLists.txt b/CMakeLists.txt index 49bbb2988c..c3e2702b96 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -63,7 +63,7 @@ function(setup_target_for_microphysics_compilation network_name output_dir) ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/primordial_chem ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/primordial_chem ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/constants - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null PARENT_SCOPE) + PARENT_SCOPE) #we need to have extern_parameters.cpp be available at configure time #the script write_probin.py writes this .cpp file so we call it here diff --git a/EOS/breakout/actual_eos.H b/EOS/breakout/actual_eos.H index e55123f6a0..342b36fba7 100644 --- a/EOS/breakout/actual_eos.H +++ b/EOS/breakout/actual_eos.H @@ -123,7 +123,7 @@ void actual_eos (I input, T& state) state.gam1 = gamma_const; // sound speed - state.cs = sqrt(gamma_const * poverrho); + state.cs = std::sqrt(gamma_const * poverrho); state.dpdr_e = poverrho; state.dpde = (gamma_const - 1.0_rt) * state.rho; diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H index 74700e64dd..1ef7779926 100644 --- a/integration/VODE/vode_dvhin.H +++ b/integration/VODE/vode_dvhin.H @@ -44,24 +44,6 @@ void dvhin (BurnT& state, DvodeT& vstate, Real& H0, int& NITER, int& IER) // Set an upper bound on h based on vstate.tout-vstate.t and the initial Y and YDOT. Real HUB = PT1 * TDIST; -#ifdef TRUE_SDC - for (int i = 1; i <= int_neqs; ++i) { - Real atol; - if (i == 1) { - atol = vstate.atol_dens; - } else if (i == NumSpec+2) { - atol = vstate.atol_enuc; - } else { - atol = vstate.atol_spec; - } - Real DELYI = PT1 * std::abs(vstate.yh(i,1)) + atol; - Real AFI = std::abs(vstate.yh(i,2)); - if (AFI * HUB > DELYI) { - HUB = DELYI / AFI; - } - } - -#else for (int i = 1; i <= int_neqs; ++i) { Real atol = i <= NumSpec ? vstate.atol_spec : vstate.atol_enuc; Real DELYI = PT1 * std::abs(vstate.yh(i,1)) + atol; @@ -70,7 +52,6 @@ void dvhin (BurnT& state, DvodeT& vstate, Real& H0, int& NITER, int& IER) HUB = DELYI / AFI; } } -#endif // Set initial guess for h as geometric mean of upper and lower bounds. int iter = 0; diff --git a/integration/VODE/vode_dvode.H b/integration/VODE/vode_dvode.H index 9466d2676b..924c010325 100644 --- a/integration/VODE/vode_dvode.H +++ b/integration/VODE/vode_dvode.H @@ -72,24 +72,12 @@ int dvode (BurnT& state, DvodeT& vstate) vstate.NQ = 1; vstate.H = 1.0_rt; -#ifdef TRUE_SDC - vstate.ewt(1) = vstate.rtol_dens * std::abs(vstate.yh(1,1)) + vstate.atol_dens; - vstate.ewt(1) = 1.0_rt / vstate.ewt(1); - - for (int i = 2; i <= NumSpec+1; ++i) { - vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec; - vstate.ewt(i) = 1.0_rt / vstate.ewt(i); - } - vstate.ewt(NumSpec+2) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+2,1)) + vstate.atol_enuc; - vstate.ewt(NumSpec+2) = 1.0_rt / vstate.ewt(NumSpec+2); -#else for (int i = 1; i <= NumSpec; ++i) { vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec; vstate.ewt(i) = 1.0_rt / vstate.ewt(i); } vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc; vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1); -#endif // Call DVHIN to set initial step size H0 to be attempted. H0 = 0.0_rt; @@ -157,24 +145,12 @@ int dvode (BurnT& state, DvodeT& vstate) } -#ifdef TRUE_SDC - vstate.ewt(1) = vstate.rtol_dens * std::abs(vstate.yh(1,1)) + vstate.atol_dens; - vstate.ewt(1) = 1.0_rt / vstate.ewt(1); - - for (int i = 2; i <= NumSpec+1; ++i) { - vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec; - vstate.ewt(i) = 1.0_rt / vstate.ewt(i); - } - vstate.ewt(NumSpec+2) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+2,1)) + vstate.atol_enuc; - vstate.ewt(NumSpec+2) = 1.0_rt / vstate.ewt(NumSpec+2); -#else for (int i = 1; i <= NumSpec; ++i) { vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec; vstate.ewt(i) = 1.0_rt / vstate.ewt(i); } vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc; vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1); -#endif } else { diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index db28b936f1..5b850225eb 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -149,7 +149,7 @@ int dvstep (BurnT& state, DvodeT& vstate) y_save(i) = vstate.y(i); } #endif -#ifdef SIMPLIFIED_SDC +#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC) Real rho_old = state.rho_orig + TOLD * state.ydot_a[SRHO]; for (int i = 1; i <= int_neqs; ++i) { y_save(i) = vstate.y(i); @@ -224,7 +224,7 @@ int dvstep (BurnT& state, DvodeT& vstate) bool valid_update = true; -#ifdef SIMPLIFIED_SDC +#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC) if (vstate.y(SEINT+1) < 0.0_rt) { valid_update = false; } @@ -271,7 +271,7 @@ int dvstep (BurnT& state, DvodeT& vstate) #endif -#ifdef SIMPLIFIED_SDC +#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC) // these are basically the same checks as with Strang // above, except now we are evolving rhoX instead of X. @@ -303,24 +303,6 @@ int dvstep (BurnT& state, DvodeT& vstate) #endif -#ifdef TRUE_SDC - // as above, we want to make sure that the mass fractions stay bounded, but we are - // evolving: - // - // y(1) = density - // y(2:2+NumSpec-1) = rho X - // y(NumSpec+2) = rho e - - if (vstate.y(1+i) < -species_failure_tolerance * vstate.y(1)) { - valid_update = false; - } - - if (vstate.y(1+i) > (1.0_rt + species_failure_tolerance) * vstate.y(1)) { - valid_update = false; - } -#endif - - } // The corrector has converged (NFLAG = 0). The local error test is diff --git a/integration/integrator_data.H b/integration/integrator_data.H index 62330ff03e..5f64ca7e07 100644 --- a/integration/integrator_data.H +++ b/integration/integrator_data.H @@ -14,7 +14,7 @@ const int INT_NEQS = SVAR_EVOLVE; #endif #ifdef TRUE_SDC -const int INT_NEQS = NumSpec + 2; +const int INT_NEQS = NumSpec + 1; #endif @@ -56,7 +56,7 @@ constexpr int integrator_neqs () #endif #ifdef TRUE_SDC - int_neqs = NumSpec + 2; + int_neqs = NumSpec + 1; #endif return int_neqs; diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H index 77d14f3231..9e7edbd4df 100644 --- a/interfaces/burn_type.H +++ b/interfaces/burn_type.H @@ -131,12 +131,6 @@ struct burn_t // dx is useful for estimating timescales for equilibriation Real dx; -#ifdef TRUE_SDC - Real E_var; - Real mom[3]; - Real f_source[NumSpec+2]; -#endif - Real mu; Real mu_e; Real y_e; diff --git a/networks/CNO_extras/actual_rhs.H b/networks/CNO_extras/actual_rhs.H index 79e1367582..03501b4031 100644 --- a/networks/CNO_extras/actual_rhs.H +++ b/networks/CNO_extras/actual_rhs.H @@ -794,8 +794,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1197,7 +1197,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H index b7aa913121..0fb4bc3422 100644 --- a/networks/ECSN/actual_rhs.H +++ b/networks/ECSN/actual_rhs.H @@ -412,8 +412,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); // Append the energy equation (this is erg/g/s) @@ -629,7 +629,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz); diff --git a/networks/ase/actual_rhs.H b/networks/ase/actual_rhs.H index 26445601af..242edac898 100644 --- a/networks/ase/actual_rhs.H +++ b/networks/ase/actual_rhs.H @@ -1093,8 +1093,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1646,7 +1646,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/general_null/network_header.template b/networks/general_null/network_header.template index 0a7c3ecec1..949863298a 100644 --- a/networks/general_null/network_header.template +++ b/networks/general_null/network_header.template @@ -23,9 +23,7 @@ constexpr int NumSpec = @@NSPEC@@; -#define NUM_EXTRA_SPECIES @@NEXTRASPEC@@ - -constexpr int NumSpecExtra = NUM_EXTRA_SPECIES; +constexpr int NumSpecExtra = @@NEXTRASPEC@@; constexpr int NumSpecTotal = NumSpec + NumSpecExtra; diff --git a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H index 6a84156b9a..5422a92f35 100644 --- a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H +++ b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H @@ -212,8 +212,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -309,7 +309,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H index cc76ef7078..be1b24da39 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_rhs.H +++ b/networks/ignition_reaclib/URCA-simple/actual_rhs.H @@ -234,8 +234,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -343,7 +343,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/nova/actual_rhs.H b/networks/nova/actual_rhs.H index a9e79d1fcc..260a47e241 100644 --- a/networks/nova/actual_rhs.H +++ b/networks/nova/actual_rhs.H @@ -487,8 +487,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -746,7 +746,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/nova2/actual_rhs.H b/networks/nova2/actual_rhs.H index 49ea50923f..9cae820c4b 100644 --- a/networks/nova2/actual_rhs.H +++ b/networks/nova2/actual_rhs.H @@ -727,8 +727,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1058,7 +1058,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/partition_test/actual_rhs.H b/networks/partition_test/actual_rhs.H index 37d0cac6c0..cbff936ed8 100644 --- a/networks/partition_test/actual_rhs.H +++ b/networks/partition_test/actual_rhs.H @@ -202,8 +202,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -335,7 +335,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/rhs.H b/networks/rhs.H index 55a0060be2..e507d786aa 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -1348,9 +1348,10 @@ void rhs (burn_t& burn_state, Array1D& ydot) // Evaluate the neutrino cooling. #ifdef NEUTRINOS + constexpr int do_derivatives{0}; Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, - sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, + sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); #else Real sneut = 0.0; #endif @@ -1495,8 +1496,9 @@ void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) // Evaluate the neutrino cooling. #ifdef NEUTRINOS Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, - sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); + constexpr int do_derivatives{1}; + sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, + sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); #else Real sneut = 0.0, dsneutdt = 0.0, dsneutdd = 0.0, dsnuda = 0.0, dsnudz = 0.0; amrex::ignore_unused(sneut, dsneutdd); diff --git a/networks/sn160/actual_rhs.H b/networks/sn160/actual_rhs.H index 8f81f200bb..647a980a70 100644 --- a/networks/sn160/actual_rhs.H +++ b/networks/sn160/actual_rhs.H @@ -10396,8 +10396,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -18028,7 +18028,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/subch_approx/actual_rhs.H b/networks/subch_approx/actual_rhs.H index a9b90e7e67..c4e168ec2c 100644 --- a/networks/subch_approx/actual_rhs.H +++ b/networks/subch_approx/actual_rhs.H @@ -1375,8 +1375,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1928,7 +1928,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/subch_full/actual_rhs.H b/networks/subch_full/actual_rhs.H index d082bd48de..c27ce0b3d7 100644 --- a/networks/subch_full/actual_rhs.H +++ b/networks/subch_full/actual_rhs.H @@ -1603,8 +1603,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -2390,7 +2390,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/networks/subch_simple/actual_rhs.H b/networks/subch_simple/actual_rhs.H index 68e84fa45f..eae7f7bbbb 100644 --- a/networks/subch_simple/actual_rhs.H +++ b/networks/subch_simple/actual_rhs.H @@ -1227,8 +1227,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); // Append the energy equation (this is erg/g/s) @@ -1723,7 +1723,8 @@ void actual_jac(const burn_t& state, MatrixType& jac) // Account for the thermal neutrino losses Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{1}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); for (int j = 1; j <= NumSpec; ++j) { Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H index 6b066a2895..9348f75f86 100644 --- a/neutrinos/sneut5.H +++ b/neutrinos/sneut5.H @@ -25,7 +25,8 @@ Real ifermi12(const Real f) // the return value Real ifermi12r; - // load the coefficients of the expansion + // load the coefficients of the expansion from Table 8 of Antia + an = 0.5e0_rt; m1 = 4; k1 = 3; @@ -60,36 +61,55 @@ Real ifermi12(const Real f) if (f < 4.0e0_rt) { - rn = f + a1(m1); + // build sum_{i=0, m1} a_i x**i. This is the numerator + // in Eq. 4 of Antia. + // + // with our 1-based indexing, this expression is + // a[1] + a[2] * f + a[3] * f**2 + ... a[m1+1] * f**m1 + // + // we do the sum starting with the largest term and working + // on a single power of f each iteration. + // + // in the starting rn here, the leading f is actually + // a1(m1+1) * f, but that element is 1 + rn = f + a1(m1); - for (int i = m1 - 1; i >= 1; --i) { - rn = rn*f + a1(i); - } + for (int i = m1 - 1; i >= 1; --i) { + rn = rn*f + a1(i); + } - den = b1(k1+1); + // now we do the denominator in Eq. 4. None of the coefficients + // are 1, so we loop over all + den = b1(k1+1); - for (int i = k1; i >= 1; --i) { - den = den*f + b1(i); - } + for (int i = k1; i >= 1; --i) { + den = den*f + b1(i); + } + + // Eq. 6 of Antia - ifermi12r = std::log(f * rn/den); + ifermi12r = std::log(f * rn/den); } else { - ff = 1.0e0_rt/std::pow(f, 1.0e0_rt/(1.0e0_rt + an)); - rn = ff + a2(m2); + ff = 1.0e0_rt/std::pow(f, 1.0e0_rt/(1.0e0_rt + an)); - for (int i = m2 - 1; i >= 1; --i) { - rn = rn*ff + a2(i); - } + // this construction is the same as above, but using the + // second set of coefficients - den = b2(k2+1); + rn = ff + a2(m2); - for (int i = k2; i >= 1; --i) { - den = den*ff + b2(i); - } + for (int i = m2 - 1; i >= 1; --i) { + rn = rn*ff + a2(i); + } + + den = b2(k2+1); + + for (int i = k2; i >= 1; --i) { + den = den*ff + b2(i); + } - ifermi12r = rn/(den*ff); + ifermi12r = rn/(den*ff); } @@ -113,7 +133,8 @@ Real zfermim12(const Real x) // return value Real zfermim12r; - // load the coefficients of the expansion + // load the coefficients of the expansion from Table 2 of Antia + m1 = 7; k1 = 7; m2 = 11; @@ -202,6 +223,8 @@ Real zfermim12(const Real x) return zfermim12r; } + +template AMREX_GPU_HOST_DEVICE AMREX_INLINE void sneut5(const Real temp, const Real den, const Real abar, const Real zbar, @@ -231,8 +254,8 @@ void sneut5(const Real temp, const Real den, a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06, c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, - dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d,f0, - f1,deni,tempi,abari,zbari,f2,f3,z,ye, + dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d{0.0},f0, + f1{0.0},deni,tempi,abari,zbari,f2,f3,z,ye, dum,dumdt,gum,gumdt,dumda,dumdz,gumda,gumdz; // pair production @@ -347,7 +370,7 @@ void sneut5(const Real temp, const Real den, t9 = temp * 1.0e-9_rt; xl = t9 * con1; xldt = 1.0e-9_rt * con1; - xlp5 = sqrt(xl); + xlp5 = std::sqrt(xl); xl2 = xl*xl; xl3 = xl2*xl; xl4 = xl3*xl; @@ -366,15 +389,18 @@ void sneut5(const Real temp, const Real den, rmda = -rm*abari; rmdz = den*abari; rmi = 1.0e0_rt/rm; -; + a0 = rm * 1.0e-9_rt; a1 = std::pow(a0, oneth); zeta = a1 * xlm1; - zetadt = -a1 * xlm2 * xldt; - a2 = oneth * a1*rmi * xlm1; - zetada = a2 * rmda; - zetadz = a2 * rmdz; -; + + if constexpr (do_derivatives) { + a2 = oneth * a1*rmi * xlm1; + zetadt = -a1 * xlm2 * xldt; + zetada = a2 * rmda; + zetadz = a2 * rmdz; + } + zeta2 = zeta * zeta; zeta3 = zeta2 * zeta; @@ -383,95 +409,120 @@ void sneut5(const Real temp, const Real den, // equation 2.8 gl = 1.0e0_rt - 13.04e0_rt*xl2 +133.5e0_rt*xl4 +1534.0e0_rt*xl6 +918.6e0_rt*xl8; - gldt = xldt*(-26.08e0_rt*xl +534.0e0_rt*xl3 +9204.0e0_rt*xl5 +7348.8e0_rt*xl7); + if constexpr (do_derivatives) { + gldt = xldt*(-26.08e0_rt*xl +534.0e0_rt*xl3 +9204.0e0_rt*xl5 +7348.8e0_rt*xl7); + } // equation 2.7 a1 = 6.002e19_rt + 2.084e20_rt*zeta + 1.872e21_rt*zeta2; - a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta; if (t9 < 10.0_rt) { - b1 = exp(-5.5924e0_rt*zeta); - b2 = -b1*5.5924e0_rt; + b1 = std::exp(-5.5924e0_rt*zeta); + if constexpr (do_derivatives) { + b2 = -b1*5.5924e0_rt; + } } else { - b1 = exp(-4.9924e0_rt*zeta); - b2 = -b1*4.9924e0_rt; + b1 = std::exp(-4.9924e0_rt*zeta); + if constexpr (do_derivatives) { + b2 = -b1*4.9924e0_rt; + } } xnum = a1 * b1; - c = a2*b1 + a1*b2; - xnumdt = c*zetadt; - xnumda = c*zetada; - xnumdz = c*zetadz; + if constexpr (do_derivatives) { + a2 = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta; + c = a2*b1 + a1*b2; + xnumdt = c*zetadt; + xnumda = c*zetada; + xnumdz = c*zetadz; + } if (t9 < 10.0_rt) { a1 = 9.383e-1_rt*xlm1 - 4.141e-1_rt*xlm2 + 5.829e-2_rt*xlm3; - a2 = -9.383e-1_rt*xlm2 + 2.0e0_rt*4.141e-1_rt*xlm3 - 3.0e0_rt*5.829e-2_rt*xlm4; + if constexpr (do_derivatives) { + a2 = -9.383e-1_rt*xlm2 + 2.0e0_rt*4.141e-1_rt*xlm3 - 3.0e0_rt*5.829e-2_rt*xlm4; + } } else { a1 = 1.2383e0_rt*xlm1 - 8.141e-1_rt*xlm2; - a2 = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3; + if constexpr (do_derivatives) { + a2 = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3; + } } - b1 = 3.0e0_rt*zeta2; - xden = zeta3 + a1; - xdendt = b1*zetadt + a2*xldt; - xdenda = b1*zetada; - xdendz = b1*zetadz; + if constexpr (do_derivatives) { + b1 = 3.0e0_rt*zeta2; + xdendt = b1*zetadt + a2*xldt; + xdenda = b1*zetada; + xdendz = b1*zetadz; + } a1 = 1.0e0_rt/xden; fpair = xnum*a1; - fpairdt = (xnumdt - fpair*xdendt)*a1; - fpairda = (xnumda - fpair*xdenda)*a1; - fpairdz = (xnumdz - fpair*xdendz)*a1; + if constexpr (do_derivatives) { + fpairdt = (xnumdt - fpair*xdendt)*a1; + fpairda = (xnumda - fpair*xdenda)*a1; + fpairdz = (xnumdz - fpair*xdendz)*a1; + } // equation 2.6 a1 = 10.7480e0_rt*xl2 + 0.3967e0_rt*xlp5 + 1.005e0_rt; - a2 = xldt*(2.0e0_rt*10.7480e0_rt*xl + 0.5e0_rt*0.3967e0_rt*xlmp5); xnum = 1.0e0_rt/a1; - xnumdt = -xnum*xnum*a2; + if constexpr (do_derivatives) { + a2 = xldt*(2.0e0_rt*10.7480e0_rt*xl + 0.5e0_rt*0.3967e0_rt*xlmp5); + xnumdt = -xnum*xnum*a2; + } a1 = 7.692e7_rt*xl3 + 9.715e6_rt*xlp5; - a2 = xldt*(3.0e0_rt*7.692e7_rt*xl2 + 0.5e0_rt*9.715e6_rt*xlmp5); - c = 1.0e0_rt/a1; b1 = 1.0e0_rt + rm*c; xden = std::pow(b1, -0.3e0_rt); - - d = -0.3e0_rt*xden/b1; - xdendt = -d*rm*c*c*a2; - xdenda = d*rmda*c; - xdendz = d*rmdz*c; + if constexpr (do_derivatives) { + a2 = xldt*(3.0e0_rt*7.692e7_rt*xl2 + 0.5e0_rt*9.715e6_rt*xlmp5); + d = -0.3e0_rt*xden/b1; + xdendt = -d*rm*c*c*a2; + xdenda = d*rmda*c; + xdendz = d*rmdz*c; + } qpair = xnum*xden; - qpairdt = xnumdt*xden + xnum*xdendt; - qpairda = xnum*xdenda; - qpairdz = xnum*xdendz; + if constexpr (do_derivatives) { + qpairdt = xnumdt*xden + xnum*xdendt; + qpairda = xnum*xdenda; + qpairdz = xnum*xdendz; + } // equation 2.5 a1 = std::exp(-2.0e0_rt*xlm1); - a2 = a1*2.0e0_rt*xlm2*xldt; spair = a1*fpair; - spairdt = a2*fpair + a1*fpairdt; - spairda = a1*fpairda; - spairdz = a1*fpairdz; + if constexpr (do_derivatives) { + a2 = a1*2.0e0_rt*xlm2*xldt; + spairdt = a2*fpair + a1*fpairdt; + spairda = a1*fpairda; + spairdz = a1*fpairdz; + } a1 = spair; spair = gl*a1; - spairdt = gl*spairdt + gldt*a1; - spairda = gl*spairda; - spairdz = gl*spairdz; + if constexpr (do_derivatives) { + spairdt = gl*spairdt + gldt*a1; + spairda = gl*spairda; + spairdz = gl*spairdz; + } a1 = tfac4*(1.0e0_rt + tfac3 * qpair); - a2 = tfac4*tfac3; a3 = spair; spair = a1*a3; - spairdt = a1*spairdt + a2*qpairdt*a3; - spairda = a1*spairda + a2*qpairda*a3; - spairdz = a1*spairdz + a2*qpairdz*a3; + if constexpr (do_derivatives) { + a2 = tfac4*tfac3; + spairdt = a1*spairdt + a2*qpairdt*a3; + spairda = a1*spairda + a2*qpairda*a3; + spairdz = a1*spairdz + a2*qpairdz*a3; + } // plasma neutrino section // for collective reactions like gamma_plasmon => nu_e + nubar_e @@ -479,19 +530,21 @@ void sneut5(const Real temp, const Real den, a1 = 1.019e-6_rt*rm; a2 = std::pow(a1, twoth); - a3 = twoth*a2/a1; b1 = std::sqrt(1.0e0_rt + a2); - b2 = 1.0e0_rt/b1; c00 = 1.0e0_rt/(temp*temp*b1); gl2 = 1.1095e11_rt * rm * c00; - gl2dt = -2.0e0_rt*gl2*tempi; - d = rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt; - gl2da = 1.1095e11_rt * (rmda*c00 - d*rmda); - gl2dz = 1.1095e11_rt * (rmdz*c00 - d*rmdz); + if constexpr (do_derivatives) { + a3 = twoth*a2/a1; + b2 = 1.0e0_rt/b1; + gl2dt = -2.0e0_rt*gl2*tempi; + d = rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt; + gl2da = 1.1095e11_rt * (rmda*c00 - d*rmda); + gl2dz = 1.1095e11_rt * (rmdz*c00 - d*rmdz); + } gl = std::sqrt(gl2); gl12 = std::sqrt(gl); @@ -502,40 +555,46 @@ void sneut5(const Real temp, const Real den, // equation 4.7 ft = 2.4e0_rt + 0.6e0_rt*gl12 + 0.51e0_rt*gl + 1.25e0_rt*gl32; gum = 1.0e0_rt/gl2; - a1 =(0.25e0_rt*0.6e0_rt*gl12 +0.5e0_rt*0.51e0_rt*gl +0.75e0_rt*1.25e0_rt*gl32)*gum; - ftdt = a1*gl2dt; - ftda = a1*gl2da; - ftdz = a1*gl2dz; + if constexpr (do_derivatives) { + a1 =(0.25e0_rt*0.6e0_rt*gl12 +0.5e0_rt*0.51e0_rt*gl +0.75e0_rt*1.25e0_rt*gl32)*gum; + ftdt = a1*gl2dt; + ftda = a1*gl2da; + ftdz = a1*gl2dz; + } // equation 4.8 a1 = 8.6e0_rt*gl2 + 1.35e0_rt*gl72; - a2 = 8.6e0_rt + 1.75e0_rt*1.35e0_rt*gl72*gum; b1 = 225.0e0_rt - 17.0e0_rt*gl + gl2; - b2 = -0.5e0_rt*17.0e0_rt*gl*gum + 1.0e0_rt; c = 1.0e0_rt/b1; fl = a1*c; - d = (a2 - fl*b2)*c; - fldt = d*gl2dt; - flda = d*gl2da; - fldz = d*gl2dz; + if constexpr (do_derivatives) { + a2 = 8.6e0_rt + 1.75e0_rt*1.35e0_rt*gl72*gum; + b2 = -0.5e0_rt*17.0e0_rt*gl*gum + 1.0e0_rt; + d = (a2 - fl*b2)*c; + fldt = d*gl2dt; + flda = d*gl2da; + fldz = d*gl2dz; + } // equation 4.9 and 4.10 cc = std::log10(2.0e0_rt*rm); xlnt = std::log10(temp); xnum = sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt); - xnumdt = -iln10*0.5e0_rt*tempi; - a2 = iln10*sixth*rmi; - xnumda = a2*rmda; - xnumdz = a2*rmdz; - xden = sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt); - xdendt = iln10*0.5e0_rt*tempi; - xdenda = a2*rmda; - xdendz = a2*rmdz; + if constexpr (do_derivatives) { + xnumdt = -iln10*0.5e0_rt*tempi; + a2 = iln10*sixth*rmi; + xnumda = a2*rmda; + xnumdz = a2*rmdz; + + xdendt = iln10*0.5e0_rt*tempi; + xdenda = a2*rmda; + xdendz = a2*rmdz; + } // equation 4.11 if (std::abs(xnum) > 0.7e0_rt || xden < 0.0e0_rt) { @@ -548,70 +607,82 @@ void sneut5(const Real temp, const Real den, } else { a1 = 0.39e0_rt - 1.25e0_rt*xnum - 0.35e0_rt*std::sin(4.5e0_rt*xnum); - a2 = -1.25e0_rt - 4.5e0_rt*0.35e0_rt*std::cos(4.5e0_rt*xnum); b1 = 0.3e0_rt * std::exp(-std::pow(4.5e0_rt*xnum + 0.9e0_rt, 2.0_rt)); - b2 = -b1*2.0e0_rt*(4.5e0_rt*xnum + 0.9e0_rt)*4.5e0_rt; c = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt*xnum); - if (c == 0.0_rt) { - dumdt = 0.0e0_rt; - dumda = 0.0e0_rt; - dumdz = 0.0e0_rt; - } else { - dumdt = xdendt + 1.25e0_rt*xnumdt; - dumda = xdenda + 1.25e0_rt*xnumda; - dumdz = xdendz + 1.25e0_rt*xnumdz; + + if constexpr (do_derivatives) { + a2 = -1.25e0_rt - 4.5e0_rt*0.35e0_rt*std::cos(4.5e0_rt*xnum); + b2 = -b1*2.0e0_rt*(4.5e0_rt*xnum + 0.9e0_rt)*4.5e0_rt; + if (c == 0.0_rt) { + dumdt = 0.0e0_rt; + dumda = 0.0e0_rt; + dumdz = 0.0e0_rt; + } else { + dumdt = xdendt + 1.25e0_rt*xnumdt; + dumda = xdenda + 1.25e0_rt*xnumda; + dumdz = xdendz + 1.25e0_rt*xnumdz; + } } d = 0.57e0_rt - 0.25e0_rt*xnum; a3 = c/d; c00 = std::exp(-a3*a3); - f1 = -c00*2.0e0_rt*a3/d; - c01 = f1*(dumdt + a3*0.25e0_rt*xnumdt); - c03 = f1*(dumda + a3*0.25e0_rt*xnumda); - c04 = f1*(dumdz + a3*0.25e0_rt*xnumdz); - fxy = 1.05e0_rt + (a1 - b1)*c00; - fxydt = (a2*xnumdt - b2*xnumdt)*c00 + (a1-b1)*c01; - fxyda = (a2*xnumda - b2*xnumda)*c00 + (a1-b1)*c03; - fxydz = (a2*xnumdz - b2*xnumdz)*c00 + (a1-b1)*c04; - + if constexpr (do_derivatives) { + f1 = -c00*2.0e0_rt*a3/d; + c01 = f1*(dumdt + a3*0.25e0_rt*xnumdt); + c03 = f1*(dumda + a3*0.25e0_rt*xnumda); + c04 = f1*(dumdz + a3*0.25e0_rt*xnumdz); + + fxydt = (a2*xnumdt - b2*xnumdt)*c00 + (a1-b1)*c01; + fxyda = (a2*xnumda - b2*xnumda)*c00 + (a1-b1)*c03; + fxydz = (a2*xnumdz - b2*xnumdz)*c00 + (a1-b1)*c04; + } } // equation 4.1 and 4.5 splas = (ft + fl) * fxy; - splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt; - splasda = (ftda + flda)*fxy + (ft+fl)*fxyda; - splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz; + if constexpr (do_derivatives) { + splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt; + splasda = (ftda + flda)*fxy + (ft+fl)*fxyda; + splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz; + } a2 = std::exp(-gl); - a3 = -0.5e0_rt*a2*gl*gum; a1 = splas; splas = a2*a1; - splasdt = a2*splasdt + a3*gl2dt*a1; - splasda = a2*splasda + a3*gl2da*a1; - splasdz = a2*splasdz + a3*gl2dz*a1; + if constexpr (do_derivatives) { + a3 = -0.5e0_rt*a2*gl*gum; + splasdt = a2*splasdt + a3*gl2dt*a1; + splasda = a2*splasda + a3*gl2da*a1; + splasdz = a2*splasdz + a3*gl2dz*a1; + } a2 = gl6; - a3 = 3.0e0_rt*gl6*gum; a1 = splas; splas = a2*a1; - splasdt = a2*splasdt + a3*gl2dt*a1; - splasda = a2*splasda + a3*gl2da*a1; - splasdz = a2*splasdz + a3*gl2dz*a1; + if constexpr (do_derivatives) { + a3 = 3.0e0_rt*gl6*gum; + splasdt = a2*splasdt + a3*gl2dt*a1; + splasda = a2*splasda + a3*gl2da*a1; + splasdz = a2*splasdz + a3*gl2dz*a1; + } a2 = 0.93153e0_rt * 3.0e21_rt * xl9; - a3 = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt; a1 = splas; splas = a2*a1; - splasdt = a2*splasdt + a3*a1; - splasda = a2*splasda; - splasdz = a2*splasdz; + if constexpr (do_derivatives) { + a3 = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt; + splasdt = a2*splasdt + a3*a1; + splasda = a2*splasda; + splasdz = a2*splasdz; + } // photoneutrino process section // for reactions like e- + gamma => e- + nu_e + nubar_e @@ -768,108 +839,136 @@ void sneut5(const Real temp, const Real den, + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2 + c03*cos3 + dd03*sin3 + c04*cos4 + dd04*sin4 + c05*cos5 + dd05*sin5 + 0.5e0_rt*c06*last; - - f0 = taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0e0_rt - + dd02*cos2*2.0e0_rt - c03*sin3*3.0e0_rt + dd03*cos3*3.0e0_rt - - c04*sin4*4.0e0_rt + dd04*cos4*4.0e0_rt - - c05*sin5*5.0e0_rt + dd05*cos5*5.0e0_rt) - - 0.5e0_rt*c06*xast*fac2*taudt; - a1 = 0.5e0_rt*c10 + c11*cos1 + dd11*sin1 + c12*cos2 + dd12*sin2 + c13*cos3 + dd13*sin3 + c14*cos4 + dd14*sin4 + c15*cos5 + dd15*sin5 + 0.5e0_rt*c16*last; - f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0e0_rt - + dd12*cos2*2.0e0_rt - c13*sin3*3.0e0_rt + dd13*cos3*3.0e0_rt - - c14*sin4*4.0e0_rt + dd14*cos4*4.0e0_rt - c15*sin5*5.0e0_rt - + dd15*cos5*5.0e0_rt) - 0.5e0_rt*c16*xast*fac2*taudt; a2 = 0.5e0_rt*c20 + c21*cos1 + dd21*sin1 + c22*cos2 + dd22*sin2 + c23*cos3 + dd23*sin3 + c24*cos4 + dd24*sin4 + c25*cos5 + dd25*sin5 + 0.5e0_rt*c26*last; - f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0e0_rt - + dd22*cos2*2.0e0_rt - c23*sin3*3.0e0_rt + dd23*cos3*3.0e0_rt - - c24*sin4*4.0e0_rt + dd24*cos4*4.0e0_rt - c25*sin5*5.0e0_rt - + dd25*cos5*5.0e0_rt) - 0.5e0_rt*c26*xast*fac2*taudt; + if constexpr (do_derivatives) { + f0 = taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0e0_rt + + dd02*cos2*2.0e0_rt - c03*sin3*3.0e0_rt + dd03*cos3*3.0e0_rt + - c04*sin4*4.0e0_rt + dd04*cos4*4.0e0_rt + - c05*sin5*5.0e0_rt + dd05*cos5*5.0e0_rt) + - 0.5e0_rt*c06*xast*fac2*taudt; + + f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0e0_rt + + dd12*cos2*2.0e0_rt - c13*sin3*3.0e0_rt + dd13*cos3*3.0e0_rt + - c14*sin4*4.0e0_rt + dd14*cos4*4.0e0_rt - c15*sin5*5.0e0_rt + + dd15*cos5*5.0e0_rt) - 0.5e0_rt*c16*xast*fac2*taudt; + + f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0e0_rt + + dd22*cos2*2.0e0_rt - c23*sin3*3.0e0_rt + dd23*cos3*3.0e0_rt + - c24*sin4*4.0e0_rt + dd24*cos4*4.0e0_rt - c25*sin5*5.0e0_rt + + dd25*cos5*5.0e0_rt) - 0.5e0_rt*c26*xast*fac2*taudt; + } + // equation 3.4 dum = a0 + a1*zeta + a2*zeta2; - dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0e0_rt*a2*zeta*zetadt; - dumda = a1*zetada + 2.0e0_rt*a2*zeta*zetada; - dumdz = a1*zetadz + 2.0e0_rt*a2*zeta*zetadz; + if constexpr (do_derivatives) { + dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0e0_rt*a2*zeta*zetadt; + dumda = a1*zetada + 2.0e0_rt*a2*zeta*zetada; + dumdz = a1*zetadz + 2.0e0_rt*a2*zeta*zetadz; + } z = std::exp(-cc*zeta); xnum = dum*z; - xnumdt = dumdt*z - dum*z*cc*zetadt; - xnumda = dumda*z - dum*z*cc*zetada; - xnumdz = dumdz*z - dum*z*cc*zetadz; + if constexpr (do_derivatives) { + xnumdt = dumdt*z - dum*z*cc*zetadt; + xnumda = dumda*z - dum*z*cc*zetada; + xnumdz = dumdz*z - dum*z*cc*zetadz; + } xden = zeta3 + 6.290e-3_rt*xlm1 + 7.483e-3_rt*xlm2 + 3.061e-4_rt*xlm3; dum = 3.0e0_rt*zeta2; - xdendt = dum*zetadt - xldt*(6.290e-3_rt*xlm2 - + 2.0e0_rt*7.483e-3_rt*xlm3 + 3.0e0_rt*3.061e-4_rt*xlm4); - xdenda = dum*zetada; - xdendz = dum*zetadz; - + if constexpr (do_derivatives) { + xdendt = dum*zetadt - xldt*(6.290e-3_rt*xlm2 + + 2.0e0_rt*7.483e-3_rt*xlm3 + 3.0e0_rt*3.061e-4_rt*xlm4); + xdenda = dum*zetada; + xdendz = dum*zetadz; + } dum = 1.0e0_rt/xden; fphot = xnum*dum; - fphotdt = (xnumdt - fphot*xdendt)*dum; - fphotda = (xnumda - fphot*xdenda)*dum; - fphotdz = (xnumdz - fphot*xdendz)*dum; + if constexpr (do_derivatives) { + fphotdt = (xnumdt - fphot*xdendt)*dum; + fphotda = (xnumda - fphot*xdenda)*dum; + fphotdz = (xnumdz - fphot*xdendz)*dum; + } // equation 3.3 a0 = 1.0e0_rt + 2.045e0_rt * xl; xnum = 0.666e0_rt*std::pow(a0, -2.066e0_rt); - xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*xldt; + if constexpr (do_derivatives) { + xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*xldt; + } dum = 1.875e8_rt*xl + 1.653e8_rt*xl2 + 8.499e8_rt*xl3 - 1.604e8_rt*xl4; - dumdt = xldt*(1.875e8_rt + 2.0e0_rt*1.653e8_rt*xl + 3.0e0_rt*8.499e8_rt*xl2 - - 4.0e0_rt*1.604e8_rt*xl3); + if constexpr (do_derivatives) { + dumdt = xldt*(1.875e8_rt + 2.0e0_rt*1.653e8_rt*xl + 3.0e0_rt*8.499e8_rt*xl2 + - 4.0e0_rt*1.604e8_rt*xl3); + } z = 1.0e0_rt/dum; xden = 1.0e0_rt + rm*z; - xdendt = -rm*z*z*dumdt; - xdenda = rmda*z; - xdendz = rmdz*z; + if constexpr (do_derivatives) { + xdendt = -rm*z*z*dumdt; + xdenda = rmda*z; + xdendz = rmdz*z; + } z = 1.0e0_rt/xden; qphot = xnum*z; - qphotdt = (xnumdt - qphot*xdendt)*z; + if constexpr (do_derivatives) { + qphotdt = (xnumdt - qphot*xdendt)*z; + } dum = -qphot*z; - qphotda = dum*xdenda; - qphotdz = dum*xdendz; + if constexpr (do_derivatives) { + qphotda = dum*xdenda; + qphotdz = dum*xdendz; + } // equation 3.2 sphot = xl5 * fphot; - sphotdt = 5.0e0_rt*xl4*xldt*fphot + xl5*fphotdt; - sphotda = xl5*fphotda; - sphotdz = xl5*fphotdz; + if constexpr (do_derivatives) { + sphotdt = 5.0e0_rt*xl4*xldt*fphot + xl5*fphotdt; + sphotda = xl5*fphotda; + sphotdz = xl5*fphotdz; + } a1 = sphot; sphot = rm*a1; - sphotdt = rm*sphotdt; - sphotda = rm*sphotda + rmda*a1; - sphotdz = rm*sphotdz + rmdz*a1; + if constexpr (do_derivatives) { + sphotdt = rm*sphotdt; + sphotda = rm*sphotda + rmda*a1; + sphotdz = rm*sphotdz + rmdz*a1; + } a1 = tfac4*(1.0e0_rt - tfac3 * qphot); - a2 = -tfac4*tfac3; a3 = sphot; sphot = a1*a3; - sphotdt = a1*sphotdt + a2*qphotdt*a3; - sphotda = a1*sphotda + a2*qphotda*a3; - sphotdz = a1*sphotdz + a2*qphotdz*a3; + if constexpr (do_derivatives) { + a2 = -tfac4*tfac3; + sphotdt = a1*sphotdt + a2*qphotdt*a3; + sphotda = a1*sphotda + a2*qphotda*a3; + sphotdz = a1*sphotdz + a2*qphotdz*a3; + } if (sphot <= 0.0_rt) { sphot = 0.0e0_rt; - sphotdt = 0.0e0_rt; - sphotda = 0.0e0_rt; - sphotdz = 0.0e0_rt; + if constexpr (do_derivatives) { + sphotdt = 0.0e0_rt; + sphotda = 0.0e0_rt; + sphotdz = 0.0e0_rt; + } } // bremsstrahlung neutrino section @@ -900,13 +999,17 @@ void sneut5(const Real temp, const Real den, // equation 5.3 dum = 7.05e6_rt * t832 + 5.12e4_rt * t83; - dumdt = (1.5e0_rt*7.05e6_rt*t812 + 3.0e0_rt*5.12e4_rt*t82)*1.0e-8_rt; + if constexpr (do_derivatives) { + dumdt = (1.5e0_rt*7.05e6_rt*t812 + 3.0e0_rt*5.12e4_rt*t82)*1.0e-8_rt; + } z = 1.0e0_rt/dum; eta = rm*z; - etadt = -rm*z*z*dumdt; - etada = rmda*z; - etadz = rmdz*z; + if constexpr (do_derivatives) { + etadt = -rm*z*z*dumdt; + etada = rmda*z; + etadz = rmdz*z; + } etam1 = 1.0e0_rt/eta; etam2 = etam1 * etam1; @@ -918,27 +1021,35 @@ void sneut5(const Real temp, const Real den, xnum = 1.0e0_rt/a0; dum = 1.0e0_rt + 1.47e0_rt*etam1 + 3.29e-2_rt*etam2; - z = -1.47e0_rt*etam2 - 2.0e0_rt*3.29e-2_rt*etam3; - dumdt = z*etadt; - dumda = z*etada; - dumdz = z*etadz; + if constexpr (do_derivatives) { + z = -1.47e0_rt*etam2 - 2.0e0_rt*3.29e-2_rt*etam3; + dumdt = z*etadt; + dumda = z*etada; + dumdz = z*etadz; + } c00 = 1.26e0_rt * (1.0e0_rt+etam1); - z = -1.26e0_rt*etam2; - c01 = z*etadt; - c03 = z*etada; - c04 = z*etadz; + if constexpr (do_derivatives) { + z = -1.26e0_rt*etam2; + c01 = z*etadt; + c03 = z*etada; + c04 = z*etadz; + } z = 1.0e0_rt/dum; xden = c00*z; - xdendt = (c01 - xden*dumdt)*z; - xdenda = (c03 - xden*dumda)*z; - xdendz = (c04 - xden*dumdz)*z; + if constexpr (do_derivatives) { + xdendt = (c01 - xden*dumdt)*z; + xdenda = (c03 - xden*dumda)*z; + xdendz = (c04 - xden*dumdz)*z; + } fbrem = xnum + xden; - fbremdt = -xnum*xnum*f0 + xdendt; - fbremda = xdenda; - fbremdz = xdendz; + if constexpr (do_derivatives) { + fbremdt = -xnum*xnum*f0 + xdendt; + fbremda = xdenda; + fbremdz = xdendz; + } // equation 5.9 a0 = 230.0e0_rt + 6.7e5_rt*t8m2 + 7.66e9_rt*t8m5; @@ -946,55 +1057,71 @@ void sneut5(const Real temp, const Real den, z = 1.0e0_rt + rm*1.0e-9_rt; dum = a0*z; - dumdt = f0*z; - z = a0*1.0e-9_rt; - dumda = z*rmda; - dumdz = z*rmdz; + if constexpr (do_derivatives) { + dumdt = f0*z; + z = a0*1.0e-9_rt; + dumda = z*rmda; + dumdz = z*rmdz; + } xnum = 1.0e0_rt/dum; - z = -xnum*xnum; - xnumdt = z*dumdt; - xnumda = z*dumda; - xnumdz = z*dumdz; + if constexpr (do_derivatives) { + z = -xnum*xnum; + xnumdt = z*dumdt; + xnumda = z*dumda; + xnumdz = z*dumdz; + } c00 = 7.75e5_rt*t832 + 247.0e0_rt * std::pow(t8, 3.85e0_rt); - dd00 = (1.5e0_rt*7.75e5_rt*t812 + 3.85e0_rt*247.0e0_rt*std::pow(t8, 2.85e0_rt))*1.0e-8_rt; - c01 = 4.07e0_rt + 0.0240e0_rt * std::pow(t8, 1.4e0_rt); - dd01 = 1.4e0_rt*0.0240e0_rt * std::pow(t8, 0.4e0_rt)*1.0e-8_rt; - c02 = 4.59e-5_rt * std::pow(t8, -0.110e0_rt); - dd02 = -0.11e0_rt*4.59e-5_rt * std::pow(t8, -1.11e0_rt)*1.0e-8_rt; + + if constexpr (do_derivatives) { + dd00 = (1.5e0_rt*7.75e5_rt*t812 + 3.85e0_rt*247.0e0_rt*std::pow(t8, 2.85e0_rt))*1.0e-8_rt; + + dd01 = 1.4e0_rt*0.0240e0_rt * std::pow(t8, 0.4e0_rt)*1.0e-8_rt; + dd02 = -0.11e0_rt*4.59e-5_rt * std::pow(t8, -1.11e0_rt)*1.0e-8_rt; + } z = std::pow(den, 0.656e0_rt); dum = c00*rmi + c01 + c02*z; - dumdt = dd00*rmi + dd01 + dd02*z; - z = -c00*rmi*rmi; - dumda = z*rmda; - dumdz = z*rmdz; + if constexpr (do_derivatives) { + dumdt = dd00*rmi + dd01 + dd02*z; + z = -c00*rmi*rmi; + dumda = z*rmda; + dumdz = z*rmdz; + } xden = 1.0e0_rt/dum; - z = -xden*xden; - xdendt = z*dumdt; - xdenda = z*dumda; - xdendz = z*dumdz; + if constexpr (do_derivatives) { + z = -xden*xden; + xdendt = z*dumdt; + xdenda = z*dumda; + xdendz = z*dumdz; + } gbrem = xnum + xden; - gbremdt = xnumdt + xdendt; - gbremda = xnumda + xdenda; - gbremdz = xnumdz + xdendz; + if constexpr (do_derivatives) { + gbremdt = xnumdt + xdendt; + gbremda = xnumda + xdenda; + gbremdz = xnumdz + xdendz; + } // equation 5.1 dum = 0.5738e0_rt*zbar*ye*t86*den; - dumdt = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt; - dumda = -dum*abari; - dumdz = 0.5738e0_rt*2.0e0_rt*ye*t86*den; + if constexpr (do_derivatives) { + dumdt = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt; + dumda = -dum*abari; + dumdz = 0.5738e0_rt*2.0e0_rt*ye*t86*den; + } z = tfac4*fbrem - tfac5*gbrem; sbrem = dum * z; - sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt); - sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda); - sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz); + if constexpr (do_derivatives) { + sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt); + sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda); + sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz); + } // liquid metal with c12 parameters (not too different for other elements) // equation 5.18 and 5.16 @@ -1078,9 +1205,11 @@ void sneut5(const Real temp, const Real den, // - 0.00069e0_rt*sin5*5.0e0_rt); dum = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*abari, oneth); - dumdt = -dum*tempi; - dumda = -oneth*dum*abari; - dumdz = 2.0e0_rt*dum*zbari; + if constexpr (do_derivatives) { + dumdt = -dum*tempi; + dumda = -oneth*dum*abari; + dumdz = 2.0e0_rt*dum*zbari; + } gm1 = 1.0e0_rt/dum; gm2 = gm1*gm1; @@ -1098,46 +1227,56 @@ void sneut5(const Real temp, const Real den, // equation 5.19 and 5.20 fliq = v*fb + (1.0e0_rt - v)*ft; - fliqdt = a0*dumdt*(fb - ft); - fliqda = a0*dumda*(fb - ft); - fliqdz = a0*dumdz*(fb - ft); + if constexpr (do_derivatives) { + fliqdt = a0*dumdt*(fb - ft); + fliqda = a0*dumda*(fb - ft); + fliqdz = a0*dumdz*(fb - ft); + } gliq = w*gb + (1.0e0_rt - w)*gt; - gliqdt = a1*dumdt*(gb - gt); - gliqda = a1*dumda*(gb - gt); - gliqdz = a1*dumdz*(gb - gt); + if constexpr (do_derivatives) { + gliqdt = a1*dumdt*(gb - gt); + gliqda = a1*dumda*(gb - gt); + gliqdz = a1*dumdz*(gb - gt); + } // equation 5.17 dum = 0.5738e0_rt*zbar*ye*t86*den; - dumdt = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt; - dumda = -dum*abari; - dumdz = 0.5738e0_rt*2.0e0_rt*ye*t86*den; + if constexpr (do_derivatives) { + dumdt = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt; + dumda = -dum*abari; + dumdz = 0.5738e0_rt*2.0e0_rt*ye*t86*den; + } z = tfac4*fliq - tfac5*gliq; sbrem = dum * z; - sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt); - sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda); - sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz); - + if constexpr (do_derivatives) { + sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt); + sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda); + sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz); + } } // recombination neutrino section // for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e // equation 6.11 solved for nu xnum = 1.10520e8_rt * den * ye /(temp*std::sqrt(temp)); - xnumdt = -1.50e0_rt*xnum*tempi; - xnumda = -xnum*abari; - xnumdz = xnum*zbari; + if constexpr (do_derivatives) { + xnumdt = -1.50e0_rt*xnum*tempi; + xnumda = -xnum*abari; + xnumdz = xnum*zbari; + } // the chemical potential nu = ifermi12(xnum); // a0 is d(nu)/d(xnum) a0 = 1.0e0_rt/(0.5e0_rt*zfermim12(nu)); - nudt = a0*xnumdt; - nuda = a0*xnumda; - nudz = a0*xnumdz; - + if constexpr (do_derivatives) { + nudt = a0*xnumdt; + nuda = a0*xnumda; + nudz = a0*xnumdz; + } nu2 = nu * nu; nu3 = nu2 * nu; @@ -1168,16 +1307,20 @@ void sneut5(const Real temp, const Real den, if (nu >= -20.0_rt && nu <= 10.0_rt) { zeta = 1.579e5_rt*zbar*zbar*tempi; - zetadt = -zeta*tempi; - zetada = 0.0e0_rt; - zetadz = 2.0e0_rt*zeta*zbari; + if constexpr (do_derivatives) { + zetadt = -zeta*tempi; + zetada = 0.0e0_rt; + zetadz = 2.0e0_rt*zeta*zbari; + } c00 = 1.0e0_rt/(1.0e0_rt + f1*nu + f2*nu2 + f3*nu3); c01 = f1 + f2*2.0e0_rt*nu + f3*3.0e0_rt*nu2; dum = zeta*c00; - dumdt = zetadt*c00 + zeta*c01*nudt; - dumda = zeta*c01*nuda; - dumdz = zetadz*c00 + zeta*c01*nudz; + if constexpr (do_derivatives) { + dumdt = zetadt*c00 + zeta*c01*nudt; + dumda = zeta*c01*nuda; + dumdz = zetadz*c00 + zeta*c01*nudz; + } z = 1.0e0_rt/dum; dd00 = std::pow(dum, -2.25_rt); @@ -1188,17 +1331,21 @@ void sneut5(const Real temp, const Real den, z = std::exp(c*nu); dd00 = b*z*(1.0e0_rt + d*dum); gum = 1.0e0_rt + dd00; - gumdt = dd00*c*nudt + b*z*d*dumdt; - gumda = dd00*c*nuda + b*z*d*dumda; - gumdz = dd00*c*nudz + b*z*d*dumdz; + if constexpr (do_derivatives) { + gumdt = dd00*c*nudt + b*z*d*dumdt; + gumda = dd00*c*nuda + b*z*d*dumda; + gumdz = dd00*c*nudz + b*z*d*dumdz; + } z = std::exp(nu); a1 = 1.0e0_rt/gum; bigj = c00 * z * a1; - bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt; - bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda; - bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz; + if constexpr (do_derivatives) { + bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt; + bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda; + bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz; + } // equation 6.5 z = std::exp(zeta + nu); @@ -1207,45 +1354,58 @@ void sneut5(const Real temp, const Real den, a2 = 1.0e0_rt/bigj; sreco = tfac6 * 2.649e-18_rt * ye * std::pow(zbar, 13.0_rt) * den * bigj*a1; - srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1); - srecoda = sreco*(-1.0e0_rt*abari + bigjda*a2 - z*(zetada+nuda)*a1); - srecodz = sreco*(14.0e0_rt*zbari + bigjdz*a2 - z*(zetadz+nudz)*a1); - + if constexpr (do_derivatives) { + srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1); + srecoda = sreco*(-1.0e0_rt*abari + bigjda*a2 - z*(zetada+nuda)*a1); + srecodz = sreco*(14.0e0_rt*zbari + bigjdz*a2 - z*(zetadz+nudz)*a1); + } } // convert from erg/cm^3/s to erg/g/s // comment these out to duplicate the itoh et al plots spair = spair*deni; - spairdt = spairdt*deni; - spairda = spairda*deni; - spairdz = spairdz*deni; + if constexpr (do_derivatives) { + spairdt = spairdt*deni; + spairda = spairda*deni; + spairdz = spairdz*deni; + } splas = splas*deni; - splasdt = splasdt*deni; - splasda = splasda*deni; - splasdz = splasdz*deni; + if constexpr (do_derivatives) { + splasdt = splasdt*deni; + splasda = splasda*deni; + splasdz = splasdz*deni; + } sphot = sphot*deni; - sphotdt = sphotdt*deni; - sphotda = sphotda*deni; - sphotdz = sphotdz*deni; + if constexpr (do_derivatives) { + sphotdt = sphotdt*deni; + sphotda = sphotda*deni; + sphotdz = sphotdz*deni; + } sbrem = sbrem*deni; - sbremdt = sbremdt*deni; - sbremda = sbremda*deni; - sbremdz = sbremdz*deni; + if constexpr (do_derivatives) { + sbremdt = sbremdt*deni; + sbremda = sbremda*deni; + sbremdz = sbremdz*deni; + } sreco = sreco*deni; - srecodt = srecodt*deni; - srecoda = srecoda*deni; - srecodz = srecodz*deni; + if constexpr (do_derivatives) { + srecodt = srecodt*deni; + srecoda = srecoda*deni; + srecodz = srecodz*deni; + } // the total neutrino loss rate snu = splas + spair + sphot + sbrem + sreco; - dsnudt = splasdt + spairdt + sphotdt + sbremdt + srecodt; - dsnuda = splasda + spairda + sphotda + sbremda + srecoda; - dsnudz = splasdz + spairdz + sphotdz + sbremdz + srecodz; + if constexpr (do_derivatives) { + dsnudt = splasdt + spairdt + sphotdt + sbremdt + srecodt; + dsnuda = splasda + spairda + sphotda + sbremda + srecoda; + dsnudz = splasdz + spairdz + sphotdz + sbremdz + srecodz; + } } diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H index 8ce878f196..1957dfe9e3 100644 --- a/nse_solver/nse_check.H +++ b/nse_solver/nse_check.H @@ -1070,7 +1070,8 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { composition(state); amrex::Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + constexpr int do_derivatives{0}; + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); #else amrex::Real sneut = 0.0_rt; #endif diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H index 55171960f0..39028c0603 100644 --- a/nse_solver/nse_solver.H +++ b/nse_solver/nse_solver.H @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -93,15 +94,6 @@ void apply_nse_exponent(T& nse_state) { #if SCREEN_METHOD != 0 - // three unit-less constants for calculating coulomb correction term - // see Calder 2007, doi:10.1086/510709 for more detail - - const amrex::Real A1 = -0.9052_rt; - const amrex::Real A2 = 0.6322_rt; - - // A3 = -0.5 * sqrt(3) - A1/sqrt(A2) - const amrex::Real A3 = 0.272433346667; - // Find n_e for original state; const amrex::Real n_e = nse_state.rho * nse_state.y_e / C::m_u; const amrex::Real Gamma_e = C::q_e * C::q_e * std::cbrt(4.0_rt * M_PI * n_e / 3.0_rt) @@ -125,12 +117,14 @@ void apply_nse_exponent(T& nse_state) { gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e; // chemical potential for coulomb correction + // see appendix of Calder 2007, doi:10.1086/510709 for more detail + + // reuse existing implementation from screening routine + Real f, df; + constexpr int do_T_derivatives = 0; + chabrier1998_helmholtz_F(gamma, 0.0_rt, f, df); - u_c = C::k_B * T_in / C::Legacy::MeV2erg * - (A1 * (std::sqrt(gamma * (A2 + gamma)) - - A2 * std::log(std::sqrt(gamma / A2) + - std::sqrt(1.0_rt + gamma / A2))) + - 2.0_rt * A3 * (std::sqrt(gamma) - std::atan(std::sqrt(gamma)))); + u_c = C::k_B * T_in / C::Legacy::MeV2erg * f; #endif // find nse mass frac diff --git a/screening/screen.H b/screening/screen.H index 40ced220bd..898772c8d8 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -613,7 +613,7 @@ void chugunov2007 (const plasma_state_t& state, template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT) +void chugunov2009_f0 (const Real gamma, const Real dlog_dT, Real& f, Real& df_dT) { // Calculates the free energy per ion in a OCP, from Chugunov and DeWitt 2009 // equation 24. @@ -742,9 +742,9 @@ void chugunov2009 (const plasma_state_t& state, // values similar to those from Chugunov 2007. Real term1, term2, term3; Real dterm1_dT = 0.0_rt, dterm2_dT = 0.0_rt, dterm3_dT = 0.0_rt; - f0(Gamma_1 / t_12, dlog_dT, term1, dterm1_dT); - f0(Gamma_2 / t_12, dlog_dT, term2, dterm2_dT); - f0(Gamma_comp / t_12, dlog_dT, term3, dterm3_dT); + chugunov2009_f0(Gamma_1 / t_12, dlog_dT, term1, dterm1_dT); + chugunov2009_f0(Gamma_2 / t_12, dlog_dT, term2, dterm2_dT); + chugunov2009_f0(Gamma_comp / t_12, dlog_dT, term3, dterm3_dT); Real h_fit = term1 + term2 - term3; Real dh_fit_dT; if constexpr (do_T_derivatives) { @@ -780,6 +780,35 @@ void chugunov2009 (const plasma_state_t& state, } } +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void chabrier1998_helmholtz_F(const Real gamma, const Real dgamma_dT, Real& f, Real& df_dT) { + // Helmholtz free energy, See Chabrier & Potekhin 1998 Eq. 28 + + // Fitted parameters, see Chabrier & Potekhin 1998 Sec.IV + constexpr Real A_1 = -0.9052_rt; + constexpr Real A_2 = 0.6322_rt; + constexpr Real A_3 = -0.5_rt * gcem::sqrt(3.0_rt) - A_1 / gcem::sqrt(A_2); + + // Precompute some expressions that are reused in the derivative + const Real sqrt_gamma = std::sqrt(gamma); + const Real sqrt_1_gamma_A2 = std::sqrt(1.0_rt + gamma/A_2); + const Real sqrt_gamma_A2_gamma = std::sqrt(gamma * (A_2 + gamma)); + const Real sqrt_gamma_A2 = std::sqrt(gamma/A_2); + + f = A_1 * (sqrt_gamma_A2_gamma - + A_2 * std::log(sqrt_gamma_A2 + sqrt_1_gamma_A2)) + + 2.0_rt * A_3 * (sqrt_gamma - std::atan(sqrt_gamma)); + + if constexpr (do_T_derivatives) { + df_dT = A_1 * ((A_2 + 2.0_rt * gamma) / (2.0_rt * sqrt_gamma_A2_gamma) + - 0.5_rt / (sqrt_gamma_A2 + sqrt_1_gamma_A2) + * (1.0_rt / sqrt_gamma_A2 + 1.0_rt / sqrt_1_gamma_A2)) + + A_3 / sqrt_gamma * (1.0_rt - 1.0_rt / (1.0_rt + gamma)); + + df_dT *= dgamma_dT; + } +} template AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -809,57 +838,22 @@ void chabrier1998 (const plasma_state_t& state, Real Gamma2 = Gamma_e * std::pow(scn_fac.z2, 5.0_rt/3.0_rt); Real Gamma12 = Gamma_e * std::pow(zcomp, 5.0_rt/3.0_rt); - [[maybe_unused]] Real Gamma1dT, Gamma2dT, Gamma12dT; + Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{}; if constexpr (do_T_derivatives) { Gamma1dT = -Gamma1 / state.temp; Gamma2dT = -Gamma2 / state.temp; Gamma12dT = -Gamma12 / state.temp; } - // Fitted parameters, see Chabrier & Potekhin 1998 Sec.IV - - const Real A_1 = -0.9052_rt; - const Real A_2 = 0.6322_rt; - const Real A_3 = -0.5_rt * std::sqrt(3.0_rt) - A_1 / std::sqrt(A_2); - - // Helmholtz free energy, See Chabrier & Potekhin 1998 Eq. 28 - - Real f1 = A_1 * (std::sqrt(Gamma1 * (A_2 + Gamma1)) - - A_2 * std::log(std::sqrt(Gamma1/A_2) + std::sqrt(1.0_rt + Gamma1/A_2))) - + 2.0_rt * A_3 * (std::sqrt(Gamma1) - std::atan(std::sqrt(Gamma1))); - - Real f2 = A_1 * (std::sqrt(Gamma2 * (A_2 + Gamma2)) - - A_2 * std::log(std::sqrt(Gamma2/A_2) + std::sqrt(1.0_rt + Gamma2/A_2))) - + 2.0_rt * A_3 * (std::sqrt(Gamma2) - std::atan(std::sqrt(Gamma2))); - - Real f12 = A_1 * (std::sqrt(Gamma12 * (A_2 + Gamma12)) - - A_2 * std::log(std::sqrt(Gamma12/A_2) + std::sqrt(1.0_rt + Gamma12/A_2))) - + 2.0_rt * A_3 * (std::sqrt(Gamma12) - std::atan(std::sqrt(Gamma12))); - [[maybe_unused]] Real f1dT, f2dT, f12dT; - - if constexpr (do_T_derivatives) { - f1dT = A_1 * ((A_2 + 2.0_rt * Gamma1) / (2.0_rt * std::sqrt(Gamma1 * (A_2 + Gamma1))) - - 0.5_rt / (std::sqrt(Gamma1 / A_2) + std::sqrt(1.0_rt + Gamma1 / A_2)) - * (1.0_rt / std::sqrt(Gamma1 / A_2) + 1.0_rt / std::sqrt(1.0_rt + Gamma1 / A_2))) - + A_3 / std::sqrt(Gamma1) * (1.0_rt - 1.0_rt / (1.0_rt + Gamma1)); - - f1dT *= Gamma1dT; + // Helmholtz free energy - f2dT = A_1 * ((A_2 + 2.0_rt * Gamma2) / (2.0_rt * std::sqrt(Gamma2 * (A_2 + Gamma2))) - - 0.5_rt / (std::sqrt(Gamma2 / A_2) + std::sqrt(1.0_rt + Gamma2 / A_2)) - * (1.0_rt / std::sqrt(Gamma2 / A_2) + 1.0_rt / std::sqrt(1.0_rt + Gamma2 / A_2))) - + A_3 / std::sqrt(Gamma2) * (1.0_rt - 1.0_rt / (1.0_rt + Gamma2)); + Real f1, f2, f12; + Real f1dT, f2dT, f12dT; - f2dT *= Gamma2dT; - - f12dT = A_1 * ((A_2 + 2.0_rt * Gamma12) / (2.0_rt * std::sqrt(Gamma12 * (A_2 + Gamma12))) - - 0.5_rt / (std::sqrt(Gamma12 / A_2) + std::sqrt(1.0_rt + Gamma12 / A_2)) - * (1.0_rt / std::sqrt(Gamma12 / A_2) + 1.0_rt / std::sqrt(1.0_rt + Gamma12 / A_2))) - + A_3 / std::sqrt(Gamma12) * (1.0_rt - 1.0_rt / (1.0_rt + Gamma12)); - - f12dT *= Gamma12dT; - } + chabrier1998_helmholtz_F(Gamma1, Gamma1dT, f1, f1dT); + chabrier1998_helmholtz_F(Gamma2, Gamma2dT, f2, f2dT); + chabrier1998_helmholtz_F(Gamma12, Gamma12dT, f12, f12dT); // Now we add quantum correction terms discussed in Alastuey 1978. // Notice in Alastuey 1978, they have a different classical term, diff --git a/unit_test/burn_cell_sdc/_parameters b/unit_test/burn_cell_sdc/_parameters index 7d0a23f6f9..3343cff5f7 100644 --- a/unit_test/burn_cell_sdc/_parameters +++ b/unit_test/burn_cell_sdc/_parameters @@ -69,3 +69,5 @@ Adv_Aux3 real 0.0e0 +mu_p real -5.0 +mu_n real -12.0 diff --git a/unit_test/burn_cell_sdc/burn_cell.H b/unit_test/burn_cell_sdc/burn_cell.H index 78705a9bdd..c83752ecfa 100644 --- a/unit_test/burn_cell_sdc/burn_cell.H +++ b/unit_test/burn_cell_sdc/burn_cell.H @@ -293,6 +293,12 @@ void burn_cell_c() burn_state.T_fixed = -1.0; +#ifdef NSE_NET + burn_state.mu_p = unit_test_rp::mu_p; + burn_state.mu_n = unit_test_rp::mu_n; + burn_state.y_e = -1.0; +#endif + std::ofstream state_over_time("state_over_time.txt"); // we will divide the total integration time into nsteps that are diff --git a/unit_test/test_neutrino_cooling/GNUmakefile b/unit_test/test_neutrino_cooling/GNUmakefile new file mode 100644 index 0000000000..024126c7f5 --- /dev/null +++ b/unit_test/test_neutrino_cooling/GNUmakefile @@ -0,0 +1,43 @@ +PRECISION = DOUBLE +PROFILE = FALSE + +DEBUG = FALSE + +DIM = 3 + +COMP = gnu + +USE_MPI = FALSE +USE_OMP = FALSE + +USE_REACT = TRUE +USE_CONDUCTIVITY = FALSE + +EBASE = main + +BL_NO_FORT = TRUE + +# define the location of the Microphysics top directory +MICROPHYSICS_HOME := ../.. + +# This sets the EOS directory +EOS_DIR := helmholtz + +# This sets the network directory +NETWORK_DIR := aprox21 + +# This isn't actually used but we need VODE to compile with CUDA +INTEGRATOR_DIR := VODE + +CONDUCTIVITY_DIR := stellar + +SCREEN_METHOD := screen5 + +EXTERN_SEARCH += . + +Bpack := ./Make.package +Blocs := . + +include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test + + diff --git a/unit_test/test_neutrino_cooling/Make.package b/unit_test/test_neutrino_cooling/Make.package new file mode 100644 index 0000000000..c2f94cfc92 --- /dev/null +++ b/unit_test/test_neutrino_cooling/Make.package @@ -0,0 +1,7 @@ +CEXE_sources += main.cpp + +CEXE_headers += test_sneut.H + +CEXE_sources += variables.cpp +CEXE_sources += neutrino_util.cpp +CEXE_headers += variables.H diff --git a/unit_test/test_neutrino_cooling/README.md b/unit_test/test_neutrino_cooling/README.md new file mode 100644 index 0000000000..b7e52c275b --- /dev/null +++ b/unit_test/test_neutrino_cooling/README.md @@ -0,0 +1,4 @@ +# test_neutrino_cooling + +Test the neutrino cooling routine, sneut5 + diff --git a/unit_test/test_neutrino_cooling/_parameters b/unit_test/test_neutrino_cooling/_parameters new file mode 100644 index 0000000000..b2a919b35e --- /dev/null +++ b/unit_test/test_neutrino_cooling/_parameters @@ -0,0 +1,11 @@ +@namespace: unit_test + +dens_min real 1.d6 +dens_max real 1.d9 +temp_min real 1.d6 +temp_max real 1.d12 + +metalicity_max real 0.1d0 + +small_temp real 1.d4 +small_dens real 1.d-4 diff --git a/unit_test/test_neutrino_cooling/ci-benchmarks/test_sneut5.out b/unit_test/test_neutrino_cooling/ci-benchmarks/test_sneut5.out new file mode 100644 index 0000000000..c314665d07 --- /dev/null +++ b/unit_test/test_neutrino_cooling/ci-benchmarks/test_sneut5.out @@ -0,0 +1,31 @@ + plotfile = test_sneut5 + time = 0 + variables minimum value maximum value + density 10 5000000000 + temperature 1000000 10000000000 + X_hydrogen-1 0.5 0.75 + X_helium-3 0 0.026315789474 + X_helium-4 0 0.25 + X_carbon-12 0 0.026315789474 + X_nitrogen-14 0 0.026315789474 + X_oxygen-16 0 0.026315789474 + X_neon-20 0 0.026315789474 + X_magnesium-24 0 0.026315789474 + X_silicon-28 0 0.026315789474 + X_sulfur-32 0 0.026315789474 + X_argon-36 0 0.026315789474 + X_calcium-40 0 0.026315789474 + X_titanium-44 0 0.026315789474 + X_chromium-48 0 0.026315789474 + X_chromium-56 0 0.026315789474 + X_iron-52 0 0.026315789474 + X_iron-54 0 0.026315789474 + X_iron-56 0 0.026315789474 + X_nickel-56 0 0.026315789474 + X_neutron 0 0.026315789474 + X_proton 0 0.026315789474 + sneut 0 3.3206297956e+23 + dsneutdt 0 2.9114557633e+14 + dsneutda -2.0028760264e+16 1.6034657416e+20 + dsneutdz -1.832532276e+20 2.2890025857e+16 + diff --git a/unit_test/test_neutrino_cooling/inputs b/unit_test/test_neutrino_cooling/inputs new file mode 100644 index 0000000000..bd2bdb9698 --- /dev/null +++ b/unit_test/test_neutrino_cooling/inputs @@ -0,0 +1,9 @@ +n_cell = 16 +max_grid_size = 32 + +unit_test.dens_min = 10.e0 +unit_test.dens_max = 5.e9 +unit_test.temp_min = 1.e6 +unit_test.temp_max = 1.e10 + +unit_test.metalicity_max = 0.5e0 diff --git a/unit_test/test_neutrino_cooling/main.cpp b/unit_test/test_neutrino_cooling/main.cpp new file mode 100644 index 0000000000..24d498f520 --- /dev/null +++ b/unit_test/test_neutrino_cooling/main.cpp @@ -0,0 +1,166 @@ +#include +#include +#include + +#include +#include +#include + + +using namespace amrex; + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +using namespace unit_test_rp; + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + + main_main(); + + amrex::Finalize(); + return 0; +} + +void main_main () +{ + + // AMREX_SPACEDIM: number of dimensions + int n_cell, max_grid_size; + + // inputs parameters + { + // ParmParse is way of reading inputs from the inputs file + ParmParse pp; + + // We need to get n_cell from the inputs file - this is the + // number of cells on each side of a square (or cubic) domain. + pp.get("n_cell", n_cell); + + // The domain is broken into boxes of size max_grid_size + max_grid_size = 32; + pp.query("max_grid_size", max_grid_size); + + } + + + Vector is_periodic(AMREX_SPACEDIM,0); + for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { + is_periodic[idim] = 1; + } + + // make BoxArray and Geometry + BoxArray ba; + Geometry geom; + { + IntVect dom_lo(AMREX_D_DECL( 0, 0, 0)); + IntVect dom_hi(AMREX_D_DECL(n_cell-1, n_cell-1, n_cell-1)); + Box domain(dom_lo, dom_hi); + + // Initialize the boxarray "ba" from the single box "bx" + ba.define(domain); + + // Break up boxarray "ba" into chunks no larger than + // "max_grid_size" along a direction + ba.maxSize(max_grid_size); + + // This defines the physical box, [0, 1] in each direction. + RealBox real_box({AMREX_D_DECL(0.0, 0.0, 0.0)}, + {AMREX_D_DECL(1.0, 1.0, 1.0)}); + + // This defines a Geometry object + geom.define(domain, &real_box, + CoordSys::cartesian, is_periodic.data()); + } + + // Nghost = number of ghost cells for each array + int Nghost = 0; + + // do the runtime parameter initializations and microphysics inits + if (ParallelDescriptor::IOProcessor()) { + std::cout << "reading extern runtime parameters ..." << std::endl; + } + + ParmParse ppa("amr"); + + init_unit_test(); + + eos_init(small_temp, small_dens); + + network_init(); + + plot_t vars; + + vars = init_variables(); + + amrex::Vector names; + get_varnames(vars, names); + + // time = starting time in the simulation + Real time = 0.0_rt; + + // How Boxes are distributed among MPI processes + DistributionMapping dm(ba); + + // we allocate our main multifabs + MultiFab state(ba, dm, vars.n_plot_comps, Nghost); + + // Initialize the state to zero; we will fill + // it in below in do_eos. + state.setVal(0.0_rt); + + // What time is it now? We'll use this to compute total run time. + Real strt_time = ParallelDescriptor::second(); + + + Real dlogrho = 0.0e0_rt; + Real dlogT = 0.0e0_rt; + Real dmetal = 0.0e0_rt; + + if (n_cell > 1) { + dlogrho = (std::log10(dens_max) - std::log10(dens_min)) / static_cast(n_cell - 1); + dlogT = (std::log10(temp_max) - std::log10(temp_min))/ static_cast(n_cell - 1); + dmetal = (metalicity_max - 0.0_rt)/ static_cast(n_cell - 1); + } + + // Initialize the state and compute the different thermodynamics + // by inverting the EOS + for ( MFIter mfi(state); mfi.isValid(); ++mfi ) + { + const Box& bx = mfi.validbox(); + + Array4 const sp = state.array(mfi); + + neut_test_C(bx, dlogrho, dlogT, dmetal, vars, sp); + + } + + // Call the timer again and compute the maximum difference between + // the start time and stop time over all processors + Real stop_time = ParallelDescriptor::second() - strt_time; + const int IOProc = ParallelDescriptor::IOProcessorNumber(); + ParallelDescriptor::ReduceRealMax(stop_time, IOProc); + + + std::string name = "test_sneut5"; + + // Write a plotfile + WriteSingleLevelPlotfile(name, state, names, geom, time, 0); + + write_job_info(name); + + // Tell the I/O Processor to write out the "run time" + amrex::Print() << "Run time = " << stop_time << std::endl; + +} diff --git a/unit_test/test_neutrino_cooling/neutrino_util.cpp b/unit_test/test_neutrino_cooling/neutrino_util.cpp new file mode 100644 index 0000000000..efa6c70bb4 --- /dev/null +++ b/unit_test/test_neutrino_cooling/neutrino_util.cpp @@ -0,0 +1,95 @@ +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include + +using namespace amrex; +using namespace unit_test_rp; + +void neut_test_C(const Box& bx, + const Real dlogrho, const Real dlogT, const Real dmetal, + const plot_t& vars, + Array4 const sp) { + + using namespace Species; + + const int ih1 = network_spec_index("hydrogen-1"); + if (ih1 < 0) amrex::Error("Error: ih1 not found"); + + const int ihe4 = network_spec_index("helium-4"); + if (ihe4 < 0) amrex::Error("Error: ihe4 not found"); + + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + + // set the composition -- approximately solar + Real metalicity = 0.0 + static_cast (k) * dmetal; + + Real xn[NumSpec]; + + // for now... the screening using 1-based indexing + Array1D ymass; + + for (int n = 0; n < NumSpec; n++) { + xn[n] = metalicity / static_cast(NumSpec - 2); + } + xn[ih1] = 0.75_rt - 0.5_rt * metalicity; + xn[ihe4] = 0.25_rt - 0.5_rt * metalicity; + + Real temp_zone = std::pow(10.0, std::log10(temp_min) + static_cast(j)*dlogT); + + Real dens_zone = std::pow(10.0, std::log10(dens_min) + static_cast(i)*dlogrho); + + // store default state + sp(i, j, k, vars.irho) = dens_zone; + sp(i, j, k, vars.itemp) = temp_zone; + for (int n = 0; n < NumSpec; n++) { + sp(i, j, k, vars.ispec+n) = xn[n]; + } + + // compute abar and zbar + + Real abar = 0.0; + for (int n = 0; n < NumSpec; n++) { + abar += xn[n] / aion[n]; + } + abar = 1.0_rt / abar; + + Real zbar = 0.0; + for (int n = 0; n < NumSpec; n++) { + zbar += zion[n] * xn[n] / aion[n]; + } + zbar *= abar; + + Real snu; + Real dsnudt; + Real dsnudd; + Real dsnuda; + Real dsnudz; + + constexpr int do_derivatives{1}; + + sneut5(temp_zone, dens_zone, abar, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); + + sp(i, j, k, vars.isneut) = snu; + sp(i, j, k, vars.isneutdt) = dsnudt; + sp(i, j, k, vars.isneutda) = dsnuda; + sp(i, j, k, vars.isneutdz) = dsnudz; + + }); + +} diff --git a/unit_test/test_neutrino_cooling/test_sneut.H b/unit_test/test_neutrino_cooling/test_sneut.H new file mode 100644 index 0000000000..89025967e3 --- /dev/null +++ b/unit_test/test_neutrino_cooling/test_sneut.H @@ -0,0 +1,14 @@ +#ifndef TEST_SCREEN_H +#define TEST_SCREEN_H + +#include +#include + +void main_main(); + +void neut_test_C(const Box& bx, + const Real dlogrho, const Real dlogT, const Real dmetal, + const plot_t& vars, + Array4 const sp); + +#endif diff --git a/unit_test/test_neutrino_cooling/variables.H b/unit_test/test_neutrino_cooling/variables.H new file mode 100644 index 0000000000..6444cf0e3d --- /dev/null +++ b/unit_test/test_neutrino_cooling/variables.H @@ -0,0 +1,37 @@ +#include +#include + +#ifndef VARIABLES_H +#define VARIABLES_H + +#include + +class plot_t { + +public: + + + int irho = -1; + int itemp = -1; + int ispec = -1; + int isneut = -1; + int isneutdt = -1; + int isneutda = -1; + int isneutdz = -1; + + int n_plot_comps = 0; + + int next_index(const int num) { + int next = n_plot_comps; + n_plot_comps += num; + return next; + } + +}; + +plot_t init_variables(); + +void get_varnames(const plot_t& p, amrex::Vector& names); + + +#endif diff --git a/unit_test/test_neutrino_cooling/variables.cpp b/unit_test/test_neutrino_cooling/variables.cpp new file mode 100644 index 0000000000..98eed26b98 --- /dev/null +++ b/unit_test/test_neutrino_cooling/variables.cpp @@ -0,0 +1,36 @@ +#include +#include + +plot_t init_variables() { + + plot_t p; + + p.irho = p.next_index(1); + p.itemp = p.next_index(1); + p.ispec = p.next_index(NumSpec); + p.isneut = p.next_index(1); + p.isneutdt = p.next_index(1); + p.isneutda = p.next_index(1); + p.isneutdz = p.next_index(1); + + return p; +} + + +void get_varnames(const plot_t& p, amrex::Vector& names) { + + names.resize(p.n_plot_comps); + + names[p.irho] = "density"; + names[p.itemp] = "temperature"; + for (int n = 0; n < NumSpec; n++) { + names[p.ispec + n] = "X_" + spec_names_cxx[n]; + } + + names[p.isneut] = "sneut"; + names[p.isneutdt] = "dsneutdt"; + names[p.isneutda] = "dsneutda"; + names[p.isneutdz] = "dsneutdz"; +} + +