From 04ff2a3e3578134387c2cf2bb740b6823e967f76 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Mon, 10 Jul 2023 13:31:54 -0400 Subject: [PATCH 1/2] update adsurl for refs.bib (#1270) --- sphinx_docs/source/refs.bib | 110 ++++++++++++++++++------------------ 1 file changed, 54 insertions(+), 56 deletions(-) diff --git a/sphinx_docs/source/refs.bib b/sphinx_docs/source/refs.bib index 015e4d2b08..d10c1b5398 100644 --- a/sphinx_docs/source/refs.bib +++ b/sphinx_docs/source/refs.bib @@ -19,7 +19,7 @@ @ARTICLE{colglaz volume = 59, pages = {264-289}, doi = {10.1016/0021-9991(85)90146-9}, - adsurl = {http://adsabs.harvard.edu/abs/1985JCoPh..59..264C}, + adsurl = {http://ui.adsabs.harvard.edu/abs/1985JCoPh..59..264C}, } @ARTICLE{castro_I, @@ -118,7 +118,7 @@ @ARTICLE{timmes:1999 volume = 125, pages = {277-294}, doi = {10.1086/313271}, - adsurl = {http://adsabs.harvard.edu/abs/1999ApJS..125..277T}, + adsurl = {http://ui.adsabs.harvard.edu/abs/1999ApJS..125..277T}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -132,7 +132,7 @@ @ARTICLE{timmes:2000 volume = 126, pages = {501-516}, doi = {10.1086/313304}, - adsurl = {http://adsabs.harvard.edu/abs/2000ApJS..126..501T}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2000ApJS..126..501T}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -168,7 +168,7 @@ @ARTICLE{millercolella:2002 volume = 183, pages = {26-82}, doi = {10.1006/jcph.2002.7158}, - adsurl = {http://adsabs.harvard.edu/abs/2002JCoPh.183...26M}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2002JCoPh.183...26M}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -192,7 +192,7 @@ @ARTICLE{chamulak:2008 volume = 677, pages = {160-168}, doi = {10.1086/528944}, - adsurl = {http://adsabs.harvard.edu/abs/2008ApJ...677..160C}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2008ApJ...677..160C}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -208,7 +208,7 @@ @ARTICLE{lowMach4 volume = 704, pages = {196-210}, doi = {10.1088/0004-637X/704/1/196}, - adsurl = {http://adsabs.harvard.edu/abs/2009ApJ...704..196Z}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2009ApJ...704..196Z}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -224,7 +224,7 @@ @ARTICLE{wdconvect eid = {8}, pages = {8}, doi = {10.1088/0004-637X/740/1/8}, - adsurl = {http://adsabs.harvard.edu/abs/2011ApJ...740....8Z}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2011ApJ...740....8Z}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -242,7 +242,7 @@ @ARTICLE{raskin:2010 volume = 724, pages = {111-125}, doi = {10.1088/0004-637X/724/1/111}, - adsurl = {http://adsabs.harvard.edu/abs/2010ApJ...724..111R}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2010ApJ...724..111R}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -261,7 +261,7 @@ @ARTICLE{xrb:III eid = {60}, pages = {60}, doi = {10.1088/0004-637X/807/1/60}, - adsurl = {http://adsabs.harvard.edu/abs/2015ApJ...807...60Z}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2015ApJ...807...60Z}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -280,26 +280,22 @@ @ARTICLE{xrb:II eid = {115}, pages = {115}, doi = {10.1088/0004-637X/788/2/115}, - adsurl = {http://adsabs.harvard.edu/abs/2014ApJ...788..115M}, + adsurl = {http://ui.adsabs.harvard.edu/abs/2014ApJ...788..115M}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @ARTICLE{wallacewoosley:1981, - author = {{Wallace}, R.~K. and {Woosley}, S.~E.}, - title = "{Explosive hydrogen burning}", - journal = {\apjs}, - keywords = {ABUNDANCE, ASTROPHYSICS, HYDROGEN, METALLIC STARS, - NUCLEAR FUSION, REACTION KINETICS, STELLAR - EVOLUTION, GAMMA RAYS, HIGH TEMPERATURE, ISOTOPES, - NEUTRON STARS, NOVAE, RESONANCE, SUPERMASSIVE STARS, - X RAY SOURCES}, - year = 1981, - month = feb, - volume = 45, - pages = {389-420}, - doi = {10.1086/190717}, - adsurl = {http://adsabs.harvard.edu/abs/1981ApJS...45..389W}, - adsnote = {Provided by the SAO/NASA Astrophysics Data System} + author = {{Wallace}, R.~K. and {Woosley}, S.~E.}, + title = "{Explosive hydrogen burning}", + journal = {\apjs}, + keywords = {Abundance, Astrophysics, Hydrogen, Metallic Stars, Nuclear Fusion, Reaction Kinetics, Stellar Evolution, Gamma Rays, High Temperature, Isotopes, Neutron Stars, Novae, Resonance, Supermassive Stars, X Ray Sources, Astrophysics}, + year = 1981, + month = feb, + volume = {45}, + pages = {389-420}, + doi = {10.1086/190717}, + adsurl = {https://ui.adsabs.harvard.edu/abs/1981ApJS...45..389W}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @article{ReacLib, @@ -324,7 +320,7 @@ @BOOK{NR title = "{Numerical recipes in FORTRAN. The art of scientific computing}", publisher = {Cambridge: University Press, |c1992, 2nd ed.}, year = 1992, - adsurl = {http://adsabs.harvard.edu/abs/1992nrfa.book.....P}, + adsurl = {http://ui.adsabs.harvard.edu/abs/1992nrfa.book.....P}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -389,43 +385,45 @@ @ARTICLE{ma:2013 adsnote = {Provided by the SAO/NASA Astrophysics Data System} } - @ARTICLE{graboske:1973, - author = {{Graboske}, H.~C. and {Dewitt}, H.~E. and {Grossman}, A.~S. and {Cooper}, M.~S.}, - title = "{Screening Factors for Nuclear Reactions. 11. Intermediate Screen-Ing and Astrophysical Applications}", - journal = {\apj}, - year = 1973, - month = apr, - volume = 181, - pages = {457-474}, - adsurl = {http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1973ApJ...181..457G&db_key=AST}, - adsnote = {Provided by the Smithsonian/NASA Astrophysics Data System} + author = {{Graboske}, H.~C. and {Dewitt}, H.~E. and {Grossman}, A.~S. and {Cooper}, M.~S.}, + title = "{Screening Factors for Nuclear Reactions. II. Intermediate Screen-Ing and Astrophysical Applications}", + journal = {\apj}, + year = 1973, + month = apr, + volume = {181}, + pages = {457-474}, + doi = {10.1086/152062}, + adsurl = {https://ui.adsabs.harvard.edu/abs/1973ApJ...181..457G}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @ARTICLE{alastuey:1978, - author = {{Alastuey}, A. and {Jancovici}, B.}, - title = "{Nuclear reaction rate enhancement in dense stellar matter}", - journal = {\apj}, - year = 1978, - month = dec, - volume = 226, - pages = {1034-1040}, - doi = {10.1086/156681}, - adsurl = {http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1978ApJ...226.1034A&db_key=AST}, - adsnote = {Provided by the Smithsonian/NASA Astrophysics Data System} + author = {{Alastuey}, A. and {Jancovici}, B.}, + title = "{Nuclear reaction rate enhancement in dense stellar matter.}", + journal = {\apj}, + keywords = {Astrophysics, Magnetohydrodynamics, Nuclear Astrophysics, Reaction Kinetics, Thermonuclear Reactions, Correlation, Nuclei (Nuclear Physics), Perturbation Theory, Potential Theory, Quantum Mechanics, Astrophysics, Nuclear Reactions:Stellar Interiors}, + year = 1978, + month = dec, + volume = {226}, + pages = {1034-1040}, + doi = {10.1086/156681}, + adsurl = {https://ui.adsabs.harvard.edu/abs/1978ApJ...226.1034A}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @ARTICLE{itoh:1979, - author = {{Itoh}, N. and {Totsuji}, H. and {Ichimaru}, S. and {Dewitt}, H.~E.}, - title = "{Enhancement of thermonuclear reaction rate due to strong screening. II - Ionic mixtures}", - journal = {\apj}, - year = 1979, - month = dec, - volume = 234, - pages = {1079-1084}, - doi = {10.1086/157590}, - adsurl = {http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1979ApJ...234.1079I&db_key=AST}, - adsnote = {Provided by the Smithsonian/NASA Astrophysics Data System} + author = {{Itoh}, N. and {Totsuji}, H. and {Ichimaru}, S. and {Dewitt}, H.~E.}, + title = "{Enhancement of thermonuclear reaction rate due to strong screening. II - Ionic mixtures}", + journal = {\apj}, + keywords = {Nonuniform Plasmas, Reaction Kinetics, Thermonuclear Reactions, Astrophysics, Binary Mixtures, Dense Plasmas, Monte Carlo Method, Astrophysics}, + year = 1979, + month = dec, + volume = {234}, + pages = {1079-1084}, + doi = {10.1086/157590}, + adsurl = {https://ui.adsabs.harvard.edu/abs/1979ApJ...234.1079I}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @ARTICLE{chugunov:2007, From ffe836cf6a5f6ce67e522cfaba9d1ae1377ca294 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Mon, 10 Jul 2023 14:24:09 -0400 Subject: [PATCH 2/2] make self-consistent nse compatible with non-neutron networks. (#1231) Due to the missing neutron rates, it is difficult to group isotopes together into a single group. This pr mainly adds another logic of defining a single_group for non-neutron networks. If isotopes fail to form a single group and there is no neutron in the network, then use a second check where we consider nse if isotopes past Z=14(Si) that have odd N or even N form two distinct groups. Most of the line changes are code formatting. --- nse_solver/nse_check.H | 312 ++++++++++++++++++++++--------------- nse_solver/nse_solver.H | 117 +++++++------- sphinx_docs/source/nse.rst | 25 +-- 3 files changed, 262 insertions(+), 192 deletions(-) diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H index eaea987bd9..ad7d3ca888 100644 --- a/nse_solver/nse_check.H +++ b/nse_solver/nse_check.H @@ -21,25 +21,23 @@ // First check to see if we're in the ballpark of nse state AMREX_GPU_HOST_DEVICE AMREX_INLINE -void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check){ +void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) { // This function gives the first estimate whether we're in the nse or not // it checks whether the molar fractions of n,p,a are approximately in NSE - int count = 0; amrex::Real r = 1.0_rt; amrex::Real r_nse = 1.0_rt; - // Check if neutron, proton, and helium are present in the network - for (int n = 0; n < NumSpec; ++n){ - if (n == NSE_INDEX::p_index || n == NSE_INDEX::n_index|| n == NSE_INDEX::he4_index){ - ++count; - } - } + nse_check = false; - if (count < 3){ - amrex::Error("Current network does not include the proton, neutron, or helium-4, thus cannot proceed with ASE."); + // raise error if no proton or helium-4 in the network. + + if (NSE_INDEX::p_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 // these two ratios are defined in the ASE paper to determine whether network is in equilibrium @@ -60,24 +58,29 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check){ // equilibrium condition: if pass proceed with ase if not proceed with regular eos integration // Eq. 14 in Kushnir paper - if (std::abs((r - r_nse) / r_nse) < 0.5_rt){ + // if there is neutron in the network + + if ((std::abs((r - r_nse) / r_nse) < 0.5_rt) && (NSE_INDEX::n_index != -1)) { nse_check = true; } - else{ - nse_check = false; + + // if there is no neutron in the network + + if ((std::abs((r - r_nse) / r_nse) < 0.25_rt) && (NSE_INDEX::n_index == -1)) { + nse_check = true; } } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void find_max_nucs(int& max_nucs, const int current_nuc_ind){ +void find_max_nucs(int& max_nucs, const int current_nuc_ind) { // Find max number of nucs possible in a reaction step given a nuclei index. // used to determine the size of reaction array max_nucs = 0; - for (int n = NumSpec-1; n > NSE_INDEX::he4_index; --n){ + 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)) @@ -86,7 +89,7 @@ void find_max_nucs(int& max_nucs, const int current_nuc_ind){ (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; } @@ -102,7 +105,7 @@ void find_max_nucs(int& max_nucs, const int current_nuc_ind){ (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; } @@ -117,7 +120,7 @@ void find_max_nucs(int& max_nucs, const int current_nuc_ind){ || (aion[n] == aion[current_nuc_ind]-3 && (zion[n] == zion[current_nuc_ind]-1 || zion[n] == zion[current_nuc_ind]-2)) - ){ + ) { ++max_nucs; } @@ -127,7 +130,7 @@ void find_max_nucs(int& max_nucs, const int current_nuc_ind){ (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; } @@ -137,7 +140,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; @@ -150,18 +153,18 @@ void get_rate_by_nuc(int& rate_index, amrex::Array1D reactants_indice // check for which rate matches these input indices // fill in the corresponding rate_index by providing array of reactants and products 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 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; } - else{ + else { found_index = false; break; } } - if (found_index){ + if (found_index) { rate_index = i; return; } @@ -173,7 +176,7 @@ 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 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 @@ -201,8 +204,8 @@ void find_fast_reaction_cycle(const int current_nuc_ind, const amrex::Array1D NSE_INDEX::he4_index; --n){ + for (int n = NumSpec-1; n > NSE_INDEX::he4_index; --n) { if ( (aion[n] == aion[current_nuc_ind]-1 && @@ -231,7 +234,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 @@ -595,14 +598,14 @@ int get_root_index(const int nuc_ind, int root_index; int scratch_ind = nuc_ind; - while(true){ + while(true) { root_index = group_ind(scratch_ind + 1); - if (root_index != scratch_ind + 1){ + if (root_index != scratch_ind + 1) { scratch_ind = root_index - 1; } - else{ + else { return root_index; } } @@ -611,7 +614,7 @@ int get_root_index(const int nuc_ind, AMREX_GPU_HOST_DEVICE AMREX_INLINE -void nse_union(const int nuc_ind_a, const int nuc_ind_b, amrex::Array1D& group_ind){ +void nse_union(const int nuc_ind_a, const int nuc_ind_b, amrex::Array1D& group_ind) { // This function joins the two group of the two nuc indices:nuc_ind_a and nuc_ind_b // The smaller group is joined to the larger group. @@ -619,7 +622,7 @@ void nse_union(const int nuc_ind_a, const int nuc_ind_b, amrex::Array1D= group_b_size){ + if (group_a_size >= group_b_size) { group_ind(root_index_b) = group_ind(root_index_a); } - - else{ + else { group_ind(root_index_a) = group_ind(root_index_b); } } AMREX_GPU_HOST_DEVICE AMREX_INLINE -int get_pair_rate_index(const int& rate_index){ +int get_pair_rate_index(const int& rate_index) { // This function finds the pair rate index given a rate. // For example given forward rate index, return reverse rate index, vice versa @@ -662,14 +664,14 @@ int get_pair_rate_index(const int& rate_index){ // if all -1 then return, since its an invalid rate - for (int n = 1; n <= 6; ++n){ - if (NSE_INDEX::rate_indices(rate_index, n) != -1){ + for (int n = 1; n <= 6; ++n) { + if (NSE_INDEX::rate_indices(rate_index, n) != -1) { valid_rate = true; break; } } - if (!valid_rate){ + if (!valid_rate) { return pair_rate_index; } @@ -698,31 +700,76 @@ int get_pair_rate_index(const int& rate_index){ AMREX_GPU_HOST_DEVICE AMREX_INLINE -bool in_single_group(const amrex::Array1D& group_ind){ +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 LIG_root_index = get_root_index(NSE_INDEX::he4_index, group_ind); + + int oddN_group = -1; + int evenN_group = -1; - for (int n = 0; n < NumSpec; ++n){ + bool in_single_group = true; - if (get_root_index(n, group_ind) == LIG_root_index){ + // Consider NSE when there is a single group with an optional LIG group + + for (int n = 0; n < NumSpec; ++n) { + if (get_root_index(n, group_ind) == LIG_root_index) { continue; } - if (nonLIG_index == -1){ + if (nonLIG_index == -1) { nonLIG_index = get_root_index(n, group_ind); + continue; } - else if (get_root_index(n, group_ind) != nonLIG_index){ - return false; + if (get_root_index(n, group_ind) != nonLIG_index) { + in_single_group = false; + break; } + + } + + // If there no neutrons are in the network and original condition failed + // Consider a looser condition by looking at nuclei heavier than Si28 + // 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; + for (int n = 0; n < NumSpec; ++n) { + + 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; + } + } + + } } - return true; + return in_single_group; } @@ -756,7 +803,10 @@ void is_fastest_rate(amrex::Array1D& merge_indices, // there can be at most 2 nonNPA isotopes in a reaction for it to be valid amrex::Array1D nonNPA_ind = {-1, -1}; - int m = 1; // m is counts the number of nonNPA isotope detected. + + // m is counts the number of nonNPA isotope detected. + + int m = 1; // find the reverse rate index of the current rate @@ -764,11 +814,11 @@ 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 <= 6; ++k) { // skip if current isotope is -1, which means null - if (NSE_INDEX::rate_indices(current_rate, k) == -1){ + if (NSE_INDEX::rate_indices(current_rate, k) == -1) { continue; } @@ -779,11 +829,11 @@ void is_fastest_rate(amrex::Array1D& merge_indices, // Determine number of nonLIG isotopes // also check whether nonLIG isotopes are already merged - if (root_index != get_root_index(NSE_INDEX::he4_index, group_ind)){ + if (root_index != get_root_index(NSE_INDEX::he4_index, group_ind)) { ++num_nonLIG; - if (root_index != nonLIG_root){ + if (root_index != nonLIG_root) { nonLIG_root = root_index; } @@ -798,17 +848,18 @@ void is_fastest_rate(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)){ + if (root_index == get_root_index(NSE_INDEX::he4_index, group_ind)) { Y_group(k) = Y(NSE_INDEX::rate_indices(current_rate, k) + 1); } - else{ - for (int n = 0; n < NumSpec; ++n){ + else { + + for (int n = 0; n < NumSpec; ++n) { // if nuclei have the same root_index, then they are the same group - if (get_root_index(n, group_ind) == root_index){ + if (get_root_index(n, group_ind) == root_index) { Y_group(k) += Y(n + 1); } @@ -820,11 +871,11 @@ void is_fastest_rate(amrex::Array1D& merge_indices, if (NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::p_index && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::n_index && NSE_INDEX::rate_indices(current_rate, k) != NSE_INDEX::he4_index - ){ + ) { // return if there are more than 2 nonNPA nuclei involved in the network - if (m > 2){ + if (m > 2) { return; } @@ -845,7 +896,7 @@ void is_fastest_rate(amrex::Array1D& merge_indices, || (NSE_INDEX::rate_indices(current_rate, 4) != -1) || (num_nonLIG > 2) || (num_nonLIG == 0) - ){ + ) { return; } @@ -857,18 +908,18 @@ void is_fastest_rate(amrex::Array1D& merge_indices, 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); - if (NSE_INDEX::rate_indices(current_rate, 2) != -1){ + if (NSE_INDEX::rate_indices(current_rate, 2) != -1) { - if (NSE_INDEX::rate_indices(current_rate, 2) == NSE_INDEX::rate_indices(current_rate, 3)){ + if (NSE_INDEX::rate_indices(current_rate, 2) == NSE_INDEX::rate_indices(current_rate, 3)) { b_f *= 0.5_rt; } b_f *= Y(NSE_INDEX::rate_indices(current_rate, 2) + 1) * state.rho; } - if (NSE_INDEX::rate_indices(current_rate, 5) != -1){ + if (NSE_INDEX::rate_indices(current_rate, 5) != -1) { - if (NSE_INDEX::rate_indices(current_rate, 5) == NSE_INDEX::rate_indices(current_rate, 6)){ + if (NSE_INDEX::rate_indices(current_rate, 5) == NSE_INDEX::rate_indices(current_rate, 6)) { b_r *= 0.5_rt; } @@ -878,13 +929,13 @@ void is_fastest_rate(amrex::Array1D& merge_indices, // Find the timescale of the rate, See Eq. 17 and Eq. 11 in Kushnir amrex::Real scratch_t = Y_group(nonNPA_ind(1)) / amrex::min(b_f, b_r); - if (nonNPA_ind(2) != -1){ + if (nonNPA_ind(2) != -1) { scratch_t = amrex::min(scratch_t, Y_group(nonNPA_ind(2)) / amrex::min(b_f, b_r)); } // return if current rate timescale is larger than previous time scale - if (fastest_t < scratch_t){ + if (fastest_t < scratch_t) { return; } @@ -899,14 +950,14 @@ void is_fastest_rate(amrex::Array1D& merge_indices, // get merge indices for nonNPA isotopes - for (int n = 1; n <= 2; ++n){ + for (int n = 1; n <= 2; ++n) { // If nonNPA index is -1, meaning null, then use LIG index - if (nonNPA_ind(n) == -1){ + if (nonNPA_ind(n) == -1) { merge_indices(n) = get_root_index(NSE_INDEX::he4_index, group_ind); } - else{ + else { merge_indices(n) = NSE_INDEX::rate_indices(current_rate, nonNPA_ind(n)); } } @@ -920,20 +971,23 @@ 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::Real& t_s) { // This function groups all the nuclei using group_ind // which contains the node # // fill in initial group_ind, group_ind go from 1 to NumSpec - for (int i = 1; i <= NumSpec; ++i){ + for (int i = 1; i <= NumSpec; ++i) { group_ind(i) = i; // let n,p,a form the same group (LIG) initially, let 1 be index of LIG if (i == NSE_INDEX::p_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; + group_ind(i) = 1; } } @@ -943,7 +997,7 @@ void nse_grouping(amrex::Array1D& group_ind, const burn_t& stat // Go through each reaction and perform grouping - while(true){ + while(true) { merge_indices(1) = -1; merge_indices(2) = -1; @@ -951,7 +1005,7 @@ void nse_grouping(amrex::Array1D& group_ind, const burn_t& stat // Find the fastest time scale of the available reaction - for (int n = 1; n <= Rates::NumRates; ++n){ + for (int n = 1; n <= Rates::NumRates; ++n) { // Check if current rate is the faster rate than previous rate @@ -960,7 +1014,7 @@ void nse_grouping(amrex::Array1D& group_ind, const burn_t& stat // if there are no reaction satisfy conditions, break - if (merge_indices(1) == -1 && merge_indices(2) == -1){ + if (merge_indices(1) == -1 && merge_indices(2) == -1) { break; } @@ -983,7 +1037,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { current_state.nse = false; - if (current_state.T < 2.0e9_rt) { + if (current_state.T < 2.5e9_rt) { return current_state.nse; } @@ -1016,7 +1070,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { // 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) { + if (nse_molar_independent || skip_molar_check) { state = nse_state; state.dx = current_state.dx; #ifndef STRANG @@ -1032,7 +1086,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { amrex::Array1D Y; amrex::Array1D ydot; - for (int n = 1; n <= NumSpec; ++n){ + for (int n = 1; n <= NumSpec; ++n) { Y(n) = state.xn[n-1] * aion_inv[n-1]; ydot(n) = 0.0_rt; } @@ -1076,27 +1130,33 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { 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; + + bool found_fast_reaction_cycle = true; + + // Additional constraint if there is neutron in the network: // Now we look through the network and see if there are fast reaction cycles // Need to separate forward and reverse rate and determine each step is fast enough. - // use vectors for now - bool found_fast_reaction_cycle = false; + if (NSE_INDEX::n_index != -1) { + found_fast_reaction_cycle = false; - // Do a reverse for loop to start from heaviest nuclei + // Do a reverse for loop to start from heaviest nuclei - for (int n = NumSpec-1; n >= 0; --n){ - if (found_fast_reaction_cycle){ - break; + for (int n = NumSpec-1; n >= 0; --n) { + if (found_fast_reaction_cycle) { + break; + } + + find_fast_reaction_cycle(n, Y, rate_eval.screened_rates, state, t_s, found_fast_reaction_cycle); } - find_fast_reaction_cycle(n, Y, rate_eval.screened_rates, state, t_s, found_fast_reaction_cycle); } // Do nse grouping if found fast reaction cycle. @@ -1108,10 +1168,10 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { // if not found_fast_reaction_cycle, state not in nse - if (!found_fast_reaction_cycle){ + if (!found_fast_reaction_cycle) { current_state.nse = false; } - else{ + else { // Do nse grouping if found_fast_reaction_cycle @@ -1121,7 +1181,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) { current_state.nse = false; - if (in_single_group(group_ind)){ + if (in_single_group(group_ind)) { current_state.nse = true; } } diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H index 87de4a7476..4f4829bf1f 100644 --- a/nse_solver/nse_solver.H +++ b/nse_solver/nse_solver.H @@ -21,12 +21,12 @@ using namespace nse_rp; AMREX_INLINE void find_rate_indices(amrex::Array1D& reactant_indices, - amrex::Array1D& product_indices, const std::string& rate_name){ + 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 <= 3; ++n) { reactant_indices(n) = -1; product_indices(n) = -1; } @@ -39,10 +39,10 @@ void find_rate_indices(amrex::Array1D& reactant_indices, size_t pos; - while (true){ + while (true) { pos = rate_name_mod.find("p_"); - if (pos == std::string::npos){ + if (pos == std::string::npos) { break; } @@ -80,7 +80,7 @@ void find_rate_indices(amrex::Array1D& reactant_indices, // Find the index of reactants and products in short_spec_names int i = 1; - while ((pos = reactant_names.find('_')) != std::string::npos){ + while ((pos = reactant_names.find('_')) != std::string::npos) { _temp = reactant_names.substr(0, pos); // aprox network need capitalize @@ -88,8 +88,8 @@ void find_rate_indices(amrex::Array1D& reactant_indices, // cap_first_letter(reactant); - for (int n = 0; n < NumSpec; ++n){ - if (short_spec_names_cxx[n] == _temp){ + for (int n = 0; n < NumSpec; ++n) { + if (short_spec_names_cxx[n] == _temp) { reactant_indices(i) = n; ++i; } @@ -99,21 +99,20 @@ void find_rate_indices(amrex::Array1D& reactant_indices, } // cap_first_letter(reactant_names); - for (int n = 0; n < NumSpec; ++n){ - if (short_spec_names_cxx[n] == reactant_names){ + for (int n = 0; n < NumSpec; ++n) { + if (short_spec_names_cxx[n] == reactant_names) { reactant_indices(i) = n; } } - i = 1; - while ((pos = product_names.find('_')) != std::string::npos){ + while ((pos = product_names.find('_')) != std::string::npos) { _temp = product_names.substr(0, pos); // cap_first_letter(product); - for (int n = 0; n < NumSpec; ++n){ - if (short_spec_names_cxx[n] == _temp){ + for (int n = 0; n < NumSpec; ++n) { + if (short_spec_names_cxx[n] == _temp) { product_indices(i) = n; ++i; } @@ -124,8 +123,8 @@ void find_rate_indices(amrex::Array1D& reactant_indices, // cap_first_letter(product_names); - for (int n = 0; n < NumSpec; ++n){ - if (short_spec_names_cxx[n] == product_names){ + for (int n = 0; n < NumSpec; ++n) { + if (short_spec_names_cxx[n] == product_names) { product_indices(i) = n; } } @@ -142,13 +141,14 @@ AMREX_INLINE void init_nse_net() { // Initialization for nse, check whether network is valid and stores indices - for (int n = 0; n < NumSpec; ++n){ + for (int n = 0; n < NumSpec; ++n) { // store index of photoionization proton by looking for Z=1 but // exclude h1 which is always at index=0 for aprox networks + #ifdef NEW_NETWORK_IMPLEMENTATION - if (zion[n] == 1 && aion[n] == 1){ - if (n == 0){ + if (zion[n] == 1 && aion[n] == 1) { + if (n == 0) { NSE_INDEX::h1_index = n; } else { @@ -156,14 +156,14 @@ void init_nse_net() { } } #else - if (zion[n] == 1 && aion[n]== 1){ + if (zion[n] == 1 && aion[n]== 1) { NSE_INDEX::p_index = n; } #endif - else if (zion[n] == 0){ + else if (zion[n] == 0) { NSE_INDEX::n_index = n; } - else if (zion[n] == 2 && aion[n] == 4){ + else if (zion[n] == 2 && aion[n] == 4) { NSE_INDEX::he4_index = n; } } @@ -171,18 +171,19 @@ void init_nse_net() { // Check if network results in singular jacobian first, require at // least one nuclei that nuclei.Z != nuclei.N. Some examples include // aprox13 and iso7. If use hybrj solver, this is no longer an issue. - if (!use_hybrid_solver){ + + if (!use_hybrid_solver) { bool singular_network = true; - for (int n = 0; n < NumSpec; ++n){ - if (n == NSE_INDEX::h1_index){ + for (int n = 0; n < NumSpec; ++n) { + if (n == NSE_INDEX::h1_index) { continue; } - if (zion[n] != aion[n] - zion[n]){ + if (zion[n] != aion[n] - zion[n]) { singular_network = false; } } - if (singular_network == true){ + if (singular_network == true) { amrex::Error("This network always results in singular jacobian matrix, thus can't find nse mass fraction using nr!"); } } @@ -197,7 +198,7 @@ void init_nse_net() { amrex::Array1D product_indices; std::string current_rate_name; - for (int n = 1; n <= Rates::NumRates; ++n){ + for (int n = 1; n <= Rates::NumRates; ++n) { // get current name of the rate current_rate_name = Rates::rate_names[n]; @@ -206,19 +207,13 @@ void init_nse_net() { 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 <= 3; ++i) { NSE_INDEX::rate_indices(n, i) = reactant_indices(i); NSE_INDEX::rate_indices(n, i+3) = product_indices(i); } } #endif - - // make sure nse_molar_independent is used if nse_skip_molar is enabled - - if (nse_skip_molar) { - nse_molar_independent = true; - } NSE_INDEX::initialized = true; } @@ -267,8 +262,8 @@ T get_nse_state(const T& state) auto tfactors = evaluate_tfactors(T_in); #endif - for (int n = 0; n < NumSpec; ++n){ - if (n == NSE_INDEX::h1_index){ + for (int n = 0; n < NumSpec; ++n) { + if (n == NSE_INDEX::h1_index) { nse_state.xn[n] = 0.0; continue; } @@ -316,8 +311,8 @@ T get_nse_state(const T& state) nse_state.y_e = 0.0_rt; - for (int n = 0; n < NumSpec; ++n){ - if (n == NSE_INDEX::h1_index){ + for (int n = 0; n < NumSpec; ++n) { + if (n == NSE_INDEX::h1_index) { continue; } @@ -351,8 +346,8 @@ void fcn(Array1D& x, Array1D& fvec, fvec(1) = -1.0_rt; - for (int n = 0; n < NumSpec; ++n){ - if (n == NSE_INDEX::h1_index){ + for (int n = 0; n < NumSpec; ++n) { + if (n == NSE_INDEX::h1_index) { continue; } @@ -392,8 +387,8 @@ void jcn(Array1D& x, Array2D& fjac, Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T; - for (int n = 0; n < NumSpec; ++n){ - if (n == NSE_INDEX::h1_index){ + for (int n = 0; n < NumSpec; ++n) { + if (n == NSE_INDEX::h1_index) { continue; } @@ -439,13 +434,13 @@ void nse_hybrid_solver(T& state, amrex::Real eps=1.0e-10_rt) { // fine tuning initial guesses - for (int i = 0; i < 20; ++i){ + for (int i = 0; i < 20; ++i) { dx = 0.5_rt; inner_x(1) = outer_x(1); inner_x(2) = outer_x(2); - for (int j = 0; j < 20; ++j){ + for (int j = 0; j < 20; ++j) { hj.x(1) = inner_x(1); hj.x(2) = inner_x(2); @@ -455,25 +450,25 @@ void nse_hybrid_solver(T& state, amrex::Real eps=1.0e-10_rt) { 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; } - if (f(1) > 0.0_rt && f(2) > 0.0_rt){ + 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){ + if (is_pos_old != is_pos_new) { dx *= 0.8_rt; } - if (is_pos_new){ + if (is_pos_new) { inner_x(1) -= dx; inner_x(2) -= dx; } @@ -492,6 +487,13 @@ void nse_hybrid_solver(T& state, amrex::Real eps=1.0e-10_rt) { // if (hj.info != 1) { // amrex::Error("failed to solve"); // } +#ifndef AMREX_USE_GPU + std::cout << "NSE solver failed with these conditions: " << std::endl; + std::cout << "Temperature: " << state.T << std::endl; + std::cout << "Density: " << state.rho << std::endl; + std::cout << "Ye: " << state.y_e << std::endl; + std::cout << "Initial mu_p and mu_n: " << state.mu_p << ", " << state.mu_n << std::endl; +#endif amrex::Error("failed to solve"); } @@ -514,7 +516,7 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { x(2) = state.mu_n; jcn(x, jac, state, flag); - fcn(x, f, state,flag); + fcn(x, f, state, flag); amrex::Real det; // store determinant for finding inverse jac amrex::Array2D inverse_jac; // store inverse jacobian @@ -522,12 +524,12 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { amrex::Real d_mu_n = std::numeric_limits::max(); // difference in chemical potential of neutron // begin newton-raphson - for (int i = 0; i < max_nse_iters; ++i){ + for (int i = 0; i < max_nse_iters; ++i) { // check if current state fulfills constraint equation if (std::abs(d_mu_p) < eps * std::abs(x(1)) && - std::abs(d_mu_n) < eps * std::abs(x(2))){ + std::abs(d_mu_n) < eps * std::abs(x(2))) { converged = true; state.mu_p = x(1); state.mu_n = x(2); @@ -540,7 +542,7 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { // if jacobians are small, then no need for scaling - if (scale_fac < 1.0e150){ + if (scale_fac < 1.0e150) { scale_fac = 1.0_rt; } @@ -550,7 +552,7 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { // check if determinant is 0 - if (det == 0.0_rt){ + if (det == 0.0_rt) { amrex::Error("Jacobian is a singular matrix! Try a different initial guess!"); } @@ -568,7 +570,7 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { // if diff goes beyond 1.0e3_rt, likely that its not making good progress.. - if (std::abs(d_mu_p) > 1.0e3_rt || std::abs(d_mu_n) > 1.0e3_rt){ + if (std::abs(d_mu_p) > 1.0e3_rt || std::abs(d_mu_n) > 1.0e3_rt) { amrex::Error("Not making good progress, breaking"); } @@ -579,7 +581,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))){ + if (std::isnan(x(1)) || std::isnan(x(2))) { amrex::Error("Nan encountered, likely due to overflow in digits or not making good progress"); } @@ -589,7 +591,7 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { fcn(x, f, state, flag); } - if (!converged){ + if (!converged) { amrex::Error("NSE solver failed to converge!"); } } @@ -598,7 +600,7 @@ void nse_nr_solver(T& state, amrex::Real eps=1.0e-10_rt) { template AMREX_GPU_HOST_DEVICE AMREX_INLINE T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt, - bool input_ye_is_valid=false){ + bool input_ye_is_valid=false) { // Check whether initialized or not if (!NSE_INDEX::initialized) { @@ -627,6 +629,7 @@ T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt, for (int n = 0; n < NumSpec; ++n) { state.y_e += state.y[SFS+n] * zion[n] * aion_inv[n] / state.rho; } + #endif } diff --git a/sphinx_docs/source/nse.rst b/sphinx_docs/source/nse.rst index e38dc36b50..fbede0efc3 100644 --- a/sphinx_docs/source/nse.rst +++ b/sphinx_docs/source/nse.rst @@ -249,7 +249,7 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. participated in this step, where :math:`b_f` and :math:`b_r` are the forward and reverse rate of the reaction, :math:`\epsilon` is a tolerance which has a default value of :math:`0.1`, and :math:`t_s` is the sound crossing time of a - simulation cell. + simulation cell. Note: we skip this check when there is no neutron involved in the network. An example of such reaction cycle would be: @@ -258,12 +258,12 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. \isotm{S}{32} (\gamma, p)(\gamma, p)(\gamma, n)(\gamma, n) \isotm{Si}{28} (\alpha, \gamma) \isotm{S}{32} - The general approach to this is to start iterations from the heavy to the light nuclei to - use them as the starting point of the cycle. Then the algorithm checks if isotopes involved - in the network can actually form a cycle using the combination reactions above. If such cycle - is formed, then we check the rates of these reactions to see if they satisfy the condition - mention previously. 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 + The general approach to this is to start iterations from the heaviest to the lightest + nuclei. Then the algorithm checks if isotopes involved in the network can actually form + a cycle using the combination reactions above. If such cycle is formed, then we check the + rates of these reactions to see if they satisfy the condition mention previously. + 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, @@ -284,7 +284,7 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. t_{i,k} < \epsilon t_s - * + * .. math:: @@ -312,10 +312,17 @@ There are 3 main criteria discussed in the :cite:`Kushnir_2020`. This means that the non-LIG nuclei have already merged. And the iteration stops once there are no reactions that can satisfy the above criteria. - At the end of the iterations, we define that the current state have reached NSE + At the end of the iterations, we define that the current state have reached NSE when there is only a single group left, or there are two groups left where 1 of them is the light-isotope-group. + When there is no neutron in the network, it can be difficult for isotopes to form + 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. + + Additional Options ------------------