diff --git a/networks/ECSN/actual_network.H b/networks/ECSN/actual_network.H index b84447dcac..191b80d76d 100644 --- a/networks/ECSN/actual_network.H +++ b/networks/ECSN/actual_network.H @@ -93,7 +93,6 @@ namespace NSE_INDEX // last index is the corresponding reverse rate index. extern AMREX_GPU_MANAGED amrex::Array2D rate_indices; - extern AMREX_GPU_MANAGED bool initialized; } #endif diff --git a/networks/ECSN/actual_network_data.cpp b/networks/ECSN/actual_network_data.cpp index bebbfe6259..a6fab43815 100644 --- a/networks/ECSN/actual_network_data.cpp +++ b/networks/ECSN/actual_network_data.cpp @@ -30,7 +30,6 @@ namespace NSE_INDEX -1, -1, 3, -1, -1, 4, 15, -1, -1, 4, -1, -1, 5, 16 }; - AMREX_GPU_MANAGED bool initialized = false; } #endif diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H index aa4cc06caf..b7aa913121 100644 --- a/networks/ECSN/actual_rhs.H +++ b/networks/ECSN/actual_rhs.H @@ -70,7 +70,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { Real ratraw, dratraw_dT; Real scor, dscor_dt; - Real scor2, dscor2_dt; + [[maybe_unused]] Real scor2, dscor2_dt; { @@ -411,9 +411,9 @@ void actual_rhs (burn_t& state, Array1D& ydot) // Get the thermal neutrino losses - Real sneut, dsneutdt, dsneutdd, snuda, snudz; + Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz); // Append the energy equation (this is erg/g/s) @@ -628,11 +628,11 @@ 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); + Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz; + 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 * snuda + (zion[j-1] - state.zbar) * state.abar * snudz); + Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz); jac.add(net_ienuc, j, -b1); } diff --git a/networks/ECSN/partition_functions.H b/networks/ECSN/partition_functions.H index 7b6a9030ed..904f1e9397 100644 --- a/networks/ECSN/partition_functions.H +++ b/networks/ECSN/partition_functions.H @@ -22,18 +22,25 @@ namespace part_fun { void interpolate_pf(const Real t9, const Real (&temp_array)[npts], const Real (&pf_array)[npts], Real& pf, Real& dpf_dT) { - // find the index of the first temperature element we are larger than + if (t9 >= temp_array[0] && t9 < temp_array[npts-1]) { - int idx = -1; + // find the largest temperature element <= t9 using a binary search - for (int i = 0; i < npts-1; ++i) { - if (t9 >= temp_array[i] && t9 < temp_array[i+1]) { - idx = i; - break; + int left = 0; + int right = npts; + + while (left < right) { + int mid = (left + right) / 2; + if (temp_array[mid] > t9) { + right = mid; + } else { + left = mid + 1; + } } - } - if (idx >= 0) { + const int idx = right - 1; + + // now we have temp_array[idx] <= t9 < temp_array[idx+1] // construct the slope -- this is (log10(pf_{i+1}) - log10(pf_i)) / (T_{i+1} - T_i) @@ -46,12 +53,12 @@ namespace part_fun { // find the derivative (with respect to T, not T9) - Real dpf_dT9 = pf * std::log(10.0_rt) * slope; + Real dpf_dT9 = pf * M_LN10 * slope; dpf_dT = dpf_dT9 / 1.e9_rt; } else { - // T < the smallest T in the partition function table + // T < the smallest T or >= the largest T in the partition function table pf = 1.0; dpf_dT = 0.0; @@ -90,48 +97,27 @@ constexpr Real get_spin_state(const int inuc) { switch (inuc) { - case H1: - spin = 2; - break; - case He4: - spin = 1; - break; - case O16: - spin = 1; - break; - case O20: - spin = 1; - break; - - case F20: - spin = 5; - break; - case Ne20: - spin = 1; - break; - case Mg24: - spin = 1; - break; - - case Al27: - spin = 6; - break; - case Si28: + case S32: spin = 1; break; + case H1: case P31: spin = 2; break; - case S32: - spin = 1; + case F20: + spin = 5; + break; + + case Al27: + spin = 6; break; diff --git a/networks/ECSN/reaclib_rates.H b/networks/ECSN/reaclib_rates.H index 02c805a332..efbfe0e818 100644 --- a/networks/ECSN/reaclib_rates.H +++ b/networks/ECSN/reaclib_rates.H @@ -910,11 +910,11 @@ fill_reaclib_rates(const tf_t& tfactors, T& rate_eval) template AMREX_GPU_HOST_DEVICE AMREX_INLINE void -fill_approx_rates([[maybe_unused]] const tf_t& tfactors, T& rate_eval) +fill_approx_rates([[maybe_unused]] const tf_t& tfactors, [[maybe_unused]] T& rate_eval) { - Real rate; - Real drate_dT; + [[maybe_unused]] Real rate{}; + [[maybe_unused]] Real drate_dT{}; } diff --git a/networks/ECSN/table_rates.H b/networks/ECSN/table_rates.H index 8df50f1fa1..ed4cbd66c5 100644 --- a/networks/ECSN/table_rates.H +++ b/networks/ECSN/table_rates.H @@ -112,6 +112,10 @@ void init_tab_info(const table_t& tf, const std::string& file, R& rhoy, T& temp, } if (1 == j) { + if (temp(i) > 33.0_rt){ + amrex::Error("init_tab_info: Table data needs to be updated"); + } + temp(i) = std::pow(10.0_rt, temp(i)); } @@ -151,11 +155,11 @@ int vector_index_lu(const int vlen, const V& vector, const Real fvar) } else { ndn = j; } - if (((nup - ndn) == 1)) { - index = ndn; - return index; - } + if ((nup - ndn) == 1) { + break; + } } + index = ndn; } return index; @@ -256,6 +260,7 @@ get_entries(const table_t& table_meta, data(itemp_hi, irhoy_hi, ivar), temp); + // NOLINTNEXTLINE(readability-suspicious-call-argument) entries(ivar) = bl_extrap(rhoy_lo, rhoy_hi, f_i, f_ip1, rhoy); }