Implement cell importance variance reduction scheme.

include/openmc/cell.h

@@ -228,6 +228,13 @@ class Cell {
   //! \return Density in [g/cm3]
   double density(int32_t instance = -1) const;
 
+  //! Get the importance of a cell instance
+  //! \param[in] index 0 for neutrons and 1 for photons.
+  //! \param[in] instance Instance index. If -1 is given, the importance for
+  //! the first instance is returned.
+  //! \return Importance
+  double importance(int index, int32_t instance = -1) const;
+
   //! Set the temperature of a cell instance
   //! \param[in] T Temperature in [K]
   //! \param[in] instance Instance index. If -1 is given, the temperature for
@@ -359,6 +366,11 @@ class Cell {
   //! \brief Unitless density multiplier(s) within this cell.
   vector<double> density_mult_;
 
+  //! \brief Importance for this cell per particle type
+  //!
+  //! May be multiple for distribcell.
+  vector<vector<double>> importance_ {{1.0}, {1.0}};
+
   //! \brief Neighboring cells in the same universe.
   NeighborList neighbors_;
 

include/openmc/geometry.h

@@ -6,6 +6,7 @@
 
 #include "openmc/array.h"
 #include "openmc/constants.h"
+#include "openmc/particle.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -52,6 +53,21 @@ bool check_cell_overlap(GeometryState& p, bool error = true);
 
 int cell_instance_at_level(const GeometryState& p, int level);
 
+//==============================================================================
+//! Get the cell importance for a particle at the specified
+//! type and universe level.
+//!
+//! \param p A particle for which to compute the instance using
+//!   its coordinates
+//! \param type ParticleType to get importance for
+//! \param level The level (zero indexed) of the geometry where the instance
+//! should be computed. \return The importance of the cell at the specified
+//! level.
+//==============================================================================
+
+double cell_importance_at_level(
+  const GeometryState& p, const ParticleType type, int level);
+
 //==============================================================================
 //! Locate a particle in the geometry tree and set its geometry data fields.
 //!

include/openmc/particle_data.h

@@ -501,6 +501,7 @@ class ParticleData : public GeometryState {
   int g_ {0};
   int g_last_;
 
+  double importance_last_ {1.0};
   double wgt_ {1.0};
   double wgt_born_ {1.0};
   double wgt_ww_born_ {-1.0};
@@ -610,6 +611,9 @@ class ParticleData : public GeometryState {
   int& g_last() { return g_last_; }
   const int& g_last() const { return g_last_; }
 
+  double& importance_last() { return importance_last_; }
+  const double& importance_last() const { return importance_last_; }
+
   // Statistic weight of particle. Setting to zero indicates that the particle
   // is dead.
   double& wgt() { return wgt_; }

include/openmc/simulation.h

@@ -22,6 +22,7 @@ constexpr int STATUS_EXIT_ON_TRIGGER {2};
 
 namespace simulation {
 
+extern bool cell_importances; //!< Use cell importances variance reduction
 extern int ct_current_file;   //!< current collision track file index
 extern "C" int current_batch; //!< current batch
 extern "C" int current_gen;   //!< current fission generation

include/openmc/source.h

@@ -91,7 +91,7 @@ class Source {
 
   // Methods for constraints
   void read_constraints(pugi::xml_node node);
-  bool satisfies_spatial_constraints(Position r) const;
+  bool satisfies_spatial_constraints(Position r, ParticleType type) const;
   bool satisfies_energy_constraints(double E) const;
   bool satisfies_time_constraints(double time) const;
 

openmc/cell.py

@@ -1,4 +1,4 @@
-from collections.abc import Iterable
+from collections.abc import Iterable, Mapping
 from math import cos, sin, pi
 from numbers import Real
 
@@ -69,6 +69,9 @@ class Cell(IDManagerMixin):
     rotation_matrix : numpy.ndarray
         The rotation matrix defined by the angles specified in the
         :attr:`Cell.rotation` property.
+    importance : dict
+        Dictionary containing importance of the cell for neutrons and photons. Defaults to 1.0.
+        Multiple importances can be given to give each distributed cell instance a unique importance.
     temperature : float or iterable of float
         Temperature of the cell in Kelvin.  Multiple temperatures can be given
         to give each distributed cell instance a unique temperature.
@@ -110,6 +113,7 @@ def __init__(self, cell_id=None, name='', fill=None, region=None):
         self.name = name
         self.fill = fill
         self.region = region
+        self._importance = {}
         self._rotation = None
         self._rotation_matrix = None
         self._temperature = None
@@ -238,6 +242,35 @@ def rotation(self, rotation):
     @property
     def rotation_matrix(self):
         return self._rotation_matrix
+
+    @property
+    def importance(self) -> dict:
+        return self._importance
+
+    @importance.setter
+    def importance(self, importance):
+        # Make sure importances are positive
+        cv.check_type('cell importance', importance, (Mapping, Iterable, Real))
+        if isinstance(importance, Mapping):
+            for key, value in importance.items():
+                cv.check_type('cell importance key', key, ('neutron', 'photon'))
+                cv.check_type('cell importance value', value, (Iterable, Real))
+                if isinstance(value, Iterable):
+                    cv.check_type('cell importance value', value, Iterable, Real)
+                    for i, val in enumerate(importance):
+                        cv.check_greater_than(f'cell importance value[{i}]', val, 0.0, True)
+                else:
+                    cv.check_greater_than('cell importance value', value, 0.0, True)
+            self._importance = importance
+
+        else:
+            if isinstance(importance, Iterable):
+                cv.check_type('cell importance', importance, Iterable, Real)
+                for i, val in enumerate(importance):
+                    cv.check_greater_than(f'cell importance[{i}]', val, 0.0, True)
+            else:
+                cv.check_greater_than('cell importance', importance, 0.0, True)
+            self._importance = {'neutron': importance, 'photon': importance}
 
     @property
     def temperature(self):
@@ -692,6 +725,14 @@ def create_surface_elements(node, element, memo=None):
             else:
                 element.set("density", str(self.density))
 
+        for key in ('neutron', 'photon'):
+            if self.importance.get(key, 1.0) != 1.0:
+                if isinstance(self.importance[key], Iterable):
+                    importance_subelement = ET.SubElement(element, f"importance_{key}")
+                    importance_subelement.text = ' '.join(str(val) for val in self.importance[key])
+                else:
+                    element.set("importance", str(self.importance[key]))
+
         if self.translation is not None:
             element.set("translation", ' '.join(map(str, self.translation)))
 
@@ -750,6 +791,16 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe):
         v = get_text(elem, 'volume')
         if v is not None:
             c.volume = float(v)
+        importance = {}
+
+        for key in ('neutron', 'photon'):
+            values = get_elem_list(elem, f"importance_{key}", float)
+            if values is not None:
+                if len(values) == 1:
+                    values = values[0]
+                importance[f"importance_{key}"] = values
+        setattr(c, "importance", importance)
+
         for key in ('temperature', 'density', 'rotation', 'translation'):
             values = get_elem_list(elem, key, float)
             if values is not None:

src/cell.cpp

@@ -22,6 +22,7 @@
 #include "openmc/material.h"
 #include "openmc/nuclide.h"
 #include "openmc/settings.h"
+#include "openmc/simulation.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -101,6 +102,15 @@ double Cell::temperature(int32_t instance) const
   }
 }
 
+double Cell::importance(int index, int32_t instance) const
+{
+  if (instance >= 0) {
+    return importance_[index][importance_[index].size() > 1 ? instance : 0];
+  } else {
+    return importance_[index][0];
+  }
+}
+
 double Cell::density_mult(int32_t instance) const
 {
   if (instance >= 0) {
@@ -337,6 +347,18 @@ void Cell::to_hdf5(hid_t cell_group) const
     write_dataset(group, "lattice", model::lattices[fill_]->id_);
   }
 
+  if (importance_[0].size() == 1) {
+    if (importance_[0][0] != 1.0)
+      write_dataset(group, "importance_neutron", importance_[0][0]);
+  } else {
+    write_dataset(group, "importance_neutron", importance_[0]);
+  }
+  if (importance_[1].size() == 1) {
+    if (importance_[1][0] != 1.0)
+      write_dataset(group, "importance_photon", importance_[1][0]);
+  } else {
+    write_dataset(group, "importance_photon", importance_[1]);
+  }
   close_group(group);
 }
 
@@ -472,6 +494,38 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
     }
   }
 
+  if (check_for_node(cell_node, "importance_neutron")) {
+    importance_[0] = get_node_array<double>(cell_node, "importance_neutron");
+    importance_[0].shrink_to_fit();
+
+    // Make sure all importancs are non-negative and greater than zero.
+    for (auto importance : importance_[0]) {
+      if (importance < 0)
+        fatal_error(fmt::format(
+          "Cell {} was specified with a negative neutron importance.", id_));
+      if (!settings::weight_windows_on) {
+        if (importance != 1.0 && importance > 0.0)
+          simulation::cell_importances = true;
+      }
+    }
+  }
+
+  if (check_for_node(cell_node, "importance_photon")) {
+    importance_[1] = get_node_array<double>(cell_node, "importance_photon");
+    importance_[1].shrink_to_fit();
+
+    // Make sure all importancs are non-negative and greater than zero.
+    for (auto importance : importance_[1]) {
+      if (importance < 0)
+        fatal_error(fmt::format(
+          "Cell {} was specified with a negative photon importance.", id_));
+      if (!settings::weight_windows_on) {
+        if (importance != 1.0 && importance > 0.0)
+          simulation::cell_importances = true;
+      }
+    }
+  }
+
   // Read the region specification.
   std::string region_spec;
   if (check_for_node(cell_node, "region")) {
@@ -1114,6 +1168,7 @@ vector<int32_t> Region::surfaces() const
 
 void read_cells(pugi::xml_node node)
 {
+  simulation::cell_importances = false;
   // Count the number of cells.
   int n_cells = 0;
   for (pugi::xml_node cell_node : node.children("cell")) {

src/geometry.cpp

@@ -97,6 +97,47 @@ int cell_instance_at_level(const GeometryState& p, int level)
   return instance;
 }
 
+double cell_importance_at_level(
+  const GeometryState& p, ParticleType type, int level)
+{
+  // throw error if the requested level is too deep for the geometry
+  if (level > model::n_coord_levels) {
+    fatal_error(
+      fmt::format("Cell importance at level {} requested, but only {} "
+                  "levels exist in the geometry.",
+        level, p.n_coord()));
+  }
+  int j = -1;
+  if (type.is_neutron())
+    j = 0;
+  else if (type.is_photon())
+    j = 1;
+  else
+    return 1.0;
+
+  // determine the cell instance
+  Cell& c {*model::cells[p.coord(level).cell()]};
+
+  // compute the cell's instance and importance
+  int instance = 0.0;
+  double importance = 1.0;
+  for (int i = 0; i < level; i++) {
+    const auto& c_i {*model::cells[p.coord(i).cell()]};
+    if (c_i.type_ == Fill::UNIVERSE) {
+      instance += c_i.offset_[c.distribcell_index_];
+    } else if (c_i.type_ == Fill::LATTICE) {
+      instance += c_i.offset_[c.distribcell_index_];
+      auto& lat {*model::lattices[p.coord(i + 1).lattice()]};
+      const auto& i_xyz {p.coord(i + 1).lattice_index()};
+      if (lat.are_valid_indices(i_xyz)) {
+        instance += lat.offset(c.distribcell_index_, i_xyz);
+      }
+    }
+    importance *= c_i.importance(j, instance);
+  }
+  return importance;
+}
+
 //==============================================================================
 
 bool find_cell_inner(

src/particle.cpp

@@ -298,6 +298,9 @@ void Particle::event_cross_surface()
   }
   n_coord_last() = n_coord();
 
+  if (simulation::cell_importances && material() != MATERIAL_VOID)
+    importance_last() = cell_importance_at_level(*this, type(), n_coord() - 1);
+
   // Set surface that particle is on and adjust coordinate levels
   surface() = boundary().surface();
   n_coord() = boundary().coord_level();
@@ -319,13 +322,47 @@ void Particle::event_cross_surface()
       add_surf_source_to_bank(*this, surf);
     }
     this->cross_surface(surf);
+    double importance = cell_importance_at_level(*this, type(), n_coord() - 1);
+    if (importance == 0.0) {
+      wgt() = 0.0;
+      return;
+    }
     // If no BC, add particle to surface source after crossing surface
     if (surf.surf_source_ && !surf.bc_) {
       add_surf_source_to_bank(*this, surf);
     }
     if (settings::weight_window_checkpoint_surface) {
       apply_weight_windows(*this);
     }
+    if (simulation::cell_importances && material() != MATERIAL_VOID) {
+      if (importance != importance_last()) {
+        if (importance < importance_last()) {
+          if (importance_last() * prn(current_seed()) < importance) {
+            wgt() *= importance_last() / importance;
+          } else {
+            wgt() = 0.;
+            return;
+          }
+        } else {
+          // do not further split the particle if above the limit
+          if (n_split() >= settings::max_history_splits)
+            return;
+          double n_s = importance / importance_last();
+          int n = static_cast<int>(n_s);
+          if (prn(current_seed()) <= n_s - n)
+            ++n;
+          n = std::min(n, settings::max_history_splits);
+          n_split() += n;
+
+          // Create secondaries and divide weight among all particles
+          for (int l = 0; l < n - 1; l++) {
+            split(wgt() / n);
+          }
+          // remaining weight is applied to current particle
+          wgt() /= n;
+        }
+      }
+    }
     event() = TallyEvent::SURFACE;
   }
   // Score cell to cell partial currents