initial implementation wip

GuySten · GuySten · commit fd474019699b · 2025-08-26T04:50:21.000-07:00
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
@@ -62,7 +62,8 @@ class Particle : public ParticleData {
   //! site may have been produced from an external source, from fission, or
   //! simply as a secondary particle.
   //! \param src Source site data
-  void from_source(const SourceSite* src);
+  //! \param particle ParticleType incident particle type
+  void from_source(const SourceSite* src, ParticleType particle);
 
   // Coarse-grained particle events
   void event_calculate_xs();
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
@@ -469,6 +469,7 @@ class ParticleData : public GeometryState {
   CacheDataMG mg_xs_cache_;
 
   ParticleType type_ {ParticleType::neutron};
+  ParticleType type_last_ {ParticleType::neutron};
 
   double E_;
   double E_last_;
@@ -566,6 +567,8 @@ class ParticleData : public GeometryState {
   // Particle type (n, p, e, gamma, etc)
   ParticleType& type() { return type_; }
   const ParticleType& type() const { return type_; }
+  ParticleType& type_last() { return type_last_; }
+  const ParticleType& type_last() const { return type_last_; }
 
   // Current particle energy, energy before collision,
   // and corresponding multigroup group indices. Energy
diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h
@@ -39,6 +39,7 @@ enum class FilterType {
   MUSURFACE,
   PARENT_NUCLIDE,
   PARTICLE,
+  PARTICLE_OUT,
   POLAR,
   SPHERICAL_HARMONICS,
   SPATIAL_LEGENDRE,
diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h
@@ -48,5 +48,45 @@ class ParticleFilter : public Filter {
   vector<ParticleType> particles_;
 };
 
+//==============================================================================
+//! Bins by type of outgoing particle (e.g. neutron, photon).
+//==============================================================================
+
+class ParticleOutFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~ParticleOutFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "particleout"; }
+  FilterType type() const override { return FilterType::PARTICLE_OUT; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const vector<ParticleType>& particles() const { return particles_; }
+
+  void set_particles(span<ParticleType> particles);
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<ParticleType> particles_;
+};
+
 } // namespace openmc
 #endif // OPENMC_TALLIES_FILTER_PARTICLE_H
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
@@ -202,6 +202,7 @@ extern std::unordered_map<int, int> tally_map;
 extern vector<unique_ptr<Tally>> tallies;
 extern vector<int> active_tallies;
 extern vector<int> active_analog_tallies;
+extern vector<int> active_particleout_analog_tallies;
 extern vector<int> active_tracklength_tallies;
 extern vector<int> active_collision_tallies;
 extern vector<int> active_meshsurf_tallies;
diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
@@ -72,14 +72,16 @@ void score_collision_tally(Particle& p);
 //! Analog tallies are triggered at every collision, not every event.
 //
 //! \param p The particle being tracked
-void score_analog_tally_ce(Particle& p);
+//! \param tallies A vector of the indices of the tallies to score to
+void score_analog_tally_ce(Particle& p, const vector<int>& tallies);
 
 //! Score tallies based on a simple count of events (for multigroup).
 //
 //! Analog tallies are triggered at every collision, not every event.
 //
 //! \param p The particle being tracked
-void score_analog_tally_mg(Particle& p);
+//! \param tallies A vector of the indices of the tallies to score to
+void score_analog_tally_mg(Particle& p, const vector<int>& tallies);
 
 //! Score tallies using a tracklength estimate of the flux.
 //
diff --git a/openmc/filter.py b/openmc/filter.py
@@ -24,9 +24,9 @@
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
-    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface',
-    'meshmaterial',
+    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'particleout', 
+    'cellinstance', 'collision', 'time', 'parentnuclide', 'weight', 'meshborn',
+    'meshsurface', 'meshmaterial',
 )
 
 _CURRENT_NAMES = (
@@ -785,6 +785,29 @@ def from_xml_element(cls, elem, **kwargs):
         filter_id = int(get_text(elem, "id"))
         bins = get_elem_list(elem, "bins", str) or []
         return cls(bins, filter_id=filter_id)
+        
+        
+class ParticleoutFilter(ParticleFilter):
+    """Bins tally events based on the outgoing particle type.
+
+    Parameters
+    ----------
+    bins : str, or sequence of str
+        The particles to tally represented as strings ('neutron', 'photon',
+        'electron', 'positron').
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : sequence of str
+        The particles to tally
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """        
 
 
 class ParentNuclideFilter(ParticleFilter):
diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py
@@ -22,9 +22,9 @@
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
     'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
     'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter',
-    'ParentNuclideFilter', 'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
-    'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter',
-    'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
+    'ParentNuclideFilter', 'ParticleFilter', 'ParticleoutFilter', 'PolarFilter',
+    'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter',
+    'UniverseFilter', 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
 ]
 
 # Tally functions
@@ -564,6 +564,10 @@ def bins(self):
         return [ParticleType(i) for i in particle_i]
 
 
+class ParticleoutFilter(ParticleFilter):
+    filter_type = 'particleout'
+
+
 class PolarFilter(Filter):
     filter_type = 'polar'
 
diff --git a/src/mesh.cpp b/src/mesh.cpp
@@ -358,7 +358,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
           site.r[ax1] = min1 + (i1 + 0.5) * d1;
           site.r[ax2] = min2 + (i2 + 0.5) * d2;
 
-          p.from_source(&site);
+          p.from_source(&site, site.particle);
 
           // Determine particle's location
           if (!exhaustive_find_cell(p)) {
diff --git a/src/particle.cpp b/src/particle.cpp
@@ -110,7 +110,7 @@ void Particle::split(double wgt)
   }
 }
 
