From 60015e5c5c34e6a606a5ddaea94b6163da4ad44f Mon Sep 17 00:00:00 2001
From: John Tramm <john.tramm@gmail.com>
Date: Tue, 26 Nov 2019 17:05:04 +0000
Subject: [PATCH] Revert "Merge branch 'event-based-device' into
 particle_owned_seeds"

This reverts commit 12463ac6b0f82045eebc1734254b734f45f7ee45, reversing
changes made to 42616905379cd2a350e85c14944b80d99a697070.
---
 examples/xml/pincell/settings.xml |   4 +-
 include/openmc/geometry.h         |   2 -
 include/openmc/particle.h         |  24 --
 include/openmc/simulation.h       |  11 +-
 src/distribution_angle.cpp        |   1 -
 src/material.cpp                  |   6 -
 src/output.cpp                    |  10 +-
 src/particle.cpp                  |  14 +-
 src/physics.cpp                   |  31 +-
 src/simulation.cpp                | 694 +-----------------------------
 src/tallies/tally.cpp             |   1 -
 11 files changed, 12 insertions(+), 786 deletions(-)

diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml
index 7b0a985f10..5845898c8d 100644
--- a/examples/xml/pincell/settings.xml
+++ b/examples/xml/pincell/settings.xml
@@ -3,8 +3,8 @@
 
   <!-- Define how many particles to run and for how many batches -->
   <run_mode>eigenvalue</run_mode>
-  <batches>10</batches>
-  <inactive>5</inactive>
+  <batches>100</batches>
+  <inactive>10</inactive>
   <particles>1000</particles>
 
   <!-- The starting source is a uniform distribution over the entire pin
diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h
index a6e1e40631..e386348c44 100644
--- a/include/openmc/geometry.h
+++ b/include/openmc/geometry.h
@@ -28,14 +28,12 @@ extern std::vector<int64_t> overlap_check_count;
 // Information about nearest boundary crossing
 //==============================================================================
 
-/*
 struct BoundaryInfo {
   double distance {INFINITY};   //!< distance to nearest boundary
   int surface_index {0}; //!< if boundary is surface, index in surfaces vector
   int coord_level;   //!< coordinate level after crossing boundary
   std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
 };
-*/
 
 //==============================================================================
 //! Check two distances by coincidence tolerance
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 760e2ad8af..6ac8aabaea 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -129,21 +129,8 @@ struct MacroXS {
   double incoherent;      //!< macroscopic incoherent xs
   double photoelectric;   //!< macroscopic photoelectric xs
   double pair_production; //!< macroscopic pair production xs
-  int i_grid;
 };
 
-//==============================================================================
-// Information about nearest boundary crossing
-//==============================================================================
-
-struct BoundaryInfo {
-  double distance {INFINITY};   //!< distance to nearest boundary
-  int surface_index {0}; //!< if boundary is surface, index in surfaces vector
-  int coord_level;   //!< coordinate level after crossing boundary
-  std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
-};
-
-
 //============================================================================
 //! State of a particle being transported through geometry
 //============================================================================
@@ -286,7 +273,6 @@ public:
 
   // Indices for various arrays
   int surface_ {0};             //!< index for surface particle is on
-  BoundaryInfo boundary_;
   int cell_born_ {-1};      //!< index for cell particle was born in
   int material_ {-1};       //!< index for current material
   int material_last_ {-1};  //!< index for last material
@@ -308,14 +294,4 @@ public:
 
 } // namespace openmc
 
-extern template class std::vector<openmc::Particle::Bank>;
-
-namespace openmc{
-
-extern std::vector<Particle> particle_bank;
-#pragma omp threadprivate(particle_bank)
-
-}
-
-
 #endif // OPENMC_PARTICLE_H
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index 7981286be2..8efac4c7f8 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -16,10 +16,6 @@ constexpr int STATUS_EXIT_NORMAL {0};
 constexpr int STATUS_EXIT_MAX_BATCH {1};
 constexpr int STATUS_EXIT_ON_TRIGGER {2};
 
-extern Particle::Bank * shared_fission_bank;
-extern int shared_fission_bank_length;
-extern int shared_fission_bank_max;
-
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -42,9 +38,7 @@ extern "C" int restart_batch;   //!< batch at which a restart job resumed
 extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
 extern "C" int total_gen;        //!< total number of generations simulated
 extern double total_weight;  //!< Total source weight in a batch
-extern int64_t thread_work_index;
-extern int64_t work_per_rank;       //!< number of particles per MPI rank
-extern int64_t work_per_thread;    //!< number of particles on a given tread
+extern int64_t work_per_rank;         //!< number of particles per MPI rank
 
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;
@@ -55,8 +49,7 @@ extern std::vector<int64_t> work_index;
 // Threadprivate variables
 extern "C" bool trace;     //!< flag to show debug information
 
-
-#pragma omp threadprivate(current_work, work_per_thread, trace)
+#pragma omp threadprivate(current_work, trace)
 
 } // namespace simulation
 
diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index f946d0602e..59263afc89 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -85,7 +85,6 @@ double AngleDistribution::sample(double E, uint64_t * prn_seeds, int stream) con
   // Sample between the ith and (i+1)th bin
   if (r > prn(prn_seeds, stream)) ++i;
 
