Implementation of forced collision variance reduction scheme

include/openmc/cell.h

@@ -228,6 +228,20 @@ class Cell {
   //! \return Density in [g/cm3]
   double density(int32_t instance = -1) const;
 
+  //! Determine forced collision flag corresponding
+  //! to a specific cell instance, taking into account
+  //! presence of distribcell
+  //! \param[in] instance of the cell
+  //! \return is forced collision enabled
+  bool forced_collision(int32_t instance) const
+  {
+    if (forced_collision_.size() > 1) {
+      return forced_collision_[instance];
+    } else {
+      return forced_collision_[0];
+    }
+  }
+
   //! 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
@@ -345,6 +359,11 @@ class Cell {
   //! \brief Index corresponding to this cell in distribcell arrays
   int distribcell_index_ {C_NONE};
 
+  //! \brief Enable forced collision within this cell.
+  //!
+  //! May be multiple for distribcell.
+  vector<bool> forced_collision_ {false};
+
   //! \brief Material(s) within this cell.
   //!
   //! May be multiple materials for distribcell.

include/openmc/simulation.h

@@ -25,6 +25,7 @@ namespace simulation {
 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
+extern bool forced_collision; //!< is forced collision used in the problem
 extern "C" bool initialized;  //!< has simulation been initialized?
 extern "C" double keff;       //!< average k over batches
 extern "C" double keff_std;   //!< standard deviation of average k

include/openmc/xml_interface.h

@@ -40,6 +40,26 @@ vector<T> get_node_array(
   return values;
 }
 
+template<>
+inline vector<bool> get_node_array<bool>(
+  pugi::xml_node node, const char* name, bool lowercase)
+{
+  // Get value of node attribute/child
+  std::string s {get_node_value(node, name, lowercase)};
+
+  // Read values one by one into vector
+  std::stringstream iss {s};
+  std::string value;
+  char first;
+  vector<bool> values;
+  while (iss >> value) {
+    first = value[0];
+    values.push_back((first == '1' || first == 't' || first == 'T' ||
+                      first == 'y' || first == 'Y'));
+  }
+  return values;
+}
+
 template<typename T>
 tensor::Tensor<T> get_node_tensor(
   pugi::xml_node node, const char* name, bool lowercase = false)

openmc/cell.py

@@ -76,6 +76,10 @@ class Cell(IDManagerMixin):
         Density of the cell in [g/cm3]. Multiple densities can be given to give
         each distributed cell instance a unique density. Densities set here will
         override the density set on materials used to fill the cell.
+    forced_collision : bool or iterable of bool
+        Forced collision flag in the cell. Multiple flags can be given to give
+        each distributed cell instance a unique forced collision flag.
+        Defaults to False.
     translation : Iterable of float
         If the cell is filled with a universe, this array specifies a vector
         that is used to translate (shift) the universe.
@@ -114,6 +118,7 @@ def __init__(self, cell_id=None, name='', fill=None, region=None):
         self._rotation_matrix = None
         self._temperature = None
         self._density = None
+        self._forced_collision = False
         self._translation = None
         self._paths = None
         self._num_instances = None
@@ -287,6 +292,17 @@ def density(self, density):
         else:
             self._density = density
 
+    @property
+    def forced_collision(self):
+        return self._forced_collision
+
+    @forced_collision.setter
+    def forced_collision(self, forced_collision):
+        cv.check_type('cell forced collision flag', forced_collision, (Iterable, bool))
+        if isinstance(forced_collision, Iterable):
+            cv.check_type('cell forced collision flag', forced_collision, Iterable, bool)
+        self._forced_collision = forced_collision
+
     @property
     def translation(self):
         return self._translation
@@ -692,6 +708,15 @@ def create_surface_elements(node, element, memo=None):
             else:
                 element.set("density", str(self.density))
 
+        if self.forced_collision:
+            if isinstance(self.density, Iterable):
+                if any(self.forced_collision):
+                    forced_collision_subelement = ET.SubElement(element, "forced_collision")
+                    forced_collision_subelement.text =  ' '.join(str(f)
+                                                                 for f in self.forced_collision)
+            else:
+                element.set("forced_collision", str(self.forced_collision))                
+
         if self.translation is not None:
             element.set("translation", ' '.join(map(str, self.translation)))
 
@@ -747,6 +772,10 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe):
             c.region = Region.from_expression(region, surfaces)
 
         # Check for other attributes
+        forced_collision = get_elem_list(elem, 'forced_collision', bool)
+        if forced_collision:
+            c.forced_collision = (forced_collision if len(forced_collision) > 1
+                                  else forced_collision[0])
         v = get_text(elem, 'volume')
         if v is not None:
             c.volume = float(v)

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 {
@@ -472,6 +473,15 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
     }
   }
 
+  if (check_for_node(cell_node, "forced_collision")) {
+    forced_collision_ = get_node_array<bool>(cell_node, "forced_collision");
+    forced_collision_.shrink_to_fit();
+    for (auto flag : forced_collision_) {
+      if (flag)
+        simulation::forced_collision = true;
+    }
+  }
+
   // Read the region specification.
   std::string region_spec;
   if (check_for_node(cell_node, "region")) {
@@ -1114,6 +1124,8 @@ vector<int32_t> Region::surfaces() const
 
 void read_cells(pugi::xml_node node)
 {
+  simulation::forced_collision = false;
+
   // Count the number of cells.
   int n_cells = 0;
   for (pugi::xml_node cell_node : node.children("cell")) {

src/particle.cpp

@@ -22,6 +22,7 @@
 #include "openmc/photon.h"
 #include "openmc/physics.h"
 #include "openmc/physics_mg.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
@@ -239,28 +240,88 @@ void Particle::event_advance()
   // Find the distance to the nearest boundary
   boundary() = distance_to_boundary(*this);
 
+  bool forced_collision = false;
+
+  double speed = this->speed();
+  double time_cutoff = settings::time_cutoff[type().transport_index()];
+  double distance_cutoff =
+    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
+
   // Sample a distance to collision
   if (type() == ParticleType::electron() ||
       type() == ParticleType::positron()) {
     collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
   } else if (macro_xs().total == 0.0) {
     collision_distance() = INFINITY;
   } else {
-    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
+    if (simulation::forced_collision &&
+        boundary().distance() <= distance_cutoff &&
+        boundary().distance() > TINY_BIT) {
+      auto c_id = coord(n_coord() - 1).cell();
+      auto& c = model::cells[c_id];
+      if (c->forced_collision(cell_instance())) {
+        if (cell_last(n_coord_last() - 1) != c_id)
+          forced_collision = true;
+      }
+    }
+    double U;
+    if (forced_collision) {
+      U = uniform_distribution(
+        std::exp(-boundary().distance() * macro_xs().total), 1, current_seed());
+    } else {
+      U = prn(current_seed());
+    }
+    collision_distance() = -std::log(U) / macro_xs().total;
   }
 
-  double speed = this->speed();
-  double time_cutoff = settings::time_cutoff[type().transport_index()];
-  double distance_cutoff =
-    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
+  double distance;
+  double dt;
+
+  if (forced_collision) {
+    distance = boundary().distance() - TINY_BIT;
+    double uncollided_wgt = wgt() * std::exp(-distance * macro_xs().total);
+    double collided_wgt = wgt() * -expm1(-distance * macro_xs().total);
+    wgt() = uncollided_wgt;
+    dt = distance / speed;
+    this->move_distance(distance);
+    this->time() += dt;
+    this->lifetime() += dt;
+
+    // Score timed track-length tallies
+    if (!model::active_timed_tracklength_tallies.empty()) {
+      score_timed_tracklength_tally(*this, distance);
+    }
+
+    // Score track-length tallies
+    if (!model::active_tracklength_tallies.empty()) {
+      score_tracklength_tally(*this, distance);
+    }
+
+    // Score track-length estimate of k-eff
+    if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
+      keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
+    }
+
+    // Score flux derivative accumulators for differential tallies.
+    if (!model::active_tallies.empty()) {
+      score_track_derivative(*this, distance);
+    }
+
+    split(uncollided_wgt);
+
+    this->move_distance(-distance);
+    this->time() -= dt;
+    this->lifetime() -= dt;
+    wgt() = collided_wgt;
+  }
 
   // Select smaller of the three distances
-  double distance =
+  distance =
     std::min({boundary().distance(), collision_distance(), distance_cutoff});
 
   // Advance particle in space and time
   this->move_distance(distance);
-  double dt = distance / speed;
+  dt = distance / speed;
   this->time() += dt;
   this->lifetime() += dt;
 

src/simulation.cpp

@@ -301,6 +301,7 @@ namespace simulation {
 int ct_current_file;
 int current_batch;
 int current_gen;
+bool forced_collision {false};
 bool initialized {false};
 double keff {1.0};
 double keff_std;

tests/unit_tests/test_forced_collision.py

@@ -0,0 +1,32 @@
+import openmc
+
+
+def test_forced_collision(run_in_tmpdir):
+    m = openmc.Material()
+    m.add_element('Li', 1.0)
+    m.set_density('g/cm3', 1.0)
+
+    surf1 = openmc.Sphere(r=99.9)
+    surf2 = openmc.Sphere(r=100, boundary_type='vacuum')
+
+    cell1 = openmc.Cell(region=-surf1)
+    cell2 = openmc.Cell(fill=m, region=-surf2 & +surf1)
+    cell2.forced_collision = True
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell1, cell2])
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        energy=openmc.stats.Discrete([5.0e6], [1.0]),
+        particle='photon',
+    )
+    model.settings.particles = 10
+    model.settings.batches = 1
+
+    tally = openmc.Tally()
+    tally.filters = [openmc.CellFilter([cell2])]
+    tally.scores = ['heating']
+    model.tallies = openmc.Tallies([tally])
+    model.run(apply_tally_results=True)
+
+    assert (tally.mean > 0.0).all(), "No heating detected"