fixed the bug -- was using wrong index in one of the event loops (should have been j, but was using i). Now to remove all the debugging print statements

src/distribution_angle.cpp

@@ -85,6 +85,7 @@ double AngleDistribution::sample(double E) const
   // Sample between the ith and (i+1)th bin
   if (r > prn()) ++i;
 
+  //assert(i >= 0 );
   // Sample i-th distribution
   double mu = distribution_[i]->sample();
 

src/physics.cpp

@@ -79,6 +79,7 @@ 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);
 
@@ -127,9 +128,12 @@ 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_) {
@@ -598,6 +602,7 @@ 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);
 
@@ -610,6 +615,7 @@ 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;
@@ -617,6 +623,7 @@ void scatter(Particle* p, int i_nuclide)
   }
 
   if (!sampled) {
+	std::cout << "we are doing an Inelastic scatter" << std::endl;
     // =======================================================================
     // INELASTIC SCATTERING
 
@@ -643,7 +650,9 @@ 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_;
   }
 
@@ -1050,28 +1059,37 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
   // sample outgoing energy and scattering cosine
   double E;
   double mu;
+  // This gives a negative sometimes for some unknown reason
   rx->products_[0].sample(E_in, E, mu);
+std::cout << "E_in = " << E_in << "E = " << E << " mu = " << mu << std::endl;
 
   // if scattering system is in center-of-mass, transfer cosine of scattering
   // 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;
 

src/simulation.cpp

@@ -183,7 +183,8 @@ void process_calculate_xs_events(int * queue, int n)
   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[i];
+		Particle * p = particles + queue[j];
       if (p->material_ == MATERIAL_VOID) continue;
 
       // If material doesn't have this nuclide, skip it
@@ -370,6 +371,7 @@ 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) {
@@ -383,6 +385,8 @@ void process_collision_events()
 
     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;
@@ -393,6 +397,8 @@ void process_collision_events()
       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
@@ -404,6 +410,7 @@ void process_collision_events()
         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;
@@ -449,12 +456,82 @@ void process_collision_events()
     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++ )
+	{
+		Particle * p = particles + Q[i];
+		Particle p_true = *p;
+		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++ )
+	{
+		Particle * p = particles + Q[i];
+		Particle p_true = *p;
+		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++ )
+	{
+		Particle * p = particles + Q[i];
+		Particle p_true = *p;
+		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++ )
+	{
+		Particle * p = particles + Q[i];
+		Particle p_true = *p;
+		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++ )
+	{
+		Particle * p = particles + Q[i];
+		Particle p_true = *p;
+		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;
@@ -491,21 +568,35 @@ void transport()
 			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 << "Performing Fuel XS Lookups..." << std::endl;
+				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);
-				//std::cout << "Done." << std::endl;
 				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();
 			}
 		}