-  //assert(i >= 0 );
   // Sample i-th distribution
   double mu = distribution_[i]->sample(prn_seeds, stream);
 
diff --git a/src/material.cpp b/src/material.cpp
index 16b9275a05..a5bc3dc138 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -751,7 +751,6 @@ void Material::calculate_xs(Particle& p) const
 
 void Material::calculate_neutron_xs(Particle& p) const
 {
-	/*
   // Find energy index on energy grid
   int neutron = static_cast<int>(Particle::Type::neutron);
   int i_grid = std::log(p.E_/data::energy_min[neutron])/simulation::log_spacing;
@@ -761,11 +760,9 @@ void Material::calculate_neutron_xs(Particle& p) const
 
   // Initialize position in i_sab_nuclides
   int j = 0;
-  */
 
   // Add contribution from each nuclide in material
   for (int i = 0; i < nuclide_.size(); ++i) {
-	  /*
     // ======================================================================
     // CHECK FOR S(A,B) TABLE
 
@@ -795,21 +792,18 @@ void Material::calculate_neutron_xs(Particle& p) const
 
     // ======================================================================
     // CALCULATE MICROSCOPIC CROSS SECTION
-	*/
 
     // Determine microscopic cross sections for this nuclide
     int i_nuclide = nuclide_[i];
 
     // Calculate microscopic cross section for this nuclide
     const auto& micro {p.neutron_xs_[i_nuclide]};
-	/*
     if (p.E_ != micro.last_E
         || p.sqrtkT_ != micro.last_sqrtkT
         || i_sab != micro.index_sab
         || sab_frac != micro.sab_frac) {
       data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p);
     }
-	*/
 
     // ======================================================================
     // ADD TO MACROSCOPIC CROSS SECTION
diff --git a/src/output.cpp b/src/output.cpp
index 29d60830eb..899bce86e2 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -67,15 +67,7 @@ void title()
     "                     ###############     %%%%%%%%%%%%%%%%\n" <<
     "                       ############     %%%%%%%%%%%%%%%\n" <<
     "                          ########     %%%%%%%%%%%%%%\n" <<
-    "                                      %%%%%%%%%%%\n" <<
-    "                              )                             (        (\n" << 
-    "                           ( /(   *   )        (     (      )\\ )     )\\ )\n" <<  
-    "         (    (   (   (    )\\())` )  /(      ( )\\    )\\    (()/( (  (()/(\n" << 
-    "         )\\   )\\  )\\  )\\  ((_)\\  ( )(_))___  )((_)((((_)(   /(_)))\\  /(_))\n" << 
-    "        ((_) ((_)((_)((_)  _((_)(_(_())|___|((_)_  )\\ _ )\\ (_)) ((_)(_))_\n" <<  
-    "        | __|\\ \\ / / | __|| \\| ||_   _|      | _ ) (_)_\\(_)/ __|| __||   \\\n" << 
-    "        | _|  \\ V /  | _| | .` |  | |        | _ \\  / _ \\  \\__ \\| _| | |) |\n" << 
-    "        |___|  \\_/   |___||_|\\_|  |_|        |___/ /_/ \\_\\ |___/|___||___/\n\n";
+    "                                      %%%%%%%%%%%\n\n";
 
   // Write version information
   std::cout <<
diff --git a/src/particle.cpp b/src/particle.cpp
index 12334f149e..98ee8920bb 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -27,14 +27,8 @@
 #include "openmc/tallies/tally_scoring.h"
 #include "openmc/track_output.h"
 
-// Explicit template instantiation definition
-template class std::vector<openmc::Particle>;
-
-
 namespace openmc {
 
-std::vector<Particle> particle_bank;
-
 //==============================================================================
 // LocalCoord implementation
 //==============================================================================
@@ -138,15 +132,14 @@ void
 Particle::transport()
 {
   // Display message if high verbosity or trace is on
-  //if (settings::verbosity >= 9 || simulation::trace) {
-    // write_message("Simulating Particle " + std::to_string(id_));
-  //}
+  if (settings::verbosity >= 9 || simulation::trace) {
+     write_message("Simulating Particle " + std::to_string(id_));
+  }
 
   // Initialize number of events to zero
   int n_event = 0;
 
   // Add paricle's starting weight to count for normalizing tallies later
-  /*
   #pragma omp atomic
   simulation::total_weight += wgt_;
 
@@ -160,7 +153,6 @@ Particle::transport()
 
   // Every particle starts with no accumulated flux derivative.
   if (!model::active_tallies.empty()) zero_flux_derivs();
-  */
 
   while (true) {
     // Set the random number stream
diff --git a/src/physics.cpp b/src/physics.cpp
index 84927c27e1..28fc4341e0 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -79,7 +79,6 @@ void collision(Particle* p)
 
 void sample_neutron_reaction(Particle* p)
 {
-	//std::cout << "particle energy e = " << p->E_ << std::endl;
   // Sample a nuclide within the material
   int i_nuclide = sample_nuclide(p);
 
@@ -128,12 +127,9 @@ void sample_neutron_reaction(Particle* p)
   }
   if (!p->alive_) return;
 
-	//std::cout << "before scatter particle energy e = " << p->E_ << " with nuclide = " << i_nuclide <<  std::endl;
   // Sample a scattering reaction and determine the secondary energy of the
   // exiting neutron
   scatter(p, i_nuclide);
-	
-  //std::cout << "after scatter particle energy e = " << p->E_ << std::endl;
 
   // Advance URR seed stream 'N' times after energy changes
   if (p->E_ != p->E_last_) {
@@ -175,7 +171,6 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 
   p->fission_ = true;
   for (int i = 0; i < nu; ++i) {
-	  /*
     // Create new bank site and get reference to last element
     bank.emplace_back();
     auto& site {bank.back()};
@@ -184,22 +179,12 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
     site.r = p->r();
     site.particle = Particle::Type::neutron;
     site.wgt = 1. / weight;
-	*/
-	  int idx;
-	  #pragma omp atomic capture
-	  idx = shared_fission_bank_length++;
-	  Particle::Bank * site = shared_fission_bank + idx;
-	  site->r = p->r();
-	  site->particle = Particle::Type::neutron;
-	  site->wgt = 1. / weight;
 
     // Sample delayed group and angle/energy for fission reaction
-
-    sample_fission_neutron(i_nuclide, rx, p->E_, site, p->prn_seeds_, p->stream_);
+    sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->prn_seeds_, p->stream_);
 
     // Set the delayed group on the particle as well
-    //p->delayed_group_ = site.delayed_group;
-    p->delayed_group_ = site->delayed_group;
+    p->delayed_group_ = site.delayed_group;
 
     // Increment the number of neutrons born delayed
     if (p->delayed_group_ > 0) {
@@ -613,7 +598,6 @@ void scatter(Particle* p, int i_nuclide)
     // Determine temperature
     double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp];
 
-	//std::cout << "we are doing an elastic scatter" << std::endl;
     // Perform collision physics for elastic scattering
     elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p);
 
@@ -626,7 +610,6 @@ void scatter(Particle* p, int i_nuclide)
     // =======================================================================
     // S(A,B) SCATTERING
 
-	//std::cout << "we are doing an SAB scatter" << std::endl;
     sab_scatter(i_nuclide, micro.index_sab, p);
 
     p->event_mt_ = ELASTIC;
@@ -634,7 +617,6 @@ void scatter(Particle* p, int i_nuclide)
   }
 
   if (!sampled) {
-	//std::cout << "we are doing an Inelastic scatter" << std::endl;
     // =======================================================================
     // INELASTIC SCATTERING
 
@@ -661,9 +643,7 @@ void scatter(Particle* p, int i_nuclide)
 
     // Perform collision physics for inelastic scattering
     const auto& rx {nuc->reactions_[i]};
-	//std::cout << "Energy before Inelastic scatter = " << p->E_ << std::endl;
     inelastic_scatter(nuc.get(), rx.get(), p);
-	//std::cout << "Energy after Inelastic scatter = " << p->E_ << std::endl;
     p->event_mt_ = rx->mt_;
   }
 
@@ -1076,29 +1056,22 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
   // angle and outgoing energy from CM to LAB
   if (rx->scatter_in_cm_) {
     double E_cm = E;
-	//std::cout << "E_cm = " << E_cm << std::endl;
 
     // determine outgoing energy in lab
     double A = nuc->awr_;
-//	std::cout << "A = " << nuc->awr_ << std::endl;
     E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm))
           / ((A + 1.0)*(A + 1.0));
-	// here's where it goes wrong
-	//std::cout << "E = " << E << std::endl;
 
     // determine outgoing angle in lab
     mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E);
-//	std::cout << "mu = " << mu << std::endl;
   }
 
   // Because of floating-point roundoff, it may be possible for mu to be
   // outside of the range [-1,1). In these cases, we just set mu to exactly -1
   // or 1
   if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
-//	std::cout << "mu = " << mu << std::endl;
 
   // Set outgoing energy and scattering angle
-//	std::cout << "FINAL E = " << E << " mu = " << mu << std::endl;
   p->E_ = E;
   p->mu_ = mu;
 
diff --git a/src/simulation.cpp b/src/simulation.cpp
index b2af14a3f7..86f20a93d3 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -2,32 +2,23 @@
 
 #include "openmc/bank.h"
 #include "openmc/capi.h"
-#include "openmc/cell.h"
 #include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
-#include "openmc/geometry.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
-#include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
 #include "openmc/particle.h"
 #include "openmc/photon.h"
-#include "openmc/physics.h"
- #include "openmc/physics_mg.h"
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
 #include "openmc/source.h"
 #include "openmc/state_point.h"
-#include "openmc/thermal.h"
 #include "openmc/timer.h"
-#include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
-#include "openmc/tallies/tally_scoring.h"
 #include "openmc/tallies/trigger.h"
-#include "openmc/track_output.h"
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -37,604 +28,6 @@
 #include <algorithm>
 #include <string>
 
-namespace openmc {
-
-	/*
-extern std::vector<Particle*> calculate_fuel_xs_queue;
-extern std::vector<Particle*> calculate_nonfuel_xs_queue;
-extern std::vector<Particle*> advance_particle_queue;
-extern std::vector<Particle*> surface_crossing_queue;
-extern std::vector<Particle*> collision_queue;
-#pragma omp threadprivate(calculate_fuel_xs_queue, calculate_nonfuel_xs_queue, advance_particle_queue, surface_crossing_queue, collision_queue)
-
-std::vector<Particle*> calculate_fuel_xs_queue;
-std::vector<Particle*> calculate_nonfuel_xs_queue;
-std::vector<Particle*> advance_particle_queue;
-std::vector<Particle*> surface_crossing_queue;
-std::vector<Particle*> collision_queue;
-*/
-
-int * calculate_fuel_xs_queue;
-int * calculate_nonfuel_xs_queue;
-int * advance_particle_queue;
-int * surface_crossing_queue;
-int * collision_queue;
-Particle * particles;
-
-int calculate_fuel_xs_queue_length    = 0;
-int calculate_nonfuel_xs_queue_length = 0;
-int advance_particle_queue_length     = 0;
-int surface_crossing_queue_length     = 0;
-int collision_queue_length            = 0;
-
-const int MAX_PARTICLES_IN_FLIGHT = 100;
-
-void init_event_queues(void)
-{
-	int n_particles = MAX_PARTICLES_IN_FLIGHT;
-	calculate_fuel_xs_queue =    new int[n_particles];
-	calculate_nonfuel_xs_queue = new int[n_particles];
-	advance_particle_queue =     new int[n_particles];
-	surface_crossing_queue =     new int[n_particles];
-	collision_queue =            new int[n_particles];
-	particles =                  new Particle[n_particles];
-}
-
-void free_event_queues(void)
-{
-	delete calculate_fuel_xs_queue;
-	delete calculate_nonfuel_xs_queue;
-	delete advance_particle_queue;
-	delete surface_crossing_queue;
-	delete collision_queue;
-	delete[] particles;
-}
-
-constexpr size_t MAX_PARTICLES_PER_THREAD {100};
-
-Particle::Bank * shared_fission_bank;
-int shared_fission_bank_length = 0;
-int shared_fission_bank_max;
-
-void init_shared_fission_bank(int max)
-{
-	shared_fission_bank_max = max;
-	shared_fission_bank = new Particle::Bank[max];
-}
-
-void free_shared_fission_bank(void)
-{
-	delete[] shared_fission_bank;
-	shared_fission_bank_length = 0;
-}
-
-// TODO: What is going on here?
-void revive_particle_from_secondary(Particle* p)
-{
-  p->from_source(&simulation::secondary_bank.back());
-  simulation::secondary_bank.pop_back();
-  // n_event = 0;
-
-  // Enter new particle in particle track file
-  if (p->write_track_) add_particle_track();
-}
-
-void dispatch_xs_event(int i)
-{
-	Particle * p = particles + i;
-	  int idx;
-  if (p->material_ == MATERIAL_VOID) {
-
-	#pragma omp atomic capture
-	idx = calculate_nonfuel_xs_queue_length++;
-	//std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl;
-    calculate_nonfuel_xs_queue[idx] = i;
-  } else {
-    if (model::materials[p->material_]->fissionable_) {
-	#pragma omp atomic capture
-	idx = calculate_fuel_xs_queue_length++;
-	//std::cout << "Dispatching particle to Fuel XS queue idx = " << idx << std::endl;
-      calculate_fuel_xs_queue[idx] = i;
-    } else {
-	#pragma omp atomic capture
-	idx = calculate_nonfuel_xs_queue_length++;
-	//std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl;
-      calculate_nonfuel_xs_queue[idx] = i;
-    }
-  }
-}
-
-void process_calculate_xs_events(int * queue, int n)
-{
-  // Save last_ members, find grid index
-  for (int i = 0; i < n; i++) {
-	  Particle *p = particles + queue[i]; 
-	  //std::cout << "particle offset = " << queue[i] << std::endl;
-    // Set the random number stream
-	// TODO: Move RNG seeds to particle storage
-    if (p->type_ == Particle::Type::neutron) {
-      prn_set_stream(STREAM_TRACKING);
-    } else {
-      prn_set_stream(STREAM_PHOTON);
-    }
-
-    // Store pre-collision particle properties
-    p->wgt_last_ = p->wgt_;
-    p->E_last_ = p->E_;
-    p->u_last_ = p->u();
-    p->r_last_ = p->r();
-
-    // If the cell hasn't been determined based on the particle's location,
-    // initiate a search for the current cell. This generally happens at the
-    // beginning of the history and again for any secondary particles
-    if (p->coord_[p->n_coord_ - 1].cell == C_NONE) {
-      if (!find_cell(p, false)) {
-        p->mark_as_lost("Could not find the cell containing particle "
-          + std::to_string(p->id_));
-        return;
-      }
-
-      // set birth cell attribute
-      if (p->cell_born_ == C_NONE) p->cell_born_ = p->coord_[p->n_coord_ - 1].cell;
-    }
-
-    // Write particle track.
-    if (p->write_track_) write_particle_track(*p);
-
-    if (settings::check_overlaps) check_cell_overlap(p);
-
-    if (settings::run_CE) {
-      if (p->material_ == p->material_last_ && p->sqrtkT_ != p->sqrtkT_last_) {
-        // Remove particle from queue
-      }
-    }
-
-    // Find energy index on energy grid
-    // TODO: Calculate this separately?
-    int neutron = static_cast<int>(Particle::Type::neutron);
-    p->macro_xs_.i_grid = std::log(p->E_/data::energy_min[neutron]) / simulation::log_spacing;
-  }
-
-  // Calculate nuclide micros
-  for (int i = 0; i < data::nuclides.size(); ++i) {
-	  // loop over particles
-    for (int j = 0; j < n; j++) {
-		//Particle * p = particles + queue[i];
-		Particle * p = particles + queue[j];
-      if (p->material_ == MATERIAL_VOID) continue;
-
-      // If material doesn't have this nuclide, skip it
-      const auto& mat {model::materials[p->material_]};
-      if (mat->mat_nuclide_index_[i] == -1) continue;
-
-      // ======================================================================
-      // CHECK FOR S(A,B) TABLE
-
-      // Check if this nuclide matches one of the S(a,b) tables specified.
-      // This relies on thermal_tables_ being sorted by .index_nuclide
-      int i_sab = C_NONE;
-      double sab_frac = 0.0;
-      for (const auto& sab : mat->thermal_tables_) {
-        if (i == sab.index_nuclide) {
-          // Get index in sab_tables
-          i_sab = sab.index_table;
-          sab_frac = sab.fraction;
-
-          // If particle energy is greater than the highest energy for the
-          // S(a,b) table, then don't use the S(a,b) table
-          //if (p->E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE;
-          if (p->E_ > data::thermal_scatt[i_sab]->energy_max_) i_sab = C_NONE;
-        }
-      }
-
-      // ======================================================================
-      // CALCULATE MICROSCOPIC CROSS SECTION
-
-      // Calculate microscopic cross section for this nuclide
-      const auto& micro {p->neutron_xs_[i]};
-      if (p->E_ != micro.last_E
-          || p->sqrtkT_ != micro.last_sqrtkT
-          || i_sab != micro.index_sab
-          || sab_frac != micro.sab_frac) {
-        data::nuclides[i]->calculate_xs(i_sab, p->macro_xs_.i_grid, sab_frac, *p);
-      }
-    }
-  }
-
-  for (int i = 0; i < n; i++) {
-	  Particle * p = particles + queue[i];
-    // Calculate microscopic and macroscopic cross sections
-    if (p->material_ != MATERIAL_VOID) {
-	  // Only works for CE, no MG support
-      //if (settings::run_CE) {
-        // If the material is the same as the last material and the
-        // temperature hasn't changed, we don't need to lookup cross
-        // sections again.
-        model::materials[p->material_]->calculate_xs(*p);
-      }
-		//else {
-      //}
-     else {
-      p->macro_xs_.total      = 0.0;
-      p->macro_xs_.absorption = 0.0;
-      p->macro_xs_.fission    = 0.0;
-      p->macro_xs_.nu_fission = 0.0;
-    }
-
-	 int idx;
-	 #pragma omp atomic capture
-	 idx = advance_particle_queue_length++;
-    advance_particle_queue[idx] = queue[i];
-  }
-}
-
-void process_advance_particle_events()
-{
-  //for (auto& p : advance_particle_queue) {
-  for (int i = 0; i < advance_particle_queue_length; i++) {
-	  Particle * p = particles + advance_particle_queue[i];
-    simulation::trace == (p->id_ == 0);
-
-    // Sample a distance to collision
-    double d_collision;
-    if (p->type_ == Particle::Type::electron ||
-        p->type_ == Particle::Type::positron) {
-      d_collision = 0.0;
-    } else if (p->macro_xs_.total == 0.0) {
-      d_collision = INFINITY;
-    } else {
-      d_collision = -std::log(prn()) / p->macro_xs_.total;
-    }
-
-    // -------------- break here? -------------------
-
-    // Find the distance to the nearest boundary
-    p->boundary_ = distance_to_boundary(p);
-
-    // Select smaller of the two distances
-    double distance;
-	int idx;
-    if (p->boundary_.distance < d_collision) {
-		#pragma omp atomic capture
-		idx = surface_crossing_queue_length++;
-      surface_crossing_queue[idx] = advance_particle_queue[i];
-      distance = p->boundary_.distance;
-    } else {
-		#pragma omp atomic capture
-		idx = collision_queue_length++;
-		collision_queue[idx] = advance_particle_queue[i];
-      distance = d_collision;
-    }
-
-    // -------------- break here? -------------------
-
-    // Advance particle
-    for (int j = 0; j < p->n_coord_; ++j) {
-      p->coord_[j].r += distance * p->coord_[j].u;
-    }
-
-    // -------------- break here? -------------------
-
-    // Score track-length tallies
-    if (!model::active_tracklength_tallies.empty()) {
-      score_tracklength_tally(p, distance);
-    }
-
-    // Score track-length estimate of k-eff
-    if (settings::run_mode == RUN_MODE_EIGENVALUE &&
-        p->type_ == Particle::Type::neutron) {
-      global_tally_tracklength += p->wgt_ * distance * p->macro_xs_.nu_fission;
-    }
-
-    // Score flux derivative accumulators for differential tallies.
-    if (!model::active_tallies.empty()) {
-      score_track_derivative(p, distance);
-    }
-  }
-
-  advance_particle_queue_length = 0;
-}
-
-void process_surface_crossing_events()
-{
-  //for (auto& p : surface_crossing_queue) {
-  for (int i = 0; i < surface_crossing_queue_length; i++) {
-	  Particle * p = particles + surface_crossing_queue[i];
-    // Set surface that particle is on and adjust coordinate levels
-    p->surface_ = p->boundary_.surface_index;
-    p->n_coord_ = p->boundary_.coord_level;
-
-    // Saving previous cell data
-    for (int j = 0; j < p->n_coord_; ++j) {
-      p->cell_last_[j] = p->coord_[j].cell;
-    }
-    p->n_coord_last_ = p->n_coord_;
-
-    if (p->boundary_.lattice_translation[0] != 0 ||
-        p->boundary_.lattice_translation[1] != 0 ||
-        p->boundary_.lattice_translation[2] != 0) {
-      // Particle crosses lattice boundary
-      cross_lattice(p, p->boundary_);
-      p->event_ = EVENT_LATTICE;
-    } else {
-      // Particle crosses surface
-      p->cross_surface();
-      p->event_ = EVENT_SURFACE;
-    }
-    // Score cell to cell partial currents
-    if (!model::active_surface_tallies.empty()) {
-      score_surface_tally(p, model::active_surface_tallies);
-    }
-
-	// TODO: Add back in support for secondary particles
-	/*
-    if (!p->alive_ && !simulation::secondary_bank.empty()) {
-      revive_particle_from_secondary(p);
-    }
-	*/
-
-    if (p->alive_)
-	{
-		dispatch_xs_event(surface_crossing_queue[i]);
-	}
-  }
-
-  surface_crossing_queue_length = 0;
-}
-
-void process_collision_events()
-{
-  //for (auto& p : collision_queue) {
-  for (int i = 0; i < collision_queue_length; i++) {
-	  Particle * p = particles + collision_queue[i];
-	  //std::cout << "Beginning collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-    // Score collision estimate of keff
-    if (settings::run_mode == RUN_MODE_EIGENVALUE &&
-        p->type_ == Particle::Type::neutron) {
-      global_tally_collision += p->wgt_ * p->macro_xs_.nu_fission
-        / p->macro_xs_.total;
-    }
-
-    // Score surface current tallies -- this has to be done before the collision
-    // since the direction of the particle will change and we need to use the
-    // pre-collision direction to figure out what mesh surfaces were crossed
-
-    if (!model::active_meshsurf_tallies.empty())
-      score_surface_tally(p, model::active_meshsurf_tallies);
-	    
-	//std::cout << "After surface tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-
-    // Clear surface component
-    p->surface_ = 0;
-
-    if (settings::run_CE) {
-      collision(p);
-    } else {
-      collision_mg(p);
-    }
-
-	//std::cout << "After collision() of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-
-    // Score collision estimator tallies -- this is done after a collision
-    // has occurred rather than before because we need information on the
-    // outgoing energy for any tallies with an outgoing energy filter
-    if (!model::active_collision_tallies.empty()) score_collision_tally(p);
-    if (!model::active_analog_tallies.empty()) {
-      if (settings::run_CE) {
-        score_analog_tally_ce(p);
-      } else {
-        score_analog_tally_mg(p);
-      }
-    }
-	//std::cout << "After analog tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-
-    // Reset banked weight during collision
-    p->n_bank_ = 0;
-    p->wgt_bank_ = 0.0;
-    for (int& v : p->n_delayed_bank_) v = 0;
-
-    // Reset fission logical
-    p->fission_ = false;
-
-    // Save coordinates for tallying purposes
-    p->r_last_current_ = p->r();
-
-    // Set last material to none since cross sections will need to be
-    // re-evaluated
-    p->material_last_ = C_NONE;
-
-    // Set all directions to base level -- right now, after a collision, only
-    // the base level directions are changed
-    for (int j = 0; j < p->n_coord_ - 1; ++j) {
-      if (p->coord_[j + 1].rotated) {
-        // If next level is rotated, apply rotation matrix
-        const auto& m {model::cells[p->coord_[j].cell]->rotation_};
-        const auto& u {p->coord_[j].u};
-        p->coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z;
-        p->coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z;
-        p->coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z;
-      } else {
-        // Otherwise, copy this level's direction
-        p->coord_[j+1].u = p->coord_[j].u;
-      }
-    }
-
-    // Score flux derivative accumulators for differential tallies.
-    if (!model::active_tallies.empty()) score_collision_derivative(p);
-
-	// TODO: Add back in secondary particle support
-	/*
-    if (!p->alive_ && !simulation::secondary_bank.empty()) {
-      revive_particle_from_secondary(p);
-    }
-	*/
-
-    if (p->alive_)
-	{
-		dispatch_xs_event(collision_queue[i]);
-	    //std::cout << "Ended collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
-		assert(std::isfinite(p->E_) );
-	}
-  }
-
-  collision_queue_length = 0;
-}
-
-void check_energies(void)
-{
-	int * Q;
-	int n;
-
-
-	Q = calculate_fuel_xs_queue;
-	n = calculate_fuel_xs_queue_length;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index xs FUEL " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = calculate_nonfuel_xs_queue;
-	n = calculate_nonfuel_xs_queue_length ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index xs Non fuel " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = advance_particle_queue;
-	n = advance_particle_queue_length     ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index advance particle " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = surface_crossing_queue;
-	n = surface_crossing_queue_length     ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index surface crossing " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-	Q = collision_queue;
-	n = collision_queue_length            ;
-	for( int i = 0; i < n; i++ )
-	{
-		if( !std::isfinite(particles[Q[i]].E_ ) )
-		{
-			std::cout << "NAN energy particle found at index collision " << Q[i] << std::endl;
-			assert(0);
-		}
-	}
-}
-
-void transport()
-{
-	int remaining_work = simulation::work_per_rank;
-	int source_offset = 0;
-
-
-
-	// Subiterations to complete sets of particles
-	while (remaining_work > 0) {
-
-		// TODO: Do this only once per PI
-		init_event_queues();
-
-		// Figure out work for this subiteration
-		int n_particles = MAX_PARTICLES_IN_FLIGHT;
-		if( n_particles > remaining_work)
-			n_particles = remaining_work;
-
-		//std::cout << "Initializing particle histories..." << std::endl;
-		// Initialize all histories
-		// TODO: Parallelize
-		for (int i = 0; i < n_particles; i++) {
-			initialize_history(particles + i, source_offset + i);
-		}
-
-		//std::cout << "Enqueing particles for XS Lookups..." << std::endl;
-		// Add all particles to advance particle queue
-		// TODO: Parallelize
-		for (int i = 0; i < n_particles; i++) {
-			dispatch_xs_event(i);
-		}
-
-		int event_kernel_executions = 0;
-		while (true) {
-			event_kernel_executions++;
-			/*
-			std::cout << "Fuel XS Lookups = " << calculate_fuel_xs_queue_length << std::endl;
-			std::cout << "Non Fuel XS Lookups = " << calculate_nonfuel_xs_queue_length << std::endl;
-			std::cout << "Advance Particles = " << advance_particle_queue_length << std::endl;
-			std::cout << "Surface Crossings = " << surface_crossing_queue_length << std::endl;
-			std::cout << "Collisions = " << collision_queue_length << std::endl;
-			*/
-			int max = std::max({calculate_fuel_xs_queue_length, calculate_nonfuel_xs_queue_length, advance_particle_queue_length, surface_crossing_queue_length, collision_queue_length});
-			check_energies();
-			if (max == 0) {
-				break;
-			} else if (max == calculate_fuel_xs_queue_length) {
-				//std::cout << "pre fuel XS check..." << std::endl;
-				//check_energies(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
-				//std::cout << "Performing Fuel XS Lookups..." << std::endl;
-				process_calculate_xs_events(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
-				calculate_fuel_xs_queue_length = 0;
-			} else if (max == calculate_nonfuel_xs_queue_length) {
-				//std::cout << "pre non fuel XS check..." << std::endl;
-				//check_energies(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
-			//	std::cout << "Performing Non Fuel XS Lookups..." << std::endl;
-				process_calculate_xs_events(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
-				calculate_nonfuel_xs_queue_length = 0;
-			} else if (max == advance_particle_queue_length) {
-				//std::cout << "pre advancing check..." << std::endl;
-				//check_energies(advance_particle_queue, advance_particle_queue_length);
-				//std::cout << "Advancing Particles..." << std::endl;
-				process_advance_particle_events();
-			} else if (max == surface_crossing_queue_length) {
-				//std::cout << "pre surface crossing check..." << std::endl;
-				//check_energies(surface_crossing_queue, surface_crossing_queue_length);
-				//std::cout << "Surface Crossings..." << std::endl;
-				process_surface_crossing_events();
-			} else if (max == collision_queue_length) {
-				//std::cout << "pre Colliding check..." << std::endl;
-				//check_energies(collision_queue, collision_queue_length);
-				//std::cout << "Colliding..." << std::endl;
-				process_collision_events();
-			}
-		}
-
-		remaining_work -= n_particles;
-		source_offset += n_particles;
-
-		// Should all be zero
-		/*
-		calculate_fuel_xs_queue_length    = 0;
-		calculate_nonfuel_xs_queue_length = 0;
-		advance_particle_queue_length     = 0;
-		surface_crossing_queue_length     = 0;
-		collision_queue_length            = 0;
-		*/
-		
-		// TODO: Do this only once per PI
-		free_event_queues();
-		std::cout << "Event kernels retired: " << event_kernel_executions << std::endl;
-	}
-	//shared_fission_bank_length = 0;
-}
-
-} // namespace openmc
-
 //==============================================================================
 // C API functions
 //==============================================================================
@@ -649,12 +42,10 @@ int openmc_run()
 
   int err = 0;
   int status = 0;
-
   while (status == 0 && err == 0) {
     err = openmc_next_batch(&status);
   }
 
-
   openmc_simulation_finalize();
   return err;
 }
@@ -665,16 +56,9 @@ int openmc_simulation_init()
 
   // Skip if simulation has already been initialized
   if (simulation::initialized) return 0;
-	
 
   // Determine how much work each process should do
   calculate_work();
-  
-  // We assume that the number of fission sites will be at most
-	// 10x the number of neutrons run on this rank
-	int shared_fission_bank_size = 10 * simulation::work_per_rank;
-	// Initializes the shared fission bank
-	init_shared_fission_bank(shared_fission_bank_size);
 
   // Allocate array for matching filter bins
   #pragma omp parallel
@@ -687,11 +71,9 @@ int openmc_simulation_init()
   allocate_banks();
 
   // Allocate tally results arrays if they're not allocated yet
-	  std::cout << "Trying to Initialize Results..." << std::endl;
   for (auto& t : model::tallies) {
     t->init_results();
   }
-	  std::cout << "Success in Initializing Results!" << std::endl;
 
   // Set up material nuclide index mapping
   for (auto& mat : model::materials) {
@@ -736,8 +118,6 @@ int openmc_simulation_finalize()
 
   // Skip if simulation was never run
   if (!simulation::initialized) return 0;
-  
-  free_shared_fission_bank();
 
   // Stop active batch timer and start finalization timer
   simulation::time_active.stop();
@@ -808,14 +188,6 @@ int openmc_next_batch(int* status)
     // ====================================================================
     // LOOP OVER PARTICLES
 
-	simulation::current_work = 1;
-
-    //#pragma omp parallel
-    {
-      transport();
-    }
-
-/*
     #pragma omp parallel for schedule(runtime)
     for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
       simulation::current_work = i_work;
@@ -827,7 +199,6 @@ int openmc_next_batch(int* status)
       // transport particle
       p.transport();
     }
-*/
 
     // Accumulate time for transport
     simulation::time_transport.stop();
@@ -884,9 +255,7 @@ int restart_batch;
 bool satisfy_triggers {false};
 int total_gen {0};
 double total_weight;
-int64_t thread_work_index;
 int64_t work_per_rank;
-int64_t work_per_thread;
 
 const RegularMesh* entropy_mesh {nullptr};
 const RegularMesh* ufs_mesh {nullptr};
@@ -1070,15 +439,7 @@ void finalize_generation()
   if (settings::run_mode == RUN_MODE_EIGENVALUE) {
 #ifdef _OPENMP
     // Join the fission bank from each thread into one global fission bank
-    //join_bank_from_threads();
-	
-	  // We need to move all the stuff from the shared_fission_bank into the real one.
-		//std::vector<Particle::Bank> shared_fission_bank_vector(shared_fission_bank, shared_fission_bank + shared_fission_bank_length);
-      //simulation::fission_bank = shared_fission_bank_vector;
-	  std::cout << "Fission bank length = " << shared_fission_bank_length << std::endl;
-	  for( int i = 0; i < shared_fission_bank_length; i++ )
-		  simulation::fission_bank.push_back(shared_fission_bank[i]);
-	  shared_fission_bank_length = 0;
+    join_bank_from_threads();
 #endif
 
     // Distribute fission bank across processors evenly
@@ -1106,10 +467,8 @@ void finalize_generation()
 
 void initialize_history(Particle* p, int64_t index_source)
 {
-	//assert(index_source - 1 >= 0 );
   // set defaults
-  //p->from_source(&simulation::source_bank[index_source - 1]);
-  p->from_source(&simulation::source_bank[index_source]);
+  p->from_source(&simulation::source_bank[index_source - 1]);
 
   // set identifier for particle
   p->id_ = simulation::work_index[mpi::rank] + index_source;
@@ -1139,29 +498,6 @@ void initialize_history(Particle* p, int64_t index_source)
       }
     }
   }
-
-// Display message if high verbosity or trace is on
-   if (settings::verbosity >= 9 || simulation::trace) {
-      write_message("Simulating Particle " + std::to_string(p->id_));
-   }
-
-    // // Initialize number of events to zero
-   // int n_event = 0;
-
-    // Add paricle's starting weight to count for normalizing tallies later
-   #pragma omp atomic
-   simulation::total_weight += p->wgt_;
-
-    // Force calculation of cross-sections by setting last energy to zero
-   if (settings::run_CE) {
-     for (auto& micro : p->neutron_xs_) micro.last_E = 0.0;
-   }
-
-    // Prepare to write out particle track.
-   if (p->write_track_) add_particle_track();
-
-    // Every particle starts with no accumulated flux derivative.
-   if (!model::active_tallies.empty()) zero_flux_derivs();
 }
 
 int overall_generation()
@@ -1192,32 +528,6 @@ void calculate_work()
     i_bank += work_i;
     simulation::work_index[i + 1] = i_bank;
   }
