From aff6316642ccb3ba11b5c474c191f85ecb1ce436 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 5 Aug 2022 11:39:48 -0500 Subject: [PATCH 1/6] Simple rejection based on a vector of cell IDs --- include/openmc/source.h | 9 +++++++-- openmc/source.py | 29 +++++++++++++++++++++++++++- src/source.cpp | 42 +++++++++++++++++++++++++++++------------ 3 files changed, 65 insertions(+), 15 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 18fddc6899d..d48b4e4ae41 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -4,6 +4,8 @@ #ifndef OPENMC_SOURCE_H #define OPENMC_SOURCE_H +#include + #include "pugixml.hpp" #include "openmc/distribution_multi.h" @@ -51,13 +53,15 @@ class Source { }; //============================================================================== -//! Source composed of independent spatial, angle, energy, and time distributions +//! Source composed of independent spatial, angle, energy, and time +//! distributions //============================================================================== class IndependentSource : public Source { public: // Constructors - IndependentSource(UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time); + IndependentSource( + UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time); explicit IndependentSource(pugi::xml_node node); //! Sample from the external source distribution @@ -82,6 +86,7 @@ class IndependentSource : public Source { UPtrAngle angle_; //!< Angular distribution UPtrDist energy_; //!< Energy distribution UPtrDist time_; //!< Time distribution + std::unordered_set cells_; //!< Cells to reject from }; //============================================================================== diff --git a/openmc/source.py b/openmc/source.py index bd4c28ba462..6ec1338a382 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -5,6 +5,7 @@ import numpy as np import h5py +import openmc import openmc.checkvalue as cv from openmc.stats.multivariate import UnitSphere, Spatial from openmc.stats.univariate import Univariate @@ -36,6 +37,8 @@ class Source: Strength of the source particle : {'neutron', 'photon'} Source particle type + cells : iterable of int or openmc.Cell + Cells to reject based on Attributes ---------- @@ -57,11 +60,14 @@ class Source: Strength of the source particle : {'neutron', 'photon'} Source particle type + cells : iterable of int or openmc.Cell + Cells to reject based on """ def __init__(self, space=None, angle=None, energy=None, time=None, filename=None, - library=None, parameters=None, strength=1.0, particle='neutron'): + library=None, parameters=None, strength=1.0, particle='neutron', + cells=None): self._space = None self._angle = None self._energy = None @@ -69,6 +75,7 @@ def __init__(self, space=None, angle=None, energy=None, time=None, filename=None self._file = None self._library = None self._parameters = None + self._cells = [] if space is not None: self.space = space @@ -86,6 +93,8 @@ def __init__(self, space=None, angle=None, energy=None, time=None, filename=None self.parameters = parameters self.strength = strength self.particle = particle + if cells is not None: + self.cells = cells @property def file(self): @@ -123,6 +132,10 @@ def strength(self): def particle(self): return self._particle + @property + def cells(self): + return self._cells + @file.setter def file(self, filename): cv.check_type('source file', filename, str) @@ -169,6 +182,13 @@ def particle(self, particle): cv.check_value('source particle', particle, ['neutron', 'photon']) self._particle = particle + @cells.setter + def cells(self, cells): + self._cells = [ + cell.id if isinstance(cell, openmc.Cell) else cell + for cell in cells + ] + def to_xml_element(self): """Return XML representation of the source @@ -196,6 +216,9 @@ def to_xml_element(self): element.append(self.energy.to_xml_element('energy')) if self.time is not None: element.append(self.time.to_xml_element('time')) + if self.cells: + cells_elem = ET.SubElement(element, "cells") + cells_elem.text = ' '.join(str(x) for x in self.cells) return element @classmethod @@ -251,6 +274,10 @@ def from_xml_element(cls, elem): if time is not None: source.time = Univariate.from_xml_element(time) + cells = elem.find('cells') + if cells is not None: + source.cells = [int(x) for x in cells.text.split()] + return source diff --git a/src/source.cpp b/src/source.cpp index efc9e4598de..17bef4377c5 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -16,8 +16,10 @@ #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/cell.h" +#include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/file_utils.h" +#include "openmc/geometry.h" #include "openmc/hdf5_interface.h" #include "openmc/material.h" #include "openmc/memory.h" @@ -151,29 +153,36 @@ IndependentSource::IndependentSource(pugi::xml_node node) double p[] {1.0}; time_ = UPtrDist {new Discrete {T, p, 1}}; } + + // Check for cells to reject from + if (check_for_node(node, "cells")) { + auto cells = get_node_array(node, "cells"); + cells_.insert(cells.cbegin(), cells.cend()); + } } } SourceSite IndependentSource::sample(uint64_t* seed) const { SourceSite site; + site.particle = particle_; // Repeat sampling source location until a good site has been found bool found = false; int n_reject = 0; static int n_accept = 0; + while (!found) { // Set particle type - site.particle = particle_; + Particle p; + p.type() = particle_; + p.u() = {0.0, 0.0, 1.0}; // Sample spatial distribution - site.r = space_->sample(seed); + p.r() = space_->sample(seed); // Now search to see if location exists in geometry - int32_t cell_index, instance; - double xyz[] {site.r.x, site.r.y, site.r.z}; - int err = openmc_find_cell(xyz, &cell_index, &instance); - found = (err != OPENMC_E_GEOMETRY); + found = exhaustive_find_cell(p); // Check if spatial site is in fissionable material if (found) { @@ -181,15 +190,22 @@ SourceSite IndependentSource::sample(uint64_t* seed) const if (space_box) { if (space_box->only_fissionable()) { // Determine material - const auto& c = model::cells[cell_index]; - auto mat_index = - c->material_.size() == 1 ? c->material_[0] : c->material_[instance]; - + auto mat_index = p.material(); if (mat_index == MATERIAL_VOID) { found = false; } else { - if (!model::materials[mat_index]->fissionable_) - found = false; + found = model::materials[mat_index]->fissionable_; + } + } + } + + // Rejection based on cells + if (!cells_.empty()) { + found = false; + for (const auto& coord : p.coord()) { + if (contains(cells_, model::cells[coord.cell]->id_)) { + found = true; + break; } } } @@ -205,6 +221,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const "definition."); } } + + site.r = p.r(); } // Sample angle From 3e5c9d83c509754be32f40fba6c3f73427b2fd92 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 1 Sep 2022 16:13:26 -0500 Subject: [PATCH 2/6] Expand source rejection to cells, materials, and universes --- include/openmc/source.h | 7 +++- openmc/source.py | 85 +++++++++++++++++++++++++++++------------ src/source.cpp | 44 ++++++++++++++++----- 3 files changed, 100 insertions(+), 36 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index d48b4e4ae41..04b9ff85638 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -80,13 +80,18 @@ class IndependentSource : public Source { Distribution* time() const { return time_.get(); } private: + // Domain types + enum class DomainType { UNIVERSE, MATERIAL, CELL }; + + // Data members ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted double strength_ {1.0}; //!< Source strength UPtrSpace space_; //!< Spatial distribution UPtrAngle angle_; //!< Angular distribution UPtrDist energy_; //!< Energy distribution UPtrDist time_; //!< Time distribution - std::unordered_set cells_; //!< Cells to reject from + DomainType domain_type_; //!< Domain type for rejection + std::unordered_set domain_ids_; //!< Domains to reject from }; //============================================================================== diff --git a/openmc/source.py b/openmc/source.py index 6ec1338a382..48e14f6f710 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -1,5 +1,7 @@ +from collections.abc import Iterable from enum import Enum from numbers import Real +import warnings from xml.etree import ElementTree as ET import numpy as np @@ -37,8 +39,9 @@ class Source: Strength of the source particle : {'neutron', 'photon'} Source particle type - cells : iterable of int or openmc.Cell - Cells to reject based on + domains : iterable of openmc.Cell, openmc.Material, or openmc.Universe + Domains to reject based on, i.e., if a sampled spatial location is not + within one of these domains, it will be rejected. Attributes ---------- @@ -60,14 +63,16 @@ class Source: Strength of the source particle : {'neutron', 'photon'} Source particle type - cells : iterable of int or openmc.Cell - Cells to reject based on + ids : Iterable of int + IDs of domains to use for rejection + domain_type : {'cell', 'material', 'universe'} + Type of domain to use for rejection """ def __init__(self, space=None, angle=None, energy=None, time=None, filename=None, library=None, parameters=None, strength=1.0, particle='neutron', - cells=None): + domains=None): self._space = None self._angle = None self._energy = None @@ -75,7 +80,6 @@ def __init__(self, space=None, angle=None, energy=None, time=None, filename=None self._file = None self._library = None self._parameters = None - self._cells = [] if space is not None: self.space = space @@ -93,8 +97,17 @@ def __init__(self, space=None, angle=None, energy=None, time=None, filename=None self.parameters = parameters self.strength = strength self.particle = particle - if cells is not None: - self.cells = cells + + self._domain_ids = [] + self._domain_type = None + if domains is not None: + if isinstance(domains[0], openmc.Cell): + self.domain_type = 'cell' + elif isinstance(domains[0], openmc.Material): + self.domain_type = 'material' + elif isinstance(domains[0], openmc.Universe): + self.domain_type = 'universe' + self.domain_ids = [d.id for d in domains] @property def file(self): @@ -133,8 +146,22 @@ def particle(self): return self._particle @property - def cells(self): - return self._cells + def domain_ids(self): + return self._domain_ids + + @property + def domain_type(self): + return self._domain_type + + @domain_ids.setter + def domain_ids(self, ids): + cv.check_type('domain IDs', ids, Iterable, Real) + self._domain_ids = ids + + @domain_type.setter + def domain_type(self, domain_type): + cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) + self._domain_type = domain_type @file.setter def file(self, filename): @@ -182,13 +209,6 @@ def particle(self, particle): cv.check_value('source particle', particle, ['neutron', 'photon']) self._particle = particle - @cells.setter - def cells(self, cells): - self._cells = [ - cell.id if isinstance(cell, openmc.Cell) else cell - for cell in cells - ] - def to_xml_element(self): """Return XML representation of the source @@ -216,9 +236,11 @@ def to_xml_element(self): element.append(self.energy.to_xml_element('energy')) if self.time is not None: element.append(self.time.to_xml_element('time')) - if self.cells: - cells_elem = ET.SubElement(element, "cells") - cells_elem.text = ' '.join(str(x) for x in self.cells) + if self.domain_ids: + dt_elem = ET.SubElement(element, "domain_type") + dt_elem.text = self.domain_type + id_elem = ET.SubElement(element, "domain_ids") + id_elem.text = ' '.join(str(uid) for uid in self.domain_ids) return element @classmethod @@ -236,7 +258,24 @@ def from_xml_element(cls, elem): Source generated from XML element """ - source = cls() + domain_type = get_text(elem, "domain_type") + if domain_type is not None: + domain_ids = [int(x) for x in get_text(elem, "domain_ids").split()] + + # Instantiate some throw-away domains that are used by the + # constructor to assign IDs + with warnings.catch_warnings(): + warnings.simplefilter('ignore', openmc.IDWarning) + if domain_type == 'cell': + domains = [openmc.Cell(uid) for uid in domain_ids] + elif domain_type == 'material': + domains = [openmc.Material(uid) for uid in domain_ids] + elif domain_type == 'universe': + domains = [openmc.Universe(uid) for uid in domain_ids] + else: + domains = None + + source = cls(domains=domains) strength = get_text(elem, 'strength') if strength is not None: @@ -274,10 +313,6 @@ def from_xml_element(cls, elem): if time is not None: source.time = Univariate.from_xml_element(time) - cells = elem.find('cells') - if cells is not None: - source.cells = [int(x) for x in cells.text.split()] - return source diff --git a/src/source.cpp b/src/source.cpp index 17bef4377c5..e443c49e6b5 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -154,10 +154,22 @@ IndependentSource::IndependentSource(pugi::xml_node node) time_ = UPtrDist {new Discrete {T, p, 1}}; } - // Check for cells to reject from - if (check_for_node(node, "cells")) { - auto cells = get_node_array(node, "cells"); - cells_.insert(cells.cbegin(), cells.cend()); + // Check for domains to reject from + if (check_for_node(node, "domain_type")) { + std::string domain_type = get_node_value(node, "domain_type"); + if (domain_type == "cell") { + domain_type_ = DomainType::CELL; + } else if (domain_type == "material") { + domain_type_ = DomainType::MATERIAL; + } else if (domain_type == "universe") { + domain_type_ = DomainType::UNIVERSE; + } else { + fatal_error(std::string( + "Unrecognized domain type for source rejection: " + domain_type)); + } + + auto ids = get_node_array(node, "domain_ids"); + domain_ids_.insert(ids.begin(), ids.end()); } } } @@ -199,13 +211,25 @@ SourceSite IndependentSource::sample(uint64_t* seed) const } } - // Rejection based on cells - if (!cells_.empty()) { + // Rejection based on cells/materials/universes + if (!domain_ids_.empty()) { found = false; - for (const auto& coord : p.coord()) { - if (contains(cells_, model::cells[coord.cell]->id_)) { - found = true; - break; + if (domain_type_ == DomainType::MATERIAL) { + auto mat_index = p.material(); + if (mat_index != MATERIAL_VOID) { + if (contains(domain_ids_, model::materials[mat_index]->id())) { + found = true; + } + } + } else { + for (const auto& coord : p.coord()) { + auto id = (domain_type_ == DomainType::CELL) + ? model::cells[coord.cell]->id_ + : model::universes[coord.universe]->id_; + if (contains(domain_ids_, id)) { + found = true; + break; + } } } } From dd7a25c570ba1741f683d4add9a2ce96e7b88959 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 26 Sep 2022 14:20:34 -0500 Subject: [PATCH 3/6] Add test for source rejection by domain --- tests/unit_tests/test_source.py | 36 +++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 76dcbcdc8eb..7558d823819 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -1,6 +1,7 @@ from math import pi import openmc +import openmc.lib import openmc.stats import numpy as np from pytest import approx @@ -96,3 +97,38 @@ def test_source_xml_roundtrip(): np.testing.assert_allclose(new_src.angle.reference_uvw, src.angle.reference_uvw) assert new_src.particle == src.particle assert new_src.strength == approx(src.strength) + + +def test_rejection(run_in_tmpdir): + # Model with two spheres inside a box + mat = openmc.Material() + mat.add_nuclide('H1', 1.0) + sph1 = openmc.Sphere(x0=3, r=1.0) + sph2 = openmc.Sphere(x0=-3, r=1.0) + cube = openmc.model.RectangularParallelepiped( + -5., 5., -5., 5., -5., 5., boundary_type='reflective' + ) + cell1 = openmc.Cell(fill=mat, region=-sph1) + cell2 = openmc.Cell(fill=mat, region=-sph2) + cell3 = openmc.Cell(region=+sph1 & +sph2 & -cube) + model = openmc.Model() + model.geometry = openmc.Geometry([cell1, cell2, cell3]) + model.settings.particles = 100 + model.settings.batches = 10 + model.settings.run_mode = 'fixed source' + + # Set up a box source with rejection on the spherical cell + space = openmc.stats.Box(*cell3.bounding_box) + model.settings.source = openmc.Source(space=space, domains=[cell1, cell2]) + + # Load up model via openmc.lib and sample source + model.export_to_xml() + openmc.lib.init() + particles = openmc.lib.sample_external_source(1000) + + # Make sure that all sampled sources are within one of the spheres + joint_region = cell1.region | cell2.region + for p in particles: + assert p.r in joint_region + + openmc.lib.finalize() From e0be9de08625b00fe0f764a505612a5268c1ecec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 28 Sep 2022 09:33:48 -0500 Subject: [PATCH 4/6] Apply @shimwell suggestions from code review Co-authored-by: Jonathan Shimwell --- tests/unit_tests/test_source.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 7558d823819..ca3ee7f551f 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -110,7 +110,8 @@ def test_rejection(run_in_tmpdir): ) cell1 = openmc.Cell(fill=mat, region=-sph1) cell2 = openmc.Cell(fill=mat, region=-sph2) - cell3 = openmc.Cell(region=+sph1 & +sph2 & -cube) + non_source_region = +sph1 & +sph2 & -cube + cell3 = openmc.Cell(region=non_source_region) model = openmc.Model() model.geometry = openmc.Geometry([cell1, cell2, cell3]) model.settings.particles = 100 @@ -130,5 +131,6 @@ def test_rejection(run_in_tmpdir): joint_region = cell1.region | cell2.region for p in particles: assert p.r in joint_region + assert p.r not in non_source_region openmc.lib.finalize() From 2b401bc1f1229b0b594bd9e60900134934e3d3ec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 28 Sep 2022 10:26:13 -0500 Subject: [PATCH 5/6] Update tests/unit_tests/test_source.py Co-authored-by: Jonathan Shimwell --- tests/unit_tests/test_source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index ca3ee7f551f..1e4766b4560 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -131,6 +131,6 @@ def test_rejection(run_in_tmpdir): joint_region = cell1.region | cell2.region for p in particles: assert p.r in joint_region - assert p.r not in non_source_region + assert p.r not in non_source_region openmc.lib.finalize() From 5d441251fc9b924ff0038290e63ffb34d5d9a560 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 30 Sep 2022 13:56:01 -0500 Subject: [PATCH 6/6] Apply @pshriwise suggestions from code review Co-authored-by: Patrick Shriwise --- src/source.cpp | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index e443c49e6b5..fc50115cd82 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -217,19 +217,14 @@ SourceSite IndependentSource::sample(uint64_t* seed) const if (domain_type_ == DomainType::MATERIAL) { auto mat_index = p.material(); if (mat_index != MATERIAL_VOID) { - if (contains(domain_ids_, model::materials[mat_index]->id())) { - found = true; - } + found = contains(domain_ids_, model::materials[mat_index]->id()); } } else { for (const auto& coord : p.coord()) { auto id = (domain_type_ == DomainType::CELL) ? model::cells[coord.cell]->id_ : model::universes[coord.universe]->id_; - if (contains(domain_ids_, id)) { - found = true; - break; - } + if (found = contains(domain_ids_, id)) break; } } }