-void Particle::from_source(const SourceSite* src)
+void Particle::from_source(const SourceSite* src, ParticleType particle)
 {
   // Reset some attributes
   clear();
@@ -141,6 +141,7 @@ void Particle::from_source(const SourceSite* src)
     E() = data::mg.energy_bin_avg_[g()];
   }
   E_last() = E();
+  type_last() = particle;
   time() = src->time;
   time_last() = src->time;
   parent_nuclide() = src->parent_nuclide;
@@ -160,6 +161,7 @@ void Particle::event_calculate_xs()
   // Store pre-collision particle properties
   wgt_last() = wgt();
   E_last() = E();
+  type_last() = type();
   u_last() = u();
   r_last() = r();
   time_last() = time();
@@ -348,9 +350,9 @@ void Particle::event_collide()
     score_collision_tally(*this);
   if (!model::active_analog_tallies.empty()) {
     if (settings::run_CE) {
-      score_analog_tally_ce(*this);
+      score_analog_tally_ce(*this, model::active_analog_tallies);
     } else {
-      score_analog_tally_mg(*this);
+      score_analog_tally_mg(*this, model::active_analog_tallies);
     }
   }
 
@@ -408,6 +410,8 @@ void Particle::event_revive_from_secondary()
     wgt() = 0.0;
   }
 
+  auto type = this->type();
+
   // Check for secondary particles if this particle is dead
   if (!alive()) {
     // Write final position for this particle
@@ -419,11 +423,20 @@ void Particle::event_revive_from_secondary()
     if (secondary_bank().empty())
       return;
 
-    from_source(&secondary_bank().back());
+    from_source(&secondary_bank().back(), type);
     secondary_bank().pop_back();
     n_event() = 0;
     bank_second_E() = 0.0;
 
+    // Score tallies affected by secondary particles
+    if (!model::active_particleout_analog_tallies.empty()) {
+      if (settings::run_CE) {
+        score_analog_tally_ce(*this, model::active_particleout_analog_tallies);
+      } else {
+        score_analog_tally_mg(*this, model::active_particleout_analog_tallies);
+      }
+    }
+
     // Subtract secondary particle energy from interim pulse-height results
     if (!model::active_pulse_height_tallies.empty() &&
         this->type() == ParticleType::photon) {
diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp
@@ -64,6 +64,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   p.r_last() = p.r();
   p.u_last() = p.u();
   p.E_last() = p.E();
+  p.type_last() = p.type();
   p.g_last() = p.g();
   p.time_last() = p.time();
 
diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp
@@ -796,7 +796,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   }
 
   site.E = 0.0;
-  this->from_source(&site);
+  this->from_source(&site, site.particle);
 
   // Locate ray
   if (lowest_coord().cell() == C_NONE) {
diff --git a/src/simulation.cpp b/src/simulation.cpp
@@ -567,7 +567,8 @@ void initialize_history(Particle& p, int64_t index_source)
   // set defaults
   if (settings::run_mode == RunMode::EIGENVALUE) {
     // set defaults for eigenvalue simulations from primary bank
-    p.from_source(&simulation::source_bank[index_source - 1]);
+    p.from_source(
+      &simulation::source_bank[index_source - 1], ParticleType::neutron);
   } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
     // initialize random number seed
     int64_t id = (simulation::total_gen + overall_generation() - 1) *
@@ -576,7 +577,7 @@ void initialize_history(Particle& p, int64_t index_source)
     uint64_t seed = init_seed(id, STREAM_SOURCE);
     // sample from external source distribution or custom library then set
     auto site = sample_external_source(&seed);
-    p.from_source(&site);
+    p.from_source(&site, site.particle);
   }
   p.current_work() = index_source;
 
diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp
@@ -146,6 +146,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<ParentNuclideFilter>(id);
   } else if (type == "particle") {
     return Filter::create<ParticleFilter>(id);
+  } else if (type == "particleout") {
+    return Filter::create<ParticleOutFilter>(id);
   } else if (type == "polar") {
     return Filter::create<PolarFilter>(id);
   } else if (type == "surface") {
diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp
@@ -6,6 +6,10 @@
 
 namespace openmc {
 
+//==============================================================================
+// ParticleFilter implementation
+//==============================================================================
+
 void ParticleFilter::from_xml(pugi::xml_node node)
 {
   auto particles = get_node_array<std::string>(node, "bins");
@@ -34,8 +38,11 @@ void ParticleFilter::set_particles(span<ParticleType> particles)
 void ParticleFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
+  // Get the pre-collision type of the particle.
+  auto type = p.type_last();
+
   for (auto i = 0; i < particles_.size(); i++) {
-    if (particles_[i] == p.type()) {
+    if (particles_[i] == type) {
       match.bins_.push_back(i);
       match.weights_.push_back(1.0);
     }
@@ -58,6 +65,66 @@ std::string ParticleFilter::text_label(int bin) const
   return fmt::format("Particle: {}", particle_type_to_str(p));
 }
 
+//==============================================================================
+// ParticleOutFilter implementation
+//==============================================================================
+
+void ParticleOutFilter::from_xml(pugi::xml_node node)
+{
+  auto particles = get_node_array<std::string>(node, "bins");
+
+  // Convert to vector of ParticleType
+  vector<ParticleType> types;
+  for (auto& p : particles) {
+    types.push_back(str_to_particle_type(p));
+  }
+  this->set_particles(types);
+}
+
+void ParticleOutFilter::set_particles(span<ParticleType> particles)
+{
+  // Clear existing particles
+  particles_.clear();
+  particles_.reserve(particles.size());
+
+  // Set particles and number of bins
+  for (auto p : particles) {
+    particles_.push_back(p);
+  }
+  n_bins_ = particles_.size();
+}
+
+void ParticleOutFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  for (auto i = 0; i < particles_.size(); i++) {
+    if (particles_[i] == p.type()) {
+      match.bins_.push_back(i);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+void ParticleOutFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  vector<std::string> particles;
+  for (auto p : particles_) {
+    particles.push_back(particle_type_to_str(p));
+  }
+  write_dataset(filter_group, "bins", particles);
+}
+
+std::string ParticleOutFilter::text_label(int bin) const
+{
+  const auto& p = particles_.at(bin);
+  return fmt::format("ParticleOut: {}", particle_type_to_str(p));
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
 extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[])
 {
   if (int err = verify_filter(idx))
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
