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 8988cb9c0..857826e09 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -373,6 +373,8 @@ Core Classes openmc.data.CoherentElastic openmc.data.FissionEnergyRelease openmc.data.DataLibrary + openmc.data.Decay + openmc.data.FissionProductYields openmc.data.WindowedMultipole 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/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/data/__init__.py b/openmc/data/__init__.py index 248569487..635c317bf 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -9,6 +9,7 @@ WMP_VERSION = b'v0.2' from .data import * from .neutron import * +from .decay import * from .reaction import * from .ace import * from .angle_distribution import * 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/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/filter.py b/openmc/filter.py index 88a88f99d..77e409c40 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/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..ce7511193 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): @@ -1526,9 +1532,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) @@ -3207,8 +3223,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 +3242,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 +3282,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 +3360,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 +3370,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 +3412,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 +3428,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 +3455,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 +3531,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 +3587,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 +3596,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 +3628,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 @@ -3637,7 +3718,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 +3768,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 +3791,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 +3857,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 +3942,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 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/tallies.py b/openmc/tallies.py index d5d6df03f..979929a26 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1320,7 +1320,8 @@ 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])) 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_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_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)