diff --git a/CMakeLists.txt b/CMakeLists.txt index 9516971b9..cfffbdccc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -361,6 +361,7 @@ list(APPEND libopenmc_SOURCES src/track_output.cpp src/urr.cpp src/volume_calc.cpp + src/weight_windows.cpp src/wmp.cpp src/xml_interface.cpp src/xsdata.cpp) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 38fc2d79f..ead440ae9 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -620,9 +620,9 @@ variable and whose sub-elements/attributes are as follows: :type: The type of the distribution. Valid options are "uniform", "discrete", "tabular", "maxwell", "watt", and "mixture". The "uniform" option produces - variates sampled from a uniform distribution over a finite interval. The - "discrete" option produces random variates that can assume a finite number - of values (i.e., a distribution characterized by a probability mass function). + variates sampled from a uniform distribution over a finite interval. The + "discrete" option produces random variates that can assume a finite number + of values (i.e., a distribution characterized by a probability mass function). The "tabular" option produces random variates sampled from a tabulated distribution where the density function is either a histogram or linearly-interpolated between tabulated points. The "watt" option produces @@ -668,8 +668,8 @@ variable and whose sub-elements/attributes are as follows: For a "mixture" distribution, this element provides a distribution and its corresponding probability. :probability: - An attribute or ``pair`` that provides the probability of a univariate distribution within a "mixture" distribution. - + An attribute or ``pair`` that provides the probability of a univariate distribution within a "mixture" distribution. + :dist: This sub-element of a ``pair`` element provides information on the corresponding univariate distribution. @@ -1007,3 +1007,55 @@ sub-elements/attributes: sample points within. *Default*: None + +---------------------------- +```` Element +---------------------------- + +The ```` element specifies all necessary parameters for +mesh-based weight windows. This element has the following +sub-elements/attributes: + + :id: + A unique integer that is used to identify the weight windows + + :mesh: + ID of a mesh that is to be used for weight windows + + *Default*: None + + :particle_type: + The particle that the weight windows will apply to (e.g., 'neutron') + + *Default*: None + + :energy_bins: + Monotonically increasing list of bounding energies in [eV] to be used for + weight windows + + *Default*: None + + :lower_ww_bounds: + Lower weight window bound for each (energy bin, mesh bin) combination. + + *Default*: None + + :upper_ww_bounds: + Upper weight window bound for each (energy bin, mesh bin) combination. + + *Default*: None + + :survival: + The ratio of survival weight and lower weight window bound. + + *Default*: 3.0 + + :max_split: + Maximum allowable number of particles when splitting + + *Default*: 10 + + :weight_cutoff: + Threshold below which particles will be terminated + + *Default*: :math:`10^{-38}` diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst index 774575b46..7c1490647 100644 --- a/docs/source/methods/neutron_physics.rst +++ b/docs/source/methods/neutron_physics.rst @@ -1627,6 +1627,8 @@ cross sections. Variance Reduction Techniques ----------------------------- +.. _survival_biasing: + Survival Biasing ---------------- @@ -1677,6 +1679,47 @@ default, the cutoff weight in OpenMC is :math:`w_c = 0.25` and the survival weight is :math:`w_s = 1.0`. These parameters vary from one Monte Carlo code to another. +Weight Windows +-------------- + +In fixed source problems, it can often be difficult to obtain sufficiently low +variance on tallies in regions that are far from the source. The `weight window +method `_ was developed to increase the +population of particles in important spatial regions and energy ranges by +controlling particle weights. Each spatial region and particle energy range is +assigned upper and lower weight bounds, :math:`w_u` and :math:`w_\ell`, +respectively. When a particle is in a given spatial region / energy range, its +weight, :math:`w`, is compared to the lower and upper bounds. If the weight of +the particle is above the upper weight bound, the particle is split into +:math:`N` particles, where + +.. math:: + :label: ww-split + + N = \min(N_{max}, \lceil w/w_u \rceil) + +and :math:`N_{max}` is a user-defined maximum number of splits. To ensure a +fair game, each of the :math:`N` particles is assigned a weight :math:`w/N`. If +the weight is below :math:`w_\ell`, it is Russian rouletted as described in +:ref:`survival_biasing` with a survival weight :math:`w_s` that is set equal to + +.. math:: + :label: ww-survival-weight + + w_s = \min(N_{max} w, f_s w_l) + +where :math:`f_s` is a user-defined survival weight ratio greater than one. + +On top of the standard weight window method described above, OpenMC implements +two additional checks intended to mitigate problems with long histories. First, +particles with a weight that falls below some very small cutoff (defaults to +:math:`10^{-38}`) are killed with no Russian rouletting. Additionally, the total +number of splits experienced by a particle is tracked and if it reaches some +maximum value, it is prohibited from splitting further. + +At present, OpenMC allows weight windows to be defined on all supported mesh +types. + .. only:: html .. rubric:: References diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 2dea91bb0..ff849b92c 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -24,6 +24,7 @@ Simulation Settings openmc.Source openmc.SourceParticle openmc.VolumeCalculation + openmc.WeightWindows openmc.Settings The following function can be used for generating a source file: diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 2596a0334..5a099f956 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -358,6 +358,10 @@ constexpr int CMFD_NOACCEL {-1}; enum class GeometryType { CSG, DAG }; +//============================================================================== +// Variance Reduction constants +constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff + } // namespace openmc #endif // OPENMC_CONSTANTS_H diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 06fccc092..7d0382385 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -305,6 +305,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 + // DagMC state variables #ifdef DAGMC moab::DagMC::RayHistory history_; @@ -426,6 +429,9 @@ public: double& collision_distance() { return collision_distance_; } int& n_event() { return n_event_; } + int n_split() const { return n_split_; } + int& n_split() { return n_split_; } + #ifdef DAGMC moab::DagMC::RayHistory& history() { return history_; } Direction& last_dir() { return last_dir_; } diff --git a/include/openmc/physics.h b/include/openmc/physics.h index b2e72e047..32a43e26d 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -89,6 +89,11 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p); void sample_secondary_photons(Particle& p, int i_nuclide); +//!Split or Roulette particles based their weight and the lower weight window +// bound. +//! \param[in] p, particle to be split or rouletted with the weight window. +void split_particle(Particle& p); + } // namespace openmc #endif // OPENMC_PHYSICS_H diff --git a/include/openmc/settings.h b/include/openmc/settings.h index f14f6c15f..50be80a50 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -55,6 +55,7 @@ extern "C" bool trigger_on; //!< tally triggers enabled? extern bool trigger_predict; //!< predict batches for triggers? extern bool ufs_on; //!< uniform fission site method on? extern bool urr_ptables_on; //!< use unresolved resonance prob. tables? +extern bool weight_windows_on; //!< are weight windows are enabled? extern bool write_all_tracks; //!< write track files for every particle? extern bool write_initial_source; //!< write out initial source file? @@ -99,6 +100,7 @@ extern std::unordered_set statepoint_batch; //!< Batches when state should be written extern std::unordered_set source_write_surf_id; //!< Surface ids where sources will be written +extern int max_splits; //!< maximum number of particle splits for weight windows extern int64_t max_surface_particles; //!< maximum number of particles to be //!< banked on surfaces per process extern TemperatureMethod diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h new file mode 100644 index 000000000..12113d528 --- /dev/null +++ b/include/openmc/weight_windows.h @@ -0,0 +1,99 @@ +#ifndef OPENMC_WEIGHT_WINDOWS_H +#define OPENMC_WEIGHT_WINDOWS_H + +#include +#include + +#include +#include + +#include "openmc/constants.h" +#include "openmc/memory.h" +#include "openmc/mesh.h" +#include "openmc/particle.h" +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +// Non-member functions +//============================================================================== + +//! Apply weight windows to a particle +//! \param[in] p Particle to apply weight windows to +void apply_weight_windows(Particle& p); + +//! Free memory associated with weight windows +void free_memory_weight_windows(); + +//============================================================================== +// Global variables +//============================================================================== + +class WeightWindows; + +namespace variance_reduction { + +extern std::unordered_map ww_map; +extern vector> weight_windows; + +} // namespace variance_reduction + +//============================================================================== +//! Individual weight window information +//============================================================================== + +struct WeightWindow { + double lower_weight {-1}; // -1 indicates invalid state + double upper_weight {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; } +}; + +//============================================================================== +//! Weight window settings +//============================================================================== + +class WeightWindows { +public: + // Constructors + WeightWindows(); + WeightWindows(pugi::xml_node node); + + // Methods + + //! Set the weight window ID + void set_id(int32_t id = -1); + + // NOTE: This is unused for now but may be used in the future + //! Write weight window settings to an HDF5 file + //! \param[in] group HDF5 group to write to + void to_hdf5(hid_t group) const; + + //! Retrieve the weight window for a particle + //! \param[in] p Particle to get weight window for + WeightWindow get_weight_window(const Particle& p) const; + + // Accessors + int32_t id() const { return id_; } + const Mesh& mesh() const { return *model::meshes[mesh_idx_]; } + +private: + // Data members + int32_t id_; //!< Unique ID + ParticleType particle_type_; //!< Particle type to apply weight windows to + vector energy_bins_; //!< Energy bins [eV] + vector lower_ww_; //!< Lower weight window bounds + vector upper_ww_; //!< Upper weight window bounds + double survival_ratio_ {3.0}; //!< Survival 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 +}; + +} // namespace openmc +#endif // OPENMC_WEIGHT_WINDOWS_H diff --git a/openmc/__init__.py b/openmc/__init__.py index d118e6a05..06a06aacc 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -11,6 +11,7 @@ from openmc.plots import * from openmc.region import * from openmc.volume import * from openmc.source import * +from openmc.weight_windows import * from openmc.settings import * from openmc.surface import * from openmc.universe import * diff --git a/openmc/mesh.py b/openmc/mesh.py index 3cd208e89..2a16be1ad 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -88,6 +88,32 @@ class MeshBase(IDManagerMixin, ABC): else: raise ValueError('Unrecognized mesh type: "' + mesh_type + '"') + @classmethod + def from_xml(cls, elem): + """Generates a mesh from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.MeshBase + an openmc mesh object + + """ + mesh_type = get_text(elem, 'type') + + if mesh_type == 'regular' or mesh_type is None: + return RegularMesh.from_xml_element(elem) + elif mesh_type == 'rectilinear': + return RectilinearMesh.from_xml_element(elem) + elif mesh_type == 'unstructured': + return UnstructuredMesh.from_xml_element(elem) + else: + raise ValueError(f'Unrecognized mesh type "{mesh_type}" found.') + class RegularMesh(MeshBase): """A regular Cartesian mesh in one, two, or three dimensions @@ -305,7 +331,6 @@ class RegularMesh(MeshBase): if self._upper_right is not None: subelement = ET.SubElement(element, "upper_right") subelement.text = ' '.join(map(str, self._upper_right)) - if self._width is not None: subelement = ET.SubElement(element, "width") subelement.text = ' '.join(map(str, self._width)) @@ -574,19 +599,19 @@ class RectilinearMesh(MeshBase): fmt = '{0: <16}{1}{2}\n' string = super().__repr__() string += fmt.format('\tDimensions', '=\t', self.n_dimension) - x_grid_str = str(self._x_grid) if not self._x_grid else len(self._x_grid) + x_grid_str = str(self._x_grid) if self._x_grid is None else len(self._x_grid) string += fmt.format('\tN X pnts:', '=\t', x_grid_str) - if self._x_grid: + if self._x_grid is not None: string += fmt.format('\tX Min:', '=\t', self._x_grid[0]) string += fmt.format('\tX Max:', '=\t', self._x_grid[-1]) - y_grid_str = str(self._y_grid) if not self._y_grid else len(self._y_grid) + y_grid_str = str(self._y_grid) if self._y_grid is None else len(self._y_grid) string += fmt.format('\tN Y pnts:', '=\t', y_grid_str) - if self._y_grid: + if self._y_grid is not None: string += fmt.format('\tY Min:', '=\t', self._y_grid[0]) string += fmt.format('\tY Max:', '=\t', self._y_grid[-1]) - z_grid_str = str(self._z_grid) if not self._z_grid else len(self._z_grid) + z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid) string += fmt.format('\tN Z pnts:', '=\t', z_grid_str) - if self._z_grid: + if self._z_grid is not None: string += fmt.format('\tZ Min:', '=\t', self._z_grid[0]) string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1]) return string @@ -603,6 +628,29 @@ class RectilinearMesh(MeshBase): return mesh + @classmethod + def from_xml_element(cls, elem): + """Generate a rectilinear mesh from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.RectilinearMesh + Rectilinear mesh object + + """ + id = int(get_text(elem, 'id')) + mesh = cls(id) + mesh.x_grid = [float(x) for x in get_text(elem, 'x_grid').split()] + mesh.y_grid = [float(y) for y in get_text(elem, 'y_grid').split()] + mesh.z_grid = [float(z) for z in get_text(elem, 'z_grid').split()] + + return mesh + def to_xml_element(self): """Return XML representation of the mesh @@ -1243,4 +1291,4 @@ class UnstructuredMesh(MeshBase): library = get_text(elem, 'library') length_multiplier = float(get_text(elem, 'length_multiplier', 1.0)) - return cls(filename, library, mesh_id, '', length_multiplier) + return cls(filename, library, mesh_id, '', length_multiplier) \ No newline at end of file diff --git a/openmc/settings.py b/openmc/settings.py index 6495acdd6..1ae159702 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -6,7 +6,7 @@ from xml.etree import ElementTree as ET from math import ceil import openmc.checkvalue as cv -from . import VolumeCalculation, Source, RegularMesh +from . import VolumeCalculation, Source, RegularMesh, WeightWindows from ._xml import clean_indentation, get_text, reorder_attributes @@ -96,6 +96,10 @@ class Settings: .. versionadded:: 0.12 max_order : None or int Maximum scattering order to apply globally when in multi-group mode. + max_splits : int + Maximum number of times a particle can split during a history + + .. versionadded:: 0.13 no_reduce : bool Indicate that all user-defined and global tallies should not be reduced across processes in a parallel calculation. @@ -197,9 +201,16 @@ class Settings: described in :ref:`verbosity`. volume_calculations : VolumeCalculation or iterable of VolumeCalculation Stochastic volume calculation specifications + weight_windows : WeightWindows iterable of WeightWindows + Weight windows to use for variance reduction + + .. versionadded:: 0.13 + weight_windows_on : bool + Whether weight windows are enabled + + .. versionadded:: 0.13 write_initial_source : bool Indicate whether to write the initial source distribution to file - """ def __init__(self): @@ -272,6 +283,9 @@ class Settings: self._event_based = None self._max_particles_in_flight = None self._write_initial_source = None + self._weight_windows = cv.CheckedList(WeightWindows, 'weight windows') + self._weight_windows_on = None + self._max_splits = None @property def run_mode(self): @@ -445,6 +459,18 @@ class Settings: def write_initial_source(self): return self._write_initial_source + @property + def weight_windows(self): + return self._weight_windows + + @property + def weight_windows_on(self): + return self._weight_windows_on + + @property + def max_splits(self): + return self._max_splits + @run_mode.setter def run_mode(self, run_mode): cv.check_value('run mode', run_mode, {x.value for x in RunMode}) @@ -835,6 +861,23 @@ class Settings: cv.check_type('write initial source', value, bool) self._write_initial_source = value + @weight_windows.setter + def weight_windows(self, value): + if not isinstance(value, MutableSequence): + value = [value] + self._weight_windows = cv.CheckedList(WeightWindows, 'weight windows', value) + + @weight_windows_on.setter + def weight_windows_on(self, value): + cv.check_type('weight windows on', value, bool) + self._weight_windows_on = value + + @max_splits.setter + def max_splits(self, value): + cv.check_type('maximum particle splits', value, Integral) + cv.check_greater_than('max particles in flight', value, 0) + self._max_splits = value + def _create_run_mode_subelement(self, root): elem = ET.SubElement(root, "run_mode") elem.text = self._run_mode.value @@ -1128,6 +1171,25 @@ class Settings: elem = ET.SubElement(root, "write_initial_source") elem.text = str(self._write_initial_source).lower() + def _create_weight_windows_subelement(self, root): + for ww in self._weight_windows: + # Add weight window information + root.append(ww.to_xml_element()) + + # See if a element already exists -- if not, add it + path = f"./mesh[@id='{ww.mesh.id}']" + if root.find(path) is None: + root.append(ww.mesh.to_xml_element()) + + if self._weight_windows_on is not None: + elem = ET.SubElement(root, "weight_windows_on") + elem.text = str(self._weight_windows_on).lower() + + def _create_max_splits_subelement(self, root): + if self._max_splits is not None: + elem = ET.SubElement(root, "max_splits") + elem.text = str(self._max_splits) + def _eigenvalue_from_xml_element(self, root): elem = root.find('eigenvalue') if elem is not None: @@ -1410,6 +1472,20 @@ class Settings: if text is not None: self.write_initial_source = text in ('true', '1') + def _weight_windows_from_xml_element(self, root): + for elem in root.findall('weight_windows'): + ww = WeightWindows.from_xml_element(elem, root) + self.weight_windows.append(ww) + + text = get_text(root, 'weight_windows_on') + if text is not None: + self.weight_windows_on = text in ('true', '1') + + def _max_splits_from_xml_element(self, root): + text = get_text(root, 'max_splits') + if text is not None: + self.max_splits = int(text) + def export_to_xml(self, path='settings.xml'): """Export simulation settings to an XML file. @@ -1464,6 +1540,8 @@ class Settings: self._create_material_cell_offsets_subelement(root_element) self._create_log_grid_bins_subelement(root_element) self._create_write_initial_source_subelement(root_element) + self._create_weight_windows_subelement(root_element) + self._create_max_splits_subelement(root_element) # Clean the indentation in the file to be user-readable clean_indentation(root_element) @@ -1538,6 +1616,8 @@ class Settings: settings._material_cell_offsets_from_xml_element(root) settings._log_grid_bins_from_xml_element(root) settings._write_initial_source_from_xml_element(root) + settings._weight_windows_from_xml_element(root) + settings._max_splits_from_xml_element(root) # TODO: Get volume calculations diff --git a/openmc/summary.py b/openmc/summary.py index 6e6412061..41bae7f26 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -12,7 +12,7 @@ _VERSION_SUMMARY = 6 class Summary: - """Summary of geometry, materials, and tallies used in a simulation. + """Summary of model used in a simulation. Attributes ---------- diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py new file mode 100644 index 000000000..ee8a4a2bf --- /dev/null +++ b/openmc/weight_windows.py @@ -0,0 +1,334 @@ +from collections.abc import Iterable +from numbers import Real, Integral + +from xml.etree import ElementTree as ET +import numpy as np + +from openmc.filter import _PARTICLES +from openmc.mesh import MeshBase +import openmc.checkvalue as cv + +from ._xml import get_text +from .mixin import IDManagerMixin + + +class WeightWindows(IDManagerMixin): + """Mesh-based weight windows + + This class enables you to specify weight window parameters that are used in + a simulation. Multiple sets of weight windows can be defined for different + meshes and different particles. An iterable of :class:`WeightWindows` + instances can be assigned to the :attr:`~openmc.Settings.weight_windows` + attribute, which is then exported to XML. + + Weight window lower/upper bounds are to be specified for each combination of + a mesh element and an energy bin. Thus the total number of bounds should be + equal to the product of the number of mesh bins and the number of energy + bins. + + .. versionadded:: 0.13 + + Parameters + ---------- + mesh : openmc.MeshBase + Mesh for the weight windows + lower_ww_bounds : Iterable of Real + A list of values for which each value is the lower bound of a weight + window + upper_ww_bounds : Iterable of Real + A list of values for which each value is the upper bound of a weight + window + upper_bound_ratio : float + Ratio of the lower to upper weight window bounds + energy_bins : Iterable of Real + A list of values for which each successive pair constitutes a range of + energies in [eV] for a single bin + particle_type : {'neutron', 'photon'} + Particle type the weight windows apply to + survival_ratio : float + Ratio of the survival weight to the lower weight window bound for + rouletting + max_split : int + Maximum allowable number of particles when splitting + weight_cutoff : float + Threshold below which particles will be terminated + id : int + Unique identifier for the weight window settings. If not + specified an identifier will automatically be assigned. + + Attributes + ---------- + id : int + Unique identifier for the weight window settings. + mesh : openmc.MeshBase + Mesh for the weight windows + particle_type : str + Particle type the weight windows apply to + energy_bins : Iterable of Real + A list of values for which each successive pair constitutes a range of + energies in [eV] for a single bin + lower_ww_bounds : Iterable of Real + A list of values for which each value is the lower bound of a weight + window + upper_ww_bounds : Iterable of Real + A list of values for which each value is the upper bound of a weight + window + survival_ratio : float + Ratio of the survival weight to the lower weight window bound for + rouletting + max_split : int + Maximum allowable number of particles when splitting + weight_cutoff : float + Threshold below which particles will be terminated + + See Also + -------- + openmc.Settings + + """ + next_id = 1 + used_ids = set() + + 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): + self.mesh = mesh + self.id = id + self.particle_type = particle_type + self.energy_bins = energy_bins + self.lower_ww_bounds = lower_ww_bounds + + cv.check_length('Lower window bounds', self.lower_ww_bounds, len(self.energy_bins)) + + if upper_ww_bounds is not None and upper_bound_ratio: + raise ValueError("Exactly one of uppwer_ww_bounds and " + "upper_bound_ratio must be present.") + + if upper_ww_bounds is None and upper_bound_ratio is None: + raise ValueError("Exactly one of uppwer_ww_bounds and " + "upper_bound_ratio must be present.") + + if upper_bound_ratio: + self.upper_ww_bounds = [ + lb * upper_bound_ratio for lb in self.lower_ww_bounds + ] + + if upper_ww_bounds is not None: + self.upper_ww_bounds = upper_ww_bounds + + if len(self.lower_ww_bounds) != len(self.upper_ww_bounds): + raise ValueError('Size of the lower and upper weight window bounds' + 'do not match') + + self.survival_ratio = survival_ratio + self.max_split = max_split + self.weight_cutoff = weight_cutoff + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tID', self._id) + string += '{: <16}=\t{}\n'.format('\tMesh:', self.mesh) + string += '{: <16}=\t{}\n'.format('\tParticle Type', self._particle_type) + string += '{: <16}=\t{}\n'.format('\tEnergy Bins', self._energy_bins) + string += '{: <16}=\t{}\n'.format('\tLower WW Bounds', self._lower_ww_bounds) + string += '{: <16}=\t{}\n'.format('\tUpper WW Bounds', self._upper_ww_bounds) + string += '{: <16}=\t{}\n'.format('\tSurvival Ratio', self._survival_ratio) + string += '{: <16}=\t{}\n'.format('\tMax Split', self._max_split) + string += '{: <16}=\t{}\n'.format('\tWeight Cutoff', self._weight_cutoff) + return string + + @property + def mesh(self): + return self._mesh + + @mesh.setter + def mesh(self, mesh): + cv.check_type('Weight window mesh', mesh, MeshBase) + self._mesh = mesh + + @property + def particle_type(self): + return self._particle_type + + @particle_type.setter + def particle_type(self, pt): + cv.check_value('Particle type', pt, _PARTICLES) + self._particle_type = pt + + @property + def energy_bins(self): + return self._energy_bins + + @energy_bins.setter + def energy_bins(self, bins): + cv.check_type('Energy bins', bins, Iterable, Real) + self._energy_bins = np.array(bins) + + @property + def lower_ww_bounds(self): + return self._lower_ww_bounds + + @lower_ww_bounds.setter + def lower_ww_bounds(self, bounds): + cv.check_type('Lower WW bounds', bounds, Iterable, Real) + self._lower_ww_bounds = np.array(bounds) + + @property + def upper_ww_bounds(self): + return self._upper_ww_bounds + + @upper_ww_bounds.setter + def upper_ww_bounds(self, bounds): + cv.check_type('Upper WW bounds', bounds, Iterable, Real) + self._upper_ww_bounds = np.array(bounds) + + @property + def survival_ratio(self): + return self._survival_ratio + + @survival_ratio.setter + def survival_ratio(self, val): + cv.check_type('Survival ratio', val, Real) + cv.check_greater_than('Survival ratio', val, 1.0, True) + self._survival_ratio = val + + @property + def max_split(self): + return self._max_split + + @max_split.setter + def max_split(self, val): + cv.check_type('Max split', val, Integral) + self._max_split = val + + @property + def weight_cutoff(self): + return self._weight_cutoff + + @weight_cutoff.setter + def weight_cutoff(self, cutoff): + cv.check_type('Weight cutoff', cutoff, Real) + cv.check_greater_than('Weight cutoff', cutoff, 0.0, True) + self._weight_cutoff = cutoff + + def to_xml_element(self): + """Return an XML representation of the weight window settings + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing the weight window information + """ + element = ET.Element('weight_windows') + + element.set('id', str(self._id)) + + subelement = ET.SubElement(element, 'mesh') + subelement.text = str(self.mesh.id) + + subelement = ET.SubElement(element, 'particle_type') + subelement.text = self.particle_type + + subelement = ET.SubElement(element, 'energy_bins') + subelement.text = ' '.join(str(e) for e in self.energy_bins) + + subelement = ET.SubElement(element, 'lower_ww_bounds') + subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds) + + subelement = ET.SubElement(element, 'upper_ww_bounds') + subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds) + + subelement = ET.SubElement(element, 'survival_ratio') + subelement.text = str(self.survival_ratio) + + subelement = ET.SubElement(element, 'max_split') + subelement.text = str(self.max_split) + + subelement = ET.SubElement(element, 'weight_cutoff') + subelement.text = str(self.weight_cutoff) + + return element + + @classmethod + def from_xml_element(cls, elem, root): + """Generate weight window settings from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + root : xml.etree.ElementTree.Element + Root element for the file where meshes can be found + + Returns + ------- + openmc.WeightWindows + Weight windows object + """ + # Get mesh for weight windows + mesh_id = int(get_text(elem, 'mesh')) + path = f"./mesh[@id='{mesh_id}']" + mesh_elem = root.find(path) + if mesh_elem is not None: + mesh = MeshBase.from_xml(mesh_elem) + + # Read all other parameters + lower_ww_bounds = [float(l) for l in get_text(elem, 'lower_ww_bounds').split()] + upper_ww_bounds = [float(u) for u in get_text(elem, 'upper_ww_bounds').split()] + 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_split = int(get_text(elem, 'max_split')) + weight_cutoff = float(get_text(elem, 'weight_cutoff')) + id = int(get_text(elem, 'id')) + + return cls( + mesh=mesh, + lower_ww_bounds=lower_ww_bounds, + upper_ww_bounds=upper_ww_bounds, + energy_bins=ebins, + particle_type=particle_type, + survival_ratio=survival_ratio, + max_split=max_split, + weight_cutoff=weight_cutoff, + id=id + ) + + @classmethod + def from_hdf5(cls, group, meshes): + """Create weight windows from HDF5 group + + Parameters + ---------- + group : h5py.Group + Group in HDF5 file + meshes : dict + Dictionary mapping IDs to mesh objects + + Returns + ------- + openmc.WeightWindows + A weight window object + """ + + id = int(group.name.split('/')[-1].lstrip('weight_windows')) + mesh_id = group['mesh'][()] + ptype = group['particle_type'][()].decode() + ebins = group['energy_bins'][()] + lower_ww_bounds = group['lower_ww_bounds'][()] + upper_ww_bounds = group['upper_ww_bounds'][()] + survival_ratio = group['survival_ratio'][()] + max_split = group['max_split'][()] + weight_cutoff = group['weight_cutoff'][()] + + return cls( + mesh=meshes[mesh_id], + lower_ww_bounds=lower_ww_bounds, + upper_ww_bounds=upper_ww_bounds, + energy_bins=ebins, + particle_type=ptype, + survival_ratio=survival_ratio, + max_split=max_split, + weight_cutoff=weight_cutoff, + id=id + ) diff --git a/src/finalize.cpp b/src/finalize.cpp index 2f6c65aaf..65880c0b8 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -25,6 +25,7 @@ #include "openmc/thermal.h" #include "openmc/timer.h" #include "openmc/volume_calc.h" +#include "openmc/weight_windows.h" #include "xtensor/xview.hpp" @@ -47,6 +48,7 @@ void free_memory() free_memory_tally(); free_memory_bank(); free_memory_plot(); + free_memory_weight_windows(); if (mpi::master) { free_memory_cmfd(); } @@ -82,6 +84,7 @@ int openmc_finalize() settings::legendre_to_tabular_points = -1; settings::material_cell_offsets = true; settings::max_particles_in_flight = 100000; + settings::max_splits = 1000; settings::n_inactive = 0; settings::n_particles = -1; settings::output_summary = true; @@ -113,6 +116,7 @@ int openmc_finalize() settings::verbosity = 7; settings::weight_cutoff = 0.25; settings::weight_survive = 1.0; + settings::weight_windows_on = false; settings::write_all_tracks = false; settings::write_initial_source = false; diff --git a/src/physics.cpp b/src/physics.cpp index bb9873bd1..196c7dc3b 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -23,6 +23,7 @@ #include "openmc/string_utils.h" #include "openmc/tallies/tally.h" #include "openmc/thermal.h" +#include "openmc/weight_windows.h" #include @@ -44,9 +45,15 @@ void collision(Particle& p) switch (p.type()) { case ParticleType::neutron: sample_neutron_reaction(p); + if (settings::weight_windows_on) { + apply_weight_windows(p); + } break; case ParticleType::photon: sample_photon_reaction(p); + if (settings::weight_windows_on) { + apply_weight_windows(p); + } break; case ParticleType::electron: sample_electron_reaction(p); diff --git a/src/settings.cpp b/src/settings.cpp index 37762517b..795b164ec 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -27,6 +27,7 @@ #include "openmc/string_utils.h" #include "openmc/tallies/trigger.h" #include "openmc/volume_calc.h" +#include "openmc/weight_windows.h" #include "openmc/xml_interface.h" namespace openmc { @@ -67,6 +68,7 @@ bool trigger_on {false}; bool trigger_predict {false}; bool ufs_on {false}; bool urr_ptables_on {true}; +bool weight_windows_on {false}; bool write_all_tracks {false}; bool write_initial_source {false}; @@ -92,6 +94,7 @@ int max_order {0}; int n_log_bins {8000}; int n_batches; int n_max_batches; +int max_splits {1000}; ResScatMethod res_scat_method {ResScatMethod::rvs}; double res_scat_energy_min {0.01}; double res_scat_energy_max {1000.0}; @@ -851,6 +854,23 @@ void read_settings_xml() if (check_for_node(root, "material_cell_offsets")) { material_cell_offsets = get_node_value_bool(root, "material_cell_offsets"); } + + // Weight window information + for (pugi::xml_node node_ww : root.children("weight_windows")) { + variance_reduction::weight_windows.emplace_back( + std::make_unique(node_ww)); + + // Enable weight windows by default if one or more are present + settings::weight_windows_on = true; + } + + if (check_for_node(root, "weight_windows_on")) { + weight_windows_on = get_node_value_bool(root, "weight_windows_on"); + } + + if (check_for_node(root, "max_splits")) { + settings::max_splits = std::stoi(get_node_value(root, "max_splits")); + } } void free_memory_settings() diff --git a/src/simulation.cpp b/src/simulation.cpp index 35d76ed87..772c90950 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -487,6 +487,9 @@ void initialize_history(Particle& p, int64_t index_source) // Reset particle event counter p.n_event() = 0; + // Reset split counter + p.n_split() = 0; + // set random number seed int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) * settings::n_particles + diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp new file mode 100644 index 000000000..c3859d7d1 --- /dev/null +++ b/src/weight_windows.cpp @@ -0,0 +1,264 @@ +#include "openmc/weight_windows.h" + +#include "openmc/error.h" +#include "openmc/file_utils.h" +#include "openmc/hdf5_interface.h" +#include "openmc/particle.h" +#include "openmc/particle_data.h" +#include "openmc/search.h" +#include "openmc/xml_interface.h" + +#include +#include + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace variance_reduction { + +std::unordered_map ww_map; +openmc::vector> weight_windows; + +} // namespace variance_reduction + +//============================================================================== +// Non-member functions +//============================================================================== + +void apply_weight_windows(Particle& p) +{ + // skip dead or no energy + if (p.E() <= 0 || !p.alive()) + return; + + bool in_domain = false; + // TODO: this is a linear search - should do something more clever + WeightWindow weight_window; + for (const auto& ww : variance_reduction::weight_windows) { + weight_window = ww->get_weight_window(p); + if (weight_window.is_valid()) + break; + } + // particle is not in any of the ww domains, do nothing + if (!weight_window.is_valid()) + return; + + // get the paramters + double weight = p.wgt(); + + // first check to see if particle should be killed for weight cutoff + if (p.wgt() < weight_window.weight_cutoff) { + p.alive() = false; + p.wgt() = 0.0; + return; + } + + // if particle's weight is above the weight window split until they are within + // the window + if (weight > weight_window.upper_weight) { + // do not further split the particle if above the limit + if (p.n_split() >= settings::max_splits) + return; + + double n_split = std::ceil(weight / weight_window.upper_weight); + double max_split = weight_window.max_split; + n_split = std::min(n_split, max_split); + + p.n_split() += n_split; + + // Create secondaries and divide weight among all particles + int i_split = std::round(n_split); + for (int l = 0; l < i_split - 1; l++) { + p.create_secondary(weight / n_split, p.u(), p.E(), p.type()); + } + // remaining weight is applied to current particle + p.wgt() = weight / n_split; + + } else if (weight <= weight_window.lower_weight) { + // if the particle weight is below the window, play Russian roulette + double weight_survive = + std::min(weight * weight_window.max_split, weight_window.survival_weight); + if (weight_survive * prn(p.current_seed()) <= weight) { + p.wgt() = weight_survive; + } else { + p.alive() = false; + p.wgt() = 0.0; + } + // else particle is in the window, continue as normal + } +} + +void free_memory_weight_windows() +{ + variance_reduction::ww_map.clear(); + variance_reduction::weight_windows.clear(); +} + +//============================================================================== +// WeightWindowSettings implementation +//============================================================================== + +WeightWindows::WeightWindows(pugi::xml_node node) +{ + // Make sure required elements are present + const vector required_elems { + "id", "particle_type", "energy_bins", "lower_ww_bounds", "upper_ww_bounds"}; + for (const auto& elem : required_elems) { + if (!check_for_node(node, elem.c_str())) { + fatal_error(fmt::format("Must specify <{}> for weight windows.", elem)); + } + } + + // Get weight windows ID + int32_t id = std::stoi(get_node_value(node, "id")); + this->set_id(id); + + // get the particle type + auto particle_type_str = std::string(get_node_value(node, "particle_type")); + particle_type_ = openmc::str_to_particle_type(particle_type_str); + + // Determine associated mesh + int32_t mesh_id = std::stoi(get_node_value(node, "mesh")); + mesh_idx_ = model::mesh_map.at(mesh_id); + + // energy bounds + energy_bins_ = get_node_array(node, "energy_bins"); + + // read the lower/upper weight bounds + lower_ww_ = get_node_array(node, "lower_ww_bounds"); + upper_ww_ = get_node_array(node, "upper_ww_bounds"); + + // get the survival value - optional + if (check_for_node(node, "survival_ratio")) { + survival_ratio_ = std::stod(get_node_value(node, "survival_ratio")); + if (survival_ratio_ <= 1) + fatal_error("Survival to lower weight window ratio must bigger than 1 " + "and less than the upper to lower weight window ratio."); + } + + // get the max split - optional + if (check_for_node(node, "max_split")) { + max_split_ = std::stod(get_node_value(node, "max_split")); + if (max_split_ <= 1) + fatal_error("max split must be larger than 1"); + } + + // weight cutoff - optional + if (check_for_node(node, "weight_cutoff")) { + weight_cutoff_ = std::stod(get_node_value(node, "weight_cutoff")); + if (weight_cutoff_ <= 0) + fatal_error("weight_cutoff must be larger than 0"); + if (weight_cutoff_ > 1) + fatal_error("weight_cutoff must be less than 1"); + } + + // make sure that the upper and lower bounds have the same size + if (upper_ww_.size() != lower_ww_.size()) { + fatal_error("The upper and lower weight window lengths do not match."); + } + + // num spatial*energy bins must match num weight bins + int num_spatial_bins = this->mesh().n_bins(); + int num_energy_bins = energy_bins_.size() - 1; + int num_weight_bins = lower_ww_.size(); + if (num_weight_bins != num_spatial_bins * num_energy_bins) { + auto err_msg = + fmt::format("In weight window domain {} the number of spatial " + "energy/spatial bins ({}) does not match the number " + "of weight bins ({})", + id_, num_energy_bins, num_weight_bins); + fatal_error(err_msg); + } +} + +void WeightWindows::set_id(int32_t id) +{ + Expects(id >= 0 || id == C_NONE); + + // Clear entry in mesh map in case one was already assigned + if (id_ != C_NONE) { + variance_reduction::ww_map.erase(id_); + id_ = C_NONE; + } + + // Ensure no other mesh has the same ID + if (variance_reduction::ww_map.find(id) != variance_reduction::ww_map.end()) { + throw std::runtime_error { + fmt::format("Two weight windows have the same ID: {}", id)}; + } + + // If no ID is specified, auto-assign the next ID in the sequence + if (id == C_NONE) { + id = 0; + for (const auto& m : variance_reduction::weight_windows) { + id = std::max(id, m->id_); + } + ++id; + } + + // Update ID and entry in the mesh map + id_ = id; + variance_reduction::ww_map[id] = + variance_reduction::weight_windows.size() - 1; +} + +WeightWindow WeightWindows::get_weight_window(const Particle& p) const +{ + // check for particle type + if (particle_type_ != p.type()) { + return {}; + } + + // Get mesh index for particle's position + const auto& mesh = this->mesh(); + int ww_index = mesh.get_bin(p.r()); + + // particle is outside the weight window mesh + if (ww_index < 0) + return {}; + + // particle energy + double E = p.E(); + + // check to make sure energy is in range, expects sorted energy values + if (E < energy_bins_.front() || E > energy_bins_.back()) + return {}; + + // get the mesh bin in energy group + int energy_bin = + lower_bound_index(energy_bins_.begin(), energy_bins_.end(), E); + + // indices now points to the correct weight for the given energy + ww_index += energy_bin * mesh.n_bins(); + + // Create individual weight window + WeightWindow ww; + ww.lower_weight = lower_ww_[ww_index]; + ww.upper_weight = upper_ww_[ww_index]; + ww.survival_weight = ww.lower_weight * survival_ratio_; + ww.max_split = max_split_; + ww.weight_cutoff = weight_cutoff_; + return ww; +} + +void WeightWindows::to_hdf5(hid_t group) const +{ + hid_t ww_group = create_group(group, fmt::format("weight_windows {}", id_)); + + write_dataset( + ww_group, "particle_type", openmc::particle_type_to_str(particle_type_)); + write_dataset(ww_group, "energy_bins", energy_bins_); + 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_split", max_split_); + write_dataset(ww_group, "weight_cutoff", weight_cutoff_); + write_dataset(ww_group, "mesh", this->mesh().id_); + + close_group(ww_group); +} + +} // namespace openmc diff --git a/tests/__init__.py b/tests/__init__.py index 5c21cd938..0228348a2 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -1,13 +1,23 @@ from contextlib import contextmanager import os +from shutil import copy import tempfile @contextmanager -def cdtemp(): - """Context manager to change to/return from a tmpdir.""" +def cdtemp(files=None): + """Context manager to change to/return from a tmpdir. + + Parameters + ---------- + files : Iterable of str or Path-like + Set of files to copy into the temporary directory + """ with tempfile.TemporaryDirectory() as tmpdir: cwd = os.getcwd() + if files: + for file in files: + copy(file, tmpdir, follow_symlinks=True) try: os.chdir(tmpdir) yield diff --git a/tests/regression_tests/weightwindows/__init__.py b/tests/regression_tests/weightwindows/__init__.py new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/tests/regression_tests/weightwindows/__init__.py @@ -0,0 +1 @@ + diff --git a/tests/regression_tests/weightwindows/inputs_true.dat b/tests/regression_tests/weightwindows/inputs_true.dat new file mode 100644 index 000000000..5459a09b9 --- /dev/null +++ b/tests/regression_tests/weightwindows/inputs_true.dat @@ -0,0 +1,94 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 200 + 2 + + + 0.001 0.001 0.001 + + + 14000000.0 1.0 + + + true + + 2 + neutron + 0.0 0.5 20000000.0 + -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 7.695031155172816e-14 -1.0 7.476545089278482e-06 -1.0 2.1612150425221703e-06 -1.0 -1.0 -1.0 2.1360401784344975e-18 -1.0 1.3455341306162382e-05 -1.0 3.353295403842037e-05 -1.0 3.0160758454323657e-07 -1.0 1.0262431550731535e-13 -1.0 1.6859219614058024e-19 -1.0 3.853097455417877e-06 -1.0 4.354946981329799e-05 -1.0 7.002861620919232e-08 -1.0 -1.0 -1.0 -1.0 -1.0 1.6205064883026195e-11 -1.0 1.7041669224797384e-09 -1.0 1.8747824667478637e-14 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 5.619835883512353e-14 -1.0 1.951356547084204e-15 -1.0 8.48176145671274e-10 -1.0 -1.0 -1.0 1.4508262073256659e-13 -1.0 4.109066205029878e-05 -1.0 0.00020573181450240786 -1.0 2.124999182004655e-05 -1.0 5.5379154134462336e-08 -1.0 1.5550455959178486e-06 -1.0 0.0007213722022489493 -1.0 0.007673145137236567 -1.0 0.0008175033526488423 -1.0 1.008394964750947e-06 -1.0 6.848642863465585e-07 -1.0 0.0010737081832502356 -1.0 0.007197662162700777 -1.0 0.0007151459162818186 -1.0 2.253382683081525e-06 -1.0 1.3181562352445092e-10 -1.0 1.93847081892941e-05 -1.0 0.00046029617115127556 -1.0 2.5629068330446215e-05 -1.0 5.504813693036933e-12 -1.0 -1.0 -1.0 3.3164825349503764e-16 -1.0 3.937160538176268e-07 -1.0 -1.0 -1.0 -1.0 -1.0 1.5172912057312398e-14 -1.0 6.315154082213375e-07 -1.0 3.855884234344845e-06 -1.0 1.7770136411255912e-05 -1.0 1.6421491993751715e-11 -1.0 5.97653814110781e-06 -1.0 0.001959411999677113 -1.0 0.013330619853379314 -1.0 0.0017950613169359904 -1.0 2.4306375693722205e-06 -1.0 6.922231352729155e-05 -1.0 0.08657573566388353 -1.0 -1.0 -1.0 0.0835648710074336 -1.0 0.00016419031150495913 -1.0 9.03108551826469e-05 -1.0 0.0848872967381878 -1.0 -1.0 -1.0 0.08288806458368345 -1.0 5.6949597066784976e-05 -1.0 8.037131813528501e-07 -1.0 0.0016637126543321873 -1.0 0.01664109232689093 -1.0 0.0018327593437608 -1.0 6.542700644645627e-07 -1.0 -1.0 -1.0 6.655926357459275e-08 -1.0 5.783974359937857e-06 -1.0 2.679340942769232e-06 -1.0 -1.0 -1.0 -1.0 -1.0 2.1040609012865134e-05 -1.0 2.9282097947816224e-05 -1.0 1.284759166582202e-05 -1.0 1.0092905083158804e-11 -1.0 1.3228013765525015e-05 -1.0 0.007136750029575816 -1.0 0.06539005235506369 -1.0 0.0064803314120165335 -1.0 2.6011787393916806e-06 -1.0 8.887234042637684e-05 -1.0 0.5 -1.0 -1.0 -1.0 0.4877856021170283 -1.0 0.00017699858357744463 -1.0 0.0001670700632870335 -1.0 0.4899999304038166 -1.0 -1.0 -1.0 0.48573842213878143 -1.0 0.00021534568699157805 -1.0 7.84289689494925e-06 -1.0 0.0063790654507599135 -1.0 0.06550426960512012 -1.0 0.006639379668946672 -1.0 7.754700849337902e-06 -1.0 2.7207433165796702e-11 -1.0 5.985108617124989e-06 -1.0 3.634644995026692e-05 -1.0 3.932389299316285e-06 -1.0 -1.0 -1.0 -1.0 -1.0 1.271128669411295e-07 -1.0 8.846493399850663e-06 -1.0 2.499931357628383e-07 -1.0 8.533030345699353e-17 -1.0 1.6699254277734109e-06 -1.0 0.0017541380096893788 -1.0 0.015148978428271613 -1.0 0.0016279794552427574 -1.0 2.380955631371226e-06 -1.0 2.1648918040467363e-05 -1.0 0.0833931797381589 -1.0 -1.0 -1.0 0.08285899874014539 -1.0 2.762955748645507e-05 -1.0 5.1161740516205715e-05 -1.0 0.08095280568305474 -1.0 -1.0 -1.0 0.08636569925236791 -1.0 3.452537179033895e-05 -1.0 3.5981031766416143e-06 -1.0 0.001570829313580841 -1.0 0.01621821782442112 -1.0 0.0015447823995470042 -1.0 7.989579285134734e-06 -1.0 -1.0 -1.0 2.594951085377588e-06 -1.0 4.323138781337963e-05 -1.0 6.948292975998466e-07 -1.0 -1.0 -1.0 -1.0 -1.0 2.4469456008493614e-09 -1.0 9.041149893307806e-07 -1.0 2.2060992580622125e-11 -1.0 -1.0 -1.0 1.681305446488884e-11 -1.0 7.091263906557291e-05 -1.0 0.00027464244062887683 -1.0 1.28584557092134e-05 -1.0 2.0295213621716706e-09 -1.0 1.2911827198500372e-08 -1.0 0.0009127050763987511 -1.0 0.007385384405891666 -1.0 0.0010019432401508087 -1.0 1.080711295768551e-05 -1.0 8.020767115181574e-07 -1.0 0.0009252532821720597 -1.0 0.007188410694198128 -1.0 0.0007245451610090619 -1.0 2.5157355790201773e-07 -1.0 1.4825219199133353e-12 -1.0 4.756872959630257e-05 -1.0 0.0002668579837809081 -1.0 2.8226276944532696e-05 -1.0 5.182269672663266e-10 -1.0 -1.0 -1.0 2.936268747135548e-14 -1.0 2.5447281241730366e-07 -1.0 3.89841667138598e-08 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 2.639268952045579e-10 -1.0 3.7051518220083886e-06 -1.0 1.1333627604823673e-09 -1.0 -1.0 -1.0 -1.0 -1.0 1.2296747260635405e-06 -1.0 7.962899789920809e-06 -1.0 8.211982683649991e-09 -1.0 -1.0 -1.0 -1.0 -1.0 2.5639656866250453e-06 -1.0 2.2487617828938326e-05 -1.0 3.2799045883372075e-06 -1.0 4.076762579823379e-17 -1.0 1.5476769237571295e-14 -1.0 1.0305672698745657e-11 -1.0 1.957092814065688e-05 -1.0 3.8523247550378815e-12 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 2.7697233171399416e-13 -1.0 -1.0 -1.0 -1.0 -1.0 + -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 7.695031155172816e-13 -10.0 7.476545089278483e-05 -10.0 2.1612150425221703e-05 -10.0 -10.0 -10.0 2.1360401784344976e-17 -10.0 0.00013455341306162382 -10.0 0.00033532954038420373 -10.0 3.0160758454323656e-06 -10.0 1.0262431550731535e-12 -10.0 1.6859219614058023e-18 -10.0 3.853097455417877e-05 -10.0 0.00043549469813297995 -10.0 7.002861620919232e-07 -10.0 -10.0 -10.0 -10.0 -10.0 1.6205064883026195e-10 -10.0 1.7041669224797384e-08 -10.0 1.8747824667478637e-13 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 5.619835883512353e-13 -10.0 1.951356547084204e-14 -10.0 8.48176145671274e-09 -10.0 -10.0 -10.0 1.4508262073256658e-12 -10.0 0.0004109066205029878 -10.0 0.0020573181450240785 -10.0 0.00021249991820046548 -10.0 5.537915413446234e-07 -10.0 1.5550455959178486e-05 -10.0 0.0072137220224894934 -10.0 0.07673145137236567 -10.0 0.008175033526488424 -10.0 1.0083949647509469e-05 -10.0 6.848642863465585e-06 -10.0 0.010737081832502356 -10.0 0.07197662162700777 -10.0 0.007151459162818186 -10.0 2.253382683081525e-05 -10.0 1.3181562352445092e-09 -10.0 0.000193847081892941 -10.0 0.004602961711512756 -10.0 0.00025629068330446217 -10.0 5.504813693036933e-11 -10.0 -10.0 -10.0 3.3164825349503764e-15 -10.0 3.937160538176268e-06 -10.0 -10.0 -10.0 -10.0 -10.0 1.51729120573124e-13 -10.0 6.3151540822133754e-06 -10.0 3.855884234344845e-05 -10.0 0.0001777013641125591 -10.0 1.6421491993751715e-10 -10.0 5.97653814110781e-05 -10.0 0.01959411999677113 -10.0 0.13330619853379314 -10.0 0.017950613169359905 -10.0 2.4306375693722206e-05 -10.0 0.0006922231352729155 -10.0 0.8657573566388354 -10.0 -10.0 -10.0 0.8356487100743359 -10.0 0.0016419031150495913 -10.0 0.000903108551826469 -10.0 0.848872967381878 -10.0 -10.0 -10.0 0.8288806458368345 -10.0 0.0005694959706678497 -10.0 8.037131813528501e-06 -10.0 0.016637126543321872 -10.0 0.1664109232689093 -10.0 0.018327593437608 -10.0 6.5427006446456266e-06 -10.0 -10.0 -10.0 6.655926357459275e-07 -10.0 5.7839743599378564e-05 -10.0 2.679340942769232e-05 -10.0 -10.0 -10.0 -10.0 -10.0 0.00021040609012865135 -10.0 0.00029282097947816227 -10.0 0.0001284759166582202 -10.0 1.0092905083158804e-10 -10.0 0.00013228013765525016 -10.0 0.07136750029575815 -10.0 0.6539005235506369 -10.0 0.06480331412016534 -10.0 2.6011787393916804e-05 -10.0 0.0008887234042637684 -10.0 5.0 -10.0 -10.0 -10.0 4.877856021170283 -10.0 0.0017699858357744464 -10.0 0.001670700632870335 -10.0 4.899999304038166 -10.0 -10.0 -10.0 4.8573842213878144 -10.0 0.0021534568699157807 -10.0 7.84289689494925e-05 -10.0 0.06379065450759913 -10.0 0.6550426960512012 -10.0 0.06639379668946672 -10.0 7.754700849337902e-05 -10.0 2.7207433165796704e-10 -10.0 5.985108617124989e-05 -10.0 0.0003634644995026692 -10.0 3.9323892993162844e-05 -10.0 -10.0 -10.0 -10.0 -10.0 1.271128669411295e-06 -10.0 8.846493399850663e-05 -10.0 2.499931357628383e-06 -10.0 8.533030345699353e-16 -10.0 1.669925427773411e-05 -10.0 0.017541380096893787 -10.0 0.15148978428271614 -10.0 0.016279794552427576 -10.0 2.380955631371226e-05 -10.0 0.00021648918040467362 -10.0 0.8339317973815891 -10.0 -10.0 -10.0 0.828589987401454 -10.0 0.0002762955748645507 -10.0 0.0005116174051620571 -10.0 0.8095280568305474 -10.0 -10.0 -10.0 0.8636569925236791 -10.0 0.0003452537179033895 -10.0 3.5981031766416144e-05 -10.0 0.015708293135808408 -10.0 0.16218217824421122 -10.0 0.015447823995470043 -10.0 7.989579285134734e-05 -10.0 -10.0 -10.0 2.594951085377588e-05 -10.0 0.0004323138781337963 -10.0 6.948292975998466e-06 -10.0 -10.0 -10.0 -10.0 -10.0 2.4469456008493615e-08 -10.0 9.041149893307805e-06 -10.0 2.2060992580622125e-10 -10.0 -10.0 -10.0 1.681305446488884e-10 -10.0 0.0007091263906557291 -10.0 0.002746424406288768 -10.0 0.000128584557092134 -10.0 2.0295213621716707e-08 -10.0 1.2911827198500372e-07 -10.0 0.009127050763987512 -10.0 0.07385384405891665 -10.0 0.010019432401508087 -10.0 0.0001080711295768551 -10.0 8.020767115181574e-06 -10.0 0.009252532821720597 -10.0 0.07188410694198127 -10.0 0.0072454516100906195 -10.0 2.515735579020177e-06 -10.0 1.4825219199133352e-11 -10.0 0.0004756872959630257 -10.0 0.0026685798378090807 -10.0 0.000282262769445327 -10.0 5.182269672663266e-09 -10.0 -10.0 -10.0 2.9362687471355483e-13 -10.0 2.5447281241730365e-06 -10.0 3.89841667138598e-07 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 2.639268952045579e-09 -10.0 3.7051518220083885e-05 -10.0 1.1333627604823672e-08 -10.0 -10.0 -10.0 -10.0 -10.0 1.2296747260635405e-05 -10.0 7.96289978992081e-05 -10.0 8.211982683649991e-08 -10.0 -10.0 -10.0 -10.0 -10.0 2.5639656866250452e-05 -10.0 0.00022487617828938325 -10.0 3.2799045883372076e-05 -10.0 4.076762579823379e-16 -10.0 1.5476769237571296e-13 -10.0 1.0305672698745658e-10 -10.0 0.0001957092814065688 -10.0 3.852324755037882e-11 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 2.7697233171399415e-12 -10.0 -10.0 -10.0 -10.0 -10.0 + 3 + 10 + 1e-38 + + + 5 6 7 + -240 -240 -240 + 240 240 240 + + + 2 + neutron + 0.0 0.5 20000000.0 + -1.0 -1.0 -1.0 5.447493368067179e-10 -1.0 4.0404626147344e-13 -1.0 9.382800796426366e-17 -1.0 -1.0 -1.0 5.752120832330886e-11 -1.0 3.520674728719278e-05 -1.0 8.76041935303527e-05 -1.0 3.136829920007009e-05 -1.0 -1.0 -1.0 4.233410073487772e-06 -1.0 0.00025519632243117644 -1.0 0.0004777967757842694 -1.0 0.00035983547999740665 -1.0 4.6258592060517866e-07 -1.0 4.587143797469485e-08 -1.0 0.00012918433677593975 -1.0 0.000469589568463265 -1.0 8.862760180332873e-05 -1.0 3.2790914744011534e-06 -1.0 1.1514579117272292e-10 -1.0 2.2950716170641463e-05 -1.0 6.92296299242814e-05 -1.0 2.1606550732170693e-05 -1.0 2.617878055106949e-16 -1.0 -1.0 -1.0 -1.0 -1.0 1.722191123879782e-08 -1.0 -1.0 -1.0 -1.0 -1.0 2.1234980618763222e-13 -1.0 1.4628847959717225e-05 -1.0 2.777791741571174e-05 -1.0 9.53380658669839e-06 -1.0 1.078687446476881e-14 -1.0 3.0571050855707084e-06 -1.0 0.0008231482263927564 -1.0 0.0038811025614225057 -1.0 0.0012297303081805393 -1.0 4.6842350882367666e-05 -1.0 0.00016757997590048827 -1.0 0.011724714376848187 -1.0 0.06827568359008944 -1.0 0.011440381012864086 -1.0 0.0003782322077390033 -1.0 5.098747893974248e-05 -1.0 0.01208261392046878 -1.0 0.06428782790571179 -1.0 0.011294708566804167 -1.0 8.168081568291806e-05 -1.0 2.34750482787598e-05 -1.0 0.0008787411098754915 -1.0 0.0051513292132259495 -1.0 0.0011028020390507682 -1.0 1.0882011018684852e-06 -1.0 1.0588257623722672e-07 -1.0 2.7044817176109556e-05 -1.0 2.6211239246796104e-05 -1.0 2.0707648529746894e-06 -1.0 -1.0 -1.0 1.3806577785135177e-06 -1.0 0.00010077009726691265 -1.0 0.00043130977711323417 -1.0 0.0001335621430311088 -1.0 3.792209399652396e-06 -1.0 0.00022057970761624768 -1.0 0.01884234453311008 -1.0 0.11510120013926417 -1.0 0.01825795827020865 -1.0 0.00020848309220162036 -1.0 0.0013332554437480327 -1.0 0.49317081643076643 -1.0 -1.0 -1.0 0.4899952354672477 -1.0 0.0014957033281168012 -1.0 0.001550036616828061 -1.0 0.4897735047310072 -1.0 -1.0 -1.0 0.5 -1.0 0.0011076796092102186 -1.0 0.0002792896409558745 -1.0 0.01721806295736377 -1.0 0.11891251244906934 -1.0 0.019040707071165303 -1.0 0.00017329528628152923 -1.0 4.670673837764874e-08 -1.0 0.00012637534695610975 -1.0 0.00027052054298385255 -1.0 0.0001497958445994315 -1.0 3.2478301680021854e-08 -1.0 7.239087311622819e-08 -1.0 0.00046693784586018913 -1.0 0.001537590300149361 -1.0 0.00023686818807101741 -1.0 8.616027294555462e-07 -1.0 0.0002362355168421952 -1.0 0.05708498338352203 -1.0 0.42372713683725016 -1.0 0.054902816619122045 -1.0 0.00030002326388383244 -1.0 0.0023295401596005313 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0027847329307350423 -1.0 0.003432919201293737 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.002672923773363389 -1.0 0.0003182108781955478 -1.0 0.05639071613735579 -1.0 0.4346984808472855 -1.0 0.05500108683629957 -1.0 0.00024175998138228048 -1.0 2.774067049392004e-07 -1.0 0.0003257646134837451 -1.0 0.0012648333883989995 -1.0 0.0002981891705231209 -1.0 4.986698316296507e-06 -1.0 3.419181071650816e-06 -1.0 0.00010814452764860558 -1.0 0.0006934466049804711 -1.0 0.00025106741008066317 -1.0 5.515522196186723e-06 -1.0 0.00014242053195529864 -1.0 0.018541112907970003 -1.0 0.12364069326350118 -1.0 0.017602971074583987 -1.0 0.00021305826462238593 -1.0 0.0015016500045374082 -1.0 0.4903495510314233 -1.0 -1.0 -1.0 0.4952432605379829 -1.0 0.0017693279949094218 -1.0 0.0009304497896205121 -1.0 0.4873897581604287 -1.0 -1.0 -1.0 0.4981497862127628 -1.0 0.0011273340972061078 -1.0 0.00014126153603265097 -1.0 0.01729816946709781 -1.0 0.11819165666519223 -1.0 0.017802624844035688 -1.0 0.0001319677450153934 -1.0 2.210707325798542e-07 -1.0 0.00017287873387910357 -1.0 0.0005197081842709918 -1.0 9.256616392487858e-05 -1.0 1.0325207105027347e-08 -1.0 2.027381220364608e-16 -1.0 1.2257360246712472e-05 -1.0 3.8561696337086434e-05 -1.0 1.8064550990294753e-05 -1.0 7.691527781829036e-17 -1.0 1.298062961617927e-05 -1.0 0.001178616843705956 -1.0 0.004441328419270369 -1.0 0.0008919466185572302 -1.0 5.272021975060514e-06 -1.0 0.00016079768514218546 -1.0 0.011128424196387618 -1.0 0.06433819434431046 -1.0 0.010955637485319249 -1.0 8.370896020975689e-05 -1.0 0.00022234471595443313 -1.0 0.01102706140784316 -1.0 0.06684484191345574 -1.0 0.011113024641218336 -1.0 0.00012558879474108074 -1.0 1.0044390728858277e-05 -1.0 0.000808100901275472 -1.0 0.005384625561469684 -1.0 0.001146714816952061 -1.0 3.5740607788414942e-06 -1.0 2.2512327083641865e-10 -1.0 6.717625881826167e-05 -1.0 6.949587369646497e-05 -1.0 2.9093345975075226e-05 -1.0 -1.0 -1.0 -1.0 -1.0 8.582011716859219e-15 -1.0 8.304845400451923e-08 -1.0 4.5757238372058234e-11 -1.0 -1.0 -1.0 1.3116807133030035e-15 -1.0 2.6866894852481906e-05 -1.0 4.3166052799897375e-05 -1.0 6.401055165502384e-06 -1.0 5.726152701849347e-16 -1.0 8.541606727879512e-07 -1.0 0.0002029396545382813 -1.0 0.0004299195570308035 -1.0 8.544159282151438e-05 -1.0 1.0995322456277382e-06 -1.0 1.289814511028015e-06 -1.0 0.00011042047857307253 -1.0 0.0006870228272606287 -1.0 0.0001399807456616496 -1.0 4.279460329782799e-09 -1.0 1.2907171695294846e-13 -1.0 1.9257248009114363e-05 -1.0 9.332371598523186e-05 -1.0 1.3106335518133802e-05 -1.0 6.574706543646946e-16 -1.0 -1.0 -1.0 3.8270536025799805e-15 -1.0 2.101790478097517e-09 -1.0 9.904481358675478e-13 -1.0 -1.0 + -10.0 -10.0 -10.0 5.447493368067179e-09 -10.0 4.0404626147344005e-12 -10.0 9.382800796426367e-16 -10.0 -10.0 -10.0 5.752120832330886e-10 -10.0 0.0003520674728719278 -10.0 0.0008760419353035269 -10.0 0.0003136829920007009 -10.0 -10.0 -10.0 4.233410073487772e-05 -10.0 0.0025519632243117644 -10.0 0.004777967757842694 -10.0 0.0035983547999740664 -10.0 4.625859206051786e-06 -10.0 4.587143797469485e-07 -10.0 0.0012918433677593976 -10.0 0.0046958956846326495 -10.0 0.0008862760180332873 -10.0 3.279091474401153e-05 -10.0 1.151457911727229e-09 -10.0 0.00022950716170641464 -10.0 0.0006922962992428139 -10.0 0.00021606550732170693 -10.0 2.6178780551069492e-15 -10.0 -10.0 -10.0 -10.0 -10.0 1.722191123879782e-07 -10.0 -10.0 -10.0 -10.0 -10.0 2.123498061876322e-12 -10.0 0.00014628847959717226 -10.0 0.0002777791741571174 -10.0 9.53380658669839e-05 -10.0 1.078687446476881e-13 -10.0 3.0571050855707086e-05 -10.0 0.008231482263927564 -10.0 0.03881102561422506 -10.0 0.012297303081805393 -10.0 0.00046842350882367664 -10.0 0.0016757997590048828 -10.0 0.11724714376848187 -10.0 0.6827568359008944 -10.0 0.11440381012864086 -10.0 0.003782322077390033 -10.0 0.0005098747893974248 -10.0 0.1208261392046878 -10.0 0.6428782790571179 -10.0 0.11294708566804167 -10.0 0.0008168081568291806 -10.0 0.000234750482787598 -10.0 0.008787411098754914 -10.0 0.0515132921322595 -10.0 0.011028020390507681 -10.0 1.0882011018684853e-05 -10.0 1.0588257623722672e-06 -10.0 0.00027044817176109554 -10.0 0.000262112392467961 -10.0 2.0707648529746893e-05 -10.0 -10.0 -10.0 1.3806577785135176e-05 -10.0 0.0010077009726691265 -10.0 0.004313097771132342 -10.0 0.001335621430311088 -10.0 3.792209399652396e-05 -10.0 0.0022057970761624767 -10.0 0.1884234453311008 -10.0 1.1510120013926417 -10.0 0.1825795827020865 -10.0 0.0020848309220162036 -10.0 0.013332554437480326 -10.0 4.931708164307665 -10.0 -10.0 -10.0 4.899952354672477 -10.0 0.014957033281168012 -10.0 0.01550036616828061 -10.0 4.897735047310072 -10.0 -10.0 -10.0 5.0 -10.0 0.011076796092102187 -10.0 0.002792896409558745 -10.0 0.17218062957363772 -10.0 1.1891251244906933 -10.0 0.19040707071165303 -10.0 0.0017329528628152924 -10.0 4.670673837764874e-07 -10.0 0.0012637534695610975 -10.0 0.0027052054298385255 -10.0 0.001497958445994315 -10.0 3.2478301680021856e-07 -10.0 7.239087311622819e-07 -10.0 0.004669378458601891 -10.0 0.01537590300149361 -10.0 0.002368681880710174 -10.0 8.61602729455546e-06 -10.0 0.002362355168421952 -10.0 0.5708498338352204 -10.0 4.237271368372502 -10.0 0.5490281661912204 -10.0 0.0030002326388383245 -10.0 0.023295401596005315 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.027847329307350423 -10.0 0.03432919201293737 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.026729237733633893 -10.0 0.0031821087819554777 -10.0 0.5639071613735579 -10.0 4.346984808472855 -10.0 0.5500108683629957 -10.0 0.0024175998138228046 -10.0 2.774067049392004e-06 -10.0 0.0032576461348374514 -10.0 0.012648333883989995 -10.0 0.002981891705231209 -10.0 4.986698316296507e-05 -10.0 3.419181071650816e-05 -10.0 0.0010814452764860557 -10.0 0.006934466049804711 -10.0 0.0025106741008066318 -10.0 5.515522196186723e-05 -10.0 0.0014242053195529865 -10.0 0.18541112907970003 -10.0 1.2364069326350118 -10.0 0.17602971074583987 -10.0 0.0021305826462238594 -10.0 0.015016500045374082 -10.0 4.903495510314233 -10.0 -10.0 -10.0 4.952432605379829 -10.0 0.01769327994909422 -10.0 0.00930449789620512 -10.0 4.873897581604287 -10.0 -10.0 -10.0 4.981497862127628 -10.0 0.011273340972061077 -10.0 0.0014126153603265096 -10.0 0.1729816946709781 -10.0 1.1819165666519222 -10.0 0.1780262484403569 -10.0 0.001319677450153934 -10.0 2.210707325798542e-06 -10.0 0.0017287873387910357 -10.0 0.0051970818427099184 -10.0 0.0009256616392487858 -10.0 1.0325207105027347e-07 -10.0 2.027381220364608e-15 -10.0 0.00012257360246712473 -10.0 0.00038561696337086434 -10.0 0.00018064550990294753 -10.0 7.6915277818290365e-16 -10.0 0.00012980629616179268 -10.0 0.01178616843705956 -10.0 0.04441328419270369 -10.0 0.008919466185572301 -10.0 5.272021975060514e-05 -10.0 0.0016079768514218546 -10.0 0.11128424196387618 -10.0 0.6433819434431045 -10.0 0.10955637485319249 -10.0 0.0008370896020975689 -10.0 0.0022234471595443312 -10.0 0.1102706140784316 -10.0 0.6684484191345574 -10.0 0.11113024641218336 -10.0 0.0012558879474108074 -10.0 0.00010044390728858278 -10.0 0.00808100901275472 -10.0 0.05384625561469684 -10.0 0.01146714816952061 -10.0 3.574060778841494e-05 -10.0 2.2512327083641864e-09 -10.0 0.0006717625881826168 -10.0 0.0006949587369646498 -10.0 0.0002909334597507523 -10.0 -10.0 -10.0 -10.0 -10.0 8.582011716859219e-14 -10.0 8.304845400451923e-07 -10.0 4.5757238372058235e-10 -10.0 -10.0 -10.0 1.3116807133030034e-14 -10.0 0.00026866894852481903 -10.0 0.00043166052799897375 -10.0 6.401055165502385e-05 -10.0 5.726152701849347e-15 -10.0 8.541606727879512e-06 -10.0 0.0020293965453828133 -10.0 0.004299195570308035 -10.0 0.0008544159282151438 -10.0 1.0995322456277382e-05 -10.0 1.289814511028015e-05 -10.0 0.0011042047857307254 -10.0 0.006870228272606287 -10.0 0.0013998074566164962 -10.0 4.279460329782799e-08 -10.0 1.2907171695294846e-12 -10.0 0.00019257248009114364 -10.0 0.0009332371598523186 -10.0 0.000131063355181338 -10.0 6.574706543646946e-15 -10.0 -10.0 -10.0 3.8270536025799805e-14 -10.0 2.101790478097517e-08 -10.0 9.904481358675478e-12 -10.0 -10.0 + 3 + 10 + 1e-38 + + 200 + + + + + 5 10 15 + -240 -240 -240 + 240 240 240 + + + 1 + + + 0.0 0.5 20000000.0 + + + neutron photon + + + 1 2 3 + flux + + diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat new file mode 100644 index 000000000..5e8a41faf --- /dev/null +++ b/tests/regression_tests/weightwindows/results_true.dat @@ -0,0 +1 @@ +c77af605fb23bc3f03bb304b82cc2d0d9efe7826167f43202fa0ecba17e6a77fc1a91bf9f4b23f175ae98c04c69067bffb96d7d68e30a126f5fa7c0080ea8367 \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/test.py b/tests/regression_tests/weightwindows/test.py new file mode 100644 index 000000000..63d773cbc --- /dev/null +++ b/tests/regression_tests/weightwindows/test.py @@ -0,0 +1,110 @@ +from copy import deepcopy +import pytest +import numpy as np + +import openmc +from openmc.stats import Discrete, Point + +from tests.testing_harness import HashedPyAPITestHarness + + +@pytest.fixture +def model(): + model = openmc.Model() + + # materials (M4 steel alloy) + m4 = openmc.Material() + m4.set_density('g/cc', 2.3) + m4.add_nuclide('H1', 0.168018676) + m4.add_nuclide("H2", 1.93244e-05) + m4.add_nuclide("O16", 0.561814465) + m4.add_nuclide("O17", 0.00021401) + m4.add_nuclide("Na23", 0.021365) + m4.add_nuclide("Al27", 0.021343) + m4.add_nuclide("Si28", 0.187439342) + m4.add_nuclide("Si29", 0.009517714) + m4.add_nuclide("Si30", 0.006273944) + m4.add_nuclide("Ca40", 0.018026179) + m4.add_nuclide("Ca42", 0.00012031) + m4.add_nuclide("Ca43", 2.51033e-05) + m4.add_nuclide("Ca44", 0.000387892) + m4.add_nuclide("Ca46", 7.438e-07) + m4.add_nuclide("Ca48", 3.47727e-05) + m4.add_nuclide("Fe54", 0.000248179) + m4.add_nuclide("Fe56", 0.003895875) + m4.add_nuclide("Fe57", 8.99727e-05) + m4.add_nuclide("Fe58", 1.19737e-05) + + s0 = openmc.Sphere(r=240) + s1 = openmc.Sphere(r=250, boundary_type='vacuum') + + c0 = openmc.Cell(fill=m4, region=-s0) + c1 = openmc.Cell(region=+s0 & -s1) + + model.geometry = openmc.Geometry([c0, c1]) + + # settings + settings = model.settings + settings.run_mode = 'fixed source' + settings.particles = 200 + settings.batches = 2 + settings.max_splits = 200 + settings.photon_transport = True + space = Point((0.001, 0.001, 0.001)) + energy = Discrete([14E6], [1.0]) + + settings.source = openmc.Source(space=space, energy=energy) + + # tally + mesh = openmc.RegularMesh() + mesh.lower_left = (-240, -240, -240) + mesh.upper_right = (240, 240, 240) + mesh.dimension = (5, 10, 15) + + mesh_filter = openmc.MeshFilter(mesh) + + e_bnds = [0.0, 0.5, 2E7] + energy_filter = openmc.EnergyFilter(e_bnds) + + particle_filter = openmc.ParticleFilter(['neutron', 'photon']) + + tally = openmc.Tally() + tally.filters = [mesh_filter, energy_filter, particle_filter] + tally.scores = ['flux'] + + model.tallies.append(tally) + + # weight windows + + # load pre-generated weight windows + # (created using the same tally as above) + ww_n_lower_bnds = np.loadtxt('ww_n.txt') + ww_p_lower_bnds = np.loadtxt('ww_p.txt') + + # create a mesh matching the one used + # to generate the weight windows + ww_mesh = openmc.RegularMesh() + ww_mesh.lower_left = (-240, -240, -240) + ww_mesh.upper_right = (240, 240, 240) + ww_mesh.dimension = (5, 6, 7) + + ww_n = openmc.WeightWindows(ww_mesh, + ww_n_lower_bnds, + None, + 10.0, + e_bnds) + + ww_p = openmc.WeightWindows(ww_mesh, + ww_p_lower_bnds, + None, + 10.0, + e_bnds) + + model.settings.weight_windows = [ww_n, ww_p] + + return model + + +def test_weightwindows(model): + test = HashedPyAPITestHarness('statepoint.2.h5', model) + test.main() diff --git a/tests/regression_tests/weightwindows/ww_n.txt b/tests/regression_tests/weightwindows/ww_n.txt new file mode 100644 index 000000000..dbb49537b --- /dev/null +++ b/tests/regression_tests/weightwindows/ww_n.txt @@ -0,0 +1,420 @@ +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +7.695031155172815952e-14 +-1.000000000000000000e+00 +7.476545089278482294e-06 +-1.000000000000000000e+00 +2.161215042522170282e-06 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.136040178434497494e-18 +-1.000000000000000000e+00 +1.345534130616238158e-05 +-1.000000000000000000e+00 +3.353295403842037074e-05 +-1.000000000000000000e+00 +3.016075845432365699e-07 +-1.000000000000000000e+00 +1.026243155073153453e-13 +-1.000000000000000000e+00 +1.685921961405802360e-19 +-1.000000000000000000e+00 +3.853097455417877125e-06 +-1.000000000000000000e+00 +4.354946981329799215e-05 +-1.000000000000000000e+00 +7.002861620919231721e-08 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.620506488302619488e-11 +-1.000000000000000000e+00 +1.704166922479738429e-09 +-1.000000000000000000e+00 +1.874782466747863741e-14 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +5.619835883512352867e-14 +-1.000000000000000000e+00 +1.951356547084204088e-15 +-1.000000000000000000e+00 +8.481761456712739755e-10 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.450826207325665860e-13 +-1.000000000000000000e+00 +4.109066205029878035e-05 +-1.000000000000000000e+00 +2.057318145024078612e-04 +-1.000000000000000000e+00 +2.124999182004654892e-05 +-1.000000000000000000e+00 +5.537915413446233643e-08 +-1.000000000000000000e+00 +1.555045595917848555e-06 +-1.000000000000000000e+00 +7.213722022489493227e-04 +-1.000000000000000000e+00 +7.673145137236566868e-03 +-1.000000000000000000e+00 +8.175033526488423340e-04 +-1.000000000000000000e+00 +1.008394964750946906e-06 +-1.000000000000000000e+00 +6.848642863465584904e-07 +-1.000000000000000000e+00 +1.073708183250235556e-03 +-1.000000000000000000e+00 +7.197662162700776967e-03 +-1.000000000000000000e+00 +7.151459162818185784e-04 +-1.000000000000000000e+00 +2.253382683081525085e-06 +-1.000000000000000000e+00 +1.318156235244509166e-10 +-1.000000000000000000e+00 +1.938470818929409920e-05 +-1.000000000000000000e+00 +4.602961711512755646e-04 +-1.000000000000000000e+00 +2.562906833044621533e-05 +-1.000000000000000000e+00 +5.504813693036932910e-12 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +3.316482534950376351e-16 +-1.000000000000000000e+00 +3.937160538176268057e-07 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.517291205731239772e-14 +-1.000000000000000000e+00 +6.315154082213375425e-07 +-1.000000000000000000e+00 +3.855884234344844889e-06 +-1.000000000000000000e+00 +1.777013641125591191e-05 +-1.000000000000000000e+00 +1.642149199375171522e-11 +-1.000000000000000000e+00 +5.976538141107809983e-06 +-1.000000000000000000e+00 +1.959411999677113103e-03 +-1.000000000000000000e+00 +1.333061985337931354e-02 +-1.000000000000000000e+00 +1.795061316935990395e-03 +-1.000000000000000000e+00 +2.430637569372220458e-06 +-1.000000000000000000e+00 +6.922231352729155390e-05 +-1.000000000000000000e+00 +8.657573566388353237e-02 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +8.356487100743359431e-02 +-1.000000000000000000e+00 +1.641903115049591277e-04 +-1.000000000000000000e+00 +9.031085518264690009e-05 +-1.000000000000000000e+00 +8.488729673818780352e-02 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +8.288806458368344621e-02 +-1.000000000000000000e+00 +5.694959706678497628e-05 +-1.000000000000000000e+00 +8.037131813528500897e-07 +-1.000000000000000000e+00 +1.663712654332187308e-03 +-1.000000000000000000e+00 +1.664109232689093068e-02 +-1.000000000000000000e+00 +1.832759343760799985e-03 +-1.000000000000000000e+00 +6.542700644645626997e-07 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +6.655926357459274907e-08 +-1.000000000000000000e+00 +5.783974359937856763e-06 +-1.000000000000000000e+00 +2.679340942769232109e-06 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.104060901286513389e-05 +-1.000000000000000000e+00 +2.928209794781622450e-05 +-1.000000000000000000e+00 +1.284759166582201940e-05 +-1.000000000000000000e+00 +1.009290508315880392e-11 +-1.000000000000000000e+00 +1.322801376552501482e-05 +-1.000000000000000000e+00 +7.136750029575815620e-03 +-1.000000000000000000e+00 +6.539005235506369085e-02 +-1.000000000000000000e+00 +6.480331412016533503e-03 +-1.000000000000000000e+00 +2.601178739391680563e-06 +-1.000000000000000000e+00 +8.887234042637683843e-05 +-1.000000000000000000e+00 +5.000000000000000000e-01 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +4.877856021170283163e-01 +-1.000000000000000000e+00 +1.769985835774446302e-04 +-1.000000000000000000e+00 +1.670700632870335029e-04 +-1.000000000000000000e+00 +4.899999304038166192e-01 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +4.857384221387814338e-01 +-1.000000000000000000e+00 +2.153456869915780514e-04 +-1.000000000000000000e+00 +7.842896894949250204e-06 +-1.000000000000000000e+00 +6.379065450759913505e-03 +-1.000000000000000000e+00 +6.550426960512012453e-02 +-1.000000000000000000e+00 +6.639379668946672301e-03 +-1.000000000000000000e+00 +7.754700849337902260e-06 +-1.000000000000000000e+00 +2.720743316579670241e-11 +-1.000000000000000000e+00 +5.985108617124988932e-06 +-1.000000000000000000e+00 +3.634644995026691954e-05 +-1.000000000000000000e+00 +3.932389299316284718e-06 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.271128669411294917e-07 +-1.000000000000000000e+00 +8.846493399850662704e-06 +-1.000000000000000000e+00 +2.499931357628383163e-07 +-1.000000000000000000e+00 +8.533030345699353466e-17 +-1.000000000000000000e+00 +1.669925427773410862e-06 +-1.000000000000000000e+00 +1.754138009689378778e-03 +-1.000000000000000000e+00 +1.514897842827161306e-02 +-1.000000000000000000e+00 +1.627979455242757438e-03 +-1.000000000000000000e+00 +2.380955631371225994e-06 +-1.000000000000000000e+00 +2.164891804046736287e-05 +-1.000000000000000000e+00 +8.339317973815890683e-02 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +8.285899874014539257e-02 +-1.000000000000000000e+00 +2.762955748645507000e-05 +-1.000000000000000000e+00 +5.116174051620571477e-05 +-1.000000000000000000e+00 +8.095280568305474045e-02 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +8.636569925236790846e-02 +-1.000000000000000000e+00 +3.452537179033895228e-05 +-1.000000000000000000e+00 +3.598103176641614321e-06 +-1.000000000000000000e+00 +1.570829313580840913e-03 +-1.000000000000000000e+00 +1.621821782442112170e-02 +-1.000000000000000000e+00 +1.544782399547004175e-03 +-1.000000000000000000e+00 +7.989579285134733983e-06 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.594951085377587857e-06 +-1.000000000000000000e+00 +4.323138781337962887e-05 +-1.000000000000000000e+00 +6.948292975998465777e-07 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.446945600849361415e-09 +-1.000000000000000000e+00 +9.041149893307805599e-07 +-1.000000000000000000e+00 +2.206099258062212484e-11 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.681305446488883933e-11 +-1.000000000000000000e+00 +7.091263906557290927e-05 +-1.000000000000000000e+00 +2.746424406288768276e-04 +-1.000000000000000000e+00 +1.285845570921339974e-05 +-1.000000000000000000e+00 +2.029521362171670568e-09 +-1.000000000000000000e+00 +1.291182719850037242e-08 +-1.000000000000000000e+00 +9.127050763987511325e-04 +-1.000000000000000000e+00 +7.385384405891665810e-03 +-1.000000000000000000e+00 +1.001943240150808675e-03 +-1.000000000000000000e+00 +1.080711295768550941e-05 +-1.000000000000000000e+00 +8.020767115181573972e-07 +-1.000000000000000000e+00 +9.252532821720597109e-04 +-1.000000000000000000e+00 +7.188410694198127913e-03 +-1.000000000000000000e+00 +7.245451610090619275e-04 +-1.000000000000000000e+00 +2.515735579020177281e-07 +-1.000000000000000000e+00 +1.482521919913335348e-12 +-1.000000000000000000e+00 +4.756872959630257218e-05 +-1.000000000000000000e+00 +2.668579837809080925e-04 +-1.000000000000000000e+00 +2.822627694453269613e-05 +-1.000000000000000000e+00 +5.182269672663266089e-10 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.936268747135548187e-14 +-1.000000000000000000e+00 +2.544728124173036570e-07 +-1.000000000000000000e+00 +3.898416671385980035e-08 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.639268952045578987e-10 +-1.000000000000000000e+00 +3.705151822008388643e-06 +-1.000000000000000000e+00 +1.133362760482367266e-09 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.229674726063540458e-06 +-1.000000000000000000e+00 +7.962899789920809096e-06 +-1.000000000000000000e+00 +8.211982683649991481e-09 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.563965686625045276e-06 +-1.000000000000000000e+00 +2.248761782893832567e-05 +-1.000000000000000000e+00 +3.279904588337207471e-06 +-1.000000000000000000e+00 +4.076762579823379031e-17 +-1.000000000000000000e+00 +1.547676923757129466e-14 +-1.000000000000000000e+00 +1.030567269874565733e-11 +-1.000000000000000000e+00 +1.957092814065688094e-05 +-1.000000000000000000e+00 +3.852324755037881526e-12 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.769723317139941641e-13 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 diff --git a/tests/regression_tests/weightwindows/ww_p.txt b/tests/regression_tests/weightwindows/ww_p.txt new file mode 100644 index 000000000..09ff01ebd --- /dev/null +++ b/tests/regression_tests/weightwindows/ww_p.txt @@ -0,0 +1,420 @@ +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +5.447493368067179064e-10 +-1.000000000000000000e+00 +4.040462614734400249e-13 +-1.000000000000000000e+00 +9.382800796426366364e-17 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +5.752120832330885919e-11 +-1.000000000000000000e+00 +3.520674728719278012e-05 +-1.000000000000000000e+00 +8.760419353035269667e-05 +-1.000000000000000000e+00 +3.136829920007008993e-05 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +4.233410073487771971e-06 +-1.000000000000000000e+00 +2.551963224311764450e-04 +-1.000000000000000000e+00 +4.777967757842694024e-04 +-1.000000000000000000e+00 +3.598354799974066488e-04 +-1.000000000000000000e+00 +4.625859206051786618e-07 +-1.000000000000000000e+00 +4.587143797469485258e-08 +-1.000000000000000000e+00 +1.291843367759397518e-04 +-1.000000000000000000e+00 +4.695895684632649786e-04 +-1.000000000000000000e+00 +8.862760180332873163e-05 +-1.000000000000000000e+00 +3.279091474401153391e-06 +-1.000000000000000000e+00 +1.151457911727229161e-10 +-1.000000000000000000e+00 +2.295071617064146262e-05 +-1.000000000000000000e+00 +6.922962992428139458e-05 +-1.000000000000000000e+00 +2.160655073217069276e-05 +-1.000000000000000000e+00 +2.617878055106949035e-16 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.722191123879781983e-08 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.123498061876322157e-13 +-1.000000000000000000e+00 +1.462884795971722476e-05 +-1.000000000000000000e+00 +2.777791741571173956e-05 +-1.000000000000000000e+00 +9.533806586698390583e-06 +-1.000000000000000000e+00 +1.078687446476880940e-14 +-1.000000000000000000e+00 +3.057105085570708401e-06 +-1.000000000000000000e+00 +8.231482263927564300e-04 +-1.000000000000000000e+00 +3.881102561422505696e-03 +-1.000000000000000000e+00 +1.229730308180539307e-03 +-1.000000000000000000e+00 +4.684235088236766580e-05 +-1.000000000000000000e+00 +1.675799759004882713e-04 +-1.000000000000000000e+00 +1.172471437684818700e-02 +-1.000000000000000000e+00 +6.827568359008943932e-02 +-1.000000000000000000e+00 +1.144038101286408600e-02 +-1.000000000000000000e+00 +3.782322077390032930e-04 +-1.000000000000000000e+00 +5.098747893974247782e-05 +-1.000000000000000000e+00 +1.208261392046878005e-02 +-1.000000000000000000e+00 +6.428782790571178907e-02 +-1.000000000000000000e+00 +1.129470856680416663e-02 +-1.000000000000000000e+00 +8.168081568291806404e-05 +-1.000000000000000000e+00 +2.347504827875979927e-05 +-1.000000000000000000e+00 +8.787411098754914583e-04 +-1.000000000000000000e+00 +5.151329213225949548e-03 +-1.000000000000000000e+00 +1.102802039050768210e-03 +-1.000000000000000000e+00 +1.088201101868485208e-06 +-1.000000000000000000e+00 +1.058825762372267240e-07 +-1.000000000000000000e+00 +2.704481717610955638e-05 +-1.000000000000000000e+00 +2.621123924679610439e-05 +-1.000000000000000000e+00 +2.070764852974689408e-06 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.380657778513517658e-06 +-1.000000000000000000e+00 +1.007700972669126494e-04 +-1.000000000000000000e+00 +4.313097771132341687e-04 +-1.000000000000000000e+00 +1.335621430311087867e-04 +-1.000000000000000000e+00 +3.792209399652395984e-06 +-1.000000000000000000e+00 +2.205797076162476780e-04 +-1.000000000000000000e+00 +1.884234453311008084e-02 +-1.000000000000000000e+00 +1.151012001392641704e-01 +-1.000000000000000000e+00 +1.825795827020864834e-02 +-1.000000000000000000e+00 +2.084830922016203648e-04 +-1.000000000000000000e+00 +1.333255443748032664e-03 +-1.000000000000000000e+00 +4.931708164307664344e-01 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +4.899952354672477139e-01 +-1.000000000000000000e+00 +1.495703328116801227e-03 +-1.000000000000000000e+00 +1.550036616828060930e-03 +-1.000000000000000000e+00 +4.897735047310072254e-01 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +5.000000000000000000e-01 +-1.000000000000000000e+00 +1.107679609210218599e-03 +-1.000000000000000000e+00 +2.792896409558744781e-04 +-1.000000000000000000e+00 +1.721806295736377085e-02 +-1.000000000000000000e+00 +1.189125124490693353e-01 +-1.000000000000000000e+00 +1.904070707116530328e-02 +-1.000000000000000000e+00 +1.732952862815292309e-04 +-1.000000000000000000e+00 +4.670673837764873805e-08 +-1.000000000000000000e+00 +1.263753469561097525e-04 +-1.000000000000000000e+00 +2.705205429838525473e-04 +-1.000000000000000000e+00 +1.497958445994315019e-04 +-1.000000000000000000e+00 +3.247830168002185438e-08 +-1.000000000000000000e+00 +7.239087311622818554e-08 +-1.000000000000000000e+00 +4.669378458601891341e-04 +-1.000000000000000000e+00 +1.537590300149361015e-03 +-1.000000000000000000e+00 +2.368681880710174145e-04 +-1.000000000000000000e+00 +8.616027294555461632e-07 +-1.000000000000000000e+00 +2.362355168421952122e-04 +-1.000000000000000000e+00 +5.708498338352203244e-02 +-1.000000000000000000e+00 +4.237271368372501623e-01 +-1.000000000000000000e+00 +5.490281661912204542e-02 +-1.000000000000000000e+00 +3.000232638838324392e-04 +-1.000000000000000000e+00 +2.329540159600531311e-03 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.784732930735042255e-03 +-1.000000000000000000e+00 +3.432919201293736875e-03 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +2.672923773363389179e-03 +-1.000000000000000000e+00 +3.182108781955477806e-04 +-1.000000000000000000e+00 +5.639071613735578692e-02 +-1.000000000000000000e+00 +4.346984808472855177e-01 +-1.000000000000000000e+00 +5.500108683629956891e-02 +-1.000000000000000000e+00 +2.417599813822804812e-04 +-1.000000000000000000e+00 +2.774067049392004042e-07 +-1.000000000000000000e+00 +3.257646134837451210e-04 +-1.000000000000000000e+00 +1.264833388398999524e-03 +-1.000000000000000000e+00 +2.981891705231208794e-04 +-1.000000000000000000e+00 +4.986698316296506713e-06 +-1.000000000000000000e+00 +3.419181071650815928e-06 +-1.000000000000000000e+00 +1.081445276486055809e-04 +-1.000000000000000000e+00 +6.934466049804711161e-04 +-1.000000000000000000e+00 +2.510674100806631662e-04 +-1.000000000000000000e+00 +5.515522196186722918e-06 +-1.000000000000000000e+00 +1.424205319552986396e-04 +-1.000000000000000000e+00 +1.854111290797000322e-02 +-1.000000000000000000e+00 +1.236406932635011752e-01 +-1.000000000000000000e+00 +1.760297107458398680e-02 +-1.000000000000000000e+00 +2.130582646223859260e-04 +-1.000000000000000000e+00 +1.501650004537408200e-03 +-1.000000000000000000e+00 +4.903495510314233030e-01 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +4.952432605379828989e-01 +-1.000000000000000000e+00 +1.769327994909421784e-03 +-1.000000000000000000e+00 +9.304497896205120647e-04 +-1.000000000000000000e+00 +4.873897581604286766e-01 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +4.981497862127627907e-01 +-1.000000000000000000e+00 +1.127334097206107765e-03 +-1.000000000000000000e+00 +1.412615360326509664e-04 +-1.000000000000000000e+00 +1.729816946709781048e-02 +-1.000000000000000000e+00 +1.181916566651922268e-01 +-1.000000000000000000e+00 +1.780262484403568810e-02 +-1.000000000000000000e+00 +1.319677450153934042e-04 +-1.000000000000000000e+00 +2.210707325798542070e-07 +-1.000000000000000000e+00 +1.728787338791035720e-04 +-1.000000000000000000e+00 +5.197081842709918220e-04 +-1.000000000000000000e+00 +9.256616392487857544e-05 +-1.000000000000000000e+00 +1.032520710502734699e-08 +-1.000000000000000000e+00 +2.027381220364608057e-16 +-1.000000000000000000e+00 +1.225736024671247218e-05 +-1.000000000000000000e+00 +3.856169633708643398e-05 +-1.000000000000000000e+00 +1.806455099029475317e-05 +-1.000000000000000000e+00 +7.691527781829036254e-17 +-1.000000000000000000e+00 +1.298062961617926916e-05 +-1.000000000000000000e+00 +1.178616843705955963e-03 +-1.000000000000000000e+00 +4.441328419270368887e-03 +-1.000000000000000000e+00 +8.919466185572301909e-04 +-1.000000000000000000e+00 +5.272021975060514019e-06 +-1.000000000000000000e+00 +1.607976851421854571e-04 +-1.000000000000000000e+00 +1.112842419638761758e-02 +-1.000000000000000000e+00 +6.433819434431045647e-02 +-1.000000000000000000e+00 +1.095563748531924904e-02 +-1.000000000000000000e+00 +8.370896020975688591e-05 +-1.000000000000000000e+00 +2.223447159544331337e-04 +-1.000000000000000000e+00 +1.102706140784315975e-02 +-1.000000000000000000e+00 +6.684484191345574366e-02 +-1.000000000000000000e+00 +1.111302464121833623e-02 +-1.000000000000000000e+00 +1.255887947410807431e-04 +-1.000000000000000000e+00 +1.004439072885827748e-05 +-1.000000000000000000e+00 +8.081009012754719698e-04 +-1.000000000000000000e+00 +5.384625561469683769e-03 +-1.000000000000000000e+00 +1.146714816952060971e-03 +-1.000000000000000000e+00 +3.574060778841494232e-06 +-1.000000000000000000e+00 +2.251232708364186507e-10 +-1.000000000000000000e+00 +6.717625881826167172e-05 +-1.000000000000000000e+00 +6.949587369646497293e-05 +-1.000000000000000000e+00 +2.909334597507522605e-05 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +8.582011716859218622e-15 +-1.000000000000000000e+00 +8.304845400451923449e-08 +-1.000000000000000000e+00 +4.575723837205823402e-11 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +1.311680713303003522e-15 +-1.000000000000000000e+00 +2.686689485248190559e-05 +-1.000000000000000000e+00 +4.316605279989737501e-05 +-1.000000000000000000e+00 +6.401055165502384025e-06 +-1.000000000000000000e+00 +5.726152701849347207e-16 +-1.000000000000000000e+00 +8.541606727879512040e-07 +-1.000000000000000000e+00 +2.029396545382813118e-04 +-1.000000000000000000e+00 +4.299195570308034968e-04 +-1.000000000000000000e+00 +8.544159282151438236e-05 +-1.000000000000000000e+00 +1.099532245627738225e-06 +-1.000000000000000000e+00 +1.289814511028014912e-06 +-1.000000000000000000e+00 +1.104204785730725317e-04 +-1.000000000000000000e+00 +6.870228272606287295e-04 +-1.000000000000000000e+00 +1.399807456616496116e-04 +-1.000000000000000000e+00 +4.279460329782799377e-09 +-1.000000000000000000e+00 +1.290717169529484602e-13 +-1.000000000000000000e+00 +1.925724800911436302e-05 +-1.000000000000000000e+00 +9.332371598523186474e-05 +-1.000000000000000000e+00 +1.310633551813380206e-05 +-1.000000000000000000e+00 +6.574706543646946099e-16 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +3.827053602579980514e-15 +-1.000000000000000000e+00 +2.101790478097516960e-09 +-1.000000000000000000e+00 +9.904481358675477728e-13 +-1.000000000000000000e+00 +-1.000000000000000000e+00 diff --git a/tests/unit_tests/weightwindows/__init__.py b/tests/unit_tests/weightwindows/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py new file mode 100644 index 000000000..aaa0bb523 --- /dev/null +++ b/tests/unit_tests/weightwindows/test.py @@ -0,0 +1,194 @@ +import os +from pathlib import Path + +import numpy as np +import pytest +from uncertainties import ufloat + +import openmc +from openmc.stats import Discrete, Point + +from tests import cdtemp + +@pytest.fixture +def model(): + openmc.reset_auto_ids() + model = openmc.Model() + + # materials (M4 steel alloy) + m4 = openmc.Material() + m4.set_density('g/cc', 2.3) + m4.add_nuclide('H1', 0.168018676) + m4.add_nuclide("H2", 1.93244e-05) + m4.add_nuclide("O16", 0.561814465) + m4.add_nuclide("O17", 0.00021401) + m4.add_nuclide("Na23", 0.021365) + m4.add_nuclide("Al27", 0.021343) + m4.add_nuclide("Si28", 0.187439342) + m4.add_nuclide("Si29", 0.009517714) + m4.add_nuclide("Si30", 0.006273944) + m4.add_nuclide("Ca40", 0.018026179) + m4.add_nuclide("Ca42", 0.00012031) + m4.add_nuclide("Ca43", 2.51033e-05) + m4.add_nuclide("Ca44", 0.000387892) + m4.add_nuclide("Ca46", 7.438e-07) + m4.add_nuclide("Ca48", 3.47727e-05) + m4.add_nuclide("Fe54", 0.000248179) + m4.add_nuclide("Fe56", 0.003895875) + m4.add_nuclide("Fe57", 8.99727e-05) + m4.add_nuclide("Fe58", 1.19737e-05) + + s0 = openmc.Sphere(r=240) + s1 = openmc.Sphere(r=250, boundary_type='vacuum') + + c0 = openmc.Cell(fill=m4, region=-s0) + c1 = openmc.Cell(region=+s0 & -s1) + + model.geometry = openmc.Geometry([c0, c1]) + + # settings + settings = model.settings + settings.run_mode = 'fixed source' + settings.particles = 500 + settings.batches = 2 + settings.max_splits = 100 + settings.photon_transport = True + space = Point((0.001, 0.001, 0.001)) + energy = Discrete([14E6], [1.0]) + + settings.source = openmc.Source(space=space, energy=energy) + + # tally + mesh = openmc.RegularMesh() + mesh.lower_left = (-240, -240, -240) + mesh.upper_right = (240, 240, 240) + mesh.dimension = (3, 5, 7) + + mesh_filter = openmc.MeshFilter(mesh) + + e_bnds = [0.0, 0.5, 2E7] + energy_filter = openmc.EnergyFilter(e_bnds) + + particle_filter = openmc.ParticleFilter(['neutron', 'photon']) + + tally = openmc.Tally() + tally.filters = [mesh_filter, energy_filter, particle_filter] + tally.scores = ['flux'] + + model.tallies.append(tally) + + return model + + +def test_weightwindows(model): + + ww_files = ('ww_n.txt', 'ww_p.txt') + cwd = Path(__file__).parent.absolute() + filepaths = [cwd / Path(f) for f in ww_files] + + with cdtemp(filepaths): + # run once with variance reduction off + model.settings.weight_windows_on = False + analog_sp = model.run() + os.rename(analog_sp, 'statepoint.analog.h5') + + # weight windows + + # load pre-generated weight windows + # (created using the same tally as above) + ww_n_lower_bnds = np.loadtxt('ww_n.txt') + ww_p_lower_bnds = np.loadtxt('ww_p.txt') + + # create a mesh matching the one used + # to generate the weight windows + ww_mesh = openmc.RegularMesh() + ww_mesh.lower_left = (-240, -240, -240) + ww_mesh.upper_right = (240, 240, 240) + ww_mesh.dimension = (5, 6, 7) + + # energy bounds matching those of the + # generated weight windows + e_bnds = [0.0, 0.5, 2E7] + + ww_n = openmc.WeightWindows(ww_mesh, + ww_n_lower_bnds, + None, + 10.0, + e_bnds, + survival_ratio=1.01) + + ww_p = openmc.WeightWindows(ww_mesh, + ww_p_lower_bnds, + None, + 10.0, + e_bnds, + survival_ratio=1.01) + + model.settings.weight_windows = [ww_n, ww_p] + + + # run again with variance reduction on + model.settings.weight_windows_on = True + ww_sp = model.run() + os.rename(ww_sp, 'statepoint.ww.h5') + + # load both statepoints and examine results + asp = openmc.StatePoint('statepoint.analog.h5') + wsp = openmc.StatePoint('statepoint.ww.h5') + + analog_tally = asp.tallies[1] + ww_tally = wsp.tallies[1] + + def compare_results(particle, analog_tally, ww_tally): + # get values from each of the tallies + an_mean = analog_tally.get_values(filters=[openmc.ParticleFilter], + filter_bins=[(particle,)]) + ww_mean = ww_tally.get_values(filters=[openmc.ParticleFilter], + filter_bins=[(particle,)]) + + # expect that more bins were scored with weight windows than + # the analog run + assert np.count_nonzero(an_mean) < np.count_nonzero(ww_mean) + + an_rel_err = analog_tally.get_values(filters=[openmc.ParticleFilter], + filter_bins=[(particle,)], + value='rel_err') + ww_rel_err = ww_tally.get_values(filters=[openmc.ParticleFilter], + filter_bins=[(particle,)], + value='rel_err') + + an_rel_err[an_mean == 0.0] = 1.0 + ww_rel_err[ww_mean == 0.0] = 1.0 + + an_avg_rel_err = np.mean(an_rel_err) + ww_avg_rel_err = np.mean(ww_rel_err) + + # expect that the average relative error in the tally + # decreases + assert an_avg_rel_err > ww_avg_rel_err + + # ensure that the value of the mesh bin containing the + # source is statistically similar in both runs + an_std_dev = analog_tally.get_values(filters=[openmc.ParticleFilter], + filter_bins=[(particle,)], + value='std_dev') + ww_std_dev = ww_tally.get_values(filters=[openmc.ParticleFilter], + filter_bins=[(particle,)], + value='std_dev') + + # index of the mesh bin containing the source for the higher + # energy group + source_bin_idx = (an_mean.shape[0]//2, 0, 0) + + an_source_bin = ufloat(an_mean[source_bin_idx], + an_std_dev[source_bin_idx]) + ww_source_bin = ufloat(ww_mean[source_bin_idx], + ww_std_dev[source_bin_idx]) + + diff = an_source_bin - ww_source_bin + + # check that values are within two combined standard deviations + assert abs(diff.nominal_value) / diff.std_dev < 2.0 + + compare_results('neutron', analog_tally, ww_tally) + compare_results('photon', analog_tally, ww_tally) diff --git a/tests/unit_tests/weightwindows/ww_n.txt b/tests/unit_tests/weightwindows/ww_n.txt new file mode 120000 index 000000000..3bd70d2c9 --- /dev/null +++ b/tests/unit_tests/weightwindows/ww_n.txt @@ -0,0 +1 @@ +../../regression_tests/weightwindows/ww_n.txt \ No newline at end of file diff --git a/tests/unit_tests/weightwindows/ww_p.txt b/tests/unit_tests/weightwindows/ww_p.txt new file mode 120000 index 000000000..8a94de709 --- /dev/null +++ b/tests/unit_tests/weightwindows/ww_p.txt @@ -0,0 +1 @@ +../../regression_tests/weightwindows/ww_p.txt \ No newline at end of file