From 0256efe673f02abacac99d297ae4cecc7a01422a Mon Sep 17 00:00:00 2001 From: ssim Date: Tue, 14 Jan 2025 13:26:15 +0100 Subject: [PATCH] update on dielectronic work - trying to fix indices for working with a superlevel --- atomic.h | 11 +++++++++-- globals.h | 1 + input.cc | 16 ++++++++++++++++ nltepop.cc | 46 +++++++++++++++++++++++++++++++++++++--------- 4 files changed, 63 insertions(+), 11 deletions(-) diff --git a/atomic.h b/atomic.h index ccf4604c4..2b00ff688 100644 --- a/atomic.h +++ b/atomic.h @@ -302,9 +302,16 @@ inline auto get_includedlevels() -> int { return includedlevels; } return globals::elements[element].ions[ion].nlevels_nlte; } +//Returns the number of autoionising levels for an ion +[[nodiscard]] inline auto get_nlevels_autoion(const int element, const int ion) -> int { + assert_testmodeonly(element < get_nelements()); + assert_testmodeonly(ion < get_nions(element)); + return globals::elements[element].ions[ion].nlevels_autoion; +} + // ion has NLTE levels, but this one is not NLTE => is in the superlevel [[nodiscard]] inline auto level_isinsuperlevel(const int element, const int ion, const int level) -> bool { - return (!is_nlte(element, ion, level) && level != 0 && (get_nlevels_nlte(element, ion) > 0)); + return (!is_nlte(element, ion, level) && level != 0 && (get_nlevels_nlte(element, ion) > 0) && level < get_nlevels(element, ion) - get_nlevels_autoion(element, ion)); } [[nodiscard]] inline auto get_nlevels_groundterm(const int element, const int ion) -> int { @@ -374,7 +381,7 @@ inline auto get_includedlevels() -> int { return includedlevels; } [[nodiscard]] inline auto ion_has_superlevel(const int element, const int ion) -> bool { assert_testmodeonly(element < get_nelements()); assert_testmodeonly(ion < get_nions(element)); - return (get_nlevels(element, ion) > get_nlevels_nlte(element, ion) + 1); + return (get_nlevels(element, ion) > get_nlevels_nlte(element, ion) + get_nlevels_autoion(element, ion) + 1); } // the number of downward bound-bound transitions from the specified level diff --git a/globals.h b/globals.h index 38d0b0b56..2e5f1642b 100644 --- a/globals.h +++ b/globals.h @@ -114,6 +114,7 @@ struct Ion { int first_nlte; // index into nlte_pops array of a grid cell int ionisinglevels; // Number of levels which have a bf-continuum int maxrecombininglevel; // level index of the highest level with a non-zero recombination rate + int nlevels_autoion; // Number of levels that can autoionise int nlevels_groundterm; int coolingoffset; int ncoolingterms; diff --git a/input.cc b/input.cc index 0a5dcfcb5..9d0baa842 100644 --- a/input.cc +++ b/input.cc @@ -937,6 +937,22 @@ void read_autoion_data() { temp_allautoion.clear(); temp_allautoion.shrink_to_fit(); + // Plan is that autoionizing levels will be explicitly included in the NLTE population solver, but that their level populations do not need to be accurately known - so if the ion has a superlevel already, then we will try to attach the autoionizing level populations to that for all purposes outside the NLTE solber. For this, the ions need to know how many autoionizing levels they have. So count those up now. + + int checkautoall = 0; + for (int element = 0; element < get_nelements(); element++) { + const int nions = get_nions(element); + for (int ion = 0; ion < nions; ion++) { + const int nlevels = get_nlevels(element, ion); + int nauto = 0; + for (int level = 0; level < nlevels; level++) { + nauto += get_nautoiondowntrans(element, ion, level); + } + globals::elements[element].ions[ion].nlevels_autoion = nauto; + checkautoall += nauto; + } + } + assert_always(checkautoall == num_autoion); } diff --git a/nltepop.cc b/nltepop.cc index 37858cf08..fbbd90728 100644 --- a/nltepop.cc +++ b/nltepop.cc @@ -42,8 +42,20 @@ auto get_nlte_vector_index(const int element, const int ion, const int level) -> // have to convert from nlte_pops index to nlte_vector index // the difference is that nlte vectors apply to a single element and include ground states // The (+ ion) term accounts for the ground state population indices that are not counted in the NLTE array + + int offset_autoion = 0; + for (int dion = 0; dion < ion; dion++) { + offset_autoion += get_nlevels_autoion(element,dion); + } + const int gs_index = - globals::elements[element].ions[ion].first_nlte - globals::elements[element].ions[0].first_nlte + ion; + globals::elements[element].ions[ion].first_nlte - globals::elements[element].ions[0].first_nlte + ion + offset_autoion; + + const int first_autoion = get_nlevels(element,ion) - get_nlevels_autoion(element,ion); + + if (ion_has_superlevel(element, ion) && level > (first_autoion - 1)) { + return (gs_index + get_nlevels_nlte(element, ion) + level - first_autoion + 2); + } // add in level or superlevel number const int level_index = gs_index + (is_nlte(element, ion, level) ? level : (get_nlevels_nlte(element, ion) + 1)); @@ -60,6 +72,7 @@ auto get_nlte_vector_index(const int element, const int ion, const int level) -> } } } + printout("index %d element %d nions %d\n", index, element, get_nions(element)); assert_always(false); return {-1, -1}; } @@ -402,7 +415,9 @@ auto get_element_superlevelpartfuncs(const int nonemptymgi, const int element) - if (ion_has_superlevel(element, ion)) { const int nlevels_nlte = get_nlevels_nlte(element, ion); const int nlevels = get_nlevels(element, ion); - for (int level = nlevels_nlte + 1; level < nlevels; level++) { + const int nlevels_autoion = get_nlevels_autoion(element, ion); + //If a superlevel exists, the levels in it will be from nlevels_nlte up to the total-numautoion + for (int level = nlevels_nlte + 1; level < (nlevels - nlevels_autoion); level++) { superlevel_partfuncs[ion] += superlevel_boltzmann(nonemptymgi, element, ion, level); } } @@ -417,11 +432,18 @@ auto get_element_superlevelpartfuncs(const int nonemptymgi, const int element) - for (int ion = 0; ion < nions; ion++) { const int nlevels_nlte = get_nlevels_nlte(element, ion); - nlte_dimension += nlevels_nlte + 1; // ground state is not counted in nlevels_nlte + // nlte_dimension += nlevels_nlte + 1; // ground state is not counted in nlevels_nlte // add super level if it exists if (ion_has_superlevel(element, ion)) { - nlte_dimension++; + nlte_dimension += nlevels_nlte + get_nlevels_autoion(element, ion) + 2; + //printout("Here 1: For element %d ion %d adding %d to nlte_dimension. \n", element, ion, nlevels_nlte + get_nlevels_autoion(element, ion) + 2); + //printout("checks: %d %d\n", nlevels_nlte, get_nlevels_autoion(element, ion)); + //if it has a superlevel then need + 2 for ground state and super and to add autionising levels + } + else { // if it doesn't have a superlevel + nlte_dimension += get_nlevels(element, ion); + //printout("Here 2: For element %d ion %d adding %d to nlte_dimension. \n", element, ion, get_nlevels(element,ion)); } } @@ -677,10 +699,10 @@ void nltepop_matrix_normalise(const int nonemptymgi, const int element, gsl_matr gsl_vector_set(pop_norm_factor_vec, column, calculate_levelpop_lte(nonemptymgi, element, ion, level)); - if ((level != 0) && (!is_nlte(element, ion, level))) { + if (level_isinsuperlevel(element, ion, level)) { // level is a superlevel, so add populations of higher levels to the norm factor for (int dummylevel = level + 1; dummylevel < get_nlevels(element, ion); dummylevel++) { - if (!is_nlte(element, ion, dummylevel)) { + if (level_isinsuperlevel(element, ion, dummylevel)) { *gsl_vector_ptr(pop_norm_factor_vec, column) += calculate_levelpop_lte(nonemptymgi, element, ion, dummylevel); } } @@ -976,7 +998,11 @@ void solve_nlte_pops_element(const int element, const int nonemptymgi, const int std::fill_n(s_renorm.begin(), nlevels_nlte + 1, 1.); for (int level = (nlevels_nlte + 1); level < nlevels; level++) { - s_renorm[level] = superlevel_boltzmann(nonemptymgi, element, ion, level) / superlevel_partfunc[ion]; + if (level_isinsuperlevel(element, ion, level)) { + s_renorm[level] = superlevel_boltzmann(nonemptymgi, element, ion, level) / superlevel_partfunc[ion]; + } else { //next clause is for any autoionising levels above the superlevel + s_renorm[level] = 1.; + } } nltepop_matrix_add_boundbound(nonemptymgi, element, ion, t_mid, s_renorm, &rate_matrix_rad_bb, &rate_matrix_coll_bb, @@ -1246,8 +1272,10 @@ void nltepop_write_to_file(const int nonemptymgi, const int timestep) { double superlevel_partfunc = 0; fprintf(nlte_file, "%d ", -1); for (int level_sl = nlevels_nlte + 1; level_sl < get_nlevels(element, ion); level_sl++) { - nnlevellte += calculate_levelpop_lte(nonemptymgi, element, ion, level_sl); - superlevel_partfunc += superlevel_boltzmann(nonemptymgi, element, ion, level_sl); + if (level_isinsuperlevel(element, ion, level_sl)) { + nnlevellte += calculate_levelpop_lte(nonemptymgi, element, ion, level_sl); + superlevel_partfunc += superlevel_boltzmann(nonemptymgi, element, ion, level_sl); + } } nnlevelnlte = slpopfactor * superlevel_partfunc;