mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge branch 'develop' into cmake
Conflicts: src/Makefile
This commit is contained in:
Paul Romano
commit
818f8d31dd
53 changed files with 892 additions and 453 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
|
||||
|
|
|
|||
|
|
@ -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 |
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
||||
----------------------------------------------
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -329,6 +329,7 @@ solver_interface.o: vector_header.o
|
|||
source.o: bank_header.o
|
||||
source.o: constants.o
|
||||
source.o: error.o
|
||||
source.o: geometry.o
|
||||
source.o: geometry_header.o
|
||||
source.o: global.o
|
||||
source.o: math.o
|
||||
|
|
|
|||
113
src/ace.F90
113
src/ace.F90
|
|
@ -115,9 +115,9 @@ contains
|
|||
! search through the list of nuclides for one which has a matching zaid
|
||||
sab => sab_tables(mat % i_sab_tables(k))
|
||||
|
||||
! Loop through nuclides and find match
|
||||
! Loop through nuclides and find match
|
||||
FIND_NUCLIDE: do j = 1, mat % n_nuclides
|
||||
if (nuclides(mat % nuclide(j)) % zaid == sab % zaid) then
|
||||
if (any(sab % zaid == nuclides(mat % nuclide(j)) % zaid)) then
|
||||
mat % i_sab_nuclides(k) = j
|
||||
exit FIND_NUCLIDE
|
||||
end if
|
||||
|
|
@ -161,16 +161,16 @@ contains
|
|||
mat % i_sab_tables(m) = temp_table
|
||||
end do SORT_SAB
|
||||
end if
|
||||
|
||||
|
||||
! Deallocate temporary arrays for names of nuclides and S(a,b) tables
|
||||
if (allocated(mat % names)) deallocate(mat % names)
|
||||
if (allocated(mat % sab_names)) deallocate(mat % sab_names)
|
||||
|
||||
end do MATERIAL_LOOP2
|
||||
|
||||
|
||||
! Avoid some valgrind leak errors
|
||||
call already_read % clear()
|
||||
|
||||
|
||||
end subroutine read_xs
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -253,7 +253,7 @@ contains
|
|||
end if
|
||||
|
||||
! Read more header and NXS and JXS
|
||||
read(UNIT=in, FMT=100) comment, mat, &
|
||||
read(UNIT=in, FMT=100) comment, mat, &
|
||||
(zaids(i), awrs(i), i=1,16), NXS, JXS
|
||||
100 format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/&
|
||||
,8I9/8I9/8I9/8I9/8I9/8I9)
|
||||
|
|
@ -277,7 +277,7 @@ contains
|
|||
ACCESS='direct', RECL=record_length)
|
||||
|
||||
! Read all header information
|
||||
read(UNIT=in, REC=location) name, awr, kT, date_, &
|
||||
read(UNIT=in, REC=location) name, awr, kT, date_, &
|
||||
comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS
|
||||
|
||||
! determine table length
|
||||
|
|
@ -333,7 +333,15 @@ contains
|
|||
sab % name = name
|
||||
sab % awr = awr
|
||||
sab % kT = kT
|
||||
sab % zaid = zaids(1)
|
||||
! Find sab % n_zaid
|
||||
do i = 1, 16
|
||||
if (zaids(i) == 0) then
|
||||
sab % n_zaid = i - 1
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
allocate(sab % zaid(sab % n_zaid))
|
||||
sab % zaid = zaids(1: sab % n_zaid)
|
||||
|
||||
call read_thermal_data(sab)
|
||||
end select
|
||||
|
|
@ -406,7 +414,7 @@ contains
|
|||
integer :: LNU ! type of nu data (polynomial or tabular)
|
||||
integer :: NC ! number of polynomial coefficients
|
||||
integer :: NR ! number of interpolation regions
|
||||
integer :: NE ! number of energies
|
||||
integer :: NE ! number of energies
|
||||
integer :: NPCR ! number of delayed neutron precursor groups
|
||||
integer :: LED ! location of energy distribution locators
|
||||
integer :: LDIS ! location of all energy distributions
|
||||
|
|
@ -809,7 +817,7 @@ contains
|
|||
|
||||
LED = JXS(10)
|
||||
|
||||
! Loop over all reactions
|
||||
! Loop over all reactions
|
||||
do i = 1, NXS(5)
|
||||
rxn => nuc % reactions(i+1) ! skip over elastic scattering
|
||||
rxn % has_energy_dist = .true.
|
||||
|
|
@ -840,7 +848,7 @@ contains
|
|||
|
||||
integer :: LDIS ! location of all energy distributions
|
||||
integer :: LNW ! location of next energy distribution if multiple
|
||||
integer :: LAW ! secondary energy distribution law
|
||||
integer :: LAW ! secondary energy distribution law
|
||||
integer :: NR ! number of interpolation regions
|
||||
integer :: NE ! number of incoming energies
|
||||
integer :: IDAT ! location of first energy distribution for given MT
|
||||
|
|
@ -1187,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)))
|
||||
|
|
@ -1213,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)
|
||||
|
|
|
|||
|
|
@ -16,7 +16,7 @@ module ace_header
|
|||
integer, allocatable :: type(:) ! type of distribution
|
||||
integer, allocatable :: location(:) ! location of each table
|
||||
real(8), allocatable :: data(:) ! angular distribution data
|
||||
|
||||
|
||||
! Type-Bound procedures
|
||||
contains
|
||||
procedure :: clear => distangle_clear ! Deallocates DistAngle
|
||||
|
|
@ -32,10 +32,10 @@ module ace_header
|
|||
type(Tab1) :: p_valid ! probability of law validity
|
||||
real(8), allocatable :: data(:) ! energy distribution data
|
||||
|
||||
! For reactions that may have multiple energy distributions such as (n.2n),
|
||||
! For reactions that may have multiple energy distributions such as (n,2n),
|
||||
! this pointer allows multiple laws to be stored
|
||||
type(DistEnergy), pointer :: next => null()
|
||||
|
||||
|
||||
! Type-Bound procedures
|
||||
contains
|
||||
procedure :: clear => distenergy_clear ! Deallocates DistEnergy
|
||||
|
|
@ -57,7 +57,7 @@ module ace_header
|
|||
logical :: has_energy_dist ! Energy distribution present?
|
||||
type(DistAngle) :: adist ! Secondary angular distribution
|
||||
type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution
|
||||
|
||||
|
||||
! Type-Bound procedures
|
||||
contains
|
||||
procedure :: clear => reaction_clear ! Deallocates Reaction
|
||||
|
|
@ -76,7 +76,7 @@ module ace_header
|
|||
logical :: multiply_smooth ! multiply by smooth cross section?
|
||||
real(8), allocatable :: energy(:) ! incident energies
|
||||
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
|
||||
|
||||
|
||||
! Type-Bound procedures
|
||||
contains
|
||||
procedure :: clear => urrdata_clear ! Deallocates UrrData
|
||||
|
|
@ -137,22 +137,37 @@ module ace_header
|
|||
! Reactions
|
||||
integer :: n_reaction ! # of reactions
|
||||
type(Reaction), pointer :: reactions(:) => null()
|
||||
|
||||
|
||||
! Type-Bound procedures
|
||||
contains
|
||||
procedure :: clear => nuclide_clear ! Deallocates Nuclide
|
||||
end type Nuclide
|
||||
|
||||
!===============================================================================
|
||||
! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
|
||||
! thermal scattering collisions which utilize a continuous secondary energy
|
||||
! representation.
|
||||
!===============================================================================
|
||||
|
||||
type DistEnergySab
|
||||
integer :: n_e_out
|
||||
real(8), allocatable :: e_out(:)
|
||||
real(8), allocatable :: e_out_pdf(:)
|
||||
real(8), allocatable :: e_out_cdf(:)
|
||||
real(8), allocatable :: mu(:,:)
|
||||
end type DistEnergySab
|
||||
|
||||
!===============================================================================
|
||||
! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
|
||||
! of light isotopes such as water, graphite, Be, etc
|
||||
!===============================================================================
|
||||
|
||||
|
||||
type SAlphaBeta
|
||||
character(10) :: name ! name of table, e.g. lwtr.10t
|
||||
integer :: zaid ! Z and A identifier, e.g. 6012 for Carbon-12
|
||||
real(8) :: awr ! weight of nucleus in neutron masses
|
||||
real(8) :: kT ! temperature in MeV (k*T)
|
||||
character(10) :: name ! name of table, e.g. lwtr.10t
|
||||
real(8) :: awr ! weight of nucleus in neutron masses
|
||||
real(8) :: kT ! temperature in MeV (k*T)
|
||||
integer :: n_zaid ! Number of valid zaids
|
||||
integer, allocatable :: zaid(:) ! List of valid Z and A identifiers, e.g. 6012
|
||||
|
||||
! threshold for S(a,b) treatment (usually ~4 eV)
|
||||
real(8) :: threshold_inelastic
|
||||
|
|
@ -162,11 +177,17 @@ module ace_header
|
|||
integer :: n_inelastic_e_in ! # of incoming E for inelastic
|
||||
integer :: n_inelastic_e_out ! # of outgoing E for inelastic
|
||||
integer :: n_inelastic_mu ! # of outgoing angles for inelastic
|
||||
integer :: secondary_mode ! secondary mode (equal/skewed)
|
||||
integer :: secondary_mode ! secondary mode (equal/skewed/continuous)
|
||||
real(8), allocatable :: inelastic_e_in(:)
|
||||
real(8), allocatable :: inelastic_sigma(:)
|
||||
real(8), allocatable :: inelastic_sigma(:)
|
||||
! The following are used only if secondary_mode is 0 or 1
|
||||
real(8), allocatable :: inelastic_e_out(:,:)
|
||||
real(8), allocatable :: inelastic_mu(:,:,:)
|
||||
! The following is used only if secondary_mode is 3
|
||||
! The different implementation is necessary because the continuous
|
||||
! representation has a variable number of outgoing energy points for each
|
||||
! incoming energy
|
||||
type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein
|
||||
|
||||
! Elastic scattering data
|
||||
integer :: elastic_mode ! elastic mode (discrete/exact)
|
||||
|
|
@ -240,125 +261,125 @@ module ace_header
|
|||
|
||||
!===============================================================================
|
||||
! DISTANGLE_CLEAR resets and deallocates data in Reaction.
|
||||
!===============================================================================
|
||||
|
||||
!===============================================================================
|
||||
|
||||
subroutine distangle_clear(this)
|
||||
|
||||
|
||||
class(DistAngle), intent(inout) :: this ! The DistAngle object to clear
|
||||
|
||||
|
||||
if (allocated(this % energy)) &
|
||||
deallocate(this % energy, this % type, this % location, this % data)
|
||||
|
||||
end subroutine distangle_clear
|
||||
|
||||
end subroutine distangle_clear
|
||||
|
||||
!===============================================================================
|
||||
! DISTENERGY_CLEAR resets and deallocates data in DistEnergy.
|
||||
!===============================================================================
|
||||
|
||||
!===============================================================================
|
||||
|
||||
recursive subroutine distenergy_clear(this)
|
||||
|
||||
|
||||
class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear
|
||||
|
||||
|
||||
! Clear p_valid
|
||||
call this % p_valid % clear()
|
||||
|
||||
|
||||
if (allocated(this % data)) &
|
||||
deallocate(this % data)
|
||||
|
||||
|
||||
if (associated(this % next)) then
|
||||
! recursively clear this item
|
||||
call this % next % clear()
|
||||
deallocate(this % next)
|
||||
end if
|
||||
|
||||
|
||||
end subroutine distenergy_clear
|
||||
|
||||
|
||||
!===============================================================================
|
||||
! REACTION_CLEAR resets and deallocates data in Reaction.
|
||||
!===============================================================================
|
||||
|
||||
!===============================================================================
|
||||
|
||||
subroutine reaction_clear(this)
|
||||
|
||||
|
||||
class(Reaction), intent(inout) :: this ! The Reaction object to clear
|
||||
|
||||
|
||||
if (allocated(this % sigma)) &
|
||||
deallocate(this % sigma)
|
||||
|
||||
|
||||
if (associated(this % edist)) then
|
||||
call this % edist % clear()
|
||||
deallocate(this % edist)
|
||||
end if
|
||||
|
||||
|
||||
call this % adist % clear()
|
||||
|
||||
end subroutine reaction_clear
|
||||
|
||||
|
||||
end subroutine reaction_clear
|
||||
|
||||
!===============================================================================
|
||||
! URRDATA_CLEAR resets and deallocates data in Reaction.
|
||||
!===============================================================================
|
||||
|
||||
!===============================================================================
|
||||
|
||||
subroutine urrdata_clear(this)
|
||||
|
||||
|
||||
class(UrrData), intent(inout) :: this ! The UrrData object to clear
|
||||
|
||||
|
||||
if (allocated(this % energy)) &
|
||||
deallocate(this % energy, this % prob)
|
||||
|
||||
end subroutine urrdata_clear
|
||||
|
||||
end subroutine urrdata_clear
|
||||
|
||||
!===============================================================================
|
||||
! NUCLIDE_CLEAR resets and deallocates data in Nuclide.
|
||||
!===============================================================================
|
||||
|
||||
!===============================================================================
|
||||
|
||||
subroutine nuclide_clear(this)
|
||||
|
||||
|
||||
class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
|
||||
|
||||
|
||||
integer :: i ! Loop counter
|
||||
|
||||
|
||||
if (allocated(this % grid_index)) &
|
||||
deallocate(this % grid_index)
|
||||
|
||||
|
||||
if (allocated(this % energy)) &
|
||||
deallocate(this % total, this % elastic, this % fission, &
|
||||
this % nu_fission, this % absorption)
|
||||
if (allocated(this % heating)) &
|
||||
deallocate(this % heating)
|
||||
|
||||
|
||||
if (allocated(this % index_fission)) &
|
||||
deallocate(this % index_fission)
|
||||
|
||||
|
||||
if (allocated(this % nu_t_data)) &
|
||||
deallocate(this % nu_t_data)
|
||||
|
||||
|
||||
if (allocated(this % nu_p_data)) &
|
||||
deallocate(this % nu_p_data)
|
||||
|
||||
|
||||
if (allocated(this % nu_d_data)) &
|
||||
deallocate(this % nu_d_data)
|
||||
|
||||
|
||||
if (allocated(this % nu_d_precursor_data)) &
|
||||
deallocate(this % nu_d_precursor_data)
|
||||
|
||||
|
||||
if (associated(this % nu_d_edist)) then
|
||||
do i = 1, size(this % nu_d_edist)
|
||||
call this % nu_d_edist(i) % clear()
|
||||
end do
|
||||
deallocate(this % nu_d_edist)
|
||||
end if
|
||||
|
||||
|
||||
if (associated(this % urr_data)) then
|
||||
call this % urr_data % clear()
|
||||
deallocate(this % urr_data)
|
||||
end if
|
||||
|
||||
|
||||
if (associated(this % reactions)) then
|
||||
do i = 1, size(this % reactions)
|
||||
call this % reactions(i) % clear()
|
||||
end do
|
||||
deallocate(this % reactions)
|
||||
end if
|
||||
|
||||
end subroutine nuclide_clear
|
||||
|
||||
end subroutine nuclide_clear
|
||||
|
||||
end module ace_header
|
||||
|
|
|
|||
|
|
@ -934,7 +934,7 @@ contains
|
|||
subroutine compute_effective_downscatter()
|
||||
|
||||
use constants, only: ZERO, CMFD_NOACCEL
|
||||
use global, only: cmfd, cmfd_downscatter
|
||||
use global, only: cmfd
|
||||
|
||||
integer :: nx ! number of mesh cells in x direction
|
||||
integer :: ny ! number of mesh cells in y direction
|
||||
|
|
|
|||
|
|
@ -117,9 +117,9 @@ contains
|
|||
|
||||
subroutine process_cmfd_options()
|
||||
|
||||
#ifdef PETSC
|
||||
use global, only: cmfd_snes_monitor, cmfd_ksp_monitor, mpi_err
|
||||
|
||||
#ifdef PETSC
|
||||
! Check for snes monitor
|
||||
if (cmfd_snes_monitor) call PetscOptionsSetValue("-snes_monitor", &
|
||||
"stdout", mpi_err)
|
||||
|
|
@ -138,9 +138,10 @@ contains
|
|||
subroutine calc_fission_source()
|
||||
|
||||
use constants, only: CMFD_NOACCEL, ZERO, TWO
|
||||
use global, only: cmfd, cmfd_coremap, master, mpi_err, entropy_on, &
|
||||
current_batch
|
||||
use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch
|
||||
|
||||
#ifdef MPI
|
||||
use global, only: mpi_err
|
||||
use mpi
|
||||
#endif
|
||||
|
||||
|
|
@ -261,13 +262,14 @@ contains
|
|||
|
||||
use constants, only: ZERO, ONE
|
||||
use error, only: warning, fatal_error
|
||||
use global, only: n_particles, meshes, source_bank, work, &
|
||||
n_user_meshes, message, cmfd, master, mpi_err
|
||||
use global, only: meshes, source_bank, work, n_user_meshes, message, &
|
||||
cmfd, master
|
||||
use mesh_header, only: StructuredMesh
|
||||
use mesh, only: count_bank_sites, get_mesh_indices
|
||||
use search, only: binary_search
|
||||
|
||||
#ifdef MPI
|
||||
use global, only: mpi_err
|
||||
use mpi
|
||||
#endif
|
||||
|
||||
|
|
|
|||
|
|
@ -20,7 +20,9 @@ contains
|
|||
|
||||
use cmfd_header, only: allocate_cmfd
|
||||
|
||||
#ifdef PETSC
|
||||
integer :: new_comm ! new mpi communicator
|
||||
#endif
|
||||
integer :: color ! color group of processor
|
||||
|
||||
! Read in cmfd input file
|
||||
|
|
@ -34,17 +36,17 @@ contains
|
|||
end if
|
||||
|
||||
! Split up procs
|
||||
# ifdef PETSC
|
||||
#ifdef PETSC
|
||||
call MPI_COMM_SPLIT(MPI_COMM_WORLD, color, 0, new_comm, mpi_err)
|
||||
# endif
|
||||
#endif
|
||||
|
||||
! assign to PETSc
|
||||
# ifdef PETSC
|
||||
#ifdef PETSC
|
||||
PETSC_COMM_WORLD = new_comm
|
||||
|
||||
! Initialize PETSc on all procs
|
||||
call PetscInitialize(PETSC_NULL_CHARACTER, mpi_err)
|
||||
# endif
|
||||
#endif
|
||||
|
||||
! Initialize timers
|
||||
call time_cmfd % reset()
|
||||
|
|
@ -575,7 +577,9 @@ contains
|
|||
end do
|
||||
|
||||
! Put cmfd tallies into active tally array and turn tallies on
|
||||
!$omp parallel
|
||||
call setup_active_cmfdtallies()
|
||||
!$omp end parallel
|
||||
tallies_on = .true.
|
||||
|
||||
end subroutine create_cmfd_tally
|
||||
|
|
|
|||
|
|
@ -102,7 +102,9 @@ contains
|
|||
subroutine init_data(adjoint)
|
||||
|
||||
use constants, only: ONE, ZERO
|
||||
#ifdef PETSC
|
||||
use global, only: cmfd_write_matrices
|
||||
#endif
|
||||
|
||||
logical, intent(in) :: adjoint ! adjoint calcualtion
|
||||
|
||||
|
|
|
|||
|
|
@ -20,7 +20,7 @@ module constants
|
|||
FILETYPE_PARTICLE_RESTART = -2
|
||||
|
||||
! ============================================================================
|
||||
! ADJUSTABLE PARAMETERS
|
||||
! ADJUSTABLE PARAMETERS
|
||||
|
||||
! NOTE: This is the only section of the constants module that should ever be
|
||||
! adjusted. Modifying constants in other sections may cause the code to fail.
|
||||
|
|
@ -50,6 +50,10 @@ module constants
|
|||
integer, parameter :: MAX_WORD_LEN = 150
|
||||
integer, parameter :: MAX_FILE_LEN = 255
|
||||
|
||||
! Maximum number of external source spatial resamples to encounter before an
|
||||
! error is thrown.
|
||||
integer, parameter :: MAX_EXTSRC_RESAMPLES = 10000
|
||||
|
||||
! ============================================================================
|
||||
! PHYSICAL CONSTANTS
|
||||
|
||||
|
|
@ -110,9 +114,9 @@ module constants
|
|||
|
||||
! Surface types
|
||||
integer, parameter :: &
|
||||
SURF_PX = 1, & ! Plane parallel to x-plane
|
||||
SURF_PY = 2, & ! Plane parallel to y-plane
|
||||
SURF_PZ = 3, & ! Plane parallel to z-plane
|
||||
SURF_PX = 1, & ! Plane parallel to x-plane
|
||||
SURF_PY = 2, & ! Plane parallel to y-plane
|
||||
SURF_PZ = 3, & ! Plane parallel to z-plane
|
||||
SURF_PLANE = 4, & ! Arbitrary plane
|
||||
SURF_CYL_X = 5, & ! Cylinder along x-axis
|
||||
SURF_CYL_Y = 6, & ! Cylinder along y-axis
|
||||
|
|
@ -143,7 +147,7 @@ module constants
|
|||
ELECTRON = 3
|
||||
|
||||
! Angular distribution type
|
||||
integer, parameter :: &
|
||||
integer, parameter :: &
|
||||
ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
|
||||
ANGLE_32_EQUI = 2, & ! 32 equiprobable bins
|
||||
ANGLE_TABULAR = 3 ! Tabular angular distribution
|
||||
|
|
@ -151,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 :: &
|
||||
|
|
@ -275,7 +280,7 @@ module constants
|
|||
SCORE_KAPPA_FISSION = -12, & ! fission energy production rate
|
||||
SCORE_CURRENT = -13, & ! partial current
|
||||
SCORE_EVENTS = -14 ! number of events
|
||||
|
||||
|
||||
! Maximum scattering order supported
|
||||
integer, parameter :: SCATT_ORDER_MAX = 10
|
||||
character(len=*), parameter :: SCATT_ORDER_MAX_PNSTR = "scatter-p10"
|
||||
|
|
@ -322,7 +327,7 @@ module constants
|
|||
|
||||
! Source angular distribution types
|
||||
integer, parameter :: &
|
||||
SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular
|
||||
SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular
|
||||
SRC_ANGLE_MONO = 2, & ! Monodirectional source
|
||||
SRC_ANGLE_TABULAR = 3 ! Tabular distribution
|
||||
|
||||
|
|
@ -332,7 +337,7 @@ module constants
|
|||
SRC_ENERGY_MAXWELL = 2, & ! Maxwell fission spectrum
|
||||
SRC_ENERGY_WATT = 3, & ! Watt fission spectrum
|
||||
SRC_ENERGY_TABULAR = 4 ! Tabular distribution
|
||||
|
||||
|
||||
! ============================================================================
|
||||
! MISCELLANEOUS CONSTANTS
|
||||
|
||||
|
|
|
|||
|
|
@ -87,7 +87,7 @@ contains
|
|||
! =======================================================================
|
||||
! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
|
||||
|
||||
if (k == 1 .and. a >= 4.0) then
|
||||
if (k == 1 .and. a >= -4.0) then
|
||||
! Since x = 0, this implies that a = -y
|
||||
F_b = F_a
|
||||
a = -y
|
||||
|
|
@ -101,7 +101,7 @@ contains
|
|||
end if
|
||||
|
||||
! =======================================================================
|
||||
! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
|
||||
! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
|
||||
|
||||
k = i
|
||||
b = ZERO
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
@ -635,8 +636,8 @@ contains
|
|||
if (check_for_node(node_sp, "source_write")) then
|
||||
call get_node_value(node_sp, "source_write", temp_str)
|
||||
call lower_case(temp_str)
|
||||
if (trim(temp_str) == 'true' .or. &
|
||||
trim(temp_str) == '1') source_separate = .true.
|
||||
if (trim(temp_str) == 'false' .or. &
|
||||
trim(temp_str) == '0') source_write = .false.
|
||||
end if
|
||||
else
|
||||
! If no <state_point> tag was present, by default write state point at
|
||||
|
|
@ -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
|
||||
|
|
|
|||
111
src/output.F90
111
src/output.F90
|
|
@ -31,6 +31,10 @@ contains
|
|||
|
||||
subroutine title()
|
||||
|
||||
#ifdef _OPENMP
|
||||
use omp_lib
|
||||
#endif
|
||||
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
|
||||
' .d88888b. 888b d888 .d8888b.', &
|
||||
' d88P" "Y88b 8888b d8888 d88P Y88b', &
|
||||
|
|
@ -46,25 +50,31 @@ contains
|
|||
|
||||
! Write version information
|
||||
write(UNIT=OUTPUT_UNIT, FMT=*) &
|
||||
' Copyright: 2011-2013 Massachusetts Institute of Technology'
|
||||
' Copyright: 2011-2013 Massachusetts Institute of Technology'
|
||||
write(UNIT=OUTPUT_UNIT, FMT=*) &
|
||||
' License: http://mit-crpg.github.io/openmc/license.html'
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') &
|
||||
' License: http://mit-crpg.github.io/openmc/license.html'
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') &
|
||||
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
|
||||
#ifdef GIT_SHA1
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",7X,A)') GIT_SHA1
|
||||
#endif
|
||||
|
||||
! Write the date and time
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') &
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",6X,A)') &
|
||||
time_stamp()
|
||||
|
||||
#ifdef MPI
|
||||
! Write number of processors
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') &
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",2X,A)') &
|
||||
trim(to_str(n_procs))
|
||||
#endif
|
||||
|
||||
#ifdef _OPENMP
|
||||
! Write number of OpenMP threads
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') &
|
||||
trim(to_str(omp_get_max_threads()))
|
||||
#endif
|
||||
|
||||
end subroutine title
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -126,7 +136,7 @@ contains
|
|||
! print header based on level
|
||||
select case (header_level)
|
||||
case (1)
|
||||
write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), &
|
||||
write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), &
|
||||
repeat('=', n) // '> ' // trim(line) // ' <' // &
|
||||
repeat('=', m), repeat('=', 75)
|
||||
case (2)
|
||||
|
|
@ -180,7 +190,7 @@ contains
|
|||
end subroutine print_usage
|
||||
|
||||
!===============================================================================
|
||||
! WRITE_MESSAGE displays an informational message to the log file and the
|
||||
! WRITE_MESSAGE displays an informational message to the log file and the
|
||||
! standard output stream.
|
||||
!===============================================================================
|
||||
|
||||
|
|
@ -215,7 +225,7 @@ contains
|
|||
! Determine last space in current line
|
||||
last_space = index(message(i_start+1:i_start+line_wrap), &
|
||||
' ', BACK=.true.)
|
||||
if (last_space == 0) then
|
||||
if (last_space == 0) then
|
||||
i_end = min(length + 1, i_start+line_wrap) - 1
|
||||
write(ou, fmt='(1X,A)') message(i_start+1:i_end)
|
||||
else
|
||||
|
|
@ -1028,8 +1038,10 @@ contains
|
|||
type(SAlphaBeta), pointer :: sab
|
||||
integer, optional :: unit
|
||||
|
||||
integer :: size_sab ! memory used by S(a,b) table
|
||||
integer :: unit_ ! unit to write to
|
||||
integer :: size_sab ! memory used by S(a,b) table
|
||||
integer :: unit_ ! unit to write to
|
||||
integer :: i ! Loop counter for parsing through sab % zaid
|
||||
integer :: char_count ! Counter for the number of characters on a line
|
||||
|
||||
! set default unit for writing information
|
||||
if (present(unit)) then
|
||||
|
|
@ -1040,7 +1052,30 @@ contains
|
|||
|
||||
! Basic S(a,b) table information
|
||||
write(unit_,*) 'S(a,b) Table ' // trim(sab % name)
|
||||
write(unit_,*) ' zaid = ' // trim(to_str(sab % zaid))
|
||||
write(unit_,'(A)',advance="no") ' zaids = '
|
||||
! Initialize the counter based on the above string
|
||||
char_count = 11
|
||||
do i = 1, sab % n_zaid
|
||||
! Deal with a line thats too long
|
||||
if (char_count >= 73) then ! 73 = 80 - (5 ZAID chars + 1 space + 1 comma)
|
||||
! End the line
|
||||
write(unit_,*) ""
|
||||
! Add 11 leading blanks
|
||||
write(unit_,'(A)', advance="no") " "
|
||||
! reset the counter to 11
|
||||
char_count = 11
|
||||
end if
|
||||
if (i < sab % n_zaid) then
|
||||
! Include a comma
|
||||
write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) // ", "
|
||||
char_count = char_count + len(trim(to_str(sab % zaid(i)))) + 2
|
||||
else
|
||||
! Don't include a comma, since we are all done
|
||||
write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i)))
|
||||
end if
|
||||
|
||||
end do
|
||||
write(unit_,*) "" ! Move to next line
|
||||
write(unit_,*) ' awr = ' // trim(to_str(sab % awr))
|
||||
write(unit_,*) ' kT = ' // trim(to_str(sab % kT))
|
||||
|
||||
|
|
@ -1242,7 +1277,7 @@ contains
|
|||
end select
|
||||
end if
|
||||
write(UNIT=ou, FMT=*)
|
||||
|
||||
|
||||
write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "========="
|
||||
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
|
||||
if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
|
||||
|
|
@ -1272,7 +1307,7 @@ contains
|
|||
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
|
||||
entropy(overall_gen)
|
||||
|
||||
if (overall_gen - n_inactive*gen_per_batch > 1) then
|
||||
if (overall_gen - n_inactive*gen_per_batch > 1) then
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
|
||||
keff, keff_std
|
||||
end if
|
||||
|
|
@ -1300,7 +1335,7 @@ contains
|
|||
entropy(current_batch*gen_per_batch)
|
||||
|
||||
! write out accumulated k-effective if after first active batch
|
||||
if (overall_gen - n_inactive*gen_per_batch > 1) then
|
||||
if (overall_gen - n_inactive*gen_per_batch > 1) then
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
|
||||
keff, keff_std
|
||||
else
|
||||
|
|
@ -1310,7 +1345,7 @@ contains
|
|||
! write out cmfd keff if it is active and other display info
|
||||
if (cmfd_on) then
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
|
||||
cmfd % k_cmfd(current_batch)
|
||||
cmfd % k_cmfd(current_batch)
|
||||
select case(trim(cmfd_display))
|
||||
case('entropy')
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
|
||||
|
|
@ -1349,7 +1384,7 @@ contains
|
|||
|
||||
! Plot id
|
||||
write(ou,100) "Plot ID:", trim(to_str(pl % id))
|
||||
|
||||
|
||||
! Plot type
|
||||
if (pl % type == PLOT_TYPE_SLICE) then
|
||||
write(ou,100) "Plot Type:", "Slice"
|
||||
|
|
@ -1387,7 +1422,7 @@ contains
|
|||
trim(to_str(pl % pixels(2)))
|
||||
else if (pl % type == PLOT_TYPE_VOXEL) then
|
||||
write(ou,100) "Voxels:", trim(to_str(pl % pixels(1))) // " " // &
|
||||
trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3)))
|
||||
trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3)))
|
||||
end if
|
||||
|
||||
write(ou,*)
|
||||
|
|
@ -1539,7 +1574,7 @@ contains
|
|||
|
||||
write(ou,100) 'Cell ID','No. Overlap Checks'
|
||||
|
||||
do i = 1, n_cells
|
||||
do i = 1, n_cells
|
||||
write(ou,101) cells(i) % id, overlap_check_cnt(i)
|
||||
if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1
|
||||
end do
|
||||
|
|
@ -1573,7 +1608,7 @@ contains
|
|||
integer :: k ! loop index for scoring bins
|
||||
integer :: n ! loop index for nuclides
|
||||
integer :: l ! loop index for user scores
|
||||
integer :: type ! type of tally filter
|
||||
integer :: type ! type of tally filter
|
||||
integer :: indent ! number of spaces to preceed output
|
||||
integer :: filter_index ! index in results array for filters
|
||||
integer :: score_index ! scoring bin index
|
||||
|
|
@ -1749,13 +1784,13 @@ contains
|
|||
score_name = 'P' // trim(to_str(t % scatt_order(k))) // &
|
||||
' Scattering Moment'
|
||||
end if
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
repeat(" ", indent), score_name, &
|
||||
to_str(t % results(score_index,filter_index) % sum), &
|
||||
trim(to_str(t % results(score_index,filter_index) % sum_sq))
|
||||
else if (t % score_bins(k) == SCORE_SCATTER_PN) then
|
||||
score_name = "Scattering Rate"
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
repeat(" ", indent), score_name, &
|
||||
to_str(t % results(score_index,filter_index) % sum), &
|
||||
trim(to_str(t % results(score_index,filter_index) % sum_sq))
|
||||
|
|
@ -1763,7 +1798,7 @@ contains
|
|||
score_index = score_index + 1
|
||||
score_name = 'P' // trim(to_str(n_order)) // &
|
||||
' Scattering Moment'
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
repeat(" ", indent), score_name, &
|
||||
to_str(t % results(score_index,filter_index) % sum), &
|
||||
trim(to_str(t % results(score_index,filter_index) % sum_sq))
|
||||
|
|
@ -1775,7 +1810,7 @@ contains
|
|||
else
|
||||
score_name = score_names(abs(t % score_bins(k)))
|
||||
end if
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
|
||||
repeat(" ", indent), score_name, &
|
||||
to_str(t % results(score_index,filter_index) % sum), &
|
||||
trim(to_str(t % results(score_index,filter_index) % sum_sq))
|
||||
|
|
@ -1813,8 +1848,8 @@ contains
|
|||
integer :: i_filter_ein ! index for incoming energy filter
|
||||
integer :: i_filter_surf ! index for surface filter
|
||||
integer :: n ! number of incoming energy bins
|
||||
integer :: len1 ! length of string
|
||||
integer :: len2 ! length of string
|
||||
integer :: len1 ! length of string
|
||||
integer :: len2 ! length of string
|
||||
integer :: filter_index ! index in results array for filters
|
||||
logical :: print_ebin ! should incoming energy bin be displayed?
|
||||
character(MAX_LINE_LEN) :: string
|
||||
|
|
@ -1864,14 +1899,14 @@ contains
|
|||
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
|
||||
matching_bins(i_filter_surf) = IN_RIGHT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Outgoing Current to Left", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
||||
matching_bins(i_filter_surf) = OUT_RIGHT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Incoming Current from Left", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
|
@ -1881,14 +1916,14 @@ contains
|
|||
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
|
||||
matching_bins(i_filter_surf) = IN_RIGHT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Incoming Current from Right", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
||||
matching_bins(i_filter_surf) = OUT_RIGHT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Outgoing Current to Right", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
|
@ -1898,14 +1933,14 @@ contains
|
|||
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
|
||||
matching_bins(i_filter_surf) = IN_FRONT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Outgoing Current to Back", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
||||
matching_bins(i_filter_surf) = OUT_FRONT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Incoming Current from Back", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
|
@ -1915,14 +1950,14 @@ contains
|
|||
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
|
||||
matching_bins(i_filter_surf) = IN_FRONT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Incoming Current from Front", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
||||
matching_bins(i_filter_surf) = OUT_FRONT
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Outgoing Current to Front", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
|
@ -1932,14 +1967,14 @@ contains
|
|||
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
|
||||
matching_bins(i_filter_surf) = IN_TOP
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Outgoing Current to Bottom", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
||||
matching_bins(i_filter_surf) = OUT_TOP
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Incoming Current from Bottom", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
|
@ -1949,14 +1984,14 @@ contains
|
|||
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
|
||||
matching_bins(i_filter_surf) = IN_TOP
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Incoming Current from Top", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
||||
matching_bins(i_filter_surf) = OUT_TOP
|
||||
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
|
||||
"Outgoing Current to Top", &
|
||||
to_str(t % results(1,filter_index) % sum), &
|
||||
trim(to_str(t % results(1,filter_index) % sum_sq))
|
||||
|
|
|
|||
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" } }? &
|
||||
|
|
|
|||
|
|
@ -3,6 +3,7 @@ module source
|
|||
use bank_header, only: Bank
|
||||
use constants
|
||||
use error, only: fatal_error
|
||||
use geometry, only: find_cell
|
||||
use geometry_header, only: BASE_UNIVERSE
|
||||
use global
|
||||
use math, only: maxwell_spectrum, watt_spectrum
|
||||
|
|
@ -73,6 +74,9 @@ contains
|
|||
real(8) :: p_max(3) ! maximum coordinates of source
|
||||
real(8) :: a ! Arbitrary parameter 'a'
|
||||
real(8) :: b ! Arbitrary parameter 'b'
|
||||
logical :: found ! Does the source particle exist within geometry?
|
||||
type(Particle) :: p ! Temporary particle for using find_cell
|
||||
integer, save :: num_resamples = 0 ! Number of resamples encountered
|
||||
|
||||
! Set weight to one by default
|
||||
site % wgt = ONE
|
||||
|
|
@ -80,11 +84,32 @@ contains
|
|||
! Sample position
|
||||
select case (external_source % type_space)
|
||||
case (SRC_SPACE_BOX)
|
||||
! Coordinates sampled uniformly over a box
|
||||
p_min = external_source % params_space(1:3)
|
||||
p_max = external_source % params_space(4:6)
|
||||
r = (/ (prn(), i = 1,3) /)
|
||||
site % xyz = p_min + r*(p_max - p_min)
|
||||
! Set particle defaults
|
||||
call p % initialize()
|
||||
! Repeat sampling source location until a good site has been found
|
||||
found = .false.
|
||||
do while (.not.found)
|
||||
! Coordinates sampled uniformly over a box
|
||||
p_min = external_source % params_space(1:3)
|
||||
p_max = external_source % params_space(4:6)
|
||||
r = (/ (prn(), i = 1,3) /)
|
||||
site % xyz = p_min + r*(p_max - p_min)
|
||||
|
||||
! Fill p with needed data
|
||||
p % coord0 % xyz = site % xyz
|
||||
p % coord0 % uvw = [ ONE, ZERO, ZERO ]
|
||||
|
||||
! Now search to see if location exists in geometry
|
||||
call find_cell(p, found)
|
||||
if (.not. found) then
|
||||
num_resamples = num_resamples + 1
|
||||
if (num_resamples == MAX_EXTSRC_RESAMPLES) then
|
||||
message = "Maximum number of external source spatial resamples &
|
||||
&reached!"
|
||||
call fatal_error()
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
||||
case (SRC_SPACE_POINT)
|
||||
! Point source
|
||||
|
|
@ -146,7 +171,7 @@ contains
|
|||
end subroutine sample_external_source
|
||||
|
||||
!===============================================================================
|
||||
! GET_SOURCE_PARTICLE returns the next source particle
|
||||
! GET_SOURCE_PARTICLE returns the next source particle
|
||||
!===============================================================================
|
||||
|
||||
subroutine get_source_particle(p, index_source)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -71,7 +71,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
|
||||
|
|
|
|||
|
|
@ -16,123 +16,161 @@ from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as Naviga
|
|||
import numpy as np
|
||||
|
||||
class AppForm(QMainWindow):
|
||||
def __init__(self, parent=None):
|
||||
def __init__(self, argv, parent=None):
|
||||
QMainWindow.__init__(self, parent)
|
||||
|
||||
# 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.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)
|
||||
|
||||
# Create the navigation toolbar, tied to the canvas
|
||||
self.mpl_toolbar = NavigationToolbar(self.canvas, self.main_frame)
|
||||
self.all_good = False
|
||||
while not self.all_good:
|
||||
if len(argv) > 1:
|
||||
cl_file = str(argv[1])
|
||||
else:
|
||||
cl_file = None
|
||||
self.get_file_data(cl_file)
|
||||
# Check that there are any mesh tallies at all
|
||||
if len(self.tally_ids) != 0:
|
||||
self.all_good = True
|
||||
else:
|
||||
# if there are not, the user will be given the choice to choose
|
||||
# another file (but only if using interactive chooser)
|
||||
if cl_file is None:
|
||||
choice = QMessageBox.critical(None, "Invalid StatePoint File",
|
||||
"File Does Not Contain Mesh " +
|
||||
"Tallies!" +
|
||||
"\nSelect Another File Or Quit",
|
||||
QMessageBox.Retry,
|
||||
QMessageBox.Abort)
|
||||
if choice == QMessageBox.Abort:
|
||||
self.all_good = False
|
||||
break
|
||||
else:
|
||||
print("Invalid StatePoint File; File Does Not Contain " +
|
||||
"Mesh Tallies!")
|
||||
self.all_good = False
|
||||
break
|
||||
|
||||
# Grid layout at bottom
|
||||
self.grid = QGridLayout()
|
||||
|
||||
# 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)
|
||||
if self.all_good:
|
||||
# Set maximum colorbar value by maximum tally data value
|
||||
self.maxvalue = self.datafile.tallies[0].results.max()
|
||||
|
||||
# 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)
|
||||
self.main_frame = QWidget()
|
||||
self.setCentralWidget(self.main_frame)
|
||||
|
||||
# Planar basis
|
||||
label_basis = QLabel("Basis:")
|
||||
self.basis = QComboBox()
|
||||
self.basis.addItems(['xy', 'yz', 'xz'])
|
||||
# Create the Figure, Canvas, and Axes
|
||||
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)
|
||||
|
||||
# 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)
|
||||
# Create the navigation toolbar, tied to the canvas
|
||||
self.mpl_toolbar = NavigationToolbar(self.canvas, self.main_frame)
|
||||
|
||||
# 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)
|
||||
|
||||
# Grid layout at bottom
|
||||
self.grid = QGridLayout()
|
||||
|
||||
self.label_filters = QLabel("Filter options:")
|
||||
# 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)
|
||||
|
||||
# Labels for all possible filters
|
||||
self.labels = {'cell': 'Cell: ', 'cellborn': 'Cell born: ',
|
||||
'surface': 'Surface: ', 'material': 'Material',
|
||||
'universe': 'Universe: ', 'energyin': 'Energy in: ',
|
||||
'energyout': 'Energy out: '}
|
||||
|
||||
# Empty reusable labels
|
||||
self.qlabels = {}
|
||||
for j in range(8):
|
||||
self.nextLabel = QLabel
|
||||
self.qlabels[j] = self.nextLabel
|
||||
|
||||
# Reusable comboboxes labelled with filter names
|
||||
self.boxes = {}
|
||||
for key in self.labels.keys():
|
||||
self.nextBox = QComboBox()
|
||||
self.connect(self.nextBox, SIGNAL('activated(int)'),
|
||||
# Tally selections
|
||||
label_tally = QLabel("Tally:")
|
||||
self.tally = QComboBox()
|
||||
# Only show options for the tallies with meshes
|
||||
self.tally.addItems([str(i + 1) for i in self.tally_ids])
|
||||
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)
|
||||
self.boxes[key] = self.nextBox
|
||||
|
||||
# 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)
|
||||
# Planar basis
|
||||
label_basis = QLabel("Basis:")
|
||||
self.basis = QComboBox()
|
||||
self.basis.addItems(['xy', 'yz', 'xz'])
|
||||
|
||||
# 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)
|
||||
# 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)
|
||||
|
||||
self._update()
|
||||
self.populate_boxes()
|
||||
self.on_draw()
|
||||
# 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)
|
||||
|
||||
def get_file_data(self):
|
||||
# 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)
|
||||
|
||||
|
||||
self.label_filters = QLabel("Filter options:")
|
||||
|
||||
# Labels for all possible filters
|
||||
self.labels = {'cell': 'Cell: ', 'cellborn': 'Cell born: ',
|
||||
'surface': 'Surface: ', 'material': 'Material',
|
||||
'universe': 'Universe: ', 'energyin': 'Energy in: ',
|
||||
'energyout': 'Energy out: '}
|
||||
|
||||
# Empty reusable labels
|
||||
self.qlabels = {}
|
||||
for j in range(8):
|
||||
self.nextLabel = QLabel
|
||||
self.qlabels[j] = self.nextLabel
|
||||
|
||||
# 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
|
||||
|
||||
# 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)
|
||||
|
||||
# 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)
|
||||
|
||||
self._update()
|
||||
self.populate_boxes()
|
||||
self.on_draw()
|
||||
|
||||
def get_file_data(self, cl_file=None):
|
||||
# Get data file name from "open file" browser
|
||||
filename = QFileDialog.getOpenFileName(self, 'Select statepoint file', '.')
|
||||
if cl_file is None:
|
||||
filename = QFileDialog.getOpenFileName(self,
|
||||
'Select statepoint file', '.')
|
||||
else:
|
||||
filename = cl_file
|
||||
|
||||
# Create StatePoint object and read in data
|
||||
self.datafile = StatePoint(str(filename))
|
||||
|
|
@ -141,32 +179,31 @@ class AppForm(QMainWindow):
|
|||
|
||||
self.setWindowTitle('Core Map Tool : ' + str(self.datafile.path))
|
||||
|
||||
# Set maximum colorbar value by maximum tally data value
|
||||
self.maxvalue = self.datafile.tallies[0].results.max()
|
||||
|
||||
self.labelList = []
|
||||
|
||||
# Read mesh dimensions
|
||||
# for mesh in self.datafile.meshes:
|
||||
# self.nx, self.ny, self.nz = mesh.dimension
|
||||
# for mesh in self.datafile.meshes:
|
||||
# self.nx, self.ny, self.nz = mesh.dimension
|
||||
|
||||
# Read filter types from statepoint file
|
||||
# Find which tallies have meshes so the rest can be ignored,
|
||||
# and for these tallies read the filter and score types
|
||||
self.tally_ids = []
|
||||
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)
|
||||
|
||||
# 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
|
||||
for itally, tally in enumerate(self.datafile.tallies):
|
||||
if 'mesh' in tally.filters:
|
||||
# Then we have a good tally, store the ID, filters and
|
||||
# scores
|
||||
self.tally_ids.append(itally)
|
||||
self.filter_types = []
|
||||
for f in tally.filters:
|
||||
self.filter_types.append(f)
|
||||
self.tally_list.append(self.filter_types)
|
||||
self.score_types = []
|
||||
for s in tally.scores:
|
||||
self.score_types.append(s)
|
||||
self.tally_scores.append(self.score_types)
|
||||
|
||||
def on_draw(self):
|
||||
""" Redraws the figure
|
||||
|
|
@ -178,70 +215,73 @@ class AppForm(QMainWindow):
|
|||
axial_level = self.axial_level.currentIndex() + 1
|
||||
is_mean = self.mean.currentIndex()
|
||||
|
||||
# get current tally index
|
||||
tally_id = self.tally_ids[self.tally.currentIndex()]
|
||||
|
||||
# Create spec_list
|
||||
spec_list = []
|
||||
for tally in self.datafile.tallies[self.tally.currentIndex()].filters.values():
|
||||
for tally in self.datafile.tallies[tally_id].filters.values():
|
||||
if tally.type == 'mesh':
|
||||
continue
|
||||
index = self.boxes[tally.type].currentIndex()
|
||||
spec_list.append((tally.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':
|
||||
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))],
|
||||
matrix[i,j] = self.datafile.get_value(tally_id,
|
||||
spec_list + [('mesh', (i + 1, j + 1, axial_level))],
|
||||
self.scoreBox.currentIndex())[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))],
|
||||
mean_val = self.datafile.get_value(tally_id,
|
||||
spec_list + [('mesh', (i + 1, j + 1, axial_level))],
|
||||
self.scoreBox.currentIndex())[0]
|
||||
if mean_val > 0.0:
|
||||
matrix[i,j] = matrix[i,j] / mean_val
|
||||
else:
|
||||
matrix[i,j] = 0.0
|
||||
|
||||
|
||||
elif self.basis.currentText() == '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))],
|
||||
matrix[i,j] = self.datafile.get_value(tally_id,
|
||||
spec_list + [('mesh', (axial_level, i + 1, j + 1))],
|
||||
self.scoreBox.currentIndex())[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))],
|
||||
mean_val = self.datafile.get_value(tally_id,
|
||||
spec_list + [('mesh', (axial_level, i + 1, j + 1))],
|
||||
self.scoreBox.currentIndex())[0]
|
||||
if mean_val > 0.0:
|
||||
matrix[i,j] = matrix[i,j] / mean_val
|
||||
else:
|
||||
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))],
|
||||
matrix[i,j] = self.datafile.get_value(tally_id,
|
||||
spec_list + [('mesh', (i + 1, axial_level, j + 1))],
|
||||
self.scoreBox.currentIndex())[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))],
|
||||
mean_val = self.datafile.get_value(tally_id,
|
||||
spec_list + [('mesh', (i + 1, axial_level, j + 1))],
|
||||
self.scoreBox.currentIndex())[0]
|
||||
if mean_val > 0.0:
|
||||
matrix[i,j] = matrix[i,j] / mean_val
|
||||
|
|
@ -255,8 +295,8 @@ class AppForm(QMainWindow):
|
|||
|
||||
# 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="nearest", origin='lower')
|
||||
self.fig.colorbar(cax)
|
||||
|
||||
self.axes.set_xticks([])
|
||||
|
|
@ -271,9 +311,11 @@ class AppForm(QMainWindow):
|
|||
'''
|
||||
# print 'Calling _update...'
|
||||
|
||||
# get current tally index
|
||||
tally_id = self.tally_ids[self.tally.currentIndex()]
|
||||
|
||||
self.mesh = self.datafile.meshes[
|
||||
self.datafile.tallies[
|
||||
self.tally.currentIndex()].filters['mesh'].bins[0] - 1]
|
||||
self.datafile.tallies[tally_id].filters['mesh'].bins[0] - 1]
|
||||
|
||||
self.nx, self.ny, self.nz = self.mesh.dimension
|
||||
|
||||
|
|
@ -289,8 +331,7 @@ class AppForm(QMainWindow):
|
|||
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()
|
||||
self.maxvalue = self.datafile.tallies[tally_id].results.max()
|
||||
# print self.maxvalue
|
||||
|
||||
# Clear and hide old filter labels
|
||||
|
|
@ -307,6 +348,9 @@ class AppForm(QMainWindow):
|
|||
def populate_boxes(self):
|
||||
# print 'Calling populate_boxes...'
|
||||
|
||||
# get current tally index
|
||||
tally_id = self.tally_ids[self.tally.currentIndex()]
|
||||
|
||||
n = 5
|
||||
labels = {'cell': 'Cell : ',
|
||||
'cellborn': 'Cell born: ',
|
||||
|
|
@ -317,8 +361,7 @@ class AppForm(QMainWindow):
|
|||
# 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]
|
||||
nextFilter = self.datafile.tallies[tally_id].filters[element]
|
||||
if element == 'mesh':
|
||||
continue
|
||||
|
||||
|
|
@ -337,7 +380,7 @@ class AppForm(QMainWindow):
|
|||
|
||||
elif element == 'energyin' or element == 'energyout':
|
||||
for i in range(nextFilter.length):
|
||||
text = (str(nextFilter.bins[i]) + ' to ' +
|
||||
text = (str(nextFilter.bins[i]) + ' to ' +
|
||||
str(nextFilter.bins[i+1]))
|
||||
combobox.addItem(text)
|
||||
|
||||
|
|
@ -351,9 +394,10 @@ class AppForm(QMainWindow):
|
|||
|
||||
def main():
|
||||
app = QApplication(sys.argv)
|
||||
form = AppForm()
|
||||
form.show()
|
||||
app.exec_()
|
||||
form = AppForm(app.arguments())
|
||||
if form.all_good:
|
||||
form.show()
|
||||
app.exec_()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
|
|
|||
|
|
@ -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]
|
||||
|
|
|
|||
|
|
@ -44,9 +44,11 @@
|
|||
<surface id="33" type="z-cylinder" coeffs="0 0 5.995" />
|
||||
|
||||
|
||||
<!-- Beryllium Reflector -->
|
||||
<!-- Bottom Beryllium Reflector -->
|
||||
<cell id="1" material="2" surfaces=" 1 -2 -30" />
|
||||
<cell id="2" material="2" surfaces="28 -29 -30" />
|
||||
|
||||
<!-- Top ZrH Reflector -->
|
||||
<cell id="2" material="7" surfaces="28 -29 -30" />
|
||||
|
||||
<!-- Duralumin base and top -->
|
||||
<cell id="3" material="4" surfaces=" 2 -4 -30 31" />
|
||||
|
|
|
|||
|
|
@ -1,7 +1,8 @@
|
|||
<?xml version="1.0"?>
|
||||
<materials>
|
||||
|
||||
<!-- pu-met-fast-021 case 2 -->
|
||||
<!-- pu-met-fast-021 case 2, modified to include Zr-H to exercise
|
||||
the capability to assign a single S(a,b) table to multiple nuclides -->
|
||||
|
||||
<default_xs>70c</default_xs>
|
||||
|
||||
|
|
@ -111,4 +112,19 @@
|
|||
<nuclide name="Cu-65" ao="2.2809e-04" />
|
||||
</material>
|
||||
|
||||
<!-- Additional ZrH to replace one of the BeO reflectors -->
|
||||
<material id="7">
|
||||
<density value="5.90168" units="g/cm3" />
|
||||
<nuclide name="H-1" ao="0.3" />
|
||||
<nuclide name="Zr-90" ao="0.15" />
|
||||
<nuclide name="Zr-91" ao="0.1" />
|
||||
<nuclide name="Zr-92" ao="0.1" />
|
||||
<nuclide name="Zr-94" ao="0.05" />
|
||||
<nuclide name="Zr-96" ao="0.05" />
|
||||
<nuclide name="U-235" ao="0.1" />
|
||||
<nuclide name="U-238" ao="0.15" />
|
||||
<sab name="h/zr" xs="10t" />
|
||||
<sab name="zr/h" xs="10t" />
|
||||
</material>
|
||||
|
||||
</materials>
|
||||
|
|
|
|||
|
|
@ -1,2 +1,2 @@
|
|||
k-combined:
|
||||
1.013747E+00 3.067630E-02
|
||||
9.761880E-01 9.170415E-03
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue