Merge remote-tracking branch 'upstream/develop' into mgxs_f_to_c

include/openmc/geometry_aux.h

@@ -7,9 +7,15 @@
 #include <cstdint>
 #include <string>
 #include <vector>
+#include <unordered_map>
 
 namespace openmc {
 
+namespace model {
+  extern std::unordered_map<int32_t, std::unordered_map<int32_t, int32_t>> universe_cell_counts;
+  extern std::unordered_map<int32_t, int32_t> universe_level_counts;
+} // namespace model
+
 void read_geometry_xml();
 
 //==============================================================================

setup.py

@@ -69,7 +69,7 @@
         'pandas', 'lxml', 'uncertainties'
     ],
     'extras_require': {
-        'test': ['pytest', 'pytest-cov'],
+        'test': ['pytest', 'pytest-cov', 'colorama'],
         'vtk': ['vtk'],
     },
 }

src/finalize.cpp

@@ -135,6 +135,10 @@ int openmc_finalize()
 
 int openmc_reset()
 {
+
+  model::universe_cell_counts.clear();
+  model::universe_level_counts.clear();
+
   for (auto& t : model::tallies) {
     t->reset();
   }

src/geometry_aux.cpp

@@ -23,6 +23,21 @@
 
 namespace openmc {
 
+namespace model {
+  std::unordered_map<int32_t, std::unordered_map<int32_t, int32_t>> universe_cell_counts;
+  std::unordered_map<int32_t, int32_t> universe_level_counts;
+} // namespace model
+
+
+// adds the cell counts of universe b to universe a
+void update_universe_cell_count(int32_t a, int32_t b) {
+  auto& universe_a_counts = model::universe_cell_counts[a];
+  const auto& universe_b_counts = model::universe_cell_counts[b];
+  for (const auto& it : universe_b_counts) {
+    universe_a_counts[it.first] += it.second;
+  }
+}
+
 void read_geometry_xml()
 {
 #ifdef DAGMC
@@ -395,19 +410,29 @@ prepare_distribcell()
 void
 count_cell_instances(int32_t univ_indx)
 {
-  for (int32_t cell_indx : model::universes[univ_indx]->cells_) {
-    Cell& c = *model::cells[cell_indx];
-    ++c.n_instances_;
 
-    if (c.type_ == FILL_UNIVERSE) {
-      // This cell contains another universe.  Recurse into that universe.
-      count_cell_instances(c.fill_);
+  const auto univ_counts = model::universe_cell_counts.find(univ_indx);
+  if (univ_counts != model::universe_cell_counts.end()) {
+    for (const auto& it : univ_counts->second) {
+      model::cells[it.first]->n_instances_ += it.second;
+    }
+  } else {
+    for (int32_t cell_indx : model::universes[univ_indx]->cells_) {
+      Cell& c = *model::cells[cell_indx];
+      ++c.n_instances_;
+      model::universe_cell_counts[univ_indx][cell_indx] += 1;
 
-    } else if (c.type_ == FILL_LATTICE) {
-      // This cell contains a lattice.  Recurse into the lattice universes.
-      Lattice& lat = *model::lattices[c.fill_];
-      for (auto it = lat.begin(); it != lat.end(); ++it) {
-        count_cell_instances(*it);
+      if (c.type_ == FILL_UNIVERSE) {
+        // This cell contains another universe.  Recurse into that universe.
+        count_cell_instances(c.fill_);
+        update_universe_cell_count(univ_indx, c.fill_);
+      } else if (c.type_ == FILL_LATTICE) {
+        // This cell contains a lattice.  Recurse into the lattice universes.
+        Lattice& lat = *model::lattices[c.fill_];
+        for (auto it = lat.begin(); it != lat.end(); ++it) {
+          count_cell_instances(*it);
+          update_universe_cell_count(univ_indx, *it);
+        }
       }
     }
   }
@@ -529,6 +554,12 @@ distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset)
 int
 maximum_levels(int32_t univ)
 {
+
+  const auto level_count = model::universe_level_counts.find(univ);
+  if (level_count != model::universe_level_counts.end()) {
+    return level_count->second;
+  }
+
   int levels_below {0};
 
   for (int32_t cell_indx : model::universes[univ]->cells_) {
@@ -546,6 +577,7 @@ maximum_levels(int32_t univ)
   }
 
   ++levels_below;
+  model::universe_level_counts[univ] = levels_below;
   return levels_below;
 }
 

src/mgxs.cpp

@@ -541,8 +541,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
   double nu_fission = xs_t->nu_fission(cache[tid].a, gin);
 
   // Find the probability of having a prompt neutron
-  double prob_prompt =
-       xs_t->prompt_nu_fission(cache[tid].a, gin);
+  double prob_prompt = xs_t->prompt_nu_fission(cache[tid].a, gin);
 
   // sample random numbers
   double xi_pd = prn() * nu_fission;
@@ -557,8 +556,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
 
     // sample the outgoing energy group
     gout = 0;
-    double prob_gout =
-         xs_t->chi_prompt(cache[tid].a, gin, gout);
+    double prob_gout = xs_t->chi_prompt(cache[tid].a, gin, gout);
     while (prob_gout < xi_gout) {
       gout++;
       prob_gout += xs_t->chi_prompt(cache[tid].a, gin, gout);
@@ -568,24 +566,20 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
     // the neutron is delayed
 
     // get the delayed group
-    dg = 0;
-    while (xi_pd >= prob_prompt) {
-      dg++;
-      prob_prompt +=
-           xs_t->delayed_nu_fission(cache[tid].a, dg, gin);
+    for (dg = 0; dg < num_delayed_groups; ++dg) {
+      prob_prompt += xs_t->delayed_nu_fission(cache[tid].a, dg, gin);
+      if (xi_pd < prob_prompt) break;
     }
 
     // adjust dg in case of round-off error
     dg = std::min(dg, num_delayed_groups - 1);
 
     // sample the outgoing energy group
     gout = 0;
-    double prob_gout =
-         xs_t->chi_delayed(cache[tid].a, dg, gin, gout);
+    double prob_gout = xs_t->chi_delayed(cache[tid].a, dg, gin, gout);
     while (prob_gout < xi_gout) {
       gout++;
-      prob_gout +=
-           xs_t->chi_delayed(cache[tid].a, dg, gin, gout);
+      prob_gout += xs_t->chi_delayed(cache[tid].a, dg, gin, gout);
     }
   }
 }

src/particle_restart.cpp

@@ -79,6 +79,9 @@ void run_particle_restart()
   int previous_run_mode;
   read_particle_restart(p, previous_run_mode);
 
+  // write track if that was requested on command line
+  if (settings::write_all_tracks) p.write_track_ = true;
+
   // Set all tallies to 0 for now (just tracking errors)
   model::tallies.clear();
 

src/tallies/tally_scoring.cpp

@@ -1720,7 +1720,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
                 model::tally_filters[i_dg_filt].get())};
               // Tally each delayed group bin individually
               for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
-                auto d = filt.groups()[d_bin];
+                auto d = filt.groups()[d_bin] - 1;
                 score = p->wgt_absorb_ * flux;
                 if (i_nuclide >= 0) {
                   score *=
@@ -1802,7 +1802,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
             model::tally_filters[i_dg_filt].get())};
           // Tally each delayed group bin individually
           for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
-            auto d = filt.groups()[d_bin];
+            auto d = filt.groups()[d_bin] - 1;
             if (i_nuclide >= 0) {
               score = flux * atom_density
                 * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION,
@@ -1842,7 +1842,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
                 model::tally_filters[i_dg_filt].get())};
               // Tally each delayed group bin individually
               for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
-                auto d = filt.groups()[d_bin];
+                auto d = filt.groups()[d_bin] - 1;
                 score = p->wgt_absorb_ * flux;
                 if (i_nuclide >= 0) {
                   score *=
@@ -1902,8 +1902,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           for (auto i = 0; i < p->n_bank_; ++i) {
             auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
             const auto& bank = simulation::fission_bank[i_bank];
-            auto g = bank.delayed_group;
-            if (g != 0) {
+            auto g = bank.delayed_group - 1;
+            if (g != -1) {
               if (i_nuclide >= 0) {
                 score += simulation::keff * atom_density * bank.wgt * flux
                   * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g,
@@ -1923,7 +1923,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
                 // Find the corresponding filter bin and then score
                 for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
                   auto d = filt.groups()[d_bin];
-                  if (d == g)
+                  if (d == g + 1)
                     score_fission_delayed_dg(i_tally, d_bin, score,
                       score_index);
                 }
@@ -1941,7 +1941,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
             model::tally_filters[i_dg_filt].get())};
           // Tally each delayed group bin individually
           for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
-            auto d = filt.groups()[d_bin];
+            auto d = filt.groups()[d_bin] - 1;
             if (i_nuclide >= 0) {
               score += atom_density * flux
                 * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE,