added a data member in particle to track # of secondary particles

include/openmc/constants.h

@@ -105,9 +105,10 @@ constexpr std::array<const char*, 39> SUBSHELLS =  {
   "Q1", "Q2", "Q3"
 };
 
-// Void material
+// Void material and nuclide
 // TODO: refactor and remove
 constexpr int MATERIAL_VOID {-1};
+constexpr int NUCLIDE_VOID  {-1};
 
 // ============================================================================
 // CROSS SECTION RELATED CONSTANTS
@@ -127,6 +128,7 @@ constexpr int TEMPERATURE_INTERPOLATION {2};
 
 // Reaction types
 // TODO: Convert to enum
+constexpr int NONE_REACTION {0};
 constexpr int TOTAL_XS {1};
 constexpr int ELASTIC {2};
 constexpr int N_NONELASTIC {3};
@@ -344,6 +346,7 @@ constexpr int ESTIMATOR_COLLISION   {3};
 // TODO: Convert to enum
 constexpr int EVENT_SURFACE {-2};
 constexpr int EVENT_LATTICE {-1};
+constexpr int EVENT_KILL    {0};
 constexpr int EVENT_SCATTER {1};
 constexpr int EVENT_ABSORB  {2};
 

include/openmc/particle.h

@@ -184,7 +184,7 @@ public:
   //! \param u Direction of the secondary particle
   //! \param E Energy of the secondary particle in [eV]
   //! \param type Particle type
-  void create_secondary(Direction u, double E, Type type) const;
+  void create_secondary(Direction u, double E, Type type);
 
   //! initialize from a source site
   //
@@ -261,6 +261,7 @@ public:
 
   // Post-collision physical data
   int n_bank_ {0};        //!< number of fission sites banked
+  int n_bank_second_ {0}; //!< number of secondary particles
   double wgt_bank_ {0.0}; //!< weight of fission sites banked
   int n_delayed_bank_[MAX_DELAYED_GROUPS];  //!< number of delayed fission
                                             //!< sites banked

src/particle.cpp

@@ -73,7 +73,7 @@ Particle::clear()
 }
 
 void
-Particle::create_secondary(Direction u, double E, Type type) const
+Particle::create_secondary(Direction u, double E, Type type)
 {
   simulation::secondary_bank.emplace_back();
 
@@ -83,6 +83,8 @@ Particle::create_secondary(Direction u, double E, Type type) const
   bank.r = this->r();
   bank.u = u;
   bank.E = settings::run_CE ? E : g_;
+
+  n_bank_second_ += 1;
 }
 
 void
@@ -157,6 +159,11 @@ Particle::transport()
     u_last_ = this->u();
     r_last_ = this->r();
 
+    // Reset event variables
+    event_ = EVENT_KILL;
+    event_nuclide_ = NUCLIDE_VOID;
+    event_mt_ = NONE_REACTION;
+
     // 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
@@ -309,6 +316,7 @@ Particle::transport()
 
       // Reset banked weight during collision
       n_bank_ = 0;
+      n_bank_second_ = 0;
       wgt_bank_ = 0.0;
       for (int& v : n_delayed_bank_) v = 0;
 

src/physics.cpp

@@ -61,7 +61,9 @@ void collision(Particle* p)
   // Display information about collision
   if (settings::verbosity >= 10 || simulation::trace) {
     std::stringstream msg;
-    if (p->type_ == Particle::Type::neutron) {
+    if (p->event_ == EVENT_KILL) {
+      msg << "    " << " Killed. Energy = " << p->E_ << " eV.";
+    } else if (p->type_ == Particle::Type::neutron) {
       msg << "    " << reaction_name(p->event_mt_) << " with " <<
         data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV.";
     } else {
@@ -226,6 +228,7 @@ void sample_photon_reaction(Particle* p)
   if (prob > cutoff) {
     double mu = element.rayleigh_scatter(alpha);
     p->u() = rotate_angle(p->u(), mu, nullptr);
+    p->event_ = EVENT_SCATTER;
     p->event_mt_ = COHERENT;
     return;
   }
@@ -268,6 +271,7 @@ void sample_photon_reaction(Particle* p)
     phi += PI;
     p->E_ = alpha_out*MASS_ELECTRON_EV;
     p->u() = rotate_angle(p->u(), mu, &phi);
+    p->event_ = EVENT_SCATTER;
     p->event_mt_ = INCOHERENT;
     return;
   }
@@ -319,6 +323,7 @@ void sample_photon_reaction(Particle* p)
         // Allow electrons to fill orbital and produce auger electrons
         // and fluorescent photons
         element.atomic_relaxation(shell, *p);
+        p->event_ = EVENT_ABSORB;
         p->event_mt_ = 533 + shell.index_subshell;
         p->alive_ = false;
         p->E_ = 0.0;
@@ -344,6 +349,7 @@ void sample_photon_reaction(Particle* p)
     u = rotate_angle(p->u(), mu_positron, nullptr);
     p->create_secondary(u, E_positron, Particle::Type::positron);
 
+    p->event_ = EVENT_ABSORB;
     p->event_mt_ = PAIR_PROD;
     p->alive_ = false;
     p->E_ = 0.0;
@@ -361,6 +367,7 @@ void sample_electron_reaction(Particle* p)
 
   p->E_ = 0.0;
   p->alive_ = false;
+  p->event_ = EVENT_ABSORB;
 }
 
 void sample_positron_reaction(Particle* p)
@@ -386,6 +393,7 @@ void sample_positron_reaction(Particle* p)
 
   p->E_ = 0.0;
   p->alive_ = false;
+  p->event_ = EVENT_ABSORB;
 }
 
 int sample_nuclide(const Particle* p)