-	#ifdef _OPENMP
-   // Determine work per thread
-   int remainder_thread = simulation::work_per_rank % omp_get_max_threads();
-
-    #pragma omp parallel
-   {
-     simulation::work_per_thread = simulation::work_per_rank / omp_get_num_threads();
-     if (omp_get_thread_num() < remainder_thread) {
-       ++simulation::work_per_thread;
-     }
-   }
-
-    int64_t work_i = 0;
-   #pragma omp parallel for ordered
-   for (int i = 0; i < omp_get_num_threads(); ++i) {
-     #pragma omp ordered
-     {
-       simulation::thread_work_index = work_i;
-       work_i += simulation::work_per_thread;
-	   /*
-       std::cout << "Thread " << omp_get_thread_num() << ": " << simulation::work_per_thread
-         << std::endl;
-		 */
-     }
-   }
- #endif
 }
 
 #ifdef OPENMC_MPI
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index bfd91df1d4..852918a0cd 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -804,7 +804,6 @@ Tally::init_triggers(pugi::xml_node node)
 void Tally::init_results()
 {
   int n_scores = scores_.size() * nuclides_.size();
-  std::cout << "TRYING TO ALLOCATE n filter bins, n scores =  " << n_filter_bins_ << " " <<  n_scores << std::endl;
   results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
 }