From 9aa6134366c802a50646d66a8f03f9e84c09679d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 21 Jan 2022 08:47:10 -0600 Subject: [PATCH] Adding otf ww adjustment using lower bound ratio. --- include/openmc/particle_data.h | 7 ++++- include/openmc/weight_windows.h | 9 +++++++ openmc/weight_windows.py | 27 +++++++++++++++++++- src/weight_windows.cpp | 22 ++++++++++++++++ tests/regression_tests/weightwindows/test.py | 6 +++-- 5 files changed, 67 insertions(+), 4 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index d31c20022e..128685fef2 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -309,7 +309,9 @@ private: int n_event_ {0}; // number of events executed in this particle's history // Weight window information - int n_split_ {0}; // Number of splits this particle has undergone + int n_split_ {0}; // Number of times this particle has been split + double ww_factor_ { + 0.0}; // Particle-specific factor for on-the-fly weight window adjustment // DagMC state variables #ifdef DAGMC @@ -435,6 +437,9 @@ public: int n_split() const { return n_split_; } int& n_split() { return n_split_; } + double ww_factor() const { return ww_factor_; } + double& ww_factor() { return ww_factor_; } + #ifdef DAGMC moab::DagMC::RayHistory& history() { return history_; } Direction& last_dir() { return last_dir_; } diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 9c42138901..c2ef29ae8f 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -52,12 +52,20 @@ extern vector> weight_windows; struct WeightWindow { double lower_weight {-1}; // -1 indicates invalid state double upper_weight {1}; + double max_lb_ratio {1}; double survival_weight {0.5}; double weight_cutoff {DEFAULT_WEIGHT_CUTOFF}; int max_split {1}; //! Whether the weight window is in a valid state bool is_valid() const { return lower_weight >= 0.0; } + + //! Adjust the weight window by a constant factor + void scale(double factor) + { + lower_weight *= factor; + upper_weight *= factor; + } }; //============================================================================== @@ -96,6 +104,7 @@ private: vector lower_ww_; //!< Lower weight window bounds vector upper_ww_; //!< Upper weight window bounds double survival_ratio_ {3.0}; //!< Survival weight ratio + double max_lb_ratio_ {1.0}; //!< Maximum lower bound to particle weight ratio double weight_cutoff_ {DEFAULT_WEIGHT_CUTOFF}; //!< Weight cutoff int max_split_ {10}; //!< Maximum value for particle splitting int32_t mesh_idx_; //!< index in meshes vector diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index cb9fc51b81..df46e318c0 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -48,6 +48,10 @@ class WeightWindows(IDManagerMixin): survival_ratio : float Ratio of the survival weight to the lower weight window bound for rouletting + max_lower_bound_ratio: float + Maximum allowed ratio of a weight window's lower bound to particle weight. + A factor will be applied to raise the weight window to the particle's weight + during transport if exceeded. Only applied if > 1.0. (Default: 1.0). max_split : int Maximum allowable number of particles when splitting weight_cutoff : float @@ -76,6 +80,9 @@ class WeightWindows(IDManagerMixin): survival_ratio : float Ratio of the survival weight to the lower weight window bound for rouletting + max_lower_bound_ratio: float + Maximum allowed ratio of a weight window's lower bound to particle + weight. (Default: 1.0) max_split : int Maximum allowable number of particles when splitting weight_cutoff : float @@ -91,7 +98,7 @@ class WeightWindows(IDManagerMixin): def __init__(self, mesh, lower_ww_bounds, upper_ww_bounds=None, upper_bound_ratio=None, energy_bins=None, particle_type='neutron', - survival_ratio=3, max_split=10, weight_cutoff=1.e-38, id=None): + survival_ratio=3, max_lb_ratio=1.0, max_split=10, weight_cutoff=1.e-38, id=None): self.mesh = mesh self.id = id self.particle_type = particle_type @@ -121,6 +128,7 @@ class WeightWindows(IDManagerMixin): 'do not match') self.survival_ratio = survival_ratio + self.max_lower_bound_ratio = max_lb_ratio self.max_split = max_split self.weight_cutoff = weight_cutoff @@ -192,6 +200,16 @@ class WeightWindows(IDManagerMixin): cv.check_greater_than('Survival ratio', val, 1.0, True) self._survival_ratio = val + @property + def max_lower_bound_ratio(self): + return self._max_lower_bound_ratio + + @max_lower_bound_ratio.setter + def max_lower_bound_ratio(self, val): + cv.check_type('Maximum lower bound ratio', val, Real) + cv.check_greater_than('Maximum lower bound ratio', val, 1.0, True) + self._max_lower_bound_ratio = val + @property def max_split(self): return self._max_split @@ -241,6 +259,9 @@ class WeightWindows(IDManagerMixin): subelement = ET.SubElement(element, 'survival_ratio') subelement.text = str(self.survival_ratio) + subelement = ET.SubElement(element, 'max_lower_bound_ratio') + subelement.text = str(self.max_lower_bound_ratio) + subelement = ET.SubElement(element, 'max_split') subelement.text = str(self.max_split) @@ -278,6 +299,7 @@ class WeightWindows(IDManagerMixin): ebins = [float(b) for b in get_text(elem, 'energy_bins').split()] particle_type = get_text(elem, 'particle_type') survival_ratio = float(get_text(elem, 'survival_ratio')) + max_lb_ratio = float(get_text(elem, 'max_lower_bound_ratio')) max_split = int(get_text(elem, 'max_split')) weight_cutoff = float(get_text(elem, 'weight_cutoff')) id = int(get_text(elem, 'id')) @@ -289,6 +311,7 @@ class WeightWindows(IDManagerMixin): energy_bins=ebins, particle_type=particle_type, survival_ratio=survival_ratio, + max_lb_ratio=max_lb_ratio, max_split=max_split, weight_cutoff=weight_cutoff, id=id @@ -318,6 +341,7 @@ class WeightWindows(IDManagerMixin): lower_ww_bounds = group['lower_ww_bounds'][()] upper_ww_bounds = group['upper_ww_bounds'][()] survival_ratio = group['survival_ratio'][()] + max_lb_ratio = group['max_lower_bound_ratio'][()] max_split = group['max_split'][()] weight_cutoff = group['weight_cutoff'][()] @@ -328,6 +352,7 @@ class WeightWindows(IDManagerMixin): energy_bins=ebins, particle_type=ptype, survival_ratio=survival_ratio, + max_lb_ratio=max_lb_ratio, max_split=max_split, weight_cutoff=weight_cutoff, id=id diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index c3859d7d15..32886436a8 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -56,6 +56,19 @@ void apply_weight_windows(Particle& p) return; } + // check if particle is far above current weight window + // only do this if the factor is not already set on the particle and a + // maximum lower bound ratio is specified + if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 && + p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) { + p.ww_factor() = + p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio); + } + + // move weight window closer to the particle weight if needed + if (p.ww_factor() > 1.0) + weight_window.scale(p.ww_factor()); + // if particle's weight is above the weight window split until they are within // the window if (weight > weight_window.upper_weight) { @@ -139,6 +152,14 @@ WeightWindows::WeightWindows(pugi::xml_node node) "and less than the upper to lower weight window ratio."); } + // get the max lower bound ratio - optional + if (check_for_node(node, "max_lower_bound_ratio")) { + max_lb_ratio_ = std::stod(get_node_value(node, "max_lower_bound_ratio")); + if (max_lb_ratio_ < 1.0) { + fatal_error("Maximum lower bound ratio must be larger than 1"); + } + } + // get the max split - optional if (check_for_node(node, "max_split")) { max_split_ = std::stod(get_node_value(node, "max_split")); @@ -254,6 +275,7 @@ void WeightWindows::to_hdf5(hid_t group) const write_dataset(ww_group, "lower_ww_bounds", lower_ww_); write_dataset(ww_group, "upper_ww_bounds", upper_ww_); write_dataset(ww_group, "survival_ratio", survival_ratio_); + write_dataset(ww_group, "max_lower_bound_ratio", max_lb_ratio_); write_dataset(ww_group, "max_split", max_split_); write_dataset(ww_group, "weight_cutoff", weight_cutoff_); write_dataset(ww_group, "mesh", this->mesh().id_); diff --git a/tests/regression_tests/weightwindows/test.py b/tests/regression_tests/weightwindows/test.py index 63d773cbcd..e9e64b0466 100644 --- a/tests/regression_tests/weightwindows/test.py +++ b/tests/regression_tests/weightwindows/test.py @@ -92,13 +92,15 @@ def model(): ww_n_lower_bnds, None, 10.0, - e_bnds) + e_bnds, + max_lb_ratio=5.0) ww_p = openmc.WeightWindows(ww_mesh, ww_p_lower_bnds, None, 10.0, - e_bnds) + e_bnds, + max_lb_ratio=5.0) model.settings.weight_windows = [ww_n, ww_p]