update on dielectronic work - trying to fix indices for working with …

…a superlevel

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

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;

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);
 }
 
 

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;