From b8f73b3496baac2fac65ae18faac1f6fec925170 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 3 Mar 2026 17:31:25 +0200 Subject: [PATCH 1/7] initial implementation --- include/openmc/cell.h | 19 ++++++++++ include/openmc/simulation.h | 1 + src/cell.cpp | 12 ++++++ src/particle.cpp | 73 +++++++++++++++++++++++++++++++++---- src/simulation.cpp | 1 + 5 files changed, 99 insertions(+), 7 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 26c34ebd41d..a35fb746f93 100644 --- a/include/openmc/cell.h +++ b/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 forced_collision_ {false}; + //! \brief Material(s) within this cell. //! //! May be multiple materials for distribcell. diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..6029793c64a 100644 --- a/include/openmc/simulation.h +++ b/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 diff --git a/src/cell.cpp b/src/cell.cpp index ebe28c3d2ce..99f226255cc 100644 --- a/src/cell.cpp +++ b/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(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 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")) { diff --git a/src/particle.cpp b/src/particle.cpp index 8d0e5b3f06e..c40cd0c5753 100644 --- a/src/particle.cpp +++ b/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,6 +240,13 @@ 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()) { @@ -246,21 +254,72 @@ void Particle::event_advance() } 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) { + auto& c = model::cells[coord(n_coord() - 1).cell()]; + if (c->forced_collision(cell_instance())) { + auto& cl = model::cells[cell_last(n_coord_last() - 1)]; + if (c != cl) + 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(); + 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(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; diff --git a/src/simulation.cpp b/src/simulation.cpp index d0d6f037f9e..910f5b5b455 100644 --- a/src/simulation.cpp +++ b/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; From 5b1ac46271a396ef9b31f4ba05c4557e2ddbcb90 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 3 Mar 2026 17:50:43 +0200 Subject: [PATCH 2/7] add python API --- openmc/cell.py | 40 +++++++++++++++++++++++++++++++--------- 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/openmc/cell.py b/openmc/cell.py index 945dff0dbd6..89fa32997a8 100644 --- a/openmc/cell.py +++ b/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,15 +772,10 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe): c.region = Region.from_expression(region, surfaces) # Check for other attributes - temperature = get_elem_list(elem, 'temperature', float) - if temperature is not None: - if len(temperature) > 1: - c.temperature = temperature - else: - c.temperature = temperature[0] - density = get_elem_list(elem, 'density', float) - if density is not None: - c.density = density if len(density) > 1 else density[0] + 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) @@ -764,6 +784,8 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe): if values is not None: if key == 'rotation' and len(values) == 9: values = np.array(values).reshape(3, 3) + elif len(values) == 1: + values = values[0] setattr(c, key, values) # Add this cell to appropriate universe From 9929ac09e021afc1bce5b2dd7fb6649d90ede137 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 3 Mar 2026 18:28:36 +0200 Subject: [PATCH 3/7] add test --- tests/unit_tests/test_forced_collision.py | 32 +++++++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 tests/unit_tests/test_forced_collision.py diff --git a/tests/unit_tests/test_forced_collision.py b/tests/unit_tests/test_forced_collision.py new file mode 100644 index 00000000000..c395b9025ea --- /dev/null +++ b/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" From e1c5a0b602d69157593f321123e21c80032f687e Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 3 Mar 2026 19:47:08 +0200 Subject: [PATCH 4/7] fix bug --- include/openmc/xml_interface.h | 20 ++++++++++++++++++++ src/particle.cpp | 2 +- 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/include/openmc/xml_interface.h b/include/openmc/xml_interface.h index 17a34e5c77a..b0c257213bc 100644 --- a/include/openmc/xml_interface.h +++ b/include/openmc/xml_interface.h @@ -40,6 +40,26 @@ vector get_node_array( return values; } +template<> +inline vector get_node_array( + 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 values; + while (iss >> value) { + first = value[0]; + values.push_back((first == '1' || first == 't' || first == 'T' || + first == 'y' || first == 'Y')); + } + return values; +} + template tensor::Tensor get_node_tensor( pugi::xml_node node, const char* name, bool lowercase = false) diff --git a/src/particle.cpp b/src/particle.cpp index c40cd0c5753..94e4201c7f0 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -305,7 +305,7 @@ void Particle::event_advance() if (!model::active_tallies.empty()) { score_track_derivative(*this, distance); } - split(wgt()); + split(uncollided_wgt); this->move_distance(-distance); this->time() -= dt; From 390974b3ca04b6e508e22068f7d9f457219cd9d7 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 3 Mar 2026 21:27:11 +0200 Subject: [PATCH 5/7] fix bug --- src/particle.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 94e4201c7f0..8db14aea7ab 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -255,11 +255,12 @@ void Particle::event_advance() collision_distance() = INFINITY; } else { if (simulation::forced_collision && - boundary().distance() <= distance_cutoff) { - auto& c = model::cells[coord(n_coord() - 1).cell()]; + 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())) { - auto& cl = model::cells[cell_last(n_coord_last() - 1)]; - if (c != cl) + if (cell_last(n_coord_last() - 1) != c_id) forced_collision = true; } } @@ -277,7 +278,7 @@ void Particle::event_advance() double dt; if (forced_collision) { - distance = boundary().distance(); + 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; @@ -305,6 +306,7 @@ void Particle::event_advance() if (!model::active_tallies.empty()) { score_track_derivative(*this, distance); } + split(uncollided_wgt); this->move_distance(-distance); @@ -410,6 +412,11 @@ void Particle::event_collide() // Clear surface component surface() = SURFACE_NONE; + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + if (settings::run_CE) { collision(*this); } else { From 195674bb97244378c4158c6234daf6da13892a1c Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 4 Mar 2026 04:49:04 +0200 Subject: [PATCH 6/7] fix test --- tests/unit_tests/test_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 3846ba4fb8d..f65b384ee56 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -265,7 +265,7 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm): # Check to see that values are assigned to the C and python representations # First python cell = model.geometry.get_all_cells()[1] - assert cell.temperature == [600.0] + assert cell.temperature == 600.0 assert cell.density == [pytest.approx(10.0, 1e-5)] assert cell.fill.get_mass_density() == pytest.approx(5.0) # Now C From 9853f1f4b56338899b0ec36d6918423bd46d2087 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Wed, 4 Mar 2026 15:11:36 +0200 Subject: [PATCH 7/7] Refactor particle.cpp by removing unnecessary loop Removed redundant loop for updating last cell coordinates. --- src/particle.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 2e04fefe305..bda1661033b 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -413,11 +413,6 @@ void Particle::event_collide() // Clear surface component surface() = SURFACE_NONE; - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); - if (settings::run_CE) { collision(*this); } else {