From 3df948b42f4f63e1ada7b9a95bdbc0b1724888be Mon Sep 17 00:00:00 2001 From: erkn Date: Mon, 11 Mar 2024 01:18:25 +0100 Subject: [PATCH 01/57] add code to reject on domains also for file sources --- include/openmc/source.h | 5 +++++ src/source.cpp | 49 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+) diff --git a/include/openmc/source.h b/include/openmc/source.h index d5ea9bd1d56..d2d80808bf0 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -111,7 +111,12 @@ class FileSource : public Source { void load_sites_from_file( const std::string& path); //!< Load source sites from file private: + // Domain types + enum class DomainType { UNIVERSE, MATERIAL, CELL }; + vector sites_; //!< Source sites from a file + DomainType domain_type_; //!< Domain type for rejection + std::unordered_set domain_ids_; //!< Domains to reject from }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index 778962667bb..6f9a6f42f70 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -291,6 +291,55 @@ FileSource::FileSource(pugi::xml_node node) } else { this->load_sites_from_file(path); } + + // 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()); + } + + //_sites contains the full set of particles - loop through to reject + + if (!domain_ids_.empty()) { + bool found = false; + for (auto site = sites_.begin(); site!=sites_.end();) { + found = false; + Particle p; + p.r() = site->r; + if (domain_type_ == DomainType::MATERIAL) { + auto mat_index = p.material(); + if (mat_index != MATERIAL_VOID) { + found = contains(domain_ids_, model::materials[mat_index]->id()); + } + } else { + for (int i = 0; i < p.n_coord(); i++) { + auto id = (domain_type_ == DomainType::CELL) + ? model::cells[p.coord(i).cell]->id_ + : model::universes[p.coord(i).universe]->id_; + if ((found = contains(domain_ids_, id))) + break; + } + } + + if (found) { + site=sites_.erase(site); + } else { + ++site; + } + } + } } FileSource::FileSource(const std::string& path) From 9e850f66d2faf4781126b3c04e70bf4b721a2e3f Mon Sep 17 00:00:00 2001 From: erkn Date: Mon, 11 Mar 2024 08:31:36 +0100 Subject: [PATCH 02/57] also add the domain restriction to python layer --- openmc/source.py | 66 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 65 insertions(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index 441bf4e76b3..9a7b4032fbc 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -657,6 +657,9 @@ class FileSource(SourceBase): Path to the source file from which sites should be sampled strength : float Strength of the source (default is 1.0) + 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 ---------- @@ -667,8 +670,16 @@ class FileSource(SourceBase): type : str Indicator of source type: 'file' + .. versionadded:: x.y.z + + 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, path: Optional[PathLike] = None, strength=1.0) -> None: + + def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None) -> None: super().__init__(strength=strength) self._path = None @@ -676,6 +687,17 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0) -> None: if path is not None: self.path = path + 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 type(self) -> str: return "file" @@ -689,6 +711,24 @@ def path(self, p: PathLike): cv.check_type('source file', p, str) self._path = p + @property + def domain_ids(self): + return self._domain_ids + + @domain_ids.setter + def domain_ids(self, ids): + cv.check_type('domain IDs', ids, Iterable, Real) + self._domain_ids = ids + + @property + def domain_type(self): + return self._domain_type + + @domain_type.setter + def domain_type(self, domain_type): + cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) + self._domain_type = domain_type + def populate_xml_element(self, element): """Add necessary file source information to an XML element @@ -700,6 +740,11 @@ def populate_xml_element(self, element): """ if self.path is not None: element.set("file", self.path) + 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) @classmethod def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: @@ -728,6 +773,25 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: if strength is not None: source.strength = float(strength) + 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) + return source From 35e160c32d82df6794c04f8b101fee510f3f1a19 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 11 Mar 2024 14:16:18 +0100 Subject: [PATCH 03/57] move rejection to the sampling step it appears that before that geometry is not initialized --- src/source.cpp | 58 +++++++++++++++++++++++--------------------------- 1 file changed, 27 insertions(+), 31 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 6f9a6f42f70..809935e3ad1 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -310,36 +310,6 @@ FileSource::FileSource(pugi::xml_node node) domain_ids_.insert(ids.begin(), ids.end()); } - //_sites contains the full set of particles - loop through to reject - - if (!domain_ids_.empty()) { - bool found = false; - for (auto site = sites_.begin(); site!=sites_.end();) { - found = false; - Particle p; - p.r() = site->r; - if (domain_type_ == DomainType::MATERIAL) { - auto mat_index = p.material(); - if (mat_index != MATERIAL_VOID) { - found = contains(domain_ids_, model::materials[mat_index]->id()); - } - } else { - for (int i = 0; i < p.n_coord(); i++) { - auto id = (domain_type_ == DomainType::CELL) - ? model::cells[p.coord(i).cell]->id_ - : model::universes[p.coord(i).universe]->id_; - if ((found = contains(domain_ids_, id))) - break; - } - } - - if (found) { - site=sites_.erase(site); - } else { - ++site; - } - } - } } FileSource::FileSource(const std::string& path) @@ -377,7 +347,33 @@ void FileSource::load_sites_from_file(const std::string& path) SourceSite FileSource::sample(uint64_t* seed) const { - size_t i_site = sites_.size() * prn(seed); + //_sites contains the full set of particles - loop through to reject + bool found = false; + size_t i_site; + while (!found) { + i_site = sites_.size() * prn(seed); + Particle p; + p.r() = sites_[i_site].r; + found = exhaustive_find_cell(p); + if (found){ + if (domain_ids_.empty()) { + if (domain_type_ == DomainType::MATERIAL) { + auto mat_index = p.material(); + if (mat_index != MATERIAL_VOID) { + found = contains(domain_ids_, model::materials[mat_index]->id()); + } + } else { + for (int i = 0; i < p.n_coord(); i++) { + auto id = (domain_type_ == DomainType::CELL) + ? model::cells[p.coord(i).cell]->id_ + : model::universes[p.coord(i).universe]->id_; + if ((found = contains(domain_ids_, id))) + break; + } + } + } + } + } return sites_[i_site]; } From 0fbc48d013516076e123e9f3b0e712f88b9994d3 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 11 Mar 2024 14:28:35 +0100 Subject: [PATCH 04/57] improve commenting and formatting --- src/source.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 809935e3ad1..90388f675da 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -309,7 +309,6 @@ FileSource::FileSource(pugi::xml_node node) auto ids = get_node_array(node, "domain_ids"); domain_ids_.insert(ids.begin(), ids.end()); } - } FileSource::FileSource(const std::string& path) @@ -347,15 +346,20 @@ void FileSource::load_sites_from_file(const std::string& path) SourceSite FileSource::sample(uint64_t* seed) const { - //_sites contains the full set of particles - loop through to reject bool found = false; size_t i_site; + while (!found) { + // Sample a particle randomly from list i_site = sites_.size() * prn(seed); Particle p; p.r() = sites_[i_site].r; + + // Reject particle if it's not in the geometry at all found = exhaustive_find_cell(p); - if (found){ + + if (found) { + // Rejection based on cells/materials/universes if (domain_ids_.empty()) { if (domain_type_ == DomainType::MATERIAL) { auto mat_index = p.material(); @@ -365,8 +369,8 @@ SourceSite FileSource::sample(uint64_t* seed) const } else { for (int i = 0; i < p.n_coord(); i++) { auto id = (domain_type_ == DomainType::CELL) - ? model::cells[p.coord(i).cell]->id_ - : model::universes[p.coord(i).universe]->id_; + ? model::cells[p.coord(i).cell]->id_ + : model::universes[p.coord(i).universe]->id_; if ((found = contains(domain_ids_, id))) break; } From 471a4aa114fb3423b86e00c2b54c8097833ed150 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 11 Mar 2024 17:30:27 +0100 Subject: [PATCH 05/57] fix flipped logic --- src/source.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/source.cpp b/src/source.cpp index 90388f675da..fb350860a84 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -360,7 +360,7 @@ SourceSite FileSource::sample(uint64_t* seed) const if (found) { // Rejection based on cells/materials/universes - if (domain_ids_.empty()) { + if (!domain_ids_.empty()) { if (domain_type_ == DomainType::MATERIAL) { auto mat_index = p.material(); if (mat_index != MATERIAL_VOID) { From 972a2810fba3c5943f3f006cea300f0005e61869 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Tue, 12 Mar 2024 14:57:21 +0100 Subject: [PATCH 06/57] add option to kill non-conforming particles add an option to, if a particle is rejected, kill that one and move to next particle instead of resampling. This is less efficient, but avoids renormalization --- include/openmc/source.h | 2 ++ src/source.cpp | 27 +++++++++++++++++++++++++-- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index d2d80808bf0..427823b95fc 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -113,10 +113,12 @@ class FileSource : public Source { private: // Domain types enum class DomainType { UNIVERSE, MATERIAL, CELL }; + enum class RejectionStrategy { KILL, RESAMPLE }; vector sites_; //!< Source sites from a file DomainType domain_type_; //!< Domain type for rejection std::unordered_set domain_ids_; //!< Domains to reject from + RejectionStrategy rejection_strategy_; //!< Procedure for rejected }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index fb350860a84..db2e1194ac4 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -309,6 +309,22 @@ FileSource::FileSource(pugi::xml_node node) auto ids = get_node_array(node, "domain_ids"); domain_ids_.insert(ids.begin(), ids.end()); } + + // Check for how to handle rejected particles + if (check_for_node(node, "rejection_strategy")) { + std::string rejection_strategy=get_node_value(node, "rejection_strategy"); + if (rejection_strategy=="kill") { + rejection_strategy_ = RejectionStrategy::KILL; + } else if (rejection_strategy=="resample") { + rejection_strategy_ = RejectionStrategy::RESAMPLE; + } else { + fatal_error(std::string( + "Unrecognized strategy source rejection: " + rejection_strategy)); + } + } else { + // Default to kill rejected particles + rejection_strategy_ = RejectionStrategy::KILL; + } } FileSource::FileSource(const std::string& path) @@ -348,12 +364,14 @@ SourceSite FileSource::sample(uint64_t* seed) const { bool found = false; size_t i_site; + SourceSite site; while (!found) { // Sample a particle randomly from list i_site = sites_.size() * prn(seed); + site = sites_[i_site]; Particle p; - p.r() = sites_[i_site].r; + p.r() = site.r; // Reject particle if it's not in the geometry at all found = exhaustive_find_cell(p); @@ -377,8 +395,13 @@ SourceSite FileSource::sample(uint64_t* seed) const } } } + if (!found && rejection_strategy_==RejectionStrategy::KILL){ + // Accept particle regardless but set wgt=0 to trigger garbage collection + found = true; + site.wgt=0.0; + } } - return sites_[i_site]; + return site; } //============================================================================== From b1537a7cee1546346afb00aebf1b89d9e18b8bbe Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Wed, 13 Mar 2024 15:56:13 +0100 Subject: [PATCH 07/57] add option to reject by hypercube of particles --- include/openmc/source.h | 4 ++++ src/source.cpp | 45 ++++++++++++++++++++++++++++++++++------- 2 files changed, 42 insertions(+), 7 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 427823b95fc..d0611eb353a 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -118,7 +118,11 @@ class FileSource : public Source { vector sites_; //!< Source sites from a file DomainType domain_type_; //!< Domain type for rejection std::unordered_set domain_ids_; //!< Domains to reject from + vector lower_left_; //!< Lower left corner cds of filter hypercube + vector upper_right_; //!< Upper right corner cds of filter hypercube RejectionStrategy rejection_strategy_; //!< Procedure for rejected + + bool inside_hypercube(SourceSite& s) const; }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index db2e1194ac4..7ed53e5f884 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -301,6 +301,8 @@ FileSource::FileSource(pugi::xml_node node) domain_type_ = DomainType::MATERIAL; } else if (domain_type == "universe") { domain_type_ = DomainType::UNIVERSE; + } else if (domain_type == "filter") { + domain_type_ = DomainType::FILTER; } else { fatal_error(std::string( "Unrecognized domain type for source rejection: " + domain_type)); @@ -309,17 +311,29 @@ FileSource::FileSource(pugi::xml_node node) auto ids = get_node_array(node, "domain_ids"); domain_ids_.insert(ids.begin(), ids.end()); } + if (check_for_node(node, "lower_left")) { + auto ids = get_node_array(node, "lower_left"); + for (auto id : ids){ + lower_left_.push_back(id); + } + } + if (check_for_node(node, "upper_right")) { + auto ids = get_node_array(node, "upper_right"); + for (auto id : ids){ + upper_right_.push_back(id); + } + } // Check for how to handle rejected particles if (check_for_node(node, "rejection_strategy")) { - std::string rejection_strategy=get_node_value(node, "rejection_strategy"); - if (rejection_strategy=="kill") { + std::string rejection_strategy = get_node_value(node, "rejection_strategy"); + if (rejection_strategy == "kill") { rejection_strategy_ = RejectionStrategy::KILL; - } else if (rejection_strategy=="resample") { + } else if (rejection_strategy == "resample") { rejection_strategy_ = RejectionStrategy::RESAMPLE; } else { fatal_error(std::string( - "Unrecognized strategy source rejection: " + rejection_strategy)); + "Unrecognized strategy source rejection: " + rejection_strategy)); } } else { // Default to kill rejected particles @@ -360,6 +374,17 @@ void FileSource::load_sites_from_file(const std::string& path) file_close(file_id); } +bool FileSource::inside_hypercube(SourceSite& s) const +{ + std::vector cds { s.r[0], s.r[1], s.r[2], s.u[0], s.u[1], s.u[2], s.E, s.time }; + for (int i = 0; i < std::max(lower_left_.size(),upper_right_.size()); i++ ){ + if ( ( i < lower_left_.size() && cds[i] < lower_left_[i] ) || + ( i>upper_right_.size() && cds[i] > upper_right_[i] ) ) + return false; + } + return true; +} + SourceSite FileSource::sample(uint64_t* seed) const { bool found = false; @@ -376,8 +401,13 @@ SourceSite FileSource::sample(uint64_t* seed) const // Reject particle if it's not in the geometry at all found = exhaustive_find_cell(p); + // If not geometry rejected possibly reject by hypercube + if (found && ( !lower_left_.empty() || !upper_right_.empty() ) ) { + found = inside_hypercube(site); + } + + // Rejection based on cells/materials/universes if (found) { - // Rejection based on cells/materials/universes if (!domain_ids_.empty()) { if (domain_type_ == DomainType::MATERIAL) { auto mat_index = p.material(); @@ -395,10 +425,11 @@ SourceSite FileSource::sample(uint64_t* seed) const } } } - if (!found && rejection_strategy_==RejectionStrategy::KILL){ + + if (!found && rejection_strategy_ == RejectionStrategy::KILL) { // Accept particle regardless but set wgt=0 to trigger garbage collection found = true; - site.wgt=0.0; + site.wgt = 0.0; } } return site; From e12732be870498d864d52149e2e48ec4b3128261 Mon Sep 17 00:00:00 2001 From: erkn Date: Thu, 14 Mar 2024 03:49:49 +0100 Subject: [PATCH 08/57] allow hypercube description in python layer --- openmc/source.py | 42 +++++++++++++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 7 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 9a7b4032fbc..c4b82eec056 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -660,6 +660,12 @@ class FileSource(SourceBase): 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. + lower_left : Iterable of double + Coordinates of the lower left corner of a phase space hypercube from which to + accept particles. The dimensions are: x [cm],y [cm],z [cm], ux [ ],uy [ ],uz[ ], E [eV], t [s] + upper_right : Iterable of double + Coordinates of the upper right corner of a phase space hypercube from which to + accept particles. The dimensions are: x [cm],y [cm],z [cm], ux [ ],uy [ ],uz[ ], E [eV], t [s] Attributes ---------- @@ -676,10 +682,14 @@ class FileSource(SourceBase): IDs of domains to use for rejection domain_type : {'cell', 'material', 'universe'} Type of domain to use for rejection + lower_left : Iterable of double + Coordinates of the lower left corner of hypercube + upper_right : Iterable of double + Coordinates of the upper right corner of hypercube """ - def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None) -> None: + def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional{Sequence{Double]] = None) -> None: super().__init__(strength=strength) self._path = None @@ -697,6 +707,10 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optio elif isinstance(domains[0], openmc.Universe): self.domain_type = 'universe' self.domain_ids = [d.id for d in domains] + if lower_left is not None: + self.lower_left = lower_left + if upper_right is not None: + self.upper_right = lower_left @property def type(self) -> str: @@ -721,13 +735,21 @@ def domain_ids(self, ids): self._domain_ids = ids @property - def domain_type(self): - return self._domain_type + def lower_left(self): + return self._lower_left + + @lower_left.setter + def lower_left(self, lower_left): + self._lower_left = lower_left + + @property + def upper_right(self): + return self._upper_right + + @upper_right.setter + def domain_type(self, upper_right): + self._upper_right = upper_right - @domain_type.setter - def domain_type(self, domain_type): - cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) - self._domain_type = domain_type def populate_xml_element(self, element): """Add necessary file source information to an XML element @@ -745,6 +767,12 @@ def populate_xml_element(self, element): 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) + if self.lower_left_ is not none: + dt_elem = et.subelement(element, "lower_left") + dt_elem.text = self.lower_left + if self.upper_right_ is not none: + dt_elem = et.subelement(element, "upper_right") + dt_elem.text = self.upper_right @classmethod def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: From 03e122ce5cbe2400acd6d7615072a5c2d3e071e7 Mon Sep 17 00:00:00 2001 From: erkn Date: Thu, 14 Mar 2024 03:57:03 +0100 Subject: [PATCH 09/57] fix typos --- openmc/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index c4b82eec056..9f43f696903 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -689,7 +689,7 @@ class FileSource(SourceBase): """ - def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional{Sequence{Double]] = None) -> None: + def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional[Sequence[Double]] = None) -> None: super().__init__(strength=strength) self._path = None From 1baf0720356e1d5d532a625d5e3a912671c753bc Mon Sep 17 00:00:00 2001 From: erkn Date: Thu, 14 Mar 2024 03:57:03 +0100 Subject: [PATCH 10/57] fix typos --- openmc/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index c4b82eec056..9f43f696903 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -689,7 +689,7 @@ class FileSource(SourceBase): """ - def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional{Sequence{Double]] = None) -> None: + def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional[Sequence[Double]] = None) -> None: super().__init__(strength=strength) self._path = None From 916315ebef63d519f059c19c95c44efa3c6dcfd7 Mon Sep 17 00:00:00 2001 From: erkn Date: Thu, 14 Mar 2024 03:59:37 +0100 Subject: [PATCH 11/57] yet another misprint --- openmc/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index 9f43f696903..a216b43b8ab 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -710,7 +710,7 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optio if lower_left is not None: self.lower_left = lower_left if upper_right is not None: - self.upper_right = lower_left + self.upper_right = upper_right @property def type(self) -> str: From 672f4e606738d7a8afb656556c34d395eb01d0c3 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 14 Mar 2024 11:58:53 +0100 Subject: [PATCH 12/57] hypercube and rejection_strategy added to python also apply a bit of formatting --- include/openmc/source.h | 4 +-- openmc/source.py | 55 ++++++++++++++++++++++++++++++----------- src/source.cpp | 15 +++++------ 3 files changed, 51 insertions(+), 23 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index d0611eb353a..54a52a1569b 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -118,8 +118,8 @@ class FileSource : public Source { vector sites_; //!< Source sites from a file DomainType domain_type_; //!< Domain type for rejection std::unordered_set domain_ids_; //!< Domains to reject from - vector lower_left_; //!< Lower left corner cds of filter hypercube - vector upper_right_; //!< Upper right corner cds of filter hypercube + vector lower_left_; //!< Lower left corner cds of filter hypercube + vector upper_right_; //!< Upper right corner cds of filter hypercube RejectionStrategy rejection_strategy_; //!< Procedure for rejected bool inside_hypercube(SourceSite& s) const; diff --git a/openmc/source.py b/openmc/source.py index a216b43b8ab..7f79d3b62e1 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -660,11 +660,11 @@ class FileSource(SourceBase): 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. - lower_left : Iterable of double - Coordinates of the lower left corner of a phase space hypercube from which to + lower_left : Iterable of float + Lower-left corner coordinates of a phase space hypercube from which to accept particles. The dimensions are: x [cm],y [cm],z [cm], ux [ ],uy [ ],uz[ ], E [eV], t [s] - upper_right : Iterable of double - Coordinates of the upper right corner of a phase space hypercube from which to + upper_right : Iterable of float + Upper-right corner coordinates of a phase space hypercube from which to accept particles. The dimensions are: x [cm],y [cm],z [cm], ux [ ],uy [ ],uz[ ], E [eV], t [s] Attributes @@ -683,13 +683,13 @@ class FileSource(SourceBase): domain_type : {'cell', 'material', 'universe'} Type of domain to use for rejection lower_left : Iterable of double - Coordinates of the lower left corner of hypercube + Lower-left corner hypercube coordinates upper_right : Iterable of double - Coordinates of the upper right corner of hypercube + Upper-right corner hypercube coordinates """ - def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional[Sequence[Double]] = None) -> None: + def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_strategy = None, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional[Sequence[Double]] = None) -> None: super().__init__(strength=strength) self._path = None @@ -711,6 +711,8 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0, domains: Optio self.lower_left = lower_left if upper_right is not None: self.upper_right = upper_right + if rejection_strategy is not None: + self.rejection_strategy = rejection_strategy @property def type(self) -> str: @@ -734,12 +736,23 @@ def domain_ids(self, ids): cv.check_type('domain IDs', ids, Iterable, Real) self._domain_ids = ids + @property + def domain_type(self): + return self._domain_type + + @domain_type.setter + def domain_type(self, domain_type): + cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) + self._domain_type = domain_type + @property def lower_left(self): return self._lower_left @lower_left.setter def lower_left(self, lower_left): + name = 'lower-left hypercube coordinates' + cv.check_type(name, lower_left, Iterable, Real) self._lower_left = lower_left @property @@ -747,9 +760,20 @@ def upper_right(self): return self._upper_right @upper_right.setter - def domain_type(self, upper_right): + def upper_right(self, upper_right): + name = 'upper-right hypercube coordinates' + cv.check_type(name, upper_right, Iterable, Real) self._upper_right = upper_right + @property + def rejection_strategy(self): + return self._rejection_strategy + + @rejection_strategy.setter + def rejection_strategy(self, rejection_strategy): + name = 'rejection strategy' + cv.check_value('rejection strategy', rejection_strategy, ('resample','kill')) + self._rejection_strategy = rejection_strategy def populate_xml_element(self, element): """Add necessary file source information to an XML element @@ -767,12 +791,15 @@ def populate_xml_element(self, element): 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) - if self.lower_left_ is not none: - dt_elem = et.subelement(element, "lower_left") - dt_elem.text = self.lower_left - if self.upper_right_ is not none: - dt_elem = et.subelement(element, "upper_right") - dt_elem.text = self.upper_right + if self._lower_left is not None: + dt_elem = ET.SubElement(element, "lower_left") + dt_elem.text = ' '.join(str(cd) for cd in self.lower_left) + if self._upper_right is not None: + dt_elem = ET.SubElement(element, "upper_right") + dt_elem.text = ' '.join(str(cd) for cd in self.upper_right) + if self.rejection_strategy is not None: + dt_elem = ET.SubElement(element, "rejection_strategy") + dt_elem.text = self.rejection_strategy @classmethod def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: diff --git a/src/source.cpp b/src/source.cpp index 7ed53e5f884..73586daa436 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -313,13 +313,13 @@ FileSource::FileSource(pugi::xml_node node) } if (check_for_node(node, "lower_left")) { auto ids = get_node_array(node, "lower_left"); - for (auto id : ids){ + for (auto id : ids) { lower_left_.push_back(id); } } if (check_for_node(node, "upper_right")) { auto ids = get_node_array(node, "upper_right"); - for (auto id : ids){ + for (auto id : ids) { upper_right_.push_back(id); } } @@ -376,10 +376,11 @@ void FileSource::load_sites_from_file(const std::string& path) bool FileSource::inside_hypercube(SourceSite& s) const { - std::vector cds { s.r[0], s.r[1], s.r[2], s.u[0], s.u[1], s.u[2], s.E, s.time }; - for (int i = 0; i < std::max(lower_left_.size(),upper_right_.size()); i++ ){ - if ( ( i < lower_left_.size() && cds[i] < lower_left_[i] ) || - ( i>upper_right_.size() && cds[i] > upper_right_[i] ) ) + std::vector cds { + s.r[0], s.r[1], s.r[2], s.u[0], s.u[1], s.u[2], s.E, s.time}; + for (int i = 0; i < std::max(lower_left_.size(), upper_right_.size()); i++) { + if ((i < lower_left_.size() && cds[i] < lower_left_[i]) || + (i > upper_right_.size() && cds[i] > upper_right_[i])) return false; } return true; @@ -402,7 +403,7 @@ SourceSite FileSource::sample(uint64_t* seed) const found = exhaustive_find_cell(p); // If not geometry rejected possibly reject by hypercube - if (found && ( !lower_left_.empty() || !upper_right_.empty() ) ) { + if (found && (!lower_left_.empty() || !upper_right_.empty())) { found = inside_hypercube(site); } From ed469451630fb3b826e21cb87a4f27f058bc8d99 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 14 Mar 2024 15:44:30 +0100 Subject: [PATCH 13/57] adjust whitespace acc. to clang-format --- include/openmc/source.h | 8 ++++---- src/source.cpp | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 54a52a1569b..6d21a2ae89c 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -115,12 +115,12 @@ class FileSource : public Source { enum class DomainType { UNIVERSE, MATERIAL, CELL }; enum class RejectionStrategy { KILL, RESAMPLE }; - vector sites_; //!< Source sites from a file - DomainType domain_type_; //!< Domain type for rejection - std::unordered_set domain_ids_; //!< Domains to reject from + vector sites_; //!< Source sites from a file + DomainType domain_type_; //!< Domain type for rejection + std::unordered_set domain_ids_; //!< Domains to reject from vector lower_left_; //!< Lower left corner cds of filter hypercube vector upper_right_; //!< Upper right corner cds of filter hypercube - RejectionStrategy rejection_strategy_; //!< Procedure for rejected + RejectionStrategy rejection_strategy_; //!< Procedure for rejected bool inside_hypercube(SourceSite& s) const; }; diff --git a/src/source.cpp b/src/source.cpp index 73586daa436..30054eab53d 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -305,7 +305,7 @@ FileSource::FileSource(pugi::xml_node node) domain_type_ = DomainType::FILTER; } else { fatal_error(std::string( - "Unrecognized domain type for source rejection: " + domain_type)); + "Unrecognized domain type for source rejection: " + domain_type)); } auto ids = get_node_array(node, "domain_ids"); From 5a892f82b486fe53931d7ec34e38cfb6c537b191 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 14 Mar 2024 16:15:41 +0100 Subject: [PATCH 14/57] remove unimplemented FILTER option --- src/source.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 30054eab53d..60a36cb2b5e 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -301,8 +301,6 @@ FileSource::FileSource(pugi::xml_node node) domain_type_ = DomainType::MATERIAL; } else if (domain_type == "universe") { domain_type_ = DomainType::UNIVERSE; - } else if (domain_type == "filter") { - domain_type_ = DomainType::FILTER; } else { fatal_error(std::string( "Unrecognized domain type for source rejection: " + domain_type)); From eed4aaf621cb270510c7c0ed4192c0cead2a46fd Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 14 Mar 2024 17:03:15 +0100 Subject: [PATCH 15/57] fix var name --- openmc/source.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 7f79d3b62e1..f2ba96632ae 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -791,10 +791,10 @@ def populate_xml_element(self, element): 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) - if self._lower_left is not None: + if self.lower_left is not None: dt_elem = ET.SubElement(element, "lower_left") dt_elem.text = ' '.join(str(cd) for cd in self.lower_left) - if self._upper_right is not None: + if self.upper_right is not None: dt_elem = ET.SubElement(element, "upper_right") dt_elem.text = ' '.join(str(cd) for cd in self.upper_right) if self.rejection_strategy is not None: From bb8e892d8ab1fdbb29a93079cc7d72994fff8391 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 14 Mar 2024 17:49:30 +0100 Subject: [PATCH 16/57] initialize attributes to None --- openmc/source.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/source.py b/openmc/source.py index f2ba96632ae..e9447b0a47f 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -693,6 +693,8 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_stra super().__init__(strength=strength) self._path = None + self._lower_left = None + self._upper_right = None if path is not None: self.path = path From 23a29b41ca5dd9d854abd0ad8d09ef8d41437d95 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 14 Mar 2024 18:38:59 +0100 Subject: [PATCH 17/57] initialize the rej. stragtegy and add doc string --- openmc/source.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openmc/source.py b/openmc/source.py index e9447b0a47f..9e0376a58d1 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -686,6 +686,8 @@ class FileSource(SourceBase): Lower-left corner hypercube coordinates upper_right : Iterable of double Upper-right corner hypercube coordinates + rejection_strategy: {'resample','kill'} + Strategy when a particle is rejected. Pick a new particle or accept and terminate. """ @@ -695,6 +697,7 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_stra self._path = None self._lower_left = None self._upper_right = None + self._rejection_strategy = None if path is not None: self.path = path From 20b27f8bbeca1a80b2e119b6bdcfa7140f6f729a Mon Sep 17 00:00:00 2001 From: church89 Date: Tue, 19 Mar 2024 16:18:16 +0100 Subject: [PATCH 18/57] correct wrong condition --- src/source.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 60a36cb2b5e..122cc7bbe51 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -377,10 +377,10 @@ bool FileSource::inside_hypercube(SourceSite& s) const std::vector cds { s.r[0], s.r[1], s.r[2], s.u[0], s.u[1], s.u[2], s.E, s.time}; for (int i = 0; i < std::max(lower_left_.size(), upper_right_.size()); i++) { - if ((i < lower_left_.size() && cds[i] < lower_left_[i]) || - (i > upper_right_.size() && cds[i] > upper_right_[i])) + if ((cds[i] < lower_left_[i]) || (cds[i] > upper_right_[i])) return false; } + return true; } From d93ac83697e7fdb4fd4cc0f2705795398bf804bd Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Tue, 19 Mar 2024 23:16:26 +0100 Subject: [PATCH 19/57] move from a single vector to discrete bounds instead of a full coordinate vector instead use lower_left upper_right for spatial coordinates, and energy_bounds and time_bounds for energy and time. This should be more workable --- include/openmc/source.h | 9 ++++-- openmc/source.py | 61 +++++++++++++++++++++++++++++++++-------- src/source.cpp | 28 +++++++++++++------ 3 files changed, 75 insertions(+), 23 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 6d21a2ae89c..d5dd317615e 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -5,6 +5,7 @@ #define OPENMC_SOURCE_H #include +#include #include "pugixml.hpp" @@ -118,11 +119,13 @@ class FileSource : public Source { vector sites_; //!< Source sites from a file DomainType domain_type_; //!< Domain type for rejection std::unordered_set domain_ids_; //!< Domains to reject from - vector lower_left_; //!< Lower left corner cds of filter hypercube - vector upper_right_; //!< Upper right corner cds of filter hypercube + std::pair time_bounds_ {-DBL_MAX,DBL_MAX}; //!< time limits + std::pair energy_bounds_ {0,DBL_MAX}; //!< energy limits + vector lower_left_; //!< Lower left corner cds of filter + vector upper_right_; //!< Upper right corner cds of filter RejectionStrategy rejection_strategy_; //!< Procedure for rejected - bool inside_hypercube(SourceSite& s) const; + bool inside_bounds(SourceSite& s) const; }; //============================================================================== diff --git a/openmc/source.py b/openmc/source.py index 9e0376a58d1..2c0cc9b89c5 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -660,12 +660,14 @@ class FileSource(SourceBase): 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. + time_bounds : Iterable of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : Iterable of float + lower, upper bounds in energy [eV] of accepted particles lower_left : Iterable of float - Lower-left corner coordinates of a phase space hypercube from which to - accept particles. The dimensions are: x [cm],y [cm],z [cm], ux [ ],uy [ ],uz[ ], E [eV], t [s] + Lower-left corner coordinates from which to accept particles upper_right : Iterable of float - Upper-right corner coordinates of a phase space hypercube from which to - accept particles. The dimensions are: x [cm],y [cm],z [cm], ux [ ],uy [ ],uz[ ], E [eV], t [s] + Upper-right corner coordinates from which to accept particles Attributes ---------- @@ -675,26 +677,31 @@ class FileSource(SourceBase): Strength of the source type : str Indicator of source type: 'file' - - .. versionadded:: x.y.z - ids : Iterable of int IDs of domains to use for rejection domain_type : {'cell', 'material', 'universe'} Type of domain to use for rejection + domain_ids : iterable of int + Ids of domains to reject based on. + time_bounds : Iterable of float + lower, upper bounds in time of accepted particles + energy_bounds : Iterable of float + lower, upper bounds in energy of accepted particles lower_left : Iterable of double - Lower-left corner hypercube coordinates + Lower-left corner coordinates upper_right : Iterable of double - Upper-right corner hypercube coordinates + Upper-right corner coordinates rejection_strategy: {'resample','kill'} Strategy when a particle is rejected. Pick a new particle or accept and terminate. """ - def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_strategy = None, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional[Sequence[Double]] = None) -> None: + def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_strategy = None, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional[Sequence[Double]] = None, time_bounds: Optional[Sequence[Double]] = None, energy_bounds: Optional[Sequence[Double]] = None) -> None: super().__init__(strength=strength) self._path = None + self._time_bounds = None + self._energy_bounds = None self._lower_left = None self._upper_right = None self._rejection_strategy = None @@ -712,6 +719,10 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_stra elif isinstance(domains[0], openmc.Universe): self.domain_type = 'universe' self.domain_ids = [d.id for d in domains] + if time_bounds is not None: + self.time_bounds = time_bounds + if energy_bounds is not None: + self.energy_bounds = energy_bounds if lower_left is not None: self.lower_left = lower_left if upper_right is not None: @@ -750,13 +761,33 @@ def domain_type(self, domain_type): cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) self._domain_type = domain_type + @property + def time_bounds(self): + return self._time_bounds + + @time_bounds.setter + def time_bounds(self, time_bounds): + name = 'time bounds' + cv.check_type(name, time_bounds, Iterable, Real) + self._time_bounds = time_bounds + + @property + def energy_bounds(self): + return self._energy_bounds + + @energy_bounds.setter + def energy_bounds(self, energy_bounds): + name = 'energy bounds' + cv.check_type(name, energy_bounds, Iterable, Real) + self._energy_bounds = energy_bounds + @property def lower_left(self): return self._lower_left @lower_left.setter def lower_left(self, lower_left): - name = 'lower-left hypercube coordinates' + name = 'lower-left coordinates' cv.check_type(name, lower_left, Iterable, Real) self._lower_left = lower_left @@ -766,7 +797,7 @@ def upper_right(self): @upper_right.setter def upper_right(self, upper_right): - name = 'upper-right hypercube coordinates' + name = 'upper-right coordinates' cv.check_type(name, upper_right, Iterable, Real) self._upper_right = upper_right @@ -796,6 +827,12 @@ def populate_xml_element(self, element): 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) + if self.time_bounds is not None: + dt_elem = ET.SubElement(element, "time_bounds") + dt_elem.text = ' '.join(str(t) for t in self.time_bounds) + if self.energy_bounds is not None: + dt_elem = ET.SubElement(element, "energy_bounds") + dt_elem.text = ' '.join(str(E) for E in self.energy_bounds) if self.lower_left is not None: dt_elem = ET.SubElement(element, "lower_left") dt_elem.text = ' '.join(str(cd) for cd in self.lower_left) diff --git a/src/source.cpp b/src/source.cpp index 122cc7bbe51..2f664e5d3ba 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -309,6 +309,15 @@ FileSource::FileSource(pugi::xml_node node) auto ids = get_node_array(node, "domain_ids"); domain_ids_.insert(ids.begin(), ids.end()); } + + if (check_for_node(node, "time_bounds")) { + auto ids = get_node_array(node, "time_bounds"); + time_bounds_=std::make_pair(ids[0],ids[1]); + } + if (check_for_node(node, "energy_bounds")) { + auto ids = get_node_array(node, "energy_bounds"); + energy_bounds_=std::make_pair(ids[0],ids[1]); + } if (check_for_node(node, "lower_left")) { auto ids = get_node_array(node, "lower_left"); for (auto id : ids) { @@ -372,15 +381,18 @@ void FileSource::load_sites_from_file(const std::string& path) file_close(file_id); } -bool FileSource::inside_hypercube(SourceSite& s) const +bool FileSource::inside_bounds(SourceSite& s) const { - std::vector cds { - s.r[0], s.r[1], s.r[2], s.u[0], s.u[1], s.u[2], s.E, s.time}; - for (int i = 0; i < std::max(lower_left_.size(), upper_right_.size()); i++) { - if ((cds[i] < lower_left_[i]) || (cds[i] > upper_right_[i])) + if (s.E < energy_bounds_.first || s.E >energy_bounds_.second) { + return false; + } else if (s.time < time_bounds_.first || s.time > time_bounds_.second) { + return false; + } else { + if (!lower_left_.empty() && (s.r[0] < lower_left_[0] || s.r[1] < lower_left_[1] || s.r[2] < lower_left_[2])) + return false; + if (!upper_right_.empty() && (s.r[0] > upper_right_[0] || s.r[1] > upper_right_[1] || s.r[2] > upper_right_[2])) return false; } - return true; } @@ -401,8 +413,8 @@ SourceSite FileSource::sample(uint64_t* seed) const found = exhaustive_find_cell(p); // If not geometry rejected possibly reject by hypercube - if (found && (!lower_left_.empty() || !upper_right_.empty())) { - found = inside_hypercube(site); + if (found) { + found = inside_bounds(site); } // Rejection based on cells/materials/universes From 86ae0640b794c6232867a67330e30a2066002dad Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 13:50:39 +0200 Subject: [PATCH 20/57] source restriction is now an intermediate class The source may now be restricted if it inherits from RestrictedSource instead of directly from Source --- include/openmc/source.h | 48 +++++--- src/source.cpp | 261 +++++++++++++++++++--------------------- 2 files changed, 153 insertions(+), 156 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index d5dd317615e..6f01720c292 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -55,12 +55,38 @@ class Source { static unique_ptr create(pugi::xml_node node); }; +//============================================================================== +//! Abstract middle layer class providing source restrictions +//============================================================================== +class RestrictedSource : public Source { +protected: + // Domain types + enum class DomainType { UNIVERSE, MATERIAL, CELL }; + enum class RejectionStrategy { KILL, RESAMPLE }; + + // members + std::unordered_set domain_ids_; //!< Domains to reject from + DomainType domain_type_; //!< Domain type for rejection + std::pair time_bounds_ {-DBL_MAX,DBL_MAX}; //!< time limits + std::pair energy_bounds_ {0,DBL_MAX}; //!< energy limits + vector lower_left_; //!< Lower left corner cds of filter + vector upper_right_; //!< Upper right corner cds of filter + RejectionStrategy rejection_strategy_; //!< Procedure for rejected + + //methods + void check_for_restriction_nodes(pugi::xml_node node); + bool inside_bounds(SourceSite& s) const; + bool inside_spatial_bounds(SourceSite& s) const; + bool inside_energy_bounds(const double E) const; + bool inside_time_bounds(const double time) const; +}; + //============================================================================== //! Source composed of independent spatial, angle, energy, and time //! distributions //============================================================================== -class IndependentSource : public Source { +class IndependentSource : public RestrictedSource { public: // Constructors IndependentSource( @@ -83,9 +109,6 @@ 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 @@ -93,15 +116,13 @@ class IndependentSource : public Source { UPtrAngle angle_; //!< Angular distribution UPtrDist energy_; //!< Energy distribution UPtrDist time_; //!< Time distribution - DomainType domain_type_; //!< Domain type for rejection - std::unordered_set domain_ids_; //!< Domains to reject from }; //============================================================================== //! Source composed of particles read from a file //============================================================================== -class FileSource : public Source { +class FileSource : public RestrictedSource { public: // Constructors explicit FileSource(pugi::xml_node node); @@ -112,20 +133,7 @@ class FileSource : public Source { void load_sites_from_file( const std::string& path); //!< Load source sites from file private: - // Domain types - enum class DomainType { UNIVERSE, MATERIAL, CELL }; - enum class RejectionStrategy { KILL, RESAMPLE }; - vector sites_; //!< Source sites from a file - DomainType domain_type_; //!< Domain type for rejection - std::unordered_set domain_ids_; //!< Domains to reject from - std::pair time_bounds_ {-DBL_MAX,DBL_MAX}; //!< time limits - std::pair energy_bounds_ {0,DBL_MAX}; //!< energy limits - vector lower_left_; //!< Lower left corner cds of filter - vector upper_right_; //!< Upper right corner cds of filter - RejectionStrategy rejection_strategy_; //!< Procedure for rejected - - bool inside_bounds(SourceSite& s) const; }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index 2f664e5d3ba..f7fccf618a7 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -79,6 +79,123 @@ unique_ptr Source::create(pugi::xml_node node) } } +//============================================================================== +// RestrictedSource functionality +//============================================================================== +void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) { + // 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()); + } + + if (check_for_node(node, "time_bounds")) { + auto ids = get_node_array(node, "time_bounds"); + time_bounds_=std::make_pair(ids[0],ids[1]); + } + if (check_for_node(node, "energy_bounds")) { + auto ids = get_node_array(node, "energy_bounds"); + energy_bounds_=std::make_pair(ids[0],ids[1]); + } + if (check_for_node(node, "lower_left")) { + auto ids = get_node_array(node, "lower_left"); + for (auto id : ids) { + lower_left_.push_back(id); + } + } + if (check_for_node(node, "upper_right")) { + auto ids = get_node_array(node, "upper_right"); + for (auto id : ids) { + upper_right_.push_back(id); + } + } + + // Check for how to handle rejected particles + if (check_for_node(node, "rejection_strategy")) { + std::string rejection_strategy = get_node_value(node, "rejection_strategy"); + if (rejection_strategy == "kill") { + rejection_strategy_ = RejectionStrategy::KILL; + } else if (rejection_strategy == "resample") { + rejection_strategy_ = RejectionStrategy::RESAMPLE; + } else { + fatal_error(std::string( + "Unrecognized strategy source rejection: " + rejection_strategy)); + } + } else { + // Default to kill rejected particles + rejection_strategy_ = RejectionStrategy::RESAMPLE; + } +} + +bool RestrictedSource::inside_bounds(SourceSite& s) const +{ + if ( inside_spatial_bounds(s) and inside_energy_bounds(s.E) and inside_time_bounds(s.time) ) + return true; + return false; +} + +bool RestrictedSource::inside_energy_bounds(double E) const +{ + if (E < energy_bounds_.first || E >energy_bounds_.second) + return false; + return true; +} + +bool RestrictedSource::inside_time_bounds(const double time) const +{ + if (time < time_bounds_.first || time > time_bounds_.second) + return false; + return true; +} + +bool RestrictedSource::inside_spatial_bounds(SourceSite& s) const +{ + bool found = false; + Particle p; + p.r() = s.r; + + // Reject particle if it's not in the geometry at all + if (not (found = exhaustive_find_cell(p)) ) + return false; + + if (!lower_left_.empty() && (s.r[0] < lower_left_[0] || s.r[1] < lower_left_[1] || s.r[2] < lower_left_[2])) + return false; + if (!upper_right_.empty() && (s.r[0] > upper_right_[0] || s.r[1] > upper_right_[1] || s.r[2] > upper_right_[2])) + return false; + + if (!domain_ids_.empty()) { + if (domain_type_ == DomainType::MATERIAL) { + auto mat_index = p.material(); + if (mat_index != MATERIAL_VOID) { + found = contains(domain_ids_, model::materials[mat_index]->id()); + } + } else { + for (int i = 0; i < p.n_coord(); i++) { + auto id = (domain_type_ == DomainType::CELL) + ? model::cells[p.coord(i).cell]->id_ + : model::universes[p.coord(i).universe]->id_; + if ((found = contains(domain_ids_, id))) + break; + } + } + } + + //particle is inside all spatial bounds + return true; +} + //============================================================================== // IndependentSource implementation //============================================================================== @@ -149,23 +266,8 @@ IndependentSource::IndependentSource(pugi::xml_node node) time_ = UPtrDist {new Discrete {T, p, 1}}; } - // 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()); - } + // Check for additional defined restrictions (from RestrictedSource) + check_for_restriction_nodes(node); } } @@ -188,8 +290,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Sample spatial distribution p.r() = space_->sample(seed); - // Now search to see if location exists in geometry - found = exhaustive_find_cell(p); + // Check if otherwise outside defined restriction bounds + found=inside_spatial_bounds(site); // Check if spatial site is in fissionable material if (found) { @@ -205,25 +307,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const } } } - - // Rejection based on cells/materials/universes - if (!domain_ids_.empty()) { - found = false; - if (domain_type_ == DomainType::MATERIAL) { - auto mat_index = p.material(); - if (mat_index != MATERIAL_VOID) { - found = contains(domain_ids_, model::materials[mat_index]->id()); - } - } else { - for (int i = 0; i < p.n_coord(); i++) { - auto id = (domain_type_ == DomainType::CELL) - ? model::cells[p.coord(i).cell]->id_ - : model::universes[p.coord(i).universe]->id_; - if ((found = contains(domain_ids_, id))) - break; - } - } - } } // Check for rejection @@ -259,7 +342,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const site.E = energy_->sample(seed); // Resample if energy falls above maximum particle energy - if (site.E < data::energy_max[p]) + if (site.E < data::energy_max[p] or (inside_energy_bounds(site.E))) break; n_reject++; @@ -291,61 +374,7 @@ FileSource::FileSource(pugi::xml_node node) } else { this->load_sites_from_file(path); } - - // 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()); - } - - if (check_for_node(node, "time_bounds")) { - auto ids = get_node_array(node, "time_bounds"); - time_bounds_=std::make_pair(ids[0],ids[1]); - } - if (check_for_node(node, "energy_bounds")) { - auto ids = get_node_array(node, "energy_bounds"); - energy_bounds_=std::make_pair(ids[0],ids[1]); - } - if (check_for_node(node, "lower_left")) { - auto ids = get_node_array(node, "lower_left"); - for (auto id : ids) { - lower_left_.push_back(id); - } - } - if (check_for_node(node, "upper_right")) { - auto ids = get_node_array(node, "upper_right"); - for (auto id : ids) { - upper_right_.push_back(id); - } - } - - // Check for how to handle rejected particles - if (check_for_node(node, "rejection_strategy")) { - std::string rejection_strategy = get_node_value(node, "rejection_strategy"); - if (rejection_strategy == "kill") { - rejection_strategy_ = RejectionStrategy::KILL; - } else if (rejection_strategy == "resample") { - rejection_strategy_ = RejectionStrategy::RESAMPLE; - } else { - fatal_error(std::string( - "Unrecognized strategy source rejection: " + rejection_strategy)); - } - } else { - // Default to kill rejected particles - rejection_strategy_ = RejectionStrategy::KILL; - } + check_for_restriction_nodes(node); } FileSource::FileSource(const std::string& path) @@ -381,21 +410,6 @@ void FileSource::load_sites_from_file(const std::string& path) file_close(file_id); } -bool FileSource::inside_bounds(SourceSite& s) const -{ - if (s.E < energy_bounds_.first || s.E >energy_bounds_.second) { - return false; - } else if (s.time < time_bounds_.first || s.time > time_bounds_.second) { - return false; - } else { - if (!lower_left_.empty() && (s.r[0] < lower_left_[0] || s.r[1] < lower_left_[1] || s.r[2] < lower_left_[2])) - return false; - if (!upper_right_.empty() && (s.r[0] > upper_right_[0] || s.r[1] > upper_right_[1] || s.r[2] > upper_right_[2])) - return false; - } - return true; -} - SourceSite FileSource::sample(uint64_t* seed) const { bool found = false; @@ -409,36 +423,11 @@ SourceSite FileSource::sample(uint64_t* seed) const Particle p; p.r() = site.r; - // Reject particle if it's not in the geometry at all - found = exhaustive_find_cell(p); - - // If not geometry rejected possibly reject by hypercube - if (found) { - found = inside_bounds(site); - } - - // Rejection based on cells/materials/universes - if (found) { - if (!domain_ids_.empty()) { - if (domain_type_ == DomainType::MATERIAL) { - auto mat_index = p.material(); - if (mat_index != MATERIAL_VOID) { - found = contains(domain_ids_, model::materials[mat_index]->id()); - } - } else { - for (int i = 0; i < p.n_coord(); i++) { - auto id = (domain_type_ == DomainType::CELL) - ? model::cells[p.coord(i).cell]->id_ - : model::universes[p.coord(i).universe]->id_; - if ((found = contains(domain_ids_, id))) - break; - } - } - } - } + found = inside_bounds(site); if (!found && rejection_strategy_ == RejectionStrategy::KILL) { // Accept particle regardless but set wgt=0 to trigger garbage collection + // I.e. kill this particle immediately and pick the next found = true; site.wgt = 0.0; } From 78a87166a71d8436ae27200efa7f8de4619c6f0a Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 13:56:44 +0200 Subject: [PATCH 21/57] apply clang-format --- include/openmc/source.h | 14 +++++++------- src/source.cpp | 26 ++++++++++++++++---------- 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 6f01720c292..d8eb2509ae1 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -4,8 +4,8 @@ #ifndef OPENMC_SOURCE_H #define OPENMC_SOURCE_H -#include #include +#include #include "pugixml.hpp" @@ -67,13 +67,13 @@ class RestrictedSource : public Source { // members std::unordered_set domain_ids_; //!< Domains to reject from DomainType domain_type_; //!< Domain type for rejection - std::pair time_bounds_ {-DBL_MAX,DBL_MAX}; //!< time limits - std::pair energy_bounds_ {0,DBL_MAX}; //!< energy limits - vector lower_left_; //!< Lower left corner cds of filter - vector upper_right_; //!< Upper right corner cds of filter + std::pair time_bounds_ {-DBL_MAX, DBL_MAX}; //!< time limits + std::pair energy_bounds_ {0, DBL_MAX}; //!< energy limits + vector lower_left_; //!< Lower left corner cds of filter + vector upper_right_; //!< Upper right corner cds of filter RejectionStrategy rejection_strategy_; //!< Procedure for rejected - //methods + // methods void check_for_restriction_nodes(pugi::xml_node node); bool inside_bounds(SourceSite& s) const; bool inside_spatial_bounds(SourceSite& s) const; @@ -133,7 +133,7 @@ class FileSource : public RestrictedSource { void load_sites_from_file( const std::string& path); //!< Load source sites from file private: - vector sites_; //!< Source sites from a file + vector sites_; //!< Source sites from a file }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index f7fccf618a7..d1caf00abc1 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -82,7 +82,8 @@ unique_ptr Source::create(pugi::xml_node node) //============================================================================== // RestrictedSource functionality //============================================================================== -void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) { +void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) +{ // Check for domains to reject from if (check_for_node(node, "domain_type")) { std::string domain_type = get_node_value(node, "domain_type"); @@ -103,11 +104,11 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) { if (check_for_node(node, "time_bounds")) { auto ids = get_node_array(node, "time_bounds"); - time_bounds_=std::make_pair(ids[0],ids[1]); + time_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "energy_bounds")) { auto ids = get_node_array(node, "energy_bounds"); - energy_bounds_=std::make_pair(ids[0],ids[1]); + energy_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "lower_left")) { auto ids = get_node_array(node, "lower_left"); @@ -141,14 +142,15 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) { bool RestrictedSource::inside_bounds(SourceSite& s) const { - if ( inside_spatial_bounds(s) and inside_energy_bounds(s.E) and inside_time_bounds(s.time) ) + if (inside_spatial_bounds(s) and inside_energy_bounds(s.E) and + inside_time_bounds(s.time)) return true; return false; } bool RestrictedSource::inside_energy_bounds(double E) const { - if (E < energy_bounds_.first || E >energy_bounds_.second) + if (E < energy_bounds_.first || E > energy_bounds_.second) return false; return true; } @@ -167,12 +169,16 @@ bool RestrictedSource::inside_spatial_bounds(SourceSite& s) const p.r() = s.r; // Reject particle if it's not in the geometry at all - if (not (found = exhaustive_find_cell(p)) ) + if (not(found = exhaustive_find_cell(p))) return false; - if (!lower_left_.empty() && (s.r[0] < lower_left_[0] || s.r[1] < lower_left_[1] || s.r[2] < lower_left_[2])) + if (!lower_left_.empty() && + (s.r[0] < lower_left_[0] || s.r[1] < lower_left_[1] || + s.r[2] < lower_left_[2])) return false; - if (!upper_right_.empty() && (s.r[0] > upper_right_[0] || s.r[1] > upper_right_[1] || s.r[2] > upper_right_[2])) + if (!upper_right_.empty() && + (s.r[0] > upper_right_[0] || s.r[1] > upper_right_[1] || + s.r[2] > upper_right_[2])) return false; if (!domain_ids_.empty()) { @@ -192,7 +198,7 @@ bool RestrictedSource::inside_spatial_bounds(SourceSite& s) const } } - //particle is inside all spatial bounds + // particle is inside all spatial bounds return true; } @@ -291,7 +297,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const p.r() = space_->sample(seed); // Check if otherwise outside defined restriction bounds - found=inside_spatial_bounds(site); + found = inside_spatial_bounds(site); // Check if spatial site is in fissionable material if (found) { From c0a2f3141523a5c1c8810b9f7c4d38a88db78d57 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 22:56:24 +0200 Subject: [PATCH 22/57] use limits instead of cfloats --- include/openmc/source.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index d8eb2509ae1..eab443da3cf 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -4,7 +4,7 @@ #ifndef OPENMC_SOURCE_H #define OPENMC_SOURCE_H -#include +#include #include #include "pugixml.hpp" @@ -67,8 +67,8 @@ class RestrictedSource : public Source { // members std::unordered_set domain_ids_; //!< Domains to reject from DomainType domain_type_; //!< Domain type for rejection - std::pair time_bounds_ {-DBL_MAX, DBL_MAX}; //!< time limits - std::pair energy_bounds_ {0, DBL_MAX}; //!< energy limits + std::pair time_bounds_ {-std::numeric_limits::max(), std::numeric_limits::max()}; //!< time limits + std::pair energy_bounds_ {0, std::numeric_limits::max()}; //!< energy limits vector lower_left_; //!< Lower left corner cds of filter vector upper_right_; //!< Upper right corner cds of filter RejectionStrategy rejection_strategy_; //!< Procedure for rejected From 6dad9f837c166568bdc469593e2283d5a81c7810 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 22:58:12 +0200 Subject: [PATCH 23/57] check bounds for pair assignement --- src/source.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index d1caf00abc1..4601e36733e 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -104,11 +104,13 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) if (check_for_node(node, "time_bounds")) { auto ids = get_node_array(node, "time_bounds"); - time_bounds_ = std::make_pair(ids[0], ids[1]); + if (id.size()==2) + time_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "energy_bounds")) { auto ids = get_node_array(node, "energy_bounds"); - energy_bounds_ = std::make_pair(ids[0], ids[1]); + if (id.size()==2) + energy_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "lower_left")) { auto ids = get_node_array(node, "lower_left"); From 6831f7d4f7bafb048daf9bea04651a70a728b624 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 22:58:46 +0200 Subject: [PATCH 24/57] default to kill - not resample --- src/source.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/source.cpp b/src/source.cpp index 4601e36733e..8b0ec0bde76 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -138,7 +138,7 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) } } else { // Default to kill rejected particles - rejection_strategy_ = RejectionStrategy::RESAMPLE; + rejection_strategy_ = RejectionStrategy::KILL; } } From b45f343ca0e2203a87d0d82307e1758aac0c7439 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 23:00:26 +0200 Subject: [PATCH 25/57] better style (dont use and...) --- src/source.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 8b0ec0bde76..371d1b96cc6 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -144,10 +144,8 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) bool RestrictedSource::inside_bounds(SourceSite& s) const { - if (inside_spatial_bounds(s) and inside_energy_bounds(s.E) and - inside_time_bounds(s.time)) - return true; - return false; + return inside_spatial_bounds(s) && inside_energy_bounds(s.E) && + inside_time_bounds(s.time); } bool RestrictedSource::inside_energy_bounds(double E) const From a0cd9ecee8649327a422fe3d80f2d9ac75a447f0 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 23:01:58 +0200 Subject: [PATCH 26/57] prefer ! to not --- src/source.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/source.cpp b/src/source.cpp index 371d1b96cc6..f079870b71a 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -169,7 +169,7 @@ bool RestrictedSource::inside_spatial_bounds(SourceSite& s) const p.r() = s.r; // Reject particle if it's not in the geometry at all - if (not(found = exhaustive_find_cell(p))) + if (!(found = exhaustive_find_cell(p))) return false; if (!lower_left_.empty() && From 4671d90908d5f667a75c8fddd06e50e97f4f67c3 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 23:15:24 +0200 Subject: [PATCH 27/57] better style --- src/source.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index f079870b71a..2dc0d1ed645 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -150,16 +150,12 @@ bool RestrictedSource::inside_bounds(SourceSite& s) const bool RestrictedSource::inside_energy_bounds(double E) const { - if (E < energy_bounds_.first || E > energy_bounds_.second) - return false; - return true; + return E > energy_bounds_.first && E < energy_bounds_.second; } bool RestrictedSource::inside_time_bounds(const double time) const { - if (time < time_bounds_.first || time > time_bounds_.second) - return false; - return true; + return time > time_bounds_.first && time < time_bounds_.second; } bool RestrictedSource::inside_spatial_bounds(SourceSite& s) const From bd8978bf6f9b5996ddac08971b198da9e1ce3ac6 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 23:21:39 +0200 Subject: [PATCH 28/57] apply clang-format --- include/openmc/source.h | 8 +++++--- src/source.cpp | 6 +++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index eab443da3cf..41386002e13 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -67,9 +67,11 @@ class RestrictedSource : public Source { // members std::unordered_set domain_ids_; //!< Domains to reject from DomainType domain_type_; //!< Domain type for rejection - std::pair time_bounds_ {-std::numeric_limits::max(), std::numeric_limits::max()}; //!< time limits - std::pair energy_bounds_ {0, std::numeric_limits::max()}; //!< energy limits - vector lower_left_; //!< Lower left corner cds of filter + std::pair time_bounds_ {-std::numeric_limits::max(), + std::numeric_limits::max()}; //!< time limits + std::pair energy_bounds_ { + 0, std::numeric_limits::max()}; //!< energy limits + vector lower_left_; //!< Lower left corner cds of filter vector upper_right_; //!< Upper right corner cds of filter RejectionStrategy rejection_strategy_; //!< Procedure for rejected diff --git a/src/source.cpp b/src/source.cpp index 2dc0d1ed645..a1be0277ef7 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -104,12 +104,12 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) if (check_for_node(node, "time_bounds")) { auto ids = get_node_array(node, "time_bounds"); - if (id.size()==2) + if (id.size() == 2) time_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "energy_bounds")) { auto ids = get_node_array(node, "energy_bounds"); - if (id.size()==2) + if (id.size() == 2) energy_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "lower_left")) { @@ -145,7 +145,7 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) bool RestrictedSource::inside_bounds(SourceSite& s) const { return inside_spatial_bounds(s) && inside_energy_bounds(s.E) && - inside_time_bounds(s.time); + inside_time_bounds(s.time); } bool RestrictedSource::inside_energy_bounds(double E) const From e342243cf36b8230a36fe3fc3ea89058cc701bea Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Mon, 8 Apr 2024 23:44:12 +0200 Subject: [PATCH 29/57] fix typo --- src/source.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index a1be0277ef7..d36ca527e5b 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -104,12 +104,12 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) if (check_for_node(node, "time_bounds")) { auto ids = get_node_array(node, "time_bounds"); - if (id.size() == 2) + if (ids.size() == 2) time_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "energy_bounds")) { auto ids = get_node_array(node, "energy_bounds"); - if (id.size() == 2) + if (ids.size() == 2) energy_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "lower_left")) { From 6ede6ef95bb3fc078e01fff234ca5ebbdc4cb3e7 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Tue, 9 Apr 2024 01:11:02 +0200 Subject: [PATCH 30/57] remove restr. source class --- include/openmc/source.h | 9 ++------- src/source.cpp | 13 +++++-------- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 41386002e13..b6b359f99c6 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -53,12 +53,7 @@ class Source { virtual double strength() const { return 1.0; } static unique_ptr create(pugi::xml_node node); -}; -//============================================================================== -//! Abstract middle layer class providing source restrictions -//============================================================================== -class RestrictedSource : public Source { protected: // Domain types enum class DomainType { UNIVERSE, MATERIAL, CELL }; @@ -88,7 +83,7 @@ class RestrictedSource : public Source { //! distributions //============================================================================== -class IndependentSource : public RestrictedSource { +class IndependentSource : public Source { public: // Constructors IndependentSource( @@ -124,7 +119,7 @@ class IndependentSource : public RestrictedSource { //! Source composed of particles read from a file //============================================================================== -class FileSource : public RestrictedSource { +class FileSource : public Source { public: // Constructors explicit FileSource(pugi::xml_node node); diff --git a/src/source.cpp b/src/source.cpp index d36ca527e5b..a5476ec4f32 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -79,10 +79,7 @@ unique_ptr Source::create(pugi::xml_node node) } } -//============================================================================== -// RestrictedSource functionality -//============================================================================== -void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) +void Source::check_for_restriction_nodes(pugi::xml_node node) { // Check for domains to reject from if (check_for_node(node, "domain_type")) { @@ -142,23 +139,23 @@ void RestrictedSource::check_for_restriction_nodes(pugi::xml_node node) } } -bool RestrictedSource::inside_bounds(SourceSite& s) const +bool Source::inside_bounds(SourceSite& s) const { return inside_spatial_bounds(s) && inside_energy_bounds(s.E) && inside_time_bounds(s.time); } -bool RestrictedSource::inside_energy_bounds(double E) const +bool Source::inside_energy_bounds(double E) const { return E > energy_bounds_.first && E < energy_bounds_.second; } -bool RestrictedSource::inside_time_bounds(const double time) const +bool Source::inside_time_bounds(const double time) const { return time > time_bounds_.first && time < time_bounds_.second; } -bool RestrictedSource::inside_spatial_bounds(SourceSite& s) const +bool Source::inside_spatial_bounds(SourceSite& s) const { bool found = false; Particle p; From 3453353cde5eccbb395e3030bdcd27847dc87e52 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 11 Apr 2024 15:11:31 +0200 Subject: [PATCH 31/57] remove useless code --- src/source.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index a5476ec4f32..ec09a10ef68 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -419,8 +419,6 @@ SourceSite FileSource::sample(uint64_t* seed) const // Sample a particle randomly from list i_site = sites_.size() * prn(seed); site = sites_[i_site]; - Particle p; - p.r() = site.r; found = inside_bounds(site); From aa8628575623b7e7ac6209154afacfbb63b40e15 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Thu, 11 Apr 2024 15:11:56 +0200 Subject: [PATCH 32/57] don't search for xs for already dead particles --- src/simulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index a7aea69424d..2771a751ac5 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -732,7 +732,7 @@ void free_memory_simulation() void transport_history_based_single_particle(Particle& p) { - while (true) { + while (p.alive()) { p.event_calculate_xs(); if (!p.alive()) break; From ca03b4d872aa6e99966c1f1131601769ca004b1c Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Tue, 16 Apr 2024 17:25:01 +0200 Subject: [PATCH 33/57] fix particle search logic to allow fissionable to work --- src/source.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index ec09a10ef68..83956825c08 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -281,19 +281,20 @@ SourceSite IndependentSource::sample(uint64_t* seed) const static int n_accept = 0; while (!found) { - // Set particle type - Particle p; - p.type() = particle_; - p.u() = {0.0, 0.0, 1.0}; // Sample spatial distribution - p.r() = space_->sample(seed); + site.r = space_->sample(seed); // Check if otherwise outside defined restriction bounds found = inside_spatial_bounds(site); // Check if spatial site is in fissionable material if (found) { + // Create a temporary particle and search for material at the site + Particle p; + p.r() = site.r; + exhaustive_find_cell(p); + auto space_box = dynamic_cast(space_.get()); if (space_box) { if (space_box->only_fissionable()) { @@ -318,8 +319,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const "definition."); } } - - site.r = p.r(); } // Sample angle From 306137d2bb896931c5def117ab4150cf733482c6 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Tue, 16 Apr 2024 17:27:05 +0200 Subject: [PATCH 34/57] revert to having resample as the default resample has been the normal behavior asn should remain so --- src/source.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 83956825c08..f8769be666c 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -134,8 +134,8 @@ void Source::check_for_restriction_nodes(pugi::xml_node node) "Unrecognized strategy source rejection: " + rejection_strategy)); } } else { - // Default to kill rejected particles - rejection_strategy_ = RejectionStrategy::KILL; + // Default to resample rejected particles + rejection_strategy_ = RejectionStrategy::RESAMPLE; } } From 45b2a0845c7d450ef407dc6f85718ba5ed52c516 Mon Sep 17 00:00:00 2001 From: Erik B Knudsen Date: Tue, 16 Apr 2024 19:09:55 +0200 Subject: [PATCH 35/57] fix domain search logic which accepted particles regardless --- src/source.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index f8769be666c..32895504e71 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -190,9 +190,8 @@ bool Source::inside_spatial_bounds(SourceSite& s) const } } } - - // particle is inside all spatial bounds - return true; + // If at this point found is false, part. is outside all cells/materials + return found; } //============================================================================== From 5b36f803fd66d5f4889bfd668e83e39fcec6e76b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 30 Apr 2024 11:30:55 -0500 Subject: [PATCH 36/57] Check for restriction nodes in MeshSource --- include/openmc/source.h | 16 ++++++++-------- src/source.cpp | 5 ++++- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index b6b359f99c6..2dfe96ec4bd 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -59,7 +59,14 @@ class Source { enum class DomainType { UNIVERSE, MATERIAL, CELL }; enum class RejectionStrategy { KILL, RESAMPLE }; - // members + // Methods + void check_for_restriction_nodes(pugi::xml_node node); + bool inside_bounds(SourceSite& s) const; + bool inside_spatial_bounds(SourceSite& s) const; + bool inside_energy_bounds(const double E) const; + bool inside_time_bounds(const double time) const; + + // Data members std::unordered_set domain_ids_; //!< Domains to reject from DomainType domain_type_; //!< Domain type for rejection std::pair time_bounds_ {-std::numeric_limits::max(), @@ -69,13 +76,6 @@ class Source { vector lower_left_; //!< Lower left corner cds of filter vector upper_right_; //!< Upper right corner cds of filter RejectionStrategy rejection_strategy_; //!< Procedure for rejected - - // methods - void check_for_restriction_nodes(pugi::xml_node node); - bool inside_bounds(SourceSite& s) const; - bool inside_spatial_bounds(SourceSite& s) const; - bool inside_energy_bounds(const double E) const; - bool inside_time_bounds(const double time) const; }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index ddba3449dcb..bc9471caa1f 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -264,7 +264,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) time_ = UPtrDist {new Discrete {T, p, 1}}; } - // Check for additional defined restrictions (from RestrictedSource) + // Check for additional defined restrictions check_for_restriction_nodes(node); } } @@ -524,6 +524,9 @@ MeshSource::MeshSource(pugi::xml_node node) } space_ = std::make_unique(mesh_idx, strengths); + + // Check for additional defined restrictions + check_for_restriction_nodes(node); } SourceSite MeshSource::sample(uint64_t* seed) const From 96e89f847eb4111f5e8ab0173a8d4c490f028093 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 30 Apr 2024 11:38:54 -0500 Subject: [PATCH 37/57] Get rid of lower_left, upper_right. Use GeomState instead of Particle --- include/openmc/source.h | 4 +--- openmc/source.py | 50 ++++++++--------------------------------- src/source.cpp | 35 ++++++----------------------- 3 files changed, 17 insertions(+), 72 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 2dfe96ec4bd..9e4d23675de 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -73,9 +73,7 @@ class Source { std::numeric_limits::max()}; //!< time limits std::pair energy_bounds_ { 0, std::numeric_limits::max()}; //!< energy limits - vector lower_left_; //!< Lower left corner cds of filter - vector upper_right_; //!< Upper right corner cds of filter - RejectionStrategy rejection_strategy_; //!< Procedure for rejected + RejectionStrategy rejection_strategy_; //!< Procedure for rejected }; //============================================================================== diff --git a/openmc/source.py b/openmc/source.py index e24ca52d572..18d78de8c6a 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -668,10 +668,6 @@ class FileSource(SourceBase): lower, upper bounds in time [s] of accepted particles energy_bounds : Iterable of float lower, upper bounds in energy [eV] of accepted particles - lower_left : Iterable of float - Lower-left corner coordinates from which to accept particles - upper_right : Iterable of float - Upper-right corner coordinates from which to accept particles Attributes ---------- @@ -691,23 +687,25 @@ class FileSource(SourceBase): lower, upper bounds in time of accepted particles energy_bounds : Iterable of float lower, upper bounds in energy of accepted particles - lower_left : Iterable of double - Lower-left corner coordinates - upper_right : Iterable of double - Upper-right corner coordinates rejection_strategy: {'resample','kill'} Strategy when a particle is rejected. Pick a new particle or accept and terminate. """ - def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_strategy = None, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, lower_left: Optional[Sequence[Double]] = None, upper_right: Optional[Sequence[Double]] = None, time_bounds: Optional[Sequence[Double]] = None, energy_bounds: Optional[Sequence[Double]] = None) -> None: + def __init__( + self, + path: Optional[PathLike] = None, + strength=1.0, + rejection_strategy = None, + domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, + time_bounds: Optional[Sequence[float]] = None, + energy_bounds: Optional[Sequence[float]] = None + ) -> None: super().__init__(strength=strength) self._path = None self._time_bounds = None self._energy_bounds = None - self._lower_left = None - self._upper_right = None self._rejection_strategy = None if path is not None: @@ -727,10 +725,6 @@ def __init__(self, path: Optional[PathLike] = None, strength=1.0, rejection_stra self.time_bounds = time_bounds if energy_bounds is not None: self.energy_bounds = energy_bounds - if lower_left is not None: - self.lower_left = lower_left - if upper_right is not None: - self.upper_right = upper_right if rejection_strategy is not None: self.rejection_strategy = rejection_strategy @@ -785,26 +779,6 @@ def energy_bounds(self, energy_bounds): cv.check_type(name, energy_bounds, Iterable, Real) self._energy_bounds = energy_bounds - @property - def lower_left(self): - return self._lower_left - - @lower_left.setter - def lower_left(self, lower_left): - name = 'lower-left coordinates' - cv.check_type(name, lower_left, Iterable, Real) - self._lower_left = lower_left - - @property - def upper_right(self): - return self._upper_right - - @upper_right.setter - def upper_right(self, upper_right): - name = 'upper-right coordinates' - cv.check_type(name, upper_right, Iterable, Real) - self._upper_right = upper_right - @property def rejection_strategy(self): return self._rejection_strategy @@ -837,12 +811,6 @@ def populate_xml_element(self, element): if self.energy_bounds is not None: dt_elem = ET.SubElement(element, "energy_bounds") dt_elem.text = ' '.join(str(E) for E in self.energy_bounds) - if self.lower_left is not None: - dt_elem = ET.SubElement(element, "lower_left") - dt_elem.text = ' '.join(str(cd) for cd in self.lower_left) - if self.upper_right is not None: - dt_elem = ET.SubElement(element, "upper_right") - dt_elem.text = ' '.join(str(cd) for cd in self.upper_right) if self.rejection_strategy is not None: dt_elem = ET.SubElement(element, "rejection_strategy") dt_elem.text = self.rejection_strategy diff --git a/src/source.cpp b/src/source.cpp index bc9471caa1f..623286e923c 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -109,18 +109,6 @@ void Source::check_for_restriction_nodes(pugi::xml_node node) if (ids.size() == 2) energy_bounds_ = std::make_pair(ids[0], ids[1]); } - if (check_for_node(node, "lower_left")) { - auto ids = get_node_array(node, "lower_left"); - for (auto id : ids) { - lower_left_.push_back(id); - } - } - if (check_for_node(node, "upper_right")) { - auto ids = get_node_array(node, "upper_right"); - for (auto id : ids) { - upper_right_.push_back(id); - } - } // Check for how to handle rejected particles if (check_for_node(node, "rejection_strategy")) { @@ -158,33 +146,24 @@ bool Source::inside_time_bounds(const double time) const bool Source::inside_spatial_bounds(SourceSite& s) const { bool found = false; - Particle p; - p.r() = s.r; + GeometryState geom_state; + geom_state.r() = s.r; // Reject particle if it's not in the geometry at all - if (!(found = exhaustive_find_cell(p))) - return false; - - if (!lower_left_.empty() && - (s.r[0] < lower_left_[0] || s.r[1] < lower_left_[1] || - s.r[2] < lower_left_[2])) - return false; - if (!upper_right_.empty() && - (s.r[0] > upper_right_[0] || s.r[1] > upper_right_[1] || - s.r[2] > upper_right_[2])) + if (!(found = exhaustive_find_cell(geom_state))) return false; if (!domain_ids_.empty()) { if (domain_type_ == DomainType::MATERIAL) { - auto mat_index = p.material(); + auto mat_index = geom_state.material(); if (mat_index != MATERIAL_VOID) { found = contains(domain_ids_, model::materials[mat_index]->id()); } } else { - for (int i = 0; i < p.n_coord(); i++) { + for (int i = 0; i < geom_state.n_coord(); i++) { auto id = (domain_type_ == DomainType::CELL) - ? model::cells[p.coord(i).cell]->id_ - : model::universes[p.coord(i).universe]->id_; + ? model::cells[geom_state.coord(i).cell]->id_ + : model::universes[geom_state.coord(i).universe]->id_; if ((found = contains(domain_ids_, id))) break; } From 44b4c20a7e7f478e80c1003bdee3d4f83aa91e96 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 30 Apr 2024 23:20:11 -0500 Subject: [PATCH 38/57] Consolidate common logic in source.py --- openmc/source.py | 453 +++++++++++++++++++++++++++-------------------- 1 file changed, 263 insertions(+), 190 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 18d78de8c6a..26717bf7cbf 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -28,6 +28,16 @@ class SourceBase(ABC): ---------- strength : float Strength of the source + 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. + time_bounds : sequence of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : sequence of float + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). Attributes ---------- @@ -35,11 +45,50 @@ class SourceBase(ABC): Indicator of source type. strength : float Strength of the source + domain_type : {'cell', 'material', 'universe'} + Type of domain to use for rejection + domain_ids : iterable of int + Ids of domains to reject based on. + time_bounds : Iterable of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : Iterable of float + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). """ - def __init__(self, strength=1.0): + def __init__( + self, + strength: Optional[float] = 1.0, + domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, + time_bounds: Optional[Sequence[float]] = None, + energy_bounds: Optional[Sequence[float]] = None, + rejection_strategy: Optional[str] = None, + ): self.strength = strength + self._time_bounds = None + self._energy_bounds = None + self._rejection_strategy = None + 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] + + if time_bounds is not None: + self.time_bounds = time_bounds + if energy_bounds is not None: + self.energy_bounds = energy_bounds + if rejection_strategy is not None: + self.rejection_strategy = rejection_strategy @property def strength(self): @@ -47,10 +96,58 @@ def strength(self): @strength.setter def strength(self, strength): - cv.check_type('source strength', strength, Real) - cv.check_greater_than('source strength', strength, 0.0, True) + cv.check_type('source strength', strength, Real, none_ok=True) + if strength is not None: + cv.check_greater_than('source strength', strength, 0.0, True) self._strength = strength + @property + def domain_ids(self): + return self._domain_ids + + @domain_ids.setter + def domain_ids(self, ids): + cv.check_type('domain IDs', ids, Iterable, Real) + self._domain_ids = ids + + @property + def domain_type(self): + return self._domain_type + + @domain_type.setter + def domain_type(self, domain_type): + cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) + self._domain_type = domain_type + + @property + def time_bounds(self): + return self._time_bounds + + @time_bounds.setter + def time_bounds(self, time_bounds): + name = 'time bounds' + cv.check_type(name, time_bounds, Iterable, Real) + self._time_bounds = tuple(time_bounds) + + @property + def energy_bounds(self): + return self._energy_bounds + + @energy_bounds.setter + def energy_bounds(self, energy_bounds): + name = 'energy bounds' + cv.check_type(name, energy_bounds, Iterable, Real) + self._energy_bounds = tuple(energy_bounds) + + @property + def rejection_strategy(self): + return self._rejection_strategy + + @rejection_strategy.setter + def rejection_strategy(self, rejection_strategy): + cv.check_value('rejection strategy', rejection_strategy, ('resample', 'kill')) + self._rejection_strategy = rejection_strategy + @abstractmethod def populate_xml_element(self, element): """Add necessary source information to an XML element @@ -73,8 +170,24 @@ def to_xml_element(self) -> ET.Element: """ element = ET.Element("source") element.set("type", self.type) - element.set("strength", str(self.strength)) + if self.strength is not None: + element.set("strength", str(self.strength)) self.populate_xml_element(element) + 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) + if self.time_bounds is not None: + dt_elem = ET.SubElement(element, "time_bounds") + dt_elem.text = ' '.join(str(t) for t in self.time_bounds) + if self.energy_bounds is not None: + dt_elem = ET.SubElement(element, "energy_bounds") + dt_elem.text = ' '.join(str(E) for E in self.energy_bounds) + if self.rejection_strategy is not None: + dt_elem = ET.SubElement(element, "rejection_strategy") + dt_elem.text = self.rejection_strategy + return element @classmethod @@ -118,6 +231,42 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: else: raise ValueError(f'Source type {source_type} is not recognized') + @staticmethod + def _get_common_arguments(elem: ET.Element) -> dict: + 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 + + time_bounds = get_text(elem, "time_bounds") + if time_bounds is not None: + time_bounds = [float(x) for x in time_bounds.split()] + + energy_bounds = get_text(elem, "energy_bounds") + if energy_bounds is not None: + energy_bounds = [float(x) for x in energy_bounds.split()] + + rejection_strategy = get_text(elem, "rejection_strategy") + + return { + 'domains': domains, + 'time_bounds': time_bounds, + 'energy_bounds': energy_bounds, + 'rejection_strategy': rejection_strategy + } + class IndependentSource(SourceBase): """Distribution of phase space coordinates for source sites. @@ -141,6 +290,13 @@ class IndependentSource(SourceBase): 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. + time_bounds : sequence of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : sequence of float + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). Attributes ---------- @@ -176,9 +332,15 @@ def __init__( time: Optional[openmc.stats.Univariate] = None, strength: float = 1.0, particle: str = 'neutron', - domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None + domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, + time_bounds: Optional[Sequence[float]] = None, + energy_bounds: Optional[Sequence[float]] = None, + rejection_strategy: Optional[str] = None, ): - super().__init__(strength) + super().__init__( + strength=strength, domains=domains, time_bounds=time_bounds, + energy_bounds=energy_bounds, rejection_strategy=rejection_strategy + ) self._space = None self._angle = None @@ -193,20 +355,8 @@ def __init__( self.energy = energy if time is not None: self.time = time - self.strength = strength self.particle = particle - 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 type(self) -> str: return 'independent' @@ -273,24 +423,6 @@ def particle(self, particle): cv.check_value('source particle', particle, ['neutron', 'photon']) self._particle = particle - @property - def domain_ids(self): - return self._domain_ids - - @domain_ids.setter - def domain_ids(self, ids): - cv.check_type('domain IDs', ids, Iterable, Real) - self._domain_ids = ids - - @property - def domain_type(self): - return self._domain_type - - @domain_type.setter - def domain_type(self, domain_type): - cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) - self._domain_type = domain_type - def populate_xml_element(self, element): """Add necessary source information to an XML element @@ -309,11 +441,6 @@ def populate_xml_element(self, element): element.append(self.energy.to_xml_element('energy')) if self.time is not None: element.append(self.time.to_xml_element('time')) - 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) @classmethod def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: @@ -333,24 +460,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: Source generated from XML element """ - 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) + kwargs = cls._get_common_arguments(elem) + source = cls(**kwargs) strength = get_text(elem, 'strength') if strength is not None: @@ -360,10 +471,6 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: if particle is not None: source.particle = particle - filename = get_text(elem, 'file') - if filename is not None: - source.file = filename - space = elem.find('space') if space is not None: source.space = Spatial.from_xml_element(space, meshes) @@ -403,6 +510,16 @@ class MeshSource(SourceBase): Sources for each element in the mesh. If spatial distributions are set on any of the source objects, they will be ignored during source site sampling. + 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. + time_bounds : sequence of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : sequence of float + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). Attributes ---------- @@ -417,9 +534,32 @@ class MeshSource(SourceBase): Strength of the source type : str Indicator of source type: 'mesh' + domain_type : {'cell', 'material', 'universe'} + Type of domain to use for rejection + domain_ids : iterable of int + Ids of domains to reject based on. + time_bounds : Iterable of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : Iterable of float + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). """ - def __init__(self, mesh: MeshBase, sources: Sequence[SourceBase]): + def __init__( + self, + mesh: MeshBase, + sources: Sequence[SourceBase], + domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, + time_bounds: Optional[Sequence[float]] = None, + energy_bounds: Optional[Sequence[float]] = None, + rejection_strategy: Optional[str] = None, + ): + super().__init__( + strength=None, domains=domains, time_bounds=time_bounds, + energy_bounds=energy_bounds, rejection_strategy=rejection_strategy + ) self.mesh = mesh self.sources = sources @@ -467,8 +607,9 @@ def sources(self, s): @strength.setter def strength(self, val): - cv.check_type('mesh source strength', val, Real) - self.set_total_strength(val) + if val is not None: + cv.check_type('mesh source strength', val, Real) + self.set_total_strength(val) def set_total_strength(self, strength: float): """Scales the element source strengths based on a desired total strength. @@ -532,7 +673,8 @@ def from_xml_element(cls, elem: ET.Element, meshes) -> openmc.MeshSource: sources = [SourceBase.from_xml_element(e) for e in elem.iterchildren('source')] sources = np.asarray(sources).reshape(mesh.dimension, order='F') - return cls(mesh, sources) + kwargs = cls._get_common_arguments(elem) + return cls(mesh, sources, **kwargs) def Source(*args, **kwargs): @@ -557,6 +699,16 @@ class CompiledSource(SourceBase): Parameters to be provided to the compiled shared library function strength : float Strength of the source + 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. + time_bounds : sequence of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : sequence of float + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). Attributes ---------- @@ -568,10 +720,33 @@ class CompiledSource(SourceBase): Strength of the source type : str Indicator of source type: 'compiled' + domain_type : {'cell', 'material', 'universe'} + Type of domain to use for rejection + domain_ids : iterable of int + Ids of domains to reject based on. + time_bounds : Iterable of float + lower, upper bounds in time [s] of accepted particles + energy_bounds : Iterable of float + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). """ - def __init__(self, library: Optional[str] = None, parameters: Optional[str] = None, strength=1.0) -> None: - super().__init__(strength=strength) + def __init__( + self, + library: Optional[str] = None, + parameters: Optional[str] = None, + strength: float = 1.0, + domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, + time_bounds: Optional[Sequence[float]] = None, + energy_bounds: Optional[Sequence[float]] = None, + rejection_strategy: Optional[str] = None, + ) -> None: + super().__init__( + strength=strength, domains=domains, time_bounds=time_bounds, + energy_bounds=energy_bounds, rejection_strategy=rejection_strategy + ) self._library = None if library is not None: @@ -635,9 +810,10 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.CompiledSource: Source generated from XML element """ - library = get_text(elem, 'library') + kwargs = cls._get_common_arguments(elem) + kwargs['library'] = get_text(elem, 'library') - source = cls(library) + source = cls(**kwargs) strength = get_text(elem, 'strength') if strength is not None: @@ -684,50 +860,31 @@ class FileSource(SourceBase): domain_ids : iterable of int Ids of domains to reject based on. time_bounds : Iterable of float - lower, upper bounds in time of accepted particles + lower, upper bounds in time [s] of accepted particles energy_bounds : Iterable of float - lower, upper bounds in energy of accepted particles - rejection_strategy: {'resample','kill'} - Strategy when a particle is rejected. Pick a new particle or accept and terminate. - + lower, upper bounds in energy [eV] of accepted particles + rejection_strategy : {'resample', 'kill'} + Strategy when a particle is rejected. Pick a new particle ('resample') + or accept and terminate ('kill'). """ def __init__( self, path: Optional[PathLike] = None, - strength=1.0, - rejection_strategy = None, + strength: float = 1.0, domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, time_bounds: Optional[Sequence[float]] = None, - energy_bounds: Optional[Sequence[float]] = None + energy_bounds: Optional[Sequence[float]] = None, + rejection_strategy = None, ) -> None: - super().__init__(strength=strength) - + super().__init__( + strength=strength, domains=domains, time_bounds=time_bounds, + energy_bounds=energy_bounds, rejection_strategy=rejection_strategy + ) self._path = None - self._time_bounds = None - self._energy_bounds = None - self._rejection_strategy = None - if path is not None: self.path = path - 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] - if time_bounds is not None: - self.time_bounds = time_bounds - if energy_bounds is not None: - self.energy_bounds = energy_bounds - if rejection_strategy is not None: - self.rejection_strategy = rejection_strategy - @property def type(self) -> str: return "file" @@ -741,54 +898,6 @@ def path(self, p: PathLike): cv.check_type('source file', p, str) self._path = p - @property - def domain_ids(self): - return self._domain_ids - - @domain_ids.setter - def domain_ids(self, ids): - cv.check_type('domain IDs', ids, Iterable, Real) - self._domain_ids = ids - - @property - def domain_type(self): - return self._domain_type - - @domain_type.setter - def domain_type(self, domain_type): - cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) - self._domain_type = domain_type - - @property - def time_bounds(self): - return self._time_bounds - - @time_bounds.setter - def time_bounds(self, time_bounds): - name = 'time bounds' - cv.check_type(name, time_bounds, Iterable, Real) - self._time_bounds = time_bounds - - @property - def energy_bounds(self): - return self._energy_bounds - - @energy_bounds.setter - def energy_bounds(self, energy_bounds): - name = 'energy bounds' - cv.check_type(name, energy_bounds, Iterable, Real) - self._energy_bounds = energy_bounds - - @property - def rejection_strategy(self): - return self._rejection_strategy - - @rejection_strategy.setter - def rejection_strategy(self, rejection_strategy): - name = 'rejection strategy' - cv.check_value('rejection strategy', rejection_strategy, ('resample','kill')) - self._rejection_strategy = rejection_strategy - def populate_xml_element(self, element): """Add necessary file source information to an XML element @@ -800,20 +909,6 @@ def populate_xml_element(self, element): """ if self.path is not None: element.set("file", self.path) - 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) - if self.time_bounds is not None: - dt_elem = ET.SubElement(element, "time_bounds") - dt_elem.text = ' '.join(str(t) for t in self.time_bounds) - if self.energy_bounds is not None: - dt_elem = ET.SubElement(element, "energy_bounds") - dt_elem.text = ' '.join(str(E) for E in self.energy_bounds) - if self.rejection_strategy is not None: - dt_elem = ET.SubElement(element, "rejection_strategy") - dt_elem.text = self.rejection_strategy @classmethod def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: @@ -833,35 +928,13 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: Source generated from XML element """ - - filename = get_text(elem, 'file') - - source = cls(filename) - + kwargs = cls._get_common_arguments(elem) + kwargs['path'] = get_text(elem, 'file') strength = get_text(elem, 'strength') if strength is not None: - source.strength = float(strength) - - 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 + kwargs['strength'] = float(strength) - source = cls(domains=domains) - - return source + return cls(**kwargs) class ParticleType(IntEnum): From cf6049927539c784f7c30089608e40a1a4299d0e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 1 May 2024 09:04:23 -0500 Subject: [PATCH 39/57] Generalize source rejection on C++ side --- include/openmc/source.h | 41 +++++++++++++++++++------------ src/source.cpp | 54 ++++++++++++++++++++++------------------- 2 files changed, 55 insertions(+), 40 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 9e4d23675de..a9160f9d2be 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -46,10 +46,8 @@ class Source { public: virtual ~Source() = default; - // Methods that must be implemented - virtual SourceSite sample(uint64_t* seed) const = 0; - // Methods that can be overridden + virtual SourceSite sample(uint64_t* seed) const; virtual double strength() const { return 1.0; } static unique_ptr create(pugi::xml_node node); @@ -59,10 +57,13 @@ class Source { enum class DomainType { UNIVERSE, MATERIAL, CELL }; enum class RejectionStrategy { KILL, RESAMPLE }; + // Methods that can be overridden + virtual SourceSite sample_no_rejection(uint64_t* seed) const = 0; + // Methods void check_for_restriction_nodes(pugi::xml_node node); bool inside_bounds(SourceSite& s) const; - bool inside_spatial_bounds(SourceSite& s) const; + bool inside_spatial_bounds(Position r) const; bool inside_energy_bounds(const double E) const; bool inside_time_bounds(const double time) const; @@ -103,6 +104,11 @@ class IndependentSource : public Source { Distribution* energy() const { return energy_.get(); } Distribution* time() const { return time_.get(); } +protected: + // Note that this method is not used since IndependentSource directly + // overrides the sample method + SourceSite sample_no_rejection(uint64_t* seed) const override { return {}; } + private: // Data members ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted @@ -124,9 +130,12 @@ class FileSource : public Source { explicit FileSource(const std::string& path); // Methods - SourceSite sample(uint64_t* seed) const override; void load_sites_from_file( const std::string& path); //!< Load source sites from file + +protected: + SourceSite sample_no_rejection(uint64_t* seed) const override; + private: vector sites_; //!< Source sites from a file }; @@ -141,16 +150,17 @@ class CompiledSourceWrapper : public Source { CompiledSourceWrapper(pugi::xml_node node); ~CompiledSourceWrapper(); + double strength() const override { return compiled_source_->strength(); } + + void setup(const std::string& path, const std::string& parameters); + +protected: // Defer implementation to custom source library - SourceSite sample(uint64_t* seed) const override + SourceSite sample_no_rejection(uint64_t* seed) const override { return compiled_source_->sample(seed); } - double strength() const override { return compiled_source_->strength(); } - - void setup(const std::string& path, const std::string& parameters); - private: void* shared_library_; //!< library from dlopen unique_ptr compiled_source_; @@ -167,11 +177,6 @@ class MeshSource : public Source { // Constructors explicit MeshSource(pugi::xml_node node); - //! Sample from the external source distribution - //! \param[inout] seed Pseudorandom seed pointer - //! \return Sampled site - SourceSite sample(uint64_t* seed) const override; - // Properties double strength() const override { return space_->total_strength(); } @@ -181,6 +186,12 @@ class MeshSource : public Source { return sources_.size() == 1 ? sources_[0] : sources_[i]; } +protected: + //! Sample from the external source distribution + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + SourceSite sample_no_rejection(uint64_t* seed) const override; + private: // Data members unique_ptr space_; //!< Mesh spatial diff --git a/src/source.cpp b/src/source.cpp index 623286e923c..1fcef600b48 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -127,9 +127,29 @@ void Source::check_for_restriction_nodes(pugi::xml_node node) } } +SourceSite Source::sample(uint64_t* seed) const +{ + bool found = false; + SourceSite site; + + while (!found) { + // Sample a particle randomly from list + site = this->sample_no_rejection(seed); + + found = inside_bounds(site); + if (!found && rejection_strategy_ == RejectionStrategy::KILL) { + // Accept particle regardless but set wgt=0 to trigger garbage collection + // I.e. kill this particle immediately and pick the next + found = true; + site.wgt = 0.0; + } + } + return site; +} + bool Source::inside_bounds(SourceSite& s) const { - return inside_spatial_bounds(s) && inside_energy_bounds(s.E) && + return inside_spatial_bounds(s.r) && inside_energy_bounds(s.E) && inside_time_bounds(s.time); } @@ -143,11 +163,11 @@ bool Source::inside_time_bounds(const double time) const return time > time_bounds_.first && time < time_bounds_.second; } -bool Source::inside_spatial_bounds(SourceSite& s) const +bool Source::inside_spatial_bounds(Position r) const { bool found = false; GeometryState geom_state; - geom_state.r() = s.r; + geom_state.r() = r; // Reject particle if it's not in the geometry at all if (!(found = exhaustive_find_cell(geom_state))) @@ -264,7 +284,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const site.r = space_->sample(seed); // Check if otherwise outside defined restriction bounds - found = inside_spatial_bounds(site); + found = inside_spatial_bounds(site.r); // Check if spatial site is in fissionable material if (found) { @@ -390,27 +410,11 @@ void FileSource::load_sites_from_file(const std::string& path) file_close(file_id); } -SourceSite FileSource::sample(uint64_t* seed) const +SourceSite FileSource::sample_no_rejection(uint64_t* seed) const { - bool found = false; - size_t i_site; - SourceSite site; - - while (!found) { - // Sample a particle randomly from list - i_site = sites_.size() * prn(seed); - site = sites_[i_site]; - - found = inside_bounds(site); - - if (!found && rejection_strategy_ == RejectionStrategy::KILL) { - // Accept particle regardless but set wgt=0 to trigger garbage collection - // I.e. kill this particle immediately and pick the next - found = true; - site.wgt = 0.0; - } - } - return site; + // Sample a particle randomly from list + size_t i_site = sites_.size() * prn(seed); + return sites_[i_site]; } //============================================================================== @@ -508,7 +512,7 @@ MeshSource::MeshSource(pugi::xml_node node) check_for_restriction_nodes(node); } -SourceSite MeshSource::sample(uint64_t* seed) const +SourceSite MeshSource::sample_no_rejection(uint64_t* seed) const { // sample location and element from mesh auto mesh_location = space_->sample_mesh(seed); From 03951509d172d3c945237b7cb17f2206e0dc951a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 1 May 2024 09:48:14 -0500 Subject: [PATCH 40/57] Preserve per-element probabilities for MeshSource --- include/openmc/distribution_spatial.h | 5 +++++ include/openmc/source.h | 11 +++++---- src/distribution_spatial.cpp | 7 +++++- src/source.cpp | 32 +++++++++++++++++++-------- 4 files changed, 39 insertions(+), 16 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index bd10a299605..28568c890d8 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -116,6 +116,11 @@ class MeshSpatial : public SpatialDistribution { //! \return Sampled element index and position within that element std::pair sample_mesh(uint64_t* seed) const; + //! Sample a mesh element + //! \param seed Pseudorandom number seed pointer + //! \return Sampled element index + int32_t sample_element_index(uint64_t* seed) const; + //! For unstructured meshes, ensure that elements are all linear tetrahedra void check_element_types() const; diff --git a/include/openmc/source.h b/include/openmc/source.h index a9160f9d2be..6b3a867d048 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -177,6 +177,11 @@ class MeshSource : public Source { // Constructors explicit MeshSource(pugi::xml_node node); + //! Sample from the external source distribution + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + SourceSite sample(uint64_t* seed) const override; + // Properties double strength() const override { return space_->total_strength(); } @@ -186,12 +191,6 @@ class MeshSource : public Source { return sources_.size() == 1 ? sources_[0] : sources_[i]; } -protected: - //! Sample from the external source distribution - //! \param[inout] seed Pseudorandom seed pointer - //! \return Sampled site - SourceSite sample_no_rejection(uint64_t* seed) const override; - private: // Data members unique_ptr space_; //!< Mesh spatial diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 8f34ea6b936..8dbd7d88706 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -281,10 +281,15 @@ void MeshSpatial::check_element_types() const } } +int32_t MeshSpatial::sample_element_index(uint64_t* seed) const +{ + return elem_idx_dist_.sample(seed); +} + std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { // Sample the CDF defined in initialization above - int32_t elem_idx = elem_idx_dist_.sample(seed); + int32_t elem_idx = this->sample_element_index(seed); return {elem_idx, mesh()->sample_element(elem_idx, seed)}; } diff --git a/src/source.cpp b/src/source.cpp index 1fcef600b48..1f4a2b33d0b 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -512,17 +512,31 @@ MeshSource::MeshSource(pugi::xml_node node) check_for_restriction_nodes(node); } -SourceSite MeshSource::sample_no_rejection(uint64_t* seed) const +SourceSite MeshSource::sample(uint64_t* seed) const { - // sample location and element from mesh - auto mesh_location = space_->sample_mesh(seed); - - // Sample source for the chosen element - int32_t element = mesh_location.first; - SourceSite site = source(element)->sample(seed); + // Sample the CDF defined in initialization above + int32_t element = space_->sample_element_index(seed); + + // Sample position and apply rejection on spatial domains + Position r; + while (true) { + r = space_->mesh()->sample_element(element, seed); + if (this->inside_spatial_bounds(r)) { + break; + } + } - // Replace spatial position with the one already sampled - site.r = mesh_location.second; + SourceSite site; + while (true) { + // Sample source for the chosen element and replace the position + site = source(element)->sample(seed); + site.r = r; + + // Apply other rejections + if (inside_energy_bounds(site.E) && inside_time_bounds(site.time)) { + break; + } + } return site; } From 71538d23e086f909d515c96f501968e1d179939c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 1 May 2024 10:10:01 -0500 Subject: [PATCH 41/57] Add missing dummy method on MeshSource --- include/openmc/source.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/openmc/source.h b/include/openmc/source.h index 6b3a867d048..a8869004423 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -191,6 +191,11 @@ class MeshSource : public Source { return sources_.size() == 1 ? sources_[0] : sources_[i]; } +protected: + // Note that this method is not used since MeshSource directly overrides the + // sample method + SourceSite sample_no_rejection(uint64_t* seed) const override { return {}; } + private: // Data members unique_ptr space_; //!< Mesh spatial From 4ca3996b26334e208e33a68953159d287636d2a9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 May 2024 09:12:04 -0500 Subject: [PATCH 42/57] Remove superfluous check in transport_history_based_single_particle --- src/simulation.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index 9827bb38282..527d98db4f1 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -761,8 +761,6 @@ void transport_history_based_single_particle(Particle& p) p.event_collide(); } p.event_revive_from_secondary(); - if (!p.alive()) - break; } p.event_death(); } From 9755ac93c72eb5beabfa5b62d2bad1cebd8320c7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 May 2024 10:31:57 -0500 Subject: [PATCH 43/57] Use single constraints argument in source classes --- include/openmc/source.h | 10 +- openmc/source.py | 392 ++++++++++++++++------------------------ src/source.cpp | 37 ++-- 3 files changed, 185 insertions(+), 254 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index a8869004423..e7e0d0598aa 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -61,11 +61,11 @@ class Source { virtual SourceSite sample_no_rejection(uint64_t* seed) const = 0; // Methods - void check_for_restriction_nodes(pugi::xml_node node); - bool inside_bounds(SourceSite& s) const; - bool inside_spatial_bounds(Position r) const; - bool inside_energy_bounds(const double E) const; - bool inside_time_bounds(const double time) const; + void check_for_constraints(pugi::xml_node node); + bool satisfies_constraints(SourceSite& s) const; + bool satisfies_spatial_constraints(Position r) const; + bool satisfies_energy_constraints(const double E) const; + bool satisfies_time_constraints(const double time) const; // Data members std::unordered_set domain_ids_; //!< Domains to reject from diff --git a/openmc/source.py b/openmc/source.py index 3a67b9c4e39..90d5ff7fe6e 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -6,7 +6,7 @@ import warnings import typing # imported separately as py3.8 requires typing.Iterable # also required to prevent typing.Union namespace overwriting Union -from typing import Optional, Sequence +from typing import Optional, Sequence, Dict, Any import lxml.etree as ET import numpy as np @@ -28,16 +28,17 @@ class SourceBase(ABC): ---------- strength : float Strength of the source - 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. - time_bounds : sequence of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : sequence of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', + the corresponding value is an iterable of :class:`openmc.Cell`, + :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled + sites must be within. For 'time_bounds' and 'energy_bounds', the + corresponding value is a sequence of floats giving the lower and upper + bounds on time in [s] or energy in [eV] that the sampled particle must + be within. The 'rejection_strategy' indicates what should happen when a + source particle is rejected: either 'resample' (pick a new particle) or + 'kill' (accept and terminate). Attributes ---------- @@ -45,50 +46,20 @@ class SourceBase(ABC): Indicator of source type. strength : float Strength of the source - domain_type : {'cell', 'material', 'universe'} - Type of domain to use for rejection - domain_ids : iterable of int - Ids of domains to reject based on. - time_bounds : Iterable of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : Iterable of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and + 'rejection_strategy'. """ def __init__( - self, - strength: Optional[float] = 1.0, - domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, - time_bounds: Optional[Sequence[float]] = None, - energy_bounds: Optional[Sequence[float]] = None, - rejection_strategy: Optional[str] = None, + self, + strength: Optional[float] = 1.0, + constraints: Optional[Dict[str, Any]] = None ): self.strength = strength - self._time_bounds = None - self._energy_bounds = None - self._rejection_strategy = None - 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] - - if time_bounds is not None: - self.time_bounds = time_bounds - if energy_bounds is not None: - self.energy_bounds = energy_bounds - if rejection_strategy is not None: - self.rejection_strategy = rejection_strategy + self.constraints = constraints @property def strength(self): @@ -102,51 +73,37 @@ def strength(self, strength): self._strength = strength @property - def domain_ids(self): - return self._domain_ids - - @domain_ids.setter - def domain_ids(self, ids): - cv.check_type('domain IDs', ids, Iterable, Real) - self._domain_ids = ids - - @property - def domain_type(self): - return self._domain_type - - @domain_type.setter - def domain_type(self, domain_type): - cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) - self._domain_type = domain_type - - @property - def time_bounds(self): - return self._time_bounds - - @time_bounds.setter - def time_bounds(self, time_bounds): - name = 'time bounds' - cv.check_type(name, time_bounds, Iterable, Real) - self._time_bounds = tuple(time_bounds) - - @property - def energy_bounds(self): - return self._energy_bounds - - @energy_bounds.setter - def energy_bounds(self, energy_bounds): - name = 'energy bounds' - cv.check_type(name, energy_bounds, Iterable, Real) - self._energy_bounds = tuple(energy_bounds) - - @property - def rejection_strategy(self): - return self._rejection_strategy - - @rejection_strategy.setter - def rejection_strategy(self, rejection_strategy): - cv.check_value('rejection strategy', rejection_strategy, ('resample', 'kill')) - self._rejection_strategy = rejection_strategy + def constraints(self) -> Dict[str, Any]: + return self._constraints + + @constraints.setter + def constraints(self, constraints: Optional[Dict[str, Any]]): + self._constraints = {} + if constraints is None: + return + + for key, value in constraints.items(): + if key == 'domains': + cv.check_type('domains', value, Iterable, + (openmc.Cell, openmc.Material, openmc.Universe)) + if isinstance(value[0], openmc.Cell): + self._constraints['domain_type'] = 'cell' + elif isinstance(value[0], openmc.Material): + self._constraints['domain_type'] = 'material' + elif isinstance(value[0], openmc.Universe): + self._constraints['domain_type'] = 'universe' + self._constraints['domain_ids'] = [d.id for d in value] + elif key == 'time_bounds': + cv.check_type('time bounds', value, Iterable, Real) + self._constraints['time_bounds'] = tuple(value) + elif key == 'energy_bounds': + cv.check_type('energy bounds', value, Iterable, Real) + self._constraints['energy_bounds'] = tuple(value) + elif key == 'rejection_strategy': + cv.check_value('rejection strategy', value, ('resample', 'kill')) + self._constraints['rejection_strategy'] = value + else: + raise ValueError('Unknown key in constraints dictionary: {key}') @abstractmethod def populate_xml_element(self, element): @@ -173,20 +130,21 @@ def to_xml_element(self) -> ET.Element: if self.strength is not None: element.set("strength", str(self.strength)) self.populate_xml_element(element) - if self.domain_ids: + constraints = self.constraints + if "domain_ids" in constraints: dt_elem = ET.SubElement(element, "domain_type") - dt_elem.text = self.domain_type + dt_elem.text = constraints["domain_type"] id_elem = ET.SubElement(element, "domain_ids") - id_elem.text = ' '.join(str(uid) for uid in self.domain_ids) - if self.time_bounds is not None: + id_elem.text = ' '.join(str(uid) for uid in constraints["domain_ids"]) + if "time_bounds" in constraints: dt_elem = ET.SubElement(element, "time_bounds") - dt_elem.text = ' '.join(str(t) for t in self.time_bounds) - if self.energy_bounds is not None: + dt_elem.text = ' '.join(str(t) for t in constraints["time_bounds"]) + if "energy_bounds" in constraints: dt_elem = ET.SubElement(element, "energy_bounds") - dt_elem.text = ' '.join(str(E) for E in self.energy_bounds) - if self.rejection_strategy is not None: + dt_elem.text = ' '.join(str(E) for E in constraints["energy_bounds"]) + if "rejection_strategy" in constraints: dt_elem = ET.SubElement(element, "rejection_strategy") - dt_elem.text = self.rejection_strategy + dt_elem.text = constraints["rejection_strategy"] return element @@ -232,7 +190,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: raise ValueError(f'Source type {source_type} is not recognized') @staticmethod - def _get_common_arguments(elem: ET.Element) -> dict: + def _get_constraints(elem: ET.Element) -> Dict[str, Any]: + constraints = {} 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()] @@ -247,25 +206,21 @@ def _get_common_arguments(elem: ET.Element) -> dict: 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 + constraints['domains'] = domains time_bounds = get_text(elem, "time_bounds") if time_bounds is not None: - time_bounds = [float(x) for x in time_bounds.split()] + constraints['time_bounds'] = [float(x) for x in time_bounds.split()] energy_bounds = get_text(elem, "energy_bounds") if energy_bounds is not None: - energy_bounds = [float(x) for x in energy_bounds.split()] + constraints['energy_bounds'] = [float(x) for x in energy_bounds.split()] rejection_strategy = get_text(elem, "rejection_strategy") + if rejection_strategy is not None: + constraints['rejection_strategy'] = rejection_strategy - return { - 'domains': domains, - 'time_bounds': time_bounds, - 'energy_bounds': energy_bounds, - 'rejection_strategy': rejection_strategy - } + return constraints class IndependentSource(SourceBase): @@ -290,13 +245,20 @@ class IndependentSource(SourceBase): 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. - time_bounds : sequence of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : sequence of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + + .. deprecated:: 0.14.1 + Use the `constraints` argument instead. + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', + the corresponding value is an iterable of :class:`openmc.Cell`, + :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled + sites must be within. For 'time_bounds' and 'energy_bounds', the + corresponding value is a sequence of floats giving the lower and upper + bounds on time in [s] or energy in [eV] that the sampled particle must + be within. The 'rejection_strategy' indicates what should happen when a + source particle is rejected: either 'resample' (pick a new particle) or + 'kill' (accept and terminate). Attributes ---------- @@ -317,10 +279,10 @@ class IndependentSource(SourceBase): particle : {'neutron', 'photon'} Source particle type - ids : Iterable of int - IDs of domains to use for rejection - domain_type : {'cell', 'material', 'universe'} - Type of domain to use for rejection + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and + 'rejection_strategy'. """ @@ -333,14 +295,14 @@ def __init__( strength: float = 1.0, particle: str = 'neutron', domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, - time_bounds: Optional[Sequence[float]] = None, - energy_bounds: Optional[Sequence[float]] = None, - rejection_strategy: Optional[str] = None, + constraints: Optional[Dict[str, Any]] = None ): - super().__init__( - strength=strength, domains=domains, time_bounds=time_bounds, - energy_bounds=energy_bounds, rejection_strategy=rejection_strategy - ) + if domains is not None: + warnings.warn("The 'domains' arguments has been replaced by the " + "'constraints' argument.", FutureWarning) + constraints = {'domains': domains} + + super().__init__(strength=strength, constraints=constraints) self._space = None self._angle = None @@ -460,8 +422,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: Source generated from XML element """ - kwargs = cls._get_common_arguments(elem) - source = cls(**kwargs) + constraints = cls._get_constraints(elem) + source = cls(constraints=constraints) strength = get_text(elem, 'strength') if strength is not None: @@ -512,16 +474,17 @@ class MeshSource(SourceBase): multidimensional array whose shape matches the mesh shape. If spatial distributions are set on any of the source objects, they will be ignored during source site sampling. - 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. - time_bounds : sequence of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : sequence of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', + the corresponding value is an iterable of :class:`openmc.Cell`, + :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled + sites must be within. For 'time_bounds' and 'energy_bounds', the + corresponding value is a sequence of floats giving the lower and upper + bounds on time in [s] or energy in [eV] that the sampled particle must + be within. The 'rejection_strategy' indicates what should happen when a + source particle is rejected: either 'resample' (pick a new particle) or + 'kill' (accept and terminate). Attributes ---------- @@ -533,32 +496,19 @@ class MeshSource(SourceBase): Strength of the source type : str Indicator of source type: 'mesh' - domain_type : {'cell', 'material', 'universe'} - Type of domain to use for rejection - domain_ids : iterable of int - Ids of domains to reject based on. - time_bounds : Iterable of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : Iterable of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and + 'rejection_strategy'. """ def __init__( self, mesh: MeshBase, sources: Sequence[SourceBase], - domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, - time_bounds: Optional[Sequence[float]] = None, - energy_bounds: Optional[Sequence[float]] = None, - rejection_strategy: Optional[str] = None, + constraints: Optional[Dict[str, Any]] = None, ): - super().__init__( - strength=None, domains=domains, time_bounds=time_bounds, - energy_bounds=energy_bounds, rejection_strategy=rejection_strategy - ) + super().__init__(strength=None, constraints=constraints) self.mesh = mesh self.sources = sources @@ -672,8 +622,8 @@ def from_xml_element(cls, elem: ET.Element, meshes) -> openmc.MeshSource: mesh = meshes[mesh_id] sources = [SourceBase.from_xml_element(e) for e in elem.iterchildren('source')] - kwargs = cls._get_common_arguments(elem) - return cls(mesh, sources, **kwargs) + constraints = cls._get_constraints(elem) + return cls(mesh, sources, constraints=constraints) def Source(*args, **kwargs): @@ -698,16 +648,17 @@ class CompiledSource(SourceBase): Parameters to be provided to the compiled shared library function strength : float Strength of the source - 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. - time_bounds : sequence of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : sequence of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', + the corresponding value is an iterable of :class:`openmc.Cell`, + :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled + sites must be within. For 'time_bounds' and 'energy_bounds', the + corresponding value is a sequence of floats giving the lower and upper + bounds on time in [s] or energy in [eV] that the sampled particle must + be within. The 'rejection_strategy' indicates what should happen when a + source particle is rejected: either 'resample' (pick a new particle) or + 'kill' (accept and terminate). Attributes ---------- @@ -719,33 +670,20 @@ class CompiledSource(SourceBase): Strength of the source type : str Indicator of source type: 'compiled' - domain_type : {'cell', 'material', 'universe'} - Type of domain to use for rejection - domain_ids : iterable of int - Ids of domains to reject based on. - time_bounds : Iterable of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : Iterable of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and + 'rejection_strategy'. """ def __init__( - self, - library: Optional[str] = None, - parameters: Optional[str] = None, - strength: float = 1.0, - domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, - time_bounds: Optional[Sequence[float]] = None, - energy_bounds: Optional[Sequence[float]] = None, - rejection_strategy: Optional[str] = None, + self, + library: Optional[str] = None, + parameters: Optional[str] = None, + strength: float = 1.0, + constraints: Optional[Dict[str, Any]] = None ) -> None: - super().__init__( - strength=strength, domains=domains, time_bounds=time_bounds, - energy_bounds=energy_bounds, rejection_strategy=rejection_strategy - ) + super().__init__(strength=strength, constraints=constraints) self._library = None if library is not None: @@ -809,7 +747,7 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.CompiledSource: Source generated from XML element """ - kwargs = cls._get_common_arguments(elem) + kwargs = {'constraints': cls._get_constraints(elem)} kwargs['library'] = get_text(elem, 'library') source = cls(**kwargs) @@ -836,13 +774,17 @@ class FileSource(SourceBase): Path to the source file from which sites should be sampled strength : float Strength of the source (default is 1.0) - 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. - time_bounds : Iterable of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : Iterable of float - lower, upper bounds in energy [eV] of accepted particles + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', + the corresponding value is an iterable of :class:`openmc.Cell`, + :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled + sites must be within. For 'time_bounds' and 'energy_bounds', the + corresponding value is a sequence of floats giving the lower and upper + bounds on time in [s] or energy in [eV] that the sampled particle must + be within. The 'rejection_strategy' indicates what should happen when a + source particle is rejected: either 'resample' (pick a new particle) or + 'kill' (accept and terminate). Attributes ---------- @@ -852,34 +794,20 @@ class FileSource(SourceBase): Strength of the source type : str Indicator of source type: 'file' - ids : Iterable of int - IDs of domains to use for rejection - domain_type : {'cell', 'material', 'universe'} - Type of domain to use for rejection - domain_ids : iterable of int - Ids of domains to reject based on. - time_bounds : Iterable of float - lower, upper bounds in time [s] of accepted particles - energy_bounds : Iterable of float - lower, upper bounds in energy [eV] of accepted particles - rejection_strategy : {'resample', 'kill'} - Strategy when a particle is rejected. Pick a new particle ('resample') - or accept and terminate ('kill'). + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and + 'rejection_strategy'. + """ def __init__( - self, - path: Optional[PathLike] = None, - strength: float = 1.0, - domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, - time_bounds: Optional[Sequence[float]] = None, - energy_bounds: Optional[Sequence[float]] = None, - rejection_strategy = None, - ) -> None: - super().__init__( - strength=strength, domains=domains, time_bounds=time_bounds, - energy_bounds=energy_bounds, rejection_strategy=rejection_strategy - ) + self, + path: Optional[PathLike] = None, + strength: float = 1.0, + constraints: Optional[Dict[str, Any]] = None + ): + super().__init__(strength=strength, constraints=constraints) self._path = None if path is not None: self.path = path @@ -927,7 +855,7 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: Source generated from XML element """ - kwargs = cls._get_common_arguments(elem) + kwargs = {'constraints': cls._get_constraints(elem)} kwargs['path'] = get_text(elem, 'file') strength = get_text(elem, 'strength') if strength is not None: diff --git a/src/source.cpp b/src/source.cpp index 1f4a2b33d0b..286240a9d8f 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -79,7 +79,7 @@ unique_ptr Source::create(pugi::xml_node node) } } -void Source::check_for_restriction_nodes(pugi::xml_node node) +void Source::check_for_constraints(pugi::xml_node node) { // Check for domains to reject from if (check_for_node(node, "domain_type")) { @@ -91,8 +91,8 @@ void Source::check_for_restriction_nodes(pugi::xml_node node) } else if (domain_type == "universe") { domain_type_ = DomainType::UNIVERSE; } else { - fatal_error(std::string( - "Unrecognized domain type for source rejection: " + domain_type)); + fatal_error( + std::string("Unrecognized domain type for constraint: " + domain_type)); } auto ids = get_node_array(node, "domain_ids"); @@ -136,7 +136,7 @@ SourceSite Source::sample(uint64_t* seed) const // Sample a particle randomly from list site = this->sample_no_rejection(seed); - found = inside_bounds(site); + found = satisfies_constraints(site); if (!found && rejection_strategy_ == RejectionStrategy::KILL) { // Accept particle regardless but set wgt=0 to trigger garbage collection // I.e. kill this particle immediately and pick the next @@ -147,23 +147,24 @@ SourceSite Source::sample(uint64_t* seed) const return site; } -bool Source::inside_bounds(SourceSite& s) const +bool Source::satisfies_constraints(SourceSite& s) const { - return inside_spatial_bounds(s.r) && inside_energy_bounds(s.E) && - inside_time_bounds(s.time); + return satisfies_spatial_constraints(s.r) && + satisfies_energy_constraints(s.E) && + satisfies_time_constraints(s.time); } -bool Source::inside_energy_bounds(double E) const +bool Source::satisfies_energy_constraints(double E) const { return E > energy_bounds_.first && E < energy_bounds_.second; } -bool Source::inside_time_bounds(const double time) const +bool Source::satisfies_time_constraints(const double time) const { return time > time_bounds_.first && time < time_bounds_.second; } -bool Source::inside_spatial_bounds(Position r) const +bool Source::satisfies_spatial_constraints(Position r) const { bool found = false; GeometryState geom_state; @@ -264,7 +265,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) } // Check for additional defined restrictions - check_for_restriction_nodes(node); + check_for_constraints(node); } } @@ -284,7 +285,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const site.r = space_->sample(seed); // Check if otherwise outside defined restriction bounds - found = inside_spatial_bounds(site.r); + found = satisfies_spatial_constraints(site.r); // Check if spatial site is in fissionable material if (found) { @@ -340,7 +341,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const site.E = energy_->sample(seed); // Resample if energy falls above maximum particle energy - if (site.E < data::energy_max[p] and (inside_energy_bounds(site.E))) + if (site.E < data::energy_max[p] and + (satisfies_energy_constraints(site.E))) break; n_reject++; @@ -374,7 +376,7 @@ FileSource::FileSource(pugi::xml_node node) } else { this->load_sites_from_file(path); } - check_for_restriction_nodes(node); + check_for_constraints(node); } FileSource::FileSource(const std::string& path) @@ -509,7 +511,7 @@ MeshSource::MeshSource(pugi::xml_node node) space_ = std::make_unique(mesh_idx, strengths); // Check for additional defined restrictions - check_for_restriction_nodes(node); + check_for_constraints(node); } SourceSite MeshSource::sample(uint64_t* seed) const @@ -521,7 +523,7 @@ SourceSite MeshSource::sample(uint64_t* seed) const Position r; while (true) { r = space_->mesh()->sample_element(element, seed); - if (this->inside_spatial_bounds(r)) { + if (this->satisfies_spatial_constraints(r)) { break; } } @@ -533,7 +535,8 @@ SourceSite MeshSource::sample(uint64_t* seed) const site.r = r; // Apply other rejections - if (inside_energy_bounds(site.E) && inside_time_bounds(site.time)) { + if (satisfies_energy_constraints(site.E) && + satisfies_time_constraints(site.time)) { break; } } From 99eb170d68bb3ca4e15318236111f43f01e5dd52 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 May 2024 12:36:24 -0500 Subject: [PATCH 44/57] Check for constraints/strength in Source constructor --- include/openmc/source.h | 19 +++++++++++-------- src/source.cpp | 37 +++++++++++++++++++------------------ 2 files changed, 30 insertions(+), 26 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index e7e0d0598aa..1b9fbcf4370 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -44,30 +44,35 @@ extern vector> external_sources; class Source { public: + // Constructors, destructors + Source() = default; + explicit Source(pugi::xml_node node); virtual ~Source() = default; // Methods that can be overridden virtual SourceSite sample(uint64_t* seed) const; - virtual double strength() const { return 1.0; } + virtual double strength() const { return strength_; } static unique_ptr create(pugi::xml_node node); protected: - // Domain types - enum class DomainType { UNIVERSE, MATERIAL, CELL }; - enum class RejectionStrategy { KILL, RESAMPLE }; - // Methods that can be overridden virtual SourceSite sample_no_rejection(uint64_t* seed) const = 0; - // Methods + // Methods for constraints void check_for_constraints(pugi::xml_node node); bool satisfies_constraints(SourceSite& s) const; bool satisfies_spatial_constraints(Position r) const; bool satisfies_energy_constraints(const double E) const; bool satisfies_time_constraints(const double time) const; +private: + // Domain types + enum class DomainType { UNIVERSE, MATERIAL, CELL }; + enum class RejectionStrategy { KILL, RESAMPLE }; + // Data members + double strength_ {1.0}; //!< Source strength std::unordered_set domain_ids_; //!< Domains to reject from DomainType domain_type_; //!< Domain type for rejection std::pair time_bounds_ {-std::numeric_limits::max(), @@ -96,7 +101,6 @@ class IndependentSource : public Source { // Properties ParticleType particle_type() const { return particle_; } - double strength() const override { return strength_; } // Make observing pointers available SpatialDistribution* space() const { return space_.get(); } @@ -112,7 +116,6 @@ class IndependentSource : public Source { private: // 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 diff --git a/src/source.cpp b/src/source.cpp index 286240a9d8f..a0caba9c339 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -47,9 +47,20 @@ vector> external_sources; } //============================================================================== -// Source create implementation +// Source implementation //============================================================================== +Source::Source(pugi::xml_node node) +{ + // Check for source strength + if (check_for_node(node, "strength")) { + strength_ = std::stod(get_node_value(node, "strength")); + } + + // Check for additional defined constraints + check_for_constraints(node); +} + unique_ptr Source::create(pugi::xml_node node) { // if the source type is present, use it to determine the type @@ -204,7 +215,7 @@ IndependentSource::IndependentSource( energy_ {std::move(energy)}, time_ {std::move(time)} {} -IndependentSource::IndependentSource(pugi::xml_node node) +IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) { // Check for particle type if (check_for_node(node, "particle")) { @@ -219,11 +230,6 @@ IndependentSource::IndependentSource(pugi::xml_node node) } } - // Check for source strength - if (check_for_node(node, "strength")) { - strength_ = std::stod(get_node_value(node, "strength")); - } - // Check for external source file if (check_for_node(node, "file")) { @@ -263,9 +269,6 @@ IndependentSource::IndependentSource(pugi::xml_node node) double p[] {1.0}; time_ = UPtrDist {new Discrete {T, p, 1}}; } - - // Check for additional defined restrictions - check_for_constraints(node); } } @@ -284,7 +287,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Sample spatial distribution site.r = space_->sample(seed); - // Check if otherwise outside defined restriction bounds + // Check if sampled position satisfies spatial constraints found = satisfies_spatial_constraints(site.r); // Check if spatial site is in fissionable material @@ -368,7 +371,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const //============================================================================== // FileSource implementation //============================================================================== -FileSource::FileSource(pugi::xml_node node) + +FileSource::FileSource(pugi::xml_node node) : Source(node) { auto path = get_node_value(node, "file", false, true); if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { @@ -376,7 +380,6 @@ FileSource::FileSource(pugi::xml_node node) } else { this->load_sites_from_file(path); } - check_for_constraints(node); } FileSource::FileSource(const std::string& path) @@ -422,7 +425,8 @@ SourceSite FileSource::sample_no_rejection(uint64_t* seed) const //============================================================================== // CompiledSourceWrapper implementation //============================================================================== -CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) + +CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) : Source(node) { // Get shared library path and parameters auto path = get_node_value(node, "library", false, true); @@ -486,7 +490,7 @@ CompiledSourceWrapper::~CompiledSourceWrapper() // MeshSource implementation //============================================================================== -MeshSource::MeshSource(pugi::xml_node node) +MeshSource::MeshSource(pugi::xml_node node) : Source(node) { int32_t mesh_id = stoi(get_node_value(node, "mesh")); int32_t mesh_idx = model::mesh_map.at(mesh_id); @@ -509,9 +513,6 @@ MeshSource::MeshSource(pugi::xml_node node) } space_ = std::make_unique(mesh_idx, strengths); - - // Check for additional defined restrictions - check_for_constraints(node); } SourceSite MeshSource::sample(uint64_t* seed) const From 31ad069aeb58fde74f4a6505b680be90a2ebae99 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 May 2024 21:36:11 -0500 Subject: [PATCH 45/57] Extend tests for source constraints --- tests/unit_tests/test_source.py | 82 ++++++++++++++++++++++++++++++--- 1 file changed, 75 insertions(+), 7 deletions(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 4783ff8f11d..729531c22c8 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -100,7 +100,8 @@ def test_source_xml_roundtrip(): assert new_src.strength == approx(src.strength) -def test_rejection(run_in_tmpdir): +@pytest.fixture +def sphere_box_model(): # Model with two spheres inside a box mat = openmc.Material() mat.add_nuclide('H1', 1.0) @@ -119,20 +120,87 @@ def test_rejection(run_in_tmpdir): model.settings.batches = 10 model.settings.run_mode = 'fixed source' + return model, cell1, cell2, cell3 + + +def test_constraints_independent(sphere_box_model, run_in_tmpdir): + model, cell1, cell2, cell3 = sphere_box_model + # Set up a box source with rejection on the spherical cell space = openmc.stats.Box(*cell3.bounding_box) - model.settings.source = openmc.IndependentSource(space=space, domains=[cell1, cell2]) + model.settings.source = openmc.IndependentSource( + space=space, constraints={'domains': [cell1, cell2]} + ) # Load up model via openmc.lib and sample source - model.export_to_xml() + model.export_to_model_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 - assert p.r not in non_source_region + assert p.r in (cell1.region | cell2.region) + assert p.r not in cell3.region + + openmc.lib.finalize() + + +def test_constraints_mesh(sphere_box_model, run_in_tmpdir): + model, cell1, cell2, cell3 = sphere_box_model + + bbox = cell3.bounding_box + mesh = openmc.RegularMesh() + mesh.lower_left = bbox.lower_left + mesh.upper_right = bbox.upper_right + mesh.dimension = (2, 1, 1) + + left_source = openmc.IndependentSource() + right_source = openmc.IndependentSource() + model.settings.source = openmc.MeshSource( + mesh, [left_source, right_source], constraints={'domains': [cell1, cell2]} + ) + + # Load up model via openmc.lib and sample source + model.export_to_model_xml() + openmc.lib.init() + particles = openmc.lib.sample_external_source(1000) + + # Make sure that all sampled sources are within one of the spheres + for p in particles: + assert p.r in (cell1.region | cell2.region) + assert p.r not in cell3.region + + openmc.lib.finalize() + + +def test_constraints_file(sphere_box_model, run_in_tmpdir): + model = sphere_box_model[0] + + # Create source file with randomly sampled source sites + rng = np.random.default_rng() + energy = rng.uniform(0., 1e6, 10_000) + time = rng.uniform(0., 1., 10_000) + particles = [openmc.SourceParticle(E=e, time=t) for e, t in zip(energy, time)] + openmc.write_source_file(particles, 'uniform_source.h5') + + # Use source file + model.settings.source = openmc.FileSource( + 'uniform_source.h5', + constraints={ + 'time_bounds': [0.25, 0.75], + 'energy_bounds': [500.e3, 1.0e6], + } + ) + + # Load up model via openmc.lib and sample source + model.export_to_model_xml() + openmc.lib.init() + particles = openmc.lib.sample_external_source(1000) + + # Make sure that all sampled sources are within energy/time bounds + for p in particles: + assert 0.25 <= p.time <= 0.75 + assert 500.e3 <= p.E <= 1.0e6 openmc.lib.finalize() @@ -156,4 +224,4 @@ def test_exceptions(): with pytest.raises(AttributeError, match=r'has no attribute \'frisbee\''): s = openmc.IndependentSource() - s.frisbee \ No newline at end of file + s.frisbee From 0e2f7473c8682999a63fdeb1ee015b1d1c977bf4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 May 2024 18:44:08 -0500 Subject: [PATCH 46/57] Respond to many @gridley comments --- include/openmc/source.h | 40 +++++++++--------- src/source.cpp | 91 ++++++++++++++++++++--------------------- 2 files changed, 65 insertions(+), 66 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 1b9fbcf4370..f7335c67e7f 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -56,19 +56,30 @@ class Source { static unique_ptr create(pugi::xml_node node); protected: - // Methods that can be overridden - virtual SourceSite sample_no_rejection(uint64_t* seed) const = 0; + //! Sample a source site (without applying constraints) + // + //! Sample from the external source distribution + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + virtual SourceSite sample_without_constraints(uint64_t* seed) const + { + return {}; + } // Methods for constraints - void check_for_constraints(pugi::xml_node node); - bool satisfies_constraints(SourceSite& s) const; + void read_constraints(pugi::xml_node node); bool satisfies_spatial_constraints(Position r) const; - bool satisfies_energy_constraints(const double E) const; - bool satisfies_time_constraints(const double time) const; + bool satisfies_energy_constraints(double E) const; + bool satisfies_time_constraints(double time) const; private: // Domain types enum class DomainType { UNIVERSE, MATERIAL, CELL }; + + // Strategy used for rejecting sites when constraints are applied. KILL means + // that sites are always accepted but if they don't satisfy constraints, they + // are given weight 0. RESAMPLE means that a new source site will be sampled + // until constraints are met. enum class RejectionStrategy { KILL, RESAMPLE }; // Data members @@ -79,7 +90,8 @@ class Source { std::numeric_limits::max()}; //!< time limits std::pair energy_bounds_ { 0, std::numeric_limits::max()}; //!< energy limits - RejectionStrategy rejection_strategy_; //!< Procedure for rejected + RejectionStrategy rejection_strategy_ { + RejectionStrategy::RESAMPLE}; //!< Procedure for rejecting }; //============================================================================== @@ -108,11 +120,6 @@ class IndependentSource : public Source { Distribution* energy() const { return energy_.get(); } Distribution* time() const { return time_.get(); } -protected: - // Note that this method is not used since IndependentSource directly - // overrides the sample method - SourceSite sample_no_rejection(uint64_t* seed) const override { return {}; } - private: // Data members ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted @@ -137,7 +144,7 @@ class FileSource : public Source { const std::string& path); //!< Load source sites from file protected: - SourceSite sample_no_rejection(uint64_t* seed) const override; + SourceSite sample_without_constraints(uint64_t* seed) const override; private: vector sites_; //!< Source sites from a file @@ -159,7 +166,7 @@ class CompiledSourceWrapper : public Source { protected: // Defer implementation to custom source library - SourceSite sample_no_rejection(uint64_t* seed) const override + SourceSite sample_without_constraints(uint64_t* seed) const override { return compiled_source_->sample(seed); } @@ -194,11 +201,6 @@ class MeshSource : public Source { return sources_.size() == 1 ? sources_[0] : sources_[i]; } -protected: - // Note that this method is not used since MeshSource directly overrides the - // sample method - SourceSite sample_no_rejection(uint64_t* seed) const override { return {}; } - private: // Data members unique_ptr space_; //!< Mesh spatial diff --git a/src/source.cpp b/src/source.cpp index a0caba9c339..7cc10fdd669 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -55,10 +55,13 @@ Source::Source(pugi::xml_node node) // Check for source strength if (check_for_node(node, "strength")) { strength_ = std::stod(get_node_value(node, "strength")); + if (strength_ < 0.0) { + fatal_error("Source strength is negative."); + } } // Check for additional defined constraints - check_for_constraints(node); + read_constraints(node); } unique_ptr Source::create(pugi::xml_node node) @@ -90,7 +93,7 @@ unique_ptr Source::create(pugi::xml_node node) } } -void Source::check_for_constraints(pugi::xml_node node) +void Source::read_constraints(pugi::xml_node node) { // Check for domains to reject from if (check_for_node(node, "domain_type")) { @@ -112,13 +115,17 @@ void Source::check_for_constraints(pugi::xml_node node) if (check_for_node(node, "time_bounds")) { auto ids = get_node_array(node, "time_bounds"); - if (ids.size() == 2) - time_bounds_ = std::make_pair(ids[0], ids[1]); + if (ids.size() != 2) { + fatal_error("Time bounds must be represented by two numbers."); + } + time_bounds_ = std::make_pair(ids[0], ids[1]); } if (check_for_node(node, "energy_bounds")) { auto ids = get_node_array(node, "energy_bounds"); - if (ids.size() == 2) - energy_bounds_ = std::make_pair(ids[0], ids[1]); + if (ids.size() != 2) { + fatal_error("Energy bounds must be represented by two numbers."); + } + energy_bounds_ = std::make_pair(ids[0], ids[1]); } // Check for how to handle rejected particles @@ -132,77 +139,70 @@ void Source::check_for_constraints(pugi::xml_node node) fatal_error(std::string( "Unrecognized strategy source rejection: " + rejection_strategy)); } - } else { - // Default to resample rejected particles - rejection_strategy_ = RejectionStrategy::RESAMPLE; } } SourceSite Source::sample(uint64_t* seed) const { - bool found = false; + bool accepted = false; SourceSite site; - while (!found) { - // Sample a particle randomly from list - site = this->sample_no_rejection(seed); - - found = satisfies_constraints(site); - if (!found && rejection_strategy_ == RejectionStrategy::KILL) { - // Accept particle regardless but set wgt=0 to trigger garbage collection - // I.e. kill this particle immediately and pick the next - found = true; + while (!accepted) { + // Sample a source site without considering constraints yet + site = this->sample_without_constraints(seed); + + // Check whether sampled site satisfies constraints + accepted = satisfies_spatial_constraints(site.r) && + satisfies_energy_constraints(site.E) && + satisfies_time_constraints(site.time); + if (!accepted && rejection_strategy_ == RejectionStrategy::KILL) { + // Accept particle but set weight to 0 so that it is killed immediately + accepted = true; site.wgt = 0.0; } } return site; } -bool Source::satisfies_constraints(SourceSite& s) const -{ - return satisfies_spatial_constraints(s.r) && - satisfies_energy_constraints(s.E) && - satisfies_time_constraints(s.time); -} - bool Source::satisfies_energy_constraints(double E) const { return E > energy_bounds_.first && E < energy_bounds_.second; } -bool Source::satisfies_time_constraints(const double time) const +bool Source::satisfies_time_constraints(double time) const { return time > time_bounds_.first && time < time_bounds_.second; } bool Source::satisfies_spatial_constraints(Position r) const { - bool found = false; GeometryState geom_state; geom_state.r() = r; // Reject particle if it's not in the geometry at all - if (!(found = exhaustive_find_cell(geom_state))) + bool found = exhaustive_find_cell(geom_state); + if (!found) return false; + // Check the geometry state against the specified constraints + bool accepted = true; if (!domain_ids_.empty()) { if (domain_type_ == DomainType::MATERIAL) { auto mat_index = geom_state.material(); if (mat_index != MATERIAL_VOID) { - found = contains(domain_ids_, model::materials[mat_index]->id()); + accepted = contains(domain_ids_, model::materials[mat_index]->id()); } } else { for (int i = 0; i < geom_state.n_coord(); i++) { auto id = (domain_type_ == DomainType::CELL) ? model::cells[geom_state.coord(i).cell]->id_ : model::universes[geom_state.coord(i).universe]->id_; - if ((found = contains(domain_ids_, id))) + if ((accepted = contains(domain_ids_, id))) break; } } } - // If at this point found is false, part. is outside all cells/materials - return found; + return accepted; } //============================================================================== @@ -277,21 +277,21 @@ 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; + // Repeat sampling source location until a good site has been accepted + bool accepted = false; int n_reject = 0; static int n_accept = 0; - while (!found) { + while (!accepted) { // Sample spatial distribution site.r = space_->sample(seed); // Check if sampled position satisfies spatial constraints - found = satisfies_spatial_constraints(site.r); + accepted = satisfies_spatial_constraints(site.r); // Check if spatial site is in fissionable material - if (found) { + if (accepted) { // Create a temporary particle and search for material at the site Particle p; p.r() = site.r; @@ -303,16 +303,16 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Determine material auto mat_index = p.material(); if (mat_index == MATERIAL_VOID) { - found = false; + accepted = false; } else { - found = model::materials[mat_index]->fissionable(); + accepted = model::materials[mat_index]->fissionable(); } } } } // Check for rejection - if (!found) { + if (!accepted) { ++n_reject; if (n_reject >= EXTSRC_REJECT_THRESHOLD && static_cast(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) { @@ -415,7 +415,7 @@ void FileSource::load_sites_from_file(const std::string& path) file_close(file_id); } -SourceSite FileSource::sample_no_rejection(uint64_t* seed) const +SourceSite FileSource::sample_without_constraints(uint64_t* seed) const { // Sample a particle randomly from list size_t i_site = sites_.size() * prn(seed); @@ -522,12 +522,9 @@ SourceSite MeshSource::sample(uint64_t* seed) const // Sample position and apply rejection on spatial domains Position r; - while (true) { + do { r = space_->mesh()->sample_element(element, seed); - if (this->satisfies_spatial_constraints(r)) { - break; - } - } + } while (!this->satisfies_spatial_constraints(r)); SourceSite site; while (true) { From cd04ecb39c4e56ccb672c38938e11c75781e38c6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 May 2024 07:26:23 -0500 Subject: [PATCH 47/57] Add fissionable source constraint --- examples/assembly/assembly.py | 9 +- .../pincell_depletion/restart_depletion.py | 5 +- examples/pincell_depletion/run_depletion.py | 5 +- include/openmc/source.h | 4 +- openmc/examples.py | 12 +- openmc/source.py | 130 ++++++++++-------- openmc/stats/multivariate.py | 11 ++ src/source.cpp | 48 ++++--- .../asymmetric_lattice/inputs_true.dat | 3 +- .../asymmetric_lattice/test.py | 6 +- .../mg_temperature_multi/inputs_true.dat | 3 +- .../mg_temperature_multi/test.py | 5 +- .../surface_tally/inputs_true.dat | 3 +- tests/regression_tests/surface_tally/test.py | 6 +- tests/unit_tests/test_lib.py | 6 +- tests/unit_tests/test_model.py | 6 +- 16 files changed, 156 insertions(+), 106 deletions(-) diff --git a/examples/assembly/assembly.py b/examples/assembly/assembly.py index a370a361144..d94cb69bfc4 100644 --- a/examples/assembly/assembly.py +++ b/examples/assembly/assembly.py @@ -112,11 +112,10 @@ def assembly_model(): model.settings.batches = 150 model.settings.inactive = 50 model.settings.particles = 1000 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - (-pitch/2, -pitch/2, -1), - (pitch/2, pitch/2, 1), - only_fissionable=True - )) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box((-pitch/2, -pitch/2, -1), (pitch/2, pitch/2, 1)), + constraints={'fissionable': True} + ) # NOTE: We never actually created a Materials object. When you export/run # using the Model object, if no materials were assigned it will look through diff --git a/examples/pincell_depletion/restart_depletion.py b/examples/pincell_depletion/restart_depletion.py index ac065d3402f..de9fc16cb1f 100644 --- a/examples/pincell_depletion/restart_depletion.py +++ b/examples/pincell_depletion/restart_depletion.py @@ -26,8 +26,9 @@ # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] -uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings.source = openmc.IndependentSource(space=uniform_dist) +uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) +settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] diff --git a/examples/pincell_depletion/run_depletion.py b/examples/pincell_depletion/run_depletion.py index e1959938f05..013c83c86a5 100644 --- a/examples/pincell_depletion/run_depletion.py +++ b/examples/pincell_depletion/run_depletion.py @@ -71,8 +71,9 @@ # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] -uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings.source = openmc.IndependentSource(space=uniform_dist) +uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) +settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] diff --git a/include/openmc/source.h b/include/openmc/source.h index f7335c67e7f..d3b5c76d8ae 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -72,7 +72,7 @@ class Source { bool satisfies_energy_constraints(double E) const; bool satisfies_time_constraints(double time) const; -private: +protected: // Domain types enum class DomainType { UNIVERSE, MATERIAL, CELL }; @@ -90,6 +90,8 @@ class Source { std::numeric_limits::max()}; //!< time limits std::pair energy_bounds_ { 0, std::numeric_limits::max()}; //!< energy limits + bool only_fissionable_ { + false}; //!< Whether site must be in fissionable material RejectionStrategy rejection_strategy_ { RejectionStrategy::RESAMPLE}; //!< Procedure for rejecting }; diff --git a/openmc/examples.py b/openmc/examples.py index cef86a47a87..fc94b8b528e 100644 --- a/openmc/examples.py +++ b/openmc/examples.py @@ -76,8 +76,10 @@ def pwr_pin_cell(): model.settings.batches = 10 model.settings.inactive = 5 model.settings.particles = 100 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True)) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]), + constraints={'fissionable': True} + ) plot = openmc.Plot.from_geometry(model.geometry) plot.pixels = (300, 300) @@ -527,8 +529,10 @@ def pwr_assembly(): model.settings.batches = 10 model.settings.inactive = 5 model.settings.particles = 100 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True)) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]), + constraints={'fissionable': True} + ) plot = openmc.Plot() plot.origin = (0.0, 0.0, 0) diff --git a/openmc/source.py b/openmc/source.py index 90d5ff7fe6e..6380cadd56e 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -30,15 +30,17 @@ class SourceBase(ABC): Strength of the source constraints : dict Constraints on sampled source particles. Valid keys include 'domains', - 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', - the corresponding value is an iterable of :class:`openmc.Cell`, - :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled - sites must be within. For 'time_bounds' and 'energy_bounds', the - corresponding value is a sequence of floats giving the lower and upper - bounds on time in [s] or energy in [eV] that the sampled particle must - be within. The 'rejection_strategy' indicates what should happen when a - source particle is rejected: either 'resample' (pick a new particle) or - 'kill' (accept and terminate). + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -48,8 +50,8 @@ class SourceBase(ABC): Strength of the source constraints : dict Constraints on sampled source particles. Valid keys include - 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and - 'rejection_strategy'. + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ @@ -99,6 +101,9 @@ def constraints(self, constraints: Optional[Dict[str, Any]]): elif key == 'energy_bounds': cv.check_type('energy bounds', value, Iterable, Real) self._constraints['energy_bounds'] = tuple(value) + elif key == 'fissionable': + cv.check_type('fissionable', value, bool) + self._constraints['fissionable'] = value elif key == 'rejection_strategy': cv.check_value('rejection strategy', value, ('resample', 'kill')) self._constraints['rejection_strategy'] = value @@ -142,6 +147,9 @@ def to_xml_element(self) -> ET.Element: if "energy_bounds" in constraints: dt_elem = ET.SubElement(element, "energy_bounds") dt_elem.text = ' '.join(str(E) for E in constraints["energy_bounds"]) + if "fissionable" in constraints: + dt_elem = ET.SubElement(element, "fissionable") + dt_elem.text = str(constraints["fissionable"]).lower() if "rejection_strategy" in constraints: dt_elem = ET.SubElement(element, "rejection_strategy") dt_elem.text = constraints["rejection_strategy"] @@ -216,6 +224,10 @@ def _get_constraints(elem: ET.Element) -> Dict[str, Any]: if energy_bounds is not None: constraints['energy_bounds'] = [float(x) for x in energy_bounds.split()] + fissionable = get_text(elem, "fissionable") + if fissionable is not None: + constraints['fissionable'] = fissionable in ('true', '1') + rejection_strategy = get_text(elem, "rejection_strategy") if rejection_strategy is not None: constraints['rejection_strategy'] = rejection_strategy @@ -250,15 +262,17 @@ class IndependentSource(SourceBase): Use the `constraints` argument instead. constraints : dict Constraints on sampled source particles. Valid keys include 'domains', - 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', - the corresponding value is an iterable of :class:`openmc.Cell`, - :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled - sites must be within. For 'time_bounds' and 'energy_bounds', the - corresponding value is a sequence of floats giving the lower and upper - bounds on time in [s] or energy in [eV] that the sampled particle must - be within. The 'rejection_strategy' indicates what should happen when a - source particle is rejected: either 'resample' (pick a new particle) or - 'kill' (accept and terminate). + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -281,8 +295,8 @@ class IndependentSource(SourceBase): Source particle type constraints : dict Constraints on sampled source particles. Valid keys include - 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and - 'rejection_strategy'. + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ @@ -476,15 +490,17 @@ class MeshSource(SourceBase): during source site sampling. constraints : dict Constraints on sampled source particles. Valid keys include 'domains', - 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', - the corresponding value is an iterable of :class:`openmc.Cell`, - :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled - sites must be within. For 'time_bounds' and 'energy_bounds', the - corresponding value is a sequence of floats giving the lower and upper - bounds on time in [s] or energy in [eV] that the sampled particle must - be within. The 'rejection_strategy' indicates what should happen when a - source particle is rejected: either 'resample' (pick a new particle) or - 'kill' (accept and terminate). + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -498,8 +514,8 @@ class MeshSource(SourceBase): Indicator of source type: 'mesh' constraints : dict Constraints on sampled source particles. Valid keys include - 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and - 'rejection_strategy'. + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ def __init__( @@ -650,15 +666,17 @@ class CompiledSource(SourceBase): Strength of the source constraints : dict Constraints on sampled source particles. Valid keys include 'domains', - 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', - the corresponding value is an iterable of :class:`openmc.Cell`, - :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled - sites must be within. For 'time_bounds' and 'energy_bounds', the - corresponding value is a sequence of floats giving the lower and upper - bounds on time in [s] or energy in [eV] that the sampled particle must - be within. The 'rejection_strategy' indicates what should happen when a - source particle is rejected: either 'resample' (pick a new particle) or - 'kill' (accept and terminate). + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -672,8 +690,8 @@ class CompiledSource(SourceBase): Indicator of source type: 'compiled' constraints : dict Constraints on sampled source particles. Valid keys include - 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and - 'rejection_strategy'. + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ def __init__( @@ -776,15 +794,17 @@ class FileSource(SourceBase): Strength of the source (default is 1.0) constraints : dict Constraints on sampled source particles. Valid keys include 'domains', - 'time_bounds', 'energy_bounds', and 'rejection_strategy'. For 'domains', - the corresponding value is an iterable of :class:`openmc.Cell`, - :class:`openmc.Material`, or :class:`openmc.Universe` for which sampled - sites must be within. For 'time_bounds' and 'energy_bounds', the - corresponding value is a sequence of floats giving the lower and upper - bounds on time in [s] or energy in [eV] that the sampled particle must - be within. The 'rejection_strategy' indicates what should happen when a - source particle is rejected: either 'resample' (pick a new particle) or - 'kill' (accept and terminate). + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -796,8 +816,8 @@ class FileSource(SourceBase): Indicator of source type: 'file' constraints : dict Constraints on sampled source particles. Valid keys include - 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', and - 'rejection_strategy'. + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 3922d601aa8..212083c3cff 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -4,6 +4,7 @@ from collections.abc import Iterable from math import cos, pi from numbers import Real +from warnings import warn import lxml.etree as ET import numpy as np @@ -768,6 +769,9 @@ class Box(Spatial): Whether spatial sites should only be accepted if they occur in fissionable materials + .. deprecated:: 0.14.1 + Use the `constraints` argument when defining a source object instead. + Attributes ---------- lower_left : Iterable of float @@ -778,6 +782,9 @@ class Box(Spatial): Whether spatial sites should only be accepted if they occur in fissionable materials + .. deprecated:: 0.14.1 + Use the `constraints` argument when defining a source object instead. + """ def __init__( @@ -818,6 +825,10 @@ def only_fissionable(self): def only_fissionable(self, only_fissionable): cv.check_type('only fissionable', only_fissionable, bool) self._only_fissionable = only_fissionable + if only_fissionable: + warn("The 'only_fissionable' has been deprecated. Use the " + "'constraints' argument when defining a source instead.", + FutureWarning) def to_xml_element(self): """Return XML representation of the box distribution diff --git a/src/source.cpp b/src/source.cpp index 7cc10fdd669..7995affabae 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -128,6 +128,10 @@ void Source::read_constraints(pugi::xml_node node) energy_bounds_ = std::make_pair(ids[0], ids[1]); } + if (check_for_node(node, "fissionable")) { + only_fissionable_ = get_node_value_bool(node, "fissionable"); + } + // Check for how to handle rejected particles if (check_for_node(node, "rejection_strategy")) { std::string rejection_strategy = get_node_value(node, "rejection_strategy"); @@ -184,7 +188,7 @@ bool Source::satisfies_spatial_constraints(Position r) const if (!found) return false; - // Check the geometry state against the specified constraints + // Check the geometry state against specified domains bool accepted = true; if (!domain_ids_.empty()) { if (domain_type_ == DomainType::MATERIAL) { @@ -202,6 +206,18 @@ bool Source::satisfies_spatial_constraints(Position r) const } } } + + // Check if spatial site is in fissionable material + if (accepted && only_fissionable_) { + // Determine material + auto mat_index = geom_state.material(); + if (mat_index == MATERIAL_VOID) { + accepted = false; + } else { + accepted = model::materials[mat_index]->fissionable(); + } + } + return accepted; } @@ -243,6 +259,15 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) space_ = UPtrSpace {new SpatialPoint()}; } + // For backwards compatibility, check for only fissionable setting on box + // source + auto space_box = dynamic_cast(space_.get()); + if (space_box) { + if (!only_fissionable_) { + only_fissionable_ = space_box->only_fissionable(); + } + } + // Determine external source angular distribution if (check_for_node(node, "angle")) { angle_ = UnitSphereDistribution::create(node.child("angle")); @@ -290,27 +315,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Check if sampled position satisfies spatial constraints accepted = satisfies_spatial_constraints(site.r); - // Check if spatial site is in fissionable material - if (accepted) { - // Create a temporary particle and search for material at the site - Particle p; - p.r() = site.r; - exhaustive_find_cell(p); - - auto space_box = dynamic_cast(space_.get()); - if (space_box) { - if (space_box->only_fissionable()) { - // Determine material - auto mat_index = p.material(); - if (mat_index == MATERIAL_VOID) { - accepted = false; - } else { - accepted = model::materials[mat_index]->fissionable(); - } - } - } - } - // Check for rejection if (!accepted) { ++n_reject; diff --git a/tests/regression_tests/asymmetric_lattice/inputs_true.dat b/tests/regression_tests/asymmetric_lattice/inputs_true.dat index 3977220da61..d3d9ea3336e 100644 --- a/tests/regression_tests/asymmetric_lattice/inputs_true.dat +++ b/tests/regression_tests/asymmetric_lattice/inputs_true.dat @@ -209,9 +209,10 @@ 10 5 - + -32 -32 0 32 32 32 + true diff --git a/tests/regression_tests/asymmetric_lattice/test.py b/tests/regression_tests/asymmetric_lattice/test.py index 0d17a6c98eb..fadf272ffd3 100644 --- a/tests/regression_tests/asymmetric_lattice/test.py +++ b/tests/regression_tests/asymmetric_lattice/test.py @@ -54,8 +54,10 @@ def __init__(self, *args, **kwargs): self._model.tallies.append(tally) # Specify summary output and correct source sampling box - self._model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-32, -32, 0], [32, 32, 32], only_fissionable = True)) + self._model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-32, -32, 0], [32, 32, 32]), + constraints={'fissionable': True} + ) def _get_results(self, hash_output=True): """Digest info in statepoint and summary and return as a string.""" diff --git a/tests/regression_tests/mg_temperature_multi/inputs_true.dat b/tests/regression_tests/mg_temperature_multi/inputs_true.dat index b2ad5063c2b..a3879251d4c 100644 --- a/tests/regression_tests/mg_temperature_multi/inputs_true.dat +++ b/tests/regression_tests/mg_temperature_multi/inputs_true.dat @@ -28,9 +28,10 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + true multi-group diff --git a/tests/regression_tests/mg_temperature_multi/test.py b/tests/regression_tests/mg_temperature_multi/test.py index 404c8c0cd05..3117e29ba03 100755 --- a/tests/regression_tests/mg_temperature_multi/test.py +++ b/tests/regression_tests/mg_temperature_multi/test.py @@ -133,8 +133,9 @@ def test_mg_temperature_multi(): # Create an initial uniform spatial source distribution over fissionable zones lower_left = (-pitch/2, -pitch/2, -1) upper_right = (pitch/2, pitch/2, 1) - uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True) - settings.source = openmc.IndependentSource(space=uniform_dist) + uniform_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) ############################################################################### # Define tallies diff --git a/tests/regression_tests/surface_tally/inputs_true.dat b/tests/regression_tests/surface_tally/inputs_true.dat index abd01266d48..a6b30aeb221 100644 --- a/tests/regression_tests/surface_tally/inputs_true.dat +++ b/tests/regression_tests/surface_tally/inputs_true.dat @@ -30,9 +30,10 @@ 10 0 - + -0.62992 -0.62992 -1 0.62992 0.62992 1 + true diff --git a/tests/regression_tests/surface_tally/test.py b/tests/regression_tests/surface_tally/test.py index 8968db9b6fb..e496ac0f65c 100644 --- a/tests/regression_tests/surface_tally/test.py +++ b/tests/regression_tests/surface_tally/test.py @@ -72,9 +72,9 @@ def __init__(self, *args, **kwargs): # Create an initial uniform spatial source distribution bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] - uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:],\ - only_fissionable=True) - settings_file.source = openmc.IndependentSource(space=uniform_dist) + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:],) + settings_file.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) self._model.settings = settings_file # Tallies file diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 1310a6a7fa0..009a9096831 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -67,8 +67,10 @@ def uo2_trigger_model(): model.settings.batches = 10 model.settings.inactive = 5 model.settings.particles = 100 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-0.5, -0.5, -1], [0.5, 0.5, 1], only_fissionable=True)) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-0.5, -0.5, -1], [0.5, 0.5, 1]), + constraints={'fissionable': True}, + ) model.settings.verbosity = 1 model.settings.keff_trigger = {'type': 'std_dev', 'threshold': 0.001} model.settings.trigger_active = True diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 16fa18d453c..404235bef1f 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -57,9 +57,9 @@ def pin_model_attributes(): # Create a uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] - uniform_dist = openmc.stats.Box( - bounds[:3], bounds[3:], only_fissionable=True) - settings.source = openmc.IndependentSource(space=uniform_dist) + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) + settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] From 91cfdcccabd53d124f716965ba8a0ab50c5183e6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 May 2024 10:33:24 -0500 Subject: [PATCH 48/57] Rename to sample / sample_with_constraints --- include/openmc/source.h | 24 +++++++++++++----------- src/source.cpp | 14 +++++++------- 2 files changed, 20 insertions(+), 18 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index d3b5c76d8ae..98e8de94fe3 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -50,22 +50,24 @@ class Source { virtual ~Source() = default; // Methods that can be overridden - virtual SourceSite sample(uint64_t* seed) const; virtual double strength() const { return strength_; } - static unique_ptr create(pugi::xml_node node); + //! Sample a source site and apply constraints + // + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + virtual SourceSite sample_with_constraints(uint64_t* seed) const; -protected: //! Sample a source site (without applying constraints) // //! Sample from the external source distribution //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - virtual SourceSite sample_without_constraints(uint64_t* seed) const - { - return {}; - } + virtual SourceSite sample(uint64_t* seed) const { return {}; } + static unique_ptr create(pugi::xml_node node); + +protected: // Methods for constraints void read_constraints(pugi::xml_node node); bool satisfies_spatial_constraints(Position r) const; @@ -111,7 +113,7 @@ class IndependentSource : public Source { //! Sample from the external source distribution //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - SourceSite sample(uint64_t* seed) const override; + SourceSite sample_with_constraints(uint64_t* seed) const override; // Properties ParticleType particle_type() const { return particle_; } @@ -146,7 +148,7 @@ class FileSource : public Source { const std::string& path); //!< Load source sites from file protected: - SourceSite sample_without_constraints(uint64_t* seed) const override; + SourceSite sample(uint64_t* seed) const override; private: vector sites_; //!< Source sites from a file @@ -168,7 +170,7 @@ class CompiledSourceWrapper : public Source { protected: // Defer implementation to custom source library - SourceSite sample_without_constraints(uint64_t* seed) const override + SourceSite sample(uint64_t* seed) const override { return compiled_source_->sample(seed); } @@ -192,7 +194,7 @@ class MeshSource : public Source { //! Sample from the external source distribution //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - SourceSite sample(uint64_t* seed) const override; + SourceSite sample_with_constraints(uint64_t* seed) const override; // Properties double strength() const override { return space_->total_strength(); } diff --git a/src/source.cpp b/src/source.cpp index 7995affabae..914e584e8e2 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -146,14 +146,14 @@ void Source::read_constraints(pugi::xml_node node) } } -SourceSite Source::sample(uint64_t* seed) const +SourceSite Source::sample_with_constraints(uint64_t* seed) const { bool accepted = false; SourceSite site; while (!accepted) { // Sample a source site without considering constraints yet - site = this->sample_without_constraints(seed); + site = this->sample(seed); // Check whether sampled site satisfies constraints accepted = satisfies_spatial_constraints(site.r) && @@ -297,7 +297,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) } } -SourceSite IndependentSource::sample(uint64_t* seed) const +SourceSite IndependentSource::sample_with_constraints(uint64_t* seed) const { SourceSite site; site.particle = particle_; @@ -419,7 +419,7 @@ void FileSource::load_sites_from_file(const std::string& path) file_close(file_id); } -SourceSite FileSource::sample_without_constraints(uint64_t* seed) const +SourceSite FileSource::sample(uint64_t* seed) const { // Sample a particle randomly from list size_t i_site = sites_.size() * prn(seed); @@ -519,7 +519,7 @@ MeshSource::MeshSource(pugi::xml_node node) : Source(node) space_ = std::make_unique(mesh_idx, strengths); } -SourceSite MeshSource::sample(uint64_t* seed) const +SourceSite MeshSource::sample_with_constraints(uint64_t* seed) const { // Sample the CDF defined in initialization above int32_t element = space_->sample_element_index(seed); @@ -533,7 +533,7 @@ SourceSite MeshSource::sample(uint64_t* seed) const SourceSite site; while (true) { // Sample source for the chosen element and replace the position - site = source(element)->sample(seed); + site = source(element)->sample_with_constraints(seed); site.r = r; // Apply other rejections @@ -596,7 +596,7 @@ SourceSite sample_external_source(uint64_t* seed) } // Sample source site from i-th source distribution - SourceSite site {model::external_sources[i]->sample(seed)}; + SourceSite site {model::external_sources[i]->sample_with_constraints(seed)}; // If running in MG, convert site.E to group if (!settings::run_CE) { From a60ff904b00cc482844186899a29f8ab7c2fbe8c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 May 2024 12:39:35 -0500 Subject: [PATCH 49/57] Make sample_with_constraints non-virtual --- include/openmc/source.h | 29 +++++++++++++++++++---------- src/source.cpp | 22 +++++++++++++--------- 2 files changed, 32 insertions(+), 19 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 98e8de94fe3..c8fa06f4503 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -56,7 +56,7 @@ class Source { // //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - virtual SourceSite sample_with_constraints(uint64_t* seed) const; + SourceSite sample_with_constraints(uint64_t* seed) const; //! Sample a source site (without applying constraints) // @@ -67,13 +67,6 @@ class Source { static unique_ptr create(pugi::xml_node node); -protected: - // Methods for constraints - void read_constraints(pugi::xml_node node); - bool satisfies_spatial_constraints(Position r) const; - bool satisfies_energy_constraints(double E) const; - bool satisfies_time_constraints(double time) const; - protected: // Domain types enum class DomainType { UNIVERSE, MATERIAL, CELL }; @@ -84,6 +77,15 @@ class Source { // until constraints are met. enum class RejectionStrategy { KILL, RESAMPLE }; + // Indicates whether derived class already handles constraints + virtual bool constraints_applied() const { return false; } + + // Methods for constraints + void read_constraints(pugi::xml_node node); + bool satisfies_spatial_constraints(Position r) const; + bool satisfies_energy_constraints(double E) const; + bool satisfies_time_constraints(double time) const; + // Data members double strength_ {1.0}; //!< Source strength std::unordered_set domain_ids_; //!< Domains to reject from @@ -113,7 +115,7 @@ class IndependentSource : public Source { //! Sample from the external source distribution //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - SourceSite sample_with_constraints(uint64_t* seed) const override; + SourceSite sample(uint64_t* seed) const override; // Properties ParticleType particle_type() const { return particle_; } @@ -124,6 +126,10 @@ class IndependentSource : public Source { Distribution* energy() const { return energy_.get(); } Distribution* time() const { return time_.get(); } +protected: + // Indicates whether derived class already handles constraints + bool constraints_applied() const override { return true; } + private: // Data members ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted @@ -194,7 +200,7 @@ class MeshSource : public Source { //! Sample from the external source distribution //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - SourceSite sample_with_constraints(uint64_t* seed) const override; + SourceSite sample(uint64_t* seed) const override; // Properties double strength() const override { return space_->total_strength(); } @@ -205,6 +211,9 @@ class MeshSource : public Source { return sources_.size() == 1 ? sources_[0] : sources_[i]; } +protected: + bool constraints_applied() const override { return true; } + private: // Data members unique_ptr space_; //!< Mesh spatial diff --git a/src/source.cpp b/src/source.cpp index 914e584e8e2..df2673251b3 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -155,14 +155,18 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const // Sample a source site without considering constraints yet site = this->sample(seed); - // Check whether sampled site satisfies constraints - accepted = satisfies_spatial_constraints(site.r) && - satisfies_energy_constraints(site.E) && - satisfies_time_constraints(site.time); - if (!accepted && rejection_strategy_ == RejectionStrategy::KILL) { - // Accept particle but set weight to 0 so that it is killed immediately + if (constraints_applied()) { accepted = true; - site.wgt = 0.0; + } else { + // Check whether sampled site satisfies constraints + accepted = satisfies_spatial_constraints(site.r) && + satisfies_energy_constraints(site.E) && + satisfies_time_constraints(site.time); + if (!accepted && rejection_strategy_ == RejectionStrategy::KILL) { + // Accept particle but set weight to 0 so that it is killed immediately + accepted = true; + site.wgt = 0.0; + } } } return site; @@ -297,7 +301,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) } } -SourceSite IndependentSource::sample_with_constraints(uint64_t* seed) const +SourceSite IndependentSource::sample(uint64_t* seed) const { SourceSite site; site.particle = particle_; @@ -519,7 +523,7 @@ MeshSource::MeshSource(pugi::xml_node node) : Source(node) space_ = std::make_unique(mesh_idx, strengths); } -SourceSite MeshSource::sample_with_constraints(uint64_t* seed) const +SourceSite MeshSource::sample(uint64_t* seed) const { // Sample the CDF defined in initialization above int32_t element = space_->sample_element_index(seed); From e3f56e7652ceb09a26b20678cdbfea4ea83fdc73 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 May 2024 07:02:05 -0500 Subject: [PATCH 50/57] Add user's guide section on source constraints --- docs/source/usersguide/settings.rst | 46 +++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index e966423f0e3..19c120c853a 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -420,6 +420,52 @@ string, which gets passed down to the ``openmc_create_source()`` function:: settings.source = openmc.CompiledSource('libsource.so', '3.5e6') +Source Constraints +------------------ + +All source classes in OpenMC have the ability to apply a set of "constraints" +that limit which sampled source sites are actually used for transport. The most +common use case is to sample source sites over some simple spatial distribution +(e.g., uniform over a box) and then only accept those that appear in a given +cell or material. This can be done with a domain constraint, which can be +specified as follows:: + + source_cell = openmc.Cell(...) + ... + + spatial_dist = openmc.stats.Box((-10., -10., -10.), (10., 10., 10.)) + source = openmc.IndependentSource( + space=spatial_dist, + constraints={'domains': [source_cell]} + ) + +For k-eigenvalue problems, a convenient constraint is available that limits +source sites to those sampled in a fissionable material:: + + source = openmc.IndependentSource( + space=spatial_dist, constraints={'fissionable': True} + ) + +Constraints can also be placed on a range of energies or times:: + + # Only use source sites between 500 keV and 1 MeV and with times under 1 sec + source = openmc.FileSource( + 'source.h5', + constraints={'energy_bounds': [500.0e3, 1.0e6], 'time_bounds': [0.0, 1.0]} + ) + +Normally, when a source site is rejected, a new one will be resampled until one +is found that meets the constraints. However, the rejection strategy can be +changed so that a rejected site will just not be simulated by specifying:: + + source = openmc.IndependentSource( + space=spatial_dist, + constraints={'domains': [cell], 'rejection_strategy': 'kill'} + ) + +In this case, the actual number of particles simulated may be less than what you +specified in :attr:`Settings.particles`. + .. _usersguide_entropy: --------------- From e468dd4ee05cbec65b6225e2abfd3d50b96bb949 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 May 2024 11:58:34 -0500 Subject: [PATCH 51/57] Put source constraints into a XML element --- docs/source/io_formats/settings.rst | 38 +++++++++++++++++++++++++-- docs/source/usersguide/settings.rst | 2 ++ openmc/source.py | 40 +++++++++++++++++------------ src/source.cpp | 8 ++++++ 4 files changed, 69 insertions(+), 19 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6318a991601..ae3266dabed 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -534,8 +534,9 @@ attributes/sub-elements: *Default*: 1.0 :type: - Indicator of source type. One of ``independent``, ``file``, ``compiled``, or ``mesh``. - The type of the source will be determined by this attribute if it is present. + Indicator of source type. One of ``independent``, ``file``, ``compiled``, or + ``mesh``. The type of the source will be determined by this attribute if it + is present. :particle: The source particle type, either ``neutron`` or ``photon``. @@ -716,6 +717,39 @@ attributes/sub-elements: mesh element and follows the format for :ref:`source_element`. The number of ```` sub-elements should correspond to the number of mesh elements. + :constraints: + This sub-element indicates the presence of constraints on sampled source + sites (see :ref:`usersguide_source_constraints` for details). It may have + the following sub-elements: + + :domain_ids: + The unique IDs of domains for which source sites must be within. + + *Default*: None + + :domain_type: + The type of each domain for source rejection ("cell", "material", or + "universe"). + + *Default*: None + + :fissionable: + A boolean indicating whether source sites must be sampled within a + material that is fissionable in order to be accepted. + + :time_bounds: + A pair of times in [s] indicating the lower and upper bound for a time + interval that source particles must be within. + + :energy_bounds: + A pair of energies in [eV] indicating the lower and upper bound for an + energy interval that source particles must be within. + + :rejection_strategy: + Either "resample", indicating that source sites should be resampled when + one is rejected, or "kill", indicating that a rejected source site is + assigned zero weight. + .. _univariate: Univariate Probability Distributions diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 19c120c853a..eb02f654c30 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -420,6 +420,8 @@ string, which gets passed down to the ``openmc_create_source()`` function:: settings.source = openmc.CompiledSource('libsource.so', '3.5e6') +.. _usersguide_source_constraints: + Source Constraints ------------------ diff --git a/openmc/source.py b/openmc/source.py index 6380cadd56e..91318f8e9ad 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -136,23 +136,25 @@ def to_xml_element(self) -> ET.Element: element.set("strength", str(self.strength)) self.populate_xml_element(element) constraints = self.constraints - if "domain_ids" in constraints: - dt_elem = ET.SubElement(element, "domain_type") - dt_elem.text = constraints["domain_type"] - id_elem = ET.SubElement(element, "domain_ids") - id_elem.text = ' '.join(str(uid) for uid in constraints["domain_ids"]) - if "time_bounds" in constraints: - dt_elem = ET.SubElement(element, "time_bounds") - dt_elem.text = ' '.join(str(t) for t in constraints["time_bounds"]) - if "energy_bounds" in constraints: - dt_elem = ET.SubElement(element, "energy_bounds") - dt_elem.text = ' '.join(str(E) for E in constraints["energy_bounds"]) - if "fissionable" in constraints: - dt_elem = ET.SubElement(element, "fissionable") - dt_elem.text = str(constraints["fissionable"]).lower() - if "rejection_strategy" in constraints: - dt_elem = ET.SubElement(element, "rejection_strategy") - dt_elem.text = constraints["rejection_strategy"] + if constraints: + constraints_elem = ET.SubElement(element, "constraints") + if "domain_ids" in constraints: + dt_elem = ET.SubElement(constraints_elem, "domain_type") + dt_elem.text = constraints["domain_type"] + id_elem = ET.SubElement(constraints_elem, "domain_ids") + id_elem.text = ' '.join(str(uid) for uid in constraints["domain_ids"]) + if "time_bounds" in constraints: + dt_elem = ET.SubElement(constraints_elem, "time_bounds") + dt_elem.text = ' '.join(str(t) for t in constraints["time_bounds"]) + if "energy_bounds" in constraints: + dt_elem = ET.SubElement(constraints_elem, "energy_bounds") + dt_elem.text = ' '.join(str(E) for E in constraints["energy_bounds"]) + if "fissionable" in constraints: + dt_elem = ET.SubElement(constraints_elem, "fissionable") + dt_elem.text = str(constraints["fissionable"]).lower() + if "rejection_strategy" in constraints: + dt_elem = ET.SubElement(constraints_elem, "rejection_strategy") + dt_elem.text = constraints["rejection_strategy"] return element @@ -199,6 +201,10 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: @staticmethod def _get_constraints(elem: ET.Element) -> Dict[str, Any]: + # Find element containing constraints + constraints_elem = elem.find("constraints") + elem = constraints_elem if constraints_elem is not None else elem + constraints = {} domain_type = get_text(elem, "domain_type") if domain_type is not None: diff --git a/src/source.cpp b/src/source.cpp index df2673251b3..cb7a545c973 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -95,6 +95,14 @@ unique_ptr Source::create(pugi::xml_node node) void Source::read_constraints(pugi::xml_node node) { + // Check for constraints node. For backwards compatibility, if no constraints + // node is given, still try searching for domain constraints from top-level + // node. + pugi::xml_node constraints_node = node.child("constraints"); + if (constraints_node) { + node = constraints_node; + } + // Check for domains to reject from if (check_for_node(node, "domain_type")) { std::string domain_type = get_node_value(node, "domain_type"); From 42d5cadf0b6ce951fa4da4c7427862e79b767b5c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 May 2024 12:26:31 -0500 Subject: [PATCH 52/57] Add comment describing Source base/derived relationship --- include/openmc/source.h | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/include/openmc/source.h b/include/openmc/source.h index c8fa06f4503..42975c8ae93 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -40,6 +40,14 @@ extern vector> external_sources; //============================================================================== //! Abstract source interface +// +//! The Source class provides the interface that must be implemented by derived +//! classes, namely the sample() method that returns a sampled source site. From +//! this base class, source rejection is handled within the +//! sample_with_constraints() method. However, note that some classes directly +//! check for constraints for efficiency reasons (like IndependentSource), in +//! which case the constraints_applied() method indicates that constraints +//! should not be checked a second time from the base class. //============================================================================== class Source { From c9d130d54a7904624f538f0b5fc1650dee4bf86d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 May 2024 21:20:18 -0500 Subject: [PATCH 53/57] Update reference inputs to reflect contraints format --- .../asymmetric_lattice/inputs_true.dat | 4 +- .../mg_temperature_multi/inputs_true.dat | 4 +- .../mgxs_library_ce_to_mg/inputs_true.dat | 5 +- .../inputs_true.dat | 5 +- .../mgxs_library_condense/inputs_true.dat | 5 +- .../mgxs_library_correction/inputs_true.dat | 5 +- .../mgxs_library_distribcell/inputs_true.dat | 1047 +++++++++-------- .../mgxs_library_hdf5/inputs_true.dat | 5 +- .../mgxs_library_histogram/inputs_true.dat | 5 +- .../mgxs_library_no_nuclides/inputs_true.dat | 5 +- .../mgxs_library_nuclides/inputs_true.dat | 5 +- .../inputs_true.dat | 5 +- .../surface_tally/inputs_true.dat | 4 +- 13 files changed, 569 insertions(+), 535 deletions(-) diff --git a/tests/regression_tests/asymmetric_lattice/inputs_true.dat b/tests/regression_tests/asymmetric_lattice/inputs_true.dat index d3d9ea3336e..302ef24c200 100644 --- a/tests/regression_tests/asymmetric_lattice/inputs_true.dat +++ b/tests/regression_tests/asymmetric_lattice/inputs_true.dat @@ -212,7 +212,9 @@ -32 -32 0 32 32 32 - true + + true + diff --git a/tests/regression_tests/mg_temperature_multi/inputs_true.dat b/tests/regression_tests/mg_temperature_multi/inputs_true.dat index a3879251d4c..0c452cd520e 100644 --- a/tests/regression_tests/mg_temperature_multi/inputs_true.dat +++ b/tests/regression_tests/mg_temperature_multi/inputs_true.dat @@ -31,7 +31,9 @@ -0.63 -0.63 -1 0.63 0.63 1 - true + + true + multi-group diff --git a/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat b/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat index d9a3c52c840..6c071263f3d 100644 --- a/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat index f84adc32967..ad3d50974b0 100644 --- a/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_condense/inputs_true.dat b/tests/regression_tests/mgxs_library_condense/inputs_true.dat index 93519cd2a2a..8606e9fdda3 100644 --- a/tests/regression_tests/mgxs_library_condense/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_condense/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_correction/inputs_true.dat b/tests/regression_tests/mgxs_library_correction/inputs_true.dat index 80cbc4bb56a..ea2d1b714a2 100644 --- a/tests/regression_tests/mgxs_library_correction/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_correction/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat b/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat index fc2c06d1cfa..bd6abbb7cb5 100644 --- a/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat @@ -1,43 +1,17 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1.26 1.26 - 17 17 - -10.71 -10.71 - + + + + + + + + + + 1.26 1.26 + 17 17 + -10.71 -10.71 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 @@ -55,487 +29,514 @@ 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 - - - - - - - - - - eigenvalue - 100 - 10 - 5 - - - -10.71 -10.71 -1 10.71 10.71 1 - - - - - - 1 - - - 0.0 20000000.0 - - - 0.0 20000000.0 - - - 1 - - - 3 - - - 1 2 3 4 5 6 - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - nu-scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - absorption - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - absorption - tracklength - - - 1 2 - total - (n,2n) - tracklength - - - 1 2 - total - (n,3n) - tracklength - - - 1 2 - total - (n,4n) - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - absorption - tracklength - - - 1 2 - total - fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - nu-fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - kappa-fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - scatter - tracklength - - - 1 2 - total - flux - analog - - - 1 2 - total - nu-scatter - analog - - - 1 2 - total - flux - analog - - - 1 2 5 30 - total - scatter - analog - - - 1 2 - total - flux - analog - - - 1 2 5 30 - total - nu-scatter - analog - - - 1 2 5 - total - nu-scatter - analog - - - 1 2 5 - total - scatter - analog - - - 1 2 - total - flux - analog - - - 1 2 5 - total - nu-fission - analog - - - 1 2 5 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - scatter - tracklength - - - 1 2 5 30 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - scatter - tracklength - - - 1 2 5 30 - total - scatter - analog - - - 1 2 5 - total - nu-scatter - analog - - - 1 2 - total - nu-fission - analog - - - 1 5 - total - nu-fission - analog - - - 1 2 - total - prompt-nu-fission - analog - - - 1 5 - total - prompt-nu-fission - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - inverse-velocity - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - prompt-nu-fission - tracklength - - - 1 2 - total - flux - analog - - - 1 2 5 - total - prompt-nu-fission - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - nu-scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 75 2 - total - delayed-nu-fission - tracklength - - - 1 75 2 - total - delayed-nu-fission - analog - - - 1 75 5 - total - delayed-nu-fission - analog - - - 1 2 - total - nu-fission - tracklength - - - 1 75 2 - total - delayed-nu-fission - tracklength - - - 1 75 - total - delayed-nu-fission - tracklength - - - 1 75 - total - decay-rate - tracklength - - - 1 2 - total - flux - analog - - - 1 75 2 5 - total - delayed-nu-fission - analog - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + + + -10.71 -10.71 -1 10.71 10.71 1 + + + + + + + 1 + + + 0.0 20000000.0 + + + 0.0 20000000.0 + + + 1 + + + 3 + + + 1 2 3 4 5 6 + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + nu-scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + absorption + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + absorption + tracklength + + + 1 2 + total + (n,2n) + tracklength + + + 1 2 + total + (n,3n) + tracklength + + + 1 2 + total + (n,4n) + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + absorption + tracklength + + + 1 2 + total + fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + nu-fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + kappa-fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + scatter + tracklength + + + 1 2 + total + flux + analog + + + 1 2 + total + nu-scatter + analog + + + 1 2 + total + flux + analog + + + 1 2 5 30 + total + scatter + analog + + + 1 2 + total + flux + analog + + + 1 2 5 30 + total + nu-scatter + analog + + + 1 2 5 + total + nu-scatter + analog + + + 1 2 5 + total + scatter + analog + + + 1 2 + total + flux + analog + + + 1 2 5 + total + nu-fission + analog + + + 1 2 5 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + scatter + tracklength + + + 1 2 5 30 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + scatter + tracklength + + + 1 2 5 30 + total + scatter + analog + + + 1 2 5 + total + nu-scatter + analog + + + 1 2 + total + nu-fission + analog + + + 1 5 + total + nu-fission + analog + + + 1 2 + total + prompt-nu-fission + analog + + + 1 5 + total + prompt-nu-fission + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + inverse-velocity + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + prompt-nu-fission + tracklength + + + 1 2 + total + flux + analog + + + 1 2 5 + total + prompt-nu-fission + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + nu-scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 75 2 + total + delayed-nu-fission + tracklength + + + 1 75 2 + total + delayed-nu-fission + analog + + + 1 75 5 + total + delayed-nu-fission + analog + + + 1 2 + total + nu-fission + tracklength + + + 1 75 2 + total + delayed-nu-fission + tracklength + + + 1 75 + total + delayed-nu-fission + tracklength + + + 1 75 + total + decay-rate + tracklength + + + 1 2 + total + flux + analog + + + 1 75 2 5 + total + delayed-nu-fission + analog + + diff --git a/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat b/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat index 93519cd2a2a..8606e9fdda3 100644 --- a/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_histogram/inputs_true.dat b/tests/regression_tests/mgxs_library_histogram/inputs_true.dat index 178a44464e9..013fcd0fa5f 100644 --- a/tests/regression_tests/mgxs_library_histogram/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_histogram/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat index 0003413acb1..5cad36b7ea3 100644 --- a/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat index 97e227177e0..9b3a22510b1 100644 --- a/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat index b73ed179e5e..1372741c947 100644 --- a/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/surface_tally/inputs_true.dat b/tests/regression_tests/surface_tally/inputs_true.dat index a6b30aeb221..192c2e89f40 100644 --- a/tests/regression_tests/surface_tally/inputs_true.dat +++ b/tests/regression_tests/surface_tally/inputs_true.dat @@ -33,7 +33,9 @@ -0.62992 -0.62992 -1 0.62992 0.62992 1 - true + + true + From bad3a01f406037a0bba45f7e7f54455de157a65e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 May 2024 22:57:02 -0500 Subject: [PATCH 54/57] Fix one more reference input file --- .../mgxs_library_distribcell/inputs_true.dat | 1050 +++++++++-------- 1 file changed, 526 insertions(+), 524 deletions(-) diff --git a/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat b/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat index bd6abbb7cb5..f1a60533491 100644 --- a/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat @@ -1,17 +1,43 @@ - - - - - - - - - - 1.26 1.26 - 17 17 - -10.71 -10.71 - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.26 1.26 + 17 17 + -10.71 -10.71 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 @@ -29,514 +55,490 @@ 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - eigenvalue - 100 - 10 - 5 - - - -10.71 -10.71 -1 10.71 10.71 1 - - - - - - - 1 - - - 0.0 20000000.0 - - - 0.0 20000000.0 - - - 1 - - - 3 - - - 1 2 3 4 5 6 - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - nu-scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - absorption - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - absorption - tracklength - - - 1 2 - total - (n,2n) - tracklength - - - 1 2 - total - (n,3n) - tracklength - - - 1 2 - total - (n,4n) - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - absorption - tracklength - - - 1 2 - total - fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - nu-fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - kappa-fission - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - scatter - tracklength - - - 1 2 - total - flux - analog - - - 1 2 - total - nu-scatter - analog - - - 1 2 - total - flux - analog - - - 1 2 5 30 - total - scatter - analog - - - 1 2 - total - flux - analog - - - 1 2 5 30 - total - nu-scatter - analog - - - 1 2 5 - total - nu-scatter - analog - - - 1 2 5 - total - scatter - analog - - - 1 2 - total - flux - analog - - - 1 2 5 - total - nu-fission - analog - - - 1 2 5 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - scatter - tracklength - - - 1 2 5 30 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - scatter - tracklength - - - 1 2 5 30 - total - scatter - analog - - - 1 2 5 - total - nu-scatter - analog - - - 1 2 - total - nu-fission - analog - - - 1 5 - total - nu-fission - analog - - - 1 2 - total - prompt-nu-fission - analog - - - 1 5 - total - prompt-nu-fission - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - inverse-velocity - tracklength - - - 1 2 - total - flux - tracklength - - - 1 2 - total - prompt-nu-fission - tracklength - - - 1 2 - total - flux - analog - - - 1 2 5 - total - prompt-nu-fission - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 2 - total - total - tracklength - - - 1 2 - total - flux - analog - - - 1 5 6 - total - nu-scatter - analog - - - 1 2 - total - flux - tracklength - - - 1 75 2 - total - delayed-nu-fission - tracklength - - - 1 75 2 - total - delayed-nu-fission - analog - - - 1 75 5 - total - delayed-nu-fission - analog - - - 1 2 - total - nu-fission - tracklength - - - 1 75 2 - total - delayed-nu-fission - tracklength - - - 1 75 - total - delayed-nu-fission - tracklength - - - 1 75 - total - decay-rate - tracklength - - - 1 2 - total - flux - analog - - - 1 75 2 5 - total - delayed-nu-fission - analog - - + + + + + + + + + + eigenvalue + 100 + 10 + 5 + + + -10.71 -10.71 -1 10.71 10.71 1 + + + true + + + + + + 1 + + + 0.0 20000000.0 + + + 0.0 20000000.0 + + + 1 + + + 3 + + + 1 2 3 4 5 6 + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + nu-scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + absorption + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + absorption + tracklength + + + 1 2 + total + (n,2n) + tracklength + + + 1 2 + total + (n,3n) + tracklength + + + 1 2 + total + (n,4n) + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + absorption + tracklength + + + 1 2 + total + fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + nu-fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + kappa-fission + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + scatter + tracklength + + + 1 2 + total + flux + analog + + + 1 2 + total + nu-scatter + analog + + + 1 2 + total + flux + analog + + + 1 2 5 30 + total + scatter + analog + + + 1 2 + total + flux + analog + + + 1 2 5 30 + total + nu-scatter + analog + + + 1 2 5 + total + nu-scatter + analog + + + 1 2 5 + total + scatter + analog + + + 1 2 + total + flux + analog + + + 1 2 5 + total + nu-fission + analog + + + 1 2 5 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + scatter + tracklength + + + 1 2 5 30 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + scatter + tracklength + + + 1 2 5 30 + total + scatter + analog + + + 1 2 5 + total + nu-scatter + analog + + + 1 2 + total + nu-fission + analog + + + 1 5 + total + nu-fission + analog + + + 1 2 + total + prompt-nu-fission + analog + + + 1 5 + total + prompt-nu-fission + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + inverse-velocity + tracklength + + + 1 2 + total + flux + tracklength + + + 1 2 + total + prompt-nu-fission + tracklength + + + 1 2 + total + flux + analog + + + 1 2 5 + total + prompt-nu-fission + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 2 + total + total + tracklength + + + 1 2 + total + flux + analog + + + 1 5 6 + total + nu-scatter + analog + + + 1 2 + total + flux + tracklength + + + 1 75 2 + total + delayed-nu-fission + tracklength + + + 1 75 2 + total + delayed-nu-fission + analog + + + 1 75 5 + total + delayed-nu-fission + analog + + + 1 2 + total + nu-fission + tracklength + + + 1 75 2 + total + delayed-nu-fission + tracklength + + + 1 75 + total + delayed-nu-fission + tracklength + + + 1 75 + total + decay-rate + tracklength + + + 1 2 + total + flux + analog + + + 1 75 2 5 + total + delayed-nu-fission + analog + + + From 42898cc60735d9793acfa90cfcd01dbd12b17b13 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 May 2024 23:24:39 -0500 Subject: [PATCH 55/57] Implement rejection fraction check in sample_with_constraints --- include/openmc/source.h | 2 +- src/source.cpp | 23 ++++++++++++++++++----- tests/unit_tests/test_source.py | 14 ++++++++++++++ 3 files changed, 33 insertions(+), 6 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 42975c8ae93..9ef8b925127 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -71,7 +71,7 @@ class Source { //! Sample from the external source distribution //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - virtual SourceSite sample(uint64_t* seed) const { return {}; } + virtual SourceSite sample(uint64_t* seed) const = 0; static unique_ptr create(pugi::xml_node node); diff --git a/src/source.cpp b/src/source.cpp index cb7a545c973..453d4befa9b 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -157,6 +157,8 @@ void Source::read_constraints(pugi::xml_node node) SourceSite Source::sample_with_constraints(uint64_t* seed) const { bool accepted = false; + static int n_reject = 0; + static int n_accept = 0; SourceSite site; while (!accepted) { @@ -170,10 +172,21 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const accepted = satisfies_spatial_constraints(site.r) && satisfies_energy_constraints(site.E) && satisfies_time_constraints(site.time); - if (!accepted && rejection_strategy_ == RejectionStrategy::KILL) { - // Accept particle but set weight to 0 so that it is killed immediately - accepted = true; - site.wgt = 0.0; + if (!accepted) { + ++n_reject; + if (n_reject >= EXTSRC_REJECT_THRESHOLD && + static_cast(n_accept) / n_reject <= + EXTSRC_REJECT_FRACTION) { + fatal_error("More than 95% of external source sites sampled were " + "rejected. Please check your source definition."); + } + + // For the "kill" strategy, accept particle but set weight to 0 so that + // it is terminated immediately + if (rejection_strategy_ == RejectionStrategy::KILL) { + accepted = true; + site.wgt = 0.0; + } } } } @@ -316,7 +329,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Repeat sampling source location until a good site has been accepted bool accepted = false; - int n_reject = 0; + static int n_reject = 0; static int n_accept = 0; while (!accepted) { diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 729531c22c8..90c9683f67d 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -205,6 +205,20 @@ def test_constraints_file(sphere_box_model, run_in_tmpdir): openmc.lib.finalize() +def test_rejection_limit(sphere_box_model, run_in_tmpdir): + model, cell1 = sphere_box_model[:2] + + # Define a point source that will get rejected 100% of the time + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point((-3., 0., 0.)), + constraints={'domains': [cell1]} + ) + + # Confirm that OpenMC doesn't run in an infinite loop + with pytest.raises(RuntimeError, match="rejected"): + model.run() + + def test_exceptions(): with pytest.raises(AttributeError, match=r'Please use the FileSource class'): From ca3ba803dcbc82de6d8566be21485ba3576acb79 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 May 2024 11:37:30 -0500 Subject: [PATCH 56/57] Fix rejection threshold check in sample_with_constraints --- src/source.cpp | 4 ++++ tests/unit_tests/test_source.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/source.cpp b/src/source.cpp index 453d4befa9b..405de998c89 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -190,6 +190,10 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const } } } + + // Increment number of accepted samples + ++n_accept; + return site; } diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 90c9683f67d..2661910b448 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -127,7 +127,7 @@ def test_constraints_independent(sphere_box_model, run_in_tmpdir): model, cell1, cell2, cell3 = sphere_box_model # Set up a box source with rejection on the spherical cell - space = openmc.stats.Box(*cell3.bounding_box) + space = openmc.stats.Box((-4., -1., -1.), (4., 1., 1.)) model.settings.source = openmc.IndependentSource( space=space, constraints={'domains': [cell1, cell2]} ) From da08d77fda1d70cd2db6813fd68578ea2aa4fc21 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 May 2024 21:28:04 -0500 Subject: [PATCH 57/57] Skip rejection limit test for MPI config --- tests/unit_tests/test_source.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 2661910b448..9a19f6f24dd 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -7,6 +7,8 @@ import pytest from pytest import approx +from tests.regression_tests import config + def test_source(): space = openmc.stats.Point() @@ -205,6 +207,7 @@ def test_constraints_file(sphere_box_model, run_in_tmpdir): openmc.lib.finalize() +@pytest.mark.skipif(config['mpi'], reason='Not compatible with MPI') def test_rejection_limit(sphere_box_model, run_in_tmpdir): model, cell1 = sphere_box_model[:2] @@ -214,9 +217,11 @@ def test_rejection_limit(sphere_box_model, run_in_tmpdir): constraints={'domains': [cell1]} ) - # Confirm that OpenMC doesn't run in an infinite loop + # Confirm that OpenMC doesn't run in an infinite loop. Note that this may + # work when running with MPI since it won't necessarily capture the error + # message correctly with pytest.raises(RuntimeError, match="rejected"): - model.run() + model.run(openmc_exec=config['exe']) def test_exceptions():