diff --git a/docs/source/conf.py b/docs/source/conf.py index 672217c84..4e21ec20b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -25,7 +25,7 @@ except ImportError: MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial', - 'h5py', 'pandas', 'opencg'] + 'h5py', 'pandas', 'uncertainties', 'opencg'] sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES) import numpy as np diff --git a/docs/source/developers.rst b/docs/source/developers.rst index d82d89744..0bdae2416 100644 --- a/docs/source/developers.rst +++ b/docs/source/developers.rst @@ -13,6 +13,7 @@ Active development of the OpenMC Monte Carlo code is currently led by: * `Jon Walsh `_ * `Sterling Harper `_ * `Will Boyd `_ +* `Samuel Shaner `_ * `Benoit Forget `_ * `Kord Smith `_ -* `Andrew Siegel `_ +* `Andrew Siegel `_ diff --git a/docs/source/devguide/workflow.rst b/docs/source/devguide/workflow.rst index 2f36436ca..a53bd114b 100644 --- a/docs/source/devguide/workflow.rst +++ b/docs/source/devguide/workflow.rst @@ -238,4 +238,4 @@ from your private repository into a public fork. .. _paid plan: https://github.com/plans .. _Bitbucket: https://bitbucket.org .. _ctest: http://www.cmake.org/cmake/help/v2.8.12/ctest.html -.. _NNDC: http://http://www.nndc.bnl.gov/endf/b7.1/acefiles.html +.. _NNDC: http://www.nndc.bnl.gov/endf/b7.1/acefiles.html diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index c2e7270b6..36c252bb1 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -384,7 +384,7 @@ x - x_0`, :math:`\bar{y} = y - y_0`, and :math:`\bar{z} = z - z_0`. We then have Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain .. math:: - :label: dist-xcylinder-2 + :label: dist-xcone-2 (v^2 + w^2 - R^2u^2) d^2 + 2 (\bar{y}v + \bar{z}w - R^2\bar{x}u) d + (\bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2) = 0 @@ -392,7 +392,7 @@ Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain Defining the terms .. math:: - :label: dist-quadric-terms + :label: dist-xcone-terms a = v^2 + w^2 - R^2u^2 @@ -896,4 +896,4 @@ Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, the gradient to the surface is .. _surfaces: http://en.wikipedia.org/wiki/Surface .. _MCNP: http://mcnp.lanl.gov .. _Serpent: http://montecarlo.vtt.fi -.. _Monte Carlo Performance benchmark: https://github.com/paulromano/benchmarks/tree/master/mc-performance/openmc +.. _Monte Carlo Performance benchmark: https://github.com/mit-crpg/benchmarks/tree/master/mc-performance/openmc diff --git a/docs/source/publications.rst b/docs/source/publications.rst index c5b29b194..097cd43d0 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -82,6 +82,10 @@ Coupling and Multi-physics Geometry -------- +- Derek M. Lax, "Memory efficient indexing algorithm for physical properties in + OpenMC," S. M. Thesis, Massachusetts Institute of Technology + (2015). ``_ + - Derek Lax, William Boyd, Nicholas Horelik, Benoit Forget, and Kord Smith, "A memory efficient algorithm for classifying unique regions in constructive solid geometries," *Proc. PHYSOR*, Kyoto, Japan, Sep. 28--Oct. 3 (2014). @@ -106,6 +110,10 @@ Miscellaneous triggers for the OpenMC Monte Carlo code," *Trans. Am. Nucl. Soc.*, **112**, 637-640 (2015). +- Kyungkwan Noh and Deokjung Lee, "Whole Core Analysis using OpenMC Monte Carlo + Code," *Trans. Kor. Nucl. Soc. Autumn Meeting*, Gyeongju, Korea, + Oct. 24-25, 2013. + - 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). @@ -139,7 +147,7 @@ Doppler Broadening - Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed multipole for cross section Doppler broadening," *J. Comput. Phys.*, **307**, 715-727 - (2016). ``_ + (2016). ``_ - Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown, "On-the-fly Doppler Broadening of Unresolved Resonance Region Cross Sections @@ -169,6 +177,11 @@ Nuclear Data - Paul K. Romano and Sterling M. Harper, "Nuclear data processing capabilities in OpenMC", *Proc. Nuclear Data*, Sep. 11-16, 2016. +- Jonathan A. Walsh, Benoit Froget, Kord S. Smith, and Forrest B. Brown, + "Neutron Cross Section Processing Methods for Improved Integral Benchmarking + of Unresolved Resonance Region Evaluations," *Eur. Phys. J. Web Conf.* + **111**, 06001 (2016). ``_ + - 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 @@ -213,10 +226,19 @@ Parallelism memory subsystem on Monte Carlo code performance," *Proc. Joint Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015). +- Hajime Fujita, Nan Dun, Aiman Fang, Zachary A. Rubinstein, Ziming Zheng, Kamil + Iskra, Jeff Hammonds, Anshu Dubey, Pavan Balaji, and Andrew A. Chien, "Using + Global View Resilience (GVR) to add Resilience to Exascale Applications," + *Proc. Supercomputing*, New Orleans, Louisiana, Nov. 16--21, 2014. + - Nicholas Horelik, Benoit Forget, Kord Smith, and Andrew Siegel, "Domain decomposition and terabyte tallies with the OpenMC Monte Carlo neutron transport code," *Proc. PHYSOR*, Kyoto Japan, Sep. 28--Oct. 3 (2014). +- John R. Tramm, Andrew R. Siegel, Tanzima Islam, and Martin Schulz, "XSBench -- + the development and verification of a performance abstraction for Monte Carlo + reactor analysis," *Proc. PHYSOR*, Kyoto, Japan, Sep 28--Oct. 3, 2014. + - Nicholas Horelik, Andrew Siegel, Benoit Forget, and Kord Smith, "Monte Carlo domain decomposition for robust nuclear reactor analysis," *Parallel Comput.*, **40**, 646--660 (2014). ``_ @@ -270,6 +292,10 @@ Parallelism Depletion --------- +- Anas Gul, K. S. Chaudri, R. Khan, and M. Azeen, "Development and verification + of LOOP: A Linkage of ORIGEN2.2 and OpenMC," *Ann. Nucl. Energy*, **99**, + 321--327 (2017). ``_ + - Kai Huang, Hongchun Wu, Yunzhao Li, and Liangzhi Cao, "Generalized depletion chain simplification based of significance analysis," *Proc. PHYSOR*, Sun Valley, Idaho, May 1-5, 2016. diff --git a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb index cd011d862..312bb8ee6 100644 --- a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb +++ b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb @@ -429,7 +429,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX////pgJFyEhJNv8RV\nUZDeAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AKHxEkM8uGp70AAAWFSURBVGje7Zs7cttADIZ9CSvX\ncrP0iCxUqbBc8Ag6xR6BhV2EvYvwFD4CCx1ABT1jMdgndpegRQnOrCbjpPlGESISC4A/gd27e8H5\n83CX3b4+iKJrRHkS4vkghMPBonRYWGwtfgD2YN+dRDUOoh6lACw0Noi9w2fESuEoAR/uVuMolX03\n9oXGT7F3eFL2iEfhUX1f4cPdL/ishs+68ai+udE4xPhexbjX2FfjGNoPj/DPNX4Tsd+EODr8FvsV\ndf1Hd9P2VvCi4+s/aXvrf+upAD+1/9GV1mkOH5X9vV6THtfvACslcaUCbESL61drBPtdI8SrFMWr\nELsXCkuFDYW75gbiP7d9Cf7bAYI/aCwUShrBvh30+lWQkzVgZ/HD4OixNCgcQpJ3BxU/Ln91elKo\nM5VEE38QtJ+Yv6cQ9xjKNYayyl8TypP8DfJnQ2H/b/N3ye9P83cT33SQv/sQh9gV7zZ/0dNj5HQa\nC5vVzv9+/WFN2w8KVaZ2BwL1+pv4g0x1QRfjq0dB4Q3kT277oP6VNL6gKxNU9a8zK+WLbi/Wwpdi\nhbboKqyxFOulHMj6v4W/AXbmUeAxrv9J/CqEBXaRKsXaodD4nsYvkT/G6H1D4SR/iPy1Roj9JsQ5\ne18/7EUHv1+Fvx/Xj5V9Ugb5K8TW4TZEEdcvoz/up0VTe9qsVIppKVX6a7D6y9ZvwEKjrtQxPtv6\nfXII9vCxKOGaIeAIfEF8IvAG8ie3vRK9rRQl+PPpSctbhfpTUCpviH+kxsZgpT91+snoX1l49KK3\niUQvICRy5aUw6l8leoVwoo3Uv1rKreF/UFLY6d9QP4L9Wf2r7EP9GOSfcsjZ56f60kz+XmVPXv+R\nuP49ff0T/53Rv6n/7m2lvXT9Wqd/VUz8hvh5M/ED6ILmt4mfHYZSaePnTWpsf/SvqV9O6dLYYClL\nEetnoH/LBLFoBvrX189uTv8++kot5vTvQD4/9jP690g9P/4z/bvo/XVG/xYoZZx+8fr3MxAtsf7t\nUOkG2JqsTtCIpgCt/qX1226KqZS7gfzJbe+c9jLrtIZ8lXD+s4umlW6AKIVrlML2/cXjgPFjlJqI\nRC+Fj0bVJe+vSh56pSdR6YkQ1ygF10Wqf0FeLta/iKn9Mv1L24ti2e+7W4n1b3T/W+L+t9H9T/Sv\nVboUmqJJon1/hZq8LnzRDlDrX1u0xRT1+6vEpomMmyYkqi95vIH8yW1PN+122KkLcNLKi/WTF01z\n/cNASrWE/l3ev6T17zX909z9X27/euK/Rf3zWP+Waf9eEv37KkWJ+rfDl6ZglNDa+cEBhwYDvkoN\nP/rX69814NaI3imq0l7OYDy/qSdDGwr7r+Y3VbzoKZr6XX2lfxfOb87qXzr+b1j/Xlp/nP6dn98M\ncdH7cn7zjPObKsYWS3Eb9w8n85smHtqQuPuZ30T2dlIT6F9xFl+n8xslegL9a4c2KRr9W4rp/GYq\numiM9Nec/j2v/yj9u1h//hv9e93vc++f63/u+rPjL3f+5Lbn1j9m/eXWf+7zh/v8+2b9e/Hzn6s/\nuPqHrb8g71n6L3f+5Lbnvn8w33+4718/+5d47//c/gO7/5E7/nPbc/tv3P4fs//I7X9y+6/fqH+v\n6j9z+9/c/ju3/8+eP+TOn9z23PkXc/7Gnf9x5483q38Xzn+582fu/Js9fy8kb/6fO39y23P3n3S8\n/S/c/Tfc/T83uX/pgv1XE/9duP+Lu/+Mvf8td/znti8kb/8ld/9nx9t/Sjw/Ltr/yt1/+337f6/b\nf0zoB3nJ/ucVc/81d/83e/957vzJbc89/8A8f8E9/5HE78XnT/4H/cs5f8Q9/8Q9f8U+/5U7f3Lb\nc88fdrzzjyvm+cuf/Uu887/c88fs88954/8vO4SjPC+2QRIAAAAldEVYdGRhdGU6Y3JlYXRlADIw\nMTYtMTAtMzFUMTI6MzY6NTEtMDU6MDD2oqCYAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTEwLTMx\nVDEyOjM2OjUxLTA1OjAwh/8YJAAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX////pgJFyEhJNv8RV\nUZDeAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+ALDQ8YOYloD4IAAAWFSURBVGje7Zs7cttADIZ9CSvX\ncrP0iCxUqbBc8Ag6xR6BhV2EvYvwFD4CCx1ABT1jMdgndpegRQnOrCbjpPlGESISC4A/gd27e8H5\n83CX3b4+iKJrRHkS4vkghMPBonRYWGwtfgD2YN+dRDUOoh6lACw0Noi9w2fESuEoAR/uVuMolX03\n9oXGT7F3eFL2iEfhUX1f4cPdL/ishs+68ai+udE4xPhexbjX2FfjGNoPj/DPNX4Tsd+EODr8FvsV\ndf1Hd9P2VvCi4+s/aXvrf+upAD+1/9GV1mkOH5X9vV6THtfvACslcaUCbESL61drBPtdI8SrFMWr\nELsXCkuFDYW75gbiP7d9Cf7bAYI/aCwUShrBvh30+lWQkzVgZ/HD4OixNCgcQpJ3BxU/Ln91elKo\nM5VEE38QtJ+Yv6cQ9xjKNYayyl8TypP8DfJnQ2H/b/N3ye9P83cT33SQv/sQh9gV7zZ/0dNj5HQa\nC5vVzv9+/WFN2w8KVaZ2BwL1+pv4g0x1QRfjq0dB4Q3kT277oP6VNL6gKxNU9a8zK+WLbi/Wwpdi\nhbboKqyxFOulHMj6v4W/AXbmUeAxrv9J/CqEBXaRKsXaodD4nsYvkT/G6H1D4SR/iPy1Roj9JsQ5\ne18/7EUHv1+Fvx/Xj5V9Ugb5K8TW4TZEEdcvoz/up0VTe9qsVIppKVX6a7D6y9ZvwEKjrtQxPtv6\nfXII9vCxKOGaIeAIfEF8IvAG8ie3vRK9rRQl+PPpSctbhfpTUCpviH+kxsZgpT91+snoX1l49KK3\niUQvICRy5aUw6l8leoVwoo3Uv1rKreF/UFLY6d9QP4L9Wf2r7EP9GOSfcsjZ56f60kz+XmVPXv+R\nuP49ff0T/53Rv6n/7m2lvXT9Wqd/VUz8hvh5M/ED6ILmt4mfHYZSaePnTWpsf/SvqV9O6dLYYClL\nEetnoH/LBLFoBvrX189uTv8++kot5vTvQD4/9jP690g9P/4z/bvo/XVG/xYoZZx+8fr3MxAtsf7t\nUOkG2JqsTtCIpgCt/qX1226KqZS7gfzJbe+c9jLrtIZ8lXD+s4umlW6AKIVrlML2/cXjgPFjlJqI\nRC+Fj0bVJe+vSh56pSdR6YkQ1ygF10Wqf0FeLta/iKn9Mv1L24ti2e+7W4n1b3T/W+L+t9H9T/Sv\nVboUmqJJon1/hZq8LnzRDlDrX1u0xRT1+6vEpomMmyYkqi95vIH8yW1PN+122KkLcNLKi/WTF01z\n/cNASrWE/l3ev6T17zX909z9X27/euK/Rf3zWP+Waf9eEv37KkWJ+rfDl6ZglNDa+cEBhwYDvkoN\nP/rX69814NaI3imq0l7OYDy/qSdDGwr7r+Y3VbzoKZr6XX2lfxfOb87qXzr+b1j/Xlp/nP6dn98M\ncdH7cn7zjPObKsYWS3Eb9w8n85smHtqQuPuZ30T2dlIT6F9xFl+n8xslegL9a4c2KRr9W4rp/GYq\numiM9Nec/j2v/yj9u1h//hv9e93vc++f63/u+rPjL3f+5Lbn1j9m/eXWf+7zh/v8+2b9e/Hzn6s/\nuPqHrb8g71n6L3f+5Lbnvn8w33+4718/+5d47//c/gO7/5E7/nPbc/tv3P4fs//I7X9y+6/fqH+v\n6j9z+9/c/ju3/8+eP+TOn9z23PkXc/7Gnf9x5483q38Xzn+582fu/Js9fy8kb/6fO39y23P3n3S8\n/S/c/Tfc/T83uX/pgv1XE/9duP+Lu/+Mvf8td/znti8kb/8ld/9nx9t/Sjw/Ltr/yt1/+337f6/b\nf0zoB3nJ/ucVc/81d/83e/957vzJbc89/8A8f8E9/5HE78XnT/4H/cs5f8Q9/8Q9f8U+/5U7f3Lb\nc88fdrzzjyvm+cuf/Uu887/c88fs88954/8vO4SjPC+2QRIAAAAldEVYdGRhdGU6Y3JlYXRlADIw\nMTYtMTEtMTNUMTU6MjQ6NTctMDU6MDCpBKn7AAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTExLTEz\nVDE1OjI0OjU3LTA1OjAw2FkRRwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] @@ -575,7 +575,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/home/romano/openmc/openmc/mgxs/library.py:370: RuntimeWarning: The P0 correction will be ignored since the scattering order 0 is greater than zero\n", + "/home/nelsonag/git/openmc/openmc/mgxs/library.py:398: RuntimeWarning: The P0 correction will be ignored since the scattering order 0 is greater than zero\n", " warn(msg, RuntimeWarning)\n" ] } @@ -731,9 +731,9 @@ " Copyright | 2011-2016 Massachusetts Institute of Technology\n", " License | http://openmc.readthedocs.io/en/latest/license.html\n", " Version | 0.8.0\n", - " Git SHA1 | da5563eddb5f2c2d6b2c9839d518de40962b78f2\n", - " Date/Time | 2016-10-31 12:36:52\n", - " OpenMP Threads | 4\n", + " Git SHA1 | 1a921e7d08fc41b72bf1dd65cd17e922222b78b1\n", + " Date/Time | 2016-11-13 15:24:57\n", + " OpenMP Threads | 8\n", "\n", " ===========================================================================\n", " ========================> INITIALIZATION <=========================\n", @@ -743,12 +743,12 @@ " Reading geometry XML file...\n", " Reading materials XML file...\n", " Reading cross sections XML file...\n", - " Reading U235 from /home/romano/openmc/scripts/nndc_hdf5/U235.h5\n", - " Reading U238 from /home/romano/openmc/scripts/nndc_hdf5/U238.h5\n", - " Reading O16 from /home/romano/openmc/scripts/nndc_hdf5/O16.h5\n", - " Reading Zr90 from /home/romano/openmc/scripts/nndc_hdf5/Zr90.h5\n", - " Reading H1 from /home/romano/openmc/scripts/nndc_hdf5/H1.h5\n", - " Reading B10 from /home/romano/openmc/scripts/nndc_hdf5/B10.h5\n", + " Reading U235 from /opt/xsdata/nndc/U235.h5\n", + " Reading U238 from /opt/xsdata/nndc/U238.h5\n", + " Reading O16 from /opt/xsdata/nndc/O16.h5\n", + " Reading Zr90 from /opt/xsdata/nndc/Zr90.h5\n", + " Reading H1 from /opt/xsdata/nndc/H1.h5\n", + " Reading B10 from /opt/xsdata/nndc/B10.h5\n", " Maximum neutron transport energy: 2.00000E+07 eV for U235\n", " Reading tallies XML file...\n", " Building neighboring cells lists for each surface...\n", @@ -819,20 +819,20 @@ "\n", " =======================> TIMING STATISTICS <=======================\n", "\n", - " Total time for initialization = 5.3479E-01 seconds\n", - " Reading cross sections = 3.8403E-01 seconds\n", - " Total time in simulation = 3.9455E+01 seconds\n", - " Time in transport only = 3.9331E+01 seconds\n", - " Time in inactive batches = 4.6776E+00 seconds\n", - " Time in active batches = 3.4778E+01 seconds\n", - " Time synchronizing fission bank = 1.0514E-02 seconds\n", - " Sampling source sites = 7.2826E-03 seconds\n", - " SEND/RECV source sites = 3.1267E-03 seconds\n", - " Time accumulating tallies = 3.7644E-04 seconds\n", - " Total time for finalization = 7.7600E-06 seconds\n", - " Total time elapsed = 4.0029E+01 seconds\n", - " Calculation Rate (inactive) = 10689.2 neutrons/second\n", - " Calculation Rate (active) = 5750.82 neutrons/second\n", + " Total time for initialization = 2.4361E-01 seconds\n", + " Reading cross sections = 1.6584E-01 seconds\n", + " Total time in simulation = 7.3766E+00 seconds\n", + " Time in transport only = 7.3482E+00 seconds\n", + " Time in inactive batches = 8.5442E-01 seconds\n", + " Time in active batches = 6.5222E+00 seconds\n", + " Time synchronizing fission bank = 4.9435E-03 seconds\n", + " Sampling source sites = 3.3815E-03 seconds\n", + " SEND/RECV source sites = 1.5249E-03 seconds\n", + " Time accumulating tallies = 9.5694E-05 seconds\n", + " Total time for finalization = 2.9260E-06 seconds\n", + " Total time elapsed = 7.6377E+00 seconds\n", + " Calculation Rate (inactive) = 58519.0 neutrons/second\n", + " Calculation Rate (active) = 30664.7 neutrons/second\n", "\n", " ============================> RESULTS <============================\n", "\n", @@ -970,7 +970,20 @@ "metadata": { "collapsed": false }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/nelsonag/git/openmc/openmc/tallies.py:1944: RuntimeWarning: invalid value encountered in true_divide\n", + " self_rel_err = data['self']['std. dev.'] / data['self']['mean']\n", + "/home/nelsonag/git/openmc/openmc/tallies.py:1945: RuntimeWarning: invalid value encountered in true_divide\n", + " other_rel_err = data['other']['std. dev.'] / data['other']['mean']\n", + "/home/nelsonag/git/openmc/openmc/tallies.py:1946: RuntimeWarning: invalid value encountered in true_divide\n", + " new_tally._mean = data['self']['mean'] / data['other']['mean']\n" + ] + } + ], "source": [ "# Create a MGXS File which can then be written to disk\n", "mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=['fuel', 'zircaloy', 'water'])\n", @@ -1107,9 +1120,9 @@ " Copyright | 2011-2016 Massachusetts Institute of Technology\n", " License | http://openmc.readthedocs.io/en/latest/license.html\n", " Version | 0.8.0\n", - " Git SHA1 | da5563eddb5f2c2d6b2c9839d518de40962b78f2\n", - " Date/Time | 2016-10-31 12:37:32\n", - " OpenMP Threads | 4\n", + " Git SHA1 | 1a921e7d08fc41b72bf1dd65cd17e922222b78b1\n", + " Date/Time | 2016-11-13 15:25:05\n", + " OpenMP Threads | 8\n", "\n", " ===========================================================================\n", " ========================> INITIALIZATION <=========================\n", @@ -1133,56 +1146,56 @@ "\n", " Bat./Gen. k Average k \n", " ========= ======== ==================== \n", - " 1/1 1.00711 \n", - " 2/1 1.01538 \n", - " 3/1 1.01664 \n", - " 4/1 1.03592 \n", - " 5/1 1.00771 \n", - " 6/1 1.00555 \n", - " 7/1 1.02573 \n", - " 8/1 1.04322 \n", - " 9/1 1.02270 \n", - " 10/1 1.02354 \n", - " 11/1 1.02023 \n", - " 12/1 1.03047 1.02535 +/- 0.00512\n", - " 13/1 1.04476 1.03182 +/- 0.00711\n", - " 14/1 1.02223 1.02942 +/- 0.00557\n", - " 15/1 1.02082 1.02770 +/- 0.00465\n", - " 16/1 1.01472 1.02554 +/- 0.00437\n", - " 17/1 1.02104 1.02489 +/- 0.00375\n", - " 18/1 1.04471 1.02737 +/- 0.00408\n", - " 19/1 1.02806 1.02745 +/- 0.00360\n", - " 20/1 1.02044 1.02675 +/- 0.00330\n", - " 21/1 1.02592 1.02667 +/- 0.00298\n", - " 22/1 1.02242 1.02632 +/- 0.00275\n", - " 23/1 0.99969 1.02427 +/- 0.00325\n", - " 24/1 1.02213 1.02412 +/- 0.00301\n", - " 25/1 1.02080 1.02390 +/- 0.00281\n", - " 26/1 1.01033 1.02305 +/- 0.00277\n", - " 27/1 1.02881 1.02339 +/- 0.00262\n", - " 28/1 1.01649 1.02300 +/- 0.00250\n", - " 29/1 1.03817 1.02380 +/- 0.00250\n", - " 30/1 1.00958 1.02309 +/- 0.00247\n", - " 31/1 1.01811 1.02285 +/- 0.00236\n", - " 32/1 1.02709 1.02305 +/- 0.00226\n", - " 33/1 1.01823 1.02284 +/- 0.00217\n", - " 34/1 1.01208 1.02239 +/- 0.00213\n", - " 35/1 1.01380 1.02204 +/- 0.00207\n", - " 36/1 1.02358 1.02210 +/- 0.00199\n", - " 37/1 1.03653 1.02264 +/- 0.00199\n", - " 38/1 1.03117 1.02294 +/- 0.00194\n", - " 39/1 1.00915 1.02247 +/- 0.00193\n", - " 40/1 1.03107 1.02275 +/- 0.00189\n", - " 41/1 1.02316 1.02277 +/- 0.00182\n", - " 42/1 1.02677 1.02289 +/- 0.00177\n", - " 43/1 0.99361 1.02200 +/- 0.00193\n", - " 44/1 1.04841 1.02278 +/- 0.00203\n", - " 45/1 0.99768 1.02206 +/- 0.00210\n", - " 46/1 1.02694 1.02220 +/- 0.00204\n", - " 47/1 1.03540 1.02256 +/- 0.00202\n", - " 48/1 1.03539 1.02289 +/- 0.00199\n", - " 49/1 1.02498 1.02295 +/- 0.00194\n", - " 50/1 1.00692 1.02255 +/- 0.00193\n", + " 1/1 0.98369 \n", + " 2/1 1.01520 \n", + " 3/1 1.03642 \n", + " 4/1 1.02658 \n", + " 5/1 1.03102 \n", + " 6/1 1.05382 \n", + " 7/1 1.01978 \n", + " 8/1 1.01753 \n", + " 9/1 1.02420 \n", + " 10/1 0.99889 \n", + " 11/1 1.04874 \n", + " 12/1 1.01382 1.03128 +/- 0.01746\n", + " 13/1 1.03987 1.03414 +/- 0.01048\n", + " 14/1 1.02282 1.03131 +/- 0.00793\n", + " 15/1 1.03282 1.03162 +/- 0.00615\n", + " 16/1 0.99669 1.02579 +/- 0.00769\n", + " 17/1 1.00052 1.02218 +/- 0.00743\n", + " 18/1 1.01124 1.02082 +/- 0.00658\n", + " 19/1 1.00629 1.01920 +/- 0.00602\n", + " 20/1 1.05322 1.02260 +/- 0.00637\n", + " 21/1 1.00763 1.02124 +/- 0.00592\n", + " 22/1 1.01841 1.02101 +/- 0.00541\n", + " 23/1 1.03430 1.02203 +/- 0.00508\n", + " 24/1 1.03064 1.02264 +/- 0.00474\n", + " 25/1 1.03272 1.02331 +/- 0.00447\n", + " 26/1 1.01226 1.02262 +/- 0.00424\n", + " 27/1 1.00883 1.02181 +/- 0.00406\n", + " 28/1 1.02712 1.02211 +/- 0.00384\n", + " 29/1 1.03146 1.02260 +/- 0.00367\n", + " 30/1 1.02964 1.02295 +/- 0.00350\n", + " 31/1 0.99832 1.02178 +/- 0.00353\n", + " 32/1 1.03420 1.02234 +/- 0.00341\n", + " 33/1 1.01860 1.02218 +/- 0.00326\n", + " 34/1 1.03328 1.02264 +/- 0.00316\n", + " 35/1 1.01865 1.02248 +/- 0.00303\n", + " 36/1 1.02643 1.02264 +/- 0.00292\n", + " 37/1 1.01070 1.02219 +/- 0.00284\n", + " 38/1 1.01871 1.02207 +/- 0.00274\n", + " 39/1 0.98827 1.02090 +/- 0.00289\n", + " 40/1 1.01740 1.02079 +/- 0.00279\n", + " 41/1 1.02920 1.02106 +/- 0.00272\n", + " 42/1 1.02496 1.02118 +/- 0.00263\n", + " 43/1 1.04288 1.02184 +/- 0.00264\n", + " 44/1 1.03749 1.02230 +/- 0.00260\n", + " 45/1 1.04338 1.02290 +/- 0.00259\n", + " 46/1 1.03146 1.02314 +/- 0.00253\n", + " 47/1 1.04668 1.02377 +/- 0.00254\n", + " 48/1 1.02707 1.02386 +/- 0.00248\n", + " 49/1 1.02589 1.02391 +/- 0.00241\n", + " 50/1 1.02100 1.02384 +/- 0.00235\n", " Creating state point statepoint.50.h5...\n", "\n", " ===========================================================================\n", @@ -1192,27 +1205,27 @@ "\n", " =======================> TIMING STATISTICS <=======================\n", "\n", - " Total time for initialization = 3.6445E-02 seconds\n", - " Reading cross sections = 4.9377E-03 seconds\n", - " Total time in simulation = 3.0983E+01 seconds\n", - " Time in transport only = 3.0902E+01 seconds\n", - " Time in inactive batches = 3.2106E+00 seconds\n", - " Time in active batches = 2.7772E+01 seconds\n", - " Time synchronizing fission bank = 9.7451E-03 seconds\n", - " Sampling source sites = 6.9236E-03 seconds\n", - " SEND/RECV source sites = 2.6796E-03 seconds\n", - " Time accumulating tallies = 3.2976E-04 seconds\n", - " Total time for finalization = 7.4870E-06 seconds\n", - " Total time elapsed = 3.1057E+01 seconds\n", - " Calculation Rate (inactive) = 15573.4 neutrons/second\n", - " Calculation Rate (active) = 7201.41 neutrons/second\n", + " Total time for initialization = 2.5229E-02 seconds\n", + " Reading cross sections = 1.1001E-02 seconds\n", + " Total time in simulation = 7.3074E+00 seconds\n", + " Time in transport only = 7.2846E+00 seconds\n", + " Time in inactive batches = 6.2995E-01 seconds\n", + " Time in active batches = 6.6774E+00 seconds\n", + " Time synchronizing fission bank = 4.8069E-03 seconds\n", + " Sampling source sites = 3.3693E-03 seconds\n", + " SEND/RECV source sites = 1.3618E-03 seconds\n", + " Time accumulating tallies = 8.0542E-05 seconds\n", + " Total time for finalization = 2.9260E-06 seconds\n", + " Total time elapsed = 7.3508E+00 seconds\n", + " Calculation Rate (inactive) = 79370.9 neutrons/second\n", + " Calculation Rate (active) = 29951.7 neutrons/second\n", "\n", " ============================> RESULTS <============================\n", "\n", - " k-effective (Collision) = 1.02346 +/- 0.00203\n", - " k-effective (Track-length) = 1.02255 +/- 0.00193\n", - " k-effective (Absorption) = 1.02775 +/- 0.00139\n", - " Combined k-effective = 1.02594 +/- 0.00120\n", + " k-effective (Collision) = 1.02315 +/- 0.00205\n", + " k-effective (Track-length) = 1.02384 +/- 0.00235\n", + " k-effective (Absorption) = 1.02372 +/- 0.00194\n", + " Combined k-effective = 1.02369 +/- 0.00173\n", " Leakage Fraction = 0.00000 +/- 0.00000\n", "\n" ] @@ -1294,8 +1307,8 @@ "output_type": "stream", "text": [ "Continuous-Energy keff = 1.024739\n", - "Multi-Group keff = 1.025941\n", - "bias [pcm]: -120.2\n" + "Multi-Group keff = 1.023689\n", + "bias [pcm]: 105.0\n" ] } ], @@ -1392,7 +1405,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 40, @@ -1401,9 +1414,9 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADDCAYAAACS2+oqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4FdXW/78rCEoLHaQrzS4IylXBggVFQbkgCAe4the8\nYntf/FlBUa/1esFyFSsiXj2iUgQbguBVIArSu7Q0AgQIJSShJuv3x0zOmTlzzqwJSTgn4/o8T57M\nXrNm7X32rFmzZ8/svYmZoSiKolR8kuJdAEVRFKVs0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP\n0ICuKIriExIuoBPRaiK6PN7l+DNDRN8R0ZBSHP82EY0syzL9GSGiIiJq5bLft9eK+uBxwsziH4AA\ngN8BHACQBeBbAF28HCvYnQDg2dLaieef+RsOA8g1/w4AWBbvcnko92gARyxlzgXw/+JdrhKUeQ+A\n+QAuLsHxPwG48wSUMw3AIQB1I+TLABQBaOHRTiGAVhY/K9G1AqAygKcArDfPcaZ57V4b73MZ5Xyq\nD5bBn9hCJ6IRAMYCeA5AQwAtALwFoJd07J+Il5k52fyrycwXlHUGRFSprG0CmGQpczIz/6sc8ihr\nJjFzMoD6AP4L4Mv4FicqDCAVwMBiARGdC6Cquc8rVMpyTIFxnQ4GUAfA6QBeB3BD1MzKx8ck1AfL\nEuFukgzjztnHRacKgNdgtNy3AngVQGVz3xUwWgUjAGSbOreb+4bCuNMdgnG3m27KUwFcZbkbfg5g\noqmzCkBHS95FMFswZtrWijHz2AhgN4CvADQ25S3NY5Oi3TkBtIZxovYB2AngM5ffH7PlZMnnbwDS\nTVtPWPYTgMcAbAKwC8AkALUjjr3TPPa/pvxvMFqAuwCMKq4vAI0A5AOoY7Hf0cyzUoyWxsdSK8Kt\nLsxznQ1gP4AVAM4uyXmwnMO7AWyA0eJ5U2gdfWxJnwWjFVvPTNcG8LVZzhxzu4m57zkAxwAUmL70\nhik/E8AsU38dgH4W+zcAWGPqZwIY4bEVlgrgCQCLLLJXADxulrdFtNYagNsAzIv0b3i4VqKU4RrT\nHxp7KOsj5vk7CKMb9iyzbHthXHO9ovmGS5nvB7DZPA//9Ho+1QdL74NSC/0SACebFRCLUQA6Azgf\nQHtze5Rl/6kAagJoAuB/ALxFRLWY+X0An8I44cnMfHMM+70ABAHUMivnLcu+mK0dIroKwAsAbgHQ\nGEAGjIApHgvgHwB+YObaAJoB+LeLrhe6AGgL4yJ7iojOMOUPALgJwGUw6mcvgHERx14O44RfR0Rn\nwfj9A2H8plrmcWDmbBgXQX/LsYNhOH9hKcoetS6IqDuArgDaMHMtM9+cyIM9nAcAuBFAJxj+09+0\n7QoRVYERTHJg1BtgBKMPATSH8SRZANNfmHkUgHkA7jP97QEiqgbjQvoERmtrAIBxRHSmae8DAEPZ\naI2dC2CuVC4LvwGoSURnEFESgFvNfKRWt8MvS3CtWLkawEJm3u5BdwCAHjCCURKAGQBmAmgAw0c/\nJaK2JShzbxiNiY4AbiaiOz2UwQ31QY8+KAX0egB2M3ORi04AwDPMnMPMOQCeAWB9mXEEwD+YuZCZ\nvweQB+CMKHZiMZ+Zf2DjdvUfGDeOYtwujgCA8cy8gpmPwmgdXUJELTzkeRRASyJqysxHmDlF0H+Y\niPYQ0V7z/wTLPgbwtGlnJYxWRHtz390ARjLzdrOMzwK4xQwAxceOZuaDzHwYhkPOYOZfmfkYjP5R\nKx/DrHvTxkAYdRaLWyPKfWoJ6uIojBv12UREzPyHeVOJxMt5eJGZDzBzJoybUgepzDAulLsA3FLs\nn8y8h5mnMfNhZs4H8CKMG2IsegJIZeaP2WAFjG6Kfub+IwDOIaKazLyfmZe72IrGf2Bc8NfCaHlt\nK+HxpaE+gB3FCSKqY57nfUR0MEL3dWbeZvrYxQCqM/PLzHyMmX8C8A0s3UceeMmsr60wnt7djlUf\nLEMflAJ6DoD6lgATjSYw7njFpJuykI2IG0IBgBpCvlZ2WLYLAJwilMdarvTihFm5OQCaejj2YRh1\ns4iIVhHRHQBARI8T0QEiyiUia0v6FWauy8x1zP93RNizOpn197cEMM105D0A1sJw0kYW/a0RvynT\n8psOwt4imQ7gLCJqCaA7gH3MvNjld34eUe4dUXSi1oV5ob8Jo/WRTUTvEFG08+rlPMSqn5hlhvE+\nZzWAC4t3EFFVInqXiNKIaB+AnwHUJqJYN/6WAC4urn8i2gvj4i+u/74wWm7pRPQTEV3sUq5ofGLa\nux3GzbbcMP2y2DebwajjxsX7mXkvM9eB0QqtEnF4TB8zSYe36yaavch4EIn6YBn6oBQYf4XxBUdv\nF50ss1DWAnptiZTkBVE0CgBUs6Std/dt1nIRUXUYTxxbYfQtItaxzLyTmYcxc1MAf4fxCNSKmV/k\n8Mub4aUsO2DcCHuYjlzs1NUjHpOtdbQdxiNn8W+qav6m4nIfBvAFjFb6YLi3zj0Rqy7MfW8y84UA\nzobx1PVwFBNu56E05doD4wnnaSIqdv6HYHRtXWQ+nhe3jIovpkh/y4TxbsJa/8nMfJ+ZxxJm7g2j\n62E6jLotSRkzYPRR9wAwNYpKPmL7r8OckFdNi29uBTAHwEVEFC2YRgYXq+1tMLoLrLSAcZ17LbP1\n+BYo5ZOJ+qB3H3QN6MycC+MlwFtEdLN59zmJiHoQ0Uum2iQAo4ioPhHVB/AkvAeSbBgvfUqC1RmX\nAQgQURIRXQ/jJWwxnwG4g4jOJ6KTYfSh/cbMmcy8G4aDDjaPvRPGixcjA6JbiKj47r0PxksTt24n\nr+WN5F0ALxQ/+hFRAyK6yeXYyQB6EdHFRFQZwNNRbP4HRouwF8ogoMeqCyK6kIg6E9FJMF6mHUL0\nOop5HkpbNmbeAKOv91FTVNMsSy4R1YWzfiL97RsA7YhosOnXlc3fdaa5HSCiZDbeQRyA8fKrpNwJ\n48VlZDcHACwH0Me8rtrAeHyPRYmuFWaeDaPr4CvzPFU2z9UlcL85LARQQESPmHVyJYxugc9KUOaH\niag2ETUH8CCc/dUlQn3Quw+KXRfMPBbGVyqjYLy5zQAwHOEXpc8BWAyguH94MYDn3UxatsfD6B/a\nQ0RTo+yXjv9fGC8V98Lop5tmKfccGDeXqTCC9+kwXjgUMxTG2/3dMN5UL7DsuwjAQiLKNX/nA8yc\n5lKmR8xH3VzzsXdnjPJGpl+HcdedRUT7AaTAeKkc9VhmXgvjC4LPYbQ6cmGck8MWnRQYTr20FA5r\nzTdWXSQDeB/GVwGpMOrxFYch+Ty41Y8X/gVgqNmYeA1G63E3jLr8LkL3dQD9iCiHiF5j5jwYXVMD\nYNTnNgAvIdwlMQRAqvnoPAzGo7AXQr+BmVOZeWm0fTC+0DgKo1txAowumqh2cHzXyl9hBIxPYFwj\nW2BcJ9YXfpE+dhRGY+AGGPX4JoAhzLzRY5kBw6eXAFgK40OGD4VyRkN90KBEPkjMpe31UOKF+ei4\nD8Zb/nSLfA6AT5n5eC4kRTluiKgIhj9uiXdZ/owk3NB/xR0i6mk+7lYHMAbAyohgfhGAC2C04hVF\n+ROhAb3icTOMx7KtMPr9Q4+ORPQRjG9aHzTf5CvKiUYf+eOIdrkoiqL4BG2hK4qi+ISTysOo+Qnh\nazBuGOOZ+eUoOvpooJQrzFzaya0cqG8riUAs3y7zLhcyRnFugDGXxDYY0+4OYOb1EXps/2S0L4wR\nr1a8TKHygKzykvwbv3ukm6hzQ+Anu2BeX+Aye5mfD44Q7Yz8bqyoM+jGD1z3z+PLRBsZ30WZYeH5\nvsDIcJn73SAPYPxy/G2iTngcngteXG16FKX0vkDLcJm7rvpRNDOfupd5QC+Zb/9ikYyC8XVvMXPk\nzP4xWtZ5cqas89X1okqbm1Y6ZNv7jkDjKWE/3bKztUMnkvYN5ZkRll3aRdTBJ7KjdG79i0O2oe9T\naDfl2VC6behLy9h8Om6oXJ7qsgp+lsv8wYRBtvS4vvMwfIr9Ov4IkYPM7VyEeniVOsX07fLocukM\nYCMzp5vftE6C8SJPUSo66ttKQlMeAb0p7HNBbEXJ5oFQlERFfVtJaMqlD907fS3bO2HMkmvFbV6p\nYiKPicJy+XHop2C0SdoiSIvIi4scshXB9RBZLpc5df9C1/35vN9DPkucsqIi4L/h/DP2/SrbWVhZ\n1tktq3jqctkXRYmLgH3hMu8MrnaoFKxNR8G6DIc8flhnkN4DYLYlvUY+fIUHv8YKWeXnPaLKgbwo\n9VZUhAPB8CBHzm3k1IlgT7KHfrfdHnRmyI6yu1GU66ywCLuD4e64k6LO8xXB7x76U06WVbBFLvPC\nYJotXVTIDtnOKLPjWn17h2NeNTvlEdCzYEzIU0wzhCf2icDa/xyEc1SrY2rjKHgYjd1BruxugfdF\nnVe+iZLXaXZZ+4B8E/qitlzm028scN2/1UMfek7tGLMUXxnOv8UNx0Q7Cw96qOOy6kPfHEPJUmcN\nA9760MuBEvi2tc98NowZdIvx8GDc3kOdT64r61wh96HXjNKHDgA1A+GFjXZ56EOv66EPPf1ND33o\nN8mOUj9KHzoA1A9cE9o+3UMfeso+D/XspQ/9iFzmvwS+jSI7zZZeh6tcbRT3oceiPLpcfgfQhoha\nkjEB/AAYE+YrSkVHfVtJaMq8hc7MhUR0H4wRi8Wfdq0r63wU5USjvq0kOuXSh87MM1GyVYkUpUKg\nvq0kMvF9KTra8inlKgLOi/i08p37ZRvZHjpmf5M/R77xUXm5yJHBkbb06uAanBuwv+CqDnkKlcLW\n8uLqO1DLdf8wvCfayNgUJe5kw1iS2uQ6+kG00/wueRbesVNGiTpYIZ+rpis3OWQFwWxUC4Tl8/+4\nxqGTcJzTNby9LwOobUnnd3XqR+LhnXfR4zeIOtxOtnM013l9TDrIGJAb9vdxDYeJdh7e6Zi51kH+\nbNn3z62+StTJ5WSH7CBXtck/Xfo/op3qt+0SdfJX1Bd16g6O8SrFwsB8+7TwSYcZt+b/ZpN9VN39\nO3QJHfqvKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP\niOvAohajwzOm5Qe3oXrAPoNaRgd5QF6Dm+UZ9r5Af1FnL9cRdXrje1s6iCACEZODPYeHRTtV68kz\n4B1i94FFbfGiaKOwp3MQR7CIEeg5JJQmLhTt8HZ5MEha39NEnam1B4k6WRktncKc+thrkTdt5xx8\n5LAjapQzoyyDdRYQ0MWS/tzD8dd5GDA3T1b55YzOos4V9JtDdlLVIKokh337zkJ5ysEHh8uD3Wiy\n7G+bPpX9bfog5+Rrv2AbLrdUSu+O3zt0ItkJedKxyZf2FXVexBOiTvVdh+yC3Em4c9cAu+wr95lN\nqzV2z0Nb6IqiKD5BA7qiKIpP0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP0ICuKIriE+L6HXpG\nUqoltRM5g1MjNFqJNnZcebqoQ3Plb1+TPhNV0CBg/+b9EHLwf7DLKtE9op2H6o8RdXKPvey6f8w8\neXHnGh13OmTHakzBsHrh72rz3pK/+e0+/CtRZzZ6iTqcJuf11NWPOWSr663BuS3CYxTeo7tFO3Hn\nQst35NvYnvbwHXrhJrmu8HSRqHJprtxm46POvDiPwTnh8QrJsRbvtrBn8imiziJ0E3VOGnSpqNN7\nl3NhloJcoPeu8OIY++vK5dlS6SxRZzg+FHW2Y6So83zW83bBnkpAlv278+8C7vVTHxdhlst+baEr\niqL4BA3oiqIoPkEDuqIoik/QgK4oiuITNKAriqL4BA3oiqIoPkEDuqIoik/QgK4oiuIT4jqwCP97\nXXh7fQ5w5nW23SljOoomkt6TBzzw6/IgjcLzRRUsoAts6dm0B9fSRJtsE7cV7dw5Th7FdGR4Fdf9\n53Vb5bofABbiYofs22oHcGOtV0Lpl+99QLTTh6aKOvfyRlGn8139RJ3PIhYMAYA8fItVuDGU3vVz\nC9FOvKlz2rbQ9pEGe1HFkt57XRPx+KQlHvx6quzXleUxMUCqMy/aD9DrYXnGmvqimRbDd4s6zd/O\nFHVyUE/UWdbgTIcsLXk/ljUILwzTBNscOpH8Zc8KUYcXy/V8YfdrRR1HtK3klA084h4brqIqAP4V\nc7+20BVFUXyCBnRFURSfoAFdURTFJ2hAVxRF8Qka0BVFUXyCBnRFURSfoAFdURTFJ2hAVxRF8Qnl\nMrCIiNIA7AdQBOAoM3eOqjjTcj/ZnwSk2e8vXZouFfMqfIhEnV/wF1Gnc95iUacr23UyOIiubB8I\ns5kGiHZeuec+UechvOmusES+Fz/WabRDtharsBfnhdIvwakTCe+R85pZ50pRpwfmijpvYIFDVogC\n1MK+sKCpaKbc8OzbbjYulHX4bnk1IqyXz8vYr+UVtEbQW05hMAgEwr7dor2c1zfLrxJ1evKPog7e\n87DKUg/ndb82h9Ehc0coTc3llcqQ4aFN+5uscnP3maLOJ5372NIpmzJxaefJNtngDPflrPJPcY93\n5TVStAjAlcy8t5zsK0q8UN9WEpby6nKhcrStKPFEfVtJWMrLMRnAbCL6nYiGllMeihIP1LeVhKW8\nuly6MPN2ImoAw/nXMfP8cspLUU4k6ttKwlIuAZ2Zt5v/dxHRNACdATidfmtfy0HOFxi8TH4xFAzK\nL0XXI0fUST0kz25X+ZSgLZ2SkuLQ+ZXSRTt5RXtEnSAF3RVSRRNY+4dzRsasFPtsd0EI+QDgfDmv\nldWzRZ29HvLagw0OWX5KxO/IjlLHm9YCm9eJ9kuLV9/O62dpvBfa/ZjT64j5BP/wUBh5MkEsXeKs\nT0deUXzN4dv7HCoOlgV3iDq5kl8DwCJZhQ87r9eUxYDxAGVA9TzkJV+uwEo5NiAo55UC+7W3ISVK\nHMiZ5JRtXAdsWg8AWC70qZR5QCeiagCSmDmPiKoD6A7gmajKzaaEt/cHgVr2L0boAjmgBwJevnJ5\nXdTpnCdHyFNqOKd2DQTssqM0Q7Szu0ieijRAzrxsLBks2ljZ6byo8rMDYXkgynS1kfBeOa+6dRqJ\nOj085DUmylcuAFA3EJ6eNH1zF9EO2pZ9b2JJfLvGl++Hto98Ng1VBv41lC5YLk+fG+jkoUDr5fOy\n44x2cl4xfM3m2y/LeSUHThV1ekp+DQB5cl7RvnIBGIHeYTk195DXcjkvHJVjDAIe8sJkh+TSQHNb\nelyG+1dyHU4hzDo19tTa5dFCbwRgGhGxaf9TZp5VDvkoyolGfVtJaMo8oDNzKoAOZW1XUeKN+raS\n6MR3xaL1Cy2JTcD2hbbdXFseEFTpc7lbptWA8aLO/BpdRZ2Td9pXLuFcBu8cYpPVbiSvXLIoSR6L\nMpO7ue5f20keMPLvXOcApmMHp2BObvjdRf/kc0U76+r0EXXOIrn/ekDRRFFn2Yu3OYUr0pGeZulm\ncT65JhznVFoT2t6ZtBUNLem/dJogHp+DF0SdD8+8V9S5EEtEHVwbZUWeHQxMsPj2G7KZ6w7IA8eK\njsmr/6wfdrqoc85EZxdp0hogKTnc312UIue16t02ok77xzaJOvyRnNfAiLFHnM4YOMM+eHLMJPcX\nCK1R23W/fk+rKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qi\nKIpPiOvAor8Xhuft2BjcgLaBGrb9E3LbijYI8sQ5L+AJUSeL5GVwejT8wJbemzwLYxp2t8newP2i\nnXxUE3XexTDX/dN/GSja6H75dIdse9U/0Dg5PJfUEaos2ilkedBEx4XywKJJte4QdYY+8b5DNje4\nC1cFxoXS1942R7SD5rJKefJ1Ua/Q9pdciH5Fr4bSnSrJg336s/vKNQDw6EPyaJ9hY2WdDrOWOWT5\nwULsDYTPe52XDot2TpkhrxD0xhh5xuF7l3wo6rx/2yCHbGHlVOQHwoOShi77VLRz/qUeBg19Japg\nV8Pqos6w2961pbcGF2BywD4v0bJH3OcpaiBMzaMtdEVRFJ+gAV1RFMUnaEBXFEXxCRrQFUVRfIIG\ndEVRFJ+gAV1RFMUnaEBXFEXxCRrQFUVRfEJcBxZ9k9QztF2QxPjDkgaAe5LfEW2MpcfljP77jajC\ng+WFYJdmXmJLB5GKAOwyXGNfgSQaXSbKK81QE/dBGpwv34uP7nf+pkkFjAH7Xwulq9SSB4O8jM9E\nnYHt5GWEKteW87rqXecgph2LGFcdsAwAeflk0U68OT9pZWg7n77GP5LCA40e5ZfE4y/CSlFn7Rh5\nZZ8anCfq1KGDDll1CqKOZUHn9MflRcD/yfKgunvxgahDL8gDi+56IOiQVV3HCPz8ayjNr8mrmeF6\n+Tqa2qCHqNMH34o6/cm+8tcCykQX2maTTb9bGDBY1X23ttAVRVF8ggZ0RVEUn6ABXVEUxSdoQFcU\nRfEJGtAVRVF8ggZ0RVEUn6ABXVEUxSdoQFcURfEJcR1Y1JkXhbYzeTOaW9IA8No8edDQhsuniDqD\nr7hJ1EnLlAdpjMNGW7oAO/BYhKzp7LminXvwtqhzc+EprvsH9fhStPEQxjhkG6rtxoJa9UPpbqPk\n1Yim3i6vCvV5G7mO+9SV8/oi568OWUqNTHDAsgTRPfKAs3izKzdcx4UHa6LAkn6+1kjx+Fu4lqiz\nH7LPnjkuXdThWc7zwlsZ/MWQsGB6fYdOJN3wk5zXl7IPTJtyvajTJ3WmQ5ZUH0hqEfZVPlPOiyaJ\nKujzr+9FnaK75LwGHcy3C/Z8jnFZt9pl9whGOrnv1ha6oiiKT9CAriiK4hM0oCuKovgEDeiKoig+\nQQO6oiiKT9CAriiK4hM0oCuKovgEDeiKoig+4bgHFhHReAA9AWQz8/mmrA6AzwG0BJAGoD8z749l\nY+qsQeHEKsLvswK2/fd2f0UsR3UqEHXO49Wizu17Joo6B/fUtaWD2YzApodtspFtR4l2crieqFMt\n/4jr/ptrTRdtLGbnKIQt2IBaaBdKdzv2m2hnTpsuos6Ah2eIOs/E9IQwd9HPDtk+OowraUtYsExe\nHQkXDJd1YlAWvt06OVzefVV3orYlfRrSxDK8SfKqVvfjTVFnwfALRJ0u1y5zyOhrgHqFB+mswnmi\nnV9xqajTce5zok6z/ltFnWdOe9ghW1l/LTaednYoPXqfHD/4PVEFK8e1EXWaQS7zJ7DHt5S6mbi0\nqf06/uesR11tnI5awMux95emhT4BwHURsscA/MjMZwCYC8DD+nCKknCobysVkuMO6Mw8H8DeCPHN\nAIqbuhMB9D5e+4oSL9S3lYpKWfehN2TmbABg5h0AGpaxfUWJF+rbSsJT3i9F5VmdFKVior6tJBxl\nPdtiNhE1YuZsIjoVwE5X7X/0DW8XFTl2b9jtfFkTycnk/vIQAL7lPFGnME+evTCYb7+GU5YCkdf1\nmlNXiXaS+YCoM+mg+/5F1VJFG/lczSFLS9luSwfXimawOrhL1Nm5TrazykMI/Cp42CFbnHLULsgI\nOg/cshbY4qEQx0+JfDuj7yOhbY7w7SzaLWZWGfL5nYKjos4O3iPqpO9wygzfDrOscRSlCDZ6+Pgg\nuFFUwZbgPlFnZZHTcTNTsux5HZLzYg/lyQjK12tdFIo6K5FpS29IcZ6bvfjBITu0NhWH16UBAHJR\n2TWP0gZ0Mv+KmQHgdhjvYW8D4P4pxpOWqW9/CgLd7G+B23XPgoSXr1xu5OWizlN7+ok6gT13R0gY\ngV5kk6xpK38N0IDlADkgd47r/kO15KlT98WYgvWCQPgrl8DK2aKdHwMNRJ1rlm0QdTbKs5Cid+Bk\nUf7AikBUHRsXlPrhs1S+3WLKP0Pb+4IzUTsQnhK2KaWJmZ8L5804kr4epqvdwHVFnS4bo0+xG+gV\n3k5ud6popwqfK+oE5n0l6vweqC3qFBWdHVV+fiAsDzz0rWiH24oqWBmoKeo0g/wJVw00d8guDdhl\n8x3v4u1cglp4Nyn6bwdK0eVCREEAKQDaEVEGEd0B4CUA1xLRHwCuNtOKUqFQ31YqKsfdQmfmWM2k\na47XpqIkAurbSkUlrisWVe6UG9ouSjuIJEsaAHoLPTYAcDXPlzNq9oyoMna7/GhKhfZ+MmoUBLWx\nX/sv9JQfergviTp0u/OdgpWhy+V8Zra/wiHbjWycB8t7h5fc8wGAzoXRu0FsVJJ/0+hCuZ9xPTm7\nkg4RIY+qh9JvtB8q2nlA1ChfnqJnQ9vzKQtdKbwaVxLLdd4X34g6n1aSHxICX8vdjbghSnlODQJt\nw77dc5jsb5nvOrsUHLwt//Yn+GtRZ9bfnV+NBjcxAr98F0rzDjmvmUny77pu9WZRh36RfbspXWxL\n16G9aEr2lY46wP29Yeso3TZWdOi/oiiKT9CAriiK4hM0oCuKovgEDeiKoig+QQO6oiiKT9CAriiK\n4hM0oCuKovgEDeiKoig+Ia4Di56u93Roe0WN9Whfb41t/zv0d9HGq0VTRJ1v/yKXZfgt8sxR/IF9\nEAAvZHDBEJts99c1RDsBfCrqzBpeyXX/u+P+JtpoHjEZEAAUIQmFlvs4z3XPBwCOXCXPK0KXyxNF\nFb0u53X7g3Mdst38IyZzeJDmoqSLRDvAeA865ceTHB5YlMvf4wfuEUqvn9ZRPL7gNrmuBv1TVAFd\nKPt1UaYzL85hcGbYt4e/M1a08/aIEaLOF2PlycIOkGxn7TunOWRZwTysDYSvv8mV5Doc7Rx75yB/\npjxobiDkyf02Y5wtvR8zMRXX22S7Ud/VxkGc4rpfW+iKoig+QQO6oiiKT9CAriiK4hM0oCuKovgE\nDeiKoig+QQO6oiiKT9CAriiK4hM0oCuKoviEuA4sWoJOoe1MHMYxSxoApk2UFwM+don8E3iKvHIJ\nnpbvbZOfvtGWXlQtC1UCTW2y23MniHam1LpF1Ekd18h1/0e4Q7Tx+7bODhnv/Rxjt90aSi+7qr1o\n57xn5RVb+Cm5jmvly4tj582JsiD1mixsmWMdAbJJtBNvNnzXIZxYvhY7aofT/f76sXj815f1EnUO\n1Y++CLjFnlywAAAJu0lEQVSVZ/CoqPPUiFccMloPUGp4UNKNY+UVlOBh7e42JJ+79lgh6swi52LK\nq2gjKlN41efR7d+RCzReHjRU/eRjos42kldOqwr7gvYHcdgh645ZrjbOQ1PXIUzaQlcURfEJGtAV\nRVF8ggZ0RVEUn6ABXVEUxSdoQFcURfEJGtAVRVF8ggZ0RVEUn6ABXVEUxSfEdWBRCi4JbR/ELmRa\n0gDAteWP/u9vKy/b8u+N8solr46+R9QZgbds6aMIol/EaIpbsuV75Lo68u9qVVjour82TRdt7Gzs\nXP1kSu1j6Nv43lC6HvJFO9lP1RZ16uyX67h9rR9EnQVnXOYUrjsKnHHYImgj2ok3bW5YGdo+sC8D\nNS3pxXSheHzL+mmiztv8mqjzZJo8iAljowwKCwaBQNi3W6G1bKeR7NcdeY2oczedI+qswnkOWXXk\noS4sKyItlQe78V75en2RHxJ1lr43RtRxjIdbl4G0ZV3tdt7pCje6X+2ehbbQFUVRfIIGdEVRFJ+g\nAV1RFMUnaEBXFEXxCRrQFUVRfIIGdEVRFJ+gAV1RFMUnaEBXFEXxCcc9sIiIxgPoCSCbmc83ZaMB\nDAWw01R7gplnxrKx441W4cSShti/u5VdoQNEptJfRZ0DbWuKOiMgDwzgAfbBM5zO4BlDbLKcSdVF\nO2dfLg/mWYEzXfcP9FA5k2igQ7aYNuMohQeJDE+VBwS90WqkqNO41nZRZ8HH14g69w5xrp6zoc4y\ntGu6M5R+y8vSOKWgLHz7bKwNbW9FFppZ0jMWDhDL8HrnYaLOfSd9KOrcc/tEUWfF+DMcsgzkYgWe\nCaWbYLdop9J2eWWfvzcbK+q88+BqUSfrjboOWT4dwXW0PpSey11EOzl1eoo6Ty79l6jz5TDZTr/1\n39oF3xDQM2IwlttyRAAghLLStNAnAHCuAwWMZeaO5l9Mh1eUBEZ9W6mQHHdAZ+b5APZG2SWP/1WU\nBEZ9W6molEcf+n1EtJyIPiAieRVbRak4qG8rCU1ZT841DsCzzMxE9ByAsQDuiqn9Yd/wdlGUiXTk\nbjsczNop6qSy3Cf3PXJFndXpbEun7AYAuywvKPcj1swWVZARdC9PGqWLNg7zKQ7ZlhR75sFd7NCJ\nZE1DD32avE/Uwa9BUWVDpWUO2faUtAhJ1ShHbjT/yo0S+faivq+GtjnStzfLk0Yt3hQ5k5MTZrk+\ng5tFlai+tjzloC1dG0fl8mz2cH43LpXt/CHbmRo84pD9vsB+7WXxLtFOHg7J5UmTyzNvfZaog20R\ndpYtcOpEe712ZC1wdB0AYPnP7lmUaUBnttXg+wC+dj3gzinh7SVBoFPEyy4PL0WrXrFF1Dmdq4k6\nPTBP1Dl/RuTNgxFoaX8KzwnIVVr//cOizopAsuv+5dRStJHHNaLKLwyEX4oGUn8R7axtda6o05jl\nl6IfH5NfZrYLRL8w2gUuCG3PHiK/CAeae9DxTkl9u/OU/wttbw0uQLNA+AXd0kXyS9ELO/9X1Plk\niFyfgdZDRJ1YvtbDIm+CAtHOvYs8nN+L5NbMXA92+gTuiyGvEtpeyw1EOzmoJ+q8uVQuz2UdPxN1\n3lgfxU7PCNlbThVYvrPocBkw65PYPX+l7XIhWPoViehUy74+AOSmnaIkJurbSoWjNJ8tBgFcCaAe\nEWUAGA2gGxF1AFAEIA3A3WVQRkU5oahvKxWV4w7ozBztOWRCKcqiKAmB+rZSUYnrikX4xtIXtIOA\n7RF9Q8fkF3b3X/6mqDPyUXnQUOeXF4k6Z79jf+lWNBk4dou9jI1Gyy9X2/60QtTpRO4vj1qx/Lbr\n+aufdwqzg5g4PhyvPp9zq2hn/uprRR3IC9HgnL8tFnXeWvGwU5gZxOyV1hj7jFMnwZjxlKWffHUR\nlq63pAfJXz8+OPE9OZNOsp2kO+UXsPSF8zrj34IYeVK4zm/qP0m0w/Pl8vznHLlPn9vJPcFNpkf5\nqnRxEPdVD5eZDsnx49b+H8nlOSiXp98N34g6jo88cgB8ESGL/torTLTvASzo0H9FURSfoAFdURTF\nJ2hAVxRF8Qka0BVFUXxC4gT0vLWyToKx7o94l+A4yK949YwtFbDMVnZVwPJvXRfvEpSczApWzwfL\nvryJE9DzK54Dra+QAb3i1TO2VMAyW9ldAcufVQHLXNFuQofKvryJE9AVRVGUUhHX79A7tg1vb94M\ntG4bodBQtnEqmsn5NJXtNIA8NwpV6hgh2Ayq1Nom6thYzquF9DEpgNOEOSaa4KDrfsBev8VsTrXX\ncztpxnwABc45vpzUkVVaQ55T5+QoVbO5EtDaKvdQyUvlOaDKlY5NwtubqwKtLWmc7MGAPMUIhDVQ\nDOT1VkBRzt3mKkBri7yVhxPc0cP1Wi1JDjkF8hQsQJS5LjdXBlpb5CRfZjjdQ0V39FCHaONBp749\nuTkbaB15DoVpntq0AGa57Cdm+eP78oCI4pOx8qeBmeMyf7n6tlLexPLtuAV0RVEUpWzRPnRFURSf\noAFdURTFJyREQCei64loPRFtIKJH410eLxBRGhGtIKJlRCTP7BUHiGg8EWUT0UqLrA4RzSKiP4jo\nh0RaSi1GeUcT0VYiWmr+XR/PMpYE9evyoaL5NXDifDvuAZ2IkgC8CWOV9XMADCQiL+/v400RgCuZ\n+QJm7hzvwsQg2ur1jwH4kZnPADAXwOMnvFSxiVZeABjLzB3Nv5knulDHg/p1uVLR/Bo4Qb4d94AO\noDOAjcyczsxHAUwCcHOcy+QFQmLUX0xirF5/M4CJ5vZEAL1PaKFciFFewLJyUAVC/bqcqGh+DZw4\n306EE9cUQKYlvdWUJToMYDYR/U5EQ+NdmBLQkJmzAYCZd8DT1/5x5z4iWk5EHyTao7QL6tcnloro\n10AZ+3YiBPSKShdm7gjgBgD3ElHXeBfoOEn071bHAWjFzB0A7AAwNs7l8Tvq1yeOMvftRAjoWQBa\nWNLNTFlCw2wsc2+uBj8NxiN2RSCbiBoBoYWPd8a5PK4w8y4OD5Z4H8BF8SxPCVC/PrFUKL8Gyse3\nEyGg/w6gDRG1JKIqAAYAmBHnMrlCRNWIqIa5XR1AdyTuKvC21eth1O3t5vZtAKaf6AIJ2MprXpzF\n9EHi1nMk6tflS0Xza+AE+HZ81xQFwMyFRHQfjCkKkgCMZ+ZEnzatEYBp5hDvkwB8ysxuUyzEhRir\n178E4EsiuhNAOoD+8SuhnRjl7UZEHWB8fZEG4O64FbAEqF+XHxXNr4ET59s69F9RFMUnJEKXi6Io\nilIGaEBXFEXxCRrQFUVRfIIGdEVRFJ+gAV1RFMUnaEBXFEXxCRrQFUVRfIIGdEVRFJ/w/wGvJaRK\nz6tq0gAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAEPCAYAAADf8cexAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XecFdX5x/HPAyIICliwBCuCiTFqBKOJBey9xZgEsP5i\n7LEQo8aKJTFRoxgTTdQYFQuIorFEwRoVuy4WFAsIYqMoShWF5fz+OHPl7uzd3XmW3Xt3me/79bov\n2LnPnXPmzpxznznTLISAiIiI5FebSldAREREKkvJgIiISM4pGRAREck5JQMiIiI5p2RAREQk55QM\niIiI5JySARERkZxTMiAiIpJzSgZERERybplJBszsf2b2RKXrIS2LmY0xs4ebYb63mtl7TT1fKS8z\nO9/MFmeMPcLMFpvZus1dr2WBme2cfF/bNPF8N0zmO7Ap55t3jUoGzKyHmV1rZhPN7Cszm5V0uieZ\nWYemrmRRuRub2eA6GmMAMjXqZYWZTU4aRanXg5WuX3Mys371LPvtRaEheTW1imxvRR1h4VVtZp+b\n2QNmttVSzHfbpG2t2JT1XRpmdnjRcpb8QTGzD5P372tkMbW2DzM708z2zxJbH4sOM7OHzWyGmX1j\nZtPMbLSZHWVmyzeyzhWX7HyVanvVZrZRUWhz3e++IvfRT3YCipd3gZm9nbSdRq9PMzvbzPZtyrp6\nLef9gJntBdwJLACGAuOA5YHtgEuB7wPHNmEdi30fGAw8AUxJvbdrM5XZkgVgLPAXwFLvfVL+6lTE\nlcDLqWmTi/6/I83TcRxB7e+8nG4BRgNtge8CJwCPm9mWIYS3GzG/7YDzgOuBuU1Wy6bxFTAQeLZ4\nopn1A7oT+6KmdBaxj7s3NX0oMCyE8E1DM0h2iv4D7AY8A1wGTANWAfoBVwNbAUc1XbXLKgAfAr+n\njr4nhPCYma2Q5ftyFRzCxOaYb9bigXnA0cTl7gocQPxdWh/4v0bO9xxim75/6avYOK5kwMzWB4YD\nk4CdQgjTi97+h5mdC+zdZLUrUQXq6NhDCIuasdyW7OMQwrBKVwJiBxhCaOqOuSFjQgh31/Vmc20X\nIYTq5pivwyshhG9HQMzsOWJHcixwSiPmV8nEpiEPAj83s5NCCMWjMQOJieBq5ahEiE91y/oDdCVx\nB+WkEMLfU+8NMbMNaWAHxszaAm1CCAvdlS2PWQ31Pc31g12hRKBgYWq5rzGzF4BDzezUEMLMSlVs\nqYQQMr+AfwDVwNYZ49sC5wITiNn7JOAPwPKpuMnAfcC2wAvEPYGJwKFFMYcTh2Wrk38L/++bvP8/\n4PGi+H5JzM+Bs4lZ7FfAo8CGJcr/d4n615hnMq0bcAMwNZnfq8BhqZhC2X1T09dLph9WNG0N4Mak\nfguIWfV/gHUzfL+TgPsyxN0EzAG+k8x7DjCduLdiqVgj/piMS5ZvKvBPoGsd62w34KWk7icl73UA\nrgJmALOTMr+TLPt5ScyOyd/7l6jvwOS9Orezou/4wAaWfQzwcGraycCbxAx/JvAi8POi9zsn9Z+c\nLNc04l74pkUxtwLvpea7IjCkaF2OB04p0SYWA1cABybf8wLgDWCXDOtyw+TzJ6Wmd06m35+avjlw\nM/B+sj4/Je79r1wUcxGl29Z3Uu3vZWA+8DlwW/H7ScxGwN0saRtTkrhOnn4mVWZ18j0tAnYveq9d\nUo9TSLUDfO1vMFBd9Hf6e1hM0jcQR4MW00DbBNYmJg0POJa1ULffJtvnBGAhsFkz9Ts3EfuBDZJt\ney7wMXBuxvo+AbzeQMzOSbnbeLYRYA9iu/0yqePbwIUl2sDAVHm7EEdh5gFfJOVslIr5Q/LZ9Ykj\nPV8msdcD7TMs9y3AzBLTr0i2mz6p6Wckdfqc2HZeAg4o0R+kt7nrimK6J+trKkv6isNL1KHefq2h\nl/cwwT7A+yGEFzLG3wAcBowgDmVvTRyC2xj4WVFcAHoRh+ZuSBb8V8CNZvZyCGE88BSxgz6RuEIL\nQ6Hji+ZRyu+JX/RlQBfiyrkV+Emq/FLSxxI7EBOEDYG/EX8sfg7cZGZdQgh/yzDPtLuJ38dVwAfA\n6sQ9hnWpfSiklHZmtmqJ6fPCkr30QDw/ZDTwPHAqseH8ltjpXFv0ueuI6+zfwF+JncWJwA/NbNuw\nZI84AN8Dbk8+fx3wTvLezcBBxMb2ArGT+i9F30kI4QkzmwIcTO3h2IOBCRm3s5VKLP/MkLQOaq/D\n44g/2MOSf1cANiNum3cmYdcD+xLX8dvEPc/tiOvpjaL5hqL5WrKM2yaffx3YE7jCzNYKIZyRquMO\nxG3nGmJHfAow0szWDSHMyrDcaRsk/36Rmr47cVsq/JD8ADgmWZbtkpgRQE/gF8BviB0kxA4FMxtM\nPIRwe7JsqxM7nq3MbIsQwlwzaw88TNzOriQmUGsTv8fOxA6qsSYTt9sBxG0YYK9kvsOTuqQ19tDQ\nIcTv6gXiNg1xx6Qwzyzz3ZPYyd/WiPJ/BbQntqmvgZnN1O8U+oRRwHPAacQf4QvMrG0I4fwM82hb\nou0tCCEUr+viNtLgNmJmmxL7g1eIQ+dfE38b6j0J0cx2Bx4A3iXugHYibhfPJNvoR0X1CcBIYt93\nBrAl8Xufmny2MepqfyclZd1KPJw+kNjO9wwhPBxCqDazQ4g7hGOI2x5J3TCzNYk/6t8QfyM+J277\nN5pZpxDCNUlcln6tfo7MdSVixnJ3xvjNkvh/pqZfSvxx7lc0bVIyrTiDXI2YOV5aNO1nFI0GlMhU\nS40MjAPaFk0/MZnH91PllxoZSM/z5OSz/VOZ3TPALJLsNim7Vj1JZejE5GQx8Nus6yE1v0nUzCaL\n9+pOL4q7MZl2VurzrwAvFv29XfL5X6bidk2m90+VXU1qbxbYIon9S2r6v5P484qm/ZGYLa+UWu/f\n0MAeStH6TWfU1RTtuQFPUzQyQBxKr2pg3rOBKxqIuQV4N7VtLgZ+l4obSdzDW7doe1mcLHdxPQvf\n29ENlFvYKzoTWJX4w7w9ca+9Gtg3FV9rb4eYbNUY4SN2ijVGA5LpPYh75aempm+aLNfvkr/7JPXa\nt776O7fvwshAb+B4YpLSPnnvDuDRom0xPTLQYPtLptUYGUimzaF0f1CoT0MjA5cncZumprdL1lnh\ntUqJun1RPD15r0n7nWRaoU8Ykoq9n9jvrtLAMj5B6X7n30UxO1PUr2fZRog7KtUU9Qn1tIGBRdPe\nII5sFPclP0zmdX3RtMIo2DWped4LfJJhm7wlWUeFddgDOD0p5+US8e1Tfy9H3Ht/KDX9K4pGA4qm\n30TcKeySmj4C+AxoV7Te6u3XGnp5ribonPw7J2P8XsQMbEhq+uXEoej0uQVvhRC+PUEohPAZcU+z\nh6OOpfw71Dy++3RSfmPmuycwNYQwvDAhmfdVxCHifs75fUX84dvBzLo2oj4Q95h2Ju7pF167EjPE\ntGtTfz9Nze/hIGKH+5iZrVp4EU9SnEsc2i82KYTwaGraHsT1/o/U9L9R+7j0UOIhhYOKpvXHt1d1\nAbWXfWo98V8C65rZFvXEzAJ+nGTlWe1JXJdXp6ZfQVyePVLTR4UQvh35CSGMJe49Z90u/0A8DDMV\neJK4Z39yCKHGCUghhK8L/zez9sn6fIG4LnpnKKcwgjcytU18Sjz0UNgmCqMJezbTFUUjgI7APskV\nD/vQuD3vcij0lekTMfcirrPCa3KJz94Vah9zbup+p1h6e/07cQ92lwyfnUTNvmdX4s5eXbJsI4WY\nn2YoHwAzWxvYBLghhPDt71MI4VXgcWr/1gRK94VrZNx2u7BkHU4A/kwcuT4wHZhqf12JJxyOIUPb\nS0Ybf0pMVJZLtb+HgZWJCQ9k69fq5TlMMDv5d6WM8YVsdELxxBDCNDP7Mnm/WKkh8S+IC7w0Piwx\nTxo53/WAUteWjyd2rullqlcI4RszO4N4CGWamT1PHOoaGkKYBmBmnYlDPgXfhBCKh6I+CyE8kaG4\nBSGEz1PT0t9vL+LGOp3aAnEvtNikEnGF9Z5+b0I6MITwjpm9RNxTvTGZPBB4PoTwfqmFKGFcCOHx\njLEAfyIO0b9i8T4BDwO3hRCeL4o5jTiS8ZGZvUw8gW1oCGFyPfNdD/gohPBVavr4oveLpbdLiA06\n63b5D+IhphWIHfJviElHDUnHcT7xEEC3orcCsVNrSE/isG6p9RFI+oUQz/D+K3FY9HAze4p4Tsmt\nxR10Y4UQPjOzR4nbR6ekTnct7XyXRj1ts7C86cs0x7DkR/Z0Sg99Ty4xrUn7nSKLqb1e33XMc17G\nvgfIvI3czpJDxJcRz/G6mzgiXdchkEJd3y3x3nhgJzNrF2qeiJn+vSn+Xfi0gUWZC+xP3AbXJo6q\nrU7cuavBzPYjHhrfnHj4pyDLCZBrEn9vjydeLZRW3Cdn6dfqlTkZCCHMMbNPiMODWRT2AutagWl1\nnZ29tGc5Z5lvXXVsSxwi9dalvvnVDAzhr8k10gcQj+9eCJxpZjuGEF4jHrc/vOgj/wN2yliPYlnO\nfm9DPI43kNLLOiP1d62Nvx51fSdDgSvN7DvEjvXHxI2/WYQQ3jKz7xL3LPcg7vmeYGbnhhD+mMQM\nN7MniVn5rsTk4Awz27/ESEiBdztd2u393aIk6L9xJ4K/mNn/ku2mYCRxePYS4nkM84jD1Q+S7T4j\nbYhtID2yUVC8JzbIzG4gdpS7EfcyzzCzH4cQ6hutyapwzsJaxGHWupKMzO1vKdXVNt8mrscfsOQc\nE5Jk/HEAMzu0jnmWalNN3u/Uo1mvKmloGwkhfGVm2xFHnPYmbncDiD9udW2Djanz0rS/RcVJUJKk\njicm6AcVTd8RuIe4zo8ljuItJF5OWnzOXF0K7fNm4jkHpbwG2fq1hnhPIHwAOMrMtg4Nn9w1mbgw\nvVhyYhlmtjpx7/MDZ9nQfDea+IJYp7T1WHLyEMRlKpUMbZz8W1imL1hyDWqx9UsVHkKYRDycUrjk\n6DXisbPDiJ34Lam6NpeJxL3MZ4uHt5w+IK73Daj53W1UOpxhxKH0AcRh4G+IQ8LNJoQwPyljhJm1\nIw7DnWtmfy4cUgohfEo8ue8aM+tGXCdnEfdUSpkMbJdc/1zcoae3jeZyEXBk8u9+8O2oQF/gzBDC\nJYVAM/teic/X1bYmEn9M3m9gZCTOJIRxxPN0/ph06k8Rr8m+MPOS1O0e4vDu1sAv64lztb8SsvYz\ndbXNh4g/NgdT+nCd12Saod8httMe1By1K7TTZtteG9pGkhGAx5PXqckl6+ebWd8QwlMlZjk5+fe7\nJd77HjAtNOPlmSGEj5MRj7PMrHcIoSp560Bi8r1H8aFqMzum1GxKTJuafL5NltHPLP1afbx3ILyU\neOLTv5If9Ros3h3tpOTPB1lymVqxU4kL/l9n2RC/mFIb+9KaSDxG/G1ylNwNap1U3IPAmmb2y6K4\ntsSTEucQj91CbEjVxI642PHUPLt2heQM22KTknm1BwghvB1CeLzoNbaRy5jFCGKCeF76DTNra2ZZ\nhpVHE9dReu/+REps8Mnx0YeAQ4md56gSx0ybjJmtkip/IXFPrg3xyoy2ZrZSKmYGcegwva6KPUg8\n1ppe7kHEbeGhpax6vZLh6euBvc1sk2RyoQNIt/NB1F4XhTPA021rZBI7uFS5he/TzDqbWbqcccln\n6/veMgvxLPVjiYc96rs5S6b2V495ZOhj6mqbIYQPiYeZ9jSzUsO74Ot7m7TfSflNib+/AR5z1C+T\nLNtIun0mCiNdJbejEK8UGAccUdx2zWxz4kjNA0tZ9Sz+Srzy4fdF0wonN387MmNmPYhXT6TV2uaS\nH/B7gF+Y2cbpD5jZakX/r7dfy7IArpGBEML7Fu8HPRwYb2bFdyDchni5y41J7OtmdjNwtJmtTNxg\ntybu7d4dQniyVBkNeJX4BZ+RnIzxNfBYcrLh0vgXcXhntJmNIJ6tegi1j3NfR7ws6yYz25Ill/j8\nhHjy1jyAEMJsM7sTOCkZvp1I3ADSN0fZiHiy3gjgLeJw7IHE40BZ9yi6m9nBJabPDSGkL9mrVwjh\nKTO7Fvi9mf2QODS3MKnnQcRjfXXe4CeZR5WZjQROSTbW54knOPUqhJT42FDi8d9AvJyoOT1u8ZLG\n54iHRDYhdpb3hhAWJHvTk5L19waxke5GPFHnpDrmCbHRPgVcYmY9WXJp4d7AZckPRHO7kvgDcQbx\nzPEvzexZ4mGnFYj3sNiDeKlhejj0lWTan5JlXwj8J4TwnsVLCy9MRq3uIx4z7UE8jPI34olsuxJH\ntu4kHt9uRxxCX0gD20wDatQzhHBLXYFFMVnbX11eAXYxs0HE72xSCOFFV63jTtD6wFVm1p+YvExP\n6rBtUp/xdX66pqbudwq+BvZI+unniSc47gn8scT5RY1VvP7q20ZGJjEXmNmPicnzB8Tj5scn/69x\nB8qU3xF/9J8zs38Tz9c4kXh5bFOMStUrOaflZuDXZtYzhDAhqc9JxN+VYcTDW8cTR8o3Sc3iFWA3\nMzuFuOMxMYTwMvHckr7Ai2Z2PXGbWYV4OeT2xO8HGujXsi6E+0X8sfwncWP7inji01PAcSSXOiRx\nbYide+GmQ5OJw5jtUvN7P6l0qctXHktN+xVxQ/qGmjcdqhHLkstsDkx9fr1kevqGHacQTyqZT0xc\ntqij/NWIycM0ltz849ASdV+VuKc9h3gJyNXEYb1vy05W6lXES01mEzfcZ9N1rmc9TErmV+r1flHc\njcS7haU/P5h4/Cs9/Ujita1zk3X7KnAxsEZD6yx5r/imQ7OIP/Q9iVnyaSXi2yXf0RekbkhVz7KX\nXL8l4p4GRhf9fUyyfqcn6/pd4iWOhcuzlicO/45Nln0WsaEemZrvLcA7qWmdiIc8PiJu728TO+vi\nmLZJvS8vUdcpwLUZ2l41cGId7w8ldvKFSxm7EzvamSy5WdBayTzOTH32XOKJjQupfdOhA4ltfHby\nepOYfPRI3u9BbBfvEROo6cAjlLgM2NHPHJ7Uo3cDcbW2RTK0v7raADH5fSLZ/r+9XI6MlxYWzceI\nOz+PJG3ha2K/8TDw6+JtnSX90qA65tVk/U5RnzCbmLCMSuI/wXfTodcaiElfWtjgNkLck7+HJTeJ\nm5Js0xuUaAPpmw7tTGzvc4l9yUigVyrmouSznVPTj0xv83Us0y3A53W815PYdq5LzfcdYl8zjriT\neRHxZNPiz36PeM5JYZsrnkc34rkVk4n9ysfEEdjDi2Lq7deyvCyZkUizSkYaqoCDQ+oWpsmQ5yfE\nDv3oStRPJE/M7EbgZyGEzg0GSy4sM48wlpajxHkQEEdeqol7l2k/Je75DG3OeomISGnupxaKZHC6\nmfUhDnstIh6L3J04BP5xIcjiI3c3Jx5KqgohjKlAXUVEck/JgDSH54gnC51DPJFnCvHY7MWpuOOI\nVxCMpfGP/hSRxtExYvmWzhkQERHJuRY/MpBc6rU7S86kFJHG6UA8e3x0aLpLx5qc2rxIk8rU7lt8\nMkDsFFrqA0lEWqODibf2banU5kWaXr3tviLJQHJXrt8Rb5jwGvGa6ZfqCJ8c/7mVeClmsUHUfigi\nxHsieQxwxgMn+g6vXHW4/4q5k866rvQbYwfBFrWX+/iLr3CXcc2Y37ri99z+Hlf82OB/iNbUMeuX\nfuP6QXBU7eXeZTv/DcYe/U+pm4DVo6FHl5TiPQL3ZB0f+HQQrFV7uX84oqE7gqeqM/4TXjvkj1D6\nYTjNztHuJwPQ61bomLrx2qRBsEGpNk+8itvD+3ifxmwDx/g2gm6H1H1vqlmDLqbLkLNqTZ/x9/SN\nUuu35omTXfEAR3KD7wPObf/At+q+QeegITBkUIkiFvvK+N+mP/Z9AFjk/Ilsj/8u7q+zWcnpjw56\nhF2G7Fprer+SF2TV7YvxfRh0yHBooN2XPRlIbql5OfFe1C8Sf9FHm9lGofSdBJNhwu9R+6mPXUtM\ng/iYb48sT3JN6e7b2nv2zvqwxyKr1FGvdl1Lvte9d607RDdsqm/ZV+ld1XBQkeVD+kZbGUwtdYtx\noFNX6Fm7viv39v4KAGMbsc69vMnACnV8oG1XWKF2fVfs7btrc1jS3Ms+9O5s97F+HTeGFVPL3bZL\n7WnfFuKs1AoNh9Qs2xkPsLZvI1i+d92X/bfpuhLL9y7Rntbs6Sqjfe+OrniA9VjV9wHnuWi96/lh\n77Ii9C7xNA1vMvBhPd9tXRZmu5Pvtzq4nt0WTWetktPbd2nPmr1rv/ddfOvvsyUPm6233VfiPgOD\niJeYDQ0hvE281/h84p0FRWTZpHYv0oKVNRlInqTUh6KHYIR4OcOjxPtsi8gyRu1epOUr98jAasTB\ntmmp6dNY8sAFEVm2qN2LtHAt5WoCo8EjrIOo/VTRdZupOi3cev0rXYPK6JfT5e7iX+4Zw55gxrAn\nakwLs+Y3VY2aSv3tftKgeI5Asfbes/6WHSv036fSVaiIAbtVugaV8f0B/vOtHh72BY8M+7LGtG9m\nzc702XInA58R70+/Rmr66tTea0gZQqNO9FsWrdeIqx+WBf1yutxd/cvdbcCOdBuwY41poWoSz/Q5\npqlq5dG4dr/BkLpPFsyhjgNymgzsXukaVMYmA37g/sxuA1ZmtwEr15j2WVU/9u1zVYOfLethghDC\nQuLjYHcuTLP44O2dqf9Z1SLSSqndi7R8lThMcAVws5m9wpJLjDoCN1WgLiJSHmr3Ii1Y2ZOBEMII\nM1sNuJA4bPgqsHsIYUa56yIi5aF2L9KyVeQEwhDCNcA1lShbRCpD7V6k5WopVxM07GiDtTLeYuyf\nJ/rmPa0RT2583ne7s73PeNxdxNm3n+2K78Q8dxnVG/puqzaVLg0HFTmaOm6pXI8pE+q4A2EddrfR\n7jLWObLu276WcsXIc9xl8Jpvu+r++gRX/Jh3dnHFd57iu3tkxR0HbOSIP805/16+9dP3jUZsZ8G3\nnd124K/dZeC7ASHLB/8tc399yK2u+Ktv891L6q99jnLFA5y85fWu+P3/8rC7jDd36uGK/w8HuMs4\nJVzpiu/24lxXfNVn2a5KqMQdCEVERKQFUTIgIiKSc0oGREREck7JgIiISM4pGRAREck5JQMiIiI5\np2RAREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzrWaBxWtefQklu/dMVPslB/6\nHnTTbf8p7vqM4Beu+C/Cyu4yDuAhV/wf3E9qgRVWnemKXxB8DyrqxZ9c8QDV+/genmSh2l1G+NRX\nxuSfre8u4+6uB7viP56yniu++0a+Bxu1n/sRs12fqLA7A6zseJjQK8/45v+7bV3hz3/+Y9/8gSce\n3ssVf8vd/gcVbYLvAVTjJv3IXQa3LnaFn3Cfbz/zin2Pc8UDhFd8dbIN/fu+m2z3vit+ztDH3GU8\nbLu54rfc6mVX/MyqlTLFaWRAREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZ\nEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIBERGRnGs1zyaYuuUnQLZnE0AP37x32MBdH3vcdz/8\nNsPcRdBtoO+ZCW3Nf3/vU1e73BU/e9ElrvjLn17kigdYsfd0V/zcq33PGQDY7fj/uOIfYV93GWGy\nr17n7fx7V/x1dowrvq21c8VX3G4GPSx7/Hu+Zw3YUF91FvzS91wOgGcOcTxbAWg3YHN3GftzgCt+\n9XU/cJdxsvO5J+33PcEVP93WcMUDPBl8z4qYPmEfdxk/t/tc8euwmruMlcIcV/xGX/h+F+YuyDZ/\njQyIiIjknJIBERGRnFMyICIiknNKBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadkQEREJOeUDIiI\niOSckgEREZGcUzIgIiKSc0oGREREcq7VPKiIgT+B1XtnCn328mxxBW2u8z1MBCD81fcQmurN3EXw\njG3hip8QernL+NU1vicofXP88q74TXd8wxUP8AK+B5BccsJJ7jIOtLtd8SeE99xlbHXkz13xwxjo\nip/x5Lqu+M4TP3PFV9zMAJ0cbfMuXzvuvuEEV/yT4VBXPEC/3/niw0nj3GXs8udHXfEbdXzHXcY5\np/seaFZ16cau+N53jnfFAyzYyxdf3da/7+vt51c/1v9bclGXc13xG67s2247d1gPuL3BOI0MiIiI\n5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIB\nERGRnFMyICIiknOt59kEz7aBDtlyl227V7lmXX2quavzFFu74rea+7K7jO2C7zMTrb+7jMuO+40r\n/lT+7ivgFX+++fs+g13xf8YXDxBm+uo1auUd3GXsyeOu+Kt4xldAd184s5zxlbY20DN7eJfu01yz\nn4LvWR53WTdXPACb+voWO6TaXcSsZ333z991m3vdZXDpYlf4c/ZrV/waP/OtO4Du5nzWxh6N2Pft\n5Ft/7U73fU8AR/EDV/xP7R5X/F42JVNcWUcGzGywmS1Ovd4qZx1EpLzU7kVavkqMDIwDdgYKKdei\nCtRBRMpL7V6kBatEMrAohDCjAuWKSOWo3Yu0YJU4gbCXmX1sZhPN7FYzW6cCdRCR8lK7F2nByp0M\nPA8cAewOHAtsADxlZp3KXA8RKR+1e5EWrqyHCUIIo4v+HGdmLwIfAL8Abqz3w9MGQdsuNad1HhBf\nIlLT/cPggeE1Js2f82VFqtLodn/tIOjUtea0HfrDjmrzIqXMH/YA84c9UGPaqFlfZfpsRS8tDCHM\nMrN3yXIB0RpDoEPv5q+UyLJg3wHxVaTjW1XM3n/LClVoiczt/pgh0EttXiSrjgP2oeOAfWpM26Nq\nCtf12bXBz1b0pkNmtiKwIfBpJeshIuWjdi/S8pT7PgOXmVlfM1vPzLYB7iFeYjSsnPUQkfJRuxdp\n+cp9mGBt4HZgVWAGMAb4cQjh8zLXQ0TKR+1epIUr9wmEOvNHJGfU7kVaPj2oSEREJOdaz4OKJr8B\nLMwUGrr6HiLU9g7/wyV69L/BFT9mxe3cZbSf7nsASdc1Gj5jNO3FNlu54keFHV3xb/U5zhUP8LfZ\nvocn/aJeMDN0AAATnklEQVSz70EfAONXPtAVv7GNd5fRf/HNrvixfzrcV8BdvvBWl/r/G1gxe3jo\n53uoTNs7giv+wV/OccUD2EW+MsInvjYP0GZrXxmP3rCfu4ywju+7fXu3P7vi131/uiseYNEDvu+q\n+g53EbS7x/fdLj7Hv/5+tLrz4VRdGw4pNteybbetrXsQERGRJqZkQEREJOeUDIiIiOSckgEREZGc\nUzIgIiKSc0oGREREck7JgIiISM4pGRAREck5JQMiIiI5p2RAREQk55QMiIiI5FyreTbBz156jW69\nsz3x9MbZvVzzNnz3nwa4mLNc8R9bd3cZe67+L1f8VZzoLmMeHV3x13K0K/7ep/wPrNut772u+G+s\nnbuM6uC7h3jvF/zPJhje5f9c8Ueddb0rftfDH3PF8+Y3sIfvIxW1rUH37PfEn335Gq7Z73jWf13x\nT9HXFQ+wzZvPuuJXOnORu4wt+z3piv97P38/sfWzr7viZ9jqrvgnN/Q9IwXAhvr67SEnn+Au47TZ\nV7viJ/zB38/vx3BXfHu+ccX3qsrWhjQyICIiknNKBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadk\nQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oGREREck7JgIiISM4pGRAREcm5VvOgoqfbbM/ybTbJ\nFHtc53+65n2Fnemv0P8ecIWHQ7I/cKWg6sOf+D6wS5W7jG1v9j28w75T7YoP8/z55sJZvu9q+S6+\nOgFcwjBX/ICN7nKX0a6rr147Xet7eBKXtPfFL7+8L77SngywYvaH0fx1jO8hWidynSt+OL4HSQGs\n+D3ntjlhsbuMwbaLK35qWNNdBtu86goftqmv3V/8pr9/3K7a912ddpq/Lxpzha9e2500xV3GHhzh\nij/206Gu+Kq5A7kgQ5xGBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadkQEREJOeUDIiIiOSckgER\nEZGcUzIgIiKSc0oGREREck7JgIiISM4pGRAREcm5VvNsgk3Cm6wc5mWKvfJp37MG3u070l2fQ/rt\n54qf/OEG7jKu4T1XfPdHHneXcRz/cMXvX93BFX/wnne64gFO5XJX/I7nOO/pD9x9RPZ73gPc0dO3\nvgEOXMVXrxGf/9RXwHG+Z3DQeaovvtLeWAQszBx+ynJbuGa/2iLf9z3gjvtc8QAM8oVfEk52F3HW\nlAdd8U+t63zmCcDlvm25zd98sz/7PV97BAi3+up01mWD3WX86cQsd/Vf4qq2/r7opBd9yz7/B77n\nJSxcIVucRgZERERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIBERGRnFMyICIi\nknNKBkRERHJOyYCIiEjOKRkQERHJuVbzbIInntsTPuudKfaE3S5zzbuTzXfXZ9MwzhV/xMyb3WV8\nNXMVV/zZvc5xl/F5WNUV33HeN674/bvc64oHeDn0ccXvuOh5dxmP9dzWFd//NP996S+Y5Ys/0p70\nfWDsXb74d6qg/4W+z1TQ1i89Q+fe2Z+n8Ej341zzf8m+dMUP6ODflu86YW9X/BnXOW/qD6x59DRX\n/NP0dZex7Y/GuuLv7burK37/Xo+44gHsYN89/f+0wPecAYDbjvHFn+SsE8D1fQ5xxZ/DH1zx+7V9\nH7i9wbgmHRkws+3N7D4z+9jMFptZrae7mNmFZvaJmc03s0fMrGdT1kFEykvtXqT1a+rDBJ2AV4ET\ngFopkpmdAfwGOAbYCpgHjDaz5Zu4HiJSPmr3Iq1ckx4mCCGMAkYBmFmp5yyeDFwUQrg/iTkMmAYc\nAIxoyrqISHmo3Yu0fmU7gdDMNgDWBB4rTAshzAZeABrxgG0RaenU7kVah3JeTbAmcQgxfbbLtOQ9\nEVn2qN2LtAIt4dJCo8RxRhFZpqndi7Qg5by0cCqxA1iDmnsJqwMNX7dy7SDo1LXmtB36w44Dmq6G\nIsuKh4bBqOE1Js2f57uUrok0ut2/89sbWK5LpxrT1uy/PWsN8F8aJ5IHC4bdy9fD768x7ZEv52X6\nbNmSgRDCJDObCuwMvA5gZp2BrYGrG5zBMUOgV7b7DIjk3p4D4qtIx3eqmN1/y7JWY2na/XevOJLO\nvTds/kqKLCM6DNifDgP2rzFt16r3uWHLnRr8bJMmA2bWCehJ3BMA6GFmmwMzQwgfAlcC55jZBGAy\ncBHwEeC/m4eItAhq9yKtX1OPDGwJPEE8FhiAy5PpNwO/CiFcamYdgWuBrsDTwJ4hBN9t7USkJVG7\nF2nlmvo+A0/SwEmJIYTzgfObslwRqRy1e5HWryVcTSAiIiIV1GoeVLTcxnOxzWdnij3AeShy5zDG\nX6G1fQ+9uOJT30OHAKy62hV/8T7+3C78rNQN4+pmRyx2xR/1qr9Oozbv5/vAn311Atiqur3vA219\n3xPAYOf6e9s2cMVftflRrvhp1fP4o+sTlTU+bMxyYdPM8T0+ftM1/ys40xW/cMezXfEAM1jN94Gj\n/dvy4ff72tiT+27tLoO+vnqNa+ur037H+NsXT/jal23l74sGXuKs1+98dQI4yrn+xu63hSt+fVsx\nU5xGBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadkQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oG\nREREck7JgIiISM4pGRAREcm5VvNsgqO7Xkf3VVfPFPtPO9Y17yGLR7rr81/n7b2PPyi4ywj/auuK\n/+z+bPegLjaQ21zxDx/vq9O11xzmigdYhw9d8eFxX50Avtmpoyve+i50l7H4r756HXHy4674F9v8\nyBXfufOrwHDXZypp9vRV4KNsbR7g1+v+y1fAx5u5wtv91zd7gOP2udkV/9Y6Q91lbDzDF7/l1y+5\ny7igo29bHryLb/72biP6x4t9dQrvu4ugzURfvRbf5O+LLj7it674NZjmiu+UcZ9fIwMiIiI5p2RA\nREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTn\nlAyIiIjknJIBERGRnGs1Dyoaz8ZMZf1MsffcPNA170U/8X8NYeRi3wfO9+ddd52/tyv+iNk3ussY\n2eUgV/yka9Zwxd/E/7niAV76ZCtX/NidNneXsemFE13x4Tzn+ga6zPM9QWbuY92cJUxwxn/mjK+w\nG5eDbu0yh1/668Gu2S/cIvu8AS4/+ixXPIDd52v33/+d/4E9jHOG9/2Bu4jB4193xc9xNsnO2/ji\nAbjV913N/9rcRXT8R7Ur3ob7+/nVbborfg4rueLbsWqmOI0MiIiI5JySARERkZxTMiAiIpJzSgZE\nRERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIBERGRnFMyICIiknOt5tkEr7Mp\n7ch2T+3Q1XcP6hN7Xequz9/ea+uKHzL4OHcZv+VqV/xB0/y53fiVfd9Vj2rfvbq72r2ueIDpa63m\nil+Vee4ypp3X1RW/8izf+gbYvMtoV/wz393eWUJPZ/xsZ3yFvWHQ0bF9Puab/ZDbfc8auPJ5/3MD\nzjv+dFf84P3+7C4D5/3w7+Gn7iK23uhVV/zbX23mir8z/NwVD3CpneuK73hhI/Z9L3a2+738RRzz\n2lDfB+b6wgd+XpUpTiMDIiIiOadkQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oGREREck7JgIiI\nSM4pGRAREck5JQMiIiI5p2RAREQk55QMiIiI5JySARERkZxrNQ8q+nzY2jCmR7bgH/rmfbf5H9wx\np9dKrvjfcrm7jNDf95CMz4d3cpfx/b6+h/y8xvdc8QO8KwMYbgNc8cdP8j9E6KoeZ7vi1+ryqbuM\nZ4bu4oo/4dDLXPFXM9AVD5854ytr/dvH06F39vgP563jmv+8Cd1c8eEDVzgAF83wPUxnfLeb3WWM\n+O9i3wcO9j9wae4iX/918eJ3XPGXXj3YFQ/Ah+f74vv4i7i7/x6u+J+/f7+7jAN73OaK34ZnXfEr\nVq3P7RnimnRkwMy2N7P7zOxjM1tsZvul3r8xmV78erAp6yAi5aV2L9L6NfVhgk7Aq8AJQF3p50PA\nGsCaycu3GygiLY3avUgr16SHCUIIo4BRAGZW14PIvw4hzGjKckWkctTuRVq/SpxAuIOZTTOzt83s\nGjNbpQJ1EJHyUrsXacHKfQLhQ8BIYBKwIfAn4EEz+0kIwX9Wi4i0Bmr3Ii1cWZOBEMKIoj/fNLM3\ngInADsAT9X74nkGwQtea03r3hz469ChS23+Ae2tMmT9/TkVq0th2P23Q5bTtWvOqnc79d6fzAN8Z\n3iJ5MXbYe4wd/l6Nact92THTZyt6aWEIYZKZfQb0pKFk4KdDYB3HdUYiuXZA8lqiY8c3mT278j+k\nWdv9GkNOpUPvjctXMZFWbosBvdhiQK8a01asWp9jtzyjwc9W9KZDZrY2sCrgv4hbRFoltXuRlqdJ\nRwbMrBMx2y+cUdzDzDYHZiavwcRjh1OTuEuAd4HRTVkPESkftXuR1q+pDxNsSRz2C8mrcNuqm4Hj\ngc2Aw4CuwCfEzuC8EMLCJq6HiJSP2r1IK9fU9xl4kvoPPVT+gKWINCm1e5HWr9U8m4CnDTrXdT+T\nlEW+q5VO7Pt3d3XOPsN3r+6tLnnRXcb3//lew0FF1hg8211Grydec8X3sSpXfI8w0RUP8Med/+iK\nv+OxX7rLGDNuV98H3nQXwSaHveyKv/q105wlXOCMn+mMr6zJ534PVt0i+wd+4Cygi/Oqxsu93zf0\nuNT33JMRAw93l8Hwj33xJ3d3F3HNc6e64m/Z/FBX/OzP13TFAzDLt/6Cv3vk+/aWK75njzfcZdx5\n1WGu+LYHLHDFD5z9aqY4PbVQREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyrnUnA1OHVboG\nFXHHyErXoEKm5XN981BOl7uUSTn+Lj7I57IPG1fpGlTGsFfK+wwvJQOtUH6TgeGVrkFljMrpcpcy\nOcffxQf5XPa8JgPDfVdxL7XWnQyIiIjIUlMyICIiknNKBkRERHKuNdyOuAMA88bXfmfRLJhd4sDK\nR74TLz6umu6v1TTfAZ0Pqz5zFzF2Tunps2bD2BJ3EQ6f+g8yLah61xU/0z5wxa8QprriAZhTx3Is\n+rLke3OrfLdtBmDiqr74Sf4ivqp62/eBd+vIzed+CeNLfSe+JwBXV88q/LeD64PlF+s3q8T3t/BL\n+LyO7cN5V15mNRxSQ/A/cXlBVYl+qz4zF9X93sIvYWapZXf2X9On+eIB3vGFVwffra+r6vlqZ31d\nx/vOLjVM8cUDTK762hW/IDjXN1D1Yenfq1lf1fHeuLGu+c/86NuVV2+7txDKe8ail5kNBG6rdD1E\nliEHhxBur3Ql6qI2L9Is6m33rSEZWBXYHZgM+J7QICLFOgDrA6NDCJ9XuC51UpsXaVKZ2n2LTwZE\nRESkeekEQhERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZEBERyblWmQyY2QlmNsnMvjKz583sR5Wu\nU3Mzs8Fmtjj1eqvS9WpqZra9md1nZh8ny7hfiZgLzewTM5tvZo+YWc9K1LUpNbTcZnZjifX/YKXq\nWwl5a/dq8zVilrk2Dy2r3be6ZMDMfglcDgwGtgBeA0ab2WoVrVh5jAPWANZMXttVtjrNohPwKnAC\nUOu6VzM7A/gNcAywFTCPuP6XL2clm0G9y514iJrrf0B5qlZ5OW73avPLbpuHFtTuW8PtiNMGAdeG\nEIYCmNmxwN7Ar4BLK1mxMlgUQphR6Uo0pxDCKGAUgJlZiZCTgYtCCPcnMYcB04ADgBHlqmdTy7Dc\nAF8v6+u/Hnlt92rzy2ibh5bV7lvVyICZtQP6AI8VpoV416RHgZ9Uql5l1CsZTppoZrea2TqVrlA5\nmdkGxMy4eP3PBl4gH+t/BzObZmZvm9k1ZrZKpStUDjlv92rz+W7zUKZ236qSAWA1oC0xKyw2jbjB\nLMueB44g3qb1WGAD4Ckz61TJSpXZmsShtDyu/4eAw4CdgNOBfsCD9exNLEvy2u7V5vPd5qGM7b41\nHiYoxaj7eMsyIYQwuujPcWb2IvAB8AvgxsrUqsXIw/ovHg5908zeACYCOwBPVKRSlbdMr3e1+Xot\n0+u+oJztvrWNDHwGVBNPpii2OrUzx2VaCGEW8C6wTJxVm9FUYieg9R/CJGJ7yMP6V7tHbT41PVfr\nvqA5232rSgZCCAuBV4CdC9OS4ZKdgWcrVa9KMLMVgQ3xPtS+FUsawlRqrv/OwNbkb/2vDaxKDta/\n2n2kNh/ltc1D87b71niY4ArgZjN7BXiReJZxR+CmSlaquZnZZcD9xGHC7sAFwCJgWCXr1dSS46E9\niXsDAD3MbHNgZgjhQ+BK4Bwzm0B8xO1FwEfAvRWobpOpb7mT12BgJLFj7AlcQtxLHF17bsuk3LV7\ntfllu81DC2v3IYRW9wKOJ24UXwHPAVtWuk5lWOZhxAbwFTAFuB3YoNL1aobl7AcsJg4LF7/+XRRz\nPvAJMD9pFD0rXe/mXG7i88hHJR3CAuB94B9At0rXu8zfUa7avdr8st3mG1r2crd7SyokIiIiOdWq\nzhkQERGRpqdkQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oGREREck7JgIiISM4pGRAREck5JQMi\nIiI5p2RAREQk55QMiIiI5Nz/A+/fWb9QS6P+AAAAAElFTkSuQmCC\n", "text/plain": [ - "" + "" ] }, "metadata": {}, diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index d66b36aa1..450cc26fa 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -374,6 +374,8 @@ Core Classes openmc.data.CoherentElastic openmc.data.FissionEnergyRelease openmc.data.DataLibrary + openmc.data.Decay + openmc.data.FissionProductYields Core Functions -------------- @@ -477,6 +479,7 @@ Functions openmc.data.endf.float_endf openmc.data.endf.get_cont_record + openmc.data.endf.get_evaluations openmc.data.endf.get_head_record openmc.data.endf.get_tab1_record openmc.data.endf.get_tab2_record diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py index ee74233c1..9cc23300d 100644 --- a/examples/python/pincell_multigroup/build-xml.py +++ b/examples/python/pincell_multigroup/build-xml.py @@ -1,3 +1,5 @@ +import numpy as np + import openmc import openmc.mgxs @@ -26,7 +28,7 @@ uo2_xsdata.set_total( 0.5644058]) uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02, 3.0020E-02, 1.1126E-01, 2.8278E-01]) -uo2_xsdata.set_scatter_matrix( +scatter_matrix = np.array( [[[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000], [0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000], [0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000], @@ -34,6 +36,8 @@ uo2_xsdata.set_scatter_matrix( [0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000], [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090], [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]]) +scatter_matrix = np.rollaxis(scatter_matrix, 0, 3) +uo2_xsdata.set_scatter_matrix(scatter_matrix) uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03, 1.85648E-02, 1.78084E-02, 8.30348E-02, 2.16004E-01]) @@ -50,7 +54,7 @@ h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435, h2o_xsdata.set_absorption([6.0105E-04, 1.5793E-05, 3.3716E-04, 1.9406E-03, 5.7416E-03, 1.5001E-02, 3.7239E-02]) -h2o_xsdata.set_scatter_matrix( +scatter_matrix = np.array( [[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000], [0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010], [0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034], @@ -58,6 +62,8 @@ h2o_xsdata.set_scatter_matrix( [0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290], [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200], [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]]) +scatter_matrix = np.rollaxis(scatter_matrix, 0, 3) +h2o_xsdata.set_scatter_matrix(scatter_matrix) mg_cross_sections_file = openmc.MGXSLibrary(groups) mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata]) diff --git a/examples/xml/pincell_multigroup/mgxs.h5 b/examples/xml/pincell_multigroup/mgxs.h5 index 1d5561b00..ae3b17e90 100644 Binary files a/examples/xml/pincell_multigroup/mgxs.h5 and b/examples/xml/pincell_multigroup/mgxs.h5 differ diff --git a/openmc/__init__.py b/openmc/__init__.py index b5530ed84..40a4dcae4 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -25,6 +25,7 @@ from openmc.statepoint import * from openmc.summary import * from openmc.particle_restart import * from openmc.mixin import * +from openmc.plotter import * try: from openmc.opencg_compatible import * diff --git a/openmc/clean_xml.py b/openmc/clean_xml.py index 564281a5c..6aaf64c66 100644 --- a/openmc/clean_xml.py +++ b/openmc/clean_xml.py @@ -65,25 +65,25 @@ def sort_xml_elements(tree): tree.extend(sorted_elements) -def clean_xml_indentation(element, level=0): +def clean_xml_indentation(element, level=0, spaces_per_level=4): """ copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint it basically walks your tree and adds spaces and newlines so the tree is printed in a nice way """ - i = "\n" + level*" " + i = "\n" + level*spaces_per_level*" " if len(element): if not element.text or not element.text.strip(): - element.text = i + " " + element.text = i + spaces_per_level*" " if not element.tail or not element.tail.strip(): element.tail = i for sub_element in element: - clean_xml_indentation(sub_element, level+1) + clean_xml_indentation(sub_element, level+1, spaces_per_level) if not sub_element.tail or not sub_element.tail.strip(): sub_element.tail = i diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 4e1fccb8b..c3df3f104 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -368,7 +368,7 @@ class CMFD(object): @cmfd_mesh.setter def cmfd_mesh(self, mesh): check_type('CMFD mesh', mesh, CMFDMesh) - self._mesh = mesh + self._cmfd_mesh = mesh @norm.setter def norm(self, norm): @@ -446,8 +446,8 @@ class CMFD(object): element.text = str(self._ktol) def _create_mesh_subelement(self): - if self._mesh is not None: - xml_element = self._mesh._get_xml_element() + if self._cmfd_mesh is not None: + xml_element = self._cmfd_mesh._get_xml_element() self._cmfd_file.append(xml_element) def _create_norm_subelement(self): diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 373136538..c361a204d 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -6,6 +6,7 @@ HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR) from .data import * from .neutron import * +from .decay import * from .reaction import * from .ace import * from .angle_distribution import * diff --git a/openmc/data/data.py b/openmc/data/data.py index a3e715812..92be7cade 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -186,3 +186,9 @@ K_BOLTZMANN = 8.6173324e-5 # Used for converting units in ACE data EV_PER_MEV = 1.0e6 + +# Avogadro's constant from CODATA 2010 +AVOGADRO = 6.02214129E23 + +# Neutron mass from CODATA 2010 in units of amu +NEUTRON_MASS = 1.008664916 diff --git a/openmc/data/decay.py b/openmc/data/decay.py new file mode 100644 index 000000000..4327b2fc3 --- /dev/null +++ b/openmc/data/decay.py @@ -0,0 +1,485 @@ +from collections import Iterable, namedtuple +from io import StringIO +from math import log +from numbers import Real +import re +from warnings import warn + +from six import string_types +import numpy as np +try: + from uncertainties import ufloat, unumpy, UFloat +except ImportError: + ufloat = UFloat = namedtuple('UFloat', ['nominal_value', 'std_dev']) + +import openmc.checkvalue as cv +from openmc.mixin import EqualityMixin +from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER +from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record + + +# Gives name and (change in A, change in Z) resulting from decay +_DECAY_MODES = { + 0: ('gamma', (0, 0)), + 1: ('beta-', (0, 1)), + 2: ('ec/beta+', (0, -1)), + 3: ('IT', (0, 0)), + 4: ('alpha', (-4, -2)), + 5: ('n', (-1, 0)), + 6: ('sf', None), + 7: ('p', (-1, -1)), + 8: ('e-', (0, 0)), + 9: ('xray', (0, 0)), + 10: ('unknown', None) +} + +_RADIATION_TYPES = { + 0: 'gamma', + 1: 'beta-', + 2: 'ec/beta+', + 4: 'alpha', + 5: 'n', + 6: 'sf', + 7: 'p', + 8: 'e-', + 9: 'xray', + 10: 'anti-neutrino', + 11: 'neutrino' +} + + +def get_decay_modes(value): + """Return sequence of decay modes given an ENDF RTYP value. + + Parameters + ---------- + value : float + ENDF definition of sequence of decay modes + + Returns + ------- + list of str + List of successive decays, e.g. ('beta-', 'neutron') + + """ + return [_DECAY_MODES[int(x)][0] for x in + str(value).strip('0').replace('.', '')] + + +class FissionProductYields(EqualityMixin): + """Independent and cumulative fission product yields. + + Parameters + ---------- + ev_or_filename : str of openmc.data.endf.Evaluation + ENDF fission product yield evaluation to read from. If given as a + string, it is assumed to be the filename for the ENDF file. + + Attributes + ---------- + cumulative : list of dict + Cumulative yields for each tabulated energy. Each item in the list is a + dictionary whose keys are nuclide names and values are cumulative + yields. The i-th dictionary corresponds to the i-th incident neutron + energy. + energies : Iterable of float or None + Energies at which fission product yields are tabulated. + independent : list of dict + Independent yields for each tabulated energy. Each item in the list is a + dictionary whose keys are nuclide names and values are independent + yields. The i-th dictionary corresponds to the i-th incident neutron + energy. + nuclide : dict + Properties of the fissioning nuclide. + + Notes + ----- + Neutron fission yields are typically not measured with a monoenergetic + source of neutrons. As such, if the fission yields are given at, e.g., + 0.0253 eV, one should interpret this as meaning that they are derived from a + typical thermal reactor flux spectrum as opposed to a monoenergetic source + at 0.0253 eV. + + """ + def __init__(self, ev_or_filename): + # Define function that can be used to read both independent and + # cumulative yields + def get_yields(file_obj): + # Determine number of energies + n_energy = get_head_record(file_obj)[2] + energies = np.zeros(n_energy) + + data = [] + for i in range(n_energy): + # Determine i-th energy and number of products + items, values = get_list_record(file_obj) + energies[i] = items[0] + n_products = items[5] + + # Get yields for i-th energy + yields = {} + for j in range(n_products): + Z, A = divmod(int(values[4*j]), 1000) + isomeric_state = int(values[4*j + 1]) + name = ATOMIC_SYMBOL[Z] + str(A) + if isomeric_state > 0: + name += '_m{}'.format(isomeric_state) + yield_j = ufloat(values[4*j + 2], values[4*j + 3]) + yields[name] = yield_j + + data.append(yields) + + return energies, data + + # Get evaluation if str is passed + if isinstance(ev_or_filename, Evaluation): + ev = ev_or_filename + else: + ev = Evaluation(ev_or_filename) + + # Assign basic nuclide properties + self.nuclide = { + 'name': ev.gnd_name, + 'atomic_number': ev.target['atomic_number'], + 'mass_number': ev.target['mass_number'], + 'isomeric_state': ev.target['isomeric_state'] + } + + # Read independent yields + if (8, 454) in ev.section: + file_obj = StringIO(ev.section[8, 454]) + self.energies, self.independent = get_yields(file_obj) + + # Read cumulative yields + if (8, 459) in ev.section: + file_obj = StringIO(ev.section[8, 459]) + energies, self.cumulative = get_yields(file_obj) + assert np.all(energies == self.energies) + + @classmethod + def from_endf(cls, ev_or_filename): + """Generate fission product yield data from an ENDF evaluation + + Parameters + ---------- + ev_or_filename : str or openmc.data.endf.Evaluation + ENDF fission product yield evaluation to read from. If given as a + string, it is assumed to be the filename for the ENDF file. + + Returns + ------- + openmc.data.FissionProductYields + Fission product yield data + + """ + return cls(ev_or_filename) + + +class DecayMode(EqualityMixin): + """Radioactive decay mode. + + Parameters + ---------- + parent : str + Parent decaying nuclide + modes : list of str + Successive decay modes + daughter_state : int + Metastable state of the daughter nuclide + energy : uncertainties.UFloat + Total decay energy in eV available in the decay process. + branching_ratio : uncertainties.UFloat + Fraction of the decay of the parent nuclide which proceeds by this mode. + + Attributes + ---------- + branching_ratio : uncertainties.UFloat + Fraction of the decay of the parent nuclide which proceeds by this mode. + daughter : str + Name of daughter nuclide produced from decay + energy : uncertainties.UFloat + Total decay energy in eV available in the decay process. + modes : list of str + Successive decay modes + parent : str + Parent decaying nuclide + + """ + + def __init__(self, parent, modes, daughter_state, energy, + branching_ratio): + self._daughter_state = daughter_state + self.parent = parent + self.modes = modes + self.energy = energy + self.branching_ratio = branching_ratio + + def __repr__(self): + return (' {}, {}>'.format( + ','.join(self.modes), self.parent, self.daughter, + self.branching_ratio)) + + @property + def branching_ratio(self): + return self._branching_ratio + + @property + def daughter(self): + # Determine atomic number and mass number of parent + symbol, A = re.match(r'([A-Zn][a-z]*)(\d+)', self.parent).groups() + A = int(A) + Z = ATOMIC_NUMBER[symbol] + + # Process changes + for mode in self.modes: + for name, changes in _DECAY_MODES.values(): + if name == mode: + if changes is not None: + delta_A, delta_Z = changes + A += delta_A + Z += delta_Z + + if self._daughter_state > 0: + return '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A, self._daughter_state) + else: + return '{}{}'.format(ATOMIC_SYMBOL[Z], A) + + @property + def energy(self): + return self._energy + + @property + def modes(self): + return self._modes + + @property + def parent(self): + return self._parent + + @branching_ratio.setter + def branching_ratio(self, branching_ratio): + cv.check_type('branching ratio', branching_ratio, UFloat) + cv.check_greater_than('branching ratio', + branching_ratio.nominal_value, 0.0, True) + if branching_ratio.nominal_value == 0.0: + warn('Decay mode {} of parent {} has a zero branching ratio.' + .format(self.modes, self.parent)) + cv.check_greater_than('branching ratio uncertainty', + branching_ratio.std_dev, 0.0, True) + self._branching_ratio = branching_ratio + + @energy.setter + def energy(self, energy): + cv.check_type('decay energy', energy, UFloat) + cv.check_greater_than('decay energy', energy.nominal_value, 0.0, True) + cv.check_greater_than('decay energy uncertainty', + energy.std_dev, 0.0, True) + self._energy = energy + + @modes.setter + def modes(self, modes): + cv.check_type('decay modes', modes, Iterable, string_types) + self._modes = modes + + @parent.setter + def parent(self, parent): + cv.check_type('parent nuclide', parent, string_types) + self._parent = parent + + +class Decay(EqualityMixin): + """Radioactive decay data. + + Parameters + ---------- + ev_or_filename : str of openmc.data.endf.Evaluation + ENDF radioactive decay data evaluation to read from. If given as a + string, it is assumed to be the filename for the ENDF file. + + Attributes + ---------- + average_energies : dict + Average decay energies in eV of each type of radiation for decay heat + applications. + decay_constant : uncertainties.UFloat + Decay constant in inverse seconds. + half_life : uncertainties.UFloat + Half-life of the decay in seconds. + modes : list + Decay mode information for each mode of decay. + nuclide : dict + Dictionary describing decaying nuclide with keys 'name', + 'excited_state', 'mass', 'stable', 'spin', and 'parity'. + spectra : dict + Resulting radiation spectra for each radiation type. + + """ + def __init__(self, ev_or_filename): + # Get evaluation if str is passed + if isinstance(ev_or_filename, Evaluation): + ev = ev_or_filename + else: + ev = Evaluation(ev_or_filename) + + file_obj = StringIO(ev.section[8, 457]) + + self.nuclide = {} + self.modes = [] + self.spectra = {} + self.average_energies = {} + + # Get head record + items = get_head_record(file_obj) + Z, A = divmod(items[0], 1000) + metastable = items[3] + self.nuclide['atomic_number'] = Z + self.nuclide['mass_number'] = A + self.nuclide['isomeric_state'] = metastable + if metastable > 0: + self.nuclide['name'] = '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A, + metastable) + else: + self.nuclide['name'] = '{}{}'.format(ATOMIC_SYMBOL[Z], A) + self.nuclide['mass'] = items[1] # AWR + self.nuclide['excited_state'] = items[2] # State of the original nuclide + self.nuclide['stable'] = (items[4] == 1) # Nucleus stability flag + + # Determine if radioactive/stable + if not self.nuclide['stable']: + NSP = items[5] # Number of radiation types + + # Half-life and decay energies + items, values = get_list_record(file_obj) + self.half_life = ufloat(items[0], items[1]) + NC = items[4]//2 + pairs = [x for x in zip(values[::2], values[1::2])] + ex = self.average_energies + ex['light'] = ufloat(*pairs[0]) + ex['electromagnetic'] = ufloat(*pairs[1]) + ex['heavy'] = ufloat(*pairs[2]) + if NC == 17: + ex['beta-'] = ufloat(*pairs[3]) + ex['beta+'] = ufloat(*pairs[4]) + ex['auger'] = ufloat(*pairs[5]) + ex['conversion'] = ufloat(*pairs[6]) + ex['gamma'] = ufloat(*pairs[7]) + ex['xray'] = ufloat(*pairs[8]) + ex['Bremsstrahlung'] = ufloat(*pairs[9]) + ex['annihilation'] = ufloat(*pairs[10]) + ex['alpha'] = ufloat(*pairs[11]) + ex['recoil'] = ufloat(*pairs[12]) + ex['SF'] = ufloat(*pairs[13]) + ex['neutron'] = ufloat(*pairs[14]) + ex['proton'] = ufloat(*pairs[15]) + ex['neutrino'] = ufloat(*pairs[16]) + + items, values = get_list_record(file_obj) + spin = items[0] + if spin == -77.777: + self.nuclide['spin'] = None + else: + self.nuclide['spin'] = spin + self.nuclide['parity'] = items[1] # Parity of the nuclide + + # Decay mode information + n_modes = items[5] # Number of decay modes + for i in range(n_modes): + decay_type = get_decay_modes(values[6*i]) + isomeric_state = int(values[6*i + 1]) + energy = ufloat(*values[6*i + 2:6*i + 4]) + branching_ratio = ufloat(*values[6*i + 4:6*(i + 1)]) + + mode = DecayMode(self.nuclide['name'], decay_type, isomeric_state, + energy, branching_ratio) + self.modes.append(mode) + + discrete_type = {0.0: None, 1.0: 'allowed', 2.0: 'first-forbidden', + 3.0: 'second-forbidden', 4.0: 'third-forbidden', + 5.0: 'fourth-forbidden', 6.0: 'fifth-forbidden'} + + # Read spectra + for i in range(NSP): + spectrum = {} + + items, values = get_list_record(file_obj) + # Decay radiation type + spectrum['type'] = _RADIATION_TYPES[items[1]] + # Continuous spectrum flag + spectrum['continuous_flag'] = {0: 'discrete', 1: 'continuous', + 2: 'both'}[items[2]] + spectrum['discrete_normalization'] = ufloat(*values[0:2]) + spectrum['energy_average'] = ufloat(*values[2:4]) + spectrum['continuous_normalization'] = ufloat(*values[4:6]) + + NER = items[5] # Number of tabulated discrete energies + + if not spectrum['continuous_flag'] == 'continuous': + # Information about discrete spectrum + spectrum['discrete'] = [] + for j in range(NER): + items, values = get_list_record(file_obj) + di = {} + di['energy'] = ufloat(*items[0:2]) + di['from_mode'] = get_decay_modes(values[0]) + di['type'] = discrete_type[values[1]] + di['intensity'] = ufloat(*values[2:4]) + if spectrum['type'] == 'ec/beta+': + di['positron_intensity'] = ufloat(*values[4:6]) + elif spectrum['type'] == 'gamma': + di['internal_pair'] = ufloat(*values[4:6]) + if len(values) >= 8: + di['total_internal_conversion'] = ufloat(*values[6:8]) + if len(values) == 12: + di['k_shell_conversion'] = ufloat(*values[8:10]) + di['l_shell_conversion'] = ufloat(*values[10:12]) + spectrum['discrete'].append(di) + + if not spectrum['continuous_flag'] == 'discrete': + # Read continuous spectrum + ci = {} + params, ci['probability'] = get_tab1_record(file_obj) + ci['type'] = get_decay_modes(params[0]) + + # Read covariance (Ek, Fk) table + LCOV = params[3] + if LCOV != 0: + items, values = get_list_record(file_obj) + ci['covariance_lb'] = items[3] + ci['covariance'] = zip(values[0::2], values[1::2]) + + spectrum['continuous'] = ci + + # Add spectrum to dictionary + self.spectra[spectrum['type']] = spectrum + + else: + items, values = get_list_record(file_obj) + items, values = get_list_record(file_obj) + self.nuclide['spin'] = items[0] + self.nuclide['parity'] = items[1] + + @property + def decay_constant(self): + if hasattr(self.half_life, 'n'): + return log(2.)/self.half_life + else: + mu, sigma = self.half_life + return ufloat(log(2.)/mu, log(2.)/mu**2*sigma) + + @classmethod + def from_endf(cls, ev_or_filename): + """Generate radioactive decay data from an ENDF evaluation + + Parameters + ---------- + ev_or_filename : str or openmc.data.endf.Evaluation + ENDF radioactive decay data evaluation to read from. If given as a + string, it is assumed to be the filename for the ENDF file. + + Returns + ------- + openmc.data.Decay + Radioactive decay data + + """ + return cls(ev_or_filename) diff --git a/openmc/data/endf.py b/openmc/data/endf.py index c0049f977..34553ad2d 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -14,9 +14,11 @@ import os from math import pi from collections import OrderedDict, Iterable +from six import string_types import numpy as np from numpy.polynomial.polynomial import Polynomial +from .data import ATOMIC_SYMBOL from .function import Tabulated1D, INTERPOLATION_SCHEME from openmc.stats.univariate import Uniform, Tabular, Legendre @@ -47,16 +49,6 @@ SUM_RULES = {1: [2, 3], _ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)') -def radiation_type(value): - p = {0: 'gamma', 1: 'beta-', 2: 'ec/beta+', 3: 'IT', - 4: 'alpha', 5: 'neutron', 6: 'sf', 7: 'proton', - 8: 'e-', 9: 'xray', 10: 'unknown'} - if value % 1.0 == 0: - return p[int(value)] - else: - return (p[int(value)], p[int(10*value % 10)]) - - def float_endf(s): """Convert string of floating point number in ENDF to float. @@ -258,14 +250,40 @@ def get_tab2_record(file_obj): return params, Tabulated2D(breakpoints, interpolation) +def get_evaluations(filename): + """Return a list of all evaluations within an ENDF file. + + Parameters + ---------- + filename : str + Path to ENDF-6 formatted file + + Returns + ------- + list + A list of :class:`openmc.data.endf.Evaluation` instances. + + """ + evaluations = [] + with open(filename, 'r') as fh: + while True: + pos = fh.tell() + line = fh.readline() + if line[66:70] == ' -1': + break + fh.seek(pos) + evaluations.append(Evaluation(fh)) + return evaluations + class Evaluation(object): """ENDF material evaluation with multiple files/sections Parameters ---------- - filename : str - Path to ENDF file to read + filename_or_obj : str or file-like + Path to ENDF file to read or an open file positioned at the start of an + ENDF material Attributes ---------- @@ -282,8 +300,11 @@ class Evaluation(object): indicator (MOD). """ - def __init__(self, filename): - fh = open(filename, 'r') + def __init__(self, filename_or_obj): + if isinstance(filename_or_obj, string_types): + fh = open(filename_or_obj, 'r') + else: + fh = filename_or_obj self.section = {} self.info = {} self.target = {} @@ -313,6 +334,7 @@ class Evaluation(object): # If end of material reached, exit loop if MAT == 0: + fh.readline() break section_data = '' @@ -396,6 +418,16 @@ class Evaluation(object): mod = 0 self.reaction_list.append((mf, mt, nc, mod)) + @property + def gnd_name(self): + symbol = ATOMIC_SYMBOL[self.target['atomic_number']] + A = self.target['mass_number'] + m = self.target['isomeric_state'] + if m > 0: + return '{}{}_m{}'.format(symbol, A, m) + else: + return '{}{}'.format(symbol, A) + class Tabulated2D(object): """Metadata for a two-dimensional function. diff --git a/openmc/data/laboratory.py b/openmc/data/laboratory.py index be449b79a..0a8908362 100644 --- a/openmc/data/laboratory.py +++ b/openmc/data/laboratory.py @@ -137,3 +137,6 @@ class LaboratoryAngleEnergy(AngleEnergy): energy_out.append(energy_out_i) return cls(tab2.breakpoints, tab2.interpolation, energy, mu, energy_out) + + def to_hdf5(self, group): + raise NotImplementedError diff --git a/openmc/data/library.py b/openmc/data/library.py index 9f6d93146..c179f78f8 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -1,10 +1,12 @@ import os import xml.etree.ElementTree as ET +from six import string_types import h5py from openmc.mixin import EqualityMixin from openmc.clean_xml import clean_xml_indentation +from openmc.checkvalue import check_type class DataLibrary(EqualityMixin): @@ -95,13 +97,14 @@ class DataLibrary(EqualityMixin): method='xml') @classmethod - def from_xml(cls, path): + def from_xml(cls, path=None): """Read cross section data library from an XML file. Parameters ---------- - path : str - Path to XML file to read. + path : str, optional + Path to XML file to read. If not provided, the + `OPENMC_CROSS_SECTIONS` environment variable will be used. Returns ------- @@ -112,6 +115,18 @@ class DataLibrary(EqualityMixin): data = cls() + # If path is None, get the cross sections from the + # OPENMC_CROSS_SECTIONS environment variable + if path is None: + path = os.environ.get('OPENMC_CROSS_SECTIONS') + + # Check to make sure there was an environmental variable. + if path is None: + raise ValueError("Either path or OPENMC_CROSS_SECTIONS " + "environmental variable must be set") + + check_type('path', path, string_types) + tree = ET.parse(path) root = tree.getroot() if root.find('directory') is not None: diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index f478aafc7..581adc417 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -511,6 +511,9 @@ class IncidentNeutron(EqualityMixin): rxs = [data[mt] for mt in SUM_RULES[mt_sum] if mt in data] if len(rxs) > 0: data.summed_reactions[mt_sum] = rx = Reaction(mt_sum) + if rx.mt == 18 and 'total_nu' in group: + tgroup = group['total_nu'] + rx.derived_products.append(Product.from_hdf5(tgroup)) for T in data.temperatures: rx.xs[T] = Sum([rx_i.xs[T] for rx_i in rxs]) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 23b864bf3..bb976f9fa 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -55,13 +55,14 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)', 195: '(n,4n2a)', 196: '(n,4npa)', 197: '(n,3p)', 198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 444: '(n,damage)', 649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)', - 849: '(n,ac)'} -REACTION_NAME.update({i: '(n,n{})'.format(i-50) for i in range(50, 91)}) -REACTION_NAME.update({i: '(n,p{})'.format(i-600) for i in range(600, 649)}) -REACTION_NAME.update({i: '(n,d{})'.format(i-650) for i in range(650, 699)}) -REACTION_NAME.update({i: '(n,t{})'.format(i-700) for i in range(700, 749)}) -REACTION_NAME.update({i: '(n,3He{})'.format(i-750) for i in range(750, 799)}) -REACTION_NAME.update({i: '(n,a{})'.format(i-800) for i in range(800, 849)}) + 849: '(n,ac)', 891: '(n,2nc)'} +REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(50, 91)}) +REACTION_NAME.update({i: '(n,p{})'.format(i - 600) for i in range(600, 649)}) +REACTION_NAME.update({i: '(n,d{})'.format(i - 650) for i in range(650, 699)}) +REACTION_NAME.update({i: '(n,t{})'.format(i - 700) for i in range(700, 749)}) +REACTION_NAME.update({i: '(n,3He{})'.format(i - 750) for i in range(750, 799)}) +REACTION_NAME.update({i: '(n,a{})'.format(i - 800) for i in range(800, 849)}) +REACTION_NAME.update({i: '(n,2n{})'.format(i - 875) for i in range(875, 891)}) def _get_products(ev, mt): @@ -817,7 +818,7 @@ class Reaction(EqualityMixin): Parameters ---------- group : h5py.Group - HDF5 group to write to + HDF5 group to read from energy : dict Dictionary whose keys are temperatures (e.g., '300K') and values are arrays of energies at which cross sections are tabulated at. diff --git a/openmc/element.py b/openmc/element.py index 6e65edab1..edfbee136 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -1,13 +1,12 @@ from collections import OrderedDict import re -import sys import os from six import string_types from xml.etree import ElementTree as ET import openmc -from openmc.checkvalue import check_type, check_length +import openmc.checkvalue as cv from openmc.data import NATURAL_ABUNDANCE, atomic_mass @@ -80,8 +79,8 @@ class Element(object): @name.setter def name(self, name): - check_type('element name', name, string_types) - check_length('element name', name, 1, 2) + cv.check_type('element name', name, string_types) + cv.check_length('element name', name, 1, 2) self._name = name @scattering.setter @@ -254,6 +253,6 @@ class Element(object): for nuclide, abundance in abundances.items(): nuc = openmc.Nuclide(nuclide) nuc.scattering = self.scattering - isotopes.append((nuc, percent*abundance, percent_type)) + isotopes.append((nuc, percent * abundance, percent_type)) return isotopes diff --git a/openmc/filter.py b/openmc/filter.py index e16499660..ac6e961cf 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1,8 +1,7 @@ -from abc import ABCMeta, abstractproperty +from abc import ABCMeta from collections import Iterable, OrderedDict import copy from numbers import Real, Integral -import sys from xml.etree import ElementTree as ET from six import add_metaclass @@ -100,9 +99,13 @@ class Filter(object): @classmethod def _recursive_subclasses(cls): """Return all subclasses and their subclasses, etc.""" - subs = cls.__subclasses__() - subsubs = [grand for s in subs for grand in s.__subclasses__()] - return subs + subsubs + all_subclasses = [] + + for subclass in cls.__subclasses__(): + all_subclasses.append(subclass) + all_subclasses.extend(subclass._recursive_subclasses()) + + return all_subclasses @classmethod def from_hdf5(cls, group, **kwargs): @@ -774,7 +777,114 @@ class MeshFilter(Filter): return df -class EnergyFilter(Filter): +class RealFilter(Filter): + """Tally modifier that describes phase-space and other characteristics + + Parameters + ---------- + bins : Iterable of Real + A grid of bin values. + + Attributes + ---------- + bins : Iterable of Real + A grid of bin values. + num_bins : Integral + The number of filter bins + stride : Integral + The number of filter, nuclide and score bins within each of this + filter's bins. + + """ + + def __gt__(self, other): + if type(self) is type(other): + # Compare largest/smallest bin edges in filters + # This logic is used when merging tallies with real filters + return self.bins[0] >= other.bins[-1] + else: + return super(RealFilter, self).__gt__(other) + + @property + def num_bins(self): + return len(self.bins) - 1 + + @num_bins.setter + def num_bins(self, num_bins): + cv.check_type('filter num_bins', num_bins, Integral) + cv.check_greater_than('filter num_bins', num_bins, 0, equality=True) + self._num_bins = num_bins + + def can_merge(self, other): + if type(self) is not type(other): + return False + + if self.bins[0] == other.bins[-1]: + # This low edge coincides with other's high edge + return True + elif self.bins[-1] == other.bins[0]: + # This high edge coincides with other's low edge + return True + else: + return False + + def merge(self, other): + if not self.can_merge(other): + msg = 'Unable to merge "{0}" with "{1}" ' \ + 'filters'.format(type(self), type(other)) + raise ValueError(msg) + + # Merge unique filter bins + merged_bins = np.concatenate((self.bins, other.bins)) + merged_bins = np.unique(merged_bins) + + # Create a new filter with these bins + return type(self)(sorted(merged_bins)) + + def is_subset(self, other): + """Determine if another filter is a subset of this filter. + + If all of the bins in the other filter are included as bins in this + filter, then it is a subset of this filter. + + Parameters + ---------- + other : openmc.Filter + The filter to query as a subset of this filter + + Returns + ------- + bool + Whether or not the other filter is a subset of this filter + + """ + + if type(self) is not type(other): + return False + elif len(self.bins) != len(other.bins): + return False + else: + return np.allclose(self.bins, other.bins) + + def get_bin_index(self, filter_bin): + i = np.where(self.bins == filter_bin[1])[0] + if len(i) == 0: + msg = 'Unable to get the bin index for Filter since "{0}" ' \ + 'is not one of the bins'.format(filter_bin) + raise ValueError(msg) + else: + return i[0] - 1 + + def get_bin(self, bin_index): + cv.check_type('bin_index', bin_index, Integral) + cv.check_greater_than('bin_index', bin_index, 0, equality=True) + cv.check_less_than('bin_index', bin_index, self.num_bins) + + # Construct 2-tuple of lower, upper bins for real-valued filters + return (self.bins[bin_index], self.bins[bin_index + 1]) + + +class EnergyFilter(RealFilter): """Bins tally events based on incident particle energy. Parameters @@ -794,23 +904,16 @@ class EnergyFilter(Filter): """ - def __gt__(self, other): - if type(self) is type(other): - # Compare largest/smallest energy bin edges in energy filters - # This logic is used when merging tallies with energy filters - return self.bins[0] >= other.bins[-1] + def get_bin_index(self, filter_bin): + # Use lower energy bound to find index for RealFilters + deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1] + min_delta = np.min(deltas) + if min_delta < 1E-3: + return deltas.argmin() - 1 else: - return super(EnergyFilter, self).__gt__(other) - - @property - def num_bins(self): - return len(self.bins) - 1 - - @num_bins.setter - def num_bins(self, num_bins): - cv.check_type('filter num_bins', num_bins, Integral) - cv.check_greater_than('filter num_bins', num_bins, 0, equality=True) - self._num_bins = num_bins + msg = 'Unable to get the bin index for Filter since "{0}" ' \ + 'is not one of the bins'.format(filter_bin) + raise ValueError(msg) def check_bins(self, bins): for edge in bins: @@ -829,62 +932,9 @@ class EnergyFilter(Filter): if bins[index] < bins[index-1]: msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ 'since they are not monotonically ' \ - 'increasing'.format(bins, self.type) + 'increasing'.format(bins, type(self)) raise ValueError(msg) - def can_merge(self, other): - if type(self) is not type(other): - return False - - if self.bins[0] == other.bins[-1]: - # This low energy edge coincides with other's high energy edge - return True - elif self.bins[-1] == other.bins[0]: - # This high energy edge coincides with other's low energy edge - return True - else: - return False - - def merge(self, other): - if not self.can_merge(other): - msg = 'Unable to merge "{0}" with "{1}" ' \ - 'filters'.format(self.type, other.type) - raise ValueError(msg) - - # Merge unique filter bins - merged_bins = np.concatenate((self.bins, other.bins)) - merged_bins = np.unique(merged_bins) - - # Create a new filter with these bins - return type(self)(sorted(merged_bins)) - - def is_subset(self, other): - if type(self) is not type(other): - return False - elif len(self.bins) != len(other.bins): - return False - else: - return np.allclose(self.bins, other.bins) - - def get_bin_index(self, filter_bin): - # Use lower energy bound to find index for energy Filters - deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1] - min_delta = np.min(deltas) - if min_delta < 1E-3: - return deltas.argmin() - 1 - else: - msg = 'Unable to get the bin index for Filter since "{0}" ' \ - 'is not one of the bins'.format(filter_bin) - raise ValueError(msg) - - def get_bin(self, bin_index): - cv.check_type('bin_index', bin_index, Integral) - cv.check_greater_than('bin_index', bin_index, 0, equality=True) - cv.check_less_than('bin_index', bin_index, self.num_bins) - - # Construct 2-tuple of lower, upper energies for energy(out) filters - return (self.bins[bin_index], self.bins[bin_index+1]) - def get_pandas_dataframe(self, data_size, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. @@ -1228,7 +1278,7 @@ class DistribcellFilter(Filter): return df -class MuFilter(Filter): +class MuFilter(RealFilter): """Bins tally events based on particle scattering angle. Parameters @@ -1256,8 +1306,82 @@ class MuFilter(Filter): """ + def check_bins(self, bins): + for edge in bins: + if not isinstance(edge, Real): + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is a non-integer or floating point ' \ + 'value'.format(edge, type(self)) + raise ValueError(msg) + elif edge < -1.: + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is less than -1'.format(edge, type(self)) + raise ValueError(msg) + elif edge > 1.: + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is greater than 1'.format(edge, type(self)) + raise ValueError(msg) -class PolarFilter(Filter): + # Check that bin edges are monotonically increasing + for index in range(1, len(bins)): + if bins[index] < bins[index-1]: + msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ + 'since they are not monotonically ' \ + 'increasing'.format(bins, type(self)) + raise ValueError(msg) + + def get_pandas_dataframe(self, data_size, **kwargs): + """Builds a Pandas DataFrame for the Filter's bins. + + This method constructs a Pandas DataFrame object for the filter with + columns annotated by filter bin information. This is a helper method + for :meth:`Tally.get_pandas_dataframe`. + + Parameters + ---------- + data_size : Integral + The total number of bins in the tally corresponding to this filter + + Returns + ------- + pandas.DataFrame + A Pandas DataFrame with one column of the lower energy bound and one + column of upper energy bound for each filter bin. The number of + rows in the DataFrame is the same as the total number of bins in the + corresponding tally, with the filter bin appropriately tiled to map + to the corresponding tally bins. + + Raises + ------ + ImportError + When Pandas is not installed + + See also + -------- + Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() + + """ + + # Initialize Pandas DataFrame + import pandas as pd + df = pd.DataFrame() + + # Extract the lower and upper energy bounds, then repeat and tile + # them as necessary to account for other filters. + lo_bins = np.repeat(self.bins[:-1], self.stride) + hi_bins = np.repeat(self.bins[1:], self.stride) + tile_factor = data_size / len(lo_bins) + lo_bins = np.tile(lo_bins, tile_factor) + hi_bins = np.tile(hi_bins, tile_factor) + + # Add the new energy columns to the DataFrame. + df.loc[:, self.short_name.lower() + ' low'] = lo_bins + df.loc[:, self.short_name.lower() + ' high'] = hi_bins + + return df + + +class PolarFilter(RealFilter): """Bins tally events based on the incident particle's direction. Parameters @@ -1285,6 +1409,30 @@ class PolarFilter(Filter): """ + def check_bins(self, bins): + for edge in bins: + if not isinstance(edge, Real): + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is a non-integer or floating point ' \ + 'value'.format(edge, type(self)) + raise ValueError(msg) + elif edge < 0.: + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is less than 0'.format(edge, type(self)) + raise ValueError(msg) + elif edge > np.pi: + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is greater than pi'.format(edge, type(self)) + raise ValueError(msg) + + # Check that bin edges are monotonically increasing + for index in range(1, len(bins)): + if bins[index] < bins[index-1]: + msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ + 'since they are not monotonically ' \ + 'increasing'.format(bins, type(self)) + raise ValueError(msg) + def get_pandas_dataframe(self, data_size, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. @@ -1330,12 +1478,13 @@ class PolarFilter(Filter): hi_bins = np.tile(hi_bins, tile_factor) # Add the new angle columns to the DataFrame. - df.loc[:, self.type + ' low'] = lo_bins + df.loc[:, 'polar low'] = lo_bins + df.loc[:, 'polar high'] = hi_bins return df -class AzimuthalFilter(Filter): +class AzimuthalFilter(RealFilter): """Bins tally events based on the incident particle's direction. Parameters @@ -1363,6 +1512,30 @@ class AzimuthalFilter(Filter): """ + def check_bins(self, bins): + for edge in bins: + if not isinstance(edge, Real): + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is a non-integer or floating point ' \ + 'value'.format(edge, type(self)) + raise ValueError(msg) + elif edge < -np.pi: + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is less than -pi'.format(edge, type(self)) + raise ValueError(msg) + elif edge > np.pi: + msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ + 'since it is greater than pi'.format(edge, type(self)) + raise ValueError(msg) + + # Check that bin edges are monotonically increasing + for index in range(1, len(bins)): + if bins[index] < bins[index-1]: + msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ + 'since they are not monotonically ' \ + 'increasing'.format(bins, type(self)) + raise ValueError(msg) + def get_pandas_dataframe(self, data_size, distribcell_paths=True): """Builds a Pandas DataFrame for the Filter's bins. @@ -1408,7 +1581,8 @@ class AzimuthalFilter(Filter): hi_bins = np.tile(hi_bins, tile_factor) # Add the new angle columns to the DataFrame. - df.loc[:, self.type + ' low'] = lo_bins + df.loc[:, 'azimuthal low'] = lo_bins + df.loc[:, 'azimuthal high'] = hi_bins return df diff --git a/openmc/macroscopic.py b/openmc/macroscopic.py index 4d3589141..f2521c4ac 100644 --- a/openmc/macroscopic.py +++ b/openmc/macroscopic.py @@ -1,5 +1,3 @@ -import sys - from six import string_types from openmc.checkvalue import check_type @@ -45,7 +43,7 @@ class Macroscopic(object): return hash((self._name)) def __repr__(self): - string = 'Nuclide - {0}\n'.format(self._name) + string = 'Macroscopic - {0}\n'.format(self._name) return string @property diff --git a/openmc/material.py b/openmc/material.py index 70a070ddd..08b001e47 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -3,9 +3,9 @@ from copy import deepcopy from numbers import Real, Integral import warnings from xml.etree import ElementTree as ET -import sys from six import string_types +import numpy as np import openmc import openmc.data @@ -595,7 +595,7 @@ class Material(object): nuclides = OrderedDict() for nuclide, density, density_type in self._nuclides: - nuclides[nuclide.name] = (nuclide, density) + nuclides[nuclide.name] = (nuclide, density, density_type) for ele, ele_pct, ele_pct_type, enr in self._elements: @@ -606,6 +606,80 @@ class Material(object): return nuclides + def get_nuclide_atom_densities(self): + """Returns all nuclides in the material and their atomic densities in + units of atom/b-cm + + Returns + ------- + nuclides : dict + Dictionary whose keys are nuclide names and values are tuples of + (nuclide, density in atom/b-cm) + + """ + + # Expand elements in to nuclides + nuclides = self.get_nuclide_densities() + + sum_density = False + if self.density_units == 'sum': + sum_density = True + density = 0. + elif self.density_units == 'macro': + density = self.density + elif self.density_units == 'g/cc' or self.density_units == 'g/cm3': + density = -self.density + elif self.density_units == 'kg/m3': + density = -0.001 * self.density + elif self.density_units == 'atom/b-cm': + density = self.density + elif self.density_units == 'atom/cm3' or self.density_units == 'atom/cc': + density = 1.E-24 * self.density + + # For ease of processing split out nuc, nuc_density, + # and nuc_density_type in to separate arrays + nucs = [] + nuc_densities = [] + nuc_density_types = [] + for nuclide in nuclides.items(): + nuc, nuc_density, nuc_density_type = nuclide[1] + nucs.append(nuc) + nuc_densities.append(nuc_density) + nuc_density_types.append(nuc_density_type) + + if sum_density: + density = np.sum(nuc_densities) + percent_in_atom = np.all(nuc_density_types == 'ao') + density_in_atom = density > 0. + sum_percent = 0. + + awrs = [] + for n, nuclide in enumerate(nuclides.items()): + awr = openmc.data.atomic_mass(nuclide[0]) + if awr is not None: + awrs.append(awr / openmc.data.NEUTRON_MASS) + else: + raise ValueError(nuclide[0] + " is invalid") + + # Now that we have the awr, lets finish calculating densities + sum_percent = np.sum(nuc_densities) + nuc_densities = nuc_densities / sum_percent + if not density_in_atom: + sum_percent = 0. + for n, nuc in enumerate(nucs): + x = nuc_densities[n] + sum_percent += x * awrs[n] + sum_percent = 1. / sum_percent + density = -density * sum_percent * \ + openmc.data.AVOGADRO / openmc.data.NEUTRON_MASS * 1.E-24 + nuc_densities = density * nuc_densities + + nuclides = OrderedDict() + for n, nuc in enumerate(nucs): + nuclides[nuc] = (nuc, nuc_densities[n]) + + return nuclides + def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") xml_element.set("name", nuclide[0].name) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 08e0bb404..7e0289064 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -59,15 +59,23 @@ class Library(object): The spatial domain(s) for which MGXS in the Library are computed correction : {'P0', None} Apply the P0 correction to scattering matrices if set to 'P0' + scatter_format : {'legendre', 'histogram'} + Representation of the angular scattering distribution (default is + 'legendre') legendre_order : int - The highest legendre moment in the scattering matrices (default is 0) + The highest Legendre moment in the scattering matrix; this is used if + :attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0) + histogram_bins : int + The number of equally-spaced bins for the histogram representation of + the angular scattering distribution; this is used if + :attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16) energy_groups : openmc.mgxs.EnergyGroups Energy group structure for energy condensation num_delayed_groups : int Number of delayed groups estimator : str or None - The tally estimator used to compute multi-group cross sections. If None, - the default for each MGXS type is used. + The tally estimator used to compute multi-group cross sections. + If None, the default for each MGXS type is used. tally_trigger : openmc.Trigger An (optional) tally precision trigger given to each tally used to compute the cross section @@ -101,7 +109,9 @@ class Library(object): self._energy_groups = None self._num_delayed_groups = 0 self._correction = 'P0' + self._scatter_format = 'legendre' self._legendre_order = 0 + self._histogram_bins = 16 self._tally_trigger = None self._all_mgxs = OrderedDict() self._sp_filename = None @@ -130,7 +140,9 @@ class Library(object): clone._domain_type = self.domain_type clone._domains = copy.deepcopy(self.domains) clone._correction = self.correction + clone._scatter_format = self.scatter_format clone._legendre_order = self.legendre_order + clone._histogram_bins = self.histogram_bins clone._energy_groups = copy.deepcopy(self.energy_groups, memo) clone._num_delayed_groups = self.num_delayed_groups clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo) @@ -209,10 +221,18 @@ class Library(object): def correction(self): return self._correction + @property + def scatter_format(self): + return self._scatter_format + @property def legendre_order(self): return self._legendre_order + @property + def histogram_bins(self): + return self._histogram_bins + @property def tally_trigger(self): return self._tally_trigger @@ -267,7 +287,7 @@ class Library(object): def by_nuclide(self, by_nuclide): cv.check_type('by_nuclide', by_nuclide, bool) - if by_nuclide == True and self.domain_type == 'mesh': + if by_nuclide and self.domain_type == 'mesh': raise ValueError('Unable to create MGXS library by nuclide with ' 'mesh domain') @@ -277,7 +297,7 @@ class Library(object): def domain_type(self, domain_type): cv.check_value('domain type', domain_type, openmc.mgxs.DOMAIN_TYPES) - if self.by_nuclide == True and domain_type == 'mesh': + if self.by_nuclide and domain_type == 'mesh': raise ValueError('Unable to create MGXS library by nuclide with ' 'mesh domain') @@ -337,26 +357,72 @@ class Library(object): def correction(self, correction): cv.check_value('correction', correction, ('P0', None)) - if correction == 'P0' and self.legendre_order > 0: - warn('The P0 correction will be ignored since the scattering ' - 'order "{}" is greater than zero'.format(self.legendre_order)) + if self.scatter_format == 'legendre': + if correction == 'P0' and self.legendre_order > 0: + msg = 'The P0 correction will be ignored since the ' \ + 'scattering order {} is greater than '\ + 'zero'.format(self.legendre_order) + warn(msg) + elif self.scatter_format == 'histogram': + msg = 'The P0 correction will be ignored since the ' \ + 'scatter format is set to histogram' + warn(msg) self._correction = correction + @scatter_format.setter + def scatter_format(self, scatter_format): + cv.check_value('scatter_format', scatter_format, + openmc.mgxs.MU_TREATMENTS) + + if scatter_format == 'histogram' and self.correction == 'P0': + msg = 'The P0 correction will be ignored since the ' \ + 'scatter format is set to histogram' + warn(msg) + self.correction = None + + self._scatter_format = scatter_format + @legendre_order.setter def legendre_order(self, legendre_order): cv.check_type('legendre_order', legendre_order, Integral) - cv.check_greater_than('legendre_order', legendre_order, 0, equality=True) + cv.check_greater_than('legendre_order', legendre_order, 0, + equality=True) cv.check_less_than('legendre_order', legendre_order, 10, equality=True) - if self.correction == 'P0' and legendre_order > 0: - msg = 'The P0 correction will be ignored since the scattering ' \ - 'order {} is greater than zero'.format(self.legendre_order) - warn(msg, RuntimeWarning) - self.correction = None + if self.scatter_format == 'legendre': + if self.correction == 'P0' and legendre_order > 0: + msg = 'The P0 correction will be ignored since the ' \ + 'scattering order {} is greater than '\ + 'zero'.format(self.legendre_order) + warn(msg, RuntimeWarning) + self.correction = None + elif self.scatter_format == 'histogram': + msg = 'The legendre order will be ignored since the ' \ + 'scatter format is set to histogram' + warn(msg) self._legendre_order = legendre_order + @histogram_bins.setter + def histogram_bins(self, histogram_bins): + cv.check_type('histogram_bins', histogram_bins, Integral) + cv.check_greater_than('histogram_bins', histogram_bins, 0) + + if self.scatter_format == 'legendre': + msg = 'The histogram bins will be ignored since the ' \ + 'scatter format is set to legendre' + warn(msg) + elif self.scatter_format == 'histogram': + if self.correction == 'P0': + msg = 'The P0 correction will be ignored since ' \ + 'a histogram representation of the scattering '\ + 'kernel is requested' + warn(msg, RuntimeWarning) + self.correction = None + + self._histogram_bins = histogram_bins + @tally_trigger.setter def tally_trigger(self, tally_trigger): cv.check_type('tally trigger', tally_trigger, openmc.Trigger) @@ -420,7 +486,7 @@ class Library(object): mgxs.delayed_groups = None else: delayed_groups \ - = list(range(1,self.num_delayed_groups+1)) + = list(range(1, self.num_delayed_groups + 1)) mgxs.delayed_groups = delayed_groups # If a tally trigger was specified, add it to the MGXS @@ -430,7 +496,9 @@ class Library(object): # Specify whether to use a transport ('P0') correction if isinstance(mgxs, openmc.mgxs.ScatterMatrixXS): mgxs.correction = self.correction + mgxs.scatter_format = self.scatter_format mgxs.legendre_order = self.legendre_order + mgxs.histogram_bins = self.histogram_bins self.all_mgxs[domain.id][mgxs_type] = mgxs @@ -809,13 +877,13 @@ class Library(object): return pickle.load(open(full_filename, 'rb')) def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', - order=None, subdomain=None): + subdomain=None): """Generates an openmc.XSdata object describing a multi-group cross section dataset for writing to an openmc.MGXSLibrary object. Note that this method does not build an XSdata object with nested temperature tables. The temperature of each - XSdata object will be left at the default value of 300K. + XSdata object will be left at the default value of 294K. Parameters ---------- @@ -830,10 +898,6 @@ class Library(object): Provide the macro or micro cross section in units of cm^-1 or barns. Defaults to 'macro'. If the Library object is not tallied by nuclide this will be set to 'macro' regardless. - order : int - Scattering order for this data entry. Default is None, - which will set the XSdata object to use the order of the - Library. subdomain : iterable of int This parameter is not used unless using a mesh domain. In that case, the subdomain is an [i,j,k] index (1-based indexing) of the @@ -863,10 +927,6 @@ class Library(object): cv.check_type('xsdata_name', xsdata_name, string_types) cv.check_type('nuclide', nuclide, string_types) cv.check_value('xs_type', xs_type, ['macro', 'micro']) - cv.check_type('order', order, (type(None), Integral)) - if order is not None: - cv.check_greater_than('order', order, 0, equality=True) - cv.check_less_than('order', order, 10, equality=True) if subdomain is not None: cv.check_iterable_type('subdomain', subdomain, Integral, max_depth=3) @@ -888,16 +948,7 @@ class Library(object): xsdata = openmc.XSdata(name, self.energy_groups) xsdata.num_delayed_groups = self.num_delayed_groups - if order is None: - # Set the order to the Library's order (the defualt behavior) - xsdata.order = self.legendre_order - else: - # Set the order of the xsdata object to the minimum of - # the provided order or the Library's order. - xsdata.order = min(order, self.legendre_order) - - # Right now only 'legendre' data and isotropic weighting is supported - self.scatter_format = 'legendre' + # Right now only isotropic weighting is supported self.representation = 'isotropic' if nuclide != 'total': @@ -1037,10 +1088,17 @@ class Library(object): # accounted for approximately by using an adjusted # absorption cross section. if 'total' in self.mgxs_types or 'transport' in self.mgxs_types: - for i in range(len(xsdata.temperatures)): - xsdata._absorption[i] \ - = np.subtract(xsdata._total[i], np.sum( - xsdata._scatter_matrix[i][0, :, :], axis=1)) + if xsdata.scatter_format == 'legendre': + for i in range(len(xsdata.temperatures)): + xsdata._absorption[i] = \ + np.subtract(xsdata._total[i], np.sum( + xsdata._scatter_matrix[i][0, :, :], axis=1)) + elif xsdata.scatter_format == 'histogram': + for i in range(len(xsdata.temperatures)): + xsdata._absorption[i] = \ + np.subtract(xsdata._total[i], np.sum(np.sum( + xsdata._scatter_matrix[i][:, :, :], axis=0), + axis=1)) return xsdata @@ -1321,7 +1379,7 @@ class Library(object): 'scattering matrix is not provided.') # Total or transport can be present, but if using # self.correction=="P0", then we should use transport. - if (((self.correction is "P0") and + if (((self.correction == "P0") and ('nu-transport' not in self.mgxs_types))): error_flag = True warn('A "nu-transport" MGXS type is required since a "P0" ' diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 1936697ee..3ad67484a 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -59,6 +59,12 @@ _DOMAINS = (openmc.Cell, openmc.Material, openmc.Mesh) +# Supported ScatterMatrixXS and NuScatterMatrixXS angular distribution types +MU_TREATMENTS = ('legendre', 'histogram') + +# Maximum Legendre order supported by OpenMC +_MAX_LEGENDRE = 10 + @add_metaclass(ABCMeta) class MGXS(object): @@ -192,7 +198,9 @@ class MGXS(object): clone._rxn_rate_tally = copy.deepcopy(self._rxn_rate_tally, memo) clone._xs_tally = copy.deepcopy(self._xs_tally, memo) clone._sparse = self.sparse + clone._loaded_sp = self._loaded_sp clone._derived = self.derived + clone._hdf5_key = self._hdf5_key clone._tallies = OrderedDict() for tally_type, tally in self.tallies.items(): @@ -325,7 +333,7 @@ class MGXS(object): @property def num_subdomains(self): - if self.domain_type.startswith('avg('): + if self.domain_type.startswith('sum('): domain_type = self.domain_type[4:-1] else: domain_type = self.domain_type @@ -789,16 +797,22 @@ class MGXS(object): if not isinstance(subdomains, string_types): cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3) + + filters.append(_DOMAIN_TO_FILTER[self.domain_type]) + subdomain_bins = [] for subdomain in subdomains: - filters.append(_DOMAIN_TO_FILTER[self.domain_type]) - filter_bins.append((subdomain,)) + subdomain_bins.append(subdomain) + filter_bins.append(tuple(subdomain_bins)) # Construct list of energy group bounds tuples for all requested groups if not isinstance(groups, string_types): cv.check_iterable_type('groups', groups, Integral) + filters.append(openmc.EnergyFilter) + energy_bins = [] for group in groups: - filters.append(openmc.EnergyFilter) - filter_bins.append((self.energy_groups.get_group_bounds(group),)) + energy_bins.append( + (self.energy_groups.get_group_bounds(group),)) + filter_bins.append(tuple(energy_bins)) # Construct a collection of the nuclides to retrieve from the xs tally if self.by_nuclide: @@ -958,30 +972,27 @@ class MGXS(object): # Construct a collection of the subdomain filter bins to average across if not isinstance(subdomains, string_types): cv.check_iterable_type('subdomains', subdomains, Integral) + subdomains = [(subdomain,) for subdomain in subdomains] + subdomains = [tuple(subdomains)] elif self.domain_type == 'distribcell': - subdomains = np.arange(self.num_subdomains) + subdomains = [i for i in range(self.num_subdomains)] + subdomains = [tuple(subdomains)] else: subdomains = None # Clone this MGXS to initialize the subdomain-averaged version avg_xs = copy.deepcopy(self) + avg_xs._rxn_rate_tally = None + avg_xs._xs_tally = None - if self.derived: - avg_xs._rxn_rate_tally = avg_xs.rxn_rate_tally.average( - filter_type=_DOMAIN_TO_FILTER[self.domain_type], - filter_bins=subdomains) - else: - avg_xs._rxn_rate_tally = None - avg_xs._xs_tally = None + # Average each of the tallies across subdomains + for tally_type, tally in avg_xs.tallies.items(): + filt_type = _DOMAIN_TO_FILTER[self.domain_type] + tally_avg = tally.summation(filter_type=filt_type, + filter_bins=subdomains) + avg_xs.tallies[tally_type] = tally_avg - # Average each of the tallies across subdomains - for tally_type, tally in avg_xs.tallies.items(): - filt_type = _DOMAIN_TO_FILTER[self.domain_type] - tally_avg = tally.average(filter_type=filt_type, - filter_bins=subdomains) - avg_xs.tallies[tally_type] = tally_avg - - avg_xs._domain_type = 'avg({0})'.format(self.domain_type) + avg_xs._domain_type = 'sum({0})'.format(self.domain_type) avg_xs.sparse = self.sparse return avg_xs @@ -1304,8 +1315,8 @@ class MGXS(object): cv.check_iterable_type('subdomains', subdomains, Integral) elif self.domain_type == 'distribcell': subdomains = np.arange(self.num_subdomains, dtype=np.int) - elif self.domain_type == 'avg(distribcell)': - domain_filter = self.xs_tally.find_filter('avg(distribcell)') + elif self.domain_type == 'sum(distribcell)': + domain_filter = self.xs_tally.find_filter('sum(distribcell)') subdomains = domain_filter.bins elif self.domain_type == 'mesh': xyz = [range(1, x+1) for x in self.domain.dimension] @@ -1526,9 +1537,19 @@ class MGXS(object): else: df = df.drop('score', axis=1) + # Determine if change-in-angle bins are included in the MGXS to + # properly tile the group boundaries + if 'mu low' in df: + # Find the length of the mu filters indirectly from the number + # of times the mu bins repeats. + num_mu = int(df.shape[0] / + df[df['mu low'] == df['mu low'][0]].shape[0]) + else: + num_mu = 1 + # Override energy groups bounds with indices all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int) - all_groups = np.repeat(all_groups, len(query_nuclides)) + all_groups = np.repeat(all_groups, len(query_nuclides) * num_mu) if 'energy low [eV]' in df and 'energyout low [eV]' in df: df.rename(columns={'energy low [eV]': 'group in'}, inplace=True) @@ -1784,17 +1805,19 @@ class MatrixMGXS(MGXS): if not isinstance(subdomains, string_types): cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3) + filters.append(_DOMAIN_TO_FILTER[self.domain_type]) + subdomain_bins = [] for subdomain in subdomains: - filters.append(_DOMAIN_TO_FILTER[self.domain_type]) - filter_bins.append((subdomain,)) + subdomain_bins.append(subdomain) + filter_bins.append(tuple(subdomain_bins)) # Construct list of energy group bounds tuples for all requested groups if not isinstance(in_groups, string_types): cv.check_iterable_type('groups', in_groups, Integral) + filters.append(openmc.EnergyFilter) for group in in_groups: - filters.append(openmc.EnergyFilter) - filter_bins.append(( - self.energy_groups.get_group_bounds(group),)) + energy_bins.append((self.energy_groups.get_group_bounds(group),)) + filter_bins.append(tuple(energy_bins)) # Construct list of energy group bounds tuples for all requested groups if not isinstance(out_groups, string_types): @@ -3207,8 +3230,8 @@ class NuScatterXS(MGXS): class ScatterMatrixXS(MatrixMGXS): - r"""A scattering matrix multi-group cross section for one or more Legendre - moments. + r"""A scattering matrix multi-group cross section with the cosine of the + change-in-angle represented as one or more Legendre moments or a histogram. This class can be used for both OpenMC input generation and tally data post-processing to compute spatially-homogenized and energy-integrated @@ -3226,7 +3249,7 @@ class ScatterMatrixXS(MatrixMGXS): For a spatial domain :math:`V`, incoming energy group :math:`[E_{g'},E_{g'-1}]`, and outgoing energy group :math:`[E_g,E_{g-1}]`, - the scattering moments are calculated as: + the Legendre scattering moments are calculated as: .. math:: @@ -3266,9 +3289,18 @@ class ScatterMatrixXS(MatrixMGXS): Attributes ---------- correction : 'P0' or None - Apply the P0 correction to scattering matrices if set to 'P0' + Apply the P0 correction to scattering matrices if set to 'P0'; this is + used only if :attr:`ScatterMatrixXS.scatter_format` is 'legendre' + scatter_format : {'legendre', or 'histogram'} + Representation of the angular scattering distribution (default is + 'legendre') legendre_order : int - The highest Legendre moment in the scattering matrix (default is 0) + The highest Legendre moment in the scattering matrix; this is used if + :attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0) + histogram_bins : int + The number of equally-spaced bins for the histogram representation of + the angular scattering distribution; this is used if + :attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16) name : str, optional Name of the multi-group cross section rxn_type : str @@ -3335,7 +3367,9 @@ class ScatterMatrixXS(MatrixMGXS): groups, by_nuclide, name) self._rxn_type = 'scatter' self._correction = 'P0' + self._scatter_format = 'legendre' self._legendre_order = 0 + self._histogram_bins = 16 self._hdf5_key = 'scatter matrix' self._estimator = 'analog' self._valid_estimators = ['analog'] @@ -3343,26 +3377,39 @@ class ScatterMatrixXS(MatrixMGXS): def __deepcopy__(self, memo): clone = super(ScatterMatrixXS, self).__deepcopy__(memo) clone._correction = self.correction + clone._scatter_format = self.scatter_format clone._legendre_order = self.legendre_order + clone._histogram_bins = self.histogram_bins return clone @property def correction(self): return self._correction + @property + def scatter_format(self): + return self._scatter_format + @property def legendre_order(self): return self._legendre_order + @property + def histogram_bins(self): + return self._histogram_bins + @property def scores(self): scores = ['flux'] - if self.correction == 'P0' and self.legendre_order == 0: - scores += ['{}-0'.format(self.rxn_type), - '{}-1'.format(self.rxn_type)] - else: - scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)] + if self.scatter_format == 'legendre': + if self.correction == 'P0' and self.legendre_order == 0: + scores += ['{}-0'.format(self.rxn_type), + '{}-1'.format(self.rxn_type)] + else: + scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)] + elif self.scatter_format == 'histogram': + scores += [self.rxn_type] return scores @@ -3372,10 +3419,15 @@ class ScatterMatrixXS(MatrixMGXS): energy = openmc.EnergyFilter(group_edges) energyout = openmc.EnergyoutFilter(group_edges) - if self.correction == 'P0' and self.legendre_order == 0: - filters = [[energy], [energy, energyout], [energyout]] - else: - filters = [[energy], [energy, energyout]] + if self.scatter_format == 'legendre': + if self.correction == 'P0' and self.legendre_order == 0: + filters = [[energy], [energy, energyout], [energyout]] + else: + filters = [[energy], [energy, energyout]] + elif self.scatter_format == 'histogram': + bins = np.linspace(-1., 1., num=self.histogram_bins + 1, + endpoint=True) + filters = [[energy], [energy, energyout, openmc.MuFilter(bins)]] return filters @@ -3383,20 +3435,24 @@ class ScatterMatrixXS(MatrixMGXS): def rxn_rate_tally(self): if self._rxn_rate_tally is None: + if self.scatter_format == 'legendre': + # If using P0 correction subtract scatter-1 from the diagonal + if self.correction == 'P0' and self.legendre_order == 0: + scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)] + scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)] + energy_filter = scatter_p0.find_filter(openmc.EnergyFilter) + energy_filter = copy.deepcopy(energy_filter) + scatter_p1 = scatter_p1.diagonalize_filter(energy_filter) + self._rxn_rate_tally = scatter_p0 - scatter_p1 - # If using P0 correction subtract scatter-1 from the diagonal - if self.correction == 'P0' and self.legendre_order == 0: - scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)] - scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)] - energy_filter = scatter_p0.find_filter(openmc.EnergyFilter) - energy_filter = copy.deepcopy(energy_filter) - scatter_p1 = scatter_p1.diagonalize_filter(energy_filter) - self._rxn_rate_tally = scatter_p0 - scatter_p1 - - # Extract scattering moment reaction rate Tally - else: - tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order) - self._rxn_rate_tally = self.tallies[tally_key] + # Extract scattering moment reaction rate Tally + else: + tally_key = '{}-P{}'.format(self.rxn_type, + self.legendre_order) + self._rxn_rate_tally = self.tallies[tally_key] + elif self.scatter_format == 'histogram': + # Extract scattering rate distribution tally + self._rxn_rate_tally = self.tallies[self.rxn_type] self._rxn_rate_tally.sparse = self.sparse @@ -3406,27 +3462,53 @@ class ScatterMatrixXS(MatrixMGXS): def correction(self, correction): cv.check_value('correction', correction, ('P0', None)) - if correction == 'P0' and self.legendre_order > 0: - msg = 'The P0 correction will be ignored since the scattering ' \ - 'order {} is greater than zero'.format(self.legendre_order) + if self.scatter_format == 'legendre': + if correction == 'P0' and self.legendre_order > 0: + msg = 'The P0 correction will be ignored since the ' \ + 'scattering order {} is greater than '\ + 'zero'.format(self.legendre_order) + warnings.warn(msg) + elif self.scatter_format == 'histogram': + msg = 'The P0 correction will be ignored since the ' \ + 'scatter format is set to histogram' warnings.warn(msg) self._correction = correction + @scatter_format.setter + def scatter_format(self, scatter_format): + cv.check_value('scatter_format', scatter_format, MU_TREATMENTS) + self._scatter_format = scatter_format + @legendre_order.setter def legendre_order(self, legendre_order): cv.check_type('legendre_order', legendre_order, Integral) - cv.check_greater_than('legendre_order', legendre_order, 0, equality=True) - cv.check_less_than('legendre_order', legendre_order, 10, equality=True) + cv.check_greater_than('legendre_order', legendre_order, 0, + equality=True) + cv.check_less_than('legendre_order', legendre_order, _MAX_LEGENDRE, + equality=True) - if self.correction == 'P0' and legendre_order > 0: - msg = 'The P0 correction will be ignored since the scattering ' \ - 'order {} is greater than zero'.format(self.legendre_order) - warnings.warn(msg, RuntimeWarning) - self.correction = None + if self.scatter_format == 'legendre': + if self.correction == 'P0' and legendre_order > 0: + msg = 'The P0 correction will be ignored since the ' \ + 'scattering order {} is greater than '\ + 'zero'.format(self.legendre_order) + warnings.warn(msg, RuntimeWarning) + self.correction = None + elif self.scatter_format == 'histogram': + msg = 'The legendre order will be ignored since the ' \ + 'scatter format is set to histogram' + warnings.warn(msg) self._legendre_order = legendre_order + @histogram_bins.setter + def histogram_bins(self, histogram_bins): + cv.check_type('histogram_bins', histogram_bins, Integral) + cv.check_greater_than('histogram_bins', histogram_bins, 0) + + self._histogram_bins = histogram_bins + def load_from_statepoint(self, statepoint): """Extracts tallies in an OpenMC StatePoint with the data needed to compute multi-group cross sections. @@ -3456,12 +3538,16 @@ class ScatterMatrixXS(MatrixMGXS): self._rxn_rate_tally = None self._loaded_sp = False - # Expand scores to match the format in the statepoint - # e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2" - if self.correction != 'P0' or self.legendre_order != 0: - tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order) - self.tallies[tally_key].scores = \ - [self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)] + if self.scatter_format == 'legendre': + # Expand scores to match the format in the statepoint + # e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2" + if self.correction != 'P0' or self.legendre_order != 0: + tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order) + self.tallies[tally_key].scores = \ + [self.rxn_type + '-{}'.format(i) + for i in range(self.legendre_order + 1)] + elif self.scatter_format == 'histogram': + self.tallies[self.rxn_type].scores = [self.rxn_type] super(ScatterMatrixXS, self).load_from_statepoint(statepoint) @@ -3508,7 +3594,7 @@ class ScatterMatrixXS(MatrixMGXS): slice_xs._xs_tally = None # Slice the Legendre order if needed - if legendre_order != 'same': + if legendre_order != 'same' and self.scatter_format == 'legendre': cv.check_type('legendre_order', legendre_order, Integral) cv.check_less_than('legendre_order', legendre_order, self.legendre_order, equality=True) @@ -3517,7 +3603,8 @@ class ScatterMatrixXS(MatrixMGXS): # Slice the scattering tally tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order) expand_scores = \ - [self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)] + [self.rxn_type + '-{}'.format(i) + for i in range(self.legendre_order + 1)] slice_xs.tallies[tally_key] = \ slice_xs.tallies[tally_key].get_slice(scores=expand_scores) @@ -3548,7 +3635,8 @@ class ScatterMatrixXS(MatrixMGXS): This method constructs a 5D NumPy array for the requested multi-group cross section data for one or more subdomains (1st dimension), energy groups in (2nd dimension), energy groups out - (3rd dimension), nuclides (4th dimension), and moments (5th dimension). + (3rd dimension), nuclides (4th dimension), and moments/histograms + (5th dimension). NOTE: The scattering moments are not multiplied by the :math:`(2l+1)/2` prefactor in the expansion of the scattering source into Legendre @@ -3618,16 +3706,21 @@ class ScatterMatrixXS(MatrixMGXS): # Construct a collection of the domain filter bins if not isinstance(subdomains, string_types): cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3) + filters.append(_DOMAIN_TO_FILTER[self.domain_type]) + subdomain_bins = [] for subdomain in subdomains: - filters.append(_DOMAIN_TO_FILTER[self.domain_type]) - filter_bins.append((subdomain,)) + subdomain_bins.append(subdomain) + filter_bins.append(tuple(subdomain_bins)) # Construct list of energy group bounds tuples for all requested groups if not isinstance(in_groups, string_types): cv.check_iterable_type('groups', in_groups, Integral) + filters.append(openmc.EnergyFilter) + energy_bins = [] for group in in_groups: - filters.append(openmc.EnergyFilter) - filter_bins.append((self.energy_groups.get_group_bounds(group),)) + energy_bins.append( + (self.energy_groups.get_group_bounds(group),)) + filter_bins.append(tuple(energy_bins)) # Construct list of energy group bounds tuples for all requested groups if not isinstance(out_groups, string_types): @@ -3637,7 +3730,7 @@ class ScatterMatrixXS(MatrixMGXS): filter_bins.append((self.energy_groups.get_group_bounds(group),)) # Construct CrossScore for requested scattering moment - if moment != 'all': + if moment != 'all' and self.scatter_format == 'legendre': cv.check_type('moment', moment, Integral) cv.check_greater_than('moment', moment, 0, equality=True) cv.check_less_than( @@ -3687,9 +3780,19 @@ class ScatterMatrixXS(MatrixMGXS): else: num_out_groups = len(out_groups) + if self.scatter_format == 'histogram': + num_mu_bins = self.histogram_bins + else: + num_mu_bins = 1 + # Reshape tally data array with separate axes for domain and energy - num_subdomains = int(xs.shape[0] / (num_in_groups * num_out_groups)) - new_shape = (num_subdomains, num_in_groups, num_out_groups) + num_subdomains = int(xs.shape[0] / + (num_mu_bins * num_in_groups * num_out_groups)) + if self.scatter_format == 'histogram': + new_shape = (num_subdomains, num_in_groups, num_out_groups, + num_mu_bins) + else: + new_shape = (num_subdomains, num_in_groups, num_out_groups) new_shape += xs.shape[1:] xs = np.reshape(xs, new_shape) @@ -3700,11 +3803,22 @@ class ScatterMatrixXS(MatrixMGXS): # Reverse data if user requested increasing energy groups since # tally data is stored in order of increasing energies if order_groups == 'increasing': - xs = xs[:, ::-1, ::-1, :] + xs = xs[:, ::-1, ::-1, ...] if squeeze: - xs = np.squeeze(xs) - xs = np.atleast_2d(xs) + # We want to squeeze out everything but the in_groups, out_groups, + # and, if needed, num_mu_bins dimension. These must not be squeezed + # so 1-group problems have the correct shape. + if self.scatter_format == 'histogram': + axes = (5, 4, 0) + else: + axes = (4, 3, 0) + # Squeeze will return a ValueError if the axis has a size greater + # than 1, so try each axis in axes one at a time, catching the + # ValueError as needed. + for axis in axes: + if xs.shape[axis] == 1: + xs = np.squeeze(xs, axis=axis) return xs @@ -3755,28 +3869,41 @@ class ScatterMatrixXS(MatrixMGXS): df = super(ScatterMatrixXS, self).get_pandas_dataframe( groups, nuclides, xs_type, distribcell_paths) - # Add a moment column to dataframe - if self.legendre_order > 0: - # Insert a column corresponding to the Legendre moments - moments = ['P{}'.format(i) for i in range(self.legendre_order+1)] - moments = np.tile(moments, int(df.shape[0] / len(moments))) - df['moment'] = moments + if self.scatter_format == 'legendre': + # Add a moment column to dataframe + if self.legendre_order > 0: + # Insert a column corresponding to the Legendre moments + moments = ['P{}'.format(i) + for i in range(self.legendre_order + 1)] + moments = np.tile(moments, int(df.shape[0] / len(moments))) + df['moment'] = moments - # Place the moment column before the mean column - columns = df.columns.tolist() - mean_index = [i for i, s in enumerate(columns) if 'mean' in s][0] - if self.domain_type == 'mesh': - df = df[columns[:mean_index] + [('moment', '')] + columns[mean_index:-1]] - else: - df = df[columns[:mean_index] + ['moment'] + columns[mean_index:-1]] + # Place the moment column before the mean column + columns = df.columns.tolist() + mean_index \ + = [i for i, s in enumerate(columns) if 'mean' in s][0] + if self.domain_type == 'mesh': + df = df[columns[:mean_index] + [('moment', '')] + + columns[mean_index:-1]] + else: + df = df[columns[:mean_index] + ['moment'] + + columns[mean_index:-1]] - # Select rows corresponding to requested scattering moment - if moment != 'all': - cv.check_type('moment', moment, Integral) - cv.check_greater_than('moment', moment, 0, equality=True) - cv.check_less_than( - 'moment', moment, self.legendre_order, equality=True) - df = df[df['moment'] == 'P{}'.format(moment)] + # Select rows corresponding to requested scattering moment + if moment != 'all': + cv.check_type('moment', moment, Integral) + cv.check_greater_than('moment', moment, 0, equality=True) + cv.check_less_than( + 'moment', moment, self.legendre_order, equality=True) + df = df[df['moment'] == 'P{}'.format(moment)] + + elif self.scatter_format == 'histogram': + # Replace the mu low and mu high columns with a single mu bin + del df['mu high'] + df.rename(columns={'mu low': 'mu bins'}, inplace=True) + bins = [i + 1 for i in range(self.histogram_bins)] + bins = np.tile(bins, int(df.shape[0] / len(bins))) + df['mu bins'] = bins return df @@ -3827,7 +3954,7 @@ class ScatterMatrixXS(MatrixMGXS): cv.check_value('xs_type', xs_type, ['macro', 'micro']) - if self.correction != 'P0': + if self.correction != 'P0' and self.scatter_format == 'legendre': rxn_type = '{0} (P{1})'.format(self.rxn_type, moment) else: rxn_type = self.rxn_type @@ -4606,16 +4733,21 @@ class Chi(MGXS): # Construct a collection of the domain filter bins if not isinstance(subdomains, string_types): cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3) + filters.append(_DOMAIN_TO_FILTER[self.domain_type]) + subdomain_bins = [] for subdomain in subdomains: - filters.append(_DOMAIN_TO_FILTER[self.domain_type]) - filter_bins.append((subdomain,)) + subdomain_bins.append(subdomain) + filter_bins.append(tuple(subdomain_bins)) # Construct list of energy group bounds tuples for all requested groups if not isinstance(groups, string_types): cv.check_iterable_type('groups', groups, Integral) + filters.append(openmc.EnergyoutFilter) + energy_bins = [] for group in groups: - filters.append(openmc.EnergyoutFilter) - filter_bins.append((self.energy_groups.get_group_bounds(group),)) + energy_bins.append( + (self.energy_groups.get_group_bounds(group),)) + filter_bins.append(tuple(energy_bins)) # If chi was computed for each nuclide in the domain if self.by_nuclide: diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 3f3a38b06..10fda4593 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1,6 +1,6 @@ from collections import Iterable from numbers import Real, Integral -import sys +import os from six import string_types import numpy as np @@ -15,7 +15,7 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \ # Supported incoming particle MGXS angular treatment representations _REPRESENTATIONS = ['isotropic', 'angle'] _SCATTER_TYPES = ['tabular', 'legendre', 'histogram'] -_XS_SHAPES = ["[Order][G][G']", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]", +_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]", "[G'][DG]", "[G][G'][DG]"] @@ -42,7 +42,7 @@ class XSdata(object): ---------- name : str Unique identifier for the xsdata object - aromic_weight_ratio : float + atomic_weight_ratio : float Atomic weight ratio of an isotope. That is, the ratio of the mass of the isotope to the mass of a single neutron. temperatures : numpy.ndarray @@ -71,50 +71,50 @@ class XSdata(object): Number of equal width angular bins that the polar angular domain is subdivided into. This only applies when :attr:`XSdata.representation` is "angle". - total : dict of numpy.ndarray + total : list of numpy.ndarray Group-wise total cross section. - absorption : dict of numpy.ndarray + absorption : list of numpy.ndarray Group-wise absorption cross section. - scatter_matrix : dict of numpy.ndarray + scatter_matrix : list of numpy.ndarray Scattering moment matrices presented with the columns representing incoming group and rows representing the outgoing group. That is, down-scatter will be above the diagonal of the resultant matrix. - multiplicity_matrix : dict of numpy.ndarray + multiplicity_matrix : list of numpy.ndarray Ratio of neutrons produced in scattering collisions to the neutrons which undergo scattering collisions; that is, the multiplicity provides the code with a scaling factor to account for neutrons produced in (n,xn) reactions. - fission : dict of numpy.ndarray + fission : list of numpy.ndarray Group-wise fission cross section. - kappa_fission : dict of numpy.ndarray + kappa_fission : list of numpy.ndarray Group-wise kappa_fission cross section. - chi : dict of numpy.ndarray + chi : list of numpy.ndarray Group-wise fission spectra ordered by increasing group index (i.e., fast to thermal). This attribute should be used if making the common approximation that the fission spectra does not depend on incoming energy. If the user does not wish to make this approximation, then this should not be provided and this information included in the :attr:`XSdata.nu_fission` attribute instead. - chi_prompt : dict of numpy.ndarray + chi_prompt : list of numpy.ndarray Group-wise prompt fission spectra ordered by increasing group index (i.e., fast to thermal). This attribute should be used if chi from prompt and delayed neutrons is being set separately. - chi_delayed : dict of numpy.ndarray + chi_delayed : list of numpy.ndarray Group-wise delayed fission spectra ordered by increasing group index (i.e., fast to thermal). This attribute should be used if chi from prompt and delayed neutrons is being set separately. - nu_fission : dict of numpy.ndarray + nu_fission : list of numpy.ndarray Group-wise fission production cross section vector (i.e., if ``chi`` is provided), or is the group-wise fission production matrix. - prompt_nu_fission : dict of numpy.ndarray + prompt_nu_fission : list of numpy.ndarray Group-wise prompt fission production cross section vector. - delayed_nu_fission : dict of numpy.ndarray + delayed_nu_fission : list of numpy.ndarray Group-wise delayed fission production cross section vector. - beta : dict of numpy.ndarray + beta : list of numpy.ndarray Delayed-group-wise delayed neutron fraction cross section vector. - decay_rate : dict of numpy.ndarray + decay_rate : list of numpy.ndarray Delayed-group-wise decay rate vector. - inverse_velocity : dict of numpy.ndarray + inverse_velocity : list of numpy.ndarray Inverse of velocity, in units of sec/cm. xs_shapes : dict of iterable of int Dictionary with keys of _XS_SHAPES and iterable of int values with the @@ -134,7 +134,7 @@ class XSdata(object): Note that some cross sections can be input in more than one shape so they are listed multiple times: - [Order][G][G']: scatter_matrix + [G][G'][Order]: scatter_matrix [G]: total, absorption, fission, kappa_fission, nu_fission, prompt_nu_fission, delayed_nu_fission, inverse_velocity @@ -304,13 +304,13 @@ class XSdata(object): self.energy_groups.num_groups, self.num_delayed_groups) - self._xs_shapes["[Order][G][G']"] \ - = (self.num_orders, self.energy_groups.num_groups, - self.energy_groups.num_groups) + self._xs_shapes["[G][G'][Order]"] \ + = (self.energy_groups.num_groups, + self.energy_groups.num_groups, self.num_orders) # If representation is by angle prepend num polar and num azim if self.representation == 'angle': - for key,shapes in self._xs_shapes.items(): + for key, shapes in self._xs_shapes.items(): self._xs_shapes[key] \ = (self.num_polar, self.num_azimuthal) + shapes @@ -338,7 +338,7 @@ class XSdata(object): def num_delayed_groups(self, num_delayed_groups): # Check validity of num_delayed_groups - check_type('num_delayed_groups', num_delayed_groups, int) + check_type('num_delayed_groups', num_delayed_groups, Integral) check_less_than('num_delayed_groups', num_delayed_groups, openmc.mgxs.MAX_DELAYED_GROUPS, equality=True) check_greater_than('num_delayed_groups', num_delayed_groups, 0, @@ -724,7 +724,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[Order][G][G']"]] + shapes = [self.xs_shapes["[G][G'][Order]"]] # Convert to a numpy array so we can easily get the shape for checking scatter = np.asarray(scatter) @@ -1521,33 +1521,41 @@ class XSdata(object): check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) - if self.scatter_format != 'legendre': - msg = 'Anisotropic scattering representations other than ' \ - 'Legendre expansions have not yet been implemented in ' \ - 'openmc.mgxs.' - raise ValueError(msg) + # Set the value of scatter_format based on the same value within + # scatter + self.scatter_format = scatter.scatter_format # If the user has not defined XSdata.order, then we will set # the order based on the data within scatter. - # Otherwise, we will check to see that XSdata.order to match + # Otherwise, we will check to see that XSdata.order matches # the order of scatter - if self.order is None: - self.order = scatter.legendre_order - else: - check_value('legendre_order', scatter.legendre_order, - [self.order]) + if self.scatter_format == 'legendre': + if self.order is None: + self.order = scatter.legendre_order + else: + check_value('legendre_order', scatter.legendre_order, + [self.order]) + elif self.scatter_format == 'histogram': + if self.order is None: + self.order = scatter.histogram_bins + else: + check_value('histogram_bins', scatter.histogram_bins, + [self.order]) i = np.where(self.temperatures == temperature)[0][0] if self.representation == 'isotropic': - # Get the scattering orders in the outermost dimension - self._scatter_matrix[i] = np.zeros((self.num_orders, - self.energy_groups.num_groups, - self.energy_groups.num_groups)) - for moment in range(self.num_orders): - self._scatter_matrix[i][moment, :, :] = \ + if self.scatter_format == 'legendre': + # Get the scattering orders in the outermost dimension + self._scatter_matrix[i] = \ + np.zeros(self.xs_shapes["[G][G'][Order]"]) + for moment in range(self.num_orders): + self._scatter_matrix[i][:, :, moment] = \ + scatter.get_xs(nuclides=nuclide, xs_type=xs_type, + moment=moment, subdomains=subdomain) + else: + self._scatter_matrix[i] = \ scatter.get_xs(nuclides=nuclide, xs_type=xs_type, - moment=moment, subdomains=subdomain) - + subdomains=subdomain) elif self.representation == 'angle': msg = 'Angular-Dependent MGXS have not yet been implemented' raise ValueError(msg) @@ -1625,6 +1633,10 @@ class XSdata(object): scatt = scatter.get_xs(nuclides=nuclide, xs_type=xs_type, moment=0, subdomains=subdomain) + if scatter.scatter_format == 'histogram': + scatt = np.sum(scatt, axis=0) + if nuscatter.scatter_format == 'histogram': + nuscatt = np.sum(nuscatt, axis=0) self._multiplicity_matrix[i] = np.divide(nuscatt, scatt) elif self.representation == 'angle': msg = 'Angular-Dependent MGXS have not yet been implemented' @@ -1641,6 +1653,7 @@ class XSdata(object): HDF5 File (a root Group) to write to """ + grp = file.create_group(self.name) if self.atomic_weight_ratio is not None: grp.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio @@ -1656,7 +1669,7 @@ class XSdata(object): if self.num_polar is not None: grp.attrs['num_polar'] = self.num_polar - grp.attrs['scatter_shape'] = np.string_("[Order][G][G']") + grp.attrs['scatter_shape'] = np.string_("[G][G'][Order]") if self.scatter_format is not None: grp.attrs['scatter_format'] = np.string_(self.scatter_format) if self.order is not None: @@ -1708,11 +1721,12 @@ class XSdata(object): (self._delayed_nu_fission[i] is None or \ self._prompt_nu_fission[i] is None): raise ValueError('nu-fission or prompt-nu-fission and ' - 'delayed-nu-fission data must be provided ' - 'when writing the HDF5 library') + 'delayed-nu-fission data must be ' + 'provided when writing the HDF5 library') if self._nu_fission[i] is not None: - xs_grp.create_dataset("nu-fission", data=self._nu_fission[i]) + xs_grp.create_dataset("nu-fission", + data=self._nu_fission[i]) if self._prompt_nu_fission[i] is not None: xs_grp.create_dataset("prompt-nu-fission", @@ -1726,7 +1740,8 @@ class XSdata(object): xs_grp.create_dataset("beta", data=self._beta[i]) if self._decay_rate[i] is not None: - xs_grp.create_dataset("decay rate", data=self._decay_rate[i]) + xs_grp.create_dataset("decay rate", + data=self._decay_rate[i]) if self._scatter_matrix[i] is None: raise ValueError('Scatter matrix must be provided when ' @@ -1735,97 +1750,88 @@ class XSdata(object): # Get the sparse scattering data to print to the library G = self.energy_groups.num_groups if self.representation == 'isotropic': - - g_out_bounds = np.zeros((G, 2), dtype=np.int) - - for g_in in range(G): - nz = np.nonzero(self._scatter_matrix[i][0, g_in, :]) - g_out_bounds[g_in, 0] = nz[0][0] - g_out_bounds[g_in, 1] = nz[0][-1] - - # Now create the flattened scatter matrix array - matrix = self._scatter_matrix[i] - flat_scatt = [] - for g_in in range(G): - for g_out in range(g_out_bounds[g_in, 0], - g_out_bounds[g_in, 1] + 1): - for l in range(len(matrix[:, g_in, g_out])): - flat_scatt.append(matrix[l, g_in, g_out]) - - # And write it. - scatt_grp = xs_grp.create_group('scatter_data') - scatt_grp.create_dataset("scatter_matrix", - data=np.array(flat_scatt)) - - # Repeat for multiplicity - if self._multiplicity_matrix[i] is not None: - - # Now create the flattened scatter matrix array - matrix = self._multiplicity_matrix[i][:, :] - flat_mult = [] - for g_in in range(G): - for g_out in range(g_out_bounds[g_in, 0], - g_out_bounds[g_in, 1] + 1): - flat_mult.append(matrix[g_in, g_out]) - - scatt_grp.create_dataset("multiplicity matrix", - data=np.array(flat_mult)) - - # And finally, adjust g_out_bounds for 1-based group counting - # and write it. - g_out_bounds[:, :] += 1 - scatt_grp.create_dataset("g_min", data=g_out_bounds[:, 0]) - scatt_grp.create_dataset("g_max", data=g_out_bounds[:, 1]) - + Np = 1 + Na = 1 elif self.representation == 'angle': Np = self.num_polar Na = self.num_azimuthal - g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int) - for p in range(Np): - for a in range(Na): - for g_in in range(G): - matrix = self._scatter_matrix[i][p, a, 0, g_in, :] - nz = np.nonzero(matrix) + g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int) + for p in range(Np): + for a in range(Na): + for g_in in range(G): + if self.scatter_format == 'legendre': + if self.representation == 'isotropic': + matrix = \ + self._scatter_matrix[i][g_in, :, 0] + elif self.representation == 'angle': + matrix = \ + self._scatter_matrix[i][p, a, g_in, :, 0] + elif self.scatter_format == 'histogram': + if self.representation == 'isotropic': + matrix = \ + np.sum(self._scatter_matrix[i][g_in, :, :], + axis=1) + elif self.representation == 'angle': + matrix = \ + np.sum(self._scatter_matrix[i][p, a, g_in, :, :], + axis=1) + nz = np.nonzero(matrix) + # It is possible that there only zeros in matrix + # and therefore nz will be empty, in that case set + # g_out_bounds to 0s + if len(nz[0]) == 0: + g_out_bounds[p, a, g_in, :] = 0 + else: g_out_bounds[p, a, g_in, 0] = nz[0][0] g_out_bounds[p, a, g_in, 1] = nz[0][-1] + # Now create the flattened scatter matrix array + flat_scatt = [] + for p in range(Np): + for a in range(Na): + if self.representation == 'isotropic': + matrix = self._scatter_matrix[i][:, :, :] + elif self.representation == 'angle': + matrix = self._scatter_matrix[i][p, a, :, :, :] + for g_in in range(G): + for g_out in range(g_out_bounds[p, a, g_in, 0], + g_out_bounds[p, a, g_in, 1] + 1): + for l in range(len(matrix[g_in, g_out, :])): + flat_scatt.append(matrix[g_in, g_out, l]) + + # And write it. + scatt_grp = xs_grp.create_group('scatter_data') + scatt_grp.create_dataset("scatter_matrix", + data=np.array(flat_scatt)) + + # Repeat for multiplicity + if self._multiplicity_matrix[i] is not None: + # Now create the flattened scatter matrix array - flat_scatt = [] + flat_mult = [] for p in range(Np): for a in range(Na): - matrix = self._scatter_matrix[i][p, a, :, :, :] + if self.representation == 'isotropic': + matrix = self._multiplicity_matrix[i][:, :] + elif self.representation == 'angle': + matrix = self._multiplicity_matrix[i][p, a, :, :] for g_in in range(G): for g_out in range(g_out_bounds[p, a, g_in, 0], g_out_bounds[p, a, g_in, 1] + 1): - for l in range(len(matrix[:, g_in, g_out])): - flat_scatt.append(matrix[l, g_in, g_out]) + flat_mult.append(matrix[g_in, g_out]) # And write it. - scatt_grp = xs_grp.create_group('scatter_data') - scatt_grp.create_dataset("scatter_matrix", - data=np.array(flat_scatt)) + scatt_grp.create_dataset("multiplicity_matrix", + data=np.array(flat_mult)) - # Repeat for multiplicity - if self._multiplicity_matrix[i] is not None: - - # Now create the flattened scatter matrix array - flat_mult = [] - for p in range(Np): - for a in range(Na): - matrix = self._multiplicity_matrix[i][p, a, :, :] - for g_in in range(G): - for g_out in range(g_out_bounds[p, a, g_in, 0], - g_out_bounds[p, a, g_in, 1] + 1): - flat_mult.append(matrix[g_in, g_out]) - - # And write it. - scatt_grp.create_dataset("multiplicity_matrix", - data=np.array(flat_mult)) - - # And finally, adjust g_out_bounds for 1-based group counting - # and write it. - g_out_bounds[:, :, :, :] += 1 + # And finally, adjust g_out_bounds for 1-based group counting + # and write it. + g_out_bounds[:, :, :, :] += 1 + if self.representation == 'isotropic': + scatt_grp.create_dataset("g_min", data=g_out_bounds[0, 0, :, 0]) + scatt_grp.create_dataset("g_max", data=g_out_bounds[0, 0, :, 1]) + elif self.representation == 'angle': scatt_grp.create_dataset("g_min", data=g_out_bounds[:, :, :, 0]) scatt_grp.create_dataset("g_max", data=g_out_bounds[:, :, :, 1]) @@ -1834,6 +1840,143 @@ class XSdata(object): xs_grp.create_dataset("inverse-velocity", data=self._inverse_velocity[i]) + @classmethod + def from_hdf5(cls, group, name, energy_groups, num_delayed_groups): + """Generate XSdata object from an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + name : str + Name of the mgxs data set. + energy_groups : openmc.mgxs.EnergyGroups + Energy group structure + num_delayed_groups : int + Number of delayed groups + + Returns + ------- + openmc.XSdata + Multi-group cross section data + + """ + + # Get a list of all the subgroups which will contain our temperature + # strings + subgroups = group.keys() + temperatures = [] + for subgroup in subgroups: + if subgroup != 'kTs': + temperatures.append(subgroup) + + # To ensure the actual floating point temperature used when creating + # the new library is consistent with that used when originally creating + # the file, get the floating point temperatures straight from the kTs + # group. + kTs_group = group['kTs'] + float_temperatures = [] + for temperature in temperatures: + kT = kTs_group[temperature].value + float_temperatures.append(kT / openmc.data.K_BOLTZMANN) + + attrs = group.attrs.keys() + if 'representation' in attrs: + representation = group.attrs['representation'].decode() + else: + representation = 'isotropic' + + data = cls(name, energy_groups, float_temperatures, representation, + num_delayed_groups) + + if 'scatter_format' in attrs: + data.scatter_format = group.attrs['scatter_format'].decode() + + # Get the remaining optional attributes + if 'atomic_weight_ratio' in attrs: + data.atomic_weight_ratio = group.attrs['atomic_weight_ratio'] + if 'order' in attrs: + data.order = group.attrs['order'] + if data.representation == 'angle': + data.num_azimuthal = group.attrs['num_azimuthal'] + data.num_polar = group.attrs['num_polar'] + + # Read the temperature-dependent datasets + for temp, float_temp in zip(temperatures, float_temperatures): + xs_types = ['total', 'absorption', 'fission', 'kappa-fission', + 'chi', 'chi-prompt', 'chi-delayed', 'nu-fission', + 'prompt-nu-fission', 'delayed-nu-fission', 'beta', + 'decay rate', 'inverse-velocity'] + + temperature_group = group[temp] + + for xs_type in xs_types: + set_func = 'set_' + xs_type.replace(' ', '_').replace('-', '_') + if xs_type in temperature_group: + getattr(data, set_func)(temperature_group[xs_type].value, + float_temp) + + scatt_group = temperature_group['scatter_data'] + + # Get scatter matrix and 'un-flatten' it + g_max = scatt_group['g_max'] + g_min = scatt_group['g_min'] + flat_scatter = scatt_group['scatter_matrix'].value + scatter_matrix = np.zeros(data.xs_shapes["[G][G'][Order]"]) + G = data.energy_groups.num_groups + if data.representation == 'isotropic': + Np = 1 + Na = 1 + elif data.representation == 'angle': + Np = data.num_polar + Na = data.num_azimuthal + flat_index = 0 + for p in range(Np): + for a in range(Na): + for g_in in range(G): + if data.representation == 'isotropic': + g_mins = g_min[g_in] + g_maxs = g_max[g_in] + elif data.representation == 'angle': + g_mins = g_min[p, a, g_in] + g_maxs = g_max[p, a, g_in] + for g_out in range(g_mins - 1, g_maxs): + for ang in range(data.num_orders): + if data.representation == 'isotropic': + scatter_matrix[g_in, g_out, ang] = \ + flat_scatter[flat_index] + elif data.representation == 'angle': + scatter_matrix[p, a, g_in, g_out, ang] = \ + flat_scatter[flat_index] + flat_index += 1 + data.set_scatter_matrix(scatter_matrix, float_temp) + + # Repeat for multiplicity + if 'multiplicity_matrix' in scatt_group: + flat_mult = scatt_group['multiplicity_matrix'].value + mult_matrix = np.zeros(data.xs_shapes["[G][G']"]) + flat_index = 0 + for p in range(Np): + for a in range(Na): + for g_in in range(G): + if data.representation == 'isotropic': + g_mins = g_min[g_in] + g_maxs = g_max[g_in] + elif data.representation == 'angle': + g_mins = g_min[p, a, g_in] + g_maxs = g_max[p, a, g_in] + for g_out in range(g_mins - 1, g_maxs): + if data.representation == 'isotropic': + mult_matrix[g_in, g_out] = \ + flat_mult[flat_index] + elif data.representation == 'angle': + mult_matrix[p, a, g_in, g_out] = \ + flat_mult[flat_index] + flat_index += 1 + data.set_multiplicity_matrix(mult_matrix, float_temp) + + return data + class MGXSLibrary(object): """Multi-Group Cross Sections file used for an OpenMC simulation. @@ -1870,14 +2013,14 @@ class MGXSLibrary(object): def num_delayed_groups(self): return self._num_delayed_groups - @property - def temperatures(self): - return self._temperatures - @property def xsdatas(self): return self._xsdatas + @property + def names(self): + return [xsdata.name for xsdata in self.xsdatas] + @energy_groups.setter def energy_groups(self, energy_groups): check_type('energy groups', energy_groups, openmc.mgxs.EnergyGroups) @@ -1885,7 +2028,7 @@ class MGXSLibrary(object): @num_delayed_groups.setter def num_delayed_groups(self, num_delayed_groups): - check_type('num_delayed_groups', num_delayed_groups, int) + check_type('num_delayed_groups', num_delayed_groups, Integral) check_greater_than('num_delayed_groups', num_delayed_groups, 0, equality=True) check_less_than('num_delayed_groups', num_delayed_groups, @@ -1914,7 +2057,7 @@ class MGXSLibrary(object): self._xsdatas.append(xsdata) def add_xsdatas(self, xsdatas): - """Add multiple xsdatas to the file. + """Add multiple XSdatas to the file. Parameters ---------- @@ -1945,6 +2088,27 @@ class MGXSLibrary(object): self._xsdatas.remove(xsdata) + def get_by_name(self, name): + """Access the XSdata objects by name + + Parameters + ---------- + name : str + Name of openmc.XSdata object to obtain + + Returns + ------- + result : openmc.XSdata or None + Provides the matching XSdata object or None, if not found + + """ + check_type("name", name, str) + result = None + for xsdata in self.xsdatas: + if name == xsdata.name: + result = xsdata + return result + def export_to_hdf5(self, filename='mgxs.h5'): """Create an hdf5 file that can be used for a simulation. @@ -1967,3 +2131,43 @@ class MGXSLibrary(object): xsdata.to_hdf5(file) file.close() + + @classmethod + def from_hdf5(cls, filename=None): + """Generate an MGXS Library from an HDF5 group or file + Parameters + ---------- + filename : str, optional + Name of HDF5 file containing MGXS data. Default is None. + If not provided, the value of the OPENMC_MG_CROSS_SECTIONS + environmental variable will be used + Returns + ------- + openmc.MGXSLibrary + Multi-group cross section data object. + """ + + # If filename is None, get the cross sections from the + # OPENMC_CROSS_SECTIONS environment variable + if filename is None: + filename = os.environ.get('OPENMC_MG_CROSS_SECTIONS') + + # Check to make sure there was an environmental variable. + if filename is None: + raise ValueError("Either path or OPENMC_MG_CROSS_SECTIONS " + "environmental variable must be set") + + check_type('filename', filename, str) + file = h5py.File(filename, 'r') + + group_structure = file.attrs['group structure'] + num_delayed_groups = file.attrs['delayed_groups'] + energy_groups = openmc.mgxs.EnergyGroups(group_structure) + data = cls(energy_groups, num_delayed_groups) + + for group_name, group in file.items(): + data.add_xsdata(openmc.XSdata.from_hdf5(group, group_name, + energy_groups, + num_delayed_groups)) + + return data diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 55afd4e9c..03062ee0e 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -1,10 +1,8 @@ -from numbers import Integral -import sys import warnings from six import string_types -from openmc.checkvalue import check_type +import openmc.checkvalue as cv class Nuclide(object): @@ -72,7 +70,7 @@ class Nuclide(object): @name.setter def name(self, name): - check_type('name', name, string_types) + cv.check_type('name', name, string_types) self._name = name if '-' in name: diff --git a/openmc/plotter.py b/openmc/plotter.py new file mode 100644 index 000000000..aa4317d30 --- /dev/null +++ b/openmc/plotter.py @@ -0,0 +1,542 @@ +from numbers import Integral, Real +from six import string_types +from itertools import chain + +import numpy as np + +import openmc.checkvalue as cv +import openmc.data + +# Supported keywords for xs plotting +PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', + 'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity', + 'slowing-down power', 'damage'] + +# Special MT values +UNITY_MT = -1 +XI_MT = -2 + +# MTs to combine to generate associated plot_types +_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt != 27] +PLOT_TYPES_MT = {'total': openmc.data.SUM_RULES[1], + 'scatter': [2] + _INELASTIC, + 'elastic': [2], + 'inelastic': _INELASTIC, + 'fission': [18], + 'absorption': [27], 'capture': [101], + 'nu-fission': [18], + 'nu-scatter': [2] + _INELASTIC, + 'unity': [UNITY_MT], + 'slowing-down power': [2] + _INELASTIC + [XI_MT], + 'damage': [444]} +# Operations to use when combining MTs the first np.add is used in reference +# to zero +PLOT_TYPES_OP = {'total': (np.add,), + 'scatter': (np.add,) * (len(PLOT_TYPES_MT['scatter']) - 1), + 'elastic': (), + 'inelastic': (np.add,) * (len(PLOT_TYPES_MT['inelastic']) - 1), + 'fission': (), 'absorption': (), + 'capture': (), 'nu-fission': (), + 'nu-scatter': (np.add,) * (len(PLOT_TYPES_MT['nu-scatter']) - 1), + 'unity': (), + 'slowing-down power': + (np.add,) * (len(PLOT_TYPES_MT['slowing-down power']) - 2) + (np.multiply,), + 'damage': ()} + +# Types of plots to plot linearly in y +PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter', + 'nu-fission / absorption', 'fission / absorption'} + + +def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, + sab_name=None, cross_sections=None, enrichment=None, **kwargs): + """Creates a figure of continuous-energy cross sections for this item + + Parameters + ---------- + this : openmc.Element, openmc.Nuclide, or openmc.Material + Object to source data from + types : Iterable of values of PLOT_TYPES + The type of cross sections to include in the plot. + divisor_types : Iterable of values of PLOT_TYPES, optional + Cross section types which will divide those produced by types + before plotting. A type of 'unity' can be used to effectively not + divide some types. + temperature : float, optional + Temperature in Kelvin to plot. If not specified, a default + temperature of 294K will be plotted. Note that the nearest + temperature in the library for each nuclide will be used as opposed + to using any interpolation. + axis : matplotlib.axes, optional + A previously generated axis to use for plotting. If not specified, + a new axis and figure will be generated. + sab_name : str, optional + Name of S(a,b) library to apply to MT=2 data when applicable; only used + for items which are instances of openmc.Element or openmc.Nuclide + cross_sections : str, optional + Location of cross_sections.xml file. Default is None. + enrichment : float, optional + Enrichment for U235 in weight percent. For example, input 4.95 for + 4.95 weight percent enriched U. Default is None. This is only used for + items which are instances of openmc.Element + **kwargs + All keyword arguments are passed to + :func:`matplotlib.pyplot.figure`. + + Returns + ------- + fig : matplotlib.figure.Figure + If axis is None, then a Matplotlib Figure of the generated + cross section will be returned. Otherwise, a value of + None will be returned as the figure and axes have already been + generated. + + """ + + from matplotlib import pyplot as plt + + if isinstance(this, openmc.Nuclide): + data_type = 'nuclide' + elif isinstance(this, openmc.Element): + data_type = 'element' + elif isinstance(this, openmc.Material): + data_type = 'material' + else: + raise TypeError("Invalid type for plotting") + + E, data = calculate_xs(this, types, temperature, sab_name, cross_sections, + enrichment) + + if divisor_types: + cv.check_length('divisor types', divisor_types, len(types), + len(types)) + Ediv, data_div = calculate_xs(this, divisor_types, temperature, + sab_name, cross_sections, enrichment) + + # Create a new union grid, interpolate data and data_div on to that + # grid, and then do the actual division + Enum = E[:] + E = np.union1d(Enum, Ediv) + data_new = np.zeros((len(types), len(E))) + + for line in range(len(types)): + data_new[line, :] = \ + np.divide(np.interp(E, Enum, data[line, :]), + np.interp(E, Ediv, data_div[line, :])) + if divisor_types[line] != 'unity': + types[line] = types[line] + ' / ' + divisor_types[line] + data = data_new + + # Generate the plot + if axis is None: + fig = plt.figure(**kwargs) + ax = fig.add_subplot(111) + else: + fig = None + ax = axis + # Set to loglog or semilogx depending on if we are plotting a data + # type which we expect to vary linearly + if set(types).issubset(PLOT_TYPES_LINEAR): + plot_func = ax.semilogx + else: + plot_func = ax.loglog + # Plot the data + for i in range(len(data)): + data[i, :] = np.nan_to_num(data[i, :]) + if np.sum(data[i, :]) > 0.: + plot_func(E, data[i, :], label=types[i]) + + ax.set_xlabel('Energy [eV]') + if divisor_types: + if data_type == 'nuclide': + ylabel = 'Nuclidic Microscopic Data' + elif data_type == 'element': + ylabel = 'Elemental Microscopic Data' + elif data_type == 'material': + ylabel = 'Macroscopic Data' + else: + if data_type == 'nuclide': + ylabel = 'Microscopic Cross Section [b]' + elif data_type == 'element': + ylabel = 'Elemental Cross Section [b]' + elif data_type == 'material': + ylabel = 'Macroscopic Cross Section [1/cm]' + ax.set_ylabel(ylabel) + ax.legend(loc='best') + # Set to the most likely expected range + ax.set_xlim((1.E-5, 20.E6)) + if this.name is not None: + ax.set_title('Cross Section for ' + this.name) + + return fig + + +def calculate_xs(this, types, temperature=294., sab_name=None, + cross_sections=None, enrichment=None): + """Calculates continuous-energy cross sections of a requested type + + Parameters + ---------- + this : openmc.Element, openmc.Nuclide, or openmc.Material + Object to source data from + types : Iterable of values of PLOT_TYPES + The type of cross sections to calculate + temperature : float, optional + Temperature in Kelvin to plot. If not specified, a default + temperature of 294K will be plotted. Note that the nearest + temperature in the library for each nuclide will be used as opposed + to using any interpolation. + sab_name : str, optional + Name of S(a,b) library to apply to MT=2 data when applicable. + cross_sections : str, optional + Location of cross_sections.xml file. Default is None. + enrichment : float, optional + Enrichment for U235 in weight percent. For example, input 4.95 for + 4.95 weight percent enriched U. Default is None + (natural composition). + + Returns + ------- + energy_grid : numpy.ndarray + Energies at which cross sections are calculated, in units of eV + data : numpy.ndarray + Cross sections calculated at the energy grid described by energy_grid + + """ + + # Check types + cv.check_type('temperature', temperature, Real) + if sab_name: + cv.check_type('sab_name', sab_name, string_types) + if enrichment: + cv.check_type('enrichment', enrichment, Real) + + if isinstance(this, openmc.Nuclide): + energy_grid, xs = _calculate_xs_nuclide(this, types, temperature, + sab_name, cross_sections) + # Convert xs (Iterable of Callable) to a grid of cross section values + # calculated on @ the points in energy_grid for consistency with the + # element and material functions. + data = np.zeros((len(types), len(energy_grid))) + for line in range(len(types)): + data[line, :] = xs[line](energy_grid) + elif isinstance(this, openmc.Element): + energy_grid, data = _calculate_xs_elem_mat(this, types, temperature, + cross_sections, sab_name, + enrichment) + elif isinstance(this, openmc.Material): + energy_grid, data = _calculate_xs_elem_mat(this, types, temperature, + cross_sections) + else: + raise TypeError("Invalid type") + + return energy_grid, data + + +def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, + cross_sections=None): + """Calculates continuous-energy cross sections of a requested type + + Parameters + ---------- + this : openmc.Nuclide + Nuclide object to source data from + types : Iterable of str or Integral + The type of cross sections to calculate; values can either be those + in openmc.PLOT_TYPES or integers which correspond to reaction + channel (MT) numbers. + temperature : float, optional + Temperature in Kelvin to plot. If not specified, a default + temperature of 294K will be plotted. Note that the nearest + temperature in the library for each nuclide will be used as opposed + to using any interpolation. + sab_name : str, optional + Name of S(a,b) library to apply to MT=2 data when applicable. + cross_sections : str, optional + Location of cross_sections.xml file. Default is None. + + Returns + ------- + energy_grid : numpy.ndarray + Energies at which cross sections are calculated, in units of eV + data : Iterable of Callable + Requested cross section functions + + """ + + # Parse the types + mts = [] + ops = [] + yields = [] + for line in types: + if line in PLOT_TYPES: + mts.append(PLOT_TYPES_MT[line]) + if line.startswith('nu'): + yields.append(True) + else: + yields.append(False) + ops.append(PLOT_TYPES_OP[line]) + else: + # Not a built-in type, we have to parse it ourselves + cv.check_type('MT in types', line, Integral) + cv.check_greater_than('MT in types', line, 0) + mts.append((line,)) + ops.append(()) + yields.append(False) + + # Load the library + library = openmc.data.DataLibrary.from_xml(cross_sections) + + # Convert temperature to format needed for access in the library + strT = "{}K".format(int(round(temperature))) + T = temperature + + # Now we can create the data sets to be plotted + energy_grid = [] + xs = [] + lib = library.get_by_material(this.name) + if lib is not None: + nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) + # Obtain the nearest temperature + if strT in nuc.temperatures: + nucT = strT + else: + data_Ts = nuc.temperatures + for t in range(len(data_Ts)): + # Take off the "K" and convert to a float + data_Ts[t] = float(data_Ts[t][:-1]) + min_delta = float('inf') + closest_t = -1 + for t in data_Ts: + if abs(data_Ts[t] - T) < min_delta: + closest_t = t + nucT = "{}K".format(int(round(data_Ts[closest_t]))) + + # Prep S(a,b) data if needed + if sab_name: + sab = openmc.data.ThermalScattering.from_hdf5(sab_name) + # Obtain the nearest temperature + if strT in sab.temperatures: + sabT = strT + else: + data_Ts = sab.temperatures + for t in range(len(data_Ts)): + # Take off the "K" and convert to a float + data_Ts[t] = float(data_Ts[t][:-1]) + min_delta = np.finfo(np.float64).max + closest_t = -1 + for t in data_Ts: + if abs(data_Ts[t] - T) < min_delta: + closest_t = t + sabT = "{}K".format(int(round(data_Ts[closest_t]))) + + # Create an energy grid composed the S(a,b) and + # the nuclide's grid + grid = nuc.energy[nucT] + sab_Emax = 0. + sab_funcs = [] + if sab.elastic_xs: + elastic = sab.elastic_xs[sabT] + if isinstance(elastic, openmc.data.CoherentElastic): + grid = np.union1d(grid, elastic.bragg_edges) + if elastic.bragg_edges[-1] > sab_Emax: + sab_Emax = elastic.bragg_edges[-1] + elif isinstance(elastic, openmc.data.Tabulated1D): + grid = np.union1d(grid, elastic.x) + if elastic.x[-1] > sab_Emax: + sab_Emax = elastic.x[-1] + sab_funcs.append(elastic) + if sab.inelastic_xs: + inelastic = sab.inelastic_xs[sabT] + grid = np.union1d(grid, inelastic.x) + if inelastic.x[-1] > sab_Emax: + sab_Emax = inelastic.x[-1] + sab_funcs.append(inelastic) + energy_grid = grid + else: + energy_grid = nuc.energy[nucT] + + for i, mt_set in enumerate(mts): + # Get the reaction xs data from the nuclide + funcs = [] + op = ops[i] + for mt in mt_set: + if mt == 2: + if sab_name: + # Then we need to do a piece-wise function of + # The S(a,b) and non-thermal data + sab_sum = openmc.data.Sum(sab_funcs) + pw_funcs = openmc.data.Regions1D( + [sab_sum, nuc[mt].xs[nucT]], + [sab_Emax]) + funcs.append(pw_funcs) + else: + funcs.append(nuc[mt].xs[nucT]) + elif mt in nuc: + if yields[i]: + # Get the total yield first if available. This will be + # used primarily for fission. + for prod in chain(nuc[mt].products, + nuc[mt].derived_products): + if prod.particle == 'neutron' and \ + prod.emission_mode == 'total': + func = openmc.data.Combination( + [nuc[mt].xs[nucT], prod.yield_], + [np.multiply]) + funcs.append(func) + break + else: + # Total doesn't exist so we have to create from + # prompt and delayed. This is used for scatter + # multiplication. + func = None + for prod in chain(nuc[mt].products, + nuc[mt].derived_products): + if prod.particle == 'neutron' and \ + prod.emission_mode != 'total': + if func: + func = openmc.data.Combination( + [prod.yield_, func], [np.add]) + else: + func = prod.yield_ + if func: + funcs.append(openmc.data.Combination( + [func, nuc[mt].xs[nucT]], [np.multiply])) + else: + # If func is still None, then there were no + # products. In that case, assume the yield is + # one as its not provided for some summed + # reactions like MT=4 + funcs.append(nuc[mt].xs[nucT]) + else: + funcs.append(nuc[mt].xs[nucT]) + elif mt == UNITY_MT: + funcs.append(lambda x: 1.) + elif mt == XI_MT: + awr = nuc.atomic_weight_ratio + alpha = ((awr - 1.) / (awr + 1.))**2 + xi = 1. + alpha * np.log(alpha) / (1. - alpha) + funcs.append(lambda x: xi) + else: + funcs.append(lambda x: 0.) + xs.append(openmc.data.Combination(funcs, op)) + else: + raise ValueError(this.name + " not in library") + + return energy_grid, xs + + +def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None, + sab_name=None, enrichment=None): + """Calculates continuous-energy cross sections of a requested type + + Parameters + ---------- + this : {openmc.Material, openmc.Element} + Object to source data from + types : Iterable of values of PLOT_TYPES + The type of cross sections to calculate + temperature : float, optional + Temperature in Kelvin to plot. If not specified, a default + temperature of 294K will be plotted. Note that the nearest + temperature in the library for each nuclide will be used as opposed + to using any interpolation. + cross_sections : str, optional + Location of cross_sections.xml file. Default is None. + sab_name : str, optional + Name of S(a,b) library to apply to MT=2 data when applicable. + enrichment : float, optional + Enrichment for U235 in weight percent. For example, input 4.95 for + 4.95 weight percent enriched U. Default is None + (natural composition). + + Returns + ------- + energy_grid : numpy.ndarray + Energies at which cross sections are calculated, in units of eV + data : numpy.ndarray + Cross sections calculated at the energy grid described by energy_grid + + """ + + if isinstance(this, openmc.Material): + if this.temperature is not None: + T = this.temperature + else: + T = temperature + else: + T = temperature + + # Load the library + library = openmc.data.DataLibrary.from_xml(cross_sections) + + if isinstance(this, openmc.Material): + # Expand elements in to nuclides with atomic densities + nuclides = this.get_nuclide_atom_densities() + # For ease of processing split out the nuclide and its fraction + nuc_fractions = {nuclide[1][0].name: nuclide[1][1] + for nuclide in nuclides.items()} + # Create a dict of [nuclide name] = nuclide object to carry forward + # with a common nuclides format between openmc.Material and + # openmc.Element objects + nuclides = {nuclide[1][0].name: nuclide[1][0] + for nuclide in nuclides.items()} + else: + # Expand elements in to nuclides with atomic densities + nuclides = this.expand(1., 'ao', enrichment=enrichment, + cross_sections=cross_sections) + # For ease of processing split out the nuclide and its fraction + nuc_fractions = {nuclide[0].name: nuclide[1] for nuclide in nuclides} + # Create a dict of [nuclide name] = nuclide object to carry forward + # with a common nuclides format between openmc.Material and + # openmc.Element objects + nuclides = {nuclide[0].name: nuclide[0] for nuclide in nuclides} + + # Identify the nuclides which have S(a,b) data + sabs = {} + for nuclide in nuclides.items(): + sabs[nuclide[0]] = None + if isinstance(this, openmc.Material): + for sab_name in this._sab: + sab = openmc.data.ThermalScattering.from_hdf5( + library.get_by_material(sab_name)['path']) + for nuc in sab.nuclides: + sabs[nuc] = library.get_by_material(sab_name)['path'] + else: + if sab_name: + sab = openmc.data.ThermalScattering.from_hdf5(sab_name) + for nuc in sab.nuclides: + sabs[nuc] = library.get_by_material(sab_name)['path'] + + # Now we can create the data sets to be plotted + xs = {} + E = [] + for nuclide in nuclides.items(): + name = nuclide[0] + nuc = nuclide[1] + sab_tab = sabs[name] + temp_E, temp_xs = calculate_xs(nuc, types, T, sab_tab, cross_sections) + E.append(temp_E) + # Since the energy grids are different, store the cross sections as + # a tabulated function so they can be calculated on any grid needed. + xs[name] = [openmc.data.Tabulated1D(temp_E, temp_xs[line]) + for line in range(len(types))] + + # Condense the data for every nuclide + # First create a union energy grid + energy_grid = E[0] + for grid in E[1:]: + energy_grid = np.union1d(energy_grid, grid) + + # Now we can combine all the nuclidic data + data = np.zeros((len(types), len(energy_grid))) + for line in range(len(types)): + if types[line] == 'unity': + data[line, :] = 1. + else: + for nuclide in nuclides.items(): + name = nuclide[0] + data[line, :] += (nuc_fractions[name] * + xs[name][line](energy_grid)) + + return energy_grid, data diff --git a/openmc/summary.py b/openmc/summary.py index 37509ec37..b283ed502 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -387,7 +387,7 @@ class Summary(object): # Set the distribcell offsets for the lattice if offsets is not None: - lattice.offsets = offsets[:, ::-1, :] + lattice.offsets = offsets # Add the Lattice to the global dictionary of all Lattices self.lattices[index] = lattice diff --git a/openmc/tallies.py b/openmc/tallies.py index fe690e2d6..35fb32175 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1317,14 +1317,15 @@ class Tally(object): # Create list of 2-tuples for energy boundary bins elif isinstance(self_filter, (openmc.EnergyFilter, - openmc.EnergyoutFilter)): + openmc.EnergyoutFilter, openmc.MuFilter, + openmc.PolarFilter, openmc.AzimuthalFilter)): bins = [] for k in range(self_filter.num_bins): bins.append((self_filter.bins[k], self_filter.bins[k+1])) # Create list of cell instance IDs for distribcell Filters elif isinstance(self_filter, openmc.DistribcellFilter): - bins = np.arange(self_filter.num_bins) + bins = [b for b in range(self_filter.num_bins)] # EnergyFunctionFilters don't have bins so just add a None elif isinstance(self_filter, openmc.EnergyFunctionFilter): @@ -2119,14 +2120,14 @@ class Tally(object): # Construct lists of tuples for the bins in each of the two filters filters = [type(filter1), type(filter2)] if isinstance(filter1, openmc.DistribcellFilter): - filter1_bins = np.arange(filter1.num_bins) + filter1_bins = [b for b in range(filter1.num_bins)] elif isinstance(filter1, openmc.EnergyFunctionFilter): filter1_bins = [None] else: filter1_bins = [filter1.get_bin(i) for i in range(filter1.num_bins)] if isinstance(filter2, openmc.DistribcellFilter): - filter2_bins = np.arange(filter2.num_bins) + filter2_bins = [b for b in range(filter2.num_bins)] elif isinstance(filter2, openmc.EnergyFunctionFilter): filter2_bins = [None] else: @@ -2975,8 +2976,16 @@ class Tally(object): # Sum across the bins in the user-specified filter for i, self_filter in enumerate(self.filters): if isinstance(self_filter, filter_type): + shape = mean.shape mean = np.take(mean, indices=bin_indices, axis=i) std_dev = np.take(std_dev, indices=bin_indices, axis=i) + + # NumPy take introduces a new dimension in output array + # for some special cases that must be removed + if len(mean.shape) > len(shape): + mean = np.squeeze(mean, axis=i) + std_dev = np.squeeze(std_dev, axis=i) + mean = np.sum(mean, axis=i, keepdims=True) std_dev = np.sum(std_dev**2, axis=i, keepdims=True) std_dev = np.sqrt(std_dev) @@ -3124,10 +3133,18 @@ class Tally(object): # Average across the bins in the user-specified filter for i, self_filter in enumerate(self.filters): if isinstance(self_filter, filter_type): + shape = mean.shape mean = np.take(mean, indices=bin_indices, axis=i) std_dev = np.take(std_dev, indices=bin_indices, axis=i) - mean = np.mean(mean, axis=i, keepdims=True) - std_dev = np.mean(std_dev**2, axis=i, keepdims=True) + + # NumPy take introduces a new dimension in output array + # for some special cases that must be removed + if len(mean.shape) > len(shape): + mean = np.squeeze(mean, axis=i) + std_dev = np.squeeze(std_dev, axis=i) + + mean = np.nanmean(mean, axis=i, keepdims=True) + std_dev = np.nanmean(std_dev**2, axis=i, keepdims=True) std_dev /= len(bin_indices) std_dev = np.sqrt(std_dev) @@ -3151,8 +3168,8 @@ class Tally(object): axis_index = self.num_filters mean = np.take(mean, indices=nuclide_bins, axis=axis_index) std_dev = np.take(std_dev, indices=nuclide_bins, axis=axis_index) - mean = np.mean(mean, axis=axis_index, keepdims=True) - std_dev = np.mean(std_dev**2, axis=axis_index, keepdims=True) + mean = np.nanmean(mean, axis=axis_index, keepdims=True) + std_dev = np.nanmean(std_dev**2, axis=axis_index, keepdims=True) std_dev /= len(nuclide_bins) std_dev = np.sqrt(std_dev) @@ -3170,8 +3187,8 @@ class Tally(object): axis_index = self.num_filters + 1 mean = np.take(mean, indices=score_bins, axis=axis_index) std_dev = np.take(std_dev, indices=score_bins, axis=axis_index) - mean = np.sum(mean, axis=axis_index, keepdims=True) - std_dev = np.sum(std_dev**2, axis=axis_index, keepdims=True) + mean = np.nanmean(mean, axis=axis_index, keepdims=True) + std_dev = np.nanmean(std_dev**2, axis=axis_index, keepdims=True) std_dev /= len(score_bins) std_dev = np.sqrt(std_dev) diff --git a/scripts/openmc-update-mgxs b/scripts/openmc-update-mgxs index cad48d0e3..54a553f77 100755 --- a/scripts/openmc-update-mgxs +++ b/scripts/openmc-update-mgxs @@ -23,12 +23,10 @@ optional arguments: from __future__ import print_function import os -from shutil import move import warnings import xml.etree.ElementTree as ET import argparse -import h5py import numpy as np import openmc.mgxs_library @@ -51,7 +49,7 @@ def parse_args(): if args['output'] == '': filename = args['input'].name - extension = filenameos.path.splitext() + extension = os.path.splitext(filename) if extension == '.xml': filename = filename[:filename.rfind('.')] + '.h5' args['output'] = filename @@ -84,6 +82,8 @@ if __name__ == '__main__': temp = tree.find('group_structure').text.strip() temp = np.array(temp.split()) group_structure = temp.astype(np.float) + # Convert from MeV to eV + group_structure *= 1.e6 energy_groups = openmc.mgxs.EnergyGroups(group_structure) temp = tree.find('inverse-velocity') if temp is not None: @@ -103,7 +103,7 @@ if __name__ == '__main__': temperature = get_data(xsdata_elem, 'kT') if temperature is not None: temperature = \ - float(temperature) / openmc.data.K_BOLTZMANN + float(temperature) / openmc.data.K_BOLTZMANN * 1.E6 else: temperature = 294. temperatures = [temperature] @@ -163,7 +163,7 @@ if __name__ == '__main__': if temp is not None: temp = np.array(temp.split(), dtype=float) total = temp.astype(np.float) - total.shape = xsd[i].vector_shape + total.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_total(total, temperature) if inverse_velocity is not None: @@ -172,41 +172,55 @@ if __name__ == '__main__': temp = get_data(xsdata_elem, 'absorption') temp = np.array(temp.split()) absorption = temp.astype(np.float) - absorption.shape = xsd[i].vector_shape + absorption.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_absorption(absorption, temperature) temp = get_data(xsdata_elem, 'scatter') temp = np.array(temp.split()) scatter = temp.astype(np.float) - scatter.shape = xsd[i].pn_matrix_shape + # This is now a flattened-array of something that started with a + # shape of [Order][G][G']; we need to unflatten and then switch the + # ordering + in_shape = (order_dim, energy_groups.num_groups, + energy_groups.num_groups) + if representation == 'angle': + in_shape = (n_pol, n_azi) + in_shape + scatter.shape = in_shape + scatter = np.swapaxes(scatter, 2, 3) + scatter = np.swapaxes(scatter, 3, 4) + else: + scatter.shape = in_shape + scatter = np.swapaxes(scatter, 0, 1) + scatter = np.swapaxes(scatter, 1, 2) + xsd[i].set_scatter_matrix(scatter, temperature) temp = get_data(xsdata_elem, 'multiplicity') if temp is not None: temp = np.array(temp.split()) multiplicity = temp.astype(np.float) - multiplicity.shape = xsd[i].matrix_shape + multiplicity.shape = xsd[i].xs_shapes["[G][G']"] xsd[i].set_multiplicity_matrix(multiplicity, temperature) temp = get_data(xsdata_elem, 'fission') if temp is not None: temp = np.array(temp.split()) fission = temp.astype(np.float) - fission.shape = xsd[i].vector_shape + fission.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_fission(fission, temperature) temp = get_data(xsdata_elem, 'kappa_fission') if temp is not None: temp = np.array(temp.split()) kappa_fission = temp.astype(np.float) - kappa_fission.shape = xsd[i].vector_shape + kappa_fission.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_kappa_fission(kappa_fission, temperature) temp = get_data(xsdata_elem, 'chi') if temp is not None: temp = np.array(temp.split()) chi = temp.astype(np.float) - chi.shape = xsd[i].vector_shape + chi.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_chi(chi, temperature) else: chi = None @@ -216,9 +230,9 @@ if __name__ == '__main__': temp = np.array(temp.split()) nu_fission = temp.astype(np.float) if chi is not None: - nu_fission.shape = xsd[i].vector_shape + nu_fission.shape = xsd[i].xs_shapes['[G]'] else: - nu_fission.shape = xsd[i].matrix_shape + nu_fission.shape = xsd[i].xs_shapes["[G][G']"] xsd[i].set_nu_fission(nu_fission, temperature) # Build library as we go, but first we have enough to initialize it diff --git a/setup.py b/setup.py index 0885c28a8..befb9c0d2 100755 --- a/setup.py +++ b/setup.py @@ -43,6 +43,7 @@ if have_setuptools: # Optional dependencies 'extras_require': { + 'decay': ['uncertainties'], 'pandas': ['pandas>=0.17.0'], 'sparse' : ['scipy'], 'vtk': ['vtk', 'silomesh'], diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index d536a587e..0fe48f7c9 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -355,7 +355,7 @@ module mgxs_header if (attribute_exists(xs_id, "scatter_shape")) then call read_attribute(temp_str, xs_id, "scatter_shape") temp_str = trim(temp_str) - if (to_lower(temp_str) /= "[order][g][g']") then + if (to_lower(temp_str) /= "[g][g'][order]") then call fatal_error("Invalid scatter_shape option!") end if end if @@ -1995,8 +1995,8 @@ module mgxs_header order_dim = order + 1 end if - ! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for passing - ! to ScattData + ! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for + ! passing to ScattData allocate(input_scatt(energy_groups, this % n_azi, this % n_pol)) index = 1 diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index 6004cbd19..684be88b8 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -274,7 +274,7 @@ contains ! Get this by summing the un-normalized P0 coefficient in matrix ! over all outgoing groups do gin = 1, groups - this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1) + this % scattxs(gin) = sum(matrix(gin) % data(:, :)) end do allocate(energy(groups)) @@ -282,7 +282,7 @@ contains ! while also normalizing matrix itself (making CDF of f(mu=1)=1) do gin = 1, groups allocate(energy(gin) % data(gmin(gin):gmax(gin))) - do gout = 1, groups + do gout = gmin(gin), gmax(gin) norm = sum(matrix(gin) % data(:, gout)) energy(gin) % data(gout) = norm if (norm /= ZERO) then diff --git a/tests/1d_mgxs.h5 b/tests/1d_mgxs.h5 index e3b85c241..8f28e02cb 100644 Binary files a/tests/1d_mgxs.h5 and b/tests/1d_mgxs.h5 differ diff --git a/tests/test_asymmetric_lattice/results_true.dat b/tests/test_asymmetric_lattice/results_true.dat index dda6f953f..34c8faa64 100644 --- a/tests/test_asymmetric_lattice/results_true.dat +++ b/tests/test_asymmetric_lattice/results_true.dat @@ -1 +1 @@ -e18c2318bab6c42a263e5079fd796b1bee609e4274884fbc7bfd8b33e59aeb5e5a167da8e093f339f766815e007bd606c002939e3f730af523cbc9cf75c53faa \ No newline at end of file +b70886031e22db9e3f0332eac703a7356504750c1e90d7083ffd16b8884d00661d0e20c6d8bead3c93369b2e7c105ca3280c7858ca6a147fa6669a5d3d530461 \ No newline at end of file diff --git a/tests/test_mg_tallies/results_true.dat b/tests/test_mg_tallies/results_true.dat index a765c0ba0..8ab726bbb 100644 --- a/tests/test_mg_tallies/results_true.dat +++ b/tests/test_mg_tallies/results_true.dat @@ -1 +1 @@ -3d2ce1b8bdd558fe9f8560e8bb91455e1fedd80a47349344399e22f5b0fac1391f1721e7f12a42bb0c458fd6c87ebf9fa38fe9539fcf69292b21cebf8e5f8989 \ No newline at end of file +864328b2c4f3c4bfa9c80c756acbedac6fbdf3d502c03c2f2a3e7461b774729c47a55341c9515eff246e1bd445485f0d0a05a5712663ee499046afdbe0c77aef \ No newline at end of file diff --git a/tests/test_mgxs_library_distribcell/results_true.dat b/tests/test_mgxs_library_distribcell/results_true.dat index 2b43005aa..905a55815 100644 --- a/tests/test_mgxs_library_distribcell/results_true.dat +++ b/tests/test_mgxs_library_distribcell/results_true.dat @@ -1,79 +1,79 @@ - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.457353 0.010474 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.405649 0.015784 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.405641 0.015787 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.066556 0.00251 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.028979 0.002712 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.037577 0.001487 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.092377 0.003628 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 7.276707e+06 287579.26286 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.390797 0.008717 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.387332 0.014241 - avg(distribcell) group in group out nuclide moment mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P0 0.387009 0.014230 -1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P1 0.047179 0.004923 -2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P2 0.015713 0.003654 -3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P3 0.005378 0.003137 - avg(distribcell) group in group out nuclide moment mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P0 0.387332 0.014241 -1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P1 0.047187 0.004933 -2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P2 0.015727 0.003654 -3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total P3 0.005387 0.003141 - avg(distribcell) group in group out nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 1.000834 0.037242 - avg(distribcell) group in group out nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.094516 0.0059 - avg(distribcell) group out nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 1.0 0.080455 - avg(distribcell) group out nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 1.0 0.080541 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 5.139437e-07 2.133314e-08 - avg(distribcell) group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.091725 0.003604 - avg(distribcell) group in group out nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.093985 0.005872 - avg(distribcell) delayedgroup group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.000021 8.253907e-07 -1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 0.000112 4.284000e-06 -2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 0.000109 4.105197e-06 -3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.000252 9.271420e-06 -4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.000112 3.888625e-06 -5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 0.000047 1.625563e-06 - avg(distribcell) delayedgroup group out nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.0 0.000000 -1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 1.0 1.414214 -2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 1.0 1.414214 -3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.0 0.000000 -4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.0 0.000000 -5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 1.0 1.414214 - avg(distribcell) delayedgroup group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.000227 0.000012 -1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 0.001209 0.000061 -2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 0.001177 0.000059 -3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.002727 0.000135 -4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.001210 0.000058 -5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 0.000504 0.000024 - avg(distribcell) delayedgroup group in nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.000000 0.000000 -1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 0.032739 0.046300 -2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 total 0.120780 0.170809 -3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.000000 0.000000 -4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.000000 0.000000 -5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 2.853000 4.034751 - avg(distribcell) delayedgroup group in group out nuclide mean std. dev. -0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 1 total 0.000000 0.000000 -1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 1 total 0.000175 0.000175 -2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 1 total 0.000178 0.000178 -3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 1 total 0.000000 0.000000 -4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 1 total 0.000000 0.000000 -5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 1 total 0.000178 0.000178 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.457353 0.010474 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.405649 0.015784 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.405641 0.015787 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.066556 0.00251 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.028979 0.002712 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.037577 0.001487 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.092377 0.003628 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 7.276707e+06 287579.26286 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.390797 0.008717 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.387332 0.014241 + sum(distribcell) group in group out nuclide moment mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P0 0.387009 0.014230 +1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P1 0.047179 0.004923 +2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P2 0.015713 0.003654 +3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P3 0.005378 0.003137 + sum(distribcell) group in group out nuclide moment mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P0 0.387332 0.014241 +1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P1 0.047187 0.004933 +2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P2 0.015727 0.003654 +3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total P3 0.005387 0.003141 + sum(distribcell) group in group out nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 1.000834 0.037242 + sum(distribcell) group in group out nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.094516 0.0059 + sum(distribcell) group out nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 1.0 0.080455 + sum(distribcell) group out nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 1.0 0.080541 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 5.139437e-07 2.133314e-08 + sum(distribcell) group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 total 0.091725 0.003604 + sum(distribcell) group in group out nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.093985 0.005872 + sum(distribcell) delayedgroup group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.000021 8.253907e-07 +1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 0.000112 4.284000e-06 +2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 0.000109 4.105197e-06 +3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.000252 9.271420e-06 +4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.000112 3.888625e-06 +5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 0.000047 1.625563e-06 + sum(distribcell) delayedgroup group out nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.0 0.000000 +1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 1.0 1.414214 +2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 1.0 1.414214 +3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.0 0.000000 +4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.0 0.000000 +5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 1.0 1.414214 + sum(distribcell) delayedgroup group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.000227 0.000012 +1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 0.001209 0.000061 +2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 0.001177 0.000059 +3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.002727 0.000135 +4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.001210 0.000058 +5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 0.000504 0.000024 + sum(distribcell) delayedgroup group in nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 total 0.000000 0.000000 +1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 total 0.032739 0.046300 +2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 total 0.120780 0.170809 +3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 total 0.000000 0.000000 +4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 total 0.000000 0.000000 +5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 total 2.853000 4.034751 + sum(distribcell) delayedgroup group in group out nuclide mean std. dev. +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 1 1 1 total 0.000000 0.000000 +1 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 2 1 1 total 0.000175 0.000175 +2 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 3 1 1 total 0.000178 0.000178 +3 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 4 1 1 total 0.000000 0.000000 +4 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 5 1 1 total 0.000000 0.000000 +5 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13... 6 1 1 total 0.000178 0.000178 diff --git a/tests/test_tallies/inputs_true.dat b/tests/test_tallies/inputs_true.dat index 69472c8d6..a8b40554b 100644 --- a/tests/test_tallies/inputs_true.dat +++ b/tests/test_tallies/inputs_true.dat @@ -1 +1 @@ -b1f616bc0e342e3886b4b27eca22c4d1e55e445979254d5ba5f9ea01048f53d0e9a01c575745ae130eb34108080d2256ab9e4cb926ea0050891a44ae37ecdd88 \ No newline at end of file +ac118c5796593df0f53491920875bbe77921882253e7f24942e365e5812d726fb3d3bab86da02b3aad93f6f8f20a090125ffa51f9fdf9aef84009702f9cc363d \ No newline at end of file diff --git a/tests/test_tallies/results_true.dat b/tests/test_tallies/results_true.dat index 639e8baf3..1d6f6bfaf 100644 --- a/tests/test_tallies/results_true.dat +++ b/tests/test_tallies/results_true.dat @@ -1 +1 @@ -ceb1420d48e097185dd0798a495676c8e968c37c6e75fd6232137a4df8f6fae33d8232b95c68f6479815a014d20ccfc7611f96b22fcb7da159b42346e00cdf4e \ No newline at end of file +a8172f7492fcc69d21f53c20523ffd0ad2882a753a3e1d263d135153cc6b3ee8d8e38f8823a4a2eb3288d730579892875d2d19f85f083dd659f91414988115ea \ No newline at end of file diff --git a/tests/test_tallies/test_tallies.py b/tests/test_tallies/test_tallies.py index 9b48729ef..033c7b812 100644 --- a/tests/test_tallies/test_tallies.py +++ b/tests/test_tallies/test_tallies.py @@ -3,12 +3,14 @@ import os import sys sys.path.insert(0, os.pardir) + from testing_harness import PyAPITestHarness from openmc.filter import * from openmc import Mesh, Tally, Tallies from openmc.source import Source from openmc.stats import Box + class TalliesTestHarness(PyAPITestHarness): def _build_inputs(self): # Build default materials/geometry @@ -21,33 +23,27 @@ class TalliesTestHarness(PyAPITestHarness): self._input_set.settings.source = Source(space=Box( [-160, -160, -183], [160, 160, 183])) - azimuthal_bins = (-3.1416, -1.8850, -0.6283, 0.6283, 1.8850, 3.1416) - azimuthal_filter1 = AzimuthalFilter(azimuthal_bins) + azimuthal_bins = (-3.14159, -1.8850, -0.6283, 0.6283, 1.8850, 3.14159) + azimuthal_filter = AzimuthalFilter(azimuthal_bins) azimuthal_tally1 = Tally() - azimuthal_tally1.filters = [azimuthal_filter1] + azimuthal_tally1.filters = [azimuthal_filter] azimuthal_tally1.scores = ['flux'] azimuthal_tally1.estimator = 'tracklength' azimuthal_tally2 = Tally() - azimuthal_tally2.filters = [azimuthal_filter1] + azimuthal_tally2.filters = [azimuthal_filter] azimuthal_tally2.scores = ['flux'] azimuthal_tally2.estimator = 'analog' - azimuthal_filter2 = AzimuthalFilter(5) - azimuthal_tally3 = Tally() - azimuthal_tally3.filters = [azimuthal_filter2] - azimuthal_tally3.scores = ['flux'] - azimuthal_tally3.estimator = 'tracklength' - mesh_2x2 = Mesh(mesh_id=1) mesh_2x2.lower_left = [-182.07, -182.07] mesh_2x2.upper_right = [182.07, 182.07] mesh_2x2.dimension = [2, 2] mesh_filter = MeshFilter(mesh_2x2) - azimuthal_tally4 = Tally() - azimuthal_tally4.filters = [azimuthal_filter2, mesh_filter] - azimuthal_tally4.scores = ['flux'] - azimuthal_tally4.estimator = 'tracklength' + azimuthal_tally3 = Tally() + azimuthal_tally3.filters = [azimuthal_filter, mesh_filter] + azimuthal_tally3.scores = ['flux'] + azimuthal_tally3.estimator = 'tracklength' cellborn_tally = Tally() cellborn_tally.filters = [CellbornFilter((10, 21, 22, 23))] @@ -76,20 +72,17 @@ class TalliesTestHarness(PyAPITestHarness): material_tally.filters = [MaterialFilter((1, 2, 3, 4))] material_tally.scores = ['total'] + mu_bins = (-1.0, -0.5, 0.0, 0.5, 1.0) + mu_filter = MuFilter(mu_bins) mu_tally1 = Tally() - mu_tally1.filters = [MuFilter((-1.0, -0.5, 0.0, 0.5, 1.0))] + mu_tally1.filters = [mu_filter] mu_tally1.scores = ['scatter', 'nu-scatter'] - mu_filter = MuFilter(5) mu_tally2 = Tally() - mu_tally2.filters = [mu_filter] + mu_tally2.filters = [mu_filter, mesh_filter] mu_tally2.scores = ['scatter', 'nu-scatter'] - mu_tally3 = Tally() - mu_tally3.filters = [mu_filter, mesh_filter] - mu_tally3.scores = ['scatter', 'nu-scatter'] - - polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.1416) + polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.14159) polar_filter = PolarFilter(polar_bins) polar_tally1 = Tally() polar_tally1.filters = [polar_filter] @@ -101,17 +94,11 @@ class TalliesTestHarness(PyAPITestHarness): polar_tally2.scores = ['flux'] polar_tally2.estimator = 'analog' - polar_filter2 = PolarFilter((5,)) polar_tally3 = Tally() - polar_tally3.filters = [polar_filter2] + polar_tally3.filters = [polar_filter, mesh_filter] polar_tally3.scores = ['flux'] polar_tally3.estimator = 'tracklength' - polar_tally4 = Tally() - polar_tally4.filters = [polar_filter2, mesh_filter] - polar_tally4.scores = ['flux'] - polar_tally4.estimator = 'tracklength' - universe_tally = Tally() universe_tally.filters = [UniverseFilter((1, 2, 3, 4, 6, 8))] universe_tally.scores = ['total'] @@ -173,10 +160,9 @@ class TalliesTestHarness(PyAPITestHarness): self._input_set.tallies = Tallies() self._input_set.tallies += ( [azimuthal_tally1, azimuthal_tally2, azimuthal_tally3, - azimuthal_tally4, cellborn_tally, dg_tally, energy_tally, - energyout_tally, transfer_tally, material_tally, mu_tally1, - mu_tally2, mu_tally3, polar_tally1, polar_tally2, polar_tally3, - polar_tally4, universe_tally]) + cellborn_tally, dg_tally, energy_tally, energyout_tally, + transfer_tally, material_tally, mu_tally1, mu_tally2, + polar_tally1, polar_tally2, polar_tally3, universe_tally]) self._input_set.tallies += score_tallies self._input_set.tallies += flux_tallies self._input_set.tallies += (scatter_tally1, scatter_tally2)