From 05cbbda50b2dc544ebd65897d6d9760ba9cfe5ac Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Fri, 15 Sep 2023 15:51:57 -0400 Subject: [PATCH 1/4] Replace tabs with spaces --- integration/VODE/vode_dvstep.H | 6 +- interfaces/burner.H | 2 +- nse_solver/make_table/burn_cell.H | 38 +- nse_solver/nse_check.H | 540 +++++++++++----------- nse_solver/nse_solver.H | 72 +-- screening/screen.H | 50 +- unit_test/burn_cell/burn_cell.H | 16 +- unit_test/burn_cell_sdc/burn_cell.H | 20 +- unit_test/test_ase/make_table/burn_cell.H | 52 +-- unit_test/test_react/main.cpp | 4 +- util/hybrj/hybrj.H | 26 +- util/hybrj/hybrj_dogleg.H | 12 +- util/hybrj/hybrj_enorm.H | 6 +- util/hybrj/hybrj_r1mpyq.H | 8 +- util/microphysics_sort.H | 6 +- 15 files changed, 429 insertions(+), 429 deletions(-) diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index 0586e3e667..e5602dd80e 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -225,9 +225,9 @@ int dvstep (BurnT& state, DvodeT& vstate) bool valid_update = true; #ifdef SDC - if (vstate.y(SEINT+1) < 0.0_rt) { - valid_update = false; - } + if (vstate.y(SEINT+1) < 0.0_rt) { + valid_update = false; + } Real rho_current = state.rho_orig + vstate.tn * state.ydot_a[SRHO]; diff --git a/interfaces/burner.H b/interfaces/burner.H index 832edede73..94a725fb8d 100644 --- a/interfaces/burner.H +++ b/interfaces/burner.H @@ -63,7 +63,7 @@ void burner (BurnT& state, Real dt) #ifdef NSE_TABLE if (in_nse(state, true) && state.success == false && dt_remaining > 0.0) { #else - if (in_nse(state, nse_skip_molar) && state.success == false && dt_remaining > 0.0) { + if (in_nse(state, nse_skip_molar) && state.success == false && dt_remaining > 0.0) { #endif #ifndef AMREX_USE_GPU diff --git a/nse_solver/make_table/burn_cell.H b/nse_solver/make_table/burn_cell.H index 4afa1e5373..ffa31b9170 100644 --- a/nse_solver/make_table/burn_cell.H +++ b/nse_solver/make_table/burn_cell.H @@ -31,32 +31,32 @@ void burn_cell_c() state.rho = rho; state.y_e = Ye; - if (state.y_e > 0.52_rt){ - state.mu_p = -1.0_rt; - state.mu_n = -16.0_rt; - } - else if (state.y_e > 0.48_rt){ - state.mu_p = -6.0_rt; - state.mu_n = -11.0_rt; - } - else if (state.y_e > 0.4_rt){ - state.mu_p = -10.0_rt; - state.mu_n = -7.0_rt; - } - else{ - state.mu_p = -18.0; - state.mu_n = -1.0; - } + if (state.y_e > 0.52_rt){ + state.mu_p = -1.0_rt; + state.mu_n = -16.0_rt; + } + else if (state.y_e > 0.48_rt){ + state.mu_p = -6.0_rt; + state.mu_n = -11.0_rt; + } + else if (state.y_e > 0.4_rt){ + state.mu_p = -10.0_rt; + state.mu_n = -7.0_rt; + } + else{ + state.mu_p = -18.0; + state.mu_n = -1.0; + } // find the nse state const bool assume_ye_is_valid = true; Real eps = 1.e-10; - use_hybrid_solver = 1; - + use_hybrid_solver = 1; + auto nse_state = get_actual_nse_state(state, eps, assume_ye_is_valid); - std::cout << std::scientific; + std::cout << std::scientific; std::cout << std::setw(20) << state.rho << " " << std::setw(20) << state.T << " " << std::fixed << std::setw(20) << state.y_e << " " diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H index 1957dfe9e3..dba9fb88a3 100644 --- a/nse_solver/nse_check.H +++ b/nse_solver/nse_check.H @@ -81,55 +81,55 @@ void find_max_nucs(int& max_nucs, const int current_nuc_ind) { for (int n = NumSpec-1; n > NSE_INDEX::he4_index; --n) { if ( - (aion[n] == aion[current_nuc_ind]-1 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]-1)) - || - (aion[n] == aion[current_nuc_ind]+3 && - (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]+4 && zion[n] == zion[current_nuc_ind]+2) - ) { + (aion[n] == aion[current_nuc_ind]-1 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]-1)) + || + (aion[n] == aion[current_nuc_ind]+3 && + (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]+4 && zion[n] == zion[current_nuc_ind]+2) + ) { ++max_nucs; } if ( - (aion[n] == aion[current_nuc_ind]+3 && - (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]+2 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 - || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]-2 && - (zion[n] == zion[current_nuc_ind]-2 || zion[n] == zion[current_nuc_ind]-1 - || zion[n] == zion[current_nuc_ind])) - ) { + (aion[n] == aion[current_nuc_ind]+3 && + (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]+2 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 + || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]-2 && + (zion[n] == zion[current_nuc_ind]-2 || zion[n] == zion[current_nuc_ind]-1 + || zion[n] == zion[current_nuc_ind])) + ) { ++max_nucs; } if ( - (aion[n] == aion[current_nuc_ind]+2 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 - || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]+1 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) - || - (aion[n] == aion[current_nuc_ind]-3 && - (zion[n] == zion[current_nuc_ind]-1 || zion[n] == zion[current_nuc_ind]-2)) - ) { + (aion[n] == aion[current_nuc_ind]+2 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 + || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]+1 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) + || + (aion[n] == aion[current_nuc_ind]-3 && + (zion[n] == zion[current_nuc_ind]-1 || zion[n] == zion[current_nuc_ind]-2)) + ) { ++max_nucs; } if ( - (aion[n] == aion[current_nuc_ind]+1 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) - || - (aion[n] == aion[current_nuc_ind]-4 && zion[n] == zion[current_nuc_ind]-2) - ) { + (aion[n] == aion[current_nuc_ind]+1 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) + || + (aion[n] == aion[current_nuc_ind]-4 && zion[n] == zion[current_nuc_ind]-2) + ) { ++max_nucs; } @@ -139,7 +139,7 @@ void find_max_nucs(int& max_nucs, const int current_nuc_ind) { AMREX_GPU_HOST_DEVICE AMREX_INLINE void get_rate_by_nuc(int& rate_index, amrex::Array1D reactants_indices, - amrex::Array1D products_indices) { + amrex::Array1D products_indices) { rate_index = -1; bool found_index = false; @@ -155,11 +155,11 @@ void get_rate_by_nuc(int& rate_index, amrex::Array1D reactants_indice for (int i = 1; i <= Rates::NumRates; ++i) { for (int j = 1; j <= 3; ++j) { if (NSE_INDEX::rate_indices(i, j) == reactants_indices(j) && NSE_INDEX::rate_indices(i, j+3) == products_indices(j)) { - found_index = true; + found_index = true; } else { - found_index = false; - break; + found_index = false; + break; } } @@ -174,8 +174,8 @@ void get_rate_by_nuc(int& rate_index, amrex::Array1D reactants_indice AMREX_GPU_HOST_DEVICE AMREX_INLINE void find_fast_reaction_cycle(const int current_nuc_ind, const amrex::Array1D& Y, - const amrex::Array1D& screened_rates, - const burn_t& state, const amrex::Real& t_s, bool& found_fast_reaction_cycle) { + const amrex::Array1D& screened_rates, + const burn_t& state, const amrex::Real& t_s, bool& found_fast_reaction_cycle) { // This function finds the fast reaction cycle of a given nuclei. // This is done to satisfy Section 4.4.3 of Kushnir @@ -226,61 +226,61 @@ void find_fast_reaction_cycle(const int current_nuc_ind, const amrex::Array1D NSE_INDEX::he4_index; --n) { if ( - (aion[n] == aion[current_nuc_ind]-1 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]-1)) - || - (aion[n] == aion[current_nuc_ind]+3 && - (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]+4 && zion[n] == zion[current_nuc_ind]+2) - ) { + (aion[n] == aion[current_nuc_ind]-1 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]-1)) + || + (aion[n] == aion[current_nuc_ind]+3 && + (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]+4 && zion[n] == zion[current_nuc_ind]+2) + ) { reactions(2, inter_nuc_ind1) = n; ++inter_nuc_ind1; } if ( - (aion[n] == aion[current_nuc_ind]+3 && - (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]+2 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 - || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]-2 && - (zion[n] == zion[current_nuc_ind]-2 || zion[n] == zion[current_nuc_ind]-1 - || zion[n] == zion[current_nuc_ind])) - ) { + (aion[n] == aion[current_nuc_ind]+3 && + (zion[n] == zion[current_nuc_ind]+1 || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]+2 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 + || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]-2 && + (zion[n] == zion[current_nuc_ind]-2 || zion[n] == zion[current_nuc_ind]-1 + || zion[n] == zion[current_nuc_ind])) + ) { reactions(3, inter_nuc_ind2) = n; ++inter_nuc_ind2; } if ( - (aion[n] == aion[current_nuc_ind]+2 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 - || zion[n] == zion[current_nuc_ind]+2)) - || - (aion[n] == aion[current_nuc_ind]+1 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) - || - (aion[n] == aion[current_nuc_ind]-3 && - (zion[n] == zion[current_nuc_ind]-1 || zion[n] == zion[current_nuc_ind]-2)) - ) { + (aion[n] == aion[current_nuc_ind]+2 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1 + || zion[n] == zion[current_nuc_ind]+2)) + || + (aion[n] == aion[current_nuc_ind]+1 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) + || + (aion[n] == aion[current_nuc_ind]-3 && + (zion[n] == zion[current_nuc_ind]-1 || zion[n] == zion[current_nuc_ind]-2)) + ) { reactions(4, inter_nuc_ind3) = n; ++inter_nuc_ind3; } if ( - (aion[n] == aion[current_nuc_ind]+1 && - (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) - || - (aion[n] == aion[current_nuc_ind]-4 && zion[n] == zion[current_nuc_ind]-2) - ) { + (aion[n] == aion[current_nuc_ind]+1 && + (zion[n] == zion[current_nuc_ind] || zion[n] == zion[current_nuc_ind]+1)) + || + (aion[n] == aion[current_nuc_ind]-4 && zion[n] == zion[current_nuc_ind]-2) + ) { reactions(5, inter_nuc_ind4) = n; - ++inter_nuc_ind4; + ++inter_nuc_ind4; } } @@ -296,7 +296,7 @@ void find_fast_reaction_cycle(const int current_nuc_ind, const amrex::Array1D& group_ind) { + const amrex::Array1D& group_ind) { // This function returns the root index of the nuclei // by providing the nuclei index [0, NumSpec-1], and group indices, group_ind @@ -698,23 +698,23 @@ bool in_single_group(const amrex::Array1D& group_ind) { if (zion[n] >= 14) { - // Get even N group index - if (evenN_group == -1 && std::fmod(aion[n] - zion[n], 2) == 0.0_rt) { - evenN_group = get_root_index(n, group_ind); - continue; - } - - // Get odd N group index - if (oddN_group == -1 && std::fmod(aion[n] - zion[n], 2) == 1.0_rt) { - oddN_group = get_root_index(n, group_ind); - continue; - } - - if ((std::fmod(aion[n] - zion[n], 2) == 0.0_rt && evenN_group != get_root_index(n, group_ind)) || - (std::fmod(aion[n] - zion[n], 2) == 1.0_rt && oddN_group != get_root_index(n, group_ind))) { - in_single_group = false; - break; - } + // Get even N group index + if (evenN_group == -1 && std::fmod(aion[n] - zion[n], 2) == 0.0_rt) { + evenN_group = get_root_index(n, group_ind); + continue; + } + + // Get odd N group index + if (oddN_group == -1 && std::fmod(aion[n] - zion[n], 2) == 1.0_rt) { + oddN_group = get_root_index(n, group_ind); + continue; + } + + if ((std::fmod(aion[n] - zion[n], 2) == 0.0_rt && evenN_group != get_root_index(n, group_ind)) || + (std::fmod(aion[n] - zion[n], 2) == 1.0_rt && oddN_group != get_root_index(n, group_ind))) { + in_single_group = false; + break; + } } } @@ -726,11 +726,11 @@ bool in_single_group(const amrex::Array1D& group_ind) { AMREX_GPU_HOST_DEVICE AMREX_INLINE void fill_merge_indices(amrex::Array1D& merge_indices, - amrex::Real& fastest_t, const int current_rate, - const amrex::Array1D& Y, const burn_t& state, - const amrex::Array1D& screened_rates, - const amrex::Array1D& group_ind, - const amrex::Real& t_s) { + amrex::Real& fastest_t, const int current_rate, + const amrex::Array1D& Y, const burn_t& state, + const amrex::Array1D& screened_rates, + const amrex::Array1D& group_ind, + const amrex::Real& t_s) { // This function fills in the merge index of the current rate. // The timescale of the current rate must be shorter (faster) than fastest_t @@ -758,15 +758,15 @@ void fill_merge_indices(amrex::Array1D& merge_indices, } if (NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::h1_index - && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::n_index - && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::he4_index - ) { + && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::n_index + && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::he4_index + ) { ++m; // skip if there are more than 2 non neutron, proton, or helium-4 in the rate if (m > 2) { - return; + return; } nonNPA_ind(m) = k; } @@ -793,7 +793,7 @@ void fill_merge_indices(amrex::Array1D& merge_indices, // return if nonLIG_root index is repeated, meaning these isotope already merged if (root_index == nonLIG_root) { - return; + return; } nonLIG_root = root_index; } @@ -810,7 +810,7 @@ void fill_merge_indices(amrex::Array1D& merge_indices, // first 3 for reactants, last 3 for products amrex::Array1D Y_group = {0.0_rt, 0.0_rt, 0.0_rt, - 0.0_rt, 0.0_rt, 0.0_rt}; + 0.0_rt, 0.0_rt, 0.0_rt}; for (int k = 2; k <= 6; ++k) { if (NSE_INDEX::rate_indices(current_rate, k) == -1) { @@ -821,16 +821,16 @@ void fill_merge_indices(amrex::Array1D& merge_indices, // let LIG group use their own Y instead of Y_group. (not sure if this is true) if (root_index == get_root_index(NSE_INDEX::he4_index, group_ind)) { - Y_group(k) = Y(NSE_INDEX::rate_indices(current_rate, k) + 1); + Y_group(k) = Y(NSE_INDEX::rate_indices(current_rate, k) + 1); } else { for (int n = 0; n < NumSpec; ++n) { - + // if nuclei have the same root_index, then they are the same group - if (root_index == get_root_index(n, group_ind)) { - Y_group(k) += Y(n + 1); - } + if (root_index == get_root_index(n, group_ind)) { + Y_group(k) += Y(n + 1); + } } } } @@ -886,10 +886,10 @@ void fill_merge_indices(amrex::Array1D& merge_indices, // If nonNPA index is -1, meaning null, then use LIG index if (nonNPA_ind(n) == -1) { - merge_indices(n) = get_root_index(NSE_INDEX::he4_index, group_ind); + merge_indices(n) = get_root_index(NSE_INDEX::he4_index, group_ind); } else { - merge_indices(n) = NSE_INDEX::rate_indices(current_rate, nonNPA_ind(n)); + merge_indices(n) = NSE_INDEX::rate_indices(current_rate, nonNPA_ind(n)); } } } @@ -899,9 +899,9 @@ void fill_merge_indices(amrex::Array1D& merge_indices, AMREX_GPU_HOST_DEVICE AMREX_INLINE void nse_grouping(amrex::Array1D& group_ind, const burn_t& state, - const amrex::Array1D& Y, - const amrex::Array1D& screened_rates, - const amrex::Real& t_s) { + const amrex::Array1D& Y, + const amrex::Array1D& screened_rates, + const amrex::Real& t_s) { // This function groups all the nuclei using group_ind // which contains the node # @@ -914,7 +914,7 @@ void nse_grouping(amrex::Array1D& group_ind, const burn_t& stat // let n,p,a form the same group (LIG) initially, let 1 be index of LIG if (i == NSE_INDEX::h1_index + 1 || i == NSE_INDEX::n_index + 1 - || i == NSE_INDEX::he4_index + 1) { + || i == NSE_INDEX::he4_index + 1) { // group_ind(i) = NSE_INDEX::he4_index; @@ -1094,7 +1094,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { for (int n = NumSpec-1; n >= 0; --n) { if (found_fast_reaction_cycle) { - break; + break; } find_fast_reaction_cycle(n, Y, rate_eval.screened_rates, state, t_s, found_fast_reaction_cycle); diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H index 39028c0603..efe699144f 100644 --- a/nse_solver/nse_solver.H +++ b/nse_solver/nse_solver.H @@ -66,10 +66,10 @@ T get_nonexponent_nse_state(const T& state) { // find nse mass frac without the exponent term. Real power = 2.0 * M_PI * network::mion(n+1) * - C::k_B * T_in / (C::hplanck * C::hplanck); + C::k_B * T_in / (C::hplanck * C::hplanck); nse_state.xn[n] = network::mion(n+1) * pf * spin / state.rho * - std::sqrt(power * power * power); + std::sqrt(power * power * power); } return nse_state; @@ -133,9 +133,9 @@ void apply_nse_exponent(T& nse_state) { // iteration will make it happy again exponent = amrex::min(500.0_rt, - (zion[n] * nse_state.mu_p + (aion[n] - zion[n]) * - nse_state.mu_n - u_c + network::bion(n+1)) / - C::k_B / T_in * C::Legacy::MeV2erg); + (zion[n] * nse_state.mu_p + (aion[n] - zion[n]) * + nse_state.mu_n - u_c + network::bion(n+1)) / + C::k_B / T_in * C::Legacy::MeV2erg); nse_state.xn[n] *= std::exp(exponent); } @@ -274,42 +274,42 @@ void nse_hybrid_solver(T& state, amrex::Real eps=1.0e-10_rt) { for (int j = 0; j < 20; ++j) { - hj.x(1) = inner_x(1); - hj.x(2) = inner_x(2); + hj.x(1) = inner_x(1); + hj.x(2) = inner_x(2); - // hybrj<2, T>(hj, state, fcn_hybrid, jcn_hybrid); - hybrj(hj, state); + // hybrj<2, T>(hj, state, fcn_hybrid, jcn_hybrid); + hybrj(hj, state); - fcn(hj.x, f, state, flag); + fcn(hj.x, f, state, flag); - if (std::abs(f(1)) < eps && std::abs(f(2)) < eps) { + if (std::abs(f(1)) < eps && std::abs(f(2)) < eps) { - state.mu_p = hj.x(1); - state.mu_n = hj.x(2); - return; - } + state.mu_p = hj.x(1); + state.mu_n = hj.x(2); + return; + } - if (f(1) > 0.0_rt && f(2) > 0.0_rt) { - is_pos_new = true; - } - else{ - is_pos_new = false; - } + if (f(1) > 0.0_rt && f(2) > 0.0_rt) { + is_pos_new = true; + } + else{ + is_pos_new = false; + } - if (is_pos_old != is_pos_new) { - dx *= 0.8_rt; - } + if (is_pos_old != is_pos_new) { + dx *= 0.8_rt; + } - if (is_pos_new) { - inner_x(1) -= dx; - inner_x(2) -= dx; - } - else{ - inner_x(1) += dx; - inner_x(2) += dx; - } + if (is_pos_new) { + inner_x(1) -= dx; + inner_x(2) -= dx; + } + else{ + inner_x(1) += dx; + inner_x(2) += dx; + } - is_pos_old = is_pos_new; + is_pos_old = is_pos_new; } outer_x(1) -= 0.5_rt; @@ -422,7 +422,7 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { // check whether solution results in nan if (std::isnan(x(1)) || std::isnan(x(2))) { - amrex::Error("Nan encountered, likely due to overflow in digits or not making good progress"); + amrex::Error("Nan encountered, likely due to overflow in digits or not making good progress"); } // update constraint @@ -488,11 +488,11 @@ T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt, for (int n = 0; n < NumSpec; ++n) { #ifdef NEW_NETWORK_IMPLEMENTATION if (n == NSE_INDEX::h1_index) { - continue; + continue; } #endif if (zion[n] != aion[n] - zion[n]) { - singular_network = false; + singular_network = false; } } diff --git a/screening/screen.H b/screening/screen.H index 898772c8d8..6ebe21717b 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -253,8 +253,8 @@ void actual_screen5 (const plasma_state_t& state, alph12dt = 0.0_rt; } - // redetermine previous factors if 3*gamma_ij/tau_ij > 1.6 - + // redetermine previous factors if 3*gamma_ij/tau_ij > 1.6 + gamef = 1.6e0_rt * tau12; if constexpr (do_T_derivatives) { gamefdt = 1.6e0_rt * tau12dt; @@ -292,8 +292,8 @@ void actual_screen5 (const plasma_state_t& state, Real gamp14 = std::pow(gamp, 0.25_rt); Real rr = 1.0_rt/gamp; - // Here we follow Eq. A9 in Wallace:1982 - // See Eq. 25 Alastuey:1978, Eq. 16 and 17 in Jancovici:1977 for reference + // Here we follow Eq. A9 in Wallace:1982 + // See Eq. 25 Alastuey:1978, Eq. 16 and 17 in Jancovici:1977 for reference Real cc = 0.896434e0_rt * gamp * scn_fac.zhat - 3.44740e0_rt * gamp14 * scn_fac.zhat2 - 0.5551e0_rt * (std::log(gamp) + scn_fac.lzav) @@ -308,14 +308,14 @@ void actual_screen5 (const plasma_state_t& state, - 0.5551e0_rt *rr * gampdt; } - // (3gamma_ij/tau_ij)^3 + // (3gamma_ij/tau_ij)^3 Real a3 = alph12 * alph12 * alph12; Real da3 = 3.0e0_rt * alph12 * alph12; - // Part of Eq. 28 in Alastuey:1978 + // Part of Eq. 28 in Alastuey:1978 qq = 0.014e0_rt + 0.0128e0_rt*alph12; - // Part of Eq. 28 in Alastuey:1978 + // Part of Eq. 28 in Alastuey:1978 rr = (5.0_rt/32.0_rt) - alph12*qq; [[maybe_unused]] Real drrdt; if constexpr (do_T_derivatives) { @@ -323,20 +323,20 @@ void actual_screen5 (const plasma_state_t& state, drrdt = -(alph12dt*qq + alph12*dqqdt); } - // Part of Eq. 28 in Alastuey:1978 + // Part of Eq. 28 in Alastuey:1978 Real ss = tau12*rr; - // Part of Eq. 31 in Alastuey:1978 + // Part of Eq. 31 in Alastuey:1978 Real tt = -0.0098e0_rt + 0.0048e0_rt*alph12; - // Part of Eq. 31 in Alastuey:1978 + // Part of Eq. 31 in Alastuey:1978 Real uu = 0.0055e0_rt + alph12*tt; - // Part of Eq. 31 in Alastuey:1978 + // Part of Eq. 31 in Alastuey:1978 Real vv = gamef * alph12 * uu; - // Exponent of Eq. 32 in Alastuey:1978, which uses Eq.28 and Eq.31 - // Strong screening factor + // Exponent of Eq. 32 in Alastuey:1978, which uses Eq.28 and Eq.31 + // Strong screening factor h12 = cc - a3 * (ss + vv); if constexpr (do_T_derivatives) { Real dssdt = tau12dt*rr + tau12*drrdt; @@ -347,8 +347,8 @@ void actual_screen5 (const plasma_state_t& state, dh12dt = dccdt - rr*alph12dt - a3*(dssdt + dvvdt); } - // See conclusion and Eq. 34 in Alastuey:1978 - // This is an extra factor to account for quantum effects + // See conclusion and Eq. 34 in Alastuey:1978 + // This is an extra factor to account for quantum effects rr = 1.0_rt - 0.0562e0_rt*a3; if constexpr (do_T_derivatives) { ss = -0.0562e0_rt*da3; @@ -358,7 +358,7 @@ void actual_screen5 (const plasma_state_t& state, Real xlgfac; [[maybe_unused]] Real dxlgfacdt; - // In extreme case, rr is 0.77, see conclusion in Alastuey:1978 + // In extreme case, rr is 0.77, see conclusion in Alastuey:1978 if (rr >= 0.77e0_rt) { xlgfac = rr; if constexpr (do_T_derivatives) { @@ -371,15 +371,15 @@ void actual_screen5 (const plasma_state_t& state, } } - // Include the extra factor that accounts for quantum effects + // Include the extra factor that accounts for quantum effects h12 = std::log(xlgfac) + h12; rr = 1.0_rt/xlgfac; if constexpr (do_T_derivatives) { dh12dt = rr*dxlgfacdt + dh12dt; } - // If gamma_ij < upper limit of intermediate regime - // then it is in the intermediate regime, else strong screening. + // If gamma_ij < upper limit of intermediate regime + // then it is in the intermediate regime, else strong screening. if (gamef <= gamefs) { Real dgamma = 1.0e0_rt/(gamefs - gamefx); @@ -389,8 +389,8 @@ void actual_screen5 (const plasma_state_t& state, vv = h12; - // Then the screening factor is a combination - // of the strong and weak screening factor. + // Then the screening factor is a combination + // of the strong and weak screening factor. h12 = h12w*rr + vv*ss; if constexpr (do_T_derivatives) { drrdt = -dgamma*gamefdt; @@ -890,21 +890,21 @@ void chabrier1998 (const plasma_state_t& state, //see Calder 2007 Eq.A8 and Alastuey1978, Eq. 24 and 31 Real quantum_corr_1 = -tau12 * (5.0_rt/32.0_rt * std::pow(b_fac, 3) - 0.014_rt * std::pow(b_fac, 4) - - 0.128_rt * std::pow(b_fac, 5)); + - 0.128_rt * std::pow(b_fac, 5)); Real quantum_corr_2 = -Gamma_eff * (0.0055_rt * std::pow(b_fac, 4) - 0.0098_rt * std::pow(b_fac, 5) - + 0.0048_rt * std::pow(b_fac, 6)); + + 0.0048_rt * std::pow(b_fac, 6)); [[maybe_unused]] Real quantum_corr_1_dT; [[maybe_unused]] Real quantum_corr_2_dT; if constexpr (do_T_derivatives) { quantum_corr_1_dT = tau12dT / tau12 * quantum_corr_1 - tau12 * b_fac_dT * (15.0_rt/32.0_rt * std::pow(b_fac, 2) - 0.014_rt * 4.0_rt * std::pow(b_fac, 3) - - 0.128_rt * 5.0_rt * std::pow(b_fac, 4)); + - 0.128_rt * 5.0_rt * std::pow(b_fac, 4)); quantum_corr_2_dT = Gamma_eff_dT / Gamma_eff * quantum_corr_2 - Gamma_eff * b_fac_dT * (0.0055_rt * 4.0_rt * std::pow(b_fac, 3) - 0.0098_rt * 5.0_rt - * std::pow(b_fac, 4) + 0.0048_rt * 6.0_rt * std::pow(b_fac, 5)); + * std::pow(b_fac, 4) + 0.0048_rt * 6.0_rt * std::pow(b_fac, 5)); } // See Calder2007 Appendix Eq. A8. diff --git a/unit_test/burn_cell/burn_cell.H b/unit_test/burn_cell/burn_cell.H index aeb5d826dc..21f1f0046f 100644 --- a/unit_test/burn_cell/burn_cell.H +++ b/unit_test/burn_cell/burn_cell.H @@ -122,7 +122,7 @@ void burn_cell_c() state_over_time << std::setw(25) << "Temperature"; for(int x = 0; x < NumSpec; ++x){ const std::string& element = short_spec_names_cxx[x]; - state_over_time << std::setw(25) << element; + state_over_time << std::setw(25) << element; } state_over_time << std::endl; state_over_time << std::setprecision(15); @@ -153,7 +153,7 @@ void burn_cell_c() Real tend = std::pow(10.0_rt, std::log10(tfirst) + dlogt * n); Real dt = tend - t; - burner(burn_state, dt); + burner(burn_state, dt); if (! burn_state.success) { amrex::Error("integration failed"); @@ -178,12 +178,12 @@ void burn_cell_c() t += dt; - state_over_time << std::setw(25) << t; - state_over_time << std::setw(25) << burn_state.T; - for (double X : burn_state.xn) { - state_over_time << std::setw(25) << X; - } - state_over_time << std::endl; + state_over_time << std::setw(25) << t; + state_over_time << std::setw(25) << burn_state.T; + for (double X : burn_state.xn) { + state_over_time << std::setw(25) << X; + } + state_over_time << std::endl; } state_over_time.close(); diff --git a/unit_test/burn_cell_sdc/burn_cell.H b/unit_test/burn_cell_sdc/burn_cell.H index c83752ecfa..6567ad8df3 100644 --- a/unit_test/burn_cell_sdc/burn_cell.H +++ b/unit_test/burn_cell_sdc/burn_cell.H @@ -324,8 +324,8 @@ void burn_cell_c() state_over_time << std::setw(12) << "Density"; state_over_time << std::setw(12) << "Temperature"; for(int x = 0; x < NumSpec; ++x){ - const std::string& element = short_spec_names_cxx[x]; - state_over_time << std::setw(12) << element; + const std::string& element = short_spec_names_cxx[x]; + state_over_time << std::setw(12) << element; } state_over_time << std::endl; @@ -355,19 +355,19 @@ void burn_cell_c() std::cout << "burning for dt = " << dt << std::endl; - burner(burn_state, dt); + burner(burn_state, dt); nstep_int += burn_state.n_step; t += dt; - state_over_time << std::setw(10) << t; - state_over_time << std::setw(12) << burn_state.rho; - state_over_time << std::setw(12) << burn_state.T; - for (int x = 0; x < NumSpec; ++x){ - state_over_time << std::setw(12) << burn_state.y[SFS+x] / burn_state.rho; - } - state_over_time << std::endl; + state_over_time << std::setw(10) << t; + state_over_time << std::setw(12) << burn_state.rho; + state_over_time << std::setw(12) << burn_state.T; + for (int x = 0; x < NumSpec; ++x){ + state_over_time << std::setw(12) << burn_state.y[SFS+x] / burn_state.rho; + } + state_over_time << std::endl; } state_over_time.close(); diff --git a/unit_test/test_ase/make_table/burn_cell.H b/unit_test/test_ase/make_table/burn_cell.H index b13295a86d..b9eb41590c 100644 --- a/unit_test/test_ase/make_table/burn_cell.H +++ b/unit_test/test_ase/make_table/burn_cell.H @@ -35,44 +35,44 @@ void burn_cell_c() state.rho = rho; state.y_e = Ye; - if (state.y_e > 0.52_rt){ - state.mu_p = -1.0_rt; - state.mu_n = -16.0_rt; - } - else if (state.y_e > 0.48_rt){ - state.mu_p = -6.0_rt; - state.mu_n = -11.0_rt; - } - else if (state.y_e > 0.4_rt){ - state.mu_p = -10.0_rt; - state.mu_n = -7.0_rt; - } - else{ - state.mu_p = -18.0; - state.mu_n = -1.0; - } + if (state.y_e > 0.52_rt){ + state.mu_p = -1.0_rt; + state.mu_n = -16.0_rt; + } + else if (state.y_e > 0.48_rt){ + state.mu_p = -6.0_rt; + state.mu_n = -11.0_rt; + } + else if (state.y_e > 0.4_rt){ + state.mu_p = -10.0_rt; + state.mu_n = -7.0_rt; + } + else{ + state.mu_p = -18.0; + state.mu_n = -1.0; + } // find the nse state const bool assume_ye_is_valid = true; Real eps = 1.e-10; - use_hybrid_solver = 1; - + use_hybrid_solver = 1; + auto nse_state = get_actual_nse_state(state, eps, assume_ye_is_valid); - for (int i = 0; i < NumSpec; ++i){ - state.xn[i] = nse_state.xn[i]; - } - - bool in_nse_state = in_nse(state); - - std::cout << std::scientific; + for (int i = 0; i < NumSpec; ++i){ + state.xn[i] = nse_state.xn[i]; + } + + bool in_nse_state = in_nse(state); + + std::cout << std::scientific; std::cout << std::setw(20) << state.rho << " " << std::setw(20) << state.T << " " << std::fixed << std::setw(20) << state.y_e << " " << std::setw(20) << state.mu_p << " " << std::setw(20) << state.mu_n << " " - << std::setw(20) << in_nse_state << std::endl; + << std::setw(20) << in_nse_state << std::endl; } } diff --git a/unit_test/test_react/main.cpp b/unit_test/test_react/main.cpp index fde5631dbe..c8788312c3 100644 --- a/unit_test/test_react/main.cpp +++ b/unit_test/test_react/main.cpp @@ -52,8 +52,8 @@ void main_main () // number of cells on each side of a square (or cubic) domain. pp.get("n_cell", n_cell); - print_every_nrhs = 0; - pp.query("print_every_nrhs", print_every_nrhs); + print_every_nrhs = 0; + pp.query("print_every_nrhs", print_every_nrhs); // The domain is broken into boxes of size max_grid_size max_grid_size = 32; diff --git a/util/hybrj/hybrj.H b/util/hybrj/hybrj.H index ce5353b977..883c258a73 100644 --- a/util/hybrj/hybrj.H +++ b/util/hybrj/hybrj.H @@ -188,18 +188,18 @@ void hybrj(hybrj_t& hj, // store the direction p and x + p. calculate the norm of p. - bool invalid = false; + bool invalid = false; - for (int i = 1; i < neqs; ++i) { - if ( amrex::isinf(hj.wa1(i)) || amrex::isnan(hj.wa1(i)) ){ - invalid = true; - } - } + for (int i = 1; i < neqs; ++i) { + if ( amrex::isinf(hj.wa1(i)) || amrex::isnan(hj.wa1(i)) ){ + invalid = true; + } + } - if (invalid) { - finished = true; - break; - } + if (invalid) { + finished = true; + break; + } for (int j = 1; j <= neqs; ++j) { hj.wa1(j) = -hj.wa1(j); @@ -344,9 +344,9 @@ void hybrj(hybrj_t& hj, // calculate the rank one modification to the jacobian // and update qtf if necessary. - if (pnorm == 0.0_rt) { - break; - } + if (pnorm == 0.0_rt) { + break; + } for (int j = 1; j <= hj.n; ++j) { Real sum = 0.0_rt; diff --git a/util/hybrj/hybrj_dogleg.H b/util/hybrj/hybrj_dogleg.H index 0bc5822b13..77a6c0d8f4 100644 --- a/util/hybrj/hybrj_dogleg.H +++ b/util/hybrj/hybrj_dogleg.H @@ -105,15 +105,15 @@ void dogleg(amrex::Array1D& r, wa1(i) += r(l) * temp; l++; - if ( amrex::isnan(wa1(i)) || amrex::isinf(wa1(i)) ) { - return; - } + if ( amrex::isnan(wa1(i)) || amrex::isinf(wa1(i)) ) { + return; + } } wa1(j) /= diag(j); - if ( amrex::isnan(wa1(j)) || amrex::isinf(wa1(j)) ) { - return; - } + if ( amrex::isnan(wa1(j)) || amrex::isinf(wa1(j)) ) { + return; + } } // calculate the norm of the scaled gradient and test for diff --git a/util/hybrj/hybrj_enorm.H b/util/hybrj/hybrj_enorm.H index 176b80a774..fa24cfc61b 100644 --- a/util/hybrj/hybrj_enorm.H +++ b/util/hybrj/hybrj_enorm.H @@ -43,9 +43,9 @@ Real enorm(const int n, amrex::Array1D& x) { for (int i = 1; i <= n; ++i) { Real xabs = std::abs(x(i)); - if ( amrex::isnan(xabs) || amrex::isinf(xabs) ) { - return 0.0_rt; - } + if ( amrex::isnan(xabs) || amrex::isinf(xabs) ) { + return 0.0_rt; + } if (xabs <= rdwarf || xabs >= agiant) { if (xabs > rdwarf) { diff --git a/util/hybrj/hybrj_r1mpyq.H b/util/hybrj/hybrj_r1mpyq.H index b31bcd1f46..1e1cade991 100644 --- a/util/hybrj/hybrj_r1mpyq.H +++ b/util/hybrj/hybrj_r1mpyq.H @@ -51,7 +51,7 @@ void r1mpyq(amrex::Array2D& a, fcos = 1.0_rt / v(j); fsin = std::sqrt(1.0_rt - fcos*fcos); } - else { + else { // if (std::abs(v(j)) <= 1.0_rt) { fsin = v(j); fcos = std::sqrt(1.0_rt - fsin*fsin); @@ -70,7 +70,7 @@ void r1mpyq(amrex::Array2D& a, fcos = 1.0_rt / w(j); fsin = std::sqrt(1.0_rt - fcos*fcos); } - else { + else { // if (std::abs(w(j)) <= 1.0_rt) { fsin = w(j); fcos = std::sqrt(1.0_rt - fsin*fsin); @@ -121,7 +121,7 @@ void r1mpyq(amrex::Array1D& a, fcos = 1.0_rt / v(j); fsin = std::sqrt(1.0_rt - fcos*fcos); } - else { + else { // if (std::abs(v(j)) <= 1.0_rt) { fsin = v(j); fcos = std::sqrt(1.0_rt - fsin*fsin); @@ -138,7 +138,7 @@ void r1mpyq(amrex::Array1D& a, fcos = 1.0_rt / w(j); fsin = std::sqrt(1.0_rt - fcos*fcos); } - else { + else { // if (std::abs(w(j)) <= 1.0_rt) { fsin = w(j); fcos = std::sqrt(1.0_rt - fsin*fsin); diff --git a/util/microphysics_sort.H b/util/microphysics_sort.H index 3949917a29..63403b7c2c 100644 --- a/util/microphysics_sort.H +++ b/util/microphysics_sort.H @@ -12,9 +12,9 @@ void sort_Array1D(amrex::Array1D& array_1d){ for (int i = l; i <= m; ++i){ for (int j = l; j <= m-i; ++j){ if (array_1d(j) > array_1d(j+1)){ - scratch = array_1d(j); - array_1d(j) = array_1d(j+1); - array_1d(j+1) = scratch; + scratch = array_1d(j); + array_1d(j) = array_1d(j+1); + array_1d(j+1) = scratch; } } } From 1f779a922a860a62a11a5511aad6653bd2c3abb8 Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Fri, 15 Sep 2023 16:42:49 -0400 Subject: [PATCH 2/4] Trim trailing whitespace --- EOS/breakout/actual_eos.H | 2 +- EOS/eos_composition.H | 2 +- EOS/helmholtz/actual_eos.H | 8 +- EOS/polytrope/actual_eos.H | 4 +- EOS/ztwd/actual_eos.H | 8 +- conductivity/stellar/actual_conductivity.H | 2 +- constants/fundamental_constants.H | 4 +- integration/VODE/actual_integrator.H | 2 +- integration/VODE/vode_dvstep.H | 2 +- integration/integrator_type.H | 2 +- interfaces/burn_type.H | 4 +- interfaces/eos_type.H | 2 +- interfaces/tfactors.H | 30 +- networks/aprox13/actual_network.H | 2 +- networks/aprox19/actual_network.H | 2 +- networks/aprox21/actual_network.H | 2 +- networks/rhs.H | 2 +- .../triple_alpha_plus_cago/actual_network.H | 2 +- nse_solver/nse_check.H | 48 +- nse_solver/nse_solver.H | 14 +- rates/aprox_rates.H | 436 +++++++++--------- rates/aprox_rates_data.H | 2 +- screening/screen.H | 54 +-- .../burn_cell_primordial_chem/burn_cell.H | 4 +- unit_test/burn_cell_sdc/burn_cell.H | 12 +- .../test_aprox_rates/aprox_rates_util.cpp | 4 +- unit_test/test_aprox_rates/variables.cpp | 2 +- unit_test/test_ase/burn_cell.H | 10 +- unit_test/test_ase/main.cpp | 2 +- unit_test/test_ase/make_table/burn_cell.H | 2 +- unit_test/test_nse/nse_example.H | 18 +- unit_test/test_screening/screening_util.cpp | 2 +- unit_test/test_sdc/react_zones.H | 2 +- unit_test/test_sdc/variables.cpp | 2 +- unit_test/write_job_info.cpp | 2 +- 35 files changed, 349 insertions(+), 349 deletions(-) diff --git a/EOS/breakout/actual_eos.H b/EOS/breakout/actual_eos.H index 342b36fba7..25067019c8 100644 --- a/EOS/breakout/actual_eos.H +++ b/EOS/breakout/actual_eos.H @@ -130,7 +130,7 @@ void actual_eos (I input, T& state) } } - // Try to avoid the expensive log function. Since we don't need entropy + // Try to avoid the expensive log function. Since we don't need entropy // in hydro solver, set it to an invalid but "nice" value for the plotfile. if constexpr (has_entropy::value) { state.s = 1.0_rt; diff --git a/EOS/eos_composition.H b/EOS/eos_composition.H index cc72c5080d..be39fa0492 100644 --- a/EOS/eos_composition.H +++ b/EOS/eos_composition.H @@ -26,7 +26,7 @@ struct eos_xderivs_t { /// The auxiliary state provides an alternate description to the composition, /// in terms of Ye, abar, and binding energy / nucleon /// -template +template AMREX_GPU_HOST_DEVICE AMREX_INLINE void set_aux_comp_from_X(state_t& state) { diff --git a/EOS/helmholtz/actual_eos.H b/EOS/helmholtz/actual_eos.H index 51e7da22d9..38ed157b8e 100644 --- a/EOS/helmholtz/actual_eos.H +++ b/EOS/helmholtz/actual_eos.H @@ -29,13 +29,13 @@ using namespace eos_rp; // derivatives), number density of electrons and positron pair (along // with their derivatives), adiabatic indices, specific heats, and // relativistically correct sound speed are also returned. -// +// // this routine assumes planckian photons, an ideal gas of ions, // and an electron-positron gas with an arbitrary degree of relativity // and degeneracy. interpolation in a table of the helmholtz free energy // is used to return the electron-positron thermodynamic quantities. // all other derivatives are analytic. -// +// // references: cox & giuli chapter 24 ; timmes & swesty apj 1999 const std::string eos_name = "helmholtz"; @@ -510,7 +510,7 @@ void apply_ions (T& state) Real sion = (pion * deni + eion) * tempi + kergavo * ytot1 * y; Real dsiondd = (dpiondd * deni - pion * deni * deni + deiondd) * tempi - kergavo * deni * ytot1; - Real dsiondt = (dpiondt * deni + deiondt) * tempi - + Real dsiondt = (dpiondt * deni + deiondt) * tempi - (pion * deni + eion) * tempi * tempi + 1.5e0_rt * kergavo * tempi * ytot1; @@ -759,7 +759,7 @@ void apply_coulomb_corrections (T& state) decoulda = s * dpcoulda; decouldz = s * dpcouldz; - s = -avo_eos * kerg / (state.abar * plasg) * + s = -avo_eos * kerg / (state.abar * plasg) * (1.5e0_rt * c2 * x - a2 * (b2 - 1.0e0_rt) * y); dscouldd = s * plasgdd; dscouldt = s * plasgdt; diff --git a/EOS/polytrope/actual_eos.H b/EOS/polytrope/actual_eos.H index 8a2b353568..8de4617635 100644 --- a/EOS/polytrope/actual_eos.H +++ b/EOS/polytrope/actual_eos.H @@ -39,7 +39,7 @@ void actual_eos_init () // Available pre-defined polytrope options: // 1: Non-relativistic, fully degenerate electron gas - // 2: Relativistic, fully degenerate electron gas + // 2: Relativistic, fully degenerate electron gas if (polytrope_type > 0) { mu_e = polytrope_mu_e; @@ -88,7 +88,7 @@ bool is_input_valid (I input) //--------------------------------------------------------------------------- -// Public interfaces +// Public interfaces //--------------------------------------------------------------------------- inline diff --git a/EOS/ztwd/actual_eos.H b/EOS/ztwd/actual_eos.H index bc9a36e29b..6a02edb802 100644 --- a/EOS/ztwd/actual_eos.H +++ b/EOS/ztwd/actual_eos.H @@ -1,10 +1,10 @@ #ifndef ACTUAL_EOS_H #define ACTUAL_EOS_H -// This is the equation of state for zero-temperature white dwarf +// This is the equation of state for zero-temperature white dwarf // matter composed of degenerate electrons: // P = A * (x * (2x**2 - 3)(x**2 + 1)**1/2 + 3 sinh**-1(x)) -// +// // where rho = B x**3 and the constants are given by: // // A = pi m_e**4 c**5 / (3 h**3) = 6.0 x 10^22 dyne cm**-2 @@ -16,7 +16,7 @@ // h = (8A / B) (1 + x**2)**(1/2) // // The internal energy is calculated using the standard relation: -// +// // h = e + P / rho #include @@ -26,7 +26,7 @@ using namespace amrex; const std::string eos_name = "ztwd"; - + const Real A = M_PI * std::pow(C::m_e, 4) * std::pow(C::c_light, 5) / (3.0_rt * std::pow(C::hplanck, 3)); const Real B2 = 8.0_rt * M_PI * std::pow(C::m_e, 3) * std::pow(C::c_light, 3) * C::m_p / (3.0_rt * std::pow(C::hplanck, 3)); const Real iter_tol = 1.e-10_rt; diff --git a/conductivity/stellar/actual_conductivity.H b/conductivity/stellar/actual_conductivity.H index 2d2540e5a0..ee308c940e 100644 --- a/conductivity/stellar/actual_conductivity.H +++ b/conductivity/stellar/actual_conductivity.H @@ -41,7 +41,7 @@ actual_conductivity (T& state) // ocond = conductive contribution to the opacity (in cm**2/g) // opac = the total opacity (in cm**2/g) // conductivity = thermal conductivity (in erg/cm/K/sec) - // + // // various physical and derived constants // con2 = con1*sqrt(4*pi*e*e/me) diff --git a/constants/fundamental_constants.H b/constants/fundamental_constants.H index c07d0beca4..2a16b6e0e0 100644 --- a/constants/fundamental_constants.H +++ b/constants/fundamental_constants.H @@ -50,10 +50,10 @@ namespace C // mass of electron constexpr amrex::Real m_e = 9.10938291e-28; // g - + // atomic mass unit constexpr amrex::Real m_u = 1.6605390666e-24; // g - + // electron charge // NIST: q_e = 1.602176565e-19 C // diff --git a/integration/VODE/actual_integrator.H b/integration/VODE/actual_integrator.H index 45cdbe324c..35fa40c950 100644 --- a/integration/VODE/actual_integrator.H +++ b/integration/VODE/actual_integrator.H @@ -93,7 +93,7 @@ void actual_integrator (BurnT& state, Real dt) integrator_to_burn(vode_state, state); #ifdef NSE - // compute the temperature based on the energy release -- we need + // compute the temperature based on the energy release -- we need // this in case we failed in our burn here because we entered NSE #ifdef AUX_THERMO diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index e5602dd80e..6a77bc94e5 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -247,7 +247,7 @@ int dvstep (BurnT& state, DvodeT& vstate) (std::abs(vstate.y(i)) > vode_increase_change_factor * std::abs(y_save(i)) || std::abs(vstate.y(i)) < vode_decrease_change_factor * std::abs(y_save(i)))) { #ifdef MICROPHYSICS_DEBUG -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU std::cout << "rejecting step based on species " << i << " from " << y_save(i) << " to " << vstate.y(i) << std::endl; #endif #endif diff --git a/integration/integrator_type.H b/integration/integrator_type.H index 510f245da6..af653507dd 100644 --- a/integration/integrator_type.H +++ b/integration/integrator_type.H @@ -18,7 +18,7 @@ void update_density_in_time(const Real time, BurnT& state) // // we are always integrating from t = 0, so there is no offset // time needed here. The indexing of ydot_a is based on - // the indices in burn_t and is 0-based + // the indices in burn_t and is 0-based state.y[SRHO] = amrex::max(state.rho_orig + state.ydot_a[SRHO] * time, EOSData::mindens); // for consistency diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H index 9e7edbd4df..c489951626 100644 --- a/interfaces/burn_type.H +++ b/interfaces/burn_type.H @@ -249,7 +249,7 @@ void eos_to_burn (const T& eos_state, BurnT& burn_state) -// Given a burn type, copy the data relevant to the eos type. +// Given a burn type, copy the data relevant to the eos type. // Note that when doing simplified SDC integration, we should // avoid using this interface because the energy includes a // contribution from the advection term. However this is useful @@ -315,7 +315,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void normalize_abundances_sdc_burn (BurnT& state) { - // Constrain the partial densities in burn_t state to sum to the + // Constrain the partial densities in burn_t state to sum to the // density. // // This is meant to be used upon exit, and we assume that diff --git a/interfaces/eos_type.H b/interfaces/eos_type.H index cdfe210a13..db5d6dd695 100644 --- a/interfaces/eos_type.H +++ b/interfaces/eos_type.H @@ -109,7 +109,7 @@ struct eos_rep_t:eos_base_t { amrex::Real abar{}; amrex::Real zbar{}; - + #ifdef NSE_NET amrex::Real mu_p{}; //chemical potential of proton amrex::Real mu_n{}; //chemical potential of neutron diff --git a/interfaces/tfactors.H b/interfaces/tfactors.H index 856795b297..983b2d43c1 100644 --- a/interfaces/tfactors.H +++ b/interfaces/tfactors.H @@ -68,10 +68,10 @@ struct tf_t { amrex::Real lnt9; }; -AMREX_GPU_HOST_DEVICE inline +AMREX_GPU_HOST_DEVICE inline tf_t get_tfactors(amrex::Real temp) { - tf_t tf; + tf_t tf; tf.temp = temp; @@ -82,13 +82,13 @@ tf_t get_tfactors(amrex::Real temp) // tf.t95 = tf.t9*tf.t94; tf.t95 = tf.t92*tf.t93; // tf.t96 = tf.t9*tf.t95; - + tf.t912 = std::sqrt(tf.t9); tf.t932 = tf.t9*tf.t912; tf.t952 = tf.t9*tf.t932; // tf.t972 = tf.t9*tf.t952; tf.t972 = tf.t92*tf.t932; - + tf.t913 = std::cbrt(tf.t9); tf.t923 = tf.t913*tf.t913; tf.t943 = tf.t9*tf.t913; @@ -100,51 +100,51 @@ tf_t get_tfactors(amrex::Real temp) // tf.t934 = tf.t914*tf.t914*tf.t914; // tf.t954 = tf.t9*tf.t914; // tf.t974 = tf.t9*tf.t934; - + // tf.t915 = std::pow(tf.t9, 0.2_rt); // tf.t935 = tf.t915*tf.t915*tf.t915; // tf.t945 = tf.t915 * tf.t935; // tf.t965 = tf.t9 * tf.t915; - + // tf.t916 = std::pow(tf.t9, 1.0_rt/6.0_rt); // tf.t976 = tf.t9 * tf.t916; // tf.t9i76 = 1.0e0_rt/tf.t976; - + // tf.t917 = std::pow(tf.t9, 1.0_rt/7.0_rt); // tf.t927 = tf.t917*tf.t917; // tf.t947 = tf.t927*tf.t927; - + // tf.t918 = std::sqrt(tf.t914); // tf.t938 = tf.t918*tf.t918*tf.t918; // tf.t958 = tf.t938*tf.t918*tf.t918; - + tf.t9i = 1.0e0_rt/tf.t9; tf.t9i2 = tf.t9i*tf.t9i; // tf.t9i3 = tf.t9i2*tf.t9i; - + tf.t9i12 = 1.0e0_rt/tf.t912; tf.t9i32 = tf.t9i*tf.t9i12; // tf.t9i52 = tf.t9i*tf.t9i32; // tf.t9i72 = tf.t9i*tf.t9i52; - + tf.t9i13 = 1.0e0_rt/tf.t913; tf.t9i23 = tf.t9i13*tf.t9i13; tf.t9i43 = tf.t9i*tf.t9i13; tf.t9i53 = tf.t9i*tf.t9i23; - + // tf.t9i14 = 1.0e0_rt/tf.t914; // tf.t9i34 = tf.t9i14*tf.t9i14*tf.t9i14; // tf.t9i54 = tf.t9i*tf.t9i14; - + // tf.t9i15 = 1.0e0_rt/tf.t915; // tf.t9i35 = tf.t9i15*tf.t9i15*tf.t9i15; // tf.t9i45 = tf.t9i15 * tf.t9i35; // tf.t9i65 = tf.t9i*tf.t9i15; - + // tf.t9i17 = 1.0e0_rt/tf.t917; // tf.t9i27 = tf.t9i17*tf.t9i17; // tf.t9i47 = tf.t9i27*tf.t9i27; - + // tf.t9i18 = 1.0e0_rt/tf.t918; // tf.t9i38 = tf.t9i18*tf.t9i18*tf.t9i18; // tf.t9i58 = tf.t9i38*tf.t9i18*tf.t9i18; diff --git a/networks/aprox13/actual_network.H b/networks/aprox13/actual_network.H index a3687a71a6..7a1cce53f4 100644 --- a/networks/aprox13/actual_network.H +++ b/networks/aprox13/actual_network.H @@ -174,7 +174,7 @@ namespace Rates { } namespace RHS { - + AMREX_GPU_HOST_DEVICE AMREX_INLINE constexpr rhs_t rhs_data (int rate) { diff --git a/networks/aprox19/actual_network.H b/networks/aprox19/actual_network.H index 55ec7d06c6..75f8bfcf22 100644 --- a/networks/aprox19/actual_network.H +++ b/networks/aprox19/actual_network.H @@ -230,7 +230,7 @@ namespace Rates { } namespace RHS { - + AMREX_GPU_HOST_DEVICE AMREX_INLINE constexpr rhs_t rhs_data (int rate) { diff --git a/networks/aprox21/actual_network.H b/networks/aprox21/actual_network.H index ab2be1cfd6..1367e6fd94 100644 --- a/networks/aprox21/actual_network.H +++ b/networks/aprox21/actual_network.H @@ -226,7 +226,7 @@ namespace Rates { } namespace RHS { - + AMREX_GPU_HOST_DEVICE AMREX_INLINE constexpr rhs_t rhs_data (int rate) { diff --git a/networks/rhs.H b/networks/rhs.H index e507d786aa..628818195b 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -78,7 +78,7 @@ constexpr int is_rate_used () // Determine the index of a given intermediate reaction. We use the // order of the original rate definitions - + // Counts up the number of intermediate reactions. An intermediate // reaction is defined as any reaction which contributes to the // construction of some other reaction. Note that an intermediate diff --git a/networks/triple_alpha_plus_cago/actual_network.H b/networks/triple_alpha_plus_cago/actual_network.H index 55ba0c0491..8d86343d8c 100644 --- a/networks/triple_alpha_plus_cago/actual_network.H +++ b/networks/triple_alpha_plus_cago/actual_network.H @@ -102,7 +102,7 @@ namespace Rates }; namespace RHS { - + AMREX_GPU_HOST_DEVICE AMREX_INLINE constexpr rhs_t rhs_data (int rate) { diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H index dba9fb88a3..15972e4ae9 100644 --- a/nse_solver/nse_check.H +++ b/nse_solver/nse_check.H @@ -35,7 +35,7 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) { if (NSE_INDEX::h1_index == -1 || NSE_INDEX::he4_index == -1) { amrex::Error("Need proton and helium-4 in the network for NSE_NET to work"); } - + // If there are proton, or helium in the network // Check if n,p,a are in equilibrium @@ -655,11 +655,11 @@ bool in_single_group(const amrex::Array1D& group_ind) { // This function checks whether all isotopes are either in the LIG group // or in another single group. - + int LIG_root_index = get_root_index(NSE_INDEX::he4_index, group_ind); int nonLIG_index = -1; - + int oddN_group = -1; int evenN_group = -1; @@ -681,7 +681,7 @@ bool in_single_group(const amrex::Array1D& group_ind) { in_single_group = false; break; } - + } // If there no neutrons are in the network and original condition failed @@ -689,7 +689,7 @@ bool in_single_group(const amrex::Array1D& group_ind) { // There seems to be two big groups after Si28 in NSE: // 1) isotopes with even N // 2) isotopes with odd N - + if (NSE_INDEX::n_index == -1 && !in_single_group) { in_single_group = true; @@ -731,7 +731,7 @@ void fill_merge_indices(amrex::Array1D& merge_indices, const amrex::Array1D& screened_rates, const amrex::Array1D& group_ind, const amrex::Real& t_s) { - + // This function fills in the merge index of the current rate. // The timescale of the current rate must be shorter (faster) than fastest_t @@ -756,15 +756,15 @@ void fill_merge_indices(amrex::Array1D& merge_indices, if (NSE_INDEX::rate_indices(current_rate, k) == -1) { continue; } - + if (NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::h1_index && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::n_index && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::he4_index ) { ++m; - + // skip if there are more than 2 non neutron, proton, or helium-4 in the rate - + if (m > 2) { return; } @@ -794,7 +794,7 @@ void fill_merge_indices(amrex::Array1D& merge_indices, if (root_index == nonLIG_root) { return; - } + } nonLIG_root = root_index; } } @@ -804,7 +804,7 @@ void fill_merge_indices(amrex::Array1D& merge_indices, if (num_nonLIG == 0 || num_nonLIG > 2) { return; } - + // Find the Y_group of the nuclei involved in the net // Y_group are the group molar fractions corresponding to the reactants and products // first 3 for reactants, last 3 for products @@ -817,9 +817,9 @@ void fill_merge_indices(amrex::Array1D& merge_indices, continue; } root_index = get_root_index(NSE_INDEX::rate_indices(current_rate, k), group_ind); - + // let LIG group use their own Y instead of Y_group. (not sure if this is true) - + if (root_index == get_root_index(NSE_INDEX::he4_index, group_ind)) { Y_group(k) = Y(NSE_INDEX::rate_indices(current_rate, k) + 1); } @@ -856,7 +856,7 @@ void fill_merge_indices(amrex::Array1D& merge_indices, } b_r *= Y(NSE_INDEX::rate_indices(current_rate, 5) + 1) * state.rho; } - + // Find the timescale of the rate, See Eq. 17 and Eq. 11 in Kushnir amrex::Real t_i = Y_group(nonNPA_ind(1)) / amrex::min(b_f, b_r); @@ -865,7 +865,7 @@ void fill_merge_indices(amrex::Array1D& merge_indices, } // return if current rate timescale is larger (slower) than previous time scale. - + if (fastest_t < t_i) { return; } @@ -921,12 +921,12 @@ void nse_grouping(amrex::Array1D& group_ind, const burn_t& stat group_ind(i) = 1; } } - + // Now do the grouping based on the timescale. - + amrex::Array1D merge_indices; amrex::Real fastest_t; - + // Go through each reaction and perform grouping while(true) { @@ -960,7 +960,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { // This function returns the boolean that tells whether we're in nse or not // Note that it only works with pynucastro network for now. - + #ifndef NEW_NETWORK_IMPLEMENTATION current_state.nse = false; @@ -987,7 +987,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { state.xn[n] = current_state.y[SFS+n] / current_state.rho; } #endif - + // Check whether state is in the ballpark of NSE if (!skip_molar_check) { @@ -1000,7 +1000,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { // A simple temperature criteria after molar fraction check for determining NSE state // By default, T_nse_net = -1.0 // So this is only enabled if the user provides value in the input file - + if (T_nse_net > 0.0_rt && T_in > T_nse_net) { current_state.nse = true; return current_state.nse; @@ -1011,7 +1011,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { // the thermodynamic condition. // Note we only do this after the first check, which should tell us whether // our current mass fractions are in the ballpark of NSE mass fractions. - + if (nse_molar_independent || skip_molar_check) { state = nse_state; state.dx = current_state.dx; @@ -1050,7 +1050,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { if (!nse_dx_independent) { burn_to_eos(state, eos_state); eos(eos_input_rt, eos_state); - + // a parameter to characterize whether a rate is fast enough t_s = state.dx / eos_state.cs; @@ -1079,7 +1079,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { // fill in energy generation rate to ydot ydot(NumSpec + 1) = enuc - sneut; - + bool found_fast_reaction_cycle = true; // Additional constraint if there is neutron in the network: diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H index efe699144f..307e3730a5 100644 --- a/nse_solver/nse_solver.H +++ b/nse_solver/nse_solver.H @@ -34,7 +34,7 @@ T get_nonexponent_nse_state(const T& state) { T nse_state = state; // set partition function and spin - + amrex::Real pf = 1.0_rt; [[maybe_unused]] amrex::Real dpf_dT; amrex::Real spin = 1.0_rt; @@ -102,7 +102,7 @@ void apply_nse_exponent(T& nse_state) { #endif amrex::Real u_c = 0.0_rt; - + for (int n = 0; n < NumSpec; ++n) { #ifdef NEW_NETWORK_IMPLEMENTATION if (n == NSE_INDEX::h1_index) { @@ -162,7 +162,7 @@ void fcn(Array1D& x, Array1D& fvec, const T& state, int& iflag) { // here state is the nse_state from get_nonexponent_nse_state amrex::ignore_unused(iflag); - + T nse_state = state; nse_state.mu_p = x(1); nse_state.mu_n = x(2); @@ -198,7 +198,7 @@ void jcn(Array1D& x, Array2D& fjac, const T& state, int& iflag) { // here state is the nse_state from get_nonexponent_nse_state amrex::ignore_unused(iflag); - + T nse_state = state; nse_state.mu_p = x(1); nse_state.mu_n = x(2); @@ -355,10 +355,10 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { // store determinant for finding inverse jac amrex::Real det; - + // store inverse jacobian amrex::Array2D inverse_jac; - + // difference in chemical potential of proton and neutron amrex::Real d_mu_p = std::numeric_limits::max(); amrex::Real d_mu_n = std::numeric_limits::max(); @@ -511,7 +511,7 @@ T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt, state.mu_p = nse_state.mu_p; state.mu_n = nse_state.mu_n; - + return nse_state; } #endif diff --git a/rates/aprox_rates.H b/rates/aprox_rates.H index 0cdc9264a7..880efdd8b4 100644 --- a/rates/aprox_rates.H +++ b/rates/aprox_rates.H @@ -8,7 +8,7 @@ using namespace amrex; using namespace network_rp; -inline +inline void rates_init() { // store rates @@ -65,10 +65,10 @@ void rates_init() } } -AMREX_GPU_HOST_DEVICE inline -void rate_c12ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_c12ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { if (use_c12ag_deboer17) { // This version computes the nuclear reaction rate for 12C(a,g)16O and its inverse @@ -209,9 +209,9 @@ void rate_c12ag(const tf_t& tf, const Real den, Real& fr, } -AMREX_GPU_HOST_DEVICE inline -void rate_triplealf(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_triplealf(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, Real& drrdt) { const Real rc28 = 0.1_rt; const Real q1 = 1.0_rt / 0.009604e0_rt; @@ -294,10 +294,10 @@ void rate_triplealf(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_c12c12(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_c12c12(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // c12 + c12 reaction // this is the C12(C12,g)Mg24 rate from Caughlin & Fowler (1988) @@ -328,10 +328,10 @@ void rate_c12c12(const tf_t& tf, const Real den, Real& fr, drrdt = 0.0e0_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_c12o16(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_c12o16(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // c12 + o16 reaction; see cf88 references 47-4 Real term, dtermdt; @@ -363,7 +363,7 @@ void rate_c12o16(const tf_t& tf, const Real den, Real& fr, zz = 1.0e0_rt/cc; term = 1.72e31_rt * t9a56 * tf.t9i32 * std::exp(-106.594_rt/t9a13) * zz; - dtermdt = term*(dt9a56/t9a56 - 1.5e0_rt*tf.t9i + dtermdt = term*(dt9a56/t9a56 - 1.5e0_rt*tf.t9i + 106.594_rt/t9a23*dt9a13 - zz*dcc); } else { @@ -380,20 +380,20 @@ void rate_c12o16(const tf_t& tf, const Real den, Real& fr, drrdt = 0.0e0_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_o16o16(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_o16o16(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // o16 + o16 // this is the O16(O16,g)S32 rate from Caughlin & Fowler (1988) const Real term = 7.10e36_rt * tf.t9i23 * - std::exp(-135.93_rt * tf.t9i13 - 0.629_rt*tf.t923 + std::exp(-135.93_rt * tf.t9i13 - 0.629_rt*tf.t923 - 0.445_rt*tf.t943 + 0.0103_rt*tf.t9*tf.t9); const Real dtermdt = -2.0_rt/3.0_rt*term*tf.t9i - + term * (1.0_rt/3.0_rt*135.93_rt*tf.t9i43 - 2.0_rt/3.0_rt*0.629_rt*tf.t9i13 + + term * (1.0_rt/3.0_rt*135.93_rt*tf.t9i43 - 2.0_rt/3.0_rt*0.629_rt*tf.t9i13 - 4.0_rt/3.0_rt*0.445_rt*tf.t913 + 0.0206_rt*tf.t9); // rates @@ -404,10 +404,10 @@ void rate_o16o16(const tf_t& tf, const Real den, Real& fr, drrdt = 0.0e0_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_o16ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_o16ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real q1 = 1.0e0_rt/2.515396e0_rt; @@ -441,10 +441,10 @@ void rate_o16ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_ne20ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_ne20ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real rc102 = 0.1_rt; const Real q1 = 1.0e0_rt/4.923961e0_rt; @@ -500,10 +500,10 @@ void rate_ne20ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_mg24ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_mg24ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real rc121 = 0.1e0_rt; @@ -546,10 +546,10 @@ void rate_mg24ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_mg24ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_mg24ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real rc148 = 0.1_rt; const Real q1 = 1.0_rt / 0.024649e0_rt; @@ -604,10 +604,10 @@ void rate_mg24ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * dtermdt * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_al27pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_al27pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // al27(p,g)si28 // champagne 1996 @@ -647,10 +647,10 @@ void rate_al27pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_al27pg_old(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_al27pg_old(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real rc147 = 0.1_rt; const Real q1 = 1.0e0_rt/0.024025e0_rt; @@ -706,10 +706,10 @@ void rate_al27pg_old(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_si28ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_si28ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // si28(a,g)s32 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -725,7 +725,7 @@ void rate_si28ag(const tf_t& tf, const Real den, Real& fr, const Real term = 4.82e22_rt * tf.t9i23 * std::exp(-61.015_rt * tf.t9i13 * aa); const Real dtermdt = term*(-2.0_rt/3.0_rt*tf.t9i + 61.015_rt*tf.t9i13*(1.0_rt/3.0_rt*tf.t9i*aa - daa)); - + // the rates fr = den * term; dfrdt = den * dtermdt * 1.0e-9_rt; @@ -737,10 +737,10 @@ void rate_si28ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_si28ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_si28ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // si28(a,p)p31 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -769,10 +769,10 @@ void rate_si28ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * dtermdt * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_p31pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_p31pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // p31(p,g)s32 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -800,10 +800,10 @@ void rate_p31pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_s32ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_s32ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // s32(a,g)ar36 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -831,10 +831,10 @@ void rate_s32ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_s32ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_s32ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // s32(a,p)cl35 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -862,15 +862,15 @@ void rate_s32ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * dtermdt * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_cl35pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_cl35pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // cl35(p,g)ar36 const Real aa = 1.0e0_rt + 1.761e-1_rt*tf.t9 - 1.322e-2_rt*tf.t92 + 5.245e-4_rt*tf.t93; const Real daa = 1.761e-1_rt - 2.0e0_rt*1.322e-2_rt*tf.t9 + 3.0e0_rt*5.245e-4_rt*tf.t92; - + const Real term = 4.48e16_rt * tf.t9i23 * std::exp(-29.483_rt * tf.t9i13 * aa); const Real dtermdt = term*(-2.0_rt/3.0_rt*tf.t9i + 29.483_rt*tf.t9i13*(1.0_rt/3.0_rt*tf.t9i*aa - daa)); @@ -885,10 +885,10 @@ void rate_cl35pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_ar36ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_ar36ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // ar36(a,g)ca40 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -916,10 +916,10 @@ void rate_ar36ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_ar36ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_ar36ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // ar36(a,p)k39 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -947,10 +947,10 @@ void rate_ar36ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * dtermdt * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_k39pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_k39pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // k39(p,g)ca40 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -978,10 +978,10 @@ void rate_k39pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_ca40ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_ca40ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // ca40(a,g)ti44 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1009,10 +1009,10 @@ void rate_ca40ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_ca40ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_ca40ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // ca40(a,p)sc43 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1040,10 +1040,10 @@ void rate_ca40ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * dtermdt * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_sc43pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_sc43pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // sc43(p,g)ca40 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1071,10 +1071,10 @@ void rate_sc43pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_ti44ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_ti44ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // ti44(a,g)cr48 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1102,10 +1102,10 @@ void rate_ti44ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_ti44ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_ti44ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // ti44(a,p)v47 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1133,10 +1133,10 @@ void rate_ti44ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * dtermdt * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_v47pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_v47pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // v47(p,g)cr48 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1164,10 +1164,10 @@ void rate_v47pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_cr48ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_cr48ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // cr48(a,g)fe52 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1195,10 +1195,10 @@ void rate_cr48ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_cr48ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_cr48ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // cr48(a,p)mn51 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1226,10 +1226,10 @@ void rate_cr48ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_mn51pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_mn51pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // mn51(p,g)fe52 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1257,10 +1257,10 @@ void rate_mn51pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe52ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe52ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe52(a,g)ni56 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1288,10 +1288,10 @@ void rate_fe52ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe52ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe52ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe52(a,p)co55 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1319,10 +1319,10 @@ void rate_fe52ap(const tf_t& tf, const Real den, Real& fr, drrdt = den * (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_co55pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_co55pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // co55(p,g)ni56 const Real z = amrex::min(tf.t9, 10.0e0_rt); @@ -1350,10 +1350,10 @@ void rate_co55pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_pp(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_pp(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // p(p,e+nu)d Real term, dtermdt; @@ -1380,10 +1380,10 @@ void rate_pp(const tf_t& tf, const Real den, Real& fr, drrdt = 0.0e0_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_png(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_png(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // p(n,g)d // smith,kawano,malany 1992 @@ -1417,10 +1417,10 @@ void rate_png(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_dpg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_dpg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // d(p,g)he3 const Real aa = 2.24e3_rt * tf.t9i23 * std::exp(-3.720_rt*tf.t9i13); @@ -1443,10 +1443,10 @@ void rate_dpg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_he3ng(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_he3ng(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // he3(n,g)he4 const Real term = 6.62_rt * (1.0e0_rt + 905.0_rt*tf.t9); @@ -1463,10 +1463,10 @@ void rate_he3ng(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_he3he3(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_he3he3(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // he3(he3,2p)he4 const Real aa = 6.04e10_rt * tf.t9i23 * std::exp(-12.276_rt*tf.t9i13); @@ -1491,10 +1491,10 @@ void rate_he3he3(const tf_t& tf, const Real den, Real& fr, drrdt = den * den * (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_he3he4(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_he3he4(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // he3(he4,g)be7 const Real aa = 1.0e0_rt + 0.0495_rt*tf.t9; @@ -1526,10 +1526,10 @@ void rate_he3he4(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_c12pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_c12pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real q1 = 1.0e0_rt/2.25e0_rt; @@ -1565,11 +1565,11 @@ void rate_c12pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_n14pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) -{ +AMREX_GPU_HOST_DEVICE inline +void rate_n14pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) +{ const Real q1 = 1.0e0_rt/10.850436e0_rt; // n14(p,g)o15 @@ -1604,10 +1604,10 @@ void rate_n14pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_n15pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_n15pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real q1 = 1.0_rt / 0.2025_rt; @@ -1646,10 +1646,10 @@ void rate_n15pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_n15pa(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_n15pa(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real theta = 0.1_rt; const Real q1 = 1.0_rt / 0.272484_rt; @@ -1692,10 +1692,10 @@ void rate_n15pa(const tf_t& tf, const Real den, Real& fr, drrdt = den * (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_o16pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_o16pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // o16(p,g)f17 const Real aa = std::exp(-0.728_rt*tf.t923); @@ -1728,10 +1728,10 @@ void rate_o16pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_n14ag(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_n14ag(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { const Real q1 = 1.0e0_rt/0.776161e0_rt; @@ -1770,10 +1770,10 @@ void rate_n14ag(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt*term + rev*dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe52ng(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe52ng(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe52(n,g)fe53 const Real tq2 = tf.t9 - 0.348e0_rt; @@ -1791,10 +1791,10 @@ void rate_fe52ng(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe53ng(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe53ng(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe53(n,g)fe54 const Real tq1 = tf.t9/0.348_rt; @@ -1816,10 +1816,10 @@ void rate_fe53ng(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe54ng(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe54ng(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe54(n,g)fe55 const Real aa = 2.307390e1_rt - 7.931795e-02_rt * tf.t9i + 7.535681e0_rt * tf.t9i13 @@ -1851,10 +1851,10 @@ void rate_fe54ng(const tf_t& tf, const Real den, Real& fr, dfrdt = dtermdt*den; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe54pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe54pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe54(p,g)co55 const Real z = amrex::min(tf.t9,10.0e0_rt); @@ -1882,10 +1882,10 @@ void rate_fe54pg(const tf_t& tf, const Real den, Real& fr, drrdt = (drevdt * term + rev * dtermdt) * 1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe54ap(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe54ap(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe54(a,p)co57 const Real aa = 3.97474900e1_rt - 6.06543100e0_rt * tf.t9i + 1.63239600e2_rt * tf.t9i13 @@ -1917,10 +1917,10 @@ void rate_fe54ap(const tf_t& tf, const Real den, Real& fr, dfrdt = drrdt*bb + rr*dbb*1.0e-9_rt; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe55ng(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe55ng(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe55(n,g)fe56 const Real aa = 1.954115e1_rt - 6.834029e-02_rt * tf.t9i + 5.379859e0_rt * tf.t9i13 @@ -1952,10 +1952,10 @@ void rate_fe55ng(const tf_t& tf, const Real den, Real& fr, dfrdt = dtermdt*den; } -AMREX_GPU_HOST_DEVICE inline -void rate_fe56pg(const tf_t& tf, const Real den, Real& fr, - Real& dfrdt, Real& rr, - Real& drrdt) +AMREX_GPU_HOST_DEVICE inline +void rate_fe56pg(const tf_t& tf, const Real den, Real& fr, + Real& dfrdt, Real& rr, + Real& drrdt) { // fe56(p,g)co57 @@ -1999,10 +1999,10 @@ void rate_fe56pg(const tf_t& tf, const Real den, Real& fr, // output: // rn56ec = ni56 electron capture rate // sn56ec = ni56 neutrino loss rate -AMREX_GPU_HOST_DEVICE inline -void langanke(const Real btemp, const Real bden, - const Real y56, const Real ye, - Real& rn56ec, Real& sn56ec) +AMREX_GPU_HOST_DEVICE inline +void langanke(const Real btemp, const Real bden, + const Real y56, const Real ye, + Real& rn56ec, Real& sn56ec) { // calculate ni56 electron capture and neutrino loss rates rn56ec = 0.0_rt; @@ -2056,10 +2056,10 @@ void langanke(const Real btemp, const Real bden, // positron capture on neutrons rnep (captures/sec/neutron), // and their associated neutrino energy loss rates // spenc (erg/sec/proton) and snepc (erg/sec/neutron) -AMREX_GPU_HOST_DEVICE inline -void ecapnuc(const Real etakep, const Real temp, - Real& rpen, Real& rnep, - Real& spenc, Real& snepc) +AMREX_GPU_HOST_DEVICE inline +void ecapnuc(const Real etakep, const Real temp, + Real& rpen, Real& rnep, + Real& spenc, Real& snepc) { const Real qn1 = -2.0716446e-06_rt; const Real ftinv = 1.0e0_rt/1083.9269e0_rt; diff --git a/rates/aprox_rates_data.H b/rates/aprox_rates_data.H index f58ddca2ae..e9a3162621 100644 --- a/rates/aprox_rates_data.H +++ b/rates/aprox_rates_data.H @@ -16,6 +16,6 @@ extern AMREX_GPU_MANAGED amrex::Array1D rfd2; extern AMREX_GPU_MANAGED amrex::Array1D tfdm; extern AMREX_GPU_MANAGED amrex::Array1D tfd0; extern AMREX_GPU_MANAGED amrex::Array1D tfd1; -extern AMREX_GPU_MANAGED amrex::Array1D tfd2; +extern AMREX_GPU_MANAGED amrex::Array1D tfd2; #endif diff --git a/screening/screen.H b/screening/screen.H index 6ebe21717b..f07fd9f4dd 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -133,7 +133,7 @@ fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Array if constexpr (do_T_derivatives) { state.taufacdt = -(1.0_rt/3.0_rt) * state.taufac * tempi; } - + Real xni = std::cbrt(rr * zbar); // Part of Eq.4 in Itoh:1979 @@ -211,11 +211,11 @@ void actual_screen5 (const plasma_state_t& state, // In Eq.4 in Itoh:1979, this term is 2*Z_1*Z_2/(Z_1^(1/3) + Z_2^(1/3)) // However here we follow Wallace:1982 Eq. A13, which is Z_1*Z_2*(2/(Z_1+Z_2))^(1/3) - + Real qq = fact * bb * scn_fac.zs13inv; // Full Equation of Wallace:1982 Eq. A13 - + Real gamef = qq * gamp; [[maybe_unused]] Real gamefdt; if constexpr (do_T_derivatives) { @@ -225,7 +225,7 @@ void actual_screen5 (const plasma_state_t& state, // Full version of Eq.6 in Itoh:1979 with extra 1/3 factor // the extra 1/3 factor is there for convenience. // tau12 = Eq.6 / 3 - + Real tau12 = state.taufac * scn_fac.aznut; [[maybe_unused]] Real tau12dt; if constexpr (do_T_derivatives) { @@ -235,7 +235,7 @@ void actual_screen5 (const plasma_state_t& state, qq = 1.0_rt/tau12; // alph12 = 3*gamma_ij/tau_ij - + Real alph12 = gamef * qq; [[maybe_unused]] Real alph12dt; if constexpr (do_T_derivatives) { @@ -270,7 +270,7 @@ void actual_screen5 (const plasma_state_t& state, // weak screening regime // Full version of Eq. 19 in Graboske:1973 by considering weak regime // and Wallace:1982 Eq. A14. Here the degeneracy factor is assumed to be 1. - + Real h12w = bb * state.qlam0z; [[maybe_unused]] Real dh12wdt; if constexpr (do_T_derivatives) { @@ -284,11 +284,11 @@ void actual_screen5 (const plasma_state_t& state, } // intermediate and strong sceening regime - + if (gamef > gamefx) { // gamma_ij^(1/4) - + Real gamp14 = std::pow(gamp, 0.25_rt); Real rr = 1.0_rt/gamp; @@ -818,7 +818,7 @@ void chabrier1998 (const plasma_state_t& state, { // Calculates screening factors based on Chabrier & Potekhin 1998, // Calder2007 and partly screen5 routine mentioned in Alastuey 1978. - + // This screening is valid for weak screening: Gamma < 0.1 // and strong screening: 1 <= Gamma <= 160 // Reference: @@ -828,18 +828,18 @@ void chabrier1998 (const plasma_state_t& state, // Alastuey 1978 // Eq. 2 in Chabrier & Potekhin 1998 - + Real Gamma_e = state.gamma_e_fac / state.temp; Real zcomp = scn_fac.z1 + scn_fac.z2; // See Calder2007 appendix Eq. A9 - + Real Gamma1 = Gamma_e * std::pow(scn_fac.z1, 5.0_rt/3.0_rt); 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); Real Gamma1dT{}, Gamma2dT{}, Gamma12dT{}; - + if constexpr (do_T_derivatives) { Gamma1dT = -Gamma1 / state.temp; Gamma2dT = -Gamma2 / state.temp; @@ -854,33 +854,33 @@ void chabrier1998 (const plasma_state_t& state, 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, // which is implemented in the strong screening limit of our screen5 routine. // See Wallace1982, Eq. A13 - + Real Gamma_eff = std::cbrt(2.0_rt) * scn_fac.z1 * scn_fac.z2 * scn_fac.zs13inv * Gamma_e; [[maybe_unused]] Real Gamma_eff_dT; if constexpr (do_T_derivatives) { Gamma_eff_dT = -Gamma_eff / state.temp; } - + // TAU/3, see Wallace1982, Eq. A2 - + Real tau12 = state.taufac * scn_fac.aznut; - + [[maybe_unused]] Real tau12dT; if constexpr (do_T_derivatives) { tau12dT = state.taufacdt * scn_fac.aznut; } // see Calder 2007 Eq. A8 - + Real b_fac = Gamma_eff / tau12; - + [[maybe_unused]] Real b_fac_dT; if constexpr (do_T_derivatives) { b_fac_dT = (Gamma_eff_dT - b_fac * tau12dT) / tau12; @@ -888,10 +888,10 @@ void chabrier1998 (const plasma_state_t& state, // Quantum correction terms (same as screen5) //see Calder 2007 Eq.A8 and Alastuey1978, Eq. 24 and 31 - + Real quantum_corr_1 = -tau12 * (5.0_rt/32.0_rt * std::pow(b_fac, 3) - 0.014_rt * std::pow(b_fac, 4) - 0.128_rt * std::pow(b_fac, 5)); - + Real quantum_corr_2 = -Gamma_eff * (0.0055_rt * std::pow(b_fac, 4) - 0.0098_rt * std::pow(b_fac, 5) + 0.0048_rt * std::pow(b_fac, 6)); @@ -906,23 +906,23 @@ void chabrier1998 (const plasma_state_t& state, * (0.0055_rt * 4.0_rt * std::pow(b_fac, 3) - 0.0098_rt * 5.0_rt * std::pow(b_fac, 4) + 0.0048_rt * 6.0_rt * std::pow(b_fac, 5)); } - + // See Calder2007 Appendix Eq. A8. // f1 + f2 - f12 gives the classical terms // The difference between this and strong screening of screen5 // is that we replaced the classical term which is f1 + f2 - f12 // using results from Chabrier&Potekhin1998. - + Real h12 = f1 + f2 - f12 + quantum_corr_1 + quantum_corr_2; - + [[maybe_unused]] Real dh12dT; if constexpr (do_T_derivatives) { dh12dT = f1dT + f2dT - f12dT + quantum_corr_1_dT + quantum_corr_2_dT; } - + Real h12_max = 300.0_rt; h12 = amrex::min(h12_max, h12); - + scor = std::exp(h12); if constexpr (do_T_derivatives) { @@ -932,7 +932,7 @@ void chabrier1998 (const plasma_state_t& state, scordt = scor * dh12dT; } } - + } template diff --git a/unit_test/burn_cell_primordial_chem/burn_cell.H b/unit_test/burn_cell_primordial_chem/burn_cell.H index ec7a439620..ba9e8388e7 100644 --- a/unit_test/burn_cell_primordial_chem/burn_cell.H +++ b/unit_test/burn_cell_primordial_chem/burn_cell.H @@ -224,9 +224,9 @@ void burn_cell_c() state_over_time << std::setw(10) << t; state_over_time << std::setw(15) << state.rho; - state_over_time << std::setw(12) << state.T; + state_over_time << std::setw(12) << state.T; for (int x = 0; x < NumSpec; ++x){ - state_over_time << std::setw(15) << state.xn[x]; + state_over_time << std::setw(15) << state.xn[x]; } state_over_time << std::endl; diff --git a/unit_test/burn_cell_sdc/burn_cell.H b/unit_test/burn_cell_sdc/burn_cell.H index 6567ad8df3..1a70dd43c7 100644 --- a/unit_test/burn_cell_sdc/burn_cell.H +++ b/unit_test/burn_cell_sdc/burn_cell.H @@ -269,7 +269,7 @@ void burn_cell_c() } #endif - // these need to be initialized + // these need to be initialized burn_state.sdc_iter = 1; burn_state.num_sdc_iters = 1; @@ -325,17 +325,17 @@ void burn_cell_c() state_over_time << std::setw(12) << "Temperature"; for(int x = 0; x < NumSpec; ++x){ const std::string& element = short_spec_names_cxx[x]; - state_over_time << std::setw(12) << element; + state_over_time << std::setw(12) << element; } state_over_time << std::endl; Real t = 0.0; state_over_time << std::setw(10) << t; - state_over_time << std::setw(12) << burn_state.rho; - state_over_time << std::setw(12) << burn_state.T; + state_over_time << std::setw(12) << burn_state.rho; + state_over_time << std::setw(12) << burn_state.T; for (int x = 0; x < NumSpec; ++x){ - state_over_time << std::setw(12) << burn_state.xn[x]; + state_over_time << std::setw(12) << burn_state.xn[x]; } state_over_time << std::endl; @@ -365,7 +365,7 @@ void burn_cell_c() state_over_time << std::setw(12) << burn_state.rho; state_over_time << std::setw(12) << burn_state.T; for (int x = 0; x < NumSpec; ++x){ - state_over_time << std::setw(12) << burn_state.y[SFS+x] / burn_state.rho; + state_over_time << std::setw(12) << burn_state.y[SFS+x] / burn_state.rho; } state_over_time << std::endl; } diff --git a/unit_test/test_aprox_rates/aprox_rates_util.cpp b/unit_test/test_aprox_rates/aprox_rates_util.cpp index 85dd357b78..5df55a7a94 100644 --- a/unit_test/test_aprox_rates/aprox_rates_util.cpp +++ b/unit_test/test_aprox_rates/aprox_rates_util.cpp @@ -139,7 +139,7 @@ void aprox_rates_test(const Box& bx, sp(i, j, k, vars.isi28ag+1) = dfrdt; sp(i, j, k, vars.isi28ag+2) = rr; sp(i, j, k, vars.isi28ag+3) = drrdt; - + rate_si28ap(tf, dens_zone, fr, dfrdt, rr, drrdt); sp(i, j, k, vars.isi28ap) = fr; @@ -415,7 +415,7 @@ void aprox_rates_test(const Box& bx, Real rn56ec; Real sn56ec; - langanke(temp_zone, dens_zone, eos_state.xn[ini56], + langanke(temp_zone, dens_zone, eos_state.xn[ini56], eos_state.y_e, rn56ec, sn56ec); sp(i, j, k, vars.ilanganke) = rn56ec; diff --git a/unit_test/test_aprox_rates/variables.cpp b/unit_test/test_aprox_rates/variables.cpp index 3910b318da..0f506c6fd8 100644 --- a/unit_test/test_aprox_rates/variables.cpp +++ b/unit_test/test_aprox_rates/variables.cpp @@ -64,7 +64,7 @@ plot_t init_variables() { p.ife55ng = p.next_index(4); p.ife56pg = p.next_index(4); - // langanke and ecapnuc are different so not included in n_tests + // langanke and ecapnuc are different so not included in n_tests p.ilanganke = p.next_index(2); p.iecapnuc = p.next_index(4); diff --git a/unit_test/test_ase/burn_cell.H b/unit_test/test_ase/burn_cell.H index 311bea92c0..b4db1497a2 100644 --- a/unit_test/test_ase/burn_cell.H +++ b/unit_test/test_ase/burn_cell.H @@ -27,9 +27,9 @@ void burn_cell_c() state.mu_n = mu_n; // set a reference cell size. - + state.dx = 1.0e6_rt; - + std::cout << "chemical potential of proton is " << mu_p << std::endl; std::cout << "chemical potential of neutron is " << mu_n << std::endl; @@ -41,7 +41,7 @@ void burn_cell_c() std::cout << "State Density (g/cm^3): " << state.rho << std::endl; std::cout << "State Temperature (K): " << state.T << std::endl; std::cout << "electron fraction is " << state.y_e << std::endl; - + std::cout << "NSE state: " << std::endl; for (int n = 0; n < NumSpec; ++n) { std::cout << short_spec_names_cxx[n] << " : " << NSE_STATE.xn[n] << std::endl; @@ -55,12 +55,12 @@ void burn_cell_c() // get eos //eos(eos_input_rt, state); - + if (in_nse(state)){ std::cout << "We're in NSE. " << std::endl; } else{ std::cout << "We're not in NSE. " << std::endl; } - + } diff --git a/unit_test/test_ase/main.cpp b/unit_test/test_ase/main.cpp index a79f2a6e2f..f172e2fe10 100644 --- a/unit_test/test_ase/main.cpp +++ b/unit_test/test_ase/main.cpp @@ -21,7 +21,7 @@ int main(int argc, char *argv[]) { ParmParse ppa("amr"); init_unit_test(); - + // C++ EOS initialization (must be done after Fortran eos_init and init_extern_parameters) eos_init(small_temp, small_dens); diff --git a/unit_test/test_ase/make_table/burn_cell.H b/unit_test/test_ase/make_table/burn_cell.H index b9eb41590c..5ef38f8ac2 100644 --- a/unit_test/test_ase/make_table/burn_cell.H +++ b/unit_test/test_ase/make_table/burn_cell.H @@ -18,7 +18,7 @@ void burn_cell_c() // set a reference cell size state.dx = 1.0e6_rt; - + Real dlogrho = (std::log10(rho_max) - std::log10(rho_min))/static_cast(nrho-1); Real dlogT = (std::log10(T_max) - std::log10(T_min))/static_cast(nT-1); Real dYe = (Ye_max - Ye_min)/(nye-1); diff --git a/unit_test/test_nse/nse_example.H b/unit_test/test_nse/nse_example.H index be91da3284..7b8246c990 100644 --- a/unit_test/test_nse/nse_example.H +++ b/unit_test/test_nse/nse_example.H @@ -42,36 +42,36 @@ void nse_example_c() state.xn[n] = massfractions[n]; } - + // normalize -- just in case - + normalize_abundances_burn(state); // compute the initial Ye - + state.y_e = y_e; - + // composition(state); - + std::cout << "electron fraction is " << state.y_e << std::endl; // set initial chemical potential of proton and neutron state.mu_p = mu_p; state.mu_n = mu_n; - + const bool assume_ye_valid = true; amrex::Real eps = 1.0e-10_rt; - + // find the nse state use_hybrid_solver = 1; - + auto NSE_STATE = get_actual_nse_state(state, eps, assume_ye_valid); std::cout << "After solving: " << std::endl; std::cout << "chemical potential of proton is " << state.mu_p << std::endl; std::cout << "chemical potential of neutron is " << state.mu_n << std::endl; - + std::cout << "NSE state: " << std::endl; for (int n = 0; n < NumSpec; ++n) { std::cout << short_spec_names_cxx[n] << " : " << NSE_STATE.xn[n] << std::endl; diff --git a/unit_test/test_screening/screening_util.cpp b/unit_test/test_screening/screening_util.cpp index 743ed83800..1c2e223904 100644 --- a/unit_test/test_screening/screening_util.cpp +++ b/unit_test/test_screening/screening_util.cpp @@ -53,7 +53,7 @@ void screen_test_C(const Box& bx, const int is32 = network_spec_index("sulfur-32"); if (is32 < 0) amrex::Error("Error: is32 not found"); - const int iar36 = network_spec_index("argon-36"); + const int iar36 = network_spec_index("argon-36"); if (iar36 < 0) amrex::Error("Error: iar36 not found"); const int ica40 = network_spec_index("calcium-40"); diff --git a/unit_test/test_sdc/react_zones.H b/unit_test/test_sdc/react_zones.H index babd078536..bea7dcfc37 100644 --- a/unit_test/test_sdc/react_zones.H +++ b/unit_test/test_sdc/react_zones.H @@ -92,7 +92,7 @@ bool do_react (const plot_t& p, const int i, const int j, const int k, burn_state.reference_time = 0.0; #endif - // these need to be initialized + // these need to be initialized burn_state.sdc_iter = 1; burn_state.num_sdc_iters = 1; diff --git a/unit_test/test_sdc/variables.cpp b/unit_test/test_sdc/variables.cpp index d5f3e1489a..8c0211454d 100644 --- a/unit_test/test_sdc/variables.cpp +++ b/unit_test/test_sdc/variables.cpp @@ -35,7 +35,7 @@ void get_varnames(const plot_t& p, amrex::Vector& names) { } names[p.irho_Hnuc] = "rho_Hnuc"; -#if NAUX_NET > 0 +#if NAUX_NET > 0 for (int n = 0; n < NumAux; n++) { names[p.iaux + n] = aux_names_cxx[n]; names[p.iaux_old + n] = "old_" + aux_names_cxx[n]; diff --git a/unit_test/write_job_info.cpp b/unit_test/write_job_info.cpp index 434fc4f900..2da75b6ee7 100644 --- a/unit_test/write_job_info.cpp +++ b/unit_test/write_job_info.cpp @@ -89,7 +89,7 @@ void write_job_info(const std::string& dir) { jobInfoFile << "COMP version: " << buildInfoGetCompVersion() << "\n"; jobInfoFile << "\n"; - + jobInfoFile << "C++ compiler: " << buildInfoGetCXXName() << "\n"; jobInfoFile << "C++ flags: " << buildInfoGetCXXFlags() << "\n"; From 053bf960922b4ae5949e8e6d11a137ab0e5c14e8 Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Thu, 5 Oct 2023 15:13:08 -0400 Subject: [PATCH 3/4] Fix more tabs found by style CI --- util/skynet/burn.py | 62 ++++----- util/skynet/compare.py | 284 ++++++++++++++++++++--------------------- 2 files changed, 173 insertions(+), 173 deletions(-) diff --git a/util/skynet/burn.py b/util/skynet/burn.py index c6154ec55e..e48f04f63a 100755 --- a/util/skynet/burn.py +++ b/util/skynet/burn.py @@ -29,7 +29,7 @@ def run_skynet(args): strongReactionLibrary = REACLIBReactionLibrary(SkyNetRoot + "/data/reaclib", ReactionType.Strong, True, LeptonMode.TreatAllAsDecayExceptLabelEC, - "Strong reactions", nuclib, opts, True) + "Strong reactions", nuclib, opts, True) symmetricFission = REACLIBReactionLibrary(SkyNetRoot + "/data/netsu_panov_symmetric_0neut", ReactionType.Strong, False, LeptonMode.TreatAllAsDecayExceptLabelEC, @@ -42,8 +42,8 @@ def run_skynet(args): # use only REACLIB weak rates weakReactionLibrary = REACLIBReactionLibrary(SkyNetRoot + "/data/reaclib", - ReactionType.Weak, False, LeptonMode.TreatAllAsDecayExceptLabelEC, - "Weak reactions", nuclib, opts, True) + ReactionType.Weak, False, LeptonMode.TreatAllAsDecayExceptLabelEC, + "Weak reactions", nuclib, opts, True) reactionLibraries = [strongReactionLibrary, symmetricFission, spontaneousFission, weakReactionLibrary] @@ -59,7 +59,7 @@ def run_skynet(args): #Y0[nuclib.NuclideIdsVsNames()["ne23"]] = 4e-4 / 23.0 # above threshold Y0[nuclib.NuclideIdsVsNames()["na23"]] = 4e-4 / 23.0 # below threshold - #XRB + #XRB #Y0[nuclib.NuclideIdsVsNames()["he4"]] = 1.0 / 4.0 #SubCh @@ -74,12 +74,12 @@ def run_skynet(args): #s = 10.0 # initial entropy in k_B / baryon #tau = 7.1 # expansion timescale in ms - time = np.linspace(0.0, 1e9) + time = np.linspace(0.0, 1e9) density_vs_time = ConstantFunction(rho) #isochoric # run through network, will save as "prefix".h5 and .log files output = net.EvolveSelfHeatingWithInitialTemperature(Y0, 0.0, 1.0E1, - T0, density_vs_time, "urca_isoc_belowthresh") + T0, density_vs_time, "urca_isoc_belowthresh") # save some additional information, not necessary as all will be # contained in h5 file @@ -87,37 +87,37 @@ def run_skynet(args): A = np.arange(len(YvsA)) np.savetxt("urca_isoc_belowthresh.txt", np.array([A, YvsA]).transpose(), - "%6i %30.20E") + "%6i %30.20E") return if __name__ == '__main__': - num_cores = multiprocessing.cpu_count() - print "Running with %i worker threads" % num_cores - pool = multiprocessing.Pool(num_cores) - num_done.value = 0 - - # use if you have trajectory files to read - #files = os.listdir(input_dir) - - # otherwise, give the initial density - #densities = [1e6, 5e6, 1e7, 5e7, 1e8, 5e8, 1e9] # multiple - densities = [1.0e9] # URCA - #densities = [2e6] #XRB Aprox13 - #densities = [1e6] #SubCh - - # number of jobs to run - size = len(densities) - args = [(r, size) for r in densities] - - # run skynet in parallel - pool.map_async(run_skynet, args) - - # done submitting jobs - pool.close() - pool.join() + num_cores = multiprocessing.cpu_count() + print "Running with %i worker threads" % num_cores + pool = multiprocessing.Pool(num_cores) + num_done.value = 0 + + # use if you have trajectory files to read + #files = os.listdir(input_dir) + + # otherwise, give the initial density + #densities = [1e6, 5e6, 1e7, 5e7, 1e8, 5e8, 1e9] # multiple + densities = [1.0e9] # URCA + #densities = [2e6] #XRB Aprox13 + #densities = [1e6] #SubCh + + # number of jobs to run + size = len(densities) + args = [(r, size) for r in densities] + + # run skynet in parallel + pool.map_async(run_skynet, args) + + # done submitting jobs + pool.close() + pool.join() # monitor progress with # $ cd work_dir diff --git a/util/skynet/compare.py b/util/skynet/compare.py index 48e6d7684f..c728082c0b 100755 --- a/util/skynet/compare.py +++ b/util/skynet/compare.py @@ -14,151 +14,151 @@ # Get Star Killer output class starkiller_data(): - # modified from Microphysics/unit_test/burn_cell/burn_cell.py - def __init__(self, inpath, runprefix='react_urca'): - files = glob.glob(inpath + runprefix + r'_[0-9]*') - data = [] - for fn in files: - d = {} - f = open(fn, 'r') - - # Get file number - fnsplit = fn.split('_') - fnum = int(fnsplit[-1]) - d['fnum'] = fnum - - # Eat lines - flines = [] - for l in f: - flines.append(l.strip()) - f.close() - - # Fill dataset - n = len(flines) - for i, xi in enumerate(flines): - if xi[-1] == ':': - # This is a field - key = xi[:-1] - xd = [] - # Collect following lines of field data - for j in range(i+1,n): - xj = flines[j] - if xj[-1] == ':': - break - else: - xd.append(xj.split()) - - # Interpret data - if (key=='nspec' or key=='neqs' or key=='i' or key=='j' or key=='k' - or key=='n_rhs' or key=='n_jac'): - # Integer - d[key] = int(xd[0][0]) - elif key=='short_spec_names': - # Strings of nuclide abbreviations - d[key] = [xj[0] for xj in xd] - elif key=='self_heat': - # Boolean - d[key] = (xd[0][0]=='T') - elif key=='xn' or key=='ydot': - # 1-D Float Array - d[key] = np.array([float(y) for y in xd[0]]) - elif key=='jac': - # 2-D Float Array - d[key] = np.array([[float(xij) for xij in xdi] for xdi in xd]) - else: - # Float - d[key] = float(xd[0][0]) - - # Store data - data.append(d) - - ## SORT file data by file number - data_sorted = sorted(data, key=lambda d: d['fnum']) - data = data_sorted - - ## INIT VARIABLES - nspec = data[0]['nspec'] - neqs = data[0]['neqs'] - self.short_spec_names = data[0]['short_spec_names'] - - # Init time, temp, ener - fnum = [] - time = [] - temp = [] - ener = [] - - # Init abundance vectors - xn = [[] for i in range(nspec)] - - # Init ydot vectors - ydot = [[] for i in range(neqs)] - - ## DATA LOOP - # Loop through data and collect - for d in data: - fnum.append(d['fnum']) - temp.append(d['T']) - ener.append(d['e']) - time.append(d['time']) - for i, xi in enumerate(d['xn']): - xn[i].append(xi) - for i, xi in enumerate(d['ydot']): - ydot[i].append(xi) - - # Convert data to numpy arrays - for i in range(nspec): - xn[i] = np.array(xn[i]) - self.X = xn - #for i in range(neqs): - # ydot[i] = np.array(ydot[i]) - # ydot[i][0] = ydot[i][1] - self.fnum = np.array(fnum) - self.temp = np.array(temp) - self.time = np.array(time) - self.dtime = np.ediff1d(time, to_begin=0) - self.ener = np.array(ener) - denerdt = np.zeros_like(ener) - denerdt[1:] = self.ener[1:]/self.dtime[1:] - self.edot = denerdt - # for now ydot is garbage might be useful later? + # modified from Microphysics/unit_test/burn_cell/burn_cell.py + def __init__(self, inpath, runprefix='react_urca'): + files = glob.glob(inpath + runprefix + r'_[0-9]*') + data = [] + for fn in files: + d = {} + f = open(fn, 'r') + + # Get file number + fnsplit = fn.split('_') + fnum = int(fnsplit[-1]) + d['fnum'] = fnum + + # Eat lines + flines = [] + for l in f: + flines.append(l.strip()) + f.close() + + # Fill dataset + n = len(flines) + for i, xi in enumerate(flines): + if xi[-1] == ':': + # This is a field + key = xi[:-1] + xd = [] + # Collect following lines of field data + for j in range(i+1,n): + xj = flines[j] + if xj[-1] == ':': + break + else: + xd.append(xj.split()) + + # Interpret data + if (key=='nspec' or key=='neqs' or key=='i' or key=='j' or key=='k' + or key=='n_rhs' or key=='n_jac'): + # Integer + d[key] = int(xd[0][0]) + elif key=='short_spec_names': + # Strings of nuclide abbreviations + d[key] = [xj[0] for xj in xd] + elif key=='self_heat': + # Boolean + d[key] = (xd[0][0]=='T') + elif key=='xn' or key=='ydot': + # 1-D Float Array + d[key] = np.array([float(y) for y in xd[0]]) + elif key=='jac': + # 2-D Float Array + d[key] = np.array([[float(xij) for xij in xdi] for xdi in xd]) + else: + # Float + d[key] = float(xd[0][0]) + + # Store data + data.append(d) + + ## SORT file data by file number + data_sorted = sorted(data, key=lambda d: d['fnum']) + data = data_sorted + + ## INIT VARIABLES + nspec = data[0]['nspec'] + neqs = data[0]['neqs'] + self.short_spec_names = data[0]['short_spec_names'] + + # Init time, temp, ener + fnum = [] + time = [] + temp = [] + ener = [] + + # Init abundance vectors + xn = [[] for i in range(nspec)] + + # Init ydot vectors + ydot = [[] for i in range(neqs)] + + ## DATA LOOP + # Loop through data and collect + for d in data: + fnum.append(d['fnum']) + temp.append(d['T']) + ener.append(d['e']) + time.append(d['time']) + for i, xi in enumerate(d['xn']): + xn[i].append(xi) + for i, xi in enumerate(d['ydot']): + ydot[i].append(xi) + + # Convert data to numpy arrays + for i in range(nspec): + xn[i] = np.array(xn[i]) + self.X = xn + #for i in range(neqs): + # ydot[i] = np.array(ydot[i]) + # ydot[i][0] = ydot[i][1] + self.fnum = np.array(fnum) + self.temp = np.array(temp) + self.time = np.array(time) + self.dtime = np.ediff1d(time, to_begin=0) + self.ener = np.array(ener) + denerdt = np.zeros_like(ener) + denerdt[1:] = self.ener[1:]/self.dtime[1:] + self.edot = denerdt + # for now ydot is garbage might be useful later? class skydata(): - def __init__(self, inpath, pressure = False): - - out = NetworkOutput.ReadFromFile(inpath) - self.time = np.array(out.Times()) - self.temp = np.array(out.TemperatureVsTime()) * 1.0E9 - self.rho = np.array(out.DensityVsTime()) - self.edot = np.array(out.HeatingRateVsTime()) - self.Y = np.array(out.YVsTime()) - self.A = np.array(out.A()) - self.Z = np.array(out.Z()) - self.X = np.array([self.A * y for y in self.Y[:]]) - - if pressure: - jmax = len(self.temp) - - # writing intermediate thermo file for EoS - # calls the Helmholtz EoS for pressure as a.out - f = open("tmp.txt", 'w') - f.write(str(jmax) + "\n") - - for i in range(jmax): - wbar = 1.0 / np.sum(self.Y[i,:]) - abar = wbar * np.sum(np.dot(self.A,self.Y[i,:])) - zbar = wbar * np.sum(np.dot(self.Z,self.Y[i,:])) - f.write("{:.15e} {:.15e} {:.15e} {:.15e}\n".format(self.temp[i], - self.rho[i], abar, zbar)) - f.close() - - process = subprocess.Popen("./a.out", shell=True) - process.wait() - self.pressure = np.loadtxt("tmp_out.txt") - os.system("rm tmp*") - - def isoidx(self, a, z): + def __init__(self, inpath, pressure = False): + + out = NetworkOutput.ReadFromFile(inpath) + self.time = np.array(out.Times()) + self.temp = np.array(out.TemperatureVsTime()) * 1.0E9 + self.rho = np.array(out.DensityVsTime()) + self.edot = np.array(out.HeatingRateVsTime()) + self.Y = np.array(out.YVsTime()) + self.A = np.array(out.A()) + self.Z = np.array(out.Z()) + self.X = np.array([self.A * y for y in self.Y[:]]) + + if pressure: + jmax = len(self.temp) + + # writing intermediate thermo file for EoS + # calls the Helmholtz EoS for pressure as a.out + f = open("tmp.txt", 'w') + f.write(str(jmax) + "\n") + + for i in range(jmax): + wbar = 1.0 / np.sum(self.Y[i,:]) + abar = wbar * np.sum(np.dot(self.A,self.Y[i,:])) + zbar = wbar * np.sum(np.dot(self.Z,self.Y[i,:])) + f.write("{:.15e} {:.15e} {:.15e} {:.15e}\n".format(self.temp[i], + self.rho[i], abar, zbar)) + f.close() + + process = subprocess.Popen("./a.out", shell=True) + process.wait() + self.pressure = np.loadtxt("tmp_out.txt") + os.system("rm tmp*") + + def isoidx(self, a, z): ''' return index of isotope in the self.z, self.z, self.isos, self.misos, self.xisos[i,:] arrays @@ -168,7 +168,7 @@ def isoidx(self, a, z): print('isotope not found') else: return wh[0] - + #sky = skydata(inpath = "./urca_abovethresh/urca_isoc_abovethresh.h5", pressure = False) #star = starkiller_data("./urca_abovethresh/", runprefix='react_urca') From 23b5b368f5359cf6036f3f61fa392518f1bda3b3 Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Thu, 5 Oct 2023 15:24:25 -0400 Subject: [PATCH 4/4] Fix trailing spaces found by style CI --- .../workflows/burn_cell_primordial_chem.yml | 2 +- .github/workflows/create_release.yml | 4 ++-- .github/workflows/find_changed_files.py | 2 +- .github/workflows/get_release_txt.py | 6 +++--- CMakeLists.txt | 18 ++++++++--------- deploy_docs_action.sh | 6 +++--- networks/ase/ase.py | 4 ++-- networks/ignition_chamulak/notes.txt | 6 +++--- networks/rprox/reaclib/make_rates_module.py | 20 +++++++++---------- sphinx_docs/notes.txt | 4 ++-- sphinx_docs/source/basics.rst | 2 +- sphinx_docs/source/integrators.rst | 8 ++++---- sphinx_docs/source/networks-overview.rst | 2 +- sphinx_docs/source/networks.rst | 4 ++-- sphinx_docs/source/nse.rst | 16 +++++++-------- sphinx_docs/source/sdc.rst | 6 +++--- sphinx_docs/source/templated_networks.rst | 4 ++-- sphinx_docs/source/unit_tests.rst | 2 +- .../test_react/test_network_combinations.py | 2 +- util/code_checker/code_checker.py | 20 +++++++++---------- 20 files changed, 69 insertions(+), 69 deletions(-) diff --git a/.github/workflows/burn_cell_primordial_chem.yml b/.github/workflows/burn_cell_primordial_chem.yml index b1c506ba90..bdc52cc939 100644 --- a/.github/workflows/burn_cell_primordial_chem.yml +++ b/.github/workflows/burn_cell_primordial_chem.yml @@ -47,7 +47,7 @@ jobs: value2=$(awk 'NR=='"$line_number"' {match($0, /[+-]?[0-9]+([.][0-9]+)?[eE]?[+-]?[0-9]+/); if (RSTART) print substr($0, RSTART, RLENGTH)}' reference_solution.out) difference=$(awk -v val1="$value1" -v val2="$value2" 'BEGIN { printf "%.2f", (val1 - val2) / val2 }') - + if (( $(echo "$difference > $threshold" | bc -l) )); then echo "Line number: $line_number" echo "Value in test.out: $value1" diff --git a/.github/workflows/create_release.yml b/.github/workflows/create_release.yml index 2e7d0fba2e..3d021ab174 100644 --- a/.github/workflows/create_release.yml +++ b/.github/workflows/create_release.yml @@ -2,8 +2,8 @@ on: push: # Sequence of patterns matched against refs/tags tags: - - '[0-9][0-9].[0-9][0-9]' - + - '[0-9][0-9].[0-9][0-9]' + name: Create Release jobs: diff --git a/.github/workflows/find_changed_files.py b/.github/workflows/find_changed_files.py index 0e7f4e3e06..8a81d06995 100644 --- a/.github/workflows/find_changed_files.py +++ b/.github/workflows/find_changed_files.py @@ -28,7 +28,7 @@ def find_files(SHAs=None): if stderr is not None: raise Exception('git diff encountered an error') - files = [f for f in stdout.decode('utf-8').strip().split('\n') + files = [f for f in stdout.decode('utf-8').strip().split('\n') if f.startswith('networks/')] print(files) diff --git a/.github/workflows/get_release_txt.py b/.github/workflows/get_release_txt.py index 50530c14ec..e14f92dc28 100644 --- a/.github/workflows/get_release_txt.py +++ b/.github/workflows/get_release_txt.py @@ -26,9 +26,9 @@ txt = txt[m.end():].strip() else: txt = "" - - # we now need to substitute characters in the string so that - # the action can deal with line breaks + + # we now need to substitute characters in the string so that + # the action can deal with line breaks txt = txt.replace('%', '%25') txt = txt.replace('\n', '%0A') txt = txt.replace('\r', '%0D') diff --git a/CMakeLists.txt b/CMakeLists.txt index c3e2702b96..d1aceaff72 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,13 +32,13 @@ function(setup_target_for_microphysics_compilation network_name output_dir) set(networkdir "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/") set(networkheadertemplatefile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/network_header.template") - set (gamma_law_dirs ${CMAKE_BINARY_DIR} ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/gcem/include + set (gamma_law_dirs ${CMAKE_BINARY_DIR} ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/gcem/include ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/gamma_law - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/constants PARENT_SCOPE) - execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${EOSparamfile} ${networkparamfile} + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${EOSparamfile} ${networkparamfile} ${VODEparamfile} ${integrationparamfile}" --use_namespace WORKING_DIRECTORY ${output_dir}/) set(gamma_law_sources ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces/eos_data.cpp @@ -57,12 +57,12 @@ function(setup_target_for_microphysics_compilation network_name output_dir) set(networkheadertemplatefile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/network_header.template") #DO NOT change the order of the directories below! - set (primordial_chem_dirs ${CMAKE_BINARY_DIR} ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/gcem/include - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/VODE ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/utils + set (primordial_chem_dirs ${CMAKE_BINARY_DIR} ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/gcem/include + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/VODE ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/utils ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces ${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/primordial_chem ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/constants PARENT_SCOPE) #we need to have extern_parameters.cpp be available at configure time @@ -72,11 +72,11 @@ function(setup_target_for_microphysics_compilation network_name output_dir) #to generate updated header files if(BUILD_UNIT_TEST) - execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${paramfile} ${EOSparamfile} + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${paramfile} ${EOSparamfile} ${networkpcparamfile} ${networkparamfile} ${VODEparamfile} ${integrationparamfile} ${unittestparamfile}" --use_namespace WORKING_DIRECTORY ${output_dir}/) else() #do not need paramfile and unittestparamfile - execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${EOSparamfile} ${networkpcparamfile} + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${EOSparamfile} ${networkpcparamfile} ${networkparamfile} ${VODEparamfile} ${integrationparamfile} " --use_namespace WORKING_DIRECTORY ${output_dir}/) endif() diff --git a/deploy_docs_action.sh b/deploy_docs_action.sh index 8e03fff881..734cbf30ce 100755 --- a/deploy_docs_action.sh +++ b/deploy_docs_action.sh @@ -9,10 +9,10 @@ TARGET_BRANCH="gh-pages" mkdir out -# if on the dev branch, use the dev_layout.html template to get the +# if on the dev branch, use the dev_layout.html template to get the # links correct if [ "$GITHUB_BRANCH" = "$DEV_BRANCH" ]; then - mv sphinx_docs/source/_templates/dev_layout.html sphinx_docs/source/_templates/layout.html + mv sphinx_docs/source/_templates/dev_layout.html sphinx_docs/source/_templates/layout.html fi # Build the Sphinx documentation @@ -24,7 +24,7 @@ mkdir -p out/docs/ if [ "$GITHUB_BRANCH" = "$MAIN_BRANCH" ]; then mkdir -p out/docs mv sphinx_docs/build/html/* out/docs -else +else mkdir -p out/docs/dev/ mv sphinx_docs/build/html/* out/docs/dev fi diff --git a/networks/ase/ase.py b/networks/ase/ase.py index 81663403ce..273bcbc310 100644 --- a/networks/ase/ase.py +++ b/networks/ase/ase.py @@ -56,14 +56,14 @@ def get_library(): print("removing: ", r) _r = subch.get_rate_by_name(r) subch.remove_rate(_r) - + # additional neutron rates to remove for r in subch.get_rates(): if (r == subch.get_rate_by_name("mg24(n,a)ne21") or r == subch.get_rate_by_name("ne21(a,n)mg24") or r == subch.get_rate_by_name("na22(n,g)na23") or r == subch.get_rate_by_name("na23(g,n)na22") or r == subch.get_rate_by_name("ne21(g,n)ne20") or r == subch.get_rate_by_name("ne20(n,g)ne21")): continue - + if pyna.Nucleus("n") in r.reactants or pyna.Nucleus("n") in r.products: print("removing neutron rates: ", r) subch.remove_rate(r) diff --git a/networks/ignition_chamulak/notes.txt b/networks/ignition_chamulak/notes.txt index d4c78368f9..08b20cda0f 100644 --- a/networks/ignition_chamulak/notes.txt +++ b/networks/ignition_chamulak/notes.txt @@ -38,12 +38,12 @@ Basic idea: --- --- this means that 60% of the time, we will create Ne20 + O16 and - 40% of the time we will create Na23 + C13. + 40% of the time we will create Na23 + C13. So the number of nuclei made from consuming 6 C12 is ~ 1.2 Ne20, 1.2 O16, 0.8 Na23, 0.8 C13. - + Doing this, the ash has an A = 18 and a Z = 8.8 We also note that since we are operating in a high temperature @@ -61,4 +61,4 @@ Basic idea: release per C12 consumed), as a function of density. - + diff --git a/networks/rprox/reaclib/make_rates_module.py b/networks/rprox/reaclib/make_rates_module.py index 12960965d1..48c4cf20b9 100644 --- a/networks/rprox/reaclib/make_rates_module.py +++ b/networks/rprox/reaclib/make_rates_module.py @@ -50,7 +50,7 @@ def modulePreamble(): def tFactorData(): """Add a routine to initialize the tfactors given the temperature. Return this text.""" - + text=""" contains @@ -74,7 +74,7 @@ def tFactorData(): def buildRateRoutine(name,aFactors): """Given a reaction name and the list of aFactors from the ReacLib rate - file, build the appropriate structure for the reaction rate as a + file, build the appropriate structure for the reaction rate as a subroutine. Return this text.""" text= """ subroutine %s(tfactors,rate,dratedt) @@ -142,7 +142,7 @@ def buildRateRoutine(name,aFactors): if not aFactors[start+4] == 0.0: str += " +ct9 &\n" if not aFactors[start+5] == 0.0: - str += " +FIVE3RD*ct953 &\n" + str += " +FIVE3RD*ct953 &\n" if not aFactors[start+6] == 0.0: str += " +(%sd0))" % (aFactors[start+6]) else: @@ -177,7 +177,7 @@ def buildRateRoutine(name,aFactors): return text def numReactantsProducts(chapter): - """Return the number of reactants and products (in particles) for a + """Return the number of reactants and products (in particles) for a reaction of ReacLib Chapter type chapter.""" r,p = 1,1 if chapter in [2,5,9,10]: p+= 1 @@ -190,11 +190,11 @@ def numReactantsProducts(chapter): return r,p def parseRateFile(file): - """Given a ReacLib v 1 rate file, file, parse it and determine an + """Given a ReacLib v 1 rate file, file, parse it and determine an appropriate name and the "a" factors that build the rate. Return the name as a string and the "a" factors as a list.""" - try: + try: fh = open(file,'r') except IOError: print "Couldn't open file %f for reading" % file @@ -207,7 +207,7 @@ def parseRateFile(file): foundLabels = False rateString = "rate" - + aFactors = [] for line in data[numHeaderLines:]: @@ -226,7 +226,7 @@ def parseRateFile(file): for term in range(nterms): start = term*lengthEntry stop = start + lengthEntry - + aFactors.append(float(line[start:stop])) if len(aFactors) % numTerms is not 0: @@ -261,8 +261,8 @@ def weakRate(name,factor): return textParameter, textSubroutine def make_rates_module(rateFiles): - """Build the F90 module containing the subroutines for the reactions - given in the list of ReacLib rate files, rateFiles, along with the helper + """Build the F90 module containing the subroutines for the reactions + given in the list of ReacLib rate files, rateFiles, along with the helper routine to set the common temperature factors. Return all this as a string for easy piping.""" diff --git a/sphinx_docs/notes.txt b/sphinx_docs/notes.txt index a15e64f36b..9b6348fcd3 100644 --- a/sphinx_docs/notes.txt +++ b/sphinx_docs/notes.txt @@ -6,10 +6,10 @@ https://stackoverflow.com/questions/11347505/what-are-some-approaches-to-outputting-a-python-data-structure-to-restructuredte - + ** index: - use + use ..index:: burn_t diff --git a/sphinx_docs/source/basics.rst b/sphinx_docs/source/basics.rst index b5c8933bef..dcf86aafcd 100644 --- a/sphinx_docs/source/basics.rst +++ b/sphinx_docs/source/basics.rst @@ -38,7 +38,7 @@ This will create an executable called ``main3d.gnu.ex``. Then you can run it as ./main3d.gnu.ex inputs_aprox21 By default, the test is built with the 21-isotope ``aprox21`` network. -Here ``inputs_aprox21`` is the inputs file that sets options. +Here ``inputs_aprox21`` is the inputs file that sets options. diff --git a/sphinx_docs/source/integrators.rst b/sphinx_docs/source/integrators.rst index 77c59767df..c06320f946 100644 --- a/sphinx_docs/source/integrators.rst +++ b/sphinx_docs/source/integrators.rst @@ -22,7 +22,7 @@ The equations we integrate to do a nuclear burn are: Here, :math:`X_k` is the mass fraction of species :math:`k`, :math:`e` is the specific nuclear energy created through reactions. Also needed are density :math:`\rho`, -temperature :math:`T`, and the specific heat. The function :math:`f` provides the energy release from reactions and can often be expressed in terms of the +temperature :math:`T`, and the specific heat. The function :math:`f` provides the energy release from reactions and can often be expressed in terms of the instantaneous reaction terms, :math:`\dot{X}_k`. As noted in the previous section, this is implemented in a network-specific manner. @@ -70,7 +70,7 @@ routine (at the moment this can be ``VODE``, ``BackwardEuler``, ``ForwardEuler`` AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void burner (burn_t& state, Real dt) -The input is a ``burn_t``. +The input is a ``burn_t``. .. note:: @@ -369,7 +369,7 @@ the allowed options are: The spectral radius is estimated by default using the power method, built into RKC. Alternately, by setting ``integrator.use_circle_theorem=1``, - the `Gershgorin circle theorem `_ + the `Gershgorin circle theorem `_ is used instead. * ``VODE``: the VODE :cite:`vode` integration package. We ported this @@ -383,7 +383,7 @@ robust. .. important:: The integrator will not abort if it encounters trouble. Instead it will - set ``burn_t burn_state.success = false`` on exit. It is up to the + set ``burn_t burn_state.success = false`` on exit. It is up to the application code to handle the failure. .. note:: diff --git a/sphinx_docs/source/networks-overview.rst b/sphinx_docs/source/networks-overview.rst index ec67cc2fd7..f4beb95e9a 100644 --- a/sphinx_docs/source/networks-overview.rst +++ b/sphinx_docs/source/networks-overview.rst @@ -42,7 +42,7 @@ for the ability to swap integrators as desired. We discuss the integrators in a later section. A network is defined by a ``.net`` file which provides a list of species -and some data about each species (its name and some isotopic data). +and some data about each species (its name and some isotopic data). An example is ``Microphysics/networks/iso7/iso7.net``: diff --git a/sphinx_docs/source/networks.rst b/sphinx_docs/source/networks.rst index 2552e0d953..065a3d0f00 100644 --- a/sphinx_docs/source/networks.rst +++ b/sphinx_docs/source/networks.rst @@ -158,7 +158,7 @@ X_a / (A_a m_u)`, our rate equation is \begin{align} \frac{dX_f}{dt} &= - \frac{r_0}{m_u} \rho X_f^2 \frac{1}{A_f} \left (\frac{T}{T_0}\right)^\nu \equiv \omegadot_f \\ - \frac{dX_a}{dt} &= \frac{1}{2}\frac{r_0}{m_u} \rho X_f^2 \frac{A_a}{A_f^2} \left (\frac{T}{T_0}\right)^\nu = \frac{r_0}{m_u} \rho X_f^2 \frac{1}{A_f} \left (\frac{T}{T_0}\right)^\nu + \frac{dX_a}{dt} &= \frac{1}{2}\frac{r_0}{m_u} \rho X_f^2 \frac{A_a}{A_f^2} \left (\frac{T}{T_0}\right)^\nu = \frac{r_0}{m_u} \rho X_f^2 \frac{1}{A_f} \left (\frac{T}{T_0}\right)^\nu \end{align} We define a new rate constant, :math:`\rt` with units of :math:`[\mathrm{s^{-1}}]` as @@ -287,6 +287,6 @@ to disable rates: :math:`\isotm{N}{13}(\alpha,p)\isotm{O}{16}` and its inverse are disabled. -Together, these parameters allow us to turn off the sequence +Together, these parameters allow us to turn off the sequence :math:`\isotm{C}{12}(p,\gamma)\isotm{N}{13}(\alpha, p)\isotm{O}{16}` that acts as a bypass for :math:`\isotm{C}{12}(\alpha, \gamma)\isotm{O}{16}`. diff --git a/sphinx_docs/source/nse.rst b/sphinx_docs/source/nse.rst index d94816ba8b..c4a207d6ce 100644 --- a/sphinx_docs/source/nse.rst +++ b/sphinx_docs/source/nse.rst @@ -254,7 +254,7 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. An example of such reaction cycle would be: .. math:: - + \isotm{S}{32} (\gamma, p)(\gamma, p)(\gamma, n)(\gamma, n) \isotm{Si}{28} (\alpha, \gamma) \isotm{S}{32} @@ -265,11 +265,11 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. If there are no isotope present in the network that would form a closed-cycle, we move on to the next nuclei. We break out of the iteration once we found a fast reaction cycle. - + * If the previous two check pass, we proceed to nuclei grouping. Initially, :math:`p`, :math:`n`, and :math:`\alpha` are grouped into a single group called the light-isotope-group, or LIG. Other isotopes belong to their own group, - which only contains themselves. During each iteration, we find all valid reaction, + which only contains themselves. During each iteration, we find all valid reaction, :math:`k`, that has the fastest time-scale, :math:`t_{i,k} = \tilde{Y}_i/\textbf{min}(b_f(k), b_r(k))`, for :math:`i` to be the isotope involved with the reaction that is different from :math:`p`, :math:`n`, and :math:`\alpha`. @@ -283,13 +283,13 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. .. math:: t_{i,k} < \epsilon t_s - + * .. math:: 2|b_f(k) - b_r(k)|/(b_f(k) + b_r(k) < \epsilon - + Here we only consider two cases of reactions: * There are exactly two isotopes involved in reaction, :math:`k`, that are not in the @@ -307,7 +307,7 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. * Reactions that involve more than 2 reactants and products * Reactions that have more than 2 non-light-isotope-group. - + * The nuclei that participate in the reaction is either in LIG or in another group. This means that the non-LIG nuclei have already merged. @@ -320,8 +320,8 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. a single group due to the missing neutron rates. Therefore, there is an alternative criteria of defining a "single group" when neutron is not present in the network: for isotopes, :math:`Z >= 14`, isotopes with odd and even :math:`N` form two - distinct groups. - + distinct groups. + Additional Options ------------------ diff --git a/sphinx_docs/source/sdc.rst b/sphinx_docs/source/sdc.rst index 8db7a2d41b..379874f75f 100644 --- a/sphinx_docs/source/sdc.rst +++ b/sphinx_docs/source/sdc.rst @@ -53,9 +53,9 @@ The reactions don’t modify the total density, :math:`\rho`, or momentum, .. math:: - \frac{d}{dt}\left ( - \begin{array}{c} \rho X_k \\ \rho e \end{array} - \right ) = + \frac{d}{dt}\left ( + \begin{array}{c} \rho X_k \\ \rho e \end{array} + \right ) = \left ( \begin{array}{c} \Adv{\rho X_k}^{n+1/2} \\ \Adv{\rho e}^{n+1/2} \\ \end{array} \right ) + diff --git a/sphinx_docs/source/templated_networks.rst b/sphinx_docs/source/templated_networks.rst index e8342dc195..c6bb55cdae 100644 --- a/sphinx_docs/source/templated_networks.rst +++ b/sphinx_docs/source/templated_networks.rst @@ -181,7 +181,7 @@ special cases (e.g., for approximate nets): This indicates that each branch happens 50% of the time. -* ``apply_identical_particle_factor`` : +* ``apply_identical_particle_factor`` : Normally for rates involving identical nuclei, we divide the rate by a factor (:math:`n!`, where `n` is the number of the same nuclei participating). This @@ -216,7 +216,7 @@ special cases (e.g., for approximate nets): To save on computation, we can create a table of reaction rates by evaluating over a grid of temperature and then interpolating - in temperature as needed. This operates only on the + in temperature as needed. This operates only on the :math:`N_A \langle \sigma v \rangle` portion of the rate. Some rates are more complex than fits into the rate tabulation diff --git a/sphinx_docs/source/unit_tests.rst b/sphinx_docs/source/unit_tests.rst index 638e739556..39b6954836 100644 --- a/sphinx_docs/source/unit_tests.rst +++ b/sphinx_docs/source/unit_tests.rst @@ -164,7 +164,7 @@ Screening Test routine, using the ``aprox21`` reaction network. This uses the same basic ideas as the tests above---a cube of data is setup and the rates are evaluated using each zone's thermodynamic -conditions. +conditions. This works for both the Fortran and C++ implementations (via ``do_cxx``). diff --git a/unit_test/test_react/test_network_combinations.py b/unit_test/test_react/test_network_combinations.py index be1f9b4dd5..eadd7fdc09 100755 --- a/unit_test/test_react/test_network_combinations.py +++ b/unit_test/test_react/test_network_combinations.py @@ -192,7 +192,7 @@ def doit(): for k, v in sorted(outcomes.items()): print("{}: {}".format(k, v)) - + if __name__ == "__main__": doit() diff --git a/util/code_checker/code_checker.py b/util/code_checker/code_checker.py index 93217100c0..9a4695d63b 100644 --- a/util/code_checker/code_checker.py +++ b/util/code_checker/code_checker.py @@ -11,7 +11,7 @@ ignore_dirs = ['tmp_build_dir', 't', 'python_library'] def find_fortran_files(): - # find Microphysics Fortran source files + # find Microphysics Fortran source files try: microphysics_home = os.environ['MICROPHYSICS_HOME'] except KeyError: @@ -31,7 +31,7 @@ def pytest_generate_tests(metafunc): def test_double_precision(filename): if any([f'/{s}/' in str(filename) for s in ignore_dirs]): - return + return with open(filename, 'r') as file_dat: @@ -43,11 +43,11 @@ def test_double_precision(filename): assert re.search(double_prec, l.split('!')[0]) is None except UnicodeDecodeError: - return + return def test_dexp(filename): if any([f'/{s}/' in str(filename) for s in ignore_dirs]): - return + return with open(filename, 'r') as file_dat: @@ -59,11 +59,11 @@ def test_dexp(filename): assert re.search(dexp, l.split('!')[0]) is None except UnicodeDecodeError: - return + return def test_dlog(filename): if any([f'/{s}/' in str(filename) for s in ignore_dirs]): - return + return with open(filename, 'r') as file_dat: @@ -77,7 +77,7 @@ def test_dlog(filename): assert re.search(dlog10, l.split('!')[0]) is None except UnicodeDecodeError: - return + return @@ -85,9 +85,9 @@ def test_check_rt(filename): """ make sure that all of the numerical constants use _rt and are defined as real(rt) """ - + if any([f'/{s}/' in str(filename) for s in ignore_dirs]): - return + return elif 'extern.F90' in str(filename): return @@ -106,4 +106,4 @@ def test_check_rt(filename): # assert 'use amrex_constants_module' not in l.split('!')[0] except UnicodeDecodeError: - return \ No newline at end of file + return \ No newline at end of file