From 62ddd33b5b753305ff89daa9c5d3eb7ec6429884 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 4 Dec 2015 21:39:00 -0500 Subject: [PATCH 01/11] add RPATH information for installed exe --- CMakeLists.txt | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 57a49508b3..5ad537e10a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -235,6 +235,14 @@ if(NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/src/xml/fox/.git) endif() add_subdirectory(src/xml/fox) +#=============================================================================== +# RPATH information +#=============================================================================== + +# add the automatically determined parts of the RPATH +# which point to directories outside the build tree to the install RPATH +set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) + #=============================================================================== # Build OpenMC executable #=============================================================================== From 42f438dcabeedc19a07b97bca5c97d7a7a3fe069 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Mon, 7 Dec 2015 22:06:56 -0500 Subject: [PATCH 02/11] shortened and fixed _align_tally_data routine --- openmc/tallies.py | 80 ++++++++++------------------------------------- 1 file changed, 17 insertions(+), 63 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 13a219dedd..ba5d19dfde 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1625,6 +1625,10 @@ class Tally(object): other_mean = copy.deepcopy(other.mean) other_std_dev = copy.deepcopy(other.std_dev) + # Initialize list of tile and repeat factors + repeat_factors = [1, 1, 1] + tile_factors = [1, 1, 1] + if self.filters != other.filters: # Determine the number of paired combinations of filter bins @@ -1634,83 +1638,33 @@ class Tally(object): # Determine the factors by which each tally operands' data arrays # must be tiled or repeated for the tally outer product - other_tile_factor = 1 - self_repeat_factor = 1 for filter in diff1: - other_tile_factor *= filter.num_bins + tile_factors[0] *= filter.num_bins for filter in diff2: - self_repeat_factor *= filter.num_bins - - # Tile / repeat the tally data for the tally outer product - self_shape = list(self_mean.shape) - other_shape = list(other_mean.shape) - self_shape[0] *= self_repeat_factor - self_mean = np.repeat(self_mean, self_repeat_factor) - self_std_dev = np.repeat(self_std_dev, self_repeat_factor) - - if self_repeat_factor == 1: - other_shape[0] *= other_tile_factor - other_mean = np.repeat(other_mean, other_tile_factor, axis=0) - other_std_dev = np.repeat(other_std_dev, other_tile_factor, - axis=0) - else: - other_mean = np.tile(other_mean, (other_tile_factor, 1, 1)) - other_std_dev = np.tile(other_std_dev, (other_tile_factor, 1, 1)) - - # NumPy repeat and tile routines return 1D flattened arrays - # Reshape arrays as 3D with filters, nuclides and scores axes - self_mean.shape = tuple(self_shape) - self_std_dev.shape = tuple(self_shape) - other_mean.shape = tuple(other_shape) - other_std_dev.shape = tuple(other_shape) + repeat_factors[0] *= filter.num_bins if self.nuclides != other.nuclides: # Determine the number of paired combinations of nuclides # between the two tallies and repeat arrays along nuclide axes - self_repeat_factor = other.num_nuclides - other_tile_factor = self.num_nuclides - - # Tile / repeat the tally data for the tally outer product - self_shape = list(self_mean.shape) - other_shape = list(other_mean.shape) - self_shape[1] *= self_repeat_factor - other_shape[1] *= other_tile_factor - self_mean = np.repeat(self_mean, self_repeat_factor, axis=1) - other_mean = np.tile(other_mean, (1, other_tile_factor, 1)) - self_std_dev = np.repeat(self_std_dev, self_repeat_factor, axis=1) - other_std_dev = np.tile(other_std_dev, (1, other_tile_factor, 1)) - - # NumPy repeat and tile routines return 1D flattened arrays - # Reshape arrays as 3D with filters, nuclides and scores axes - self_mean.shape = tuple(self_shape) - self_std_dev.shape = tuple(self_shape) - other_mean.shape = tuple(other_shape) - other_std_dev.shape = tuple(other_shape) + repeat_factors[1] = other.num_nuclides + tile_factors[1] = self.num_nuclides if self.scores != other.scores: # Determine the number of paired combinations of score bins # between the two tallies and repeat arrays along score axes - self_repeat_factor = other.num_score_bins - other_tile_factor = self.num_score_bins + repeat_factors[2] = other.num_score_bins + tile_factors[2] = self.num_score_bins - # Tile / repeat the tally data for the tally outer product - self_shape = list(self_mean.shape) - other_shape = list(other_mean.shape) - self_shape[2] *= self_repeat_factor - other_shape[2] *= other_tile_factor - self_mean = np.repeat(self_mean, self_repeat_factor, axis=2) - other_mean = np.tile(other_mean, (1, 1, other_tile_factor)) - self_std_dev = np.repeat(self_std_dev, self_repeat_factor, axis=2) - other_std_dev = np.tile(other_std_dev, (1, 1, other_tile_factor)) + # Repeat the self tally + for i in range(3): + self_mean = np.repeat(self_mean, repeat_factors[i], axis=i) + self_std_dev = np.repeat(self_std_dev, repeat_factors[i], axis=i) - # NumPy repeat and tile routines return 1D flattened arrays - # Reshape arrays as 3D with filters, nuclides and scores axes - self_mean.shape = tuple(self_shape) - self_std_dev.shape = tuple(self_shape) - other_mean.shape = tuple(other_shape) - other_std_dev.shape = tuple(other_shape) + # Tile the other tally + other_mean = np.tile(other_mean, tile_factors) + other_std_dev = np.tile(other_std_dev, tile_factors) data = {} data['self'] = {} From 36f633c453205f863ba2881692c3b0eb45d5fdc6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 30 Oct 2015 08:13:20 -0500 Subject: [PATCH 03/11] Incremented version. Added release notes for 0.7.1. --- docs/source/conf.py | 2 +- docs/source/releasenotes.rst | 83 +++++++++++++++++++++++------------- setup.py | 2 +- src/constants.F90 | 2 +- 4 files changed, 56 insertions(+), 33 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 79a7604e39..118ff2c03b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -55,7 +55,7 @@ copyright = u'2011-2015, Massachusetts Institute of Technology' # The short X.Y version. version = "0.7" # The full version, including alpha/beta/rc tags. -release = "0.7.0" +release = "0.7.1" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/docs/source/releasenotes.rst b/docs/source/releasenotes.rst index dffc17201f..7320781e2b 100644 --- a/docs/source/releasenotes.rst +++ b/docs/source/releasenotes.rst @@ -1,9 +1,30 @@ .. _releasenotes: ============================== -Release Notes for OpenMC 0.7.0 +Release Notes for OpenMC 0.7.1 ============================== +This release of OpenMC provides some substantial improvements over version +0.7.0. Non-simple cell regions can now be defined through the ``|`` (union) and +``~`` (complement) operators. Similar changes in the Python API also allow +complex cell regions to be defined. A true secondary particle bank now exists; +this is crucial for photon transport (to be added in the next minor release). A +rich API for multi-group cross section generation has been added via the +``openmc.mgxs`` Python module. + +Various improvements to tallies have also been made. It is now possible to +explicitly specify that a collision estimator be used in a tally. A new +``delayedgroup`` filter and ``delayed-nu-fission`` score allow a user to obtain +delayed fission neutron production rates filtered by delayed group. Finally, the +new ``inverse-velocity`` score may be useful for calculating kinetics +parameters. + +.. caution:: In previous versions, depending on how OpenMC was compiled binary + output was either given in HDF5 or a flat binary format. With this + version, all binary output is now HDF5 which means you **must** + have HDF5 in order to install OpenMC. Please consult the user's + guide for instructions on how to compile with HDF5. + ------------------- System Requirements ------------------- @@ -17,36 +38,41 @@ the problem at hand (mostly on the number of nuclides in the problem). New Features ------------ -- Complete Python API -- Python 3 compatability for all scripts -- All scripts consistently named openmc-* and installed together -- New 'distribcell' tally filter for repeated cells -- Ability to specify outer lattice universe -- XML input validation utility (openmc-validate-xml) -- Support for hexagonal lattices -- Material union energy grid method -- Tally triggers -- Remove dependence on PETSc -- Significant OpenMP performance improvements -- Support for Fortran 2008 MPI interface -- Use of Travis CI for continuous integration -- Simplifications and improvements to test suite +- Support for complex cell regions (union and complement operators) +- Generic quadric surface type +- Improved handling of secondary particles +- Binary output is now solely HDF5 +- ``openmc.mgxs`` Python module enabling multi-group cross section generation +- Collision estimator for tallies +- Delayed fission neutron production tallies with ability to filter by delayed + group +- Inverse velocity tally score +- Performance improvements for binary search +- Performance improvements for reaction rate tallies --------- Bug Fixes --------- -- b5f712_: Fix bug in spherical harmonics tallies -- e6675b_: Ensure all constants are double precision -- 04e2c1_: Fix potential bug in sample_nuclide routine -- 6121d9_: Fix bugs related to particle track files -- 2f0e89_: Fixes for nuclide specification in tallies +- 299322_: Bug with material filter when void material present +- d74840_: Fix triggers on tallies with multiple filters +- c29a81_: Correctly handle maximum transport energy +- 3edc23_: Fixes in the nu-scatter score +- 629e3b_: Assume unspecified surface coefficients are zero in Python API +- 5dbe8b_: Fix energy filters for openmc-plot-mesh-tally +- ff66f4_: Fixes in the openmc-plot-mesh-tally script +- 441fd4_: Fix bug in kappa-fission score +- 7e5974_: Allow fixed source simulations from Python API -.. _b5f712: https://github.com/mit-crpg/openmc/commit/b5f712 -.. _e6675b: https://github.com/mit-crpg/openmc/commit/e6675b -.. _04e2c1: https://github.com/mit-crpg/openmc/commit/04e2c1 -.. _6121d9: https://github.com/mit-crpg/openmc/commit/6121d9 -.. _2f0e89: https://github.com/mit-crpg/openmc/commit/2f0e89 +.. _299322: https://github.com/mit-crpg/openmc/commit/299322 +.. _d74840: https://github.com/mit-crpg/openmc/commit/d74840 +.. _c29a81: https://github.com/mit-crpg/openmc/commit/c29a81 +.. _3edc23: https://github.com/mit-crpg/openmc/commit/3edc23 +.. _629e3b: https://github.com/mit-crpg/openmc/commit/629e3b +.. _5dbe8b: https://github.com/mit-crpg/openmc/commit/5dbe8b +.. _ff66f4: https://github.com/mit-crpg/openmc/commit/ff66f4 +.. _441fd4: https://github.com/mit-crpg/openmc/commit/441fd4 +.. _7e5974: https://github.com/mit-crpg/openmc/commit/7e5974 ------------ Contributors @@ -55,13 +81,10 @@ Contributors This release contains new contributions from the following people: - `Will Boyd `_ -- `Matt Ellis `_ - `Sterling Harper `_ -- `Bryan Herman `_ -- `Nicholas Horelik `_ - `Colin Josey `_ -- `William Lyu `_ - `Adam Nelson `_ - `Paul Romano `_ -- `Anthony Scopatz `_ +- `Kelly Rowland `_ +- `Sam Shaner `_ - `Jon Walsh `_ diff --git a/setup.py b/setup.py index 0c3d5c116f..907d80e31b 100644 --- a/setup.py +++ b/setup.py @@ -10,7 +10,7 @@ except ImportError: have_setuptools = False kwargs = {'name': 'openmc', - 'version': '0.7.0', + 'version': '0.7.1', 'packages': ['openmc', 'openmc.mgxs'], 'scripts': glob.glob('scripts/openmc-*'), diff --git a/src/constants.F90 b/src/constants.F90 index 375c517e76..ba77f35aba 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -8,7 +8,7 @@ module constants ! OpenMC major, minor, and release numbers integer, parameter :: VERSION_MAJOR = 0 integer, parameter :: VERSION_MINOR = 7 - integer, parameter :: VERSION_RELEASE = 0 + integer, parameter :: VERSION_RELEASE = 1 ! Revision numbers for binary files integer, parameter :: REVISION_STATEPOINT = 14 From 6800f8ed219f63871527b8c32450e1ffe3c10e4c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Dec 2015 14:02:10 -0600 Subject: [PATCH 04/11] Update list of publications in documentation --- docs/source/publications.rst | 63 ++++++++++++++++++++++++++++++++---- 1 file changed, 56 insertions(+), 7 deletions(-) diff --git a/docs/source/publications.rst b/docs/source/publications.rst index 79d089605c..369b9d9775 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -26,6 +26,10 @@ Overviews Benchmarking ------------ +- Khurrum S. Chaudri and Sikander M. Mirza, "Burnup dependent Monte Carlo + neutron physics calculations of IAEA MTR benchmark," *Prog. Nucl. Energy*, + **81**, 43-52 (2015). ``_ + - Daniel J. Kelly, Brian N. Aviles, Paul K. Romano, Bryan R. Herman, Nicholas E. Horelik, and Benoit Forget, "Analysis of select BEAVRS PWR benchmark cycle 1 results using MC21 and OpenMC," *Proc. PHYSOR*, Kyoto, @@ -57,13 +61,8 @@ Coupling and Multi-physics - Bryan R. Herman, Benoit Forget, and Kord Smith, "Progress toward Monte Carlo-thermal hydraulic coupling using low-order nonlinear diffusion - acceleration methods." In press, *Ann. Nucl. Energy*, - (2014). ``_ - -- Adam G. Nelson and William R. Martin, "Improved Convergence of Monte Carlo - Generated Multi-Group Scattering Moments," *Proc. Int. Conf. Mathematics and - Computational Methods Applied to Nuclear Science and Engineering*, Sun Valley, - Idaho, May 5--9 (2013). + acceleration methods." *Ann. Nucl. Energy*, **84**, 63-72 + (2015). ``_ - Bryan R. Herman, Benoit Forget, and Kord Smith, "Utilizing CMFD in OpenMC to Estimate Dominance Ratio and Adjoint," *Trans. Am. Nucl. Soc.*, **109**, @@ -81,19 +80,65 @@ Geometry Miscellaneous ------------- +- William Boyd, Sterling Harper, and Paul K. Romano, "Equipping OpenMC for the + big data era," Accepted, *PHYSOR 2016*, Sun Valley, Idaho, May 1-5, 2016. + +- Qicang Shen, William Boyd, Benoit Forget, and Kord Smith, "Tally precision + triggers for the OpenMC Monte Carlo code," *Trans. Am. Nucl. Soc.*, **112**, + 637-640 (2015). + - Timothy P. Burke, Brian C. Kiedrowski, and William R. Martin, "Flux and Reaction Rate Kernel Density Estimators in OpenMC," *Trans. Am. Nucl. Soc.*, **109**, 683-686 (2013). +------------------------------------ +Multi-group Cross Section Generation +------------------------------------ + +- Adam G. Nelson and William R. Martin, "Improved Monte Carlo tallying of + multi-group scattering moments using the NDPP code," *Trans. Am. Nucl. Soc.*, + **113**, 645-648 (2015) + +- Adam G. Nelson and William R. Martin, "Improved Monte Carlo tallying of + multi-group scattering moment matrices," *Trans. Am. Nucl. Soc.*, **110**, + 217-220 (2014). + +- Adam G. Nelson and William R. Martin, "Improved Convergence of Monte Carlo + Generated Multi-Group Scattering Moments," *Proc. Int. Conf. Mathematics and + Computational Methods Applied to Nuclear Science and Engineering*, Sun Valley, + Idaho, May 5--9 (2013). + ------------ Nuclear Data ------------ +- Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed multipole + for cross section Doppler broadening," *J. Comput. Phys.*, In Press + (2016). ``_ + +- Colin Josey, Benoit Forget, and Kord Smith, "Windowed multipole sensitivity to + target accuracy of the optimization procedure," *J. Nucl. Sci. Technol.*, + **52**, 987-992 (2015). ``_ + +- Jonathan A. Walsh, Paul K. Romano, Benoit Forget, and Kord S. Smith, + "Optimizations of the energy grid search algorithm in continuous-energy Monte + Carlo particle transport codes", *Comput. Phys. Commun.*, **196**, 134-142 + (2015). ``_ + - Jonathan A. Walsh, Benoit Forget, Kord S. Smith, Brian C. Kiedrowski, and Forrest B. Brown, "Direct, on-the-fly calculation of unresolved resonance region cross sections in Monte Carlo simulations," *Proc. Joint Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015). +- Amanda L. Lund, Andrew R. Siegel, Benoit Forget, Colin Josey, and + Paul K. Romano, "Using fractional cascading to accelerate cross section + lookups in Monte Carlo particle transport calculations," *Proc. Joint + Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015). + +- Ronald O. Rahaman, Andrew R. Siegel, and Paul K. Romano, "Monte Carlo + performance analysis for varying cross section parameter regimes," + *Proc. Joint Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015). + - Paul K. Romano and Timothy H. Trumbull, "Comparison of algorithms for Doppler broadening pointwise tabulated cross sections," *Ann. Nucl. Energy*, **75**, 358--364 (2015). ``_ @@ -114,6 +159,10 @@ Nuclear Data Parallelism ----------- +- Paul K. Romano, John R. Tramm, and Andrew R. Siegel, "Efficacy of hardware + threading for Monte Carlo particle transport calculations on multi- and + many-core systems," Accepted, *PHYSOR 2016*, Sun Valley, Idaho, May 1-5, 2016. + - David Ozog, Allen D. Malony, and Andrew R. Siegel, "A performance analysis of SIMD algorithms for Monte Carlo simulations of nuclear reactor cores," *Proc. IEEE Int. Parallel and Distributed Processing Symposium*, Hyderabad, From 9cde4ce12c10b6acacacbee310cd08f19623d955 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 7 Dec 2015 22:16:15 -0600 Subject: [PATCH 05/11] Clarify installation instructions for HDF5 pre-req and Ubuntu PPA --- docs/source/usersguide/install.rst | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index dcf990bdad..e3f0df5e91 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -8,7 +8,7 @@ Installation and Configuration Installing on Ubuntu with PPA ----------------------------- -For users with Ubuntu 11.10 or later, a binary package for OpenMC is available +For users with Ubuntu 15.04 or later, a binary package for OpenMC is available through a Personal Package Archive (PPA) and can be installed through the APT package manager. First, add the following PPA to the repository sources: @@ -28,6 +28,9 @@ Now OpenMC should be recognized within the repository and can be installed: sudo apt-get install openmc +Binary packages from this PPA may exist for earlier versions of Ubuntu, but they +are no longer supported. + -------------------- Building from Source -------------------- @@ -74,6 +77,12 @@ Prerequisites You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial. + .. important:: + + OpenMC uses various parts of the HDF5 Fortran 2003 API; as such you + must include ``--enable-fortran2003`` or else OpenMC will not be able + to compile. + On Debian derivatives, HDF5 and/or parallel HDF5 can be installed through the APT package manager: From f84a4ee4309d76200981c92896f28e8fa1bb8ad3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 8 Dec 2015 07:19:31 -0600 Subject: [PATCH 06/11] Avoid DeprecationWarning from xml.etree.Element.getchildren() --- openmc/clean_xml.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/clean_xml.py b/openmc/clean_xml.py index aefd30ac7b..564281a5cc 100644 --- a/openmc/clean_xml.py +++ b/openmc/clean_xml.py @@ -1,7 +1,7 @@ def sort_xml_elements(tree): # Retrieve all children of the root XML node in the tree - elements = tree.getchildren() + elements = list(tree) # Initialize empty lists for the sorted and comment elements sorted_elements = [] From b17a8e07e3ae085653b4c57ff2bee6e941dc28e0 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 8 Dec 2015 13:19:02 -0500 Subject: [PATCH 07/11] added test for tally arithmetic --- tests/test_tally_arithmetic/geometry.xml | 8 ++ tests/test_tally_arithmetic/materials.xml | 11 +++ tests/test_tally_arithmetic/results_true.dat | 49 ++++++++++ tests/test_tally_arithmetic/settings.xml | 18 ++++ .../test_tally_arithmetic.py | 91 +++++++++++++++++++ 5 files changed, 177 insertions(+) create mode 100644 tests/test_tally_arithmetic/geometry.xml create mode 100644 tests/test_tally_arithmetic/materials.xml create mode 100644 tests/test_tally_arithmetic/results_true.dat create mode 100644 tests/test_tally_arithmetic/settings.xml create mode 100644 tests/test_tally_arithmetic/test_tally_arithmetic.py diff --git a/tests/test_tally_arithmetic/geometry.xml b/tests/test_tally_arithmetic/geometry.xml new file mode 100644 index 0000000000..bc56030e18 --- /dev/null +++ b/tests/test_tally_arithmetic/geometry.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/tests/test_tally_arithmetic/materials.xml b/tests/test_tally_arithmetic/materials.xml new file mode 100644 index 0000000000..e7947a92da --- /dev/null +++ b/tests/test_tally_arithmetic/materials.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/tests/test_tally_arithmetic/results_true.dat b/tests/test_tally_arithmetic/results_true.dat new file mode 100644 index 0000000000..85183e3e8f --- /dev/null +++ b/tests/test_tally_arithmetic/results_true.dat @@ -0,0 +1,49 @@ +Tally + ID = 10000 + Name = (tally 1 + tally 2) + Filters = + energy [ 0. 20.] + (cell + material) (array([1]), array([1])) + Nuclides = (U-235 + U-235) (U-235 + Pu-239) (U-238 + U-235) (U-238 + Pu-239) + Scores = [(fission + fission), (fission + absorption), (nu-fission + fission), (nu-fission + absorption)] + Estimator = tracklength +[[[ 0.07510122 0.07839105 0.13713377 0.1404236 ] + [ 0.0916553 0.09321683 0.15368785 0.15524938] + [ 0.04596277 0.0492526 0.06126275 0.06455259] + [ 0.06251685 0.06407838 0.07781683 0.07937836]]]Tally + ID = 10001 + Name = (tally 1 - tally 2) + Filters = + energy [ 0. 20.] + (cell - material) (array([1]), array([1])) + Nuclides = (U-235 - U-235) (U-235 - Pu-239) (U-238 - U-235) (U-238 - Pu-239) + Scores = [(fission - fission), (fission - absorption), (nu-fission - fission), (nu-fission - absorption)] + Estimator = tracklength +[[[ 0. -0.00328983 0.06203255 0.05874271] + [-0.01655408 -0.01811561 0.04547847 0.04391694] + [-0.02913845 -0.03242829 -0.01383847 -0.0171283 ] + [-0.04569253 -0.04725406 -0.03039255 -0.03195408]]]Tally + ID = 10002 + Name = (tally 1 * tally 2) + Filters = + energy [ 0. 20.] + (cell * material) (array([1]), array([1])) + Nuclides = (U-235 * U-235) (U-235 * Pu-239) (U-238 * U-235) (U-238 * Pu-239) + Scores = [(fission * fission), (fission * absorption), (nu-fission * fission), (nu-fission * absorption)] + Estimator = tracklength +[[[ 0.00141005 0.00153358 0.00373941 0.00406702] + [ 0.00203166 0.0020903 0.00538792 0.00554342] + [ 0.00031588 0.00034356 0.00089041 0.00096841] + [ 0.00045514 0.00046827 0.00128294 0.00131997]]]Tally + ID = 10003 + Name = (tally 1 / tally 2) + Filters = + energy [ 0. 20.] + (cell / material) (array([1]), array([1])) + Nuclides = (U-235 / U-235) (U-235 / Pu-239) (U-238 / U-235) (U-238 / Pu-239) + Scores = [(fission / fission), (fission / absorption), (nu-fission / fission), (nu-fission / absorption)] + Estimator = tracklength +[[[ 1. 0.91944666 2.65197177 2.4383466 ] + [ 0.69403611 0.67456727 1.84056418 1.78893335] + [ 0.22402189 0.20597618 0.63147153 0.58060439] + [ 0.15547928 0.15111784 0.43826405 0.42597002]]] \ No newline at end of file diff --git a/tests/test_tally_arithmetic/settings.xml b/tests/test_tally_arithmetic/settings.xml new file mode 100644 index 0000000000..a69dde686a --- /dev/null +++ b/tests/test_tally_arithmetic/settings.xml @@ -0,0 +1,18 @@ + + + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + + + diff --git a/tests/test_tally_arithmetic/test_tally_arithmetic.py b/tests/test_tally_arithmetic/test_tally_arithmetic.py new file mode 100644 index 0000000000..ee1a7a6481 --- /dev/null +++ b/tests/test_tally_arithmetic/test_tally_arithmetic.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python + +import os +import sys +import glob +import hashlib +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness +import openmc + + +class TallyArithmeticTestHarness(TestHarness): + def _build_inputs(self): + + u235 = openmc.Nuclide('U-235') + u238 = openmc.Nuclide('U-238') + pu239 = openmc.Nuclide('Pu-239') + + # Instantiate energy filter + energy_filter = openmc.Filter(type='energy', bins=[0., 20.]) + + # Create tallies + tally_1 = openmc.Tally(name='tally 1') + tally_1.add_filter(openmc.Filter(type='cell', bins=[1])) + tally_1.add_filter(energy_filter) + tally_1.add_score('fission') + tally_1.add_score('nu-fission') + tally_1.add_nuclide(u235) + tally_1.add_nuclide(u238) + + tally_2 = openmc.Tally(name='tally 2') + tally_2.add_filter(openmc.Filter(type='material', bins=[1])) + tally_2.add_filter(energy_filter) + tally_2.add_score('fission') + tally_2.add_score('absorption') + tally_2.add_nuclide(u235) + tally_2.add_nuclide(pu239) + + # Export tallies to file + tallies_file = openmc.TalliesFile() + tallies_file.add_tally(tally_1) + tallies_file.add_tally(tally_2) + tallies_file.export_to_xml() + + def _get_results(self, hash_output=False): + """Digest info in the statepoint and return as a string.""" + + # Read the statepoint file. + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] + sp = openmc.StatePoint(statepoint) + + # Read the summary file. + summary = glob.glob(os.path.join(os.getcwd(), 'summary.h5'))[0] + su = openmc.Summary(summary) + sp.link_with_summary(su) + + # Load the tallies + tally_1 = sp.get_tally(name='tally 1') + tally_2 = sp.get_tally(name='tally 2') + + # Perform all the tally arithmetic operations and output results + outstr = '' + tally_3 = tally_1 + tally_2 + outstr += tally_3.__repr__() + outstr += str(tally_3.mean) + + tally_3 = tally_1 - tally_2 + outstr += tally_3.__repr__() + outstr += str(tally_3.mean) + + tally_3 = tally_1 * tally_2 + outstr += tally_3.__repr__() + outstr += str(tally_3.mean) + + tally_3 = tally_1 / tally_2 + outstr += tally_3.__repr__() + outstr += str(tally_3.mean) + + print(outstr) + + # Hash the results if necessary + if hash_output: + sha512 = hashlib.sha512() + sha512.update(outstr.encode('utf-8')) + outstr = sha512.hexdigest() + + return outstr + +if __name__ == '__main__': + harness = TallyArithmeticTestHarness('statepoint.10.*', True) + harness.main() From b81fa69db3907d129163eff8a6b83dc8d7b38761 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 8 Dec 2015 13:27:45 -0500 Subject: [PATCH 08/11] added tallies.xml file to tally arithmetic test --- tests/test_tally_arithmetic/tallies.xml | 15 +++++++++ .../test_tally_arithmetic.py | 32 ------------------- 2 files changed, 15 insertions(+), 32 deletions(-) create mode 100644 tests/test_tally_arithmetic/tallies.xml diff --git a/tests/test_tally_arithmetic/tallies.xml b/tests/test_tally_arithmetic/tallies.xml new file mode 100644 index 0000000000..2a2eb48361 --- /dev/null +++ b/tests/test_tally_arithmetic/tallies.xml @@ -0,0 +1,15 @@ + + + + + + U-235 U-238 + fission nu-fission + + + + + U-235 Pu-239 + fission absorption + + diff --git a/tests/test_tally_arithmetic/test_tally_arithmetic.py b/tests/test_tally_arithmetic/test_tally_arithmetic.py index ee1a7a6481..f71c960e92 100644 --- a/tests/test_tally_arithmetic/test_tally_arithmetic.py +++ b/tests/test_tally_arithmetic/test_tally_arithmetic.py @@ -10,38 +10,6 @@ import openmc class TallyArithmeticTestHarness(TestHarness): - def _build_inputs(self): - - u235 = openmc.Nuclide('U-235') - u238 = openmc.Nuclide('U-238') - pu239 = openmc.Nuclide('Pu-239') - - # Instantiate energy filter - energy_filter = openmc.Filter(type='energy', bins=[0., 20.]) - - # Create tallies - tally_1 = openmc.Tally(name='tally 1') - tally_1.add_filter(openmc.Filter(type='cell', bins=[1])) - tally_1.add_filter(energy_filter) - tally_1.add_score('fission') - tally_1.add_score('nu-fission') - tally_1.add_nuclide(u235) - tally_1.add_nuclide(u238) - - tally_2 = openmc.Tally(name='tally 2') - tally_2.add_filter(openmc.Filter(type='material', bins=[1])) - tally_2.add_filter(energy_filter) - tally_2.add_score('fission') - tally_2.add_score('absorption') - tally_2.add_nuclide(u235) - tally_2.add_nuclide(pu239) - - # Export tallies to file - tallies_file = openmc.TalliesFile() - tallies_file.add_tally(tally_1) - tallies_file.add_tally(tally_2) - tallies_file.export_to_xml() - def _get_results(self, hash_output=False): """Digest info in the statepoint and return as a string.""" From 4888662b9cf543d31203ddf82961055442b1b9f5 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 9 Dec 2015 10:41:19 -0500 Subject: [PATCH 09/11] fixed tally arithmetic to work with unique filters in both tally operands --- openmc/mgxs/mgxs.py | 12 +- openmc/tallies.py | 152 +++++++++--------- tests/test_tally_arithmetic/results_true.dat | 12 +- .../test_tally_arithmetic.py | 10 +- 4 files changed, 98 insertions(+), 88 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 635c822e31..597e28e292 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -47,7 +47,7 @@ _DOMAINS = [openmc.Cell, class MGXS(object): - """An abstract multi-group cross section for some energy group structure + """An abstract multi-group cross section for some energy group structure within some spatial domain. This class can be used for both OpenMC input generation and tally data @@ -95,7 +95,7 @@ class MGXS(object): num_sumbdomains : Integral The number of subdomains is unity for 'material', 'cell' and 'universe' domain types. When the This is equal to the number of cell instances - for 'distribcell' domain types (it is equal to unity prior to loading + for 'distribcell' domain types (it is equal to unity prior to loading tally data from a statepoint file). num_nuclides : Integral The number of nuclides for which the multi-group cross section is @@ -1420,8 +1420,8 @@ class AbsorptionXS(MGXS): class CaptureXS(MGXS): """A capture multi-group cross section. - - The Neutron capture reaction rate is defined as the difference between + + The Neutron capture reaction rate is defined as the difference between OpenMC's 'absorption' and 'fission' reaction rate score types. This includes not only radiative capture, but all forms of neutron disappearance aside from fission (e.g., MT > 100). @@ -1723,8 +1723,8 @@ class ScatterMatrixXS(MGXS): if self.correction == 'P0': scatter_p1 = self.tallies['scatter-P1'] scatter_p1 = scatter_p1.get_slice(scores=['scatter-P1']) - energy_filter = openmc.Filter(type='energy') - energy_filter.bins = self.energy_groups.group_edges + energy_filter = copy.deepcopy(self.tallies['scatter']. + find_filter('energy')) scatter_p1 = scatter_p1.diagonalize_filter(energy_filter) rxn_tally = self.tallies['scatter'] - scatter_p1 else: diff --git a/openmc/tallies.py b/openmc/tallies.py index ba5d19dfde..3eec325ab9 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1200,12 +1200,12 @@ class Tally(object): # Tile the nuclide bins into a DataFrame column nuclides = np.repeat(nuclides, len(self.scores)) tile_factor = data_size / len(nuclides) - df['nuclide'] = np.tile(nuclides, tile_factor) + df['nuclide'] = np.tile(nuclides, int(tile_factor)) # Include column for scores if user requested it if scores: tile_factor = data_size / len(self.scores) - df['score'] = np.tile(self.scores, tile_factor) + df['score'] = np.tile(self.scores, int(tile_factor)) # Append columns with mean, std. dev. for each tally bin df['mean'] = self.mean.ravel() @@ -1425,25 +1425,30 @@ class Tally(object): # Pickle the Tally results to a file pickle.dump(tally_results, open(filename, 'wb')) - def _outer_product(self, other, binary_op): + def _hybrid_product(self, other, binary_op): """Combines filters, scores and nuclides with another tally. - This is a helper method for the tally arithmetic methods. The filters, - scores and nuclides from both tallies are enumerated into all possible - combinations and expressed as CrossFilter, CrossScore and - CrossNuclide objects in the new derived tally. + This is a helper method for the tally arithmetic methods. It is called a + "hybrid product" because it performs a combination of a tensor + (or Kronecker) product and entrywise (or Hadamard) product. The filters, + nuclides, and scores from both tallies are combined using an entrywise + (or Hadamard) product on matching filters. If all nuclides are identical + in the two tallies, the entrywise product is performed across nuclides; + else the tensor product is performed. If all scores are identical in the + two tallies, the entrywise product is performed across scores; else the + tensor product is performed. Parameters ---------- other : Tally - The tally on the right hand side of the outer product + The tally on the right hand side of the hybrid product binary_op : {'+', '-', '*', '/', '^'} - The binary operation in the outer product + The binary operation in the hybrid product Returns ------- Tally - A new Tally that is the outer product with this one. + A new Tally that is the hybrid product with this one. Raises ------ @@ -1473,25 +1478,6 @@ class Tally(object): self_copy = copy.deepcopy(self) other_copy = copy.deepcopy(other) - # Find any shared filters between the two tallies - filter_intersect = [] - for filter in self_copy.filters: - if filter in other_copy.filters: - filter_intersect.append(filter) - - # Align the shared filters in successive order - for i, filter in enumerate(filter_intersect): - self_index = self_copy.filters.index(filter) - other_index = other_copy.filters.index(filter) - - # If necessary, swap self filter - if self_index != i: - self_copy.swap_filters(filter, self_copy.filters[i], inplace=True) - - # If necessary, swap other filter - if other_index != i: - other_copy.swap_filters(filter, other_copy.filters[i], inplace=True) - data = self_copy._align_tally_data(other_copy) if binary_op == '+': @@ -1535,7 +1521,7 @@ class Tally(object): for self_filter in self_copy.filters: new_tally.add_filter(self_filter) - # Generate filter "outer products" for non-identical filters + # Generate filter entrywise product for non-identical filters else: # Find the common longest sequence of shared filters @@ -1600,7 +1586,7 @@ class Tally(object): This is a helper method to construct a dict of dicts of the "aligned" data arrays from each tally for tally arithmetic. The method analyzes - the filters, scores and nuclides in both tally's and determines how to + the filters, scores and nuclides in both tallies and determines how to appropriately align the data for vectorized arithmetic. For example, if the two tallies have different filters, this method will use NumPy 'tile' and 'repeat' operations to the new data arrays such that all @@ -1620,51 +1606,71 @@ class Tally(object): """ + # Get the set of filters that each tally is missing + other_missing_filters = set(self.filters).difference(set(other.filters)) + self_missing_filters = set(other.filters).difference(set(self.filters)) + + # Add other_missing_filters to other + for filter in other_missing_filters: + filter = copy.deepcopy(filter) + repeat_factor = filter.num_bins + other._mean = np.repeat(other.mean, repeat_factor, axis=0) + other.sum = np.repeat(other.sum, repeat_factor, axis=0) + other._std_dev = np.repeat(other.std_dev, repeat_factor, axis=0) + other.sum_sq = np.repeat(other.sum_sq, repeat_factor, axis=0) + other.add_filter(filter) + + # Correct the stride for other filters + stride = other.num_nuclides * other.num_score_bins + for filter in reversed(other.filters): + filter.stride = stride + stride *= filter.num_bins + + # Add self_missing_filters to self + for filter in self_missing_filters: + filter = copy.deepcopy(filter) + repeat_factor = filter.num_bins + self._mean = np.repeat(self.mean, repeat_factor, axis=0) + self.sum = np.repeat(self.sum, repeat_factor, axis=0) + self._std_dev = np.repeat(self.std_dev, repeat_factor, axis=0) + self.sum_sq = np.repeat(self.sum_sq, repeat_factor, axis=0) + self.add_filter(filter) + + # Correct the stride for self filters + stride = self.num_nuclides * self.num_score_bins + for filter in reversed(self.filters): + filter.stride = stride + stride *= filter.num_bins + + # Align other filters with self filters + for i, filter in enumerate(self.filters): + other_index = other.filters.index(filter) + + # If necessary, swap other filter + if other_index != i: + other.swap_filters(filter, other.filters[i], inplace=True) + + # Deep copy the mean and std dev data self_mean = copy.deepcopy(self.mean) self_std_dev = copy.deepcopy(self.std_dev) other_mean = copy.deepcopy(other.mean) other_std_dev = copy.deepcopy(other.std_dev) - # Initialize list of tile and repeat factors - repeat_factors = [1, 1, 1] - tile_factors = [1, 1, 1] - - if self.filters != other.filters: - - # Determine the number of paired combinations of filter bins - # between the two tallies and repeat arrays along filter axes - diff1 = list(set(self.filters).difference(set(other.filters))) - diff2 = list(set(other.filters).difference(set(self.filters))) - - # Determine the factors by which each tally operands' data arrays - # must be tiled or repeated for the tally outer product - for filter in diff1: - tile_factors[0] *= filter.num_bins - for filter in diff2: - repeat_factors[0] *= filter.num_bins - + # If the tallies do not have identical sets of nuclides, tile and repeat + # to perform cross product of data for each nuclide if self.nuclides != other.nuclides: + self_mean = np.repeat(self_mean, other.num_nuclides, axis=1) + self_std_dev = np.repeat(self_std_dev, other.num_nuclides, axis=1) + other_mean = np.tile(other_mean, (1, self.num_nuclides, 1)) + other_std_dev = np.tile(other_std_dev, (1, self.num_nuclides, 1)) - # Determine the number of paired combinations of nuclides - # between the two tallies and repeat arrays along nuclide axes - repeat_factors[1] = other.num_nuclides - tile_factors[1] = self.num_nuclides - + # If the tallies do not have identical sets of scores, tile and repeat + # to perform cross product of data for each score if self.scores != other.scores: - - # Determine the number of paired combinations of score bins - # between the two tallies and repeat arrays along score axes - repeat_factors[2] = other.num_score_bins - tile_factors[2] = self.num_score_bins - - # Repeat the self tally - for i in range(3): - self_mean = np.repeat(self_mean, repeat_factors[i], axis=i) - self_std_dev = np.repeat(self_std_dev, repeat_factors[i], axis=i) - - # Tile the other tally - other_mean = np.tile(other_mean, tile_factors) - other_std_dev = np.tile(other_std_dev, tile_factors) + self_mean = np.repeat(self_mean, other.num_score_bins, axis=2) + self_std_dev = np.repeat(self_std_dev, other.num_score_bins, axis=2) + other_mean = np.tile(other_mean, (1, 1, self.num_score_bins)) + other_std_dev = np.tile(other_std_dev, (1, 1, self.num_score_bins)) data = {} data['self'] = {} @@ -1845,7 +1851,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='+') + new_tally = self._hybrid_product(other, binary_op='+') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -1915,7 +1921,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='-') + new_tally = self._hybrid_product(other, binary_op='-') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -1985,7 +1991,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='*') + new_tally = self._hybrid_product(other, binary_op='*') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -2055,7 +2061,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='/') + new_tally = self._hybrid_product(other, binary_op='/') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -2128,7 +2134,7 @@ class Tally(object): raise ValueError(msg) if isinstance(power, Tally): - new_tally = self._outer_product(power, binary_op='^') + new_tally = self._hybrid_product(power, binary_op='^') elif isinstance(power, Real): new_tally = Tally(name='derived') diff --git a/tests/test_tally_arithmetic/results_true.dat b/tests/test_tally_arithmetic/results_true.dat index 85183e3e8f..23dd38cfb6 100644 --- a/tests/test_tally_arithmetic/results_true.dat +++ b/tests/test_tally_arithmetic/results_true.dat @@ -2,8 +2,9 @@ Tally ID = 10000 Name = (tally 1 + tally 2) Filters = + cell [1] energy [ 0. 20.] - (cell + material) (array([1]), array([1])) + material [1] Nuclides = (U-235 + U-235) (U-235 + Pu-239) (U-238 + U-235) (U-238 + Pu-239) Scores = [(fission + fission), (fission + absorption), (nu-fission + fission), (nu-fission + absorption)] Estimator = tracklength @@ -14,8 +15,9 @@ Tally ID = 10001 Name = (tally 1 - tally 2) Filters = + cell [1] energy [ 0. 20.] - (cell - material) (array([1]), array([1])) + material [1] Nuclides = (U-235 - U-235) (U-235 - Pu-239) (U-238 - U-235) (U-238 - Pu-239) Scores = [(fission - fission), (fission - absorption), (nu-fission - fission), (nu-fission - absorption)] Estimator = tracklength @@ -26,8 +28,9 @@ Tally ID = 10002 Name = (tally 1 * tally 2) Filters = + cell [1] energy [ 0. 20.] - (cell * material) (array([1]), array([1])) + material [1] Nuclides = (U-235 * U-235) (U-235 * Pu-239) (U-238 * U-235) (U-238 * Pu-239) Scores = [(fission * fission), (fission * absorption), (nu-fission * fission), (nu-fission * absorption)] Estimator = tracklength @@ -38,8 +41,9 @@ Tally ID = 10003 Name = (tally 1 / tally 2) Filters = + cell [1] energy [ 0. 20.] - (cell / material) (array([1]), array([1])) + material [1] Nuclides = (U-235 / U-235) (U-235 / Pu-239) (U-238 / U-235) (U-238 / Pu-239) Scores = [(fission / fission), (fission / absorption), (nu-fission / fission), (nu-fission / absorption)] Estimator = tracklength diff --git a/tests/test_tally_arithmetic/test_tally_arithmetic.py b/tests/test_tally_arithmetic/test_tally_arithmetic.py index f71c960e92..2f49061f96 100644 --- a/tests/test_tally_arithmetic/test_tally_arithmetic.py +++ b/tests/test_tally_arithmetic/test_tally_arithmetic.py @@ -29,19 +29,19 @@ class TallyArithmeticTestHarness(TestHarness): # Perform all the tally arithmetic operations and output results outstr = '' tally_3 = tally_1 + tally_2 - outstr += tally_3.__repr__() + outstr += repr(tally_3) outstr += str(tally_3.mean) tally_3 = tally_1 - tally_2 - outstr += tally_3.__repr__() + outstr += repr(tally_3) outstr += str(tally_3.mean) tally_3 = tally_1 * tally_2 - outstr += tally_3.__repr__() + outstr += repr(tally_3) outstr += str(tally_3.mean) tally_3 = tally_1 / tally_2 - outstr += tally_3.__repr__() + outstr += repr(tally_3) outstr += str(tally_3.mean) print(outstr) @@ -55,5 +55,5 @@ class TallyArithmeticTestHarness(TestHarness): return outstr if __name__ == '__main__': - harness = TallyArithmeticTestHarness('statepoint.10.*', True) + harness = TallyArithmeticTestHarness('statepoint.10.h5', True) harness.main() From 4d3046a5976998b7371ebd5de856213e2972333a Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 10 Dec 2015 00:53:01 -0500 Subject: [PATCH 10/11] added capability to perform tensor or entrywise products across nuclides and scores in tally arithmetic --- openmc/tallies.py | 395 +++++++++++++++++++++++++++++++++++----------- 1 file changed, 303 insertions(+), 92 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 3eec325ab9..5656987cbf 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -24,6 +24,7 @@ if sys.version_info[0] >= 3: # "Static" variable for auto-generated Tally IDs AUTO_TALLY_ID = 10000 +_PRODUCT_TYPES = ['tensor', 'entrywise'] def reset_auto_tally_id(): global AUTO_TALLY_ID @@ -1425,18 +1426,20 @@ class Tally(object): # Pickle the Tally results to a file pickle.dump(tally_results, open(filename, 'wb')) - def _hybrid_product(self, other, binary_op): + def hybrid_product(self, other, binary_op, filter_product='None', + nuclide_product='None', score_product='None'): """Combines filters, scores and nuclides with another tally. This is a helper method for the tally arithmetic methods. It is called a - "hybrid product" because it performs a combination of a tensor - (or Kronecker) product and entrywise (or Hadamard) product. The filters, + "hybrid product" because it performs a combination of tensor + (or Kronecker) and entrywise (or Hadamard) products. The filters, nuclides, and scores from both tallies are combined using an entrywise - (or Hadamard) product on matching filters. If all nuclides are identical - in the two tallies, the entrywise product is performed across nuclides; - else the tensor product is performed. If all scores are identical in the - two tallies, the entrywise product is performed across scores; else the - tensor product is performed. + (or Hadamard) product on matching filters. By default, if all nuclides + are identical in the two tallies, the entrywise product is performed + across nuclides; else the tensor product is performed. By default, if all + scores are identical in the two tallies, the entrywise product is + performed across scores; else the tensor product is performed. Users can + also call the method explicitly and specify the desired product. Parameters ---------- @@ -1444,6 +1447,21 @@ class Tally(object): The tally on the right hand side of the hybrid product binary_op : {'+', '-', '*', '/', '^'} The binary operation in the hybrid product + filter_product : str, optional + The type of product (tensor or entrywise) to be performed between + filter data. The default is the entrywise product. Currently only + the entrywise product is supported since a tally cannot contain + two of the same tallies. + nuclide_product : str, optional + The type of product (tensor or entrywise) to be performed between + nuclide data. The default is the entrywise product if all nuclides + between the two tallies are the same; otherwise the default is + the tensor product. + score_product : str, optional + The type of product (tensor or entrywise) to be performed between + score data. The default is the entrywise product if all scores + between the two tallies are the same; otherwise the default is + the tensor product. Returns ------- @@ -1458,6 +1476,33 @@ class Tally(object): """ + # Set default value for filter product if it was not set + if filter_product == 'None': + filter_product = 'entrywise' + elif filter_product == 'tensor': + msg = 'Unable to perform Tally arithmetic with a tensor product' \ + 'for the filter data as this not currently supported.' + raise ValueError(msg) + + # Set default value for nuclide product if it was not set + if nuclide_product == 'None': + if self.nuclides == other.nuclides: + nuclide_product = 'entrywise' + else: + nuclide_product = 'tensor' + + # Set default value for score product if it was not set + if score_product == 'None': + if self.scores == other.scores: + score_product = 'entrywise' + else: + score_product = 'tensor' + + # Check product types + cv.check_value('filter product', filter_product, _PRODUCT_TYPES) + cv.check_value('nuclide product', nuclide_product, _PRODUCT_TYPES) + cv.check_value('score product', score_product, _PRODUCT_TYPES) + # Check that results have been read if not other.derived and other.sum is None: msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ @@ -1473,13 +1518,23 @@ class Tally(object): new_name = '({0} {1} {2})'.format(self.name, binary_op, other.name) new_tally.name = new_name + # Query the mean and std dev so the tally data is read in from file + # if it has not present + self.mean + self.std_dev + other.mean + other.std_dev + # Create copies of self and other tallies to rearrange for tally # arithmetic self_copy = copy.deepcopy(self) other_copy = copy.deepcopy(other) - data = self_copy._align_tally_data(other_copy) + # Align the tally data based on desired hybrid product + data = self_copy._align_tally_data(other_copy, filter_product, + nuclide_product, score_product) + # Perform tally arithmetic operation if binary_op == '+': new_tally._mean = data['self']['mean'] + data['other']['mean'] new_tally._std_dev = np.sqrt(data['self']['std. dev.']**2 + @@ -1509,6 +1564,13 @@ class Tally(object): new_tally._std_dev = np.abs(new_tally.mean) * \ np.sqrt(first_term**2 + second_term**2) + # Convert any infs and nans to zero + new_tally._mean[np.isinf(new_tally._mean)] = 0 + new_tally._mean = np.nan_to_num(new_tally._mean) + new_tally._std_dev[np.isinf(new_tally._std_dev)] = 0 + new_tally._std_dev = np.nan_to_num(new_tally._std_dev) + + # Set tally attributes if self_copy.estimator == other_copy.estimator: new_tally.estimator = self_copy.estimator if self_copy.with_summary and other_copy.with_summary: @@ -1516,43 +1578,28 @@ class Tally(object): if self_copy.num_realizations == other_copy.num_realizations: new_tally.num_realizations = self_copy.num_realizations - # If filters are identical, simply reuse them in derived tally - if self_copy.filters == other_copy.filters: + # Add filters to the new tally + if filter_product == 'entrywise': for self_filter in self_copy.filters: - new_tally.add_filter(self_filter) - - # Generate filter entrywise product for non-identical filters + new_tally.filters.append(self_filter) else: + all_filters = [self_copy.filters, other_copy.filters] + for self_filter, other_filter in itertools.product(*all_filters): + new_filter = CrossFilter(self_filter, other_filter, binary_op) + new_tally.add_filter(new_filter) - # Find the common longest sequence of shared filters - match = 0 - for self_filter, other_filter in zip(self_copy.filters, other_copy.filters): - if self_filter == other_filter: - match += 1 - else: - break + # Add nuclides to the new tally + if nuclide_product == 'entrywise': + for self_nuclide in self_copy.nuclides: + new_tally.nuclides.append(self_nuclide) + else: + all_nuclides = [self_copy.nuclides, other_copy.nuclides] + for self_nuclide, other_nuclide in itertools.product(*all_nuclides): + new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op) + new_tally.add_nuclide(new_nuclide) - match_filters = self_copy.filters[:match] - cross_filters = [self_copy.filters[match:], other_copy.filters[match:]] - - # Simply reuse shared filters in derived tally - for filter in match_filters: - new_tally.add_filter(filter) - - # Use cross filters to combine non-shared filters in derived tally - if len(self_copy.filters) != match and len(other_copy.filters) == match: - for filter in cross_filters[0]: - new_tally.add_filter(filter) - elif len(self_copy.filters) == match and len(other_copy.filters) != match: - for filter in cross_filters[1]: - new_tally.add_filter(filter) - else: - for self_filter, other_filter in itertools.product(*cross_filters): - new_filter = CrossFilter(self_filter, other_filter, binary_op) - new_tally.add_filter(new_filter) - - # Generate score "outer products" - if self_copy.scores == other_copy.scores: + # Add scores to the new tally + if score_product == 'entrywise': new_tally.num_score_bins = self_copy.num_score_bins for self_score in self_copy.scores: new_tally.add_score(self_score) @@ -1563,16 +1610,6 @@ class Tally(object): new_score = CrossScore(self_score, other_score, binary_op) new_tally.add_score(new_score) - # Generate nuclide "outer products" - if self_copy.nuclides == other_copy.nuclides: - for self_nuclide in self_copy.nuclides: - new_tally.nuclides.append(self_nuclide) - else: - all_nuclides = [self_copy.nuclides, other_copy.nuclides] - for self_nuclide, other_nuclide in itertools.product(*all_nuclides): - new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op) - new_tally.add_nuclide(new_nuclide) - # Correct each Filter's stride stride = new_tally.num_nuclides * new_tally.num_score_bins for filter in reversed(new_tally.filters): @@ -1581,7 +1618,8 @@ class Tally(object): return new_tally - def _align_tally_data(self, other): + def _align_tally_data(self, other, filter_product, nuclide_product, + score_product): """Aligns data from two tallies for tally arithmetic. This is a helper method to construct a dict of dicts of the "aligned" @@ -1597,6 +1635,15 @@ class Tally(object): ---------- other : Tally The tally to outer product with this tally + filter_product : str + The type of product (tensor or entry) to be performed between filter + data. + nuclide_product : str + The type of product (tensor or entry) to be performed between nuclide + data. + score_product : str + The type of product (tensor or entry) to be performed between score + data. Returns ------- @@ -1610,38 +1657,22 @@ class Tally(object): other_missing_filters = set(self.filters).difference(set(other.filters)) self_missing_filters = set(other.filters).difference(set(self.filters)) - # Add other_missing_filters to other + # Add filters present in self but not in other to other for filter in other_missing_filters: filter = copy.deepcopy(filter) repeat_factor = filter.num_bins other._mean = np.repeat(other.mean, repeat_factor, axis=0) - other.sum = np.repeat(other.sum, repeat_factor, axis=0) other._std_dev = np.repeat(other.std_dev, repeat_factor, axis=0) - other.sum_sq = np.repeat(other.sum_sq, repeat_factor, axis=0) other.add_filter(filter) - # Correct the stride for other filters - stride = other.num_nuclides * other.num_score_bins - for filter in reversed(other.filters): - filter.stride = stride - stride *= filter.num_bins - - # Add self_missing_filters to self + # Add filters present in other but not in self to self for filter in self_missing_filters: filter = copy.deepcopy(filter) repeat_factor = filter.num_bins self._mean = np.repeat(self.mean, repeat_factor, axis=0) - self.sum = np.repeat(self.sum, repeat_factor, axis=0) self._std_dev = np.repeat(self.std_dev, repeat_factor, axis=0) - self.sum_sq = np.repeat(self.sum_sq, repeat_factor, axis=0) self.add_filter(filter) - # Correct the stride for self filters - stride = self.num_nuclides * self.num_score_bins - for filter in reversed(self.filters): - filter.stride = stride - stride *= filter.num_bins - # Align other filters with self filters for i, filter in enumerate(self.filters): other_index = other.filters.index(filter) @@ -1650,28 +1681,98 @@ class Tally(object): if other_index != i: other.swap_filters(filter, other.filters[i], inplace=True) + # Repeat and tile the data by nuclide in preparation for performing + # the tensor product across nuclides. + if nuclide_product == 'tensor': + self._mean = np.repeat(self.mean, other.num_nuclides, axis=1) + self._std_dev = np.repeat(self.std_dev, other.num_nuclides, axis=1) + other._mean = np.tile(other.mean, (1, self.num_nuclides, 1)) + other._std_dev = np.tile(other.std_dev, (1, self.num_nuclides, 1)) + + # Add nuclides to each tally such that each tally contains the complete + # set of nuclides necessary to perform an entrywise product. New nuclides + # added to a tally will have all their scores set to zero. + else: + + # Get the set of nuclides that each tally is missing + other_missing_nuclides = set(self.nuclides).difference(set(other.nuclides)) + self_missing_nuclides = set(other.nuclides).difference(set(self.nuclides)) + + # Add nuclides present in self but not in other to other + for nuclide in other_missing_nuclides: + other._mean = np.insert(other.mean, other.num_nuclides, 0, axis=1) + other._std_dev = np.insert(other.std_dev, other.num_nuclides, 0, axis=1) + other.add_nuclide(nuclide) + + # Add nuclides present in other but not in self to self + for nuclide in self_missing_nuclides: + self._mean = np.insert(self.mean, self.num_nuclides, 0, axis=1) + self._std_dev = np.insert(self.std_dev, self.num_nuclides, 0, axis=1) + self.add_nuclide(nuclide) + + # Align other nuclides with self nuclides + for i, nuclide in enumerate(self.nuclides): + other_index = other.nuclides.index(nuclide) + + # If necessary, swap other nuclide + if other_index != i: + other.swap_nuclides(nuclide, other.nuclides[i]) + + # Repeat and tile the data by score in preparation for performing + # the tensor product across scores. + if score_product == 'tensor': + self._mean = np.repeat(self.mean, other.num_score_bins, axis=2) + self._std_dev = np.repeat(self.std_dev, other.num_score_bins, axis=2) + other._mean = np.tile(other.mean, (1, 1, self.num_score_bins)) + other._std_dev = np.tile(other.std_dev, (1, 1, self.num_score_bins)) + + # Add scores to each tally such that each tally contains the complete set + # of scores necessary to perform an entrywise product. New scores added + # to a tally will be set to zero. + else: + + # Get the set of scores that each tally is missing + other_missing_scores = set(self.scores).difference(set(other.scores)) + self_missing_scores = set(other.scores).difference(set(self.scores)) + + # Add scores present in self but not in other to other + for score in other_missing_scores: + other._mean = np.insert(other.mean, other.num_score_bins, 0, axis=2) + other._std_dev = np.insert(other.std_dev, other.num_score_bins, 0, axis=2) + other.add_score(score) + + # Add scores present in other but not in self to self + for score in self_missing_scores: + self._mean = np.insert(self.mean, self.num_score_bins, 0, axis=2) + self._std_dev = np.insert(self.std_dev, self.num_score_bins, 0, axis=2) + self.add_score(score) + + # Align other scores with self scores + for i, score in enumerate(self.scores): + other_index = other.scores.index(score) + + # If necessary, swap other score + if other_index != i: + other.swap_scores(score, other.scores[i]) + + # Correct the stride for other filters + stride = other.num_nuclides * other.num_score_bins + for filter in reversed(other.filters): + filter.stride = stride + stride *= filter.num_bins + + # Correct the stride for self filters + stride = self.num_nuclides * self.num_score_bins + for filter in reversed(self.filters): + filter.stride = stride + stride *= filter.num_bins + # Deep copy the mean and std dev data self_mean = copy.deepcopy(self.mean) self_std_dev = copy.deepcopy(self.std_dev) other_mean = copy.deepcopy(other.mean) other_std_dev = copy.deepcopy(other.std_dev) - # If the tallies do not have identical sets of nuclides, tile and repeat - # to perform cross product of data for each nuclide - if self.nuclides != other.nuclides: - self_mean = np.repeat(self_mean, other.num_nuclides, axis=1) - self_std_dev = np.repeat(self_std_dev, other.num_nuclides, axis=1) - other_mean = np.tile(other_mean, (1, self.num_nuclides, 1)) - other_std_dev = np.tile(other_std_dev, (1, self.num_nuclides, 1)) - - # If the tallies do not have identical sets of scores, tile and repeat - # to perform cross product of data for each score - if self.scores != other.scores: - self_mean = np.repeat(self_mean, other.num_score_bins, axis=2) - self_std_dev = np.repeat(self_std_dev, other.num_score_bins, axis=2) - other_mean = np.tile(other_mean, (1, 1, self.num_score_bins)) - other_std_dev = np.tile(other_std_dev, (1, 1, self.num_score_bins)) - data = {} data['self'] = {} data['other'] = {} @@ -1808,6 +1909,116 @@ class Tally(object): if not inplace: return swap_tally + def swap_nuclides(self, nuclide1, nuclide2): + """Reverse the ordering of two nuclides in this tally + + This is a helper method for tally arithmetic which helps align the data + in two tallies with shared nuclides. This method reverses the order of + the two nuclides in place. + + Parameters + ---------- + nuclide1 : Nuclide + The nuclide to swap with nuclide2 + + nuclide2 : Nuclide + The nuclide to swap with nuclide1 + + Raises + ------ + ValueError + If this is a derived tally or this method is called before the tally + is populated with data by the StatePoint.read_results() method. + + """ + + # Check that results have been read + if not self.derived and self.sum is None: + msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ + 'since it does not contain any results.'.format(self.id) + raise ValueError(msg) + + # Swap the nuclides in the Tally + nuclide1_index = self.nuclides.index(nuclide1) + nuclide2_index = self.nuclides.index(nuclide2) + self.nuclides[nuclide1_index] = nuclide2 + self.nuclides[nuclide2_index] = nuclide1 + + # Copy the tally data + self_mean = copy.deepcopy(self.mean) + self_std_dev = copy.deepcopy(self.std_dev) + + # Swap nuclide 1 in place of nuclide 2 + self._mean = np.delete(self.mean, nuclide2_index, axis=1) + self._mean = np.insert(self.mean, nuclide2_index, + self_mean[:,nuclide1_index,:], axis=1) + self._std_dev = np.delete(self.std_dev, nuclide2_index, axis=1) + self._std_dev = np.insert(self.std_dev, nuclide2_index, + self_mean[:,nuclide1_index,:], axis=1) + + # Swap nuclide 2 in place of nuclide 1 + self._mean = np.delete(self.mean, nuclide1_index, axis=1) + self._mean = np.insert(self.mean, nuclide1_index, + self_mean[:,nuclide2_index,:], axis=1) + self._std_dev = np.delete(self.std_dev, nuclide1_index, axis=1) + self._std_dev = np.insert(self.std_dev, nuclide1_index, + self_mean[:,nuclide2_index,:], axis=1) + + def swap_scores(self, score1, score2): + """Reverse the ordering of two scores in this tally + + This is a helper method for tally arithmetic which helps align the data + in two tallies with shared scores. This method copies reverses the order + of the two scores in place. + + Parameters + ---------- + score1 : Score + The score to swap with score2 + + score2 : Score + The score to swap with score1 + + Raises + ------ + ValueError + If this is a derived tally or this method is called before the tally + is populated with data by the StatePoint.read_results() method. + + """ + + # Check that results have been read + if not self.derived and self.sum is None: + msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ + 'since it does not contain any results.'.format(self.id) + raise ValueError(msg) + + # Swap the scores in the Tally + score1_index = self.scores.index(score1) + score2_index = self.scores.index(score2) + self.scores[score1_index] = score2 + self.scores[score2_index] = score1 + + # Copy the tally data + self_mean = copy.deepcopy(self.mean) + self_std_dev = copy.deepcopy(self.std_dev) + + # Swap score 1 in place of score 2 + self._mean = np.delete(self.mean, score2_index, axis=2) + self._mean = np.insert(self.mean, score2_index, + self_mean[:,:,score1_index], axis=2) + self._std_dev = np.delete(self.std_dev, score2_index, axis=2) + self._std_dev = np.insert(self.std_dev, score2_index, + self_mean[:,:,score1_index], axis=2) + + # Swap score 2 in place of score 1 + self._mean = np.delete(self.mean, score1_index, axis=2) + self._mean = np.insert(self.mean, score1_index, + self_mean[:,:,score2_index], axis=2) + self._std_dev = np.delete(self.std_dev, score1_index, axis=2) + self._std_dev = np.insert(self.std_dev, score1_index, + self_mean[:,:,score2_index], axis=2) + def __add__(self, other): """Adds this tally to another tally or scalar value. @@ -1851,7 +2062,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._hybrid_product(other, binary_op='+') + new_tally = self.hybrid_product(other, binary_op='+') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -1921,7 +2132,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._hybrid_product(other, binary_op='-') + new_tally = self.hybrid_product(other, binary_op='-') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -1991,7 +2202,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._hybrid_product(other, binary_op='*') + new_tally = self.hybrid_product(other, binary_op='*') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -2061,7 +2272,7 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._hybrid_product(other, binary_op='/') + new_tally = self.hybrid_product(other, binary_op='/') elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -2134,7 +2345,7 @@ class Tally(object): raise ValueError(msg) if isinstance(power, Tally): - new_tally = self._hybrid_product(power, binary_op='^') + new_tally = self.hybrid_product(power, binary_op='^') elif isinstance(power, Real): new_tally = Tally(name='derived') From cae38c30bded40da7b963b0f4068ecb59eb4ce14 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 10 Dec 2015 01:03:17 -0500 Subject: [PATCH 11/11] fixed typos and removed unnecessary repeat_factor variable --- openmc/tallies.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 5656987cbf..8b4aa89ddd 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1519,7 +1519,7 @@ class Tally(object): new_tally.name = new_name # Query the mean and std dev so the tally data is read in from file - # if it has not present + # if it has not already been read in. self.mean self.std_dev other.mean @@ -1636,14 +1636,14 @@ class Tally(object): other : Tally The tally to outer product with this tally filter_product : str - The type of product (tensor or entry) to be performed between filter - data. + The type of product (tensor or entrywise) to be performed between + filter data. nuclide_product : str - The type of product (tensor or entry) to be performed between nuclide - data. + The type of product (tensor or entrywise) to be performed between + nuclide data. score_product : str - The type of product (tensor or entry) to be performed between score - data. + The type of product (tensor or entrywise) to be performed between + score data. Returns ------- @@ -1660,17 +1660,15 @@ class Tally(object): # Add filters present in self but not in other to other for filter in other_missing_filters: filter = copy.deepcopy(filter) - repeat_factor = filter.num_bins - other._mean = np.repeat(other.mean, repeat_factor, axis=0) - other._std_dev = np.repeat(other.std_dev, repeat_factor, axis=0) + other._mean = np.repeat(other.mean, filter.num_bins, axis=0) + other._std_dev = np.repeat(other.std_dev, filter.num_bins, axis=0) other.add_filter(filter) # Add filters present in other but not in self to self for filter in self_missing_filters: filter = copy.deepcopy(filter) - repeat_factor = filter.num_bins - self._mean = np.repeat(self.mean, repeat_factor, axis=0) - self._std_dev = np.repeat(self.std_dev, repeat_factor, axis=0) + self._mean = np.repeat(self.mean, filter.num_bins, axis=0) + self._std_dev = np.repeat(self.std_dev, filter.num_bins, axis=0) self.add_filter(filter) # Align other filters with self filters