diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H index e3dfeba212..8510ed6de4 100644 --- a/nse_solver/nse_check.H +++ b/nse_solver/nse_check.H @@ -137,8 +137,8 @@ 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){ +void get_rate_by_nuc(int& rate_index, amrex::Array1D reactants_indices, + amrex::Array1D products_indices){ // fill in the corresponding rate_index by providing array of reactants and products indices @@ -148,7 +148,7 @@ void get_rate_by_nuc(int& rate_index, amrex::Array1D reactants_indice // if all -1 then return, since its an invalid rate - for (int n = 1; n <= 6; ++n){ + for (int n = 1; n <= 8; ++n){ if (NSE_INDEX::rate_indices(rate_index, n) != -1){ valid_rate = true; break; @@ -163,14 +163,14 @@ void get_rate_by_nuc(int& rate_index, amrex::Array1D reactants_indice // sort indices just to make sure - sort_Array1D<1, 3>(reactants_indices); - sort_Array1D<1, 3>(products_indices); + sort_Array1D<1, 4>(reactants_indices); + sort_Array1D<1, 4>(products_indices); // check for which rate matches these input indices 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)){ + for (int j = 1; j <= 4; ++j){ + if (NSE_INDEX::rate_indices(i, j) == reactants_indices(j) && NSE_INDEX::rate_indices(i, j+4) == products_indices(j)){ found_index = true; } else{ @@ -210,8 +210,8 @@ void find_fast_reaction_cycle(const int current_nuc_ind, const amrex::Array1D reaction_scratch; - amrex::Array1D product_scratch; + amrex::Array1D reaction_scratch; + amrex::Array1D product_scratch; int rate_index; // store index of the rate bool have_required_nucs; // variable used to check we have required nuc @@ -374,11 +374,13 @@ void find_fast_reaction_cycle(const int current_nuc_ind, const amrex::Array1D& merge_indices, // first 3 is 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}; + 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}; // Here Y_nonNPA represent 2 group isotopes that are not n,p,a // there can be at most 2 nonNPA isotopes in a reaction for it to be valid @@ -760,7 +770,7 @@ void is_fastest_rate(amrex::Array1D& merge_indices, // Check if there are > 2 isotopes or 0 isotopes in LIG (light isotope group) in this reaction - for (int k = 1; k <= 6; ++k){ + for (int k = 1; k <= 8; ++k){ // skip if current isotope is -1, which means null @@ -837,8 +847,8 @@ void is_fastest_rate(amrex::Array1D& merge_indices, // 3) number of nonLIG isotopes are more than 2 or none. if ( (pair_rate_index == -1) - || (NSE_INDEX::rate_indices(current_rate, 1) != -1) - || (NSE_INDEX::rate_indices(current_rate, 4) != -1) + || (NSE_INDEX::rate_indices(current_rate, 2) != -1) + || (NSE_INDEX::rate_indices(current_rate, 6) != -1) || (num_nonLIG > 2) || (num_nonLIG == 0) ){ @@ -850,25 +860,25 @@ void is_fastest_rate(amrex::Array1D& merge_indices, amrex::Real b_f; amrex::Real b_r; - b_f = screened_rates(current_rate) * Y(NSE_INDEX::rate_indices(current_rate, 3) + 1); - b_r = screened_rates(pair_rate_index) * Y(NSE_INDEX::rate_indices(current_rate, 6) + 1); + b_f = screened_rates(current_rate) * Y(NSE_INDEX::rate_indices(current_rate, 4) + 1); + b_r = screened_rates(pair_rate_index) * Y(NSE_INDEX::rate_indices(current_rate, 8) + 1); - if (NSE_INDEX::rate_indices(current_rate, 2) != -1){ + if (NSE_INDEX::rate_indices(current_rate, 3) != -1){ - if (NSE_INDEX::rate_indices(current_rate, 2) == NSE_INDEX::rate_indices(current_rate, 3)){ + if (NSE_INDEX::rate_indices(current_rate, 3) == NSE_INDEX::rate_indices(current_rate, 4)){ b_f *= 0.5_rt; } - b_f *= Y(NSE_INDEX::rate_indices(current_rate, 2) + 1) * state.rho; + b_f *= Y(NSE_INDEX::rate_indices(current_rate, 3) + 1) * state.rho; } - if (NSE_INDEX::rate_indices(current_rate, 5) != -1){ + if (NSE_INDEX::rate_indices(current_rate, 7) != -1){ - if (NSE_INDEX::rate_indices(current_rate, 5) == NSE_INDEX::rate_indices(current_rate, 6)){ + if (NSE_INDEX::rate_indices(current_rate, 7) == NSE_INDEX::rate_indices(current_rate, 8)){ b_r *= 0.5_rt; } - b_r *= Y(NSE_INDEX::rate_indices(current_rate, 5) + 1) * state.rho; + b_r *= Y(NSE_INDEX::rate_indices(current_rate, 7) + 1) * state.rho; } // Find the timescale of the rate, See Eq. 17 and Eq. 11 in Kushnir @@ -1002,34 +1012,160 @@ bool in_nse(T& state){ ydot(n) = 0.0_rt; } - rate_t rate_eval; - constexpr int do_T_derivatives = 0; - evaluate_rates(state, rate_eval); +#ifdef NEW_NETWORK_IMPLEMENTATION - eos_t eos_state; // need eos_state for speed of sound - amrex::Real t_s; // a parameter to characterize whether a rate is fast enough + //amrex::Array1D rates + + static_assert(nrhs == neqs || nrhs == 2 * neqs); - // Here we need to find t_s, sound crossing timescale for a single zone + rhs_state_t rhs_state; - if constexpr (std::is_same::value){ + rhs_state.rho = state.rho; + rhs_state.eta = state.eta; + rhs_state.y_e = state.y_e; - // if state is already in eos, then assume eos() is called + // Convert X to Y. + for (int n = 1; n <= NumSpec; ++n) { + rhs_state.y(n) = state.xn[n-1] * aion_inv[n-1]; + } - // see Section 3 in Kushnir - t_s = state.dx / state.cs; + // Set up the state data, which is the same for all screening factors. + fill_plasma_state(rhs_state.pstate, state.T, state.rho, rhs_state.y); + // Initialize the rate temperature term. + rhs_state.tf = get_tfactors(state.T); + if (use_tables) { + rhs_state.tab.initialize(state.T); } - else{ - // if state is burn_t then convert to eos_t and call eos() + // Initialize the RHS terms. + for (int n = 1; n <= NumSpec + 1; ++n) { + ydot(n) = 0.0; + } - burn_to_eos(state, eos_state); - eos(eos_input_rt, eos_state); - t_s = state.dx / eos_state.cs; + // Count up number of intermediate rates (rates that are used in any other reaction). + constexpr int num_intermediate = num_intermediate_reactions(); + + // We cannot have a zero-sized array, so just set the array size to 1 in that case. + constexpr int intermediate_array_size = num_intermediate > 0 ? num_intermediate : 1; + + // Define forward and reverse (and d/dT) rate arrays. + Array1D intermediate_rates; + + // Fill all intermediate rates first. + constexpr_for<1, Rates::NumRates+1>([&] (auto n) + { + constexpr int rate = n; + + constexpr int index = locate_intermediate_rate_index(rate); + + if constexpr (index >= 1) { + construct_rate(rhs_state, intermediate_rates(index)); + } + }); + + // Loop over all rates, and then loop over all species, and if the + // rate affects that given species, add its contribution to the RHS. + constexpr_for<1, Rates::NumRates+1>([&] (auto n1) + { + constexpr int rate = n1; + + rate_t rates; + + // We only need to compute the rate at this point if it's not intermediate. If it + // is intermediate, retrieve it from the cached array. + + constexpr int index = locate_intermediate_rate_index(rate); + if constexpr (index < 0) { + construct_rate(rhs_state, rates); + } + else { + rates = intermediate_rates(index); + } + + // Locate all intermediate rates needed to augment this reaction. + // To keep the problem bounded we assume that there are no more than + // three intermediate reactions needed. + + rate_t rates1, rates2, rates3; + + fill_additional_rates(intermediate_rates, rates1, rates2, rates3); + + // Perform rate postprocessing, using additional reactions as inputs. + // If there is no postprocessing for this rate, this will be a no-op. + + postprocess_rate(rhs_state, rates, rates1, rates2, rates3); + + constexpr_for<1, NumSpec+1>([&] (auto n2) + { + constexpr int species = n2; + + if constexpr (is_rate_used()) { + constexpr int use_T_derivatives = 0; + auto [forward_term, reverse_term] = rhs_term(burn_state, rates); + + if constexpr (nrhs == 2 * neqs) { + if (forward_term >= 0.0_rt) { + ydot(2 * species - 1) += forward_term; + ydot(2 * species ) += -reverse_term; + } + else { + ydot(2 * species - 1) += reverse_term; + ydot(2 * species ) += -forward_term; + } + } + else { + ydot(species) += forward_term + reverse_term; + } + } + }); + }); + + // Evaluate the neutrino cooling. +#ifdef NEUTRINOS + Real sneut, dsneutdt, dsneutdd, snuda, snudz; + sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, + sneut, dsneutdt, dsneutdd, snuda, snudz); +#else + Real sneut = 0.0; +#endif + + // Compute the energy RHS term. + if constexpr (nrhs == 2 * neqs) { + ydot(2 * net_ienuc - 1) = 0.0; + ydot(2 * net_ienuc ) = sneut; } + else { + ydot(net_ienuc) = -sneut; + } + + constexpr_for<1, NumSpec+1>([&] (auto n) + { + constexpr int species = n; + if constexpr (nrhs == 2 * neqs) { + ydot(2 * net_ienuc - 1) += ener_gener_rate(ydot(2 * species - 1)); + ydot(2 * net_ienuc ) += ener_gener_rate(ydot(2 * species )); + } + else { + ydot(net_ienuc) += ener_gener_rate(ydot(species)); + } + }); + // Find ydot + + actual_rhs(state, ydot); + +#else + // Find rates + + rate_t rate_eval; + constexpr int do_T_derivatives = 0; + evaluate_rates(state, rate_eval); + + // Find ydot + rhs_nuc(state, ydot, Y, rate_eval.screened_rates); // Find energy generation rate @@ -1041,13 +1177,36 @@ bool in_nse(T& state){ amrex::Real sneut, dsneutdt, dsneutdd, snuda, snudz; sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); #else - amrex::Real sneut = 0.0_rt, dsneutdt =0.0_rt , dsneutdd = 0.0_rt, snuda = 0.0_rt, snudz = 0.0_rt; + amrex::Real sneut = 0.0_rt; #endif // fill in energy generation rate to ydot ydot(NumSpec + 1) = enuc - sneut; +#endif + + eos_t eos_state; // need eos_state for speed of sound + amrex::Real t_s; // a parameter to characterize whether a rate is fast enough + + // Here we need to find t_s, sound crossing timescale for a single zone + if constexpr (std::is_same::value){ + + // if state is already in eos, then assume eos() is called + + // see Section 3 in Kushnir + t_s = state.dx / state.cs; + + } + else{ + + // if state is burn_t then convert to eos_t and call eos() + + burn_to_eos(state, eos_state); + eos(eos_input_rt, eos_state); + t_s = state.dx / eos_state.cs; + } + // Now we look through the network and see if there are fast reaction cycles // Need to separate forward and reverse rate and detemine each step is fast enough. // use vectors for now diff --git a/nse_solver/nse_index.H b/nse_solver/nse_index.H index b182cd7282..f89c28d132 100644 --- a/nse_solver/nse_index.H +++ b/nse_solver/nse_index.H @@ -18,7 +18,7 @@ namespace NSE_INDEX extern AMREX_GPU_MANAGED int he4_index; // array for storing indices to access different rates - extern AMREX_GPU_MANAGED amrex::Array2D rate_indices; + extern AMREX_GPU_MANAGED amrex::Array2D rate_indices; } #endif diff --git a/nse_solver/nse_index.cpp b/nse_solver/nse_index.cpp index 19efe39038..ad4cfc321d 100644 --- a/nse_solver/nse_index.cpp +++ b/nse_solver/nse_index.cpp @@ -7,4 +7,4 @@ AMREX_GPU_MANAGED int NSE_INDEX::p_index {-1}; AMREX_GPU_MANAGED int NSE_INDEX::h1_index {-1}; AMREX_GPU_MANAGED int NSE_INDEX::n_index {-1}; AMREX_GPU_MANAGED int NSE_INDEX::he4_index {-1}; -AMREX_GPU_MANAGED amrex::Array2D NSE_INDEX::rate_indices; +AMREX_GPU_MANAGED amrex::Array2D NSE_INDEX::rate_indices; diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H index 5af83309c6..3509c6ec91 100644 --- a/nse_solver/nse_solver.H +++ b/nse_solver/nse_solver.H @@ -16,15 +16,18 @@ #include #include #include +#include + +#ifndef NEW_NETWORK_IMPLEMENTATION AMREX_INLINE -void find_rate_indices(amrex::Array1D& reactant_indices, - amrex::Array1D& product_indices, const std::string& rate_name){ +void find_rate_indices(amrex::Array1D& reactant_indices, + amrex::Array1D& product_indices, const std::string& rate_name){ // fill the corresponding reactant and product indices of a given rate name // used during initialization - for (int n = 1; n <= 3; ++n){ + for (int n = 1; n <= 4; ++n){ reactant_indices(n) = -1; product_indices(n) = -1; } @@ -71,14 +74,12 @@ void find_rate_indices(amrex::Array1D& reactant_indices, return; } - std::string _temp; - // Find the index of reactants and products in short_spec_names int i = 1; while ((pos = reactant_names.find('_')) != std::string::npos){ - _temp = reactant_names.substr(0, pos); + std::string _temp = reactant_names.substr(0, pos); // aprox network need capitalize // but we only work with pynucastro networks for now. @@ -105,7 +106,7 @@ void find_rate_indices(amrex::Array1D& reactant_indices, i = 1; while ((pos = product_names.find('_')) != std::string::npos){ - _temp = product_names.substr(0, pos); + std::string _temp = product_names.substr(0, pos); // cap_first_letter(product); for (int n = 0; n < NumSpec; ++n){ @@ -128,11 +129,11 @@ void find_rate_indices(amrex::Array1D& reactant_indices, // Sort rate_indices. - sort_Array1D<1, 3>(reactant_indices); - sort_Array1D<1, 3>(product_indices); + sort_Array1D<1, 4>(reactant_indices); + sort_Array1D<1, 4>(product_indices); } - +#endif AMREX_INLINE void init_nse_net() { @@ -145,10 +146,10 @@ void init_nse_net() { #ifdef NEW_NETWORK_IMPLEMENTATION if (zion[n] == 1 && aion[n] == 1){ if (n == 0){ - NSE_INDEX::h1_index = n; + NSE_INDEX::h1_index = n; } else { - NSE_INDEX::p_index = n; + NSE_INDEX::p_index = n; // photodisintegration proton } } #else @@ -183,28 +184,87 @@ void init_nse_net() { } } -#ifndef NEW_NETWORK_IMPLEMENTATION +#ifdef NEW_NETWORK_IMPLEMENTATION + + // Fill in rate_indices for aprox networks + + for (int n = 1; n <= Rates::NumRates; ++n){ + + RHS::rhs_t data = RHS::rhs_data(n); + + amrex::Array1D reactant_indices{-1, -1, -1 ,-1}; + amrex::Array1D product_indices{-1, -1, -1, -1}; + + int m = 1; + + for (int i = 0; i < data.number_A; ++i){ + reactant_indices(m) = data.species_A - 1; + m += 1; + } + + for (int i = 0; i < data.number_B; ++i){ + reactant_indices(m) = data.species_B - 1; + m += 1; + } + + for (int i = 0; i < data.number_C; ++i){ + reactant_indices(m) = data.species_C - 1; + m += 1; + } + + m = 1; + + for (int i = 0; i < data.number_D; ++i){ + product_indices(m) = data.species_D - 1; + m += 1; + } + + for (int i = 0; i < data.number_E; ++i){ + product_indices(m) = data.species_E - 1; + m += 1; + } + + for (int i = 0; i < data.number_F; ++i){ + product_indices(m) = data.species_F - 1; + m += 1; + } + + // Sort rate_indices. + + sort_Array1D<1, 4>(reactant_indices); + sort_Array1D<1, 4>(product_indices); + + // fill in rate_indices + + for (int i = 1; i <= 4; ++i){ + NSE_INDEX::rate_indices(n, i) = reactant_indices(i); + NSE_INDEX::rate_indices(n, i+4) = product_indices(i); + } + + } + +#else + // rate names initialization for non aprox networks... - // Fill in rate_indices + // Fill in rate_indices // used for referencing rate by nuclei indices - amrex::Array1D reactant_indices; - amrex::Array1D product_indices; - std::string current_rate_name; + amrex::Array1D reactant_indices; + amrex::Array1D product_indices; for (int n = 1; n <= Rates::NumRates; ++n){ // get current name of the rate - current_rate_name = Rates::rate_names[n]; + std::string current_rate_name = Rates::rate_names[n]; // fill in reactant_indices and product indices find_rate_indices(reactant_indices, product_indices, current_rate_name); // fill in rate_indices - for (int i = 1; i <= 3; ++i){ + for (int i = 1; i <= 4; ++i){ NSE_INDEX::rate_indices(n, i) = reactant_indices(i); - NSE_INDEX::rate_indices(n, i+3) = product_indices(i); + NSE_INDEX::rate_indices(n, i+4) = product_indices(i); } }