diff --git a/.gitignore b/.gitignore index 023c4f2c1a..9f720f9894 100644 --- a/.gitignore +++ b/.gitignore @@ -19,6 +19,7 @@ src/openmc # Documentation builds docs/build +docs/source/_images/*.pdf # xml-fortran reader src/xml-fortran/xmlreader diff --git a/LICENSE b/LICENSE index f8d9cd7c04..0b9c351371 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2011-2013 Massachusetts Institute of Technology +Copyright (c) 2011-2014 Massachusetts Institute of Technology Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/docs/Makefile b/docs/Makefile index 8909256f3d..89c71dc074 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -6,13 +6,19 @@ SPHINXOPTS = SPHINXBUILD = sphinx-build PAPER = BUILDDIR = build +IMAGEDIR = source/_images # Internal variables. PAPEROPT_a4 = -D latex_paper_size=a4 PAPEROPT_letter = -D latex_paper_size=letter ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source -.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest +# SVG to PDF conversion +SVG2PDF = inkscape +PDFS = $(patsubst %.svg,%.pdf,$(wildcard $(IMAGEDIR)/*.svg)) + + +.PHONY: help images clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest help: @echo "Please use \`make ' where is one of" @@ -33,8 +39,16 @@ help: @echo " linkcheck to check all external links for integrity" @echo " doctest to run all doctests embedded in the documentation (if enabled)" +# Pattern rule for converting SVG to PDF +%.pdf: %.svg + $(SVG2PDF) -f $< -A $@ + +# Rule to build PDFs +images: $(PDFS) + clean: -rm -rf $(BUILDDIR)/* + -rm $(PDFS) html: $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html @@ -91,14 +105,14 @@ epub: @echo @echo "Build finished. The epub file is in $(BUILDDIR)/epub." -latex: +latex: images $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex @echo @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex." @echo "Run \`make' in that directory to run these through (pdf)latex" \ "(use \`make latexpdf' here to do that automatically)." -latexpdf: +latexpdf: images $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex @echo "Running LaTeX files through pdflatex..." make -C $(BUILDDIR)/latex all-pdf diff --git a/docs/img/3dcore.png b/docs/source/_images/3dcore.png similarity index 100% rename from docs/img/3dcore.png rename to docs/source/_images/3dcore.png diff --git a/docs/img/3dgeomplot.png b/docs/source/_images/3dgeomplot.png similarity index 100% rename from docs/img/3dgeomplot.png rename to docs/source/_images/3dgeomplot.png diff --git a/docs/source/_images/Tracks.png b/docs/source/_images/Tracks.png new file mode 100644 index 0000000000..39c83cd595 Binary files /dev/null and b/docs/source/_images/Tracks.png differ diff --git a/docs/img/atr.png b/docs/source/_images/atr.png similarity index 100% rename from docs/img/atr.png rename to docs/source/_images/atr.png diff --git a/docs/img/fluxplot.png b/docs/source/_images/fluxplot.png similarity index 100% rename from docs/img/fluxplot.png rename to docs/source/_images/fluxplot.png diff --git a/docs/img/fork.png b/docs/source/_images/fork.png similarity index 100% rename from docs/img/fork.png rename to docs/source/_images/fork.png diff --git a/docs/img/halfspace.svg b/docs/source/_images/halfspace.svg similarity index 100% rename from docs/img/halfspace.svg rename to docs/source/_images/halfspace.svg diff --git a/docs/img/master-slave.png b/docs/source/_images/master-slave.png similarity index 100% rename from docs/img/master-slave.png rename to docs/source/_images/master-slave.png diff --git a/docs/img/nearest-neighbor-example.png b/docs/source/_images/nearest-neighbor-example.png similarity index 100% rename from docs/img/nearest-neighbor-example.png rename to docs/source/_images/nearest-neighbor-example.png diff --git a/docs/img/nearest-neighbor.png b/docs/source/_images/nearest-neighbor.png similarity index 100% rename from docs/img/nearest-neighbor.png rename to docs/source/_images/nearest-neighbor.png diff --git a/docs/img/openmc.png b/docs/source/_images/openmc.png similarity index 100% rename from docs/img/openmc.png rename to docs/source/_images/openmc.png diff --git a/docs/img/plotmeshtally.png b/docs/source/_images/plotmeshtally.png similarity index 100% rename from docs/img/plotmeshtally.png rename to docs/source/_images/plotmeshtally.png diff --git a/docs/img/pullrequest.png b/docs/source/_images/pullrequest.png similarity index 100% rename from docs/img/pullrequest.png rename to docs/source/_images/pullrequest.png diff --git a/docs/img/union.svg b/docs/source/_images/union.svg similarity index 100% rename from docs/img/union.svg rename to docs/source/_images/union.svg diff --git a/docs/img/uniongrid.svg b/docs/source/_images/uniongrid.svg similarity index 100% rename from docs/img/uniongrid.svg rename to docs/source/_images/uniongrid.svg diff --git a/docs/source/conf.py b/docs/source/conf.py index f3553c33e1..7fb6226318 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -39,7 +39,7 @@ master_doc = 'index' # General information about the project. project = u'OpenMC' -copyright = u'2011-2013, Massachusetts Institute of Technology' +copyright = u'2011-2014, Massachusetts Institute of Technology' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -48,7 +48,7 @@ copyright = u'2011-2013, Massachusetts Institute of Technology' # The short X.Y version. version = "0.5" # The full version, including alpha/beta/rc tags. -release = "0.5.2" +release = "0.5.3" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. @@ -121,7 +121,7 @@ html_title = "OpenMC Documentation" # The name of an image file (relative to this directory) to place at the top # of the sidebar. -html_logo = '../img/openmc.png' +html_logo = '_images/openmc.png' # The name of an image file (within the static path) to use as favicon of the # docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 @@ -188,6 +188,8 @@ latex_documents = [ u'Massachusetts Institute of Technology', 'manual'), ] +latex_elements = {'preamble': '\\usepackage{enumitem}\\setlistdepth{9}'} + # The name of an image file (relative to this directory) to place at the top of # the title page. #latex_logo = None diff --git a/docs/source/devguide/index.rst b/docs/source/devguide/index.rst index 1eb7436d03..1ceba324c1 100644 --- a/docs/source/devguide/index.rst +++ b/docs/source/devguide/index.rst @@ -15,6 +15,6 @@ as debugging. structures styleguide workflow - xml-fortran + xml-parsing statepoint voxel diff --git a/docs/source/devguide/workflow.rst b/docs/source/devguide/workflow.rst index 9562c9584e..f3793628af 100644 --- a/docs/source/devguide/workflow.rst +++ b/docs/source/devguide/workflow.rst @@ -60,7 +60,7 @@ features and bug fixes. The general steps for contributing are as follows: repository with the same name under your personal account. As such, you can commit to it as you please without disrupting other developers. - .. image:: ../../img/fork.png + .. image:: ../_images/fork.png 2. Clone your fork of OpenMC and create a branch that branches off of *develop*: @@ -77,7 +77,7 @@ features and bug fixes. The general steps for contributing are as follows: 4. Issue a pull request from GitHub and select the *develop* branch of mit-crpg/openmc as the target. - .. image:: ../../img/pullrequest.png + .. image:: ../_images/pullrequest.png At a minimum, you should describe what the changes you've made are and why you are making them. If the changes are related to an oustanding issue, make diff --git a/docs/source/devguide/xml-fortran.rst b/docs/source/devguide/xml-fortran.rst deleted file mode 100644 index 456cd662f1..0000000000 --- a/docs/source/devguide/xml-fortran.rst +++ /dev/null @@ -1,40 +0,0 @@ -.. _devguide_xml-fortran: - -========================= -xml-fortran Input Parsing -========================= - -OpenMC relies on the xml-fortran package for reading and intrepreting the XML -input files for geometry, materials, settings, tallies, etc. The use of an XML -format makes writing input files considerably more flexible than would otherwise -be possible. - -With the xml-fortran package, extending the user input files to include new tags -is fairly straightforward. A "template" file exists for each diferent type of -input file that tells xml-fortran what to expect in a file. These template files -can be found in the src/templates directory. The steps for modifying/adding -input are as follows: - -1. Add a ````` tag to the desired template file, -e.g. src/templates/geometry_t.xml. See the `xml-fortran documentation`_ for a -description of the acceptable fields. - -2. In the input_xml module, any input given in your new tag will be read -automatically through a call to, e.g. read_xml_file_geometry_t. Whatever -variable name you specified should have the data available. - -3. Add code in the appropriate subroutine to check the variable for any possible -errors. - -4. Add a variable in OpenMC to copy the temporary variable into if there are no -errors. - -A set of `RELAX NG`_ schemata exists that enables real-time validation of input -files when using the GNU Emacs text editor. You should also modify the RELAX NG -schema for the template you changed (e.g. src/templates/geometry.rnc) so that -those who use Emacs can confirm whether their input is valid before they -run. You will need to be familiar with RELAX NG `compact syntax`_. - -.. _xml-fortran documentation: http://xml-fortran.sourceforge.net/documentation.html -.. _RELAX NG: http://relaxng.org/ -.. _compact syntax: http://relaxng.org/compact-tutorial-20030326.html diff --git a/docs/source/devguide/xml-parsing.rst b/docs/source/devguide/xml-parsing.rst new file mode 100644 index 0000000000..8310d53e4e --- /dev/null +++ b/docs/source/devguide/xml-parsing.rst @@ -0,0 +1,38 @@ +.. _devguide_xml-parsing: + +================= +XML Input Parsing +================= + +OpenMC relies on the FoX_ Fortran XML library for reading and intrepreting the +XML input files for geometry, materials, settings, tallies, etc. The use of an +XML format makes writing input files considerably more flexible than would +otherwise be possible. + +With the FoX library, extending the user input files to include new tags is +fairly straightforward. The steps for modifying/adding input are as follows: + +1. Add appropriate calls to procedures from the `xml_interface module`_, such as +``check_for_node``, ``get_node_value``, and ``get_node_array``. All input +reading is performed in the `input_xml module`_. + +2. Make sure that your input can be categorized as one of the datatypes from +`XML Schema Part 2`_ and that parsing of the data appropriately reflects +this. For example, for a boolean_ value, true can be represented either by "true" +or by "1". + +3. Add code to check the variable for any possible errors. + +A set of `RELAX NG`_ schemata exists that enables real-time validation of input +files when using the GNU Emacs text editor. You should also modify the RELAX NG +schema for the file you changed (e.g. src/relaxng/geometry.rnc) so that +those who use Emacs can confirm whether their input is valid before they +run. You will need to be familiar with RELAX NG `compact syntax`_. + +.. _FoX: https://github.com/andreww/fox +.. _xml_interface module: https://github.com/mit-crpg/openmc/blob/develop/src/xml_interface.F90 +.. _input_xml module: https://github.com/mit-crpg/openmc/blob/develop/src/input_xml.F90 +.. _XML Schema Part 2: http://www.w3.org/TR/xmlschema-2/ +.. _boolean: http://www.w3.org/TR/xmlschema-2/#boolean +.. _RELAX NG: http://relaxng.org/ +.. _compact syntax: http://relaxng.org/compact-tutorial-20030326.html diff --git a/docs/source/license.rst b/docs/source/license.rst index a7d7f29760..e7f4b3a69d 100644 --- a/docs/source/license.rst +++ b/docs/source/license.rst @@ -4,7 +4,7 @@ License Agreement ================= -Copyright © 2011-2013 Massachusetts Institute of Technology +Copyright © 2011-2014 Massachusetts Institute of Technology Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index c564896c79..a193dde96d 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -80,7 +80,7 @@ dashed box would need to be stored on a per-nuclide basis, and the union grid would need to be stored once. This method is also referred to as *double indexing* and is available as an option in Serpent (see paper by Leppanen_). -.. figure:: ../../img/uniongrid.svg +.. figure:: ../_images/uniongrid.* :width: 600px :align: center :figclass: align-center diff --git a/docs/source/methods/eigenvalue.rst b/docs/source/methods/eigenvalue.rst index a5ba0bf459..fe99ba22ec 100644 --- a/docs/source/methods/eigenvalue.rst +++ b/docs/source/methods/eigenvalue.rst @@ -108,6 +108,40 @@ at plots of :math:`k_{eff}` and the Shannon entropy. A number of methods have been proposed (see e.g. [Romano]_, [Ueki]_), but each of these is not without problems. +--------------------------- +Uniform Fission Site Method +--------------------------- + +Generally speaking, the variance of a Monte Carlo tally will be inversely +proportional to the number of events that score to the tally. In a reactor +problem, this implies that regions with low relative power density will have +higher variance that regions with high relative power density. One method to +circumvent the uneven distribution of relative errors is the uniform fission +site (UFS) method introduced by [Sutton]_. In this method, the portion of the +problem containing fissionable material is subdivided into a number of cells +(typically using a structured mesh). Rather than producing + +.. math:: + + m = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t} + +fission sites at each collision where :math:`w` is the weight of the neutron, +:math:`k` is the previous-generation estimate of the neutron multiplication +factor, :math:`\nu\Sigma_f` is the neutron production cross section, and +:math:`\Sigma_t` is the total cross section, in the UFS method we produce + +.. math:: + + m_{UFS} = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t} \frac{v_i}{s_i} + +fission sites at each collision where :math:`v_i` is the fraction of the total +volume occupied by cell :math:`i` and :math:`s_i` is the fraction of the fission +source contained in cell :math:`i`. To ensure that no bias is introduced, the +weight of each fission site stored in the fission bank is :math:`s_i/v_i` rather +than unity. By ensuring that the expected number of fission sites in each mesh +cell is constant, the collision density across all cells, and hence the variance +of tallies, is more uniform than it would be otherwise. + .. _Shannon entropy: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-06-3737_entropy.pdf .. [Lieberoth] J. Lieberoth, "A Monte Carlo Technique to Solve the Static @@ -119,5 +153,9 @@ problems. *Proc. International Conference on Mathematics, Computational Methods, and Reactor Physics*, Saratoga Springs, New York (2009). +.. [Sutton] Daniel J. Kelly, Thomas M. Sutton, and Stephen C. Wilson, "MC21 + Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark + Problem," *Proc. PHYSOR 2012*, Knoxville, Tennessee, Apr. 15--20 (2012). + .. [Ueki] Taro Ueki, "On-the-Fly Judgments of Monte Carlo Fission Source Convergence," *Trans. Am. Nucl. Soc.*, **98**, 512 (2008). diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index d4824f5d1c..1515ffa891 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -45,7 +45,7 @@ surface by a combination of the unique ID of the surface and a positive/negative sign. The following illustration shows an example of an ellipse with unique ID 1 dividing space into two half-spaces. -.. figure:: ../../img/halfspace.svg +.. figure:: ../_images/halfspace.* :align: center :figclass: align-center @@ -60,7 +60,7 @@ half-space references whose intersection defines the region. The region is then assigned a material defined elsewhere. The following illustration shows an example of a cell defined as the intersection of an ellipse and two planes. -.. figure:: ../../img/union.svg +.. figure:: ../_images/union.* :align: center :figclass: align-center diff --git a/docs/source/methods/parallelization.rst b/docs/source/methods/parallelization.rst index 7f0b9bd6eb..66ead52a2a 100644 --- a/docs/source/methods/parallelization.rst +++ b/docs/source/methods/parallelization.rst @@ -65,7 +65,7 @@ in the case of an eigenvalue calculation). This idea is illustrated in .. _figure-master-slave: -.. figure:: ../../img/master-slave.png +.. figure:: ../_images/master-slave.png :align: center :figclass: align-center @@ -122,7 +122,7 @@ needed. This concept is illustrated in :ref:`Figure 2 .. _figure-nearest-neighbor: -.. figure:: ../../img/nearest-neighbor.png +.. figure:: ../_images/nearest-neighbor.png :align: center :figclass: align-center @@ -203,7 +203,7 @@ communicated between adjacent nodes. .. _figure-neighbor-example: -.. figure:: ../../img/nearest-neighbor-example.png +.. figure:: ../_images/nearest-neighbor-example.png :align: center :figclass: align-center diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index c79144dfe3..54dc913455 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -790,6 +790,7 @@ outgoing angle is \mu = \frac{1}{A} \ln \left ( \xi_4 e^A + (1 - \xi_4) e^{-A} \right ). +.. _ace-law-61: ACE Law 61 - Correlated Energy and Angle Distribution +++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -952,7 +953,7 @@ as v_n \bar{\sigma} (v_n, T) = \int d\mathbf{v}_T v_r \sigma(v_r) M (\mathbf{v}_T) - + where :math:`v_n` is the magnitude of the velocity of the neutron, :math:`\bar{\sigma}` is an effective cross section, :math:`T` is the temperature of the target material, :math:`\mathbf{v}_T` is the velocity of the target @@ -1321,7 +1322,7 @@ given analytically by \mu = 1 - \frac{E_i}{E} -where :math:`E_i` is the energy of the Bragg edge that scattered the neutron. +where :math:`E_i` is the energy of the Bragg edge that scattered the neutron. Outgoing Angle for Incoherent Elastic Scattering ------------------------------------------------ @@ -1348,18 +1349,24 @@ where the interpolation factor is defined as Outgoing Energy and Angle for Inelastic Scattering -------------------------------------------------- -On each |sab| table, there is a correlated angle-energy secondary distribution -for neutron thermal inelastic scattering. While the documentation for the ACE -format implies that there are a series of equiprobable outgoing energies, the -outgoing energies may have non-uniform probability distribution. In particular, -if the thermal data were processed with :math:`iwt = 0` in NJOY, then the first -and last outgoing energies have a relative probability of 1, the second and -second to last energies have a relative probability of 4, and all other energies -have a relative probability of 10. The procedure to determine the outgoing -energy and angle is as such. First, the interpolation factor is determined from -equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is sampled -either from a uniform distribution or from the aforementioned skewed -distribution. The outgoing energy is then interpolated between values +Each |sab| table provides a correlated angle-energy secondary distribution for +neutron thermal inelastic scattering. There are three representations used +in the ACE thermal scattering data: equiprobable discrete outgoing +energies, non-uniform yet still discrete outgoing energies, and continuous +outgoing energies with corresponding probability and cumulative distribution +functions provided in tabular format. These three representations all +represent the angular distribution in a common format, using a series of +discrete equiprobable outgoing cosines. + +Equi-Probable Outgoing Energies ++++++++++++++++++++++++++++++++ + +If the thermal data was processed with :math:`iwt = 1` in NJOY, then the +outgoing energy spectra is represented in the ACE data as a set of discrete and +equiprobable outgoing energies. The procedure to determine the outgoing energy +and angle is as such. First, the interpolation factor is determined from +equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is +sampled from a uniform distribution and then interpolated between values corresponding to neighboring incoming energies: .. math:: @@ -1380,6 +1387,37 @@ uniformly and then the final cosine is interpolated on the incoming energy grid: where :math:`\mu_{i,j,k}` is the k-th outgoing cosine corresponding to the j-th outgoing energy and the i-th incoming energy. +Skewed Equi-Probable Outgoing Energies +++++++++++++++++++++++++++++++++++++++ + +If the thermal data was processed with :math:`iwt=0` in NJOY, then the +outgoing energy spectra is represented in the ACE data according to the +following: the first and last outgoing energies have a relative probability of +1, the second and second-to-last energies have a relative probability of 4, and +all other energies have a relative probability of 10. The procedure to +determine the outgoing energy and angle is similar to the method discussed +above, except that the sampled probability distribution is now skewed +accordingly. + +Continuous Outgoing Energies +++++++++++++++++++++++++++++ + +If the thermal data was processed with :math:`iwt=2` in NJOY, then the +outgoing energy spectra is represented by a continuous outgoing energy spectra +in tabular form with linear-linear interpolation. The sampling of the outgoing +energy portion of this format is very similar to :ref:`ACE Law 61`, +but the sampling of the correlated angle is performed as it was in the other +two representations discussed in this sub-section. In the Law 61 algorithm, +we found an interpolation factor :math:`f`, statistically sampled an incoming +energy bin :math:`\ell`, and sampled an outgoing energy bin :math:`j` based on +the tabulated cumulative distribution function. Once the outgoing energy has +been determined with equation :eq:`ace-law-4-energy`, we then need to decide +which angular distribution data to use. Like the linear-linear interpolation +case in Law 61, the angular distribution closest to the sampled value of the +cumulative distribution function for the outgoing energy is utilized. The +actual algorithm utilized to sample the outgoing angle is shown in equation +:eq:`inelastic-angle`. + .. _probability_tables: ---------------------------------------------- diff --git a/docs/source/publications.rst b/docs/source/publications.rst index d757b9c3ce..df7ca639bc 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -4,6 +4,53 @@ Publications ============ +- Benoit Forget, Sheng Xu, and Kord Smith, "Direct Doppler broadening in Monte + Carlo simulations using the multipole representation," *Ann. Nucl. Energy*, + **64**, 78--85 (2014). ``_ + +- Andrew Siegel, Kord Smith, Kyle Felker, Paul Romano, Benoit Forget, and Peter + Beckman, "Improved cache performance in Monte Carlo transport calculations + using energy banding," *Comput. Phys. Commun.* + (2013). ``_ + +- Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Validation of OpenMC + Reactor Physics Simulations with the B&W 1810 Series Benchmarks," + *Trans. Am. Nucl. Soc.*, **109**, 1301--1304 (2013). + +- Bryan R. Herman, Benoit Forget, and Kord Smith, "Utilizing CMFD in OpenMC to + Estimate Dominance Ratio and Adjoint," *Trans. Am. Nucl. Soc.*, **109**, + 1389-1392 (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). + +- Paul K. Romano, Benoit Forget, Kord Smith, and Andrew Siegel, "On the use of + tally servers in Monte Carlo simulations of light-water reactors," + *Proc. Joint International Conference on Supercomputing in Nuclear + Applications and Monte Carlo*, Paris, France, Oct. 27--31 (2013). + +- Paul K. Romano, Nicholas E. Horelik, Bryan R. Herman, Adam G. Nelson, Benoit + Forget, and Kord Smith, "OpenMC: A State-of-the-Art Monte Carlo Code for + Research and Development," *Proc. Joint International Conference on + Supercomputing in Nuclear Applications and Monte Carlo*, Paris, France, + Oct. 27--31 (2013). + +- Kyle G. Felker, Andrew R. Siegel, Kord S. Smith, Paul K. Romano, and Benoit + Forget, "The energy band memory server algorithm for parallel Monte Carlo + calculations," *Proc. Joint International Conference on Supercomputing in + Nuclear Applications and Monte Carlo*, Paris, France, Oct. 27--31 (2013). + +- John R. Tramm and Andrew R. Siegel, "Memory Bottlenecks and Memory Contention + in Multi-Core Monte Carlo Transport Codes," *Proc. Joint International + Conference on Supercomputing in Nuclear Applications and Monte Carlo*, Paris, + France, Oct. 27--31 (2013). + +- Andrew R. Siegel, Kord Smith, Paul K. Romano, Benoit Forget, and Kyle Felker, + "Multi-core performance studies of a Monte Carlo neutron transport code," + *Int. J. High Perform. Comput. Appl.*, **28** (1), 87--96 + (2014). ``_ + - Paul K. Romano, Andrew R. Siegel, Benoit Forget, and Kord Smith, "Data decomposition of Monte Carlo particle transport simulations via tally servers," *J. Comput. Phys.*, **252**, 20--36 diff --git a/docs/source/releasenotes/index.rst b/docs/source/releasenotes/index.rst index 876d953ff8..9164fa5943 100644 --- a/docs/source/releasenotes/index.rst +++ b/docs/source/releasenotes/index.rst @@ -10,6 +10,7 @@ bugs fixed, and known issues for each successive release. .. toctree:: :maxdepth: 1 + notes_0.5.3 notes_0.5.2 notes_0.5.1 notes_0.5.0 diff --git a/docs/source/releasenotes/notes_0.5.3.rst b/docs/source/releasenotes/notes_0.5.3.rst index fb7b1aa79a..6d93f9aece 100644 --- a/docs/source/releasenotes/notes_0.5.3.rst +++ b/docs/source/releasenotes/notes_0.5.3.rst @@ -4,10 +4,6 @@ Release Notes for OpenMC 0.5.3 ============================== -.. note:: - These release notes are for an upcoming release of OpenMC and are still - subject to change. - ------------------- System Requirements ------------------- @@ -21,12 +17,20 @@ the problem at hand (mostly on the number of nuclides in the problem). New Features ------------ -- Special run mode --tallies removed. +- Output interface enhanced to allow multiple files handles to be opened +- Particle restart file linked to output interface - Particle restarts and state point restarts are both identified with the -r command line flag. -- New regression test suite. -- All memory leaks fixed. -- Shared-memory parallelism with OpenMP. +- Particle instance no longer global, passed to all physics routines +- Physics routines refactored to rely less on global memory, more arguments + passed in +- CMFD routines refactored and now can compute dominance ratio on the fly +- PETSc 3.4.2 or higher must be used and compiled with fortran datatype support +- Memory leaks fixed except for ones from xml-fortran package +- Test suite enhanced to test output with different compiler options +- Description of OpenMC development workflow added +- OpenMP shared-memory parallelism added +- Special run mode --tallies removed. --------- Bug Fixes @@ -35,7 +39,7 @@ Bug Fixes - 2b1e8a_: Normalize direction vector after reflecting particle. - 5853d2_: Set blank default for cross section listing alias. - e178c7_: Fix infinite loop with words greater than 80 characters in write_message. -- c18a6e_: Chcek for valid secondary mode on S(a,b) tables. +- c18a6e_: Check for valid secondary mode on S(a,b) tables. - 82c456_: Fix bug where last process could have zero particles. .. _2b1e8a: https://github.com/mit-crpg/openmc/commit/2b1e8a diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 7a61ba2b5d..de6d2577a0 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -405,6 +405,15 @@ integers: the batch number, generation number, and particle number. *Default*: None +.. _track: + +```` Element +------------------- + +The ```` element specifies particles for which OpenMC will output binary files describing particle position at every step of its transport. This element should be followed by triplets of integers. Each triplet describes one particle. The integers in each triplet specify the batch number, generation number, and particle number, respectively. + + *Default*: None + ```` Element ------------------------ @@ -643,9 +652,10 @@ Each ```` element can have the following attributes or sub-elements: --------------------- The ```` can be used to represent repeating structures (e.g. fuel pins -in an assembly) or other geometry which naturally fits into a two-dimensional -structured mesh. Each cell within the lattice is filled with a specified -universe. A ```` accepts the following attributes or sub-elements: +in an assembly) or other geometry which naturally fits into a two- or +three-dimensional structured mesh. Each cell within the lattice is filled with a +specified universe. A ```` accepts the following attributes or +sub-elements: :id: A unique integer that can be used to identify the surface. @@ -657,18 +667,19 @@ universe. A ```` accepts the following attributes or sub-elements: *Default*: rectangular :dimension: - Two integers representing the number of lattice cells in the x- and y- - directions, respectively. + Two or three integers representing the number of lattice cells in the x- and + y- (and z-) directions, respectively. *Default*: None :lower_left: - The coordinates of the lower-left corner of the lattice. + The coordinates of the lower-left corner of the lattice. If the lattice is + two-dimensional, only the x- and y-coordinates are specified. *Default*: None :width: - The width of the lattice cell in the x- and y- directions. + The width of the lattice cell in the x- and y- (and z-) directions. *Default*: None @@ -1203,13 +1214,9 @@ with "false". ```` Element ------------------ -If a structured mesh is desired as a filter for a tally, it must be specified in -a separate element with the tag name ````. This element has the following +The CMFD mesh is a structured Cartesian mesh. This element has the following attributes/sub-elements: - :type: - The type of structured mesh. Only "rectangular" is currently supported. - :lower_left: The lower-left corner of the structured mesh. If only two coordinate are given, it is assumed that the mesh is an x-y mesh. diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 84f7633c2f..164b0151be 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -78,7 +78,7 @@ Prerequisites ./configure --prefix=/opt/hdf5/1.8.11-gnu --enable-fortran \ --enable-fortran2003 --enable-parallel - You may omit '--enable-parallel' if you want to compile HDF5_ in serial. + You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial. * PETSc_ for CMFD acceleration @@ -93,7 +93,7 @@ Prerequisites --with-fortran-datatypes The BLAS/LAPACK library is not required to be downloaded and can be linked - explicitly (e.g., Intel MLK library). + explicitly (e.g., Intel MKL library). * git_ version control software for obtaining source code @@ -358,6 +358,7 @@ OpenMC accepts the following command line flags: -r, --restart file Restart a previous run from a state point or a particle restart file -s, --threads N Run with *N* OpenMP threads +-t, --track Write tracks for all particles -v, --version Show version information ----------------------------------------------------- diff --git a/docs/source/usersguide/processing.rst b/docs/source/usersguide/processing.rst index 7105e9f7d1..fc426ba2f4 100644 --- a/docs/source/usersguide/processing.rst +++ b/docs/source/usersguide/processing.rst @@ -38,7 +38,7 @@ running OpenMC with the -plot or -p command-line option (See Plotting in 2D -------------- -.. image:: ../../img/atr.png +.. image:: ../_images/atr.png :height: 200px After running OpenMC to obtain PPM files, images should be saved to another @@ -58,7 +58,7 @@ Ubuntu: ``sudo apt-get install imagemagick``). Images are then converted like: Plotting in 3D -------------- -.. image:: ../../img/3dgeomplot.png +.. image:: ../_images/3dgeomplot.png :height: 200px The binary VOXEL files output by OpenMC can not be viewed directly by any @@ -162,7 +162,7 @@ tasks will be described here in the following sections. Plotting in 2D -------------- -.. image:: ../../img/plotmeshtally.png +.. image:: ../_images/plotmeshtally.png :height: 200px For simple viewing of 2D slices of a mesh plot, the utility plot_mesh_tally.py @@ -170,7 +170,7 @@ is provided. This utility provides an interactive GUI to explore and plot mesh tallies for any scores and filter bins. It requires statepoint.py, as well as `PyQt `_. -.. image:: ../../img/fluxplot.png +.. image:: ../_images/fluxplot.png :height: 200px Alternatively, the user can write their own Python script to manipulate the data @@ -249,7 +249,7 @@ two heatmaps in the previous figure. Plotting in 3D -------------- -.. image:: ../../img/3dcore.png +.. image:: ../_images/3dcore.png :height: 200px As with 3D plots of the geometry, meshtally data needs to be put into a standard @@ -353,9 +353,89 @@ file. Note that the data contained in the output from ``StatePoint.extract_result`` is already in a Numpy array that can be reshaped and dumped to MATLAB in one step. +---------------------------- +Particle Track Visualization +---------------------------- +.. image:: ../_images/Tracks.png + :height: 200px +OpenMC can dump particle tracks—the position of particles as they are +transported through the geometry. There are two ways to make OpenMC output +tracks: all particle tracks through a commandline argument or specific particle +tracks through settings.xml. +Running OpenMC with the argument "-t", "-track", or "--track" will cause a track +file to be created for every particle transported in the code. +The settings.xml file can dictate that specific particle tracks are output. +These particles are specified withen a ''track'' element. The ''track'' element +should contain triplets of integers specifying the batch, generation, and +particle numbers, respectively. For example, to output the tracks for particles +3 and 4 of batch 1 and generation 2 the settings.xml file should contain: +.. code-block:: xml + + 1 2 3 + 1 2 4 + + +After running OpenMC, the directory should contain a file of the form +"track_(batch #)_(generation #)_(particle #).(binary or h5)" for each particle +tracked. These track files can be converted into VTK poly data files with the +"track.py" utility. The usage of track.py is of the form "track.py [-o OUT] IN" +where OUT is the optional output filename and IN is one or more filenames +describing track files. The default output name is "track.pvtp". A common +usage of track.py is "track.py track*.binary" which will use the data from all +binary track files in the directory to write a "track.pvtp" VTK output file. +The .pvtp file can then be read and plotted by 3d visualization programs such as +Paraview. + +---------------------- +Source Site Processing +---------------------- + +For eigenvalue problems, OpenMC will store information on the fission source +sites in the statepoint file by default. For each source site, the weight, +position, sampled direction, and sampled energy are stored. To extract this data +from a statepoint file, the statepoint.py Python module can be used. Below is an +example of an interactive ipython session using the statepoint.py Python module: + +.. code-block:: python + + In [1]: import statepoint + + In [2]: sp = statepoint.StatePoint('statepoint.100.h5') + + In [3]: sp.read_source() + + In [4]: len(sp.source) + Out[4]: 1000 + + In [5]: sp.source[0:10] + Out[5]: + [, + , + , + , + , + , + , + , + , + ] + + In [6]: site = sp.source[0] + + In [7]: site.weight + Out[7]: 1.0 + + In [8]: site.xyz + Out[8]: array([ 2.21980946, -8.92686048, 87.93720485]) + + In [9]: site.uvw + Out[9]: array([ 0.06740523, 0.50612814, 0.85982024]) + + In [10]: site.E + Out[10]: 0.93292326356564159 diff --git a/docs/source/usersguide/troubleshoot.rst b/docs/source/usersguide/troubleshoot.rst index 49cf41c8de..a13272aa65 100644 --- a/docs/source/usersguide/troubleshoot.rst +++ b/docs/source/usersguide/troubleshoot.rst @@ -19,25 +19,6 @@ you are using a compiler that does not support type-bound procedures from Fortran 2003. This affects any version of gfortran prior to 4.6. Downloading and installing the latest gfortran_ compiler should resolve this problem. -Fatal Error: Wrong module version '4' (expected '9') for file 'xml_data_cmfd_t.mod' opened at (1) -************************************************************************************************* - -The `.mod` modules files that are created by gfortran are versioned and -sometimes are usually not backwards compatible. If gfortran is upgraded and the -modules files for xml-fortran source files are not deleted, this error may -occur. To fix this, clear out all module and object files with :program:`make -distclean` and then recompiling. - -Fatal Error: File 'xml_data_cmfd_t.mod' opened at (1) is not a GFORTRAN module file -*********************************************************************************** - -When OpenMC compiles, the first thing it needs to do is compile source in the -xml-fortran subdirectory. If you compiled everything with a compiler other than -gfortran, performed a :program:`make clean`, and then tried to :program:`make` -with gfortran, the xml-fortran modules would have been compiled with a different -compiler. To fix this, try clearing out all module and object files with -:program:`make distclean` and then recompiling. - gfortran: unrecognized option '-cpp' ************************************ diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index fb0fdd7177..db91d7a8f0 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -27,6 +27,9 @@ Restart a previous run from a state point or a particle restart file named .BI \-s " N" "\fR,\fP \-\-threads" " N" Use \fIN\fP OpenMP threads. .TP +.B "\-t\fR, \fP\-\-track" +Write tracks for all particles. +.TP .B "\-v\fR, \fP\-\-version" Show version information. .TP @@ -43,7 +46,7 @@ to locate ACE format cross section libraries if the user has not specified the tag in .I settings.xml\fP. .SH LICENSE -Copyright \(co 2011-2013 Massachusetts Institute of Technology. +Copyright \(co 2011-2014 Massachusetts Institute of Technology. .PP Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index b1c339921e..579bf26fb1 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -329,6 +329,7 @@ solver_interface.o: vector_header.o source.o: bank_header.o source.o: constants.o source.o: error.o +source.o: geometry.o source.o: geometry_header.o source.o: global.o source.o: math.o @@ -370,6 +371,11 @@ tally_initialize.o: tally_header.o timer_header.o: constants.o +track_output.o: global.o +track_output.o: output_interface.o +track_output.o: particle_header.o +track_output.o: string.o + tracking.o: cross_section.o tracking.o: error.o tracking.o: geometry.o @@ -381,10 +387,10 @@ tracking.o: physics.o tracking.o: random_lcg.o tracking.o: string.o tracking.o: tally.o +tracking.o: track_output.o vector_header.o: constants.o xml_interface.o: constants.o xml_interface.o: error.o xml_interface.o: global.o - diff --git a/src/Makefile b/src/Makefile index 03c370e07a..ad2c36730f 100644 --- a/src/Makefile +++ b/src/Makefile @@ -181,17 +181,17 @@ endif ifeq ($(OPENMP),yes) ifeq ($(COMPILER),intel) - F90FLAGS += -openmp -DOPENMP + F90FLAGS += -openmp LDFLAGS += -openmp endif ifeq ($(COMPILER),gnu) - F90FLAGS += -fopenmp -DOPENMP + F90FLAGS += -fopenmp LDFLAGS += -fopenmp endif ifeq ($(COMPILER),ibm) - F90FLAGS += -qsmp=omp -WF,-DOPENMP + F90FLAGS += -qsmp=omp LDFLAGS += -qsmp=omp endif endif @@ -218,7 +218,7 @@ endif ifeq ($(MACHINE),bluegene) F90 = /bgsys/drivers/ppcfloor/comm/xl/bin/mpixlf2003 - F90FLAGS = -WF,-DNO_F2008,-DMPI -O3 + F90FLAGS = -WF,-DNO_F2008,-DMPI,-DRESTRICTED_ASSOCIATED_BUG -O3 LDFLAGS = -lmpich.cnkf90 endif @@ -233,7 +233,7 @@ endif ifeq ($(MACHINE),bluegeneq) F90 = mpixlf2003 - F90FLAGS = -WF,-DNO_F2008,-DMPI -O5 + F90FLAGS = -WF,-DNO_F2008,-DMPI,-DRESTRICTED_ASSOCIATED_BUG -O5 endif #=============================================================================== diff --git a/src/ace.F90 b/src/ace.F90 index c37a546c01..aa6de5ab64 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -115,9 +115,9 @@ contains ! search through the list of nuclides for one which has a matching zaid sab => sab_tables(mat % i_sab_tables(k)) - ! Loop through nuclides and find match + ! Loop through nuclides and find match FIND_NUCLIDE: do j = 1, mat % n_nuclides - if (nuclides(mat % nuclide(j)) % zaid == sab % zaid) then + if (any(sab % zaid == nuclides(mat % nuclide(j)) % zaid)) then mat % i_sab_nuclides(k) = j exit FIND_NUCLIDE end if @@ -161,16 +161,16 @@ contains mat % i_sab_tables(m) = temp_table end do SORT_SAB end if - + ! Deallocate temporary arrays for names of nuclides and S(a,b) tables if (allocated(mat % names)) deallocate(mat % names) if (allocated(mat % sab_names)) deallocate(mat % sab_names) end do MATERIAL_LOOP2 - + ! Avoid some valgrind leak errors call already_read % clear() - + end subroutine read_xs !=============================================================================== @@ -244,8 +244,16 @@ contains ! Read first line of header read(UNIT=in, FMT='(A10,2G12.0,1X,A10)') name, awr, kT, date_ + ! Check that correct xs was found -- if cross_sections.xml is broken, the + ! location of the table may be wrong + if(adjustl(name) /= adjustl(listing % name)) then + message = "XS listing entry " // trim(listing % name) // " did not & + &match ACE data, " // trim(name) // " found instead." + call fatal_error() + end if + ! Read more header and NXS and JXS - read(UNIT=in, FMT=100) comment, mat, & + read(UNIT=in, FMT=100) comment, mat, & (zaids(i), awrs(i), i=1,16), NXS, JXS 100 format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/& ,8I9/8I9/8I9/8I9/8I9/8I9) @@ -269,7 +277,7 @@ contains ACCESS='direct', RECL=record_length) ! Read all header information - read(UNIT=in, REC=location) name, awr, kT, date_, & + read(UNIT=in, REC=location) name, awr, kT, date_, & comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS ! determine table length @@ -325,7 +333,15 @@ contains sab % name = name sab % awr = awr sab % kT = kT - sab % zaid = zaids(1) + ! Find sab % n_zaid + do i = 1, 16 + if (zaids(i) == 0) then + sab % n_zaid = i - 1 + exit + end if + end do + allocate(sab % zaid(sab % n_zaid)) + sab % zaid = zaids(1: sab % n_zaid) call read_thermal_data(sab) end select @@ -398,7 +414,7 @@ contains integer :: LNU ! type of nu data (polynomial or tabular) integer :: NC ! number of polynomial coefficients integer :: NR ! number of interpolation regions - integer :: NE ! number of energies + integer :: NE ! number of energies integer :: NPCR ! number of delayed neutron precursor groups integer :: LED ! location of energy distribution locators integer :: LDIS ! location of all energy distributions @@ -801,7 +817,7 @@ contains LED = JXS(10) - ! Loop over all reactions + ! Loop over all reactions do i = 1, NXS(5) rxn => nuc % reactions(i+1) ! skip over elastic scattering rxn % has_energy_dist = .true. @@ -832,7 +848,7 @@ contains integer :: LDIS ! location of all energy distributions integer :: LNW ! location of next energy distribution if multiple - integer :: LAW ! secondary energy distribution law + integer :: LAW ! secondary energy distribution law integer :: NR ! number of interpolation regions integer :: NE ! number of incoming energies integer :: IDAT ! location of first energy distribution for given MT @@ -1179,15 +1195,10 @@ contains integer :: NE_out ! number of outgoing energies integer :: NMU ! number of outgoing angles integer :: JXS4 ! location of elastic energy table + integer(8), allocatable :: LOCC(:) ! Location of inelastic data - ! read secondary energy mode for inelastic scattering and check + ! read secondary energy mode for inelastic scattering table % secondary_mode = NXS(7) - if (table % secondary_mode /= SAB_SECONDARY_EQUAL .and. & - table % secondary_mode /= SAB_SECONDARY_SKEWED) then - message = "Unsupported secondary mode on S(a,b) table " // & - trim(adjustl(table % name)) // ": " // to_str(table % secondary_mode) - call fatal_error() - end if ! read number of inelastic energies and allocate arrays NE_in = int(XSS(JXS(1))) @@ -1205,29 +1216,67 @@ contains ! allocate space for outgoing energy/angle for inelastic ! scattering - NE_out = NXS(4) - NMU = NXS(3) + 1 - table % n_inelastic_e_out = NE_out - table % n_inelastic_mu = NMU - allocate(table % inelastic_e_out(NE_out, NE_in)) - allocate(table % inelastic_mu(NMU, NE_out, NE_in)) + if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. & + table % secondary_mode == SAB_SECONDARY_SKEWED) then + NMU = NXS(3) + 1 + table % n_inelastic_mu = NMU + NE_out = NXS(4) + table % n_inelastic_e_out = NE_out + allocate(table % inelastic_e_out(NE_out, NE_in)) + allocate(table % inelastic_mu(NMU, NE_out, NE_in)) + else if (table % secondary_mode == SAB_SECONDARY_CONT) then + NMU = NXS(3) - 1 + table % n_inelastic_mu = NMU + allocate(table % inelastic_data(NE_in)) + allocate(LOCC(NE_in)) + ! NE_out will be determined later + end if ! read outgoing energy/angle distribution for inelastic scattering - lc = JXS(3) - 1 - do i = 1, NE_in - do j = 1, NE_out - ! read outgoing energy - table % inelastic_e_out(j,i) = XSS(lc + 1) + if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. & + table % secondary_mode == SAB_SECONDARY_SKEWED) then + lc = JXS(3) - 1 + do i = 1, NE_in + do j = 1, NE_out + ! read outgoing energy + table % inelastic_e_out(j,i) = XSS(lc + 1) - ! read outgoing angles for this outgoing energy - do k = 1, NMU - table % inelastic_mu(k,j,i) = XSS(lc + 1 + k) + ! read outgoing angles for this outgoing energy + do k = 1, NMU + table % inelastic_mu(k,j,i) = XSS(lc + 1 + k) + end do + + ! advance pointer + lc = lc + 1 + NMU end do - - ! advance pointer - lc = lc + 1 + NMU end do - end do + else if (table % secondary_mode == SAB_SECONDARY_CONT) then + ! Get the location pointers to each Ein's DistEnergySAB data + LOCC = get_int(NE_in) + ! Get the number of outgoing energies and allocate space accordingly + do i = 1, NE_in + NE_out = int(XSS(XSS_index + i - 1)) + table % inelastic_data(i) % n_e_out = NE_out + allocate(table % inelastic_data(i) % e_out (NE_out)) + allocate(table % inelastic_data(i) % e_out_pdf (NE_out)) + allocate(table % inelastic_data(i) % e_out_cdf (NE_out)) + allocate(table % inelastic_data(i) % mu (NMU, NE_out)) + end do + + ! Now we can fill the inelastic_data(i) attributes + do i = 1, NE_in + XSS_index = LOCC(i) + NE_out = table % inelastic_data(i) % n_e_out + do j = 1, NE_out + table % inelastic_data(i) % e_out(j) = XSS(XSS_index + 1) + table % inelastic_data(i) % e_out_pdf(j) = XSS(XSS_index + 2) + table % inelastic_data(i) % e_out_cdf(j) = XSS(XSS_index + 3) + table % inelastic_data(i) % mu(:, j) = & + XSS(XSS_index + 4: XSS_index + 4 + NMU - 1) + XSS_index = XSS_index + 4 + NMU - 1 + end do + end do + end if ! read number of elastic energies and allocate arrays JXS4 = JXS(4) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index d075c22889..e5bf1bb3a4 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -16,7 +16,7 @@ module ace_header integer, allocatable :: type(:) ! type of distribution integer, allocatable :: location(:) ! location of each table real(8), allocatable :: data(:) ! angular distribution data - + ! Type-Bound procedures contains procedure :: clear => distangle_clear ! Deallocates DistAngle @@ -32,10 +32,10 @@ module ace_header type(Tab1) :: p_valid ! probability of law validity real(8), allocatable :: data(:) ! energy distribution data - ! For reactions that may have multiple energy distributions such as (n.2n), + ! For reactions that may have multiple energy distributions such as (n,2n), ! this pointer allows multiple laws to be stored type(DistEnergy), pointer :: next => null() - + ! Type-Bound procedures contains procedure :: clear => distenergy_clear ! Deallocates DistEnergy @@ -57,7 +57,7 @@ module ace_header logical :: has_energy_dist ! Energy distribution present? type(DistAngle) :: adist ! Secondary angular distribution type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution - + ! Type-Bound procedures contains procedure :: clear => reaction_clear ! Deallocates Reaction @@ -76,7 +76,7 @@ module ace_header logical :: multiply_smooth ! multiply by smooth cross section? real(8), allocatable :: energy(:) ! incident energies real(8), allocatable :: prob(:,:,:) ! actual probabibility tables - + ! Type-Bound procedures contains procedure :: clear => urrdata_clear ! Deallocates UrrData @@ -137,22 +137,37 @@ module ace_header ! Reactions integer :: n_reaction ! # of reactions type(Reaction), pointer :: reactions(:) => null() - + ! Type-Bound procedures contains procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide +!=============================================================================== +! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic +! thermal scattering collisions which utilize a continuous secondary energy +! representation. +!=============================================================================== + + type DistEnergySab + integer :: n_e_out + real(8), allocatable :: e_out(:) + real(8), allocatable :: e_out_pdf(:) + real(8), allocatable :: e_out_cdf(:) + real(8), allocatable :: mu(:,:) + end type DistEnergySab + !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc !=============================================================================== - + type SAlphaBeta - character(10) :: name ! name of table, e.g. lwtr.10t - integer :: zaid ! Z and A identifier, e.g. 6012 for Carbon-12 - real(8) :: awr ! weight of nucleus in neutron masses - real(8) :: kT ! temperature in MeV (k*T) + character(10) :: name ! name of table, e.g. lwtr.10t + real(8) :: awr ! weight of nucleus in neutron masses + real(8) :: kT ! temperature in MeV (k*T) + integer :: n_zaid ! Number of valid zaids + integer, allocatable :: zaid(:) ! List of valid Z and A identifiers, e.g. 6012 ! threshold for S(a,b) treatment (usually ~4 eV) real(8) :: threshold_inelastic @@ -162,11 +177,17 @@ module ace_header integer :: n_inelastic_e_in ! # of incoming E for inelastic integer :: n_inelastic_e_out ! # of outgoing E for inelastic integer :: n_inelastic_mu ! # of outgoing angles for inelastic - integer :: secondary_mode ! secondary mode (equal/skewed) + integer :: secondary_mode ! secondary mode (equal/skewed/continuous) real(8), allocatable :: inelastic_e_in(:) - real(8), allocatable :: inelastic_sigma(:) + real(8), allocatable :: inelastic_sigma(:) + ! The following are used only if secondary_mode is 0 or 1 real(8), allocatable :: inelastic_e_out(:,:) real(8), allocatable :: inelastic_mu(:,:,:) + ! The following is used only if secondary_mode is 3 + ! The different implementation is necessary because the continuous + ! representation has a variable number of outgoing energy points for each + ! incoming energy + type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein ! Elastic scattering data integer :: elastic_mode ! elastic mode (discrete/exact) @@ -214,8 +235,9 @@ module ace_header real(8) :: kappa_fission ! microscopic energy-released from fission ! Information for S(a,b) use - integer :: index_sab ! index in sab_tables (zero means no table) - real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table + integer :: index_sab ! index in sab_tables (zero means no table) + integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide + real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table ! Information for URR probability table use logical :: use_ptable ! in URR range with probability tables? @@ -239,125 +261,125 @@ module ace_header !=============================================================================== ! DISTANGLE_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - +!=============================================================================== + subroutine distangle_clear(this) - + class(DistAngle), intent(inout) :: this ! The DistAngle object to clear - + if (allocated(this % energy)) & deallocate(this % energy, this % type, this % location, this % data) - - end subroutine distangle_clear + + end subroutine distangle_clear !=============================================================================== ! DISTENERGY_CLEAR resets and deallocates data in DistEnergy. -!=============================================================================== - +!=============================================================================== + recursive subroutine distenergy_clear(this) - + class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear - + ! Clear p_valid call this % p_valid % clear() - + if (allocated(this % data)) & deallocate(this % data) - + if (associated(this % next)) then ! recursively clear this item call this % next % clear() deallocate(this % next) end if - + end subroutine distenergy_clear - + !=============================================================================== ! REACTION_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - +!=============================================================================== + subroutine reaction_clear(this) - + class(Reaction), intent(inout) :: this ! The Reaction object to clear - + if (allocated(this % sigma)) & deallocate(this % sigma) - + if (associated(this % edist)) then call this % edist % clear() deallocate(this % edist) end if - + call this % adist % clear() - - end subroutine reaction_clear - + + end subroutine reaction_clear + !=============================================================================== ! URRDATA_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - +!=============================================================================== + subroutine urrdata_clear(this) - + class(UrrData), intent(inout) :: this ! The UrrData object to clear - + if (allocated(this % energy)) & deallocate(this % energy, this % prob) - - end subroutine urrdata_clear + + end subroutine urrdata_clear !=============================================================================== ! NUCLIDE_CLEAR resets and deallocates data in Nuclide. -!=============================================================================== - +!=============================================================================== + subroutine nuclide_clear(this) - + class(Nuclide), intent(inout) :: this ! The Nuclide object to clear - + integer :: i ! Loop counter - + if (allocated(this % grid_index)) & deallocate(this % grid_index) - + if (allocated(this % energy)) & deallocate(this % total, this % elastic, this % fission, & this % nu_fission, this % absorption) if (allocated(this % heating)) & deallocate(this % heating) - + if (allocated(this % index_fission)) & deallocate(this % index_fission) - + if (allocated(this % nu_t_data)) & deallocate(this % nu_t_data) - + if (allocated(this % nu_p_data)) & deallocate(this % nu_p_data) - + if (allocated(this % nu_d_data)) & deallocate(this % nu_d_data) - + if (allocated(this % nu_d_precursor_data)) & deallocate(this % nu_d_precursor_data) - + if (associated(this % nu_d_edist)) then do i = 1, size(this % nu_d_edist) call this % nu_d_edist(i) % clear() end do deallocate(this % nu_d_edist) end if - + if (associated(this % urr_data)) then call this % urr_data % clear() deallocate(this % urr_data) end if - + if (associated(this % reactions)) then do i = 1, size(this % reactions) call this % reactions(i) % clear() end do deallocate(this % reactions) end if - - end subroutine nuclide_clear + + end subroutine nuclide_clear end module ace_header diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index d0e80aaa90..8dd0fa35ff 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -934,7 +934,7 @@ contains subroutine compute_effective_downscatter() use constants, only: ZERO, CMFD_NOACCEL - use global, only: cmfd, cmfd_downscatter + use global, only: cmfd integer :: nx ! number of mesh cells in x direction integer :: ny ! number of mesh cells in y direction diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 7a35a12144..6234d9fcf9 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -117,9 +117,9 @@ contains subroutine process_cmfd_options() +#ifdef PETSC use global, only: cmfd_snes_monitor, cmfd_ksp_monitor, mpi_err -#ifdef PETSC ! Check for snes monitor if (cmfd_snes_monitor) call PetscOptionsSetValue("-snes_monitor", & "stdout", mpi_err) @@ -138,9 +138,10 @@ contains subroutine calc_fission_source() use constants, only: CMFD_NOACCEL, ZERO, TWO - use global, only: cmfd, cmfd_coremap, master, mpi_err, entropy_on, & - current_batch + use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch + #ifdef MPI + use global, only: mpi_err use mpi #endif @@ -261,13 +262,14 @@ contains use constants, only: ZERO, ONE use error, only: warning, fatal_error - use global, only: n_particles, meshes, source_bank, work, & - n_user_meshes, message, cmfd, master, mpi_err + use global, only: meshes, source_bank, work, n_user_meshes, message, & + cmfd, master use mesh_header, only: StructuredMesh use mesh, only: count_bank_sites, get_mesh_indices use search, only: binary_search #ifdef MPI + use global, only: mpi_err use mpi #endif diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index ac3d59824e..6f259f9497 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -20,7 +20,9 @@ contains use cmfd_header, only: allocate_cmfd +#ifdef PETSC integer :: new_comm ! new mpi communicator +#endif integer :: color ! color group of processor ! Read in cmfd input file @@ -34,17 +36,17 @@ contains end if ! Split up procs -# ifdef PETSC +#ifdef PETSC call MPI_COMM_SPLIT(MPI_COMM_WORLD, color, 0, new_comm, mpi_err) -# endif +#endif ! assign to PETSc -# ifdef PETSC +#ifdef PETSC PETSC_COMM_WORLD = new_comm ! Initialize PETSc on all procs call PetscInitialize(PETSC_NULL_CHARACTER, mpi_err) -# endif +#endif ! Initialize timers call time_cmfd % reset() @@ -575,7 +577,9 @@ contains end do ! Put cmfd tallies into active tally array and turn tallies on +!$omp parallel call setup_active_cmfdtallies() +!$omp end parallel tallies_on = .true. end subroutine create_cmfd_tally diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90 index 6e9917ac5b..69be017644 100644 --- a/src/cmfd_loss_operator.F90 +++ b/src/cmfd_loss_operator.F90 @@ -26,8 +26,10 @@ contains integer :: nnz ! number of nonzeros in matrix integer :: n_i ! number of interior cells integer :: n_c ! number of corner cells - integer :: n_s ! number side cells + integer :: n_s ! number of side cells + integer :: n_e ! number of edge cells integer :: nz_c ! number of non-zero corner cells + integer :: nz_e ! number of non-zero edge cells integer :: nz_s ! number of non-zero side cells integer :: nz_i ! number of non-zero interior cells @@ -48,13 +50,16 @@ contains if (cmfd_coremap) then nnz = preallocate_loss_matrix(nx, ny, nz, ng, n) else ! structured Cartesian grid - n_c = 4 ! define # of corners - n_s = 2*(nx + ny) - 8 ! define # of sides - n_i = nx*ny - (n_c + n_S) ! define # of interiors - nz_c = ng*n_c*(3 + ng - 1) ! define # nonzero corners - nz_s = ng*n_s*(4 + ng - 1) ! define # nonzero sides - nz_i = ng*n_i*(5 + ng - 1) ! define # nonzero interiors - nnz = nz_c + nz_s + nz_i + n_c = 8 ! define # of corners + n_e = 4*(nx - 2) + 4*(ny - 2) + 4*(nz - 2) ! define # of edges + n_s = 2*(nx - 2)*(ny - 2) + 2*(nx - 2)*(nz - 2) & + + 2*(ny - 2)*(nz - 2) ! define # of sides + n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors + nz_c = ng*n_c*(4 + ng - 1) ! define # nonzero corners + nz_e = ng*n_e*(5 + ng - 1) ! define # nonzero edges + nz_s = ng*n_s*(6 + ng - 1) ! define # nonzero sides + nz_i = ng*n_i*(7 + ng - 1) ! define # nonzero interiors + nnz = nz_c + nz_e + nz_s + nz_i end if ! Configure loss matrix diff --git a/src/cmfd_power_solver.F90 b/src/cmfd_power_solver.F90 index c1309b72cd..f2761d0a1f 100644 --- a/src/cmfd_power_solver.F90 +++ b/src/cmfd_power_solver.F90 @@ -102,7 +102,9 @@ contains subroutine init_data(adjoint) use constants, only: ONE, ZERO +#ifdef PETSC use global, only: cmfd_write_matrices +#endif logical, intent(in) :: adjoint ! adjoint calcualtion diff --git a/src/constants.F90 b/src/constants.F90 index 6997ec319c..05327153d0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -21,7 +21,7 @@ module constants FILETYPE_SOURCE = -3 ! ============================================================================ - ! ADJUSTABLE PARAMETERS + ! ADJUSTABLE PARAMETERS ! NOTE: This is the only section of the constants module that should ever be ! adjusted. Modifying constants in other sections may cause the code to fail. @@ -51,6 +51,10 @@ module constants integer, parameter :: MAX_WORD_LEN = 150 integer, parameter :: MAX_FILE_LEN = 255 + ! Maximum number of external source spatial resamples to encounter before an + ! error is thrown. + integer, parameter :: MAX_EXTSRC_RESAMPLES = 10000 + ! ============================================================================ ! PHYSICAL CONSTANTS @@ -111,9 +115,9 @@ module constants ! Surface types integer, parameter :: & - SURF_PX = 1, & ! Plane parallel to x-plane - SURF_PY = 2, & ! Plane parallel to y-plane - SURF_PZ = 3, & ! Plane parallel to z-plane + SURF_PX = 1, & ! Plane parallel to x-plane + SURF_PY = 2, & ! Plane parallel to y-plane + SURF_PZ = 3, & ! Plane parallel to z-plane SURF_PLANE = 4, & ! Arbitrary plane SURF_CYL_X = 5, & ! Cylinder along x-axis SURF_CYL_Y = 6, & ! Cylinder along y-axis @@ -144,7 +148,7 @@ module constants ELECTRON = 3 ! Angular distribution type - integer, parameter :: & + integer, parameter :: & ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution ANGLE_32_EQUI = 2, & ! 32 equiprobable bins ANGLE_TABULAR = 3 ! Tabular angular distribution @@ -152,7 +156,8 @@ module constants ! Secondary energy mode for S(a,b) inelastic scattering integer, parameter :: & SAB_SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins - SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins + SAB_SECONDARY_SKEWED = 1, & ! Skewed outgoing energy bins + SAB_SECONDARY_CONT = 2 ! Continuous, linear-linear interpolation ! Elastic mode for S(a,b) elastic scattering integer, parameter :: & @@ -276,7 +281,7 @@ module constants SCORE_KAPPA_FISSION = -12, & ! fission energy production rate SCORE_CURRENT = -13, & ! partial current SCORE_EVENTS = -14 ! number of events - + ! Maximum scattering order supported integer, parameter :: SCATT_ORDER_MAX = 10 character(len=*), parameter :: SCATT_ORDER_MAX_PNSTR = "scatter-p10" @@ -323,7 +328,7 @@ module constants ! Source angular distribution types integer, parameter :: & - SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular + SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular SRC_ANGLE_MONO = 2, & ! Monodirectional source SRC_ANGLE_TABULAR = 3 ! Tabular distribution @@ -333,7 +338,7 @@ module constants SRC_ENERGY_MAXWELL = 2, & ! Maxwell fission spectrum SRC_ENERGY_WATT = 3, & ! Watt fission spectrum SRC_ENERGY_TABULAR = 4 ! Tabular distribution - + ! ============================================================================ ! MISCELLANEOUS CONSTANTS diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 4c9ffb8c88..299033590d 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -93,6 +93,8 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E) then call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then + call calculate_nuclide_xs(i_nuclide, i_sab, p % E) end if ! ======================================================================== @@ -235,13 +237,8 @@ contains end if end if - ! Set last evaluated energy -- if we're in S(a,b) region, force - ! re-calculation of cross-section - if (i_sab == 0) then - micro_xs(i_nuclide) % last_E = E - else - micro_xs(i_nuclide) % last_E = ZERO - end if + micro_xs(i_nuclide) % last_E = E + micro_xs(i_nuclide) % last_index_sab = i_sab end subroutine calculate_nuclide_xs diff --git a/src/doppler.F90 b/src/doppler.F90 index 12fb34bbfe..e888eb2b18 100644 --- a/src/doppler.F90 +++ b/src/doppler.F90 @@ -87,7 +87,7 @@ contains ! ======================================================================= ! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE - if (k == 1 .and. a >= 4.0) then + if (k == 1 .and. a >= -4.0) then ! Since x = 0, this implies that a = -y F_b = F_a a = -y @@ -101,7 +101,7 @@ contains end if ! ======================================================================= - ! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4 + ! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4 k = i b = ZERO diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 62ce1c47be..4d32025699 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -166,7 +166,7 @@ contains subroutine finalize_generation() -#ifdef OPENMP +#ifdef _OPENMP ! Join the fission bank from each thread into one global fission bank call join_bank_from_threads() #endif @@ -827,7 +827,7 @@ contains end subroutine replay_batch_history -#ifdef OPENMP +#ifdef _OPENMP !=============================================================================== ! JOIN_BANK_FROM_THREADS !=============================================================================== diff --git a/src/global.F90 b/src/global.F90 index f5ce1dbcca..c27f34f9e8 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -158,7 +158,7 @@ module global ! Source and fission bank type(Bank), allocatable, target :: source_bank(:) type(Bank), allocatable, target :: fission_bank(:) -#ifdef OPENMP +#ifdef _OPENMP type(Bank), allocatable, target :: master_fission_bank(:) #endif integer(8) :: n_bank ! # of sites in fission bank @@ -206,7 +206,7 @@ module global integer :: MPI_BANK ! MPI datatype for fission bank integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult -#ifdef OPENMP +#ifdef _OPENMP integer :: n_threads = NONE ! number of OpenMP threads integer :: thread_id ! ID of a given thread #endif @@ -285,6 +285,10 @@ module global integer :: trace_gen integer(8) :: trace_particle + ! Particle tracks + logical :: write_all_tracks = .false. + integer, allocatable :: track_identifiers(:,:) + ! Particle restart run logical :: particle_restart_run = .false. @@ -440,7 +444,7 @@ contains !$omp parallel if (allocated(fission_bank)) deallocate(fission_bank) !$omp end parallel -#ifdef OPENMP +#ifdef _OPENMP if (allocated(master_fission_bank)) deallocate(master_fission_bank) #endif if (allocated(source_bank)) deallocate(source_bank) @@ -457,6 +461,9 @@ contains call active_tracklength_tallies % clear() call active_current_tallies % clear() call active_tallies % clear() + + ! Deallocate track_identifiers + if (allocated(track_identifiers)) deallocate(track_identifiers) ! Deallocate dictionaries call cell_dict % clear() diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index 915e6af070..c9e367200c 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -100,6 +100,7 @@ contains integer :: i, j, k, m integer :: n_x, n_y, n_z + integer :: length(3) integer, allocatable :: lattice_universes(:,:,:) type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() @@ -312,8 +313,8 @@ contains end do end do end do - call su % write_data(lattice_universes, "universes", & - length=(/n_x, n_y, n_z/), & + length = [n_x, n_y, n_z] + call su % write_data(lattice_universes, "universes", length=length, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) deallocate(lattice_universes) diff --git a/src/initialize.F90 b/src/initialize.F90 index 6b2023a9b6..8a4b22e90f 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -26,7 +26,7 @@ module initialize use mpi #endif -#ifdef OPENMP +#ifdef _OPENMP use omp_lib #endif @@ -408,7 +408,7 @@ contains ! Read number of threads i = i + 1 -#ifdef OPENMP +#ifdef _OPENMP ! Read and set number of OpenMP threads n_threads = str_to_int(argv(i)) if (n_threads < 1) then @@ -430,6 +430,9 @@ contains case ('-eps_tol', '-ksp_gmres_restart') ! Handle options that would be based to PETSC i = i + 1 + case ('-t', '-track', '--track') + write_all_tracks = .true. + i = i + 1 case default message = "Unknown command line option: " // argv(i) call fatal_error() @@ -863,14 +866,15 @@ contains call fatal_error() end if -#ifdef OPENMP +#ifdef _OPENMP ! If OpenMP is being used, each thread needs its own private fission ! bank. Since the private fission banks need to be combined at the end of a ! generation, there is also a 'master_fission_bank' that is used to collect ! the sites from each thread. + n_threads = omp_get_max_threads() + !$omp parallel - n_threads = omp_get_num_threads() thread_id = omp_get_thread_num() if (thread_id == 0) then diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cc12d0c41d..9fffea5d31 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -56,6 +56,7 @@ contains integer :: temp_int_array3(3) integer, allocatable :: temp_int_array(:) integer(8) :: temp_long + integer :: n_tracks logical :: file_exists logical :: check character(MAX_FILE_LEN) :: env_variable @@ -228,7 +229,7 @@ contains ! Number of OpenMP threads if (check_for_node(doc, "threads")) then -#ifdef OPENMP +#ifdef _OPENMP if (n_threads == NONE) then call get_node_value(doc, "threads", n_threads) if (n_threads < 1) then @@ -460,6 +461,25 @@ contains trace_particle = int(temp_int_array3(3), 8) end if + ! Particle tracks + if (check_for_node(doc, "track")) then + ! Make sure that there are three values per particle + n_tracks = get_arraysize_integer(doc, "track") + if (mod(n_tracks, 3) /= 0) then + message = "Number of integers specified in 'track' is not divisible & + &by 3. Please provide 3 integers per particle to be tracked." + call fatal_error() + end if + + ! Allocate space and get list of tracks + allocate(temp_int_array(n_tracks)) + call get_node_array(doc, "track", temp_int_array) + + ! Reshape into track_identifiers + allocate(track_identifiers(3, n_tracks/3)) + track_identifiers = reshape(temp_int_array, [3, n_tracks/3]) + end if + ! Shannon Entropy mesh if (check_for_node(doc, "entropy")) then @@ -2986,7 +3006,19 @@ contains end if ! set filetype, record length, and number of entries - listing % filetype = filetype + if (check_for_node(node_ace, "filetype")) then + temp_str = '' + call get_node_value(node_ace, "filetype", temp_str) + if (temp_str == 'ascii') then + listing % filetype = ASCII + else if (temp_str == 'binary') then + listing % filetype = BINARY + end if + else + listing % filetype = filetype + end if + + ! Set record length and entries for binary files if (filetype == BINARY) then listing % recl = recl listing % entries = entries diff --git a/src/output.F90 b/src/output.F90 index b1290f057e..ee0a8391c2 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -31,6 +31,10 @@ contains subroutine title() +#ifdef _OPENMP + use omp_lib +#endif + write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & ' .d88888b. 888b d888 .d8888b.', & ' d88P" "Y88b 8888b d8888 d88P Y88b', & @@ -46,25 +50,31 @@ contains ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Copyright: 2011-2013 Massachusetts Institute of Technology' + ' Copyright: 2011-2014 Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT=*) & - ' License: http://mit-crpg.github.io/openmc/license.html' - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & + ' License: http://mit-crpg.github.io/openmc/license.html' + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') & VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1 + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",7X,A)') GIT_SHA1 #endif ! Write the date and time - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') & + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",6X,A)') & time_stamp() #ifdef MPI ! Write number of processors - write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') & + write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",2X,A)') & trim(to_str(n_procs)) #endif +#ifdef _OPENMP + ! Write number of OpenMP threads + write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') & + trim(to_str(omp_get_max_threads())) +#endif + end subroutine title !=============================================================================== @@ -126,7 +136,7 @@ contains ! print header based on level select case (header_level) case (1) - write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & + write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & repeat('=', n) // '> ' // trim(line) // ' <' // & repeat('=', m), repeat('=', 75) case (2) @@ -172,6 +182,7 @@ contains write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point' write(OUTPUT_UNIT,*) ' or a particle restart file' write(OUTPUT_UNIT,*) ' -s, --threads Number of OpenMP threads' + write(OUTPUT_UNIT,*) ' -t, --track Write tracks for all particles' write(OUTPUT_UNIT,*) ' -v, --version Show version information' write(OUTPUT_UNIT,*) ' -?, --help Show this message' end if @@ -179,7 +190,7 @@ contains end subroutine print_usage !=============================================================================== -! WRITE_MESSAGE displays an informational message to the log file and the +! WRITE_MESSAGE displays an informational message to the log file and the ! standard output stream. !=============================================================================== @@ -214,7 +225,7 @@ contains ! Determine last space in current line last_space = index(message(i_start+1:i_start+line_wrap), & ' ', BACK=.true.) - if (last_space == 0) then + if (last_space == 0) then i_end = min(length + 1, i_start+line_wrap) - 1 write(ou, fmt='(1X,A)') message(i_start+1:i_end) else @@ -1027,8 +1038,10 @@ contains type(SAlphaBeta), pointer :: sab integer, optional :: unit - integer :: size_sab ! memory used by S(a,b) table - integer :: unit_ ! unit to write to + integer :: size_sab ! memory used by S(a,b) table + integer :: unit_ ! unit to write to + integer :: i ! Loop counter for parsing through sab % zaid + integer :: char_count ! Counter for the number of characters on a line ! set default unit for writing information if (present(unit)) then @@ -1039,7 +1052,30 @@ contains ! Basic S(a,b) table information write(unit_,*) 'S(a,b) Table ' // trim(sab % name) - write(unit_,*) ' zaid = ' // trim(to_str(sab % zaid)) + write(unit_,'(A)',advance="no") ' zaids = ' + ! Initialize the counter based on the above string + char_count = 11 + do i = 1, sab % n_zaid + ! Deal with a line thats too long + if (char_count >= 73) then ! 73 = 80 - (5 ZAID chars + 1 space + 1 comma) + ! End the line + write(unit_,*) "" + ! Add 11 leading blanks + write(unit_,'(A)', advance="no") " " + ! reset the counter to 11 + char_count = 11 + end if + if (i < sab % n_zaid) then + ! Include a comma + write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) // ", " + char_count = char_count + len(trim(to_str(sab % zaid(i)))) + 2 + else + ! Don't include a comma, since we are all done + write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) + end if + + end do + write(unit_,*) "" ! Move to next line write(unit_,*) ' awr = ' // trim(to_str(sab % awr)) write(unit_,*) ' kT = ' // trim(to_str(sab % kT)) @@ -1241,7 +1277,7 @@ contains end select end if write(UNIT=ou, FMT=*) - + write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "=========" write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========" if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========" @@ -1271,7 +1307,7 @@ contains if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & entropy(overall_gen) - if (overall_gen - n_inactive*gen_per_batch > 1) then + if (overall_gen - n_inactive*gen_per_batch > 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & keff, keff_std end if @@ -1299,7 +1335,7 @@ contains entropy(current_batch*gen_per_batch) ! write out accumulated k-effective if after first active batch - if (overall_gen - n_inactive*gen_per_batch > 1) then + if (overall_gen - n_inactive*gen_per_batch > 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & keff, keff_std else @@ -1309,7 +1345,7 @@ contains ! write out cmfd keff if it is active and other display info if (cmfd_on) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % k_cmfd(current_batch) + cmfd % k_cmfd(current_batch) select case(trim(cmfd_display)) case('entropy') write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & @@ -1348,7 +1384,7 @@ contains ! Plot id write(ou,100) "Plot ID:", trim(to_str(pl % id)) - + ! Plot type if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Plot Type:", "Slice" @@ -1386,7 +1422,7 @@ contains trim(to_str(pl % pixels(2))) else if (pl % type == PLOT_TYPE_VOXEL) then write(ou,100) "Voxels:", trim(to_str(pl % pixels(1))) // " " // & - trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3))) + trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3))) end if write(ou,*) @@ -1538,7 +1574,7 @@ contains write(ou,100) 'Cell ID','No. Overlap Checks' - do i = 1, n_cells + do i = 1, n_cells write(ou,101) cells(i) % id, overlap_check_cnt(i) if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1 end do @@ -1572,7 +1608,7 @@ contains integer :: k ! loop index for scoring bins integer :: n ! loop index for nuclides integer :: l ! loop index for user scores - integer :: type ! type of tally filter + integer :: type ! type of tally filter integer :: indent ! number of spaces to preceed output integer :: filter_index ! index in results array for filters integer :: score_index ! scoring bin index @@ -1748,13 +1784,13 @@ contains score_name = 'P' // trim(to_str(t % scatt_order(k))) // & ' Scattering Moment' end if - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) else if (t % score_bins(k) == SCORE_SCATTER_PN) then score_name = "Scattering Rate" - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) @@ -1762,7 +1798,7 @@ contains score_index = score_index + 1 score_name = 'P' // trim(to_str(n_order)) // & ' Scattering Moment' - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) @@ -1774,7 +1810,7 @@ contains else score_name = score_names(abs(t % score_bins(k))) end if - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) @@ -1812,8 +1848,8 @@ contains integer :: i_filter_ein ! index for incoming energy filter integer :: i_filter_surf ! index for surface filter integer :: n ! number of incoming energy bins - integer :: len1 ! length of string - integer :: len2 ! length of string + integer :: len1 ! length of string + integer :: len2 ! length of string integer :: filter_index ! index in results array for filters logical :: print_ebin ! should incoming energy bin be displayed? character(MAX_LINE_LEN) :: string @@ -1863,14 +1899,14 @@ contains mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1880,14 +1916,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1897,14 +1933,14 @@ contains mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1914,14 +1950,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Front", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Front", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1931,14 +1967,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) matching_bins(i_filter_surf) = IN_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1948,14 +1984,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Top", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Top", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 201de518eb..e1c696ff16 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -76,6 +76,9 @@ module particle_header ! Statistical data integer :: n_collision ! # of collisions + ! Track output + logical :: write_track = .false. + contains procedure :: initialize => initialize_particle procedure :: clear => clear_particle diff --git a/src/physics.F90 b/src/physics.F90 index fb1f891943..2d3cfeb39f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -47,7 +47,7 @@ contains if (verbosity >= 10 .or. trace) then message = " " // trim(reaction_name(p % event_MT)) // " with " // & trim(adjustl(nuclides(p % event_nuclide) % name)) // & - ". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV." + ". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV." call write_message() end if @@ -158,7 +158,7 @@ contains call fatal_error() end if - ! Find atom density + ! Find atom density i_nuclide = mat % nuclide(i) atom_density = mat % atom_density(i) @@ -209,7 +209,7 @@ contains i_reaction = nuc % index_fission(1) return end if - + ! Get grid index and interpolatoin factor and sample fission cdf i_grid = micro_xs(i_nuclide) % index_grid f = micro_xs(i_nuclide) % interp_factor @@ -226,7 +226,7 @@ contains if (i_grid < rxn % threshold) cycle ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + f*(rxn%sigma(i_grid - rxn%threshold + 2))) ! Create fission bank sites if fission occus @@ -382,7 +382,7 @@ contains if (i_grid < rxn % threshold) cycle ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + f*(rxn%sigma(i_grid - rxn%threshold + 2))) end do @@ -490,13 +490,24 @@ contains integer :: k ! outgoing cosine bin integer :: n_energy_out ! number of outgoing energy bins real(8) :: f ! interpolation factor - real(8) :: r ! used for skewed sampling + real(8) :: r ! used for skewed sampling & continuous real(8) :: E_ij ! outgoing energy j for E_in(i) real(8) :: E_i1j ! outgoing energy j for E_in(i+1) real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j) real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j) real(8) :: prob ! probability for sampling Bragg edge type(SAlphaBeta), pointer, save :: sab => null() + ! Following are needed only for SAB_SECONDARY_CONT scattering + integer :: l ! sampled incoming E bin (is i or i + 1) + real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_J ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_J ! endpoints interpolated between i and i+1 + real(8) :: E_l_j, E_l_j1 ! adjacent E on outgoing grid l + real(8) :: p_l_j, p_l_j1 ! adjacent p on outgoing grid l + real(8) :: c_j, c_j1 ! cumulative probability + real(8) :: frac ! interpolation factor on outgoing energy + real(8) :: r1 ! RNG for outgoing energy + !$omp threadprivate(sab) ! Get pointer to S(a,b) table @@ -513,7 +524,7 @@ contains f = ZERO else i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E) - f = (E - sab%elastic_e_in(i)) / & + f = (E - sab%elastic_e_in(i)) / & (sab%elastic_e_in(i+1) - sab%elastic_e_in(i)) end if @@ -554,8 +565,7 @@ contains ! Outgoing energy is same as incoming energy -- no need to do anything else - ! Determine number of outgoing energy and angle bins - n_energy_out = sab % n_inelastic_e_out + ! Perform inelastic calculations ! Get index and interpolation factor for inelastic grid if (E < sab % inelastic_e_in(1)) then @@ -563,61 +573,150 @@ contains f = ZERO else i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) - f = (E - sab%inelastic_e_in(i)) / & + f = (E - sab%inelastic_e_in(i)) / & (sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i)) end if ! Now that we have an incoming energy bin, we need to determine the ! outgoing energy bin. This will depend on the "secondary energy ! mode". If the mode is 0, then the outgoing energy bin is chosen from a - ! set of equally-likely bins. However, if the mode is 1, then the first + ! set of equally-likely bins. If the mode is 1, then the first ! two and last two bins are skewed to have lower probabilities than the ! other bins (0.1 for the first and last bins and 0.4 for the second and - ! second to last bins, relative to a normal bin probability of 1) + ! second to last bins, relative to a normal bin probability of 1). + ! Finally, if the mode is 2, then a continuous distribution (with + ! accompanying PDF and CDF is utilized) - if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then - ! All bins equally likely - j = 1 + int(prn() * n_energy_out) - elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then - r = prn() * (n_energy_out - 3) - if (r > ONE) then - ! equally likely N-4 middle bins - j = int(r) + 2 - elseif (r > 0.6) then - ! second to last bin has relative probability of 0.4 - j = n_energy_out - 1 - elseif (r > 0.5) then - ! last bin has relative probability of 0.1 - j = n_energy_out - elseif (r > 0.1) then - ! second bin has relative probability of 0.4 - j = 2 - else - ! first bin has relative probability of 0.1 - j = 1 + if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. & + (sab % secondary_mode == SAB_SECONDARY_SKEWED)) then + if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then + ! All bins equally likely + + j = 1 + int(prn() * sab % n_inelastic_e_out) + elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then + ! Distribution skewed away from edge points + + ! Determine number of outgoing energy and angle bins + n_energy_out = sab % n_inelastic_e_out + + r = prn() * (n_energy_out - 3) + if (r > ONE) then + ! equally likely N-4 middle bins + j = int(r) + 2 + elseif (r > 0.6) then + ! second to last bin has relative probability of 0.4 + j = n_energy_out - 1 + elseif (r > 0.5) then + ! last bin has relative probability of 0.1 + j = n_energy_out + elseif (r > 0.1) then + ! second bin has relative probability of 0.4 + j = 2 + else + ! first bin has relative probability of 0.1 + j = 1 + end if end if + + ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) + E_ij = sab % inelastic_e_out(j,i) + E_i1j = sab % inelastic_e_out(j,i+1) + + ! Outgoing energy + E = (1 - f)*E_ij + f*E_i1j + + ! Sample outgoing cosine bin + k = 1 + int(prn() * sab % n_inelastic_mu) + + ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) + mu_ijk = sab % inelastic_mu(k,j,i) + mu_i1jk = sab % inelastic_mu(k,j,i+1) + + ! Cosine of angle between incoming and outgoing neutron + mu = (1 - f)*mu_ijk + f*mu_i1jk + + else if (sab % secondary_mode == SAB_SECONDARY_CONT) then + ! Continuous secondary energy - this is to be similar to + ! Law 61 interpolation on outgoing energy + + ! Sample between ith and (i+1)th bin + r = prn() + if (f > r) then + l = i + 1 + else + l = i + end if + + ! Determine endpoints on grid i + n_energy_out = sab % inelastic_data(i) % n_e_out + E_i_1 = sab % inelastic_data(i) % e_out(1) + E_i_J = sab % inelastic_data(i) % e_out(n_energy_out) + + ! Determine endpoints on grid i + 1 + n_energy_out = sab % inelastic_data(i + 1) % n_e_out + E_i1_1 = sab % inelastic_data(i + 1) % e_out(1) + E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out) + + E_1 = E_i_1 + f * (E_i1_1 - E_i_1) + E_J = E_i_J + f * (E_i1_J - E_i_J) + + ! Determine outgoing energy bin + ! (First reset n_energy_out to the right value) + n_energy_out = sab % inelastic_data(l) % n_e_out + r1 = prn() + c_j = sab % inelastic_data(l) % e_out_cdf(1) + do j = 1, n_energy_out - 1 + c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1) + if (r1 < c_j1) exit + c_j = c_j1 + end do + + ! check to make sure k is <= n_energy_out - 1 + j = min(j, n_energy_out - 1) + + ! Get the data to interpolate between + E_l_j = sab % inelastic_data(l) % e_out(j) + p_l_j = sab % inelastic_data(l) % e_out_pdf(j) + + ! Next part assumes linear-linear interpolation in standard + E_l_j1 = sab % inelastic_data(l) % e_out(j + 1) + p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1) + + ! Find secondary energy (variable E) + frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j) + if (frac == ZERO) then + E = E_l_j + (r1 - c_j) / p_l_j + else + E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + & + TWO * frac * (r1 - c_j))) - p_l_j) / frac + end if + + ! Now interpolate between incident energy bins i and i + 1 + if (l == i) then + E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1) + else + E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1) + end if + + ! Find angular distribution for closest outgoing energy bin + if (r1 - c_j < c_j1 - r1) then + j = j + else + j = j + 1 + end if + + ! Sample outgoing cosine bin + k = 1 + int(prn() * sab % n_inelastic_mu) + + ! Will use mu from the randomly chosen incoming and closest outgoing + ! energy bins + mu = sab % inelastic_data(l) % mu(k, j) + else message = "Invalid secondary energy mode on S(a,b) table " // & trim(sab % name) - end if - - ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) - E_ij = sab % inelastic_e_out(j,i) - E_i1j = sab % inelastic_e_out(j,i+1) - - ! Outgoing energy - E = (1 - f)*E_ij + f*E_i1j - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_inelastic_mu) - - ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) - mu_ijk = sab % inelastic_mu(k,j,i) - mu_i1jk = sab % inelastic_mu(k,j,i+1) - - ! Cosine of angle between incoming and outgoing neutron - mu = (1 - f)*mu_ijk + f*mu_i1jk - end if + end if ! (inelastic secondary energy treatment) + end if ! (elastic or inelastic) ! change direction of particle uvw = rotate_angle(uvw, mu) @@ -974,7 +1073,7 @@ contains E_cm = E ! determine outgoing energy in lab - E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & + E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & / ((A+ONE)*(A+ONE)) ! determine outgoing angle in lab @@ -1037,7 +1136,7 @@ contains r = ONE else i = binary_search(rxn % adist % energy, n, E) - r = (E - rxn % adist % energy(i)) / & + r = (E - rxn % adist % energy(i)) / & (rxn % adist % energy(i+1) - rxn % adist % energy(i)) end if @@ -1166,7 +1265,7 @@ contains end if end function rotate_angle - + !=============================================================================== ! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary ! neutron from a collision or for a prompt/delayed fission neutron @@ -1322,7 +1421,7 @@ contains ! ======================================================================= ! CONTINUOUS TABULAR DISTRIBUTION - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR == 1) then @@ -1348,7 +1447,7 @@ contains r = ONE else i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & + r = (E_in - edist%data(lc+i)) / & (edist%data(lc+i+1) - edist%data(lc+i)) end if @@ -1452,7 +1551,7 @@ contains ! ======================================================================= ! MAXWELL FISSION SPECTRUM - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) @@ -1484,7 +1583,7 @@ contains ! ======================================================================= ! EVAPORATION SPECTRUM - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) @@ -1565,7 +1664,7 @@ contains call fatal_error() end if - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR > 0) then @@ -1587,7 +1686,7 @@ contains r = ONE else i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & + r = (E_in - edist%data(lc+i)) / & (edist%data(lc+i+1) - edist%data(lc+i)) end if @@ -1714,11 +1813,11 @@ contains if (.not. present(mu_out)) then ! call write_particle_restart() - message = "Law 44 called without giving mu_out as argument." + message = "Law 61 called without giving mu_out as argument." call fatal_error() end if - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR > 0) then @@ -1740,7 +1839,7 @@ contains r = ONE else i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & + r = (E_in - edist%data(lc+i)) / & (edist%data(lc+i+1) - edist%data(lc+i)) end if diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 83ff6409f9..e3a9885af1 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -5,8 +5,13 @@ module random_lcg private save + ! Random number streams + integer, parameter :: N_STREAMS = 2 + integer, parameter :: STREAM_TRACKING = 1 + integer, parameter :: STREAM_TALLIES = 2 + integer(8) :: prn_seed0 ! original seed - integer(8) :: prn_seed ! current seed + integer(8) :: prn_seed(N_STREAMS) ! current seed integer(8) :: prn_mult ! multiplication factor, g integer(8) :: prn_add ! additive factor, c integer :: prn_bits ! number of bits, M @@ -14,6 +19,7 @@ module random_lcg integer(8) :: prn_mask ! 2^M - 1 integer(8) :: prn_stride ! stride between particles real(8) :: prn_norm ! 2^(-M) + integer :: stream ! current RNG stream !$omp threadprivate(prn_seed) @@ -21,6 +27,8 @@ module random_lcg public :: initialize_prng public :: set_particle_seed public :: prn_skip + public :: prn_set_stream + public :: STREAM_TRACKING, STREAM_TALLIES contains @@ -35,12 +43,12 @@ contains ! This algorithm uses bit-masking to find the next integer(8) value to be ! used to calculate the random number - prn_seed = iand(prn_mult*prn_seed + prn_add, prn_mask) + prn_seed(stream) = iand(prn_mult*prn_seed(stream) + prn_add, prn_mask) ! Once the integer is calculated, we just need to divide by 2**m, ! represented here as multiplying by a pre-calculated factor - pseudo_rn = prn_seed * prn_norm + pseudo_rn = prn_seed(stream) * prn_norm end function prn @@ -53,8 +61,13 @@ contains use global, only: seed + integer :: i + + stream = STREAM_TRACKING prn_seed0 = seed - prn_seed = prn_seed + do i = 1, N_STREAMS + prn_seed(i) = prn_seed0 + i - 1 + end do prn_mult = 2806196910506780709_8 prn_add = 1_8 prn_bits = 63 @@ -74,7 +87,11 @@ contains integer(8), intent(in) :: id - prn_seed = prn_skip_ahead(id*prn_stride, prn_seed0) + integer :: i + + do i = 1, N_STREAMS + prn_seed(i) = prn_skip_ahead(id*prn_stride, prn_seed0 + i - 1) + end do end subroutine set_particle_seed @@ -86,7 +103,7 @@ contains integer(8), intent(in) :: n ! number of seeds to skip - prn_seed = prn_skip_ahead(n, prn_seed) + prn_seed(stream) = prn_skip_ahead(n, prn_seed(stream)) end subroutine prn_skip @@ -151,4 +168,18 @@ contains end function prn_skip_ahead +!=============================================================================== +! PRN_SET_STREAM changes the random number stream. If random numbers are needed +! in routines not used directly for tracking (e.g. physics), this allows the +! numbers to be generated without affecting reproducibility of the physics. +!=============================================================================== + + subroutine prn_set_stream(i) + + integer, intent(in) :: i + + stream = i + + end subroutine prn_set_stream + end module random_lcg diff --git a/src/relaxng/cross_sections.rnc b/src/relaxng/cross_sections.rnc index e82bb986fb..75c0ea29c2 100644 --- a/src/relaxng/cross_sections.rnc +++ b/src/relaxng/cross_sections.rnc @@ -10,7 +10,9 @@ element cross_sections { (element temperature { xsd:double } | attribute temperature { xsd:double }) & (element path { xsd:string { maxLength = "255" } } | attribute path { xsd:string { maxLength = "255" } }) & - (element location { xsd:int } | attribute location { xsd:int })? + (element location { xsd:int } | attribute location { xsd:int })? & + (element filetype { ( "ascii" | "binary" ) } | + attribute filetype { ( "ascii" | "binary" ) })? }* & element directory { xsd:string { maxLength = "255" } }? & diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 99f64d1b4a..0e4942e357 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -116,6 +116,8 @@ element settings { element trace { list { xsd:positiveInteger+ } }? & + element track { list { xsd:positiveInteger+ } }? & + element verbosity { xsd:positiveInteger }? & element uniform_fs{ diff --git a/src/source.F90 b/src/source.F90 index 67d8ec8013..231a23cc5c 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -3,6 +3,7 @@ module source use bank_header, only: Bank use constants use error, only: fatal_error + use geometry, only: find_cell use geometry_header, only: BASE_UNIVERSE use global use math, only: maxwell_spectrum, watt_spectrum @@ -73,6 +74,9 @@ contains real(8) :: p_max(3) ! maximum coordinates of source real(8) :: a ! Arbitrary parameter 'a' real(8) :: b ! Arbitrary parameter 'b' + logical :: found ! Does the source particle exist within geometry? + type(Particle) :: p ! Temporary particle for using find_cell + integer, save :: num_resamples = 0 ! Number of resamples encountered ! Set weight to one by default site % wgt = ONE @@ -80,11 +84,32 @@ contains ! Sample position select case (external_source % type_space) case (SRC_SPACE_BOX) - ! Coordinates sampled uniformly over a box - p_min = external_source % params_space(1:3) - p_max = external_source % params_space(4:6) - r = (/ (prn(), i = 1,3) /) - site % xyz = p_min + r*(p_max - p_min) + ! Set particle defaults + call p % initialize() + ! Repeat sampling source location until a good site has been found + found = .false. + do while (.not.found) + ! Coordinates sampled uniformly over a box + p_min = external_source % params_space(1:3) + p_max = external_source % params_space(4:6) + r = (/ (prn(), i = 1,3) /) + site % xyz = p_min + r*(p_max - p_min) + + ! Fill p with needed data + p % coord0 % xyz = site % xyz + p % coord0 % uvw = [ ONE, ZERO, ZERO ] + + ! Now search to see if location exists in geometry + call find_cell(p, found) + if (.not. found) then + num_resamples = num_resamples + 1 + if (num_resamples == MAX_EXTSRC_RESAMPLES) then + message = "Maximum number of external source spatial resamples & + &reached!" + call fatal_error() + end if + end if + end do case (SRC_SPACE_POINT) ! Point source @@ -146,7 +171,7 @@ contains end subroutine sample_external_source !=============================================================================== -! GET_SOURCE_PARTICLE returns the next source particle +! GET_SOURCE_PARTICLE returns the next source particle !=============================================================================== subroutine get_source_particle(p, index_source) @@ -155,6 +180,7 @@ contains integer(8), intent(in) :: index_source integer(8) :: particle_seed ! unique index for particle + integer :: i type(Bank), pointer, save :: src => null() !$omp threadprivate(src) @@ -177,6 +203,21 @@ contains if (current_batch == trace_batch .and. current_gen == trace_gen .and. & p % id == trace_particle) trace = .true. + ! Set particle track. + p % write_track = .false. + if (write_all_tracks) then + p % write_track = .true. + else if (allocated(track_identifiers)) then + do i=1, size(track_identifiers(1,:)) + if (current_batch == track_identifiers(1,i) .and. & + ¤t_gen == track_identifiers(2,i) .and. & + &p % id == track_identifiers(3,i)) then + p % write_track = .true. + exit + end if + end do + end if + end subroutine get_source_particle !=============================================================================== diff --git a/src/state_point.F90 b/src/state_point.F90 index 306fb8f37f..046d748be5 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -515,6 +515,7 @@ contains character(19) :: current_time integer :: i integer :: j + integer :: length(4) integer :: int_array(3) integer, allocatable :: temp_array(:) real(8) :: real_array(3) @@ -593,10 +594,9 @@ contains call sp % read_data(cmfd % indices, "indicies", length=4, group="cmfd") call sp % read_data(cmfd % k_cmfd, "k_cmfd", length=restart_batch, & group="cmfd") + length = cmfd % indices([4,1,2,3]) call sp % read_data(cmfd % cmfd_src, "cmfd_src", & - length=(/cmfd % indices(4), cmfd % indices(1), & - cmfd % indices(2), cmfd % indices(3)/), & - group="cmfd") + length=length, group="cmfd") call sp % read_data(cmfd % entropy, "cmfd_entropy", & length=restart_batch, group="cmfd") call sp % read_data(cmfd % balance, "cmfd_balance", & diff --git a/src/track_output.F90 b/src/track_output.F90 new file mode 100644 index 0000000000..79a2809066 --- /dev/null +++ b/src/track_output.F90 @@ -0,0 +1,79 @@ +!=============================================================================== +! TRACK_OUTPUT handles output of particle tracks (the paths taken by particles +! as they are transported through the geometry). +!=============================================================================== + +module track_output + + use global + use output_interface, only: BinaryOutput + use particle_header, only: Particle + use string, only: to_str + + implicit none + + integer, private :: n_tracks ! total number of tracks + real(8), private, allocatable :: coords(:,:) ! track coordinates +!$omp threadprivate(n_tracks, coords) + +contains + +!=============================================================================== +! INITIALIZE_PARTICLE_TRACK +!=============================================================================== + + subroutine initialize_particle_track() + n_tracks = 0 + end subroutine initialize_particle_track + +!=============================================================================== +! WRITE_PARTICLE_TRACK copies particle position to an array. +!=============================================================================== + + subroutine write_particle_track(p) + type(Particle), intent(in) :: p + + real(8), allocatable :: new_coords(:, :) + + ! Add another column to coords + n_tracks = n_tracks + 1 + if (allocated(coords)) then + allocate(new_coords(3, n_tracks)) + new_coords(:, 1:n_tracks-1) = coords + call move_alloc(FROM=new_coords, TO=coords) + else + allocate(coords(3,1)) + end if + + ! Write current coordinates into the newest column. + coords(:, n_tracks) = p % coord0 % xyz + end subroutine write_particle_track + +!=============================================================================== +! FINALIZE_PARTICLE_TRACK writes the particle track array to disk. +!=============================================================================== + + subroutine finalize_particle_track(p) + type(Particle), intent(in) :: p + + integer :: length(2) + character(MAX_FILE_LEN) :: fname + type(BinaryOutput) :: binout + +#ifdef HDF5 + fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) & + // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) & + // '.h5' +#else + fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) & + // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) & + // '.binary' +#endif + call binout % file_create(fname) + length = [3, n_tracks] + call binout % write_data(coords, 'coordinates', length=length) + call binout % file_close() + deallocate(coords) + end subroutine finalize_particle_track + +end module track_output diff --git a/src/tracking.F90 b/src/tracking.F90 index 04175b243f..c7ee31f0b2 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -13,6 +13,8 @@ module tracking use string, only: to_str use tally, only: score_analog_tally, score_tracklength_tally, & score_surface_current + use track_output, only: initialize_particle_track, write_particle_track, & + finalize_particle_track contains @@ -67,8 +69,16 @@ contains ! Force calculation of cross-sections by setting last energy to zero micro_xs % last_E = ZERO + ! Prepare to write out particle track. + if (p % write_track) then + call initialize_particle_track() + endif + do while (p % alive) + ! Write particle track. + if (p % write_track) call write_particle_track(p) + if (check_overlaps) call check_cell_overlap(p) ! Calculate microscopic and macroscopic cross sections -- note: if the @@ -196,6 +206,12 @@ contains end do + ! Finish particle track output. + if (p % write_track) then + call write_particle_track(p) + call finalize_particle_track(p) + endif + end subroutine transport end module tracking diff --git a/src/utils/convert_xsdir.py b/src/utils/convert_xsdir.py index fa29f46eca..edee1a66c7 100755 --- a/src/utils/convert_xsdir.py +++ b/src/utils/convert_xsdir.py @@ -35,8 +35,12 @@ class Xsdir(object): words = line.split() if words: if words[0].lower().startswith('datapath'): - index = line.index('=') - self.datapath = line[index+1:].strip() + if '=' in words[0]: + index = line.index('=') + self.datapath = line[index+1:].strip() + else: + if len(line.strip()) > 8: + self.datapath = line[8:].strip() else: self.f.seek(0) @@ -71,7 +75,7 @@ class Xsdir(object): # Handle continuation lines while words[-1] == '+': extraWords = self.f.readline().split() - words = words + extraWords + words = words[:-1] + extraWords assert len(words) >= 7 # Create XsdirTable object and add to line diff --git a/src/utils/plot_mesh_tally.py b/src/utils/plot_mesh_tally.py index 04b6b4592a..fb5b8cb839 100755 --- a/src/utils/plot_mesh_tally.py +++ b/src/utils/plot_mesh_tally.py @@ -1,262 +1,279 @@ #!/usr/bin/env python2 -'''Python script to plot tally data generated by OpenMC.''' +"""Python script to plot tally data generated by OpenMC.""" +import os import sys -from statepoint import * -# Color intensity dependent on individual score? - -from PyQt4.QtCore import * -from PyQt4.QtGui import * -import matplotlib.pyplot as plt +from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg +from matplotlib.backends.backend_tkagg import NavigationToolbar2TkAgg from matplotlib.figure import Figure -from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas -from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar +import matplotlib.pyplot as plt import numpy as np -class AppForm(QMainWindow): - def __init__(self, parent=None): - QMainWindow.__init__(self, parent) +from statepoint import * + +if sys.version_info[0] < 3: + import Tkinter as tk +else: + import tkinter as tk +import tkFileDialog +import tkFont +import tkMessageBox +import ttk + + +class MeshPlotter(tk.Frame): + def __init__(self, parent, filename): + tk.Frame.__init__(self, parent) + + self.labels = {'cell': 'Cell:', 'cellborn': 'Cell born:', + 'surface': 'Surface:', 'material': 'Material:', + 'universe': 'Universe:', 'energyin': 'Energy in:', + 'energyout': 'Energy out:'} + + self.filterBoxes = {} # Read data from source or leakage fraction file - self.get_file_data() - self.main_frame = QWidget() - self.setCentralWidget(self.main_frame) - - # Create the Figure, Canvas, and Axes + self.get_file_data(filename) + + # Set up top-level window + top = self.winfo_toplevel() + top.title('Mesh Tally Plotter: ' + filename) + top.rowconfigure(0, weight=1) + top.columnconfigure(0, weight=1) + self.grid(sticky=tk.W+tk.N) + + # Create widgets and draw to screen + self.create_widgets() + self.update() + + def create_widgets(self): + figureFrame = tk.Frame(self) + figureFrame.grid(row=0, column=0) + + # Create the Figure and Canvas self.dpi = 100 - self.fig = Figure((5.0, 15.0), dpi=self.dpi) - self.canvas = FigureCanvas(self.fig) - self.canvas.setParent(self.main_frame) - self.axes = self.fig.add_subplot(111) - + self.fig = Figure((5.0, 5.0), dpi=self.dpi) + self.canvas = FigureCanvasTkAgg(self.fig, master=figureFrame) + self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1) + # Create the navigation toolbar, tied to the canvas - self.mpl_toolbar = NavigationToolbar(self.canvas, self.main_frame) + self.mpl_toolbar = NavigationToolbar2TkAgg(self.canvas, figureFrame) + self.mpl_toolbar.update() + self.canvas._tkcanvas.pack(side=tk.TOP, fill=tk.BOTH, expand=1) - # Grid layout at bottom - self.grid = QGridLayout() + # Create frame for comboboxes + self.selectFrame = tk.Frame(self) + self.selectFrame.grid(row=1, column=0, sticky=tk.W+tk.E) - # Overall layout - self.vbox = QVBoxLayout() - self.vbox.addWidget(self.canvas) - self.vbox.addWidget(self.mpl_toolbar) - self.vbox.addLayout(self.grid) - self.main_frame.setLayout(self.vbox) + # Tally selection + labelTally = tk.Label(self.selectFrame, text='Tally:') + labelTally.grid(row=0, column=0, sticky=tk.W) + self.tallyBox = ttk.Combobox(self.selectFrame, state='readonly') + self.tallyBox['values'] = [self.datafile.tallies[i].id + for i in self.meshTallies] + self.tallyBox.current(0) + self.tallyBox.grid(row=0, column=1, sticky=tk.W+tk.E) + self.tallyBox.bind('<>', self.update) - # Tally selections - label_tally = QLabel("Tally:") - self.tally = QComboBox() - self.tally.addItems([(str(i + 1)) for i in range(self.n_tallies)]) - self.connect(self.tally, SIGNAL('activated(int)'), - self._update) - self.connect(self.tally, SIGNAL('activated(int)'), - self.populate_boxes) - self.connect(self.tally, SIGNAL('activated(int)'), - self.on_draw) + # Planar basis selection + labelBasis = tk.Label(self.selectFrame, text='Basis:') + labelBasis.grid(row=1, column=0, sticky=tk.W) + self.basisBox = ttk.Combobox(self.selectFrame, state='readonly') + self.basisBox['values'] = ('xy', 'yz', 'xz') + self.basisBox.current(0) + self.basisBox.grid(row=1, column=1, sticky=tk.W+tk.E) + self.basisBox.bind('<>', self.update) - # Planar basis - label_basis = QLabel("Basis:") - self.basis = QComboBox() - self.basis.addItems(['xy', 'yz', 'xz']) + # Axial level + labelAxial = tk.Label(self.selectFrame, text='Axial level:') + labelAxial.grid(row=2, column=0, sticky=tk.W) + self.axialBox = ttk.Combobox(self.selectFrame, state='readonly') + self.axialBox.grid(row=2, column=1, sticky=tk.W+tk.E) + self.axialBox.bind('<>', self.redraw) - # Update window when 'Basis' selection is changed - self.connect(self.basis, SIGNAL('activated(int)'), - self._update) - self.connect(self.basis, SIGNAL('activated(int)'), - self.populate_boxes) - self.connect(self.basis, SIGNAL('activated(int)'), - self.on_draw) + # Option for mean/uncertainty + labelMean = tk.Label(self.selectFrame, text='Mean/Uncertainty:') + labelMean.grid(row=3, column=0, sticky=tk.W) + self.meanBox = ttk.Combobox(self.selectFrame, state='readonly') + self.meanBox['values'] = ('Mean', 'Absolute uncertainty', + 'Relative uncertainty') + self.meanBox.current(0) + self.meanBox.grid(row=3, column=1, sticky=tk.W+tk.E) + self.meanBox.bind('<>', self.update) - # Axial level within selected basis - label_axial_level = QLabel("Axial Level:") - self.axial_level = QComboBox() - self.connect(self.axial_level, SIGNAL('activated(int)'), - self.on_draw) - - # Add Option to plot mean or uncertainty - label_mean = QLabel("Mean or Uncertainty:") - self.mean = QComboBox() - self.mean.addItems(['Mean','Absolute Uncertainty', - 'Relative Uncertainty']) - - # Update window when mean selection is changed - self.connect(self.mean, SIGNAL('activated(int)'), - self.on_draw) - + # Scores + labelScore = tk.Label(self.selectFrame, text='Score:') + labelScore.grid(row=4, column=0, sticky=tk.W) + self.scoreBox = ttk.Combobox(self.selectFrame, state='readonly') + self.scoreBox.grid(row=4, column=1, sticky=tk.W+tk.E) + self.scoreBox.bind('<>', self.redraw) - self.label_filters = QLabel("Filter options:") + # Filter label + font = tkFont.Font(weight='bold') + labelFilters = tk.Label(self.selectFrame, text='Filters:', font=font) + labelFilters.grid(row=5, column=0, sticky=tk.W) - # Labels for all possible filters - self.labels = {'cell': 'Cell: ', 'cellborn': 'Cell born: ', - 'surface': 'Surface: ', 'material': 'Material', - 'universe': 'Universe: ', 'energyin': 'Energy in: ', - 'energyout': 'Energy out: '} + def update(self, event=None): + if not event: + widget = None + else: + widget = event.widget - # Empty reusable labels - self.qlabels = {} - for j in range(8): - self.nextLabel = QLabel - self.qlabels[j] = self.nextLabel + tally_id = self.meshTallies[self.tallyBox.current()] + selectedTally = self.datafile.tallies[tally_id] - # Reusable comboboxes labelled with filter names - self.boxes = {} - for key in self.labels.keys(): - self.nextBox = QComboBox() - self.connect(self.nextBox, SIGNAL('activated(int)'), - self.on_draw) - self.boxes[key] = self.nextBox + # Get mesh for selected tally + self.mesh = self.datafile.meshes[ + selectedTally.filters['mesh'].bins[0] - 1] - # Combobox to select among scores - self.score_label = QLabel("Score:") - self.scoreBox = QComboBox() - for item in self.tally_scores[0]: - self.scoreBox.addItems(str(item)) - self.connect(self.scoreBox, SIGNAL('activated(int)'), - self.on_draw) + # Get mesh dimensions + self.nx, self.ny, self.nz = self.mesh.dimension - # Fill layout - self.grid.addWidget(label_tally, 0, 0) - self.grid.addWidget(self.tally, 0, 1) - self.grid.addWidget(label_basis, 1, 0) - self.grid.addWidget(self.basis, 1, 1) - self.grid.addWidget(label_axial_level, 2, 0) - self.grid.addWidget(self.axial_level, 2, 1) - self.grid.addWidget(label_mean, 3, 0) - self.grid.addWidget(self.mean, 3, 1) - self.grid.addWidget(self.label_filters, 4, 0) + # Repopulate comboboxes baesd on current basis selection + text = self.basisBox['values'][self.basisBox.current()] + if text == 'xy': + self.axialBox['values'] = [str(i+1) for i in range(self.nz)] + elif text == 'yz': + self.axialBox['values'] = [str(i+1) for i in range(self.nx)] + else: + self.axialBox['values'] = [str(i+1) for i in range(self.ny)] + self.axialBox.current(0) - self._update() - self.populate_boxes() - self.on_draw() + # If update() was called by a change in the basis combobox, we don't + # need to repopulate the filters + if widget == self.basisBox: + self.redraw() + return - def get_file_data(self): - # Get data file name from "open file" browser - filename = QFileDialog.getOpenFileName(self, 'Select statepoint file', '.') + # Update scores + self.scoreBox['values'] = selectedTally.scores + self.scoreBox.current(0) - # Create StatePoint object and read in data - self.datafile = StatePoint(str(filename)) - self.datafile.read_results() - self.datafile.generate_stdev() + # Remove any filter labels/comboboxes that exist + for row in range(6, self.selectFrame.grid_size()[1]): + for w in self.selectFrame.grid_slaves(row=row): + w.grid_forget() + w.destroy() - self.setWindowTitle('Core Map Tool : ' + str(self.datafile.path)) + # create a label/combobox for each filter in selected tally + count = 0 + for filterType in selectedTally.filters: + if filterType == 'mesh': + continue + count += 1 - # Set maximum colorbar value by maximum tally data value - self.maxvalue = self.datafile.tallies[0].results.max() + # Create label and combobox for this filter + label = tk.Label(self.selectFrame, text=self.labels[filterType]) + label.grid(row=count+6, column=0, sticky=tk.W) + combobox = ttk.Combobox(self.selectFrame, state='readonly') + self.filterBoxes[filterType] = combobox - self.labelList = [] + # Set combobox items + f = selectedTally.filters[filterType] + if filterType in ['energyin', 'energyout']: + combobox['values'] = ['{0} to {1}'.format(*f.bins[i:i+2]) + for i in range(f.length)] + else: + combobox['values'] = [str(i) for i in f.bins] - # Read mesh dimensions -# for mesh in self.datafile.meshes: -# self.nx, self.ny, self.nz = mesh.dimension + combobox.current(0) + combobox.grid(row=count+6, column=1, sticky=tk.W+tk.E) + combobox.bind('<>', self.redraw) - # Read filter types from statepoint file - self.n_tallies = len(self.datafile.tallies) - self.tally_list = [] - for tally in self.datafile.tallies: - self.filter_types = [] - for f in tally.filters: - self.filter_types.append(f) - self.tally_list.append(self.filter_types) + # If There are no filters, leave a 'None available' message + if count == 0: + count += 1 + label = tk.Label(self.selectFrame, text="None Available") + label.grid(row=count+6, column=0, sticky=tk.W) - # Read score types from statepoint file - self.tally_scores = [] - for tally in self.datafile.tallies: - self.score_types = [] - for s in tally.scores: - self.score_types.append(s) - self.tally_scores.append(self.score_types) -# print 'self.tally_scores = ', self.tally_scores + self.redraw() - def on_draw(self): - """ Redraws the figure - """ + def redraw(self, event=None): + basis = self.basisBox.current() + 1 + axial_level = self.axialBox.current() + 1 + is_mean = self.meanBox.current() -# print 'Calling on_draw...' - # Get selected basis, axial_level and stage - basis = self.basis.currentIndex() + 1 - axial_level = self.axial_level.currentIndex() + 1 - is_mean = self.mean.currentIndex() + # Get selected tally + tally_id = self.meshTallies[self.tallyBox.current()] + selectedTally = self.datafile.tallies[tally_id] # Create spec_list spec_list = [] - for tally in self.datafile.tallies[self.tally.currentIndex()].filters.values(): - if tally.type == 'mesh': + for f in selectedTally.filters.values(): + if f.type == 'mesh': continue - index = self.boxes[tally.type].currentIndex() - spec_list.append((tally.type, index)) - + index = self.filterBoxes[f.type].current() + spec_list.append((f.type, index)) + # Take is_mean and convert it to an index of the score score_loc = is_mean if score_loc > 1: score_loc = 1 - - if self.basis.currentText() == 'xy': + + text = self.basisBox['values'][self.basisBox.current()] + if text == 'xy': matrix = np.zeros((self.nx, self.ny)) for i in range(self.nx): for j in range(self.ny): - matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, j, axial_level))], - self.scoreBox.currentIndex())[score_loc] - # Calculate relative uncertainty from absolute, if - # requested + matrix[i, j] = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, j + 1, axial_level))], + self.scoreBox.current())[score_loc] + # Calculate relative uncertainty from absolute, if requested if is_mean == 2: # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, j, axial_level))], - self.scoreBox.currentIndex())[0] + mean_val = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, j + 1, axial_level))], + self.scoreBox.current())[0] if mean_val > 0.0: - matrix[i,j] = matrix[i,j] / mean_val + matrix[i, j] = matrix[i, j] / mean_val else: - matrix[i,j] = 0.0 - - elif self.basis.currentText() == 'yz': + matrix[i, j] = 0.0 + + elif text == 'yz': matrix = np.zeros((self.ny, self.nz)) for i in range(self.ny): for j in range(self.nz): - matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (axial_level, i, j))], - self.scoreBox.currentIndex())[score_loc] - # Calculate relative uncertainty from absolute, if - # requested + matrix[i, j] = self.datafile.get_value(tally_id, + spec_list + [('mesh', (axial_level, i + 1, j + 1))], + self.scoreBox.current())[score_loc] + # Calculate relative uncertainty from absolute, if requested if is_mean == 2: # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (axial_level, i, j))], - self.scoreBox.currentIndex())[0] + mean_val = self.datafile.get_value(tally_id, + spec_list + [('mesh', (axial_level, i + 1, j + 1))], + self.scoreBox.current())[0] if mean_val > 0.0: - matrix[i,j] = matrix[i,j] / mean_val + matrix[i, j] = matrix[i, j] / mean_val else: - matrix[i,j] = 0.0 - + matrix[i, j] = 0.0 + else: matrix = np.zeros((self.nx, self.nz)) for i in range(self.nx): for j in range(self.nz): - matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, axial_level, j))], - self.scoreBox.currentIndex())[score_loc] - # Calculate relative uncertainty from absolute, if - # requested + matrix[i, j] = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, axial_level, j + 1))], + self.scoreBox.current())[score_loc] + # Calculate relative uncertainty from absolute, if requested if is_mean == 2: # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, axial_level, j))], - self.scoreBox.currentIndex())[0] + mean_val = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, axial_level, j + 1))], + self.scoreBox.current())[0] if mean_val > 0.0: - matrix[i,j] = matrix[i,j] / mean_val + matrix[i, j] = matrix[i, j] / mean_val else: - matrix[i,j] = 0.0 - -# print spec_list + matrix[i, j] = 0.0 # Clear the figure self.fig.clear() # Make figure, set up color bar self.axes = self.fig.add_subplot(111) - cax = self.axes.imshow(matrix, vmin=0.0, vmax=matrix.max(), - interpolation="nearest") + cax = self.axes.imshow(matrix.transpose(), vmin=0.0, vmax=matrix.max(), + interpolation='none', origin='lower') self.fig.colorbar(cax) self.axes.set_xticks([]) @@ -266,95 +283,43 @@ class AppForm(QMainWindow): # Draw canvas self.canvas.draw() - def _update(self): - '''Updates widget to display new relevant comboboxes and figure data - ''' -# print 'Calling _update...' + def get_file_data(self, filename): + # Create StatePoint object and read in data + self.datafile = StatePoint(filename) + self.datafile.read_results() + self.datafile.generate_stdev() - self.mesh = self.datafile.meshes[ - self.datafile.tallies[ - self.tally.currentIndex()].filters['mesh'].bins[0] - 1] + # Find which tallies are mesh tallies + self.meshTallies = [] + for itally, tally in enumerate(self.datafile.tallies): + if 'mesh' in tally.filters: + self.meshTallies.append(itally) - self.nx, self.ny, self.nz = self.mesh.dimension - - # Clear axial level combobox - self.axial_level.clear() - - # Repopulate axial level combobox based on current basis selection - if (self.basis.currentText() == 'xy'): - self.axial_level.addItems([str(i+1) for i in range(self.nz)]) - elif (self.basis.currentText() == 'yz'): - self.axial_level.addItems([str(i+1) for i in range(self.nx)]) - else: - self.axial_level.addItems([str(i+1) for i in range(self.ny)]) - - # Determine maximum value from current tally data set - self.maxvalue = self.datafile.tallies[ - self.tally.currentIndex()].results.max() -# print self.maxvalue - - # Clear and hide old filter labels - for item in self.labelList: - item.clear() - - # Clear and hide old filter boxes - for j in self.labels: - self.boxes[j].clear() - self.boxes[j].setParent(None) - - self.update() - - def populate_boxes(self): -# print 'Calling populate_boxes...' - - n = 5 - labels = {'cell': 'Cell : ', - 'cellborn': 'Cell born: ', - 'surface': 'Surface: ', - 'material': 'Material: ', - 'universe': 'Universe: '} - - # For each filter in newly-selected tally, name a label and fill the - # relevant combobox with options - for element in self.tally_list[self.tally.currentIndex()]: - nextFilter = self.datafile.tallies[ - self.tally.currentIndex()].filters[element] - if element == 'mesh': - continue - - label = QLabel(self.labels[element]) - self.labelList.append(label) - combobox = self.boxes[element] - self.grid.addWidget(label, n, 0) - self.grid.addWidget(combobox, n, 1) - n += 1 - -# print element - if element in ['cell', 'cellborn', 'surface', 'material', 'universe']: - combobox.addItems([str(i) for i in nextFilter.bins]) -# for i in nextFilter.bins: -# print i - - elif element == 'energyin' or element == 'energyout': - for i in range(nextFilter.length): - text = (str(nextFilter.bins[i]) + ' to ' + - str(nextFilter.bins[i+1])) - combobox.addItem(text) - - self.scoreBox.clear() - for item in self.tally_scores[self.tally.currentIndex()]: - self.scoreBox.addItem(str(item)) - self.grid.addWidget(self.score_label, n, 0) - self.grid.addWidget(self.scoreBox, n, 1) + if not self.meshTallies: + tkMessageBox.showerror("Invalid StatePoint File", + "File does not contain mesh tallies!") + sys.exit(1) +if __name__ == '__main__': + # Hide root window + root = tk.Tk() + root.withdraw() -def main(): - app = QApplication(sys.argv) - form = AppForm() - form.show() - app.exec_() + # If no filename given as command-line argument, open file dialog + if len(sys.argv) < 2: + filename = tkFileDialog.askopenfilename(title='Select statepoint file', + initialdir='.') + else: + filename = sys.argv[1] + if filename: + # Check to make sure file exists + if not os.path.isfile(filename): + tkMessageBox.showerror("File not found", + "Could not find regular file: " + filename) + sys.exit(1) -if __name__ == "__main__": - main() + app = MeshPlotter(root, filename) + root.deiconify() + root.mainloop() diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index ec9d9edf78..430abee8b8 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -1,7 +1,6 @@ #!/usr/bin/env python2 import struct -from math import sqrt from collections import OrderedDict import numpy as np @@ -10,10 +9,10 @@ import scipy.stats filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', 5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'} -score_types = {-1: 'flux', +score_types = {-1: 'flux', -2: 'total', -3: 'scatter', - -4: 'nu-scatter', + -4: 'nu-scatter', -5: 'scatter-n', -6: 'scatter-pn', -7: 'transport', @@ -276,7 +275,7 @@ class StatePoint(object): f.bins = self._get_int(path=base+'bins') else: f.bins = self._get_int(f.length, path=base+'bins') - + base = 'tallies/tally' + str(i+1) + '/' # Read nuclide bins @@ -379,7 +378,7 @@ class StatePoint(object): Calculates the sample mean and standard deviation of the mean for each tally bin. """ - + # Determine number of realizations n = self.n_realizations @@ -387,14 +386,14 @@ class StatePoint(object): for i in range(len(self.global_tallies)): # Get sum and sum of squares s, s2 = self.global_tallies[i] - + # Calculate sample mean and replace value s /= n self.global_tallies[i,0] = s # Calculate standard deviation if s != 0.0: - self.global_tallies[i,1] = t_value*sqrt((s2/n - s*s)/(n-1)) + self.global_tallies[i,1] = t_value*np.sqrt((s2/n - s*s)/(n-1)) # Regular tallies for t in self.tallies: @@ -402,14 +401,14 @@ class StatePoint(object): for j in range(t.results.shape[1]): # Get sum and sum of squares s, s2 = t.results[i,j] - + # Calculate sample mean and replace value s /= n t.results[i,j,0] = s # Calculate standard deviation if s != 0.0: - t.results[i,j,1] = t_value*sqrt((s2/n - s*s)/(n-1)) + t.results[i,j,1] = t_value*np.sqrt((s2/n - s*s)/(n-1)) def get_value(self, tally_index, spec_list, score_index): """Returns a tally score given a list of filters to satisfy. @@ -458,7 +457,7 @@ class StatePoint(object): filter_index += value*t.filters[f_type].stride else: filter_index += f_index*t.filters[f_type].stride - + # Return the desired result from Tally.results. This could be the sum and # sum of squares, or it could be mean and stdev if self.generate_stdev() # has been called already. @@ -531,7 +530,7 @@ class StatePoint(object): for i in range(n_filters): # compute indices for filter combination - filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % + filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1 # append in dictionary bin with filter @@ -544,14 +543,14 @@ class StatePoint(object): dims.reverse() dims = np.asarray(dims) if score_str == 'current': - dims += 1 - meshmax[1:4] = dims + dims += 1 + meshmax[1:4] = dims mesh_bins = np.zeros((n_bins,3)) - mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % + mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1 - mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % + mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1 - mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % + mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1 data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1], mesh_bins[:,2])}) @@ -576,7 +575,7 @@ class StatePoint(object): def _get_data(self, n, typeCode, size): return list(struct.unpack('={0}{1}'.format(n,typeCode), self._f.read(n*size))) - + def _get_int(self, n=1, path=None): if self._hdf5: return [int(v) for v in self._f[path].value] diff --git a/src/utils/track.py b/src/utils/track.py new file mode 100755 index 0000000000..a2dca34e14 --- /dev/null +++ b/src/utils/track.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python2 +"""Convert binary particle track to VTK poly data. + +Usage information can be obtained by running 'track.py --help': + + usage: track.py [-h] [-o OUT] IN [IN ...] + + Convert particle track file to a .pvtp file. + + positional arguments: + IN Input particle track data filename(s). + + optional arguments: + -h, --help show this help message and exit + -o OUT, --out OUT Output VTK poly data filename. + +""" + +import os +import argparse +import struct +import vtk + + +def _parse_args(): + # Create argument parser. + parser = argparse.ArgumentParser( + description='Convert particle track file to a .pvtp file.') + parser.add_argument('input', metavar='IN', type=str, nargs='+', + help='Input particle track data filename(s).') + parser.add_argument('-o', '--out', metavar='OUT', type=str, dest='out', + help='Output VTK poly data filename.') + + # Parse and return commandline arguments. + return parser.parse_args() + + +def main(): + # Parse commandline arguments. + args = _parse_args() + + # Check input file extensions. + for fname in args.input: + if not (fname.endswith('.h5') or fname.endswith('.binary')): + raise ValueError("Input file names must either end with '.h5' or" + "'.binary'.") + + # Make sure that the output filename ends with '.pvtp'. + if not args.out: + args.out = 'tracks.pvtp' + elif os.path.splitext(args.out)[1] != '.pvtp': + args.out = ''.join([args.out, '.pvtp']) + + # Import HDF library if HDF files are present + for fname in args.input: + if fname.endswith('.h5'): + import h5py + break + + # Initialize data arrays and offset. + points = vtk.vtkPoints() + cells = vtk.vtkCellArray() + point_offset = 0 + for fname in args.input: + # Write coordinate values to points array. + if fname.endswith('.binary'): + track = open(fname, 'rb').read() + coords = [struct.unpack("ddd", track[24*i : 24*(i+1)]) + for i in range(len(track)/24)] + n_points = len(coords) + for triplet in coords: + points.InsertNextPoint(triplet) + else: + coords = h5py.File(fname).get('coordinates') + n_points = coords.shape[0] + for i in range(n_points): + points.InsertNextPoint(coords[i,:]) + + # Create VTK line and assign points to line. + line = vtk.vtkPolyLine() + line.GetPointIds().SetNumberOfIds(n_points) + for i in range(n_points): + line.GetPointIds().SetId(i, point_offset+i) + + cells.InsertNextCell(line) + point_offset += n_points + data = vtk.vtkPolyData() + data.SetPoints(points) + data.SetLines(cells) + + writer = vtk.vtkXMLPPolyDataWriter() + writer.SetInput(data) + writer.SetFileName(args.out) + writer.Write() + + + +if __name__ == '__main__': + main() diff --git a/src/xml/dom/FoX_dom.F90 b/src/xml/dom/FoX_dom.F90 index 6e1ac36a7b..4e2fc1e9cc 100644 --- a/src/xml/dom/FoX_dom.F90 +++ b/src/xml/dom/FoX_dom.F90 @@ -74,6 +74,7 @@ module FoX_dom public :: createAttribute public :: createEntityReference public :: getElementsByTagName + public :: getChildrenByTagName public :: getElementById public :: importNode diff --git a/src/xml/dom/m_dom_dom.F90 b/src/xml/dom/m_dom_dom.F90 index 01f6265b07..ec5b788d9b 100644 --- a/src/xml/dom/m_dom_dom.F90 +++ b/src/xml/dom/m_dom_dom.F90 @@ -349,6 +349,7 @@ module m_dom_dom public :: createEntityReference public :: createEmptyEntityReference public :: getElementsByTagName + public :: getChildrenByTagName public :: importNode public :: createElementNS public :: createAttributeNS @@ -6908,6 +6909,166 @@ endif end function getElementsByTagName + function getChildrenByTagName(doc, tagName, name, ex)result(list) + type(DOMException), intent(out), optional :: ex + type(Node), pointer :: doc + character(len=*), intent(in), optional :: tagName, name + type(NodeList), pointer :: list + + type(NodeListPtr), pointer :: nll(:), temp_nll(:) + type(Node), pointer :: arg, this, treeroot + logical :: doneChildren, doneAttributes, allElements + integer :: i, i_tree + + if (.not.associated(doc)) then + if (getFoX_checks().or.FoX_NODE_IS_NULL<200) then + call throw_exception(FoX_NODE_IS_NULL, "getElementsByTagName", ex) + if (present(ex)) then + if (inException(ex)) then + return + endif + endif +endif + + endif + + if (doc%nodeType==DOCUMENT_NODE) then + if (present(name).or..not.present(tagName)) then + if (getFoX_checks().or.FoX_INVALID_NODE<200) then + call throw_exception(FoX_INVALID_NODE, "getElementsByTagName", ex) + if (present(ex)) then + if (inException(ex)) then + return + endif + endif +endif + + endif + elseif (doc%nodeType==ELEMENT_NODE) then + if (present(name).or..not.present(tagName)) then + if (getFoX_checks().or.FoX_INVALID_NODE<200) then + call throw_exception(FoX_INVALID_NODE, "getElementsByTagName", ex) + if (present(ex)) then + if (inException(ex)) then + return + endif + endif +endif + + endif + else + if (getFoX_checks().or.FoX_INVALID_NODE<200) then + call throw_exception(FoX_INVALID_NODE, "getElementsByTagName", ex) + if (present(ex)) then + if (inException(ex)) then + return + endif + endif +endif + + endif + + if (doc%nodeType==DOCUMENT_NODE) then + arg => getDocumentElement(doc) + else + arg => doc + endif + + allocate(list) + allocate(list%nodes(0)) + list%element => doc + if (present(name)) list%nodeName => vs_str_alloc(name) + if (present(tagName)) list%nodeName => vs_str_alloc(tagName) + + allElements = (str_vs(list%nodeName)=="*") + + if (doc%nodeType==DOCUMENT_NODE) then + nll => doc%docExtras%nodelists + elseif (doc%nodeType==ELEMENT_NODE) then + nll => doc%ownerDocument%docExtras%nodelists + endif + allocate(temp_nll(size(nll)+1)) + do i = 1, size(nll) + temp_nll(i)%this => nll(i)%this + enddo + temp_nll(i)%this => list + deallocate(nll) + if (doc%nodeType==DOCUMENT_NODE) then + doc%docExtras%nodelists => temp_nll + elseif (doc%nodeType==ELEMENT_NODE) then + doc%ownerDocument%docExtras%nodelists => temp_nll + endif + + treeroot => arg + + i_tree = 0 + doneChildren = .false. + doneAttributes = .false. + this => treeroot + do + if (.not.doneChildren.and..not.(getNodeType(this)==ELEMENT_NODE.and.doneAttributes)) then + if (this%nodeType==ELEMENT_NODE) then + if ((allElements .or. str_vs(this%nodeName)==tagName) & + .and..not.(getNodeType(doc)==ELEMENT_NODE.and.associated(this, arg))) & + call append(list, this) + doneAttributes = .true. + endif + + else + if (getNodeType(this)==ELEMENT_NODE.and..not.doneChildren) then + doneAttributes = .true. + else + + endif + endif + + + if (.not.doneChildren) then + if (getNodeType(this)==ELEMENT_NODE.and..not.doneAttributes) then + if (getLength(getAttributes(this))>0) then + this => item(getAttributes(this), 0) + else + doneAttributes = .true. + endif + elseif (hasChildNodes(this) .and. .not. associated(getParentNode(this), treeroot)) then + this => getFirstChild(this) + doneChildren = .false. + doneAttributes = .false. + else + doneChildren = .true. + doneAttributes = .false. + endif + + else ! if doneChildren + + if (associated(this, treeroot)) exit + if (getNodeType(this)==ATTRIBUTE_NODE) then + if (i_tree item(getAttributes(getOwnerElement(this)), i_tree) + doneChildren = .false. + else + i_tree= 0 + this => getOwnerElement(this) + doneAttributes = .true. + doneChildren = .false. + endif + elseif (associated(getNextSibling(this))) then + + this => getNextSibling(this) + doneChildren = .false. + doneAttributes = .false. + else + this => getParentNode(this) + endif + endif + + enddo + + + + end function getChildrenByTagName + function importNode(doc , arg, deep , ex)result(np) type(DOMException), intent(out), optional :: ex type(Node), pointer :: doc diff --git a/src/xml_interface.F90 b/src/xml_interface.F90 index 4752c08734..870b4288be 100644 --- a/src/xml_interface.F90 +++ b/src/xml_interface.F90 @@ -72,7 +72,7 @@ contains ! node name. This should only be used for checking a single occurance of a ! sub-element node. To check for sub-element nodes that repeat, use ! get_node_list and get_list_size. This is to minimize number of calls -! to getElementsByTagName. +! to getChildrenByTagName. !=============================================================================== function check_for_node(ptr, node_name) result(found) @@ -94,7 +94,7 @@ contains if (associated(temp_ptr)) return ! Check for a sub-element - elem_list => getElementsByTagName(ptr, trim(node_name)) + elem_list => getChildrenByTagName(ptr, trim(node_name)) ! Get the length of the list if (getLength(elem_list) == 0) then @@ -124,7 +124,7 @@ contains found_ = .false. ! Check for a sub-element - elem_list => getElementsByTagName(in_ptr, trim(node_name)) + elem_list => getChildrenByTagName(in_ptr, trim(node_name)) ! Get the length of the list if (getLength(elem_list) == 0) return @@ -149,7 +149,7 @@ contains type(NodeList), pointer, intent(out) :: out_ptr ! Check for a sub-element - out_ptr => getElementsByTagName(in_ptr, trim(node_name)) + out_ptr => getChildrenByTagName(in_ptr, trim(node_name)) end subroutine get_node_list @@ -525,7 +525,7 @@ contains if (associated(out_ptr)) return ! Check for a sub-element - elem_list => getElementsByTagName(in_ptr, trim(node_name)) + elem_list => getChildrenByTagName(in_ptr, trim(node_name)) ! Get the length of the list if (getLength(elem_list) == 0) then diff --git a/tests/test_salphabeta_multiple/geometry.xml b/tests/test_salphabeta_multiple/geometry.xml index 2d427d8cb3..63f69f7438 100644 --- a/tests/test_salphabeta_multiple/geometry.xml +++ b/tests/test_salphabeta_multiple/geometry.xml @@ -44,9 +44,11 @@ - + - + + + diff --git a/tests/test_salphabeta_multiple/materials.xml b/tests/test_salphabeta_multiple/materials.xml index 97ff56c05f..6e975bde97 100644 --- a/tests/test_salphabeta_multiple/materials.xml +++ b/tests/test_salphabeta_multiple/materials.xml @@ -1,7 +1,8 @@ - + 70c @@ -111,4 +112,19 @@ + + + + + + + + + + + + + + + diff --git a/tests/test_salphabeta_multiple/results_true.dat b/tests/test_salphabeta_multiple/results_true.dat index f26a309aae..bbd61cdb53 100644 --- a/tests/test_salphabeta_multiple/results_true.dat +++ b/tests/test_salphabeta_multiple/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.013747E+00 3.067630E-02 +9.761880E-01 9.170415E-03 diff --git a/tests/test_track_output/geometry.xml b/tests/test_track_output/geometry.xml new file mode 100644 index 0000000000..f3566ab17f --- /dev/null +++ b/tests/test_track_output/geometry.xml @@ -0,0 +1,305 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 + 2 2 2 3 2 2 2 2 2 2 2 3 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 3 2 2 2 2 2 2 2 3 2 2 2 + 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 + 2 2 2 1 2 2 2 2 2 2 2 1 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 3 2 2 2 2 2 2 2 3 2 2 2 + 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + + + + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 + 2 2 2 1 2 2 2 2 2 2 2 1 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 3 2 2 2 2 2 2 2 1 2 2 2 + 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + + + + + + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + + + + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 + 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 + 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 + 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + + + + + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + + + + + + + + + + -12.2682 -12.2682 + 1.63576 1.63576 + + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 + + + + + + + + + + + + + + + + + + + -85.8774 -85.8774 + 24.5364 24.5364 + + 999 999 130 140 150 999 999 + 999 220 230 240 250 260 999 + 130 320 777 222 333 360 150 + 410 240 444 111 666 240 470 + 510 520 888 555 998 560 570 + 999 620 630 240 650 660 999 + 999 999 510 740 570 999 999 + + + + + + + + diff --git a/tests/test_track_output/materials.xml b/tests/test_track_output/materials.xml new file mode 100644 index 0000000000..bbf880d960 --- /dev/null +++ b/tests/test_track_output/materials.xml @@ -0,0 +1,95 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_track_output/settings.xml b/tests/test_track_output/settings.xml new file mode 100644 index 0000000000..ef74341dbb --- /dev/null +++ b/tests/test_track_output/settings.xml @@ -0,0 +1,30 @@ + + + + + + + + + + 2 + 0 + 100 + + + + + + + + + -1 -1 -1 1 1 1 + + + + + 1 1 1 + 1 1 2 + + + diff --git a/tests/test_track_output/test_track_output.py b/tests/test_track_output/test_track_output.py new file mode 100644 index 0000000000..0e7674154c --- /dev/null +++ b/tests/test_track_output/test_track_output.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob + +from nose.plugins.skip import SkipTest + +from nose_mpi import NoseMPI + +pwd = os.path.dirname(__file__) + + +def setup(): + os.putenv('PWD', pwd) + os.chdir(pwd) + + +def test_run(): + openmc_path = pwd + '/../../src/openmc' + if int(NoseMPI.mpi_np) > 0: + proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE) + returncode = proc.wait() + print(proc.communicate()[0]) + assert returncode == 0 + + +def test_created_outputs(): + outputs = [glob.glob(''.join((pwd, '/track_1_1_1.*')))] + outputs.append(glob.glob(''.join((pwd, '/track_1_1_2.*')))) + for files in outputs: + assert len(files) == 1 + assert files[0].endswith('binary') or files[0].endswith('h5') + + +def test_outputs(): + # If vtk python module is not available, we can't run track.py so skip this + # test + try: + import vtk + except ImportError: + raise SkipTest + + call(['../../src/utils/track.py', '-o', 'poly'] + + glob.glob(''.join((pwd, '/track*')))) + poly = ''.join((pwd, '/poly.pvtp')) + assert os.path.isfile(poly) + metric = ''.join((pwd, '/true_poly.pvtp')) + compare = filecmp.cmp(poly, metric) + if not compare: + os.rename('poly.pvtp', 'error_poly.pvtp') + assert compare + + +def teardown(): + temp_files = glob.glob(''.join((pwd, '/statepoint*'))) + temp_files = temp_files + glob.glob(''.join((pwd, '/track*'))) + temp_files = temp_files + glob.glob(''.join((pwd, '/poly*'))) + for f in temp_files: + if os.path.exists(f): + os.remove(f) diff --git a/tests/test_track_output/true_poly.pvtp b/tests/test_track_output/true_poly.pvtp new file mode 100644 index 0000000000..6ded87a0ec --- /dev/null +++ b/tests/test_track_output/true_poly.pvtp @@ -0,0 +1,9 @@ + + + + + + + + +