mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
Merge branch 'develop' into nndc-data
This commit is contained in:
Paul Romano
commit
f3bfaeb2b6
53 changed files with 904 additions and 465 deletions
1
.gitignore
vendored
1
.gitignore
vendored
|
|
@ -19,6 +19,7 @@ src/openmc
|
|||
|
||||
# Documentation builds
|
||||
docs/build
|
||||
docs/source/_images/*.pdf
|
||||
|
||||
# xml-fortran reader
|
||||
src/xml-fortran/xmlreader
|
||||
|
|
|
|||
2
LICENSE
2
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
|
||||
|
|
|
|||
|
|
@ -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 <target>' where <target> 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
|
||||
|
|
|
|||
|
Before Width: | Height: | Size: 278 KiB After Width: | Height: | Size: 278 KiB |
|
Before Width: | Height: | Size: 80 KiB After Width: | Height: | Size: 80 KiB |
|
Before Width: | Height: | Size: 92 KiB After Width: | Height: | Size: 92 KiB |
|
Before Width: | Height: | Size: 460 KiB After Width: | Height: | Size: 460 KiB |
|
Before Width: | Height: | Size: 26 KiB After Width: | Height: | Size: 26 KiB |
|
Before Width: | Height: | Size: 8.4 KiB After Width: | Height: | Size: 8.4 KiB |
|
Before Width: | Height: | Size: 2 KiB After Width: | Height: | Size: 2 KiB |
|
Before Width: | Height: | Size: 6.7 KiB After Width: | Height: | Size: 6.7 KiB |
|
Before Width: | Height: | Size: 34 KiB After Width: | Height: | Size: 34 KiB |
|
Before Width: | Height: | Size: 4 KiB After Width: | Height: | Size: 4 KiB |
|
Before Width: | Height: | Size: 8.7 KiB After Width: | Height: | Size: 8.7 KiB |
|
Before Width: | Height: | Size: 71 KiB After Width: | Height: | Size: 71 KiB |
|
Before Width: | Height: | Size: 17 KiB After Width: | Height: | Size: 17 KiB |
|
Before Width: | Height: | Size: 4.9 KiB After Width: | Height: | Size: 4.9 KiB |
|
Before Width: | Height: | Size: 35 KiB After Width: | Height: | Size: 35 KiB |
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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).
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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<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:
|
||||
|
||||
----------------------------------------------
|
||||
|
|
|
|||
|
|
@ -48,8 +48,8 @@ Publications
|
|||
|
||||
- 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.*
|
||||
(2013). `<http://dx.doi.org/10.1177/1094342013492179>`_
|
||||
*Int. J. High Perform. Comput. Appl.*, **28** (1), 87--96
|
||||
(2014). `<http://dx.doi.org/10.1177/1094342013492179>`_
|
||||
|
||||
- Paul K. Romano, Andrew R. Siegel, Benoit Forget, and Kord Smith, "Data
|
||||
decomposition of Monte Carlo particle transport simulations via tally
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <http://www.riverbankcomputing.com/software/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
|
||||
|
|
@ -357,7 +357,7 @@ and dumped to MATLAB in one step.
|
|||
Particle Track Visualization
|
||||
----------------------------
|
||||
|
||||
.. image:: ../../img/Tracks.png
|
||||
.. image:: ../_images/Tracks.png
|
||||
:height: 200px
|
||||
|
||||
OpenMC can dump particle tracks—the position of particles as they are
|
||||
|
|
|
|||
|
|
@ -46,7 +46,7 @@ to locate ACE format cross section libraries if the user has not specified the
|
|||
<cross_sections> 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
|
||||
|
|
|
|||
10
src/Makefile
10
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
|
||||
|
||||
#===============================================================================
|
||||
|
|
|
|||
83
src/ace.F90
83
src/ace.F90
|
|
@ -1195,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)))
|
||||
|
|
@ -1221,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)
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ 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()
|
||||
|
||||
|
|
@ -143,6 +143,20 @@ module ace_header
|
|||
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
|
||||
|
|
@ -163,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(:)
|
||||
! 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)
|
||||
|
|
|
|||
|
|
@ -577,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
|
||||
|
|
|
|||
|
|
@ -155,7 +155,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 :: &
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -821,7 +821,7 @@ contains
|
|||
|
||||
end subroutine replay_batch_history
|
||||
|
||||
#ifdef OPENMP
|
||||
#ifdef _OPENMP
|
||||
!===============================================================================
|
||||
! JOIN_BANK_FROM_THREADS
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -205,7 +205,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
|
||||
|
|
@ -438,7 +438,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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@ module initialize
|
|||
use mpi
|
||||
#endif
|
||||
|
||||
#ifdef OPENMP
|
||||
#ifdef _OPENMP
|
||||
use omp_lib
|
||||
#endif
|
||||
|
||||
|
|
@ -372,7 +372,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
|
||||
|
|
@ -830,14 +830,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
|
||||
|
|
|
|||
|
|
@ -228,7 +228,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
|
||||
|
|
@ -474,7 +474,8 @@ contains
|
|||
allocate(temp_int_array(n_tracks))
|
||||
call get_node_array(doc, "track", temp_int_array)
|
||||
|
||||
! Reshape into track_identifiers -- note automatic array allocation
|
||||
! Reshape into track_identifiers
|
||||
allocate(track_identifiers(3, n_tracks/3))
|
||||
track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
|
||||
end if
|
||||
|
||||
|
|
@ -2938,7 +2939,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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
225
src/physics.F90
225
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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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" } }? &
|
||||
|
|
|
|||
|
|
@ -97,6 +97,7 @@ contains
|
|||
|
||||
! 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)
|
||||
|
|
|
|||
|
|
@ -464,6 +464,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)
|
||||
|
|
@ -542,10 +543,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", &
|
||||
|
|
|
|||
|
|
@ -56,6 +56,7 @@ contains
|
|||
subroutine finalize_particle_track(p)
|
||||
type(Particle), intent(in) :: p
|
||||
|
||||
integer :: length(2)
|
||||
character(MAX_FILE_LEN) :: fname
|
||||
type(BinaryOutput) :: binout
|
||||
|
||||
|
|
@ -69,7 +70,8 @@ contains
|
|||
// '.binary'
|
||||
#endif
|
||||
call binout % file_create(fname)
|
||||
call binout % write_data(coords, 'coordinates', length=(/3, n_tracks/))
|
||||
length = [3, n_tracks]
|
||||
call binout % write_data(coords, 'coordinates', length=length)
|
||||
call binout % file_close()
|
||||
deallocate(coords)
|
||||
end subroutine finalize_particle_track
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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('<<ComboboxSelected>>', 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('<<ComboboxSelected>>', 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('<<ComboboxSelected>>', 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('<<ComboboxSelected>>', 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('<<ComboboxSelected>>', 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('<<ComboboxSelected>>', 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()
|
||||
|
|
|
|||
|
|
@ -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]
|
||||
|
|
|
|||
|
|
@ -74,6 +74,7 @@ module FoX_dom
|
|||
public :: createAttribute
|
||||
public :: createEntityReference
|
||||
public :: getElementsByTagName
|
||||
public :: getChildrenByTagName
|
||||
public :: getElementById
|
||||
|
||||
public :: importNode
|
||||
|
|
|
|||
|
|
@ -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<getLength(getAttributes(getOwnerElement(this)))-1) then
|
||||
i_tree= i_tree+ 1
|
||||
this => 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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue