From 05cbbda50b2dc544ebd65897d6d9760ba9cfe5ac Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Fri, 15 Sep 2023 15:51:57 -0400 Subject: [PATCH] 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; } } }