crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

merged Paul's latest changes and made additional changes to handle batches

Browse source
This commit is contained in:
Bryan Herman 2012-03-27 11:47:19 -07:00
commit 599ddf9f8e
56 changed files with 1762 additions and 668 deletions
Download patch file Download diff file Expand all files Collapse all files

18
LICENSE
View file

@ -0,0 +1,18 @@
Copyright (c) 2011-2012 Massachusetts Institute of Technology
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

1
docs/source/index.rst
View file

@ -25,6 +25,7 @@ Contents
:maxdepth: 1
install
releasenotes/index
methods/index
usersguide/index
developers

15
docs/source/releasenotes/index.rst Normal file
View file

@ -0,0 +1,15 @@
.. _releasenotes:
=============
Release Notes
=============
The release notes for OpenMC give a list of system requirements, new features,
bugs fixed, and known issues for each successive release.
.. toctree::
:maxdepth: 1
notes_0.4.2
notes_0.4.1
notes_0.4.0

36
docs/source/releasenotes/notes_0.4.0.rst Normal file
View file

@ -0,0 +1,36 @@
.. _notes_0.4.0:
==============================
Release Notes for OpenMC 0.4.0
==============================
-------------------
System Requirements
-------------------
There are no special requirements for running the OpenMC code. As of this
release, OpenMC has been tested on a variety of Linux distributions as well as
Mac OS X. However, it has not been tested yet on any releases of Microsoft
Windows. Memory requirements will vary depending on the size of the problem at
hand (mostly on the number of nuclides in the problem).
------------
New Features
------------
- The probability table method for treatment of energy self-shielding in the
unresolved resonance range has been implemented and is now turned on by
default.
- Calculation of Shannon entropy for assessing convergence of the fission source
distribution.
- Ability to compile with the PGI Fortran compiler.
- Ability to run on IBM BlueGene/P machines.
- Completely rewrote how nested universes are handled. Geometry is now much more
robust.
---------
Bug Fixes
---------
- Many geometry errors have been fixed. The Monte Carlo performance benchmark
can now be successfully run in OpenMC.

55
docs/source/releasenotes/notes_0.4.1.rst Normal file
View file

@ -0,0 +1,55 @@
.. _notes_0.4.1:
==============================
Release Notes for OpenMC 0.4.1
==============================
-------------------
System Requirements
-------------------
There are no special requirements for running the OpenMC code. As of this
release, OpenMC has been tested on a variety of Linux distributions as well as
Mac OS X. However, it has not been tested yet on any releases of Microsoft
Windows. Memory requirements will vary depending on the size of the problem at
hand (mostly on the number of nuclides in the problem).
------------
New Features
------------
- A batching method has been implemented so that statistics can be calculated
based on multiple generations instead of a single generation. This can help to
overcome problems with underpredicted variance in problems where there is
correlation between successive fission source iterations.
- Users now have the option to select a non-unionized energy grid for problems
with many nuclides where the use of a unionized grid is not feasible.
- Improved plotting capability (Nick Horelik). The plotting input is now in
``plots.xml`` instead of ``plot.xml``.
- Added multiple estimators for k-effective and added a global tally for
leakage.
- Moved cross section-related output into cross_sections.out.
- Improved timing capabilities.
- Can now use more than 2**31 - 1 particles per generation.
- Improved fission bank synchronization method. This also necessitated changing
the source bank to be of type Bank rather than of type Particle.
- Added HDF5 output (not complete yet).
- Major changes to tally implementation.
---------
Bug Fixes
---------
- `b206a8`_: Fixed subtle error in the sampling of energy distributions.
- `800742`_: Fixed error in sampling of angle and rotating angles.
- `a07c08`_: Fixed bug in linear-linear interpolation during sampling energy.
- `a75283`_: Fixed energy and energyout tally filters to support many bins.
- `95cfac`_: Fixed error in cell neighbor searches.
- `83a803`_: Fixed bug related to probability tables.
.. _b206a8: https://github.com/mit-crpg/openmc/commit/b206a8
.. _800742: https://github.com/mit-crpg/openmc/commit/800742
.. _a07c08: https://github.com/mit-crpg/openmc/commit/a07c08
.. _a75283: https://github.com/mit-crpg/openmc/commit/a75283
.. _95cfac: https://github.com/mit-crpg/openmc/commit/95cfac
.. _83a803: https://github.com/mit-crpg/openmc/commit/83a803

37
docs/source/releasenotes/notes_0.4.2.rst Normal file
View file

@ -0,0 +1,37 @@
.. _notes_0.4.2:
==============================
Release Notes for OpenMC 0.4.2
==============================
.. note::
These release notes are for an upcoming release of OpenMC and are still
subject to change.
-------------------
System Requirements
-------------------
There are no special requirements for running the OpenMC code. As of this
release, OpenMC has been tested on a variety of Linux distributions as well as
Mac OS X. However, it has not been tested yet on any versions of Microsoft
Windows. Memory requirements will vary depending on the size of the problem at
hand (mostly on the number of nuclides in the problem).
------------
New Features
------------
- Ability to specify 'point' source type.
- Ability to change random number seed.
- User's can now specify units='sum' on a <density> tag. This tells the code
that the total material density is the sum of the atom fractions listed for
each nuclide on the material.
---------
Bug Fixes
---------
- `2a5b9c`_: Fixed regression in plotting.
.. _2a5b9c: https://github.com/mit-crpg/openmc/commit/2a5b9c

162
docs/source/usersguide/input.rst
View file

@ -307,16 +307,23 @@ settings.xml file.
The ``criticality`` element indicates that a criticality calculation should be
performed. It has the following attributes/sub-elements:
:cycles:
The number of total fission source iterations.
:batches:
The total number of batches, where each batch corresponds to multiple
fission source iterations. Batching is done to eliminate correlation between
realizations of random variables.
*Default*: None
:generations_per_batch:
The number of total fission source iterations per batch.
*Default*: 1
:inactive:
The number of inactive fission source iterations. In general, the starting
cycles in a criticality calculation can not be used to contribute to tallies
since the fission source distribution and eigenvalue are generally not
converged immediately
The number of inactive batches. In general, the starting cycles in a
criticality calculation can not be used to contribute to tallies since the
fission source distribution and eigenvalue are generally not converged
immediately.
*Default*: None
@ -353,6 +360,18 @@ default. This element has the following attributes/sub-elements:
*Default*: 1.0
``energy_grid`` Element
-----------------------
The ``energy_grid`` element determines the treatment of the energy grid during a
simulation. Setting this element to "nuclide" will cause OpenMC to use a
nuclide's energy grid when determining what points to interpolate between for
determining cross sections (i.e. non-unionized energy grid). To use a unionized
energy grid, set this element to "union". Note that the unionized energy grid
treatment is slightly different than that employed in Serpent.
*Default*: union
``entropy`` Element
-------------------
@ -385,6 +404,14 @@ sub-elements and can be set to either "off" or "on".
*Default*: on
``seed`` Element
----------------
The ``seed`` element is used to set the seed used for the linear congruential
pseudo-random number generator.
*Default*: 1
``source`` Element
------------------
@ -573,40 +600,115 @@ tallies. This element should be followed by "yes" or "no"
Geometry Plotting Specification -- plot.xml
-------------------------------------------
A rudimentary plotting capability is available in OpenMC by specifying a
plot.xml file and subsequently running with the command-line flag ``-plot``. The
root element of the plot.xml is simply ``<plot>`` and four sub-elements can be
defined to configure the plotting range and resolution.
A basic 2D plotting capability is available in OpenMC by creating a
plots.xml file and subsequently running with the command-line flag ``-plot``. The
root element of the plot.xml is simply ``<plots>`` and any number output
figures can be defined with ``<plot>`` sub-elements.
``origin`` Element
``plot`` Element
------------------
The ``origin`` element has no attributes/sub-elements and indicates the
Cartesian coordinates of the center of the plot.
Each plot must contain a combination of the following attributes or sub-elements:
*Default*: None
:id:
The unique ``id`` of the plot.
``width`` Element
-----------------
*Default*: None - Required entry
The ``width`` element has no attributes/sub-elements and indicates the width of
the plot in each of the basis directions.
:filename:
Filename for the output plot file.
*Default*: None
*Default*: "plot"
``basis`` Element
-----------------
:color:
Keyword for plot coloring. This can only be either ``cell`` or ``mat``,
which colors regions by cells and materials, respectively.
The ``basis`` element has no attributes/sub-elements and indicates the specified
basis for plotting. The only option option currently accepted is "xy".
*Default*: ``cell``
*Default*: xy
:origin:
Specifies the XYZ coordinate of the center of the plot. Should be 3 floats
separated by spaces.
``pixel`` Element
-----------------
*Default*: None - Required entry
The ``pixel`` element has no attributes/sub-elements and indicates the distance
between horizontal rays sent through the geometry to record surface crossings. A
smaller ``pixel`` will result in a higher-resolution plot.
:width:
Specifies the width of the plot along each of the basis directions.
Should be 2 or 3 floats separated by spaces for 2D plots and 3D plots,
respectively.
*Default*: 0.01
*Default*: None - Required entry
:type:
Keyword for type of plot to be produced. Currently only ``slice`` plots are
implemented, which create 2D pixel maps saved in the PPM file format. PPM
files can be displayed in most viewers (e.g. the default Gnome viewer,
IrfanView, etc.).
.. note:: Since the PPM format is saved without any kind of compression,
the resulting file sizes can be quite large. Saving the image in
the PNG format can often times reduce the file size by orders of
magnitude without any loss of image quality.
*Default*: "slice"
``plot`` elements of ``type`` ``slice`` also contain the following attributes or
sub-elements:
:basis:
Keyword specifying the plane of the plot for ``slice`` type plots. Can be
one of: ``xy``, ``xz``, ``yz``.
*Default*: ``xy``
:pixels:
Specifies the number of pixes to be used along each of the basis directions for
``slice`` plots. Should be 2 integers separated by spaces.
.. warning:: The ``pixels`` input determines the output file size. For the PPM
format, 10 million pixels will result in a file just under 30 MB in
size.
.. warning:: If the aspect ratio defined in ``pixels`` does not match the aspect
ratio defined in ``width`` the plot may appear stretched or squeezed.
.. warning:: Geometry features along a basis direction smaller than ``width``/``pixels``
along that basis direction may not appear in the plot.
*Default*: None - Required entry for ``slice`` plots
:background:
Specifies the RGB color of the regions where no OpenMC cell can be found. Should
be 3 integers deparated by spaces.
*Default*: 0 0 0 (white)
:col_spec:
Any number of this optional tag may be included in each ``plot`` element, which can
override the default random colors for cells or materials. Each ``col_spec``
element must contain ``id`` and ``rgb`` sub-elements.
:id:
Specifies the cell or material unique id for the color specification.
:rgb:
Specifies the custom color for the cell or material. Should be 3 intergers separated
by spaces.
*Default*: None
:mask:
The special ``mask`` sub-element allows for the selective plotting of *only*
user-specified cells or materials. Only one ``mask`` element is allowed per ``plot``
element, and it must contain as atributes or sub-elements a background masking color and
a list of cells or materials to plot:
:components:
List of unique ``id`` numbers of the cells or materials to plot. Should be any number
of integers separated by spaces.
:background:
Color to apply to all cells or materials not in the ``components`` list of cells or
materials to plot. This overrides any ``col_spec`` color specifications.
*Default*: None

11
examples/basic/materials.xml
View file

@ -1,17 +1,18 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<default_xs>70c</default_xs>
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="03c" ao="1.0" />
<nuclide name="U-235" ao="1.0" />
</material>
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H-1" xs="03c" ao="2.0" />
<nuclide name="O-16" xs="03c" ao="1.0" />
<sab name="lwe7" xs="00t" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<sab name="lwtr" xs="10t" />
</material>
</materials>

7
examples/basic/settings.xml
View file

@ -1,14 +1,9 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Path for cross section library -->
<cross_sections>
/home/paulromano/openmc/cross_sections_serpent.xml
</cross_sections>
<!-- Parameters for criticality calculation -->
<criticality>
<cycles>15</cycles>
<batches>15</batches>
<inactive>5</inactive>
<particles>10000</particles>
</criticality>

28
examples/basic/tallies.xml Normal file
View file

@ -0,0 +1,28 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filters>
<cell>100</cell>
</filters>
<scores>total scatter nu-scatter absorption fission nu-fission</scores>
</tally>
<tally id="2">
<filters>
<cell>100</cell>
<energy>0 20.0</energy>
</filters>
<scores>total scatter nu-scatter absorption fission nu-fission</scores>
</tally>
<tally id="3">
<filters>
<cell>100</cell>
<energy>0 20.0</energy>
<energyout>0 20.0</energyout>
</filters>
<scores>scatter nu-scatter nu-fission</scores>
</tally>
</tallies>

4
examples/lattice/nested/geometry.xml
View file

@ -14,7 +14,7 @@
<lattice id="4">
<type>rectangular</type>
<dimension>2 2</dimension>
<origin>-1.0 -1.0</origin>
<lower_left>-1.0 -1.0</lower_left>
<width>1.0 1.0</width>
<universes>
1 2
@ -26,7 +26,7 @@
<lattice id="6">
<type>rectangular</type>
<dimension>2 2</dimension>
<origin>-2.0 -2.0</origin>
<lower_left>-2.0 -2.0</lower_left>
<width>2.0 2.0</width>
<universes>
5 5

6
examples/lattice/nested/plot.xml
View file

@ -1,6 +0,0 @@
<?xml version="1.0"?>
<plot>
<origin>0. 0. 0.</origin>
<width>3.999 3.999</width>
<pixel>0.01</pixel>
</plot>

10
examples/lattice/nested/plots.xml Normal file
View file

@ -0,0 +1,10 @@
<?xml version="1.0"?>
<plots>
<plot id="1" color="mat">
<origin>0. 0. 0.</origin>
<width>4.0 4.0</width>
<pixels>400 400</pixels>
</plot>
</plots>

2
examples/lattice/nested/settings.xml
View file

@ -3,7 +3,7 @@
<!-- Parameters for criticality calculation -->
<criticality>
<cycles>20</cycles>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
</criticality>

6
examples/lattice/nested/tallies.xml
View file

@ -4,7 +4,7 @@
<mesh id="1">
<type>rectangular</type>
<dimension>4 4</dimension>
<origin>-2.0 -2.0</origin>
<lower_left>-2.0 -2.0</lower_left>
<width>1.0 1.0</width>
</mesh>
@ -13,7 +13,7 @@
<mesh>1</mesh>
</filters>
<macros>total</macros>
<scores>total</scores>
</tally>
</tallies>
</tallies>

2
examples/lattice/simple/geometry.xml
View file

@ -12,7 +12,7 @@
<lattice id="5">
<type>rectangular</type>
<dimension>4 4</dimension>
<origin>-2.0 -2.0</origin>
<lower_left>-2.0 -2.0</lower_left>
<width>1.0 1.0</width>
<universes>
1 2 1 2

6
examples/lattice/simple/plot.xml
View file

@ -1,6 +0,0 @@
<?xml version="1.0"?>
<plot>
<origin>0. 0. 0.</origin>
<width>3.999 3.999</width>
<pixel>0.01</pixel>
</plot>

10
examples/lattice/simple/plots.xml Normal file
View file

@ -0,0 +1,10 @@
<?xml version="1.0"?>
<plots>
<plot id="1" color="mat">
<origin>0. 0. 0.</origin>
<width>4.0 4.0</width>
<pixels>400 400</pixels>
</plot>
</plots>

2
examples/lattice/simple/settings.xml
View file

@ -3,7 +3,7 @@
<!-- Parameters for criticality calculation -->
<criticality>
<cycles>20</cycles>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
</criticality>

6
examples/lattice/simple/tallies.xml
View file

@ -4,7 +4,7 @@
<mesh id="1">
<type>rectangular</type>
<dimension>4 4</dimension>
<origin>-2.0 -2.0</origin>
<lower_left>-2.0 -2.0</lower_left>
<width>1.0 1.0</width>
</mesh>
@ -13,7 +13,7 @@
<mesh>1</mesh>
</filters>
<macros>total</macros>
<scores>total</scores>
</tally>
</tallies>
</tallies>

5
examples/reflective/materials.xml
View file

@ -1,10 +1,11 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<default_xs>70c</default_xs>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="03c" ao="1.0" />
<nuclide name="U-235" ao="1.0" />
</material>
</materials>

7
examples/reflective/settings.xml
View file

@ -1,14 +1,9 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Path for cross section library -->
<cross_sections>
/home/paulromano/openmc/cross_sections_serpent.xml
</cross_sections>
<!-- Parameters for criticality calculation -->
<criticality>
<cycles>500</cycles>
<batches>500</batches>
<inactive>10</inactive>
<particles>10000</particles>
</criticality>

8
src/DEPENDENCIES
View file

@ -86,6 +86,7 @@ energy_grid.o: output.o
finalize.o: global.o
finalize.o: hdf5_interface.o
finalize.o: output.o
finalize.o: source.o
finalize.o: tally.o
finalize.o: timing.o
@ -145,6 +146,7 @@ initialize.o: random_lcg.o
initialize.o: source.o
initialize.o: string.o
initialize.o: tally.o
initialize.o: tally_header.o
initialize.o: timing.o
input_xml.o: cmfd_input.o
@ -194,12 +196,13 @@ main.o: intercycle.o
main.o: output.o
main.o: physics.o
main.o: plot.o
main.o: random_lcg.o
main.o: source.o
main.o: string.o
main.o: tally.o
main.o: timing.o
mesh.o: constants.o
mesh.o: global.o
mesh.o: mesh_header.o
mesh.o: particle_header.o
@ -226,6 +229,7 @@ physics.o: geometry.o
physics.o: geometry_header.o
physics.o: global.o
physics.o: interpolation.o
physics.o: material_header.o
physics.o: output.o
physics.o: particle_header.o
physics.o: random_lcg.o
@ -233,6 +237,8 @@ physics.o: search.o
physics.o: string.o
physics.o: tally.o
plot_header.o: constants.o
plot.o: constants.o
plot.o: error.o
plot.o: geometry.o

8
src/Makefile
View file

@ -23,6 +23,12 @@ OPTIMIZE = no
USE_MPI = no
USE_HDF5 = yes
#===============================================================================
# Add git SHA-1 hash
#===============================================================================
GIT_SHA1 = $(shell git log -1 | head -n 1 | awk '{print $$2}')
#===============================================================================
# GNU Fortran compiler options
#===============================================================================
@ -233,7 +239,7 @@ neat:
.PHONY: all xml-fortran clean neat distclean
%.o: %.F90
$(F90) $(F90FLAGS) -Ixml-fortran -Ixml-fortran/templates -c $<
$(F90) $(F90FLAGS) -DGIT_SHA1="\"$(GIT_SHA1)\"" -Ixml-fortran -Ixml-fortran/templates -c $<
#===============================================================================
# Dependencies

7
src/ace.F90
View file

@ -22,8 +22,6 @@ module ace
real(8), allocatable :: XSS(:) ! Cross section data
integer :: XSS_index ! current index in XSS data
type(DictionaryCI), pointer :: already_read => null()
private :: NXS
private :: JXS
private :: XSS
@ -47,6 +45,7 @@ contains
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
type(SAB_Table), pointer :: sab => null()
type(DictionaryCI), pointer :: already_read => null()
! ==========================================================================
! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES
@ -154,6 +153,8 @@ contains
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
! Read the ACE table into the appropriate entry on the nuclides
! array
call read_ace_table(index_nuclides, index_list)
! Print out information on table to cross_sections.out file
@ -172,6 +173,8 @@ contains
index_list = dict_get_key(xs_listing_dict, name)
index_sab = dict_get_key(sab_dict, name)
! Read the ACE table into the appropriate entry on the sab_tables
! array
call read_ace_table(index_sab, index_list)
call dict_add_key(already_read, name, 0)

5
src/bank_header.F90
View file

@ -9,6 +9,11 @@ module bank_header
!===============================================================================
type Bank
! The 'sequence' attribute is used here to ensure that the data listed
! appears in the given order. This is important for MPI purposes when bank
! sites are sent from one processor to another.
sequence
integer(8) :: id ! Unique ID
real(8) :: xyz(3) ! location of bank particle
real(8) :: uvw(3) ! diretional cosines

32
src/constants.F90
View file

@ -7,7 +7,7 @@ module constants
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 4
integer, parameter :: VERSION_RELEASE = 0
integer, parameter :: VERSION_RELEASE = 1
! ============================================================================
! ADJUSTABLE PARAMETERS
@ -243,7 +243,7 @@ module constants
EVENT_FISSION = 3
! Tally score type
integer, parameter :: N_SCORE_TYPES = 15
integer, parameter :: N_SCORE_TYPES = 16
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
@ -252,14 +252,15 @@ module constants
SCORE_SCATTER_1 = -5, & ! first scattering moment
SCORE_SCATTER_2 = -6, & ! second scattering moment
SCORE_SCATTER_3 = -7, & ! third scattering moment
SCORE_N_1N = -8, & ! (n,1n) rate
SCORE_N_2N = -9, & ! (n,2n) rate
SCORE_N_3N = -10, & ! (n,3n) rate
SCORE_N_4N = -11, & ! (n,4n) rate
SCORE_ABSORPTION = -12, & ! absorption rate
SCORE_FISSION = -13, & ! fission rate
SCORE_NU_FISSION = -14, & ! neutron production rate
SCORE_CURRENT = -15 ! partial current
SCORE_DIFFUSION = -8, & ! diffusion coefficient
SCORE_N_1N = -9, & ! (n,1n) rate
SCORE_N_2N = -10, & ! (n,2n) rate
SCORE_N_3N = -11, & ! (n,3n) rate
SCORE_N_4N = -12, & ! (n,4n) rate
SCORE_ABSORPTION = -13, & ! absorption rate
SCORE_FISSION = -14, & ! fission rate
SCORE_NU_FISSION = -15, & ! neutron production rate
SCORE_CURRENT = -16 ! partial current
! Tally map bin finding
integer, parameter :: NO_BIN_FOUND = -1
@ -312,11 +313,17 @@ module constants
integer, parameter :: ERROR_INT = -huge(0)
real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8
! Energy grid methods
integer, parameter :: &
GRID_NUCLIDE = 1, & ! non-unionized energy grid (MCNP)
GRID_UNION = 2, & ! union grid with pointers
GRID_LETHARGY = 3 ! lethargy mapping (MC21)
! Source types
integer, parameter :: &
SRC_BOX = 1, & ! Source in a rectangular prism
SRC_CELL = 2, & ! Source in a cell
SRC_SURFACE = 3 ! Source on a surface
SRC_POINT = 2, & ! Source at a single point
SRC_FILE = 3 ! Source from a file
integer, parameter :: &
PROB_SOURCE = 1, & ! External source problem
@ -327,5 +334,6 @@ module constants
integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file
integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file
integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file
integer, parameter :: UNIT_SOURCE = 15 ! unit # for writing source file
end module constants

29
src/cross_section.F90
View file

@ -37,7 +37,7 @@ contains
mat => materials(p % material)
! Find energy index on unionized grid
call find_energy_index()
if (grid_method == GRID_UNION) call find_energy_index()
! Check if there's an S(a,b) table for this material
if (mat % has_sab_table) then
@ -107,11 +107,28 @@ contains
! Set pointer to nuclide
nuc => nuclides(index_nuclide)
! TODO: If not using unionized energy grid, we need to find the index on the
! nuclide energy grid using lethargy mapping or whatever other technique
! Determine index on nuclide energy grid
select case (grid_method)
case (GRID_UNION)
! If we're using the unionized grid with pointers, finding the index on
! the nuclide energy grid is as simple as looking up the pointer
! get index on nuclide energy grid
IE = nuc % grid_index(p % IE)
IE = nuc % grid_index(p % IE)
case (GRID_NUCLIDE)
! If we're not using the unionized grid, we have to do a binary search on
! the nuclide energy grid in order to determine which points to
! interpolate between
if (p % E < nuc % energy(1)) then
IE = 1
elseif (p % E > nuc % energy(nuc % n_grid)) then
IE = nuc % n_grid - 1
else
IE = binary_search(nuc % energy, nuc % n_grid, p % E)
end if
end select
! check for rare case where two energy points are the same
if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1
@ -171,7 +188,7 @@ contains
end if
! Set last evaluated energy
micro_xs(index_nuclide) % last_E = p % E
if (index_sab == 0) micro_xs(index_nuclide) % last_E = p % E
end subroutine calculate_nuclide_xs

2
src/energy_grid.F90
View file

@ -37,7 +37,7 @@ contains
! Create dictionary for keeping track of cross sections already added
call dict_create(already_added)
! Loop over all files
! Loop over all materials
do i = 1, n_materials
mat => materials(i)

7
src/error.F90
View file

@ -76,9 +76,10 @@ contains
write(ERROR_UNIT,*)
! Write information on current cycle and particle
if (current_cycle > 0) then
write(ERROR_UNIT,'(1X,A,I11) ') 'Cycle: ', current_cycle
write(ERROR_UNIT,'(1X,A,I11)') 'Particle: ', p % id
if (current_batch > 0) then
write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch
write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen
write(ERROR_UNIT,'(1X,A,I12)') 'Particle: ', p % id
write(ERROR_UNIT,*)
end if

6
src/finalize.F90
View file

@ -2,6 +2,7 @@ module finalize
use global
use output, only: print_runtime
use source, only: write_source_binary
use tally, only: write_tallies, tally_statistics
use timing, only: timer_start, timer_stop
@ -23,10 +24,13 @@ contains
! Start finalization timer
call timer_start(time_finalize)
! Calculate statistics for tallies and write to tallies.out
if (.not. plotting) then
! Calculate statistics for tallies and write to tallies.out
call tally_statistics()
if (master) call write_tallies()
! Write out binary source
if (write_source) call write_source_binary()
end if
! stop timers and show timing statistics

6
src/geometry.F90
View file

@ -136,6 +136,12 @@ contains
! Set cell on this level
p % coord % cell = index_cell
! Show cell information on trace
if (verbosity >= 10 .or. trace) then
message = " Entering cell " // trim(to_str(c % id))
call write_message()
end if
if (c % type == CELL_NORMAL) then
! =================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH

1
src/geometry_header.F90
View file

@ -64,7 +64,6 @@ module geometry_header
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
integer, allocatable :: rgb(:) ! plotting color
end type Cell
! array index of universe 0

30
src/global.F90
View file

@ -83,7 +83,8 @@ module global
type(DictionaryCI), pointer :: xs_listing_dict => null()
! Unionized energy grid
integer :: n_grid ! number of points on unionized grid
integer :: grid_method ! how to treat the energy grid
integer :: n_grid ! number of points on unionized grid
real(8), allocatable :: e_grid(:) ! energies on unionized grid
! Unreoslved resonance probablity tables
@ -126,11 +127,13 @@ module global
! ============================================================================
! CRITICALITY SIMULATION VARIABLES
integer(8) :: n_particles = 10000 ! # of particles per cycle
integer :: n_cycles = 500 ! # of cycles
integer :: n_inactive = 50 ! # of inactive cycles
integer :: n_active ! # of active cycles
integer :: current_cycle = 0 ! current cycle
integer(8) :: n_particles ! # of particles per generation
integer :: n_batches ! # of batches
integer :: n_inactive ! # of inactive batches
integer :: n_active ! # of active batches
integer :: gen_per_batch = 1 ! # of generations per batch
integer :: current_batch = 0 ! current batch
integer :: current_gen = 0 ! current generation
! Random Number seed
integer(8) :: seed = 1_8
@ -147,16 +150,19 @@ module global
integer(8) :: work ! number of particles per processor
integer(8) :: maxwork ! maximum number of particles per processor
! cycle keff
! single-genreation keff
real(8) :: keff = ONE
real(8) :: keff_std
! Shannon entropy
logical :: entropy_on = .false.
real(8) :: entropy ! value of shannon entropy
real(8), allocatable :: entropy_p(:) ! fraction of source sites in each cell
real(8) :: entropy ! value of shannon entropy
real(8), allocatable :: entropy_p(:,:,:) ! % of source sites in each cell
type(StructuredMesh), pointer :: entropy_mesh
! Write source at end of simulation
logical :: write_source = .false.
! ============================================================================
! PARALLEL PROCESSING VARIABLES
@ -212,6 +218,9 @@ module global
! Message used in message/warning/fatal_error
character(MAX_LINE_LEN) :: message
! Random number seed
integer(8) :: seed = 1_8
! Problem type
integer :: problem_type = PROB_CRITICALITY
@ -221,7 +230,8 @@ module global
! Trace for single particle
logical :: trace
integer :: trace_cycle
integer :: trace_batch
integer :: trace_gen
integer(8) :: trace_particle
! ============================================================================

17
src/hdf5_interface.F90
View file

@ -101,19 +101,22 @@ contains
call hdf5_make_double(hdf5_output_file, "n_particles", real(n_particles,8))
! Use H5LT interface to write n_cycles, n_inactive, and n_active
call hdf5_make_integer(hdf5_output_file, "n_cycles", n_cycles)
call hdf5_make_integer(hdf5_output_file, "n_batches", n_batches)
call hdf5_make_integer(hdf5_output_file, "n_inactive", n_inactive)
call hdf5_make_integer(hdf5_output_file, "n_active", n_active)
call hdf5_make_integer(hdf5_output_file, "gen_per_batch", gen_per_batch)
! Add description of each variable
call h5ltset_attribute_string_f(hdf5_output_file, "n_particles", &
"description", "Number of particles per cycle", hdf5_err)
call h5ltset_attribute_string_f(hdf5_output_file, "n_cycles", &
"description", "Total number of cycles", hdf5_err)
call h5ltset_attribute_string_f(hdf5_output_file, "n_batches", &
"description", "Total number of batches", hdf5_err)
call h5ltset_attribute_string_f(hdf5_output_file, "n_inactive", &
"description", "Number of inactive cycles", hdf5_err)
call h5ltset_attribute_string_f(hdf5_output_file, "n_active", &
"description", "Number of active cycles", hdf5_err)
call h5ltset_attribute_string_f(hdf5_output_file, "gen_per_batch", &
"description", "Number of generations per batch", hdf5_err)
end if
call hdf5_write_geometry()
@ -812,7 +815,7 @@ contains
call h5ltset_attribute_string_f(timing_group, "time_transport", &
"description", "Time in transport only (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_intercycle", &
"description", "Total time between cycles (s)", hdf5_err)
"description", "Total time between generations (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_tallies", &
"description", "Time between cycles accumulating tallies (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_sample", &
@ -820,16 +823,16 @@ contains
call h5ltset_attribute_string_f(timing_group, "time_sendrecv", &
"description", "Time between cycles SEND/RECVing source sites (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_inactive", &
"description", "Total time in inactive cycles (s)", hdf5_err)
"description", "Total time in inactive batches (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_active", &
"description", "Total time in active cycles (s)", hdf5_err)
"description", "Total time in active batches (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_finalize", &
"description", "Total time for finalization (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_total", &
"description", "Total time elapsed (s)", hdf5_err)
! Write calculation rate
total_particles = n_particles * n_cycles
total_particles = n_particles * n_batches * gen_per_batch
speed = real(total_particles) / (time_inactive % elapsed + &
time_active % elapsed)
call hdf5_make_double(timing_group, "neutrons_per_second", speed)

28
src/initialize.F90
View file

@ -3,23 +3,24 @@ module initialize
use ace, only: read_xs
use bank_header, only: Bank
use constants
use datatypes, only: dict_create, dict_add_key, dict_get_key, &
dict_has_key, dict_keys
use datatypes, only: dict_create, dict_get_key, dict_has_key, &
dict_keys
use datatypes_header, only: ListKeyValueII, DictionaryII
use energy_grid, only: unionized_grid
use error, only: fatal_error
use geometry, only: neighbor_lists
use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE
use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE
use global
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict
cells_in_univ_dict, read_plots_xml
use output, only: title, header, print_summary, print_geometry, &
print_plot, create_summary_file, &
create_xs_summary_file
use random_lcg, only: initialize_prng
use source, only: initialize_source
use string, only: to_str, starts_with, ends_with, lower_case
use tally, only: create_tally_map, TallyObject
use tally, only: create_tally_map
use tally_header, only: TallyObject
use timing, only: timer_start, timer_stop
#ifdef MPI
@ -33,8 +34,6 @@ module initialize
implicit none
type(DictionaryII), pointer :: build_dict => null()
contains
!===============================================================================
@ -81,6 +80,10 @@ contains
! Initialize random number generator
call initialize_prng()
! Read plots.xml if it exists -- this has to be done separate from the other
! XML files because we need the PRNG to be initialized first
if (plotting) call read_plots_xml()
! Set up universe structures
call prepare_universes()
@ -106,9 +109,11 @@ contains
call timer_stop(time_read_xs)
! Construct unionized energy grid from cross-sections
call timer_start(time_unionize)
call unionized_grid()
call timer_stop(time_unionize)
if (grid_method == GRID_UNION) then
call timer_start(time_unionize)
call unionized_grid()
call timer_stop(time_unionize)
end if
! Create tally map
call create_tally_map()
@ -282,9 +287,6 @@ contains
! Create special dictionary used in input_xml
call dict_create(cells_in_univ_dict)
! Create special dictionary for building universes
call dict_create(build_dict)
end subroutine create_dictionaries

287
src/input_xml.F90
View file

@ -2,8 +2,7 @@ module input_xml
use cmfd_input, only: read_cmfd_xml
use constants
use datatypes, only: dict_create, dict_add_key, dict_has_key, &
dict_get_key
use datatypes, only: dict_add_key, dict_has_key, dict_get_key
use error, only: fatal_error, warning
use geometry_header, only: Cell, Surface, Lattice
use global
@ -34,7 +33,6 @@ contains
call read_materials_xml()
call read_tallies_xml()
call read_cmfd_xml()
if (plotting) call read_plots_xml()
end subroutine read_input_xml
@ -69,6 +67,10 @@ contains
! Initialize XML scalar variables
cross_sections_ = ""
verbosity_ = 0
energy_grid_ = "union"
seed_ = 0_8
write_source_ = ""
source_ % type = ""
! Parse settings.xml file
call read_xml_file_settings_t(filename)
@ -93,7 +95,7 @@ contains
end if
! Criticality information
if (criticality % cycles > 0) then
if (criticality % batches > 0) then
problem_type = PROB_CRITICALITY
! Check number of particles
@ -104,41 +106,71 @@ contains
n_particles = str_to_int(criticality % particles)
! Copy cycle information
n_cycles = criticality % cycles
n_inactive = criticality % inactive
n_active = n_cycles - n_inactive
n_batches = criticality % batches
n_inactive = criticality % inactive
n_active = n_batches - n_inactive
gen_per_batch = criticality % generations_per_batch
else
message = "Need to specify number of batches with <batches> tag."
call fatal_error()
end if
<<<<<<< HEAD:src/input_xml.F90
! Random number seed
if (dble(criticality % seed) > 0.0_8) seed = dble(criticality % seed)
=======
! Copy random number seed if specified
if (seed_ > 0) seed = seed_
! Energy grid methods
select case (energy_grid_)
case ('nuclide')
grid_method = GRID_NUCLIDE
case ('union')
grid_method = GRID_UNION
case ('lethargy')
message = "Lethargy mapped energy grid not yet supported."
call fatal_error()
case default
message = "Unknown energy grid method: " // energy_grid_
call fatal_error()
end select
>>>>>>> master:src/input_xml.F90
! Verbosity
if (verbosity_ > 0) verbosity = verbosity_
if (associated(source_ % coeffs)) then
if (len(source_ % type) > 0) then
! Determine external source type
type = source_ % type
call lower_case(type)
select case (trim(type))
select case (type)
case ('box')
external_source % type = SRC_BOX
coeffs_reqd = 6
case ('point')
external_source % type = SRC_POINT
coeffs_reqd = 3
case ('file')
external_source % type = SRC_FILE
case default
message = "Invalid source type: " // trim(source_ % type)
call fatal_error()
end select
! Coefficients for external surface
n = size(source_ % coeffs)
if (n < coeffs_reqd) then
message = "Not enough coefficients specified for external source."
call fatal_error()
elseif (n > coeffs_reqd) then
message = "Too many coefficients specified for external source."
call fatal_error()
else
allocate(external_source % values(n))
external_source % values = source_ % coeffs
if (type /= 'file') then
n = size(source_ % coeffs)
if (n < coeffs_reqd) then
message = "Not enough coefficients specified for external source."
call fatal_error()
elseif (n > coeffs_reqd) then
message = "Too many coefficients specified for external source."
call fatal_error()
else
allocate(external_source % values(n))
external_source % values = source_ % coeffs
end if
end if
end if
@ -156,8 +188,9 @@ contains
! Particle trace
if (associated(trace_)) then
trace_cycle = trace_(1)
trace_particle = trace_(2)
trace_batch = trace_(1)
trace_gen = trace_(2)
trace_particle = trace_(3)
end if
! Shannon Entropy mesh
@ -210,6 +243,9 @@ contains
entropy_on = .true.
end if
! Check if the user has specified to write binary source file
if (trim(write_source_) == 'on') write_source = .true.
end subroutine read_settings_xml
!===============================================================================
@ -266,14 +302,6 @@ contains
c % material = cell_(i) % material
c % fill = cell_(i) % fill
! Set plot color
if (plotting) then
allocate(c % rgb(3))
c % rgb(1) = prn()*255
c % rgb(2) = prn()*255
c % rgb(3) = prn()*255
end if
! Check to make sure that either material or fill was specified
if (c % material == 0 .and. c % fill == 0) then
message = "Neither material nor fill was specified for cell " // &
@ -505,6 +533,7 @@ contains
integer :: n ! number of nuclides
real(8) :: val ! value entered for density
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
character(3) :: default_xs ! default xs identifier (e.g. 70c)
character(12) :: name ! name of nuclide
character(MAX_WORD_LEN) :: units ! units on density
@ -544,40 +573,43 @@ contains
! Copy material id
m % id = material_(i) % id
! Set plot color
if (plotting) then
allocate(m % rgb(3))
m % rgb(1) = prn()*255
m % rgb(2) = prn()*255
m % rgb(3) = prn()*255
end if
! Copy density -- the default value for the units is given in the
! material_t.xml file and doesn't need to be specified here, hence case
! default results in an error.
val = material_(i) % density % value
if (val <= ZERO) then
message = "Need to specify a positive density on material " // &
trim(to_str(m % id)) // "."
call fatal_error()
end if
! Copy value and units
val = material_(i) % density % value
units = material_(i) % density % units
call lower_case(units)
select case(trim(units))
case ('g/cc', 'g/cm3')
m % density = -val
case ('kg/m3')
m % density = -0.001 * val
case ('atom/b-cm')
m % density = val
case ('atom/cm3', 'atom/cc')
m % density = 1.0e-24 * val
case default
message = "Unkwown units '" // trim(material_(i) % density % units) &
// "' specified on material " // trim(to_str(m % id))
call fatal_error()
end select
if (units == 'sum') then
! If the user gave the units as 'sum', then the total density of the
! material is taken to be the sum of the atom fractions listed on the
! nuclides
sum_density = .true.
else
! Check for erroneous density
sum_density = .false.
if (val <= ZERO) then
message = "Need to specify a positive density on material " // &
trim(to_str(m % id)) // "."
call fatal_error()
end if
! Adjust material density based on specified units
call lower_case(units)
select case(trim(units))
case ('g/cc', 'g/cm3')
m % density = -val
case ('kg/m3')
m % density = -0.001 * val
case ('atom/b-cm')
m % density = val
case ('atom/cm3', 'atom/cc')
m % density = 1.0e-24 * val
case default
message = "Unkwown units '" // trim(material_(i) % density % units) &
// "' specified on material " // trim(to_str(m % id))
call fatal_error()
end select
end if
! Check to ensure material has at least one nuclide
if (.not. associated(material_(i) % nuclides)) then
@ -652,6 +684,9 @@ contains
end if
end do
! Determine density if it is a sum value
if (sum_density) m % density = sum(m % atom_percent)
! Add material to dictionary
call dict_add_key(material_dict, m % id, i)
@ -986,6 +1021,11 @@ contains
case ('scatter-3')
t % score_bins(j) % scalar = SCORE_SCATTER_3
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('diffusion')
t % score_bins(j) % scalar = SCORE_DIFFUSION
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n1n')
@ -1116,10 +1156,11 @@ contains
use xml_data_plots_t
integer i
integer i, j
integer n_cols, col_id
logical :: file_exists ! does plots.xml file exist?
character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
type(Plot), pointer :: pl => null()
type(Plot), pointer :: pl => null()
! Check if plots.xml exists
filename = trim(path_input) // "plots.xml"
@ -1145,26 +1186,33 @@ contains
! Copy data into plots
pl % id = plot_(i) % id
pl % aspect = plot_(i) % aspect
! Set output file path
pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
"_" // trim(plot_(i) % filename) // ".ppm"
! Copy plot pixel size
if (size(plot_(i) % pixels) == 2) then
pl % pixels = plot_(i) % pixels
else
message = "<pixels> must be length 2 in plot " // to_str(i)
message = "<pixels> must be length 2 in plot " // to_str(pl % id)
call fatal_error()
end if
select case (plot_(i) % color)
case ("cell")
pl % color = PLOT_COLOR_CELLS
case ("mat")
pl % color = PLOT_COLOR_MATS
case default
message = "Unsupported plot color '" // plot_(i) % color &
// "' in plot " // trim(to_str(i))
call fatal_error()
end select
! Copy plot background color
if (associated(plot_(i) % background)) then
if (size(plot_(i) % background) == 3) then
pl % not_found % rgb = plot_(i) % background
else
message = "Bad background RGB " &
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
else
pl % not_found % rgb = (/ 255, 255, 255 /)
end if
! Copy plot type
select case (plot_(i) % type)
case ("slice")
pl % type = PLOT_TYPE_SLICE
@ -1172,10 +1220,11 @@ contains
! pl % type = PLOT_TYPE_POINTS
case default
message = "Unsupported plot type '" // plot_(i) % type &
// "' in plot " // trim(to_str(i))
// "' in plot " // trim(to_str(pl % id))
call fatal_error()
end select
! Copy plot basis
select case (plot_(i) % basis)
case ("xy")
pl % basis = PLOT_BASIS_XY
@ -1185,7 +1234,7 @@ contains
pl % basis = PLOT_BASIS_YZ
case default
message = "Unsupported plot basis '" // plot_(i) % basis &
// "' in plot " // trim(to_str(i))
// "' in plot " // trim(to_str(pl % id))
call fatal_error()
end select
@ -1194,7 +1243,7 @@ contains
pl % origin = plot_(i) % origin
else
message = "Origin must be length 3 " &
// "in plot " // trim(to_str(i))
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
@ -1206,10 +1255,88 @@ contains
pl % width(2) = plot_(i) % width(2)
else
message = "Bad plot width " &
// "in plot " // trim(to_str(i))
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
! Copy plot color type and initialize all colors randomly
select case (plot_(i) % color)
case ("cell")
pl % color_by = PLOT_COLOR_CELLS
allocate(pl % colors(n_cells))
do j = 1, n_cells
pl % colors(j) % rgb(1) = prn()*255
pl % colors(j) % rgb(2) = prn()*255
pl % colors(j) % rgb(3) = prn()*255
end do
case ("mat", "material")
pl % color_by = PLOT_COLOR_MATS
allocate(pl % colors(n_materials))
do j = 1, n_materials
pl % colors(j) % rgb(1) = prn()*255
pl % colors(j) % rgb(2) = prn()*255
pl % colors(j) % rgb(3) = prn()*255
end do
case default
message = "Unsupported plot color type '" // plot_(i) % color &
// "' in plot " // trim(to_str(pl % id))
call fatal_error()
end select
! Copy user specified colors
if (associated(plot_(i) % col_spec_)) then
n_cols = size(plot_(i) % col_spec_)
do j = 1, n_cols
if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then
message = "Bad RGB " &
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
col_id = plot_(i) % col_spec_(j) % id
if (pl % color_by == PLOT_COLOR_CELLS) then
if (dict_has_key(cell_dict, col_id)) then
pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb
else
message = "Could not find cell " // trim(to_str(col_id)) // &
" specified in plot " // trim(to_str(pl % id))
call fatal_error()
end if
else if (pl % color_by == PLOT_COLOR_MATS) then
if (dict_has_key(material_dict, col_id)) then
pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb
else
message = "Could not find material " // trim(to_str(col_id)) // &
" specified in plot " // trim(to_str(pl % id))
call fatal_error()
end if
end if
end do
end if
! Alter colors based on mask information
if (associated(plot_(i) % mask_)) then
if (size(plot_(i) % mask_) > 1) then
message = "Mutliple masks" // &
" specified in plot " // trim(to_str(pl % id))
call fatal_error()
else if (.not. size(plot_(i) % mask_) == 0) then
do j=1,size(pl % colors)
if (.not. any(j .eq. plot_(i) % mask_(1) % components)) then
pl % colors(j) % rgb = plot_(i) % mask_(1) % background
end if
end do
end if
end if
end do

167
src/intercycle.F90
View file

@ -4,13 +4,13 @@ module intercycle
use error, only: fatal_error, warning
use global
use mesh, only: get_mesh_bin
use mesh, only: count_fission_sites
use mesh_header, only: StructuredMesh
use output, only: write_message
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use string, only: to_str
use tally, only: accumulate_cycle_estimate
use tally, only: accumulate_batch_estimate
use tally_header, only: TallyObject
use timing, only: timer_start, timer_stop
@ -22,8 +22,8 @@ contains
!===============================================================================
! SYNCHRONIZE_BANK samples source sites from the fission sites that were
! accumulated during the cycle. This routine is what allows this Monte Carlo to
! scale to large numbers of processors where other codes cannot.
! accumulated during the generation. This routine is what allows this Monte
! Carlo to scale to large numbers of processors where other codes cannot.
!===============================================================================
subroutine synchronize_bank()
@ -73,7 +73,7 @@ contains
total = n_bank
#endif
! If there are not that many particles per cycle, it's possible that no
! If there are not that many particles per generation, it's possible that no
! fission sites were created at all on a single processor. Rather than add
! extra logic to treat this circumstance, we really want to ensure the user
! runs enough particles to avoid this in the first place.
@ -83,12 +83,12 @@ contains
call fatal_error()
end if
! Make sure all processors start at the same point for random sampling by
! using the current_cycle as the starting seed. Then skip ahead in the
! sequence using the starting index in the 'global' fission bank for each
! processor Skip ahead however many random numbers are needed
! Make sure all processors start at the same point for random sampling. Then
! skip ahead in the sequence using the starting index in the 'global'
! fission bank for each processor.
call set_particle_seed(int(current_cycle,8))
call set_particle_seed(int((current_batch - 1)*gen_per_batch + &
current_gen,8))
call prn_skip(start)
! Determine how many fission sites we need to sample from the source bank
@ -257,20 +257,22 @@ contains
! Send we initiated a series of asynchronous ISENDs and IRECVs, now we have
! to ensure that the data has actually been communicated before moving on to
! the next cycle
! the next generation
call MPI_WAITALL(n_request, request, MPI_STATUSES_IGNORE, mpi_err)
! Deallocate space for bank_position on the last cycle
if (current_cycle == n_cycles) deallocate(bank_position)
! Deallocate space for bank_position on the very last generation
if (current_batch == n_batches .and. current_gen == gen_per_batch) &
deallocate(bank_position)
#else
source_bank = temp_sites(1:n_particles)
#endif
call timer_stop(time_ic_sendrecv)
! Deallocate space for the temporary source bank on the last cycle
if (current_cycle == n_cycles) deallocate(temp_sites)
! Deallocate space for the temporary source bank on the last generation
if (current_batch == n_batches .and. current_gen == gen_per_batch) &
deallocate(temp_sites)
end subroutine synchronize_bank
@ -281,12 +283,10 @@ contains
subroutine shannon_entropy()
integer :: i ! index for bank sites
integer :: i, j, k ! index for bank sites
integer :: n ! # of boxes in each dimension
integer :: bin ! index in entropy_p
integer(8) :: total_bank ! total # of fission bank sites
integer, save :: n_box ! total # of boxes on mesh
logical :: outside_box ! were there sites outside entropy box?
real(8) :: total ! total weight of fission bank sites
logical :: sites_outside ! were there sites outside entropy box?
type(StructuredMesh), pointer :: m => null()
! Get pointer to entropy mesh
@ -310,96 +310,61 @@ contains
m % dimension = n
end if
! Determine total number of mesh boxes
n_box = product(m % dimension)
! allocate and determine width
allocate(m % width(3))
m % width = (m % upper_right - m % lower_left) / m % dimension
! allocate p
allocate(entropy_p(n_box))
allocate(entropy_p(m % dimension(1), m % dimension(2), &
m % dimension(3)))
end if
! initialize p
entropy_p = ZERO
outside_box = .false.
! loop over fission sites and count how many are in each mesh box
FISSION_SITES: do i = 1, int(n_bank,4)
! determine scoring bin for entropy mesh
call get_mesh_bin(m, fission_bank(i) % xyz, bin)
! if outside mesh, skip particle
if (bin == NO_BIN_FOUND) then
outside_box = .true.
cycle
end if
! add to appropriate mesh box
entropy_p(bin) = entropy_p(bin) + 1
end do FISSION_SITES
! count number of fission sites over mesh
call count_fission_sites(m, entropy_p, total, sites_outside)
! display warning message if there were sites outside entropy box
if (outside_box) then
if (sites_outside) then
message = "Fission source site(s) outside of entropy box."
call warning()
end if
#ifdef MPI
! collect values from all processors
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
else
call MPI_REDUCE(entropy_p, entropy_p, n_box, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
end if
! determine total number of bank sites
call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
#else
total_bank = n_bank
#endif
! sum values to obtain shannon entropy
if (master) then
! Normalize to number of bank sites
entropy_p = entropy_p / total_bank
! Normalize to total weight of bank sites
entropy_p = entropy_p / total
entropy = 0
do i = 1, n_box
if (entropy_p(i) > 0) then
entropy = entropy - entropy_p(i) * log(entropy_p(i))/log(2.0)
end if
do i = 1, m % dimension(1)
do j = 1, m % dimension(2)
do k = 1, m % dimension(3)
if (entropy_p(i,j,k) > ZERO) then
entropy = entropy - entropy_p(i,j,k) * &
log(entropy_p(i,j,k))/log(2.0)
end if
end do
end do
end do
end if
end subroutine shannon_entropy
!===============================================================================
! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the
! mean and standard deviation of the mean for active cycles and displays them
! CALCULATE_KEFF calculates the single batch estimate of keff as well as the
! mean and standard deviation of the mean for active batches and displays them
!===============================================================================
subroutine calculate_keff()
integer :: n ! active cycle number
real(8) :: k_cycle ! single cycle estimate of keff
integer :: n ! active batch number
real(8) :: k_batch ! single batch estimate of keff
#ifdef MPI
real(8) :: global_temp(N_GLOBAL_TALLIES)
#endif
message = "Calculate cycle keff..."
message = "Calculate batch keff..."
call write_message(8)
! Since the creation of bank sites was originally weighted by the last
! cycle keff, we need to multiply by that keff to get the current cycle's
! value
#ifdef MPI
global_tallies(K_ANALOG) % value = n_bank * keff
! Copy global tallies into array to be reduced
global_temp = global_tallies(:) % value
@ -416,59 +381,65 @@ contains
! Reset value on other processors
global_tallies(:) % value = ZERO
end if
#else
global_tallies(K_ANALOG) % value = n_bank * keff
#endif
! Collect statistics and print output
if (master) then
k_cycle = global_tallies(K_ANALOG) % value/n_particles
k_batch = global_tallies(K_ANALOG) % value/(n_particles*gen_per_batch)
if (current_cycle > n_inactive) then
! Active cycle number
n = current_cycle - n_inactive
if (tallies_on) then
! Active batch number
n = current_batch - n_inactive
! Accumulate single cycle realizations of k
call accumulate_cycle_estimate(global_tallies)
! Accumulate single batch realizations of k
call accumulate_batch_estimate(global_tallies)
! Determine mean and standard deviation of mean
! Determine sample mean of keff
keff = global_tallies(K_ANALOG) % sum/n
keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - keff*keff)/n)
! Display output for this cycle
if (current_cycle > n_inactive + 1) then
! Display output for this batch
if (current_batch > n_inactive + 1) then
! Determine standard deviation of the sample mean of keff
keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - &
keff*keff)/(n - 1))
if (entropy_on) then
if (cmfd_on) then
write(UNIT=OUTPUT_UNIT, FMT=104) current_cycle, k_cycle, &
write(UNIT=OUTPUT_UNIT, FMT=104) current_batch, k_batch, &
entropy, keff, keff_std, cmfd%keff, cmfd%entropy
else
write(UNIT=OUTPUT_UNIT, FMT=103) current_cycle, k_cycle, &
write(UNIT=OUTPUT_UNIT, FMT=103) current_batch, k_batch, &
entropy, keff, keff_std
end if
else
if (cmfd_on) then
write(UNIT=OUTPUT_UNIT, FMT=105) current_cycle, k_cycle, &
write(UNIT=OUTPUT_UNIT, FMT=105) current_batch, k_batch, &
keff, keff_std, cmfd%keff
else
write(UNIT=OUTPUT_UNIT, FMT=101) current_cycle, k_cycle, &
write(UNIT=OUTPUT_UNIT, FMT=101) current_batch, k_batch, &
keff, keff_std
end if
end if
else
if (entropy_on) then
write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy
write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy
else
write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle
write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch
end if
end if
else
! Display output for inactive cycle
! Display output for inactive batch
if (entropy_on) then
write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy
write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy
else
write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle
write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch
end if
keff = k_cycle
! Set keff
keff = k_batch
! Reset tally values
global_tallies(:) % value = ZERO
end if
end if

85
src/main.F90
View file

@ -9,10 +9,9 @@ program main
use output, only: write_message, header
use plotter, only: run_plot
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: get_source_particle
use string, only: to_str
use tally, only: synchronize_tallies
use tally, only: synchronize_tallies, add_to_score
use timing, only: timer_start, timer_stop
#ifdef MPI
@ -55,25 +54,26 @@ contains
! Display column titles
if (entropy_on) then
<<<<<<< HEAD:src/main.F90
if (cmfd_on) then
message = " Cycle k(cycle) Entropy Average k CMFD k CMFD Ent."
message = " Cycle k(batch) Entropy Average k CMFD k CMFD Ent."
call write_message(1)
message = " ===== ======== ======= =================== ========== ========="
call write_message(1)
else
message = " Cycle k(cycle) Entropy Average k"
message = " Cycle k(batch) Entropy Average k"
call write_message(1)
message = " ===== ======== ======= ==================="
call write_message(1)
end if
else
if (cmfd_on) then
message = " Cycle k(cycle) Average k CMFD k"
message = " Cycle k(batch) Average k CMFD k"
call write_message(1)
message = " ===== ======== =================== =========="
call write_message(1)
else
message = " Cycle k(cycle) Average k"
message = " Cycle k(batch) Average k"
call write_message(1)
message = " ===== ======== ==================="
call write_message(1)
@ -81,39 +81,56 @@ contains
end if
! ==========================================================================
! LOOP OVER CYCLES
CYCLE_LOOP: do current_cycle = 1, n_cycles
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_batches
message = "Simulating cycle " // trim(to_str(current_cycle)) // "..."
message = "Simulating batch " // trim(to_str(current_batch)) // "..."
call write_message(8)
! Set all tallies to zero
n_bank = 0
! =======================================================================
! LOOP OVER HISTORIES
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
! Start timer for transport
call timer_start(time_transport)
! Set all tallies to zero
n_bank = 0
HISTORY_LOOP: do i = 1, work
! ====================================================================
! LOOP OVER HISTORIES
! grab source particle from bank
call get_source_particle(i)
! Start timer for transport
call timer_start(time_transport)
! transport particle
call transport()
HISTORY_LOOP: do i = 1, work
end do HISTORY_LOOP
! grab source particle from bank
call get_source_particle(i)
! Accumulate time for transport
call timer_stop(time_transport)
! transport particle
call transport()
! =======================================================================
! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC
end do HISTORY_LOOP
! Start timer for inter-cycle synchronization
call timer_start(time_intercycle)
! Accumulate time for transport
call timer_stop(time_transport)
! ====================================================================
! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC
! Start timer for inter-cycle synchronization
call timer_start(time_intercycle)
! Distribute fission bank across processors evenly
call synchronize_bank()
! Add to analog estimate of keff -- since the creation of bank sites
! was originally weighted by the last cycle keff, we need to multiply
! by that keff to get the current cycle's value
call add_to_score(global_tallies(K_ANALOG), n_bank * keff)
! Stop timer for inter-cycle synchronization
call timer_stop(time_intercycle)
end do GENERATION_LOOP
! Collect tallies
if (tallies_on) then
@ -125,30 +142,22 @@ contains
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Distribute fission bank across processors evenly
call synchronize_bank()
! Collect results and statistics
call calculate_keff()
! print cycle information
! run cmfd
if (current_cycle > n_inactive .and. cmfd_on) then
if (current_batch > n_inactive .and. cmfd_on) then
call execute_cmfd()
end if
! Turn tallies on once inactive cycles are complete
if (current_cycle == n_inactive) then
if (current_batch == n_inactive) then
tallies_on = .true.
call timer_stop(time_inactive)
call timer_start(time_active)
end if
! Stop timer for inter-cycle synchronization
call timer_stop(time_intercycle)
end do CYCLE_LOOP
end do BATCH_LOOP
call timer_stop(time_active)

2
src/material_header.F90
View file

@ -21,8 +21,6 @@ module material_header
character(12) :: sab_name ! name of S(a,b) table
integer :: sab_table = 0 ! index in sab_tables
integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table
integer, allocatable :: rgb(:) ! plotting color
end type Material
end module material_header

175
src/mesh.F90
View file

@ -1,9 +1,14 @@
module mesh
use constants
use global
use mesh_header
use particle_header, only: Particle
#ifdef MPI
use mpi
#endif
implicit none
contains
@ -137,4 +142,174 @@ contains
end subroutine bin_to_mesh_indices
!===============================================================================
! COUNT_FISSION_SITES determines the number of fission bank sites in each cell
! of a given mesh. This can be used for a variety of purposes (Shannon entropy,
! CMFD, uniform fission source weighting)
!===============================================================================
subroutine count_fission_sites(m, count, total, sites_outside)
type(StructuredMesh), pointer :: m ! mesh to count sites
real(8), intent(out) :: count(:,:,:) ! weight of sites in each cell
real(8), intent(out) :: total ! total weight of sites
logical, optional :: sites_outside ! were there sites outside mesh?
integer :: i ! loop index for local fission sites
integer :: ijk(3) ! indices on mesh
real(8) :: weight ! accumulated weight of sites
logical :: in_mesh ! was single site outside mesh?
logical :: outside ! was any site outside mesh?
#ifdef MPI
integer :: n ! total size of count variable
#endif
! initialize variables
count = ZERO
weight = ZERO
outside = .false.
! loop over fission sites and count how many are in each mesh box
FISSION_SITES: do i = 1, int(n_bank,4)
! determine scoring bin for entropy mesh
call get_mesh_indices(m, fission_bank(i) % xyz, ijk, in_mesh)
! if outside mesh, skip particle
if (.not. in_mesh) then
outside = .true.
cycle
end if
! add weight
weight = weight + ONE
! add to appropriate mesh box
! TODO: if tracking weight through bank, add weight instead
count(ijk(1),ijk(2),ijk(3)) = count(ijk(1),ijk(2),ijk(3)) + 1
end do FISSION_SITES
#ifdef MPI
! determine total number of mesh cells
n = size(count,1) * size(count,2) * size(count,3)
! collect values from all processors
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, count, n, MPI_REAL8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
else
call MPI_REDUCE(count, count, n, MPI_REAL8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
end if
! Check if there were sites outside the mesh for any processor
if (present(sites_outside)) then
call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, &
MPI_COMM_WORLD, mpi_err)
end if
! determine total weight of bank sites
call MPI_REDUCE(weight, total, 1, MPI_REAL8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
#else
total = weight
sites_outside = outside
#endif
end subroutine count_fission_sites
!===============================================================================
! MESH_INTERSECT determines if a line between xyz0 and xyz1 intersects the outer
! boundary of the given mesh. This is important for determining whether a track
! will score to a mesh tally.
!===============================================================================
function mesh_intersects(m, xyz0, xyz1) result(intersects)
type(StructuredMesh), pointer :: m
real(8), intent(in) :: xyz0(3)
real(8), intent(in) :: xyz1(3)
logical :: intersects
real(8) :: x0, y0, z0 ! track start point
real(8) :: x1, y1, z1 ! track end point
real(8) :: xi, yi, zi ! track intersection point with mesh
real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh
real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
y0 = xyz0(2)
z0 = xyz0(3)
! Copy coordinates of ending point
x1 = xyz1(1)
y1 = xyz1(2)
z1 = xyz1(3)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
ym0 = m % lower_left(2)
zm0 = m % lower_left(3)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
ym1 = m % upper_right(2)
zm1 = m % upper_right(3)
! Check if line intersects bottom surface -- calculate the intersection
! point (y,z)
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects left surface -- calculate the intersection point
! (x,z)
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects front surface -- calculate the intersection point
! (x,y)
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
! Check if line intersects top surface -- calculate the intersection
! point (y,z)
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects right surface -- calculate the intersection point
! (x,z)
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects back surface -- calculate the intersection point
! (x,y)
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end function mesh_intersects
end module mesh

14
src/mesh_header.F90
View file

@ -8,13 +8,13 @@ module mesh_header
!===============================================================================
type StructuredMesh
integer :: id
integer :: type
integer :: n_dimension
integer, allocatable :: dimension(:)
real(8), allocatable :: lower_left(:)
real(8), allocatable :: upper_right(:)
real(8), allocatable :: width(:)
integer :: id ! user-specified id
integer :: type ! rectangular, hexagonal
integer :: n_dimension ! rank of mesh
integer, allocatable :: dimension(:) ! number of cells in each direction
real(8), allocatable :: lower_left(:) ! lower-left corner of mesh
real(8), allocatable :: upper_right(:) ! upper-right corner of mesh
real(8), allocatable :: width(:) ! width of each mesh cell
end type StructuredMesh
end module mesh_header

217
src/output.F90
View file

@ -29,8 +29,8 @@ contains
subroutine title()
character(10) :: today_date
character(8) :: today_time
character(10) :: date_
character(10) :: time_
write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
' .d88888b. 888b d888 .d8888b.', &
@ -50,20 +50,28 @@ contains
' Developed At: Massachusetts Institute of Technology'
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') &
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1
#endif
! Write the date and time
call get_today(today_date, today_time)
call date_and_time(DATE=date_, TIME=time_)
date_ = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8)
time_ = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') &
trim(today_date), trim(today_time)
trim(date_), trim(time_)
! Write information to summary file
call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1)
write(UNIT=UNIT_SUMMARY, FMT=*) &
"Copyright: 2011 Massachusetts Institute of Technology"
"Copyright: 2011-2012 Massachusetts Institute of Technology"
write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') &
"Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1
#endif
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') &
trim(today_date), trim(today_time)
trim(date_), trim(time_)
! Write information on number of processors
#ifdef MPI
@ -86,8 +94,8 @@ contains
integer, optional :: unit ! unit to write to
integer, optional :: level ! specified header level
integer :: n
integer :: m
integer :: n ! number of = signs on left
integer :: m ! number of = signs on right
integer :: unit_ ! unit to write to
integer :: header_level ! actual header level
character(MAX_LINE_LEN) :: line
@ -137,14 +145,15 @@ contains
subroutine write_message(level)
integer, optional :: level
integer, optional :: level ! verbosity level
integer :: n_lines
integer :: i
integer :: n_lines ! number of lines needed
integer :: i ! index for lines
! Only allow master to print to screen
if (.not. master .and. present(level)) return
! TODO: Take care of line wrapping so words don't get cut off
if (.not. present(level) .or. level <= verbosity) then
n_lines = (len_trim(message)-1)/79 + 1
do i = 1, n_lines
@ -154,55 +163,13 @@ contains
end subroutine write_message
!===============================================================================
! GET_TODAY determines the date and time at which the program began execution
! and returns it in a readable format
!===============================================================================
subroutine get_today(today_date, today_time)
character(10), intent(out) :: today_date
character(8), intent(out) :: today_time
integer :: val(8)
character(8) :: date_
character(10) :: time_
character(5) :: zone
call date_and_time(date_, time_, zone, val)
! val(1) = year (YYYY)
! val(2) = month (MM)
! val(3) = day (DD)
! val(4) = timezone
! val(5) = hours (HH)
! val(6) = minutes (MM)
! val(7) = seconds (SS)
! val(8) = milliseconds
if (val(2) < 10) then
if (val(3) < 10) then
today_date = date_(6:6) // "/" // date_(8:8) // "/" // date_(1:4)
else
today_date = date_(6:6) // "/" // date_(7:8) // "/" // date_(1:4)
end if
else
if (val(3) < 10) then
today_date = date_(5:6) // "/" // date_(8:8) // "/" // date_(1:4)
else
today_date = date_(5:6) // "/" // date_(7:8) // "/" // date_(1:4)
end if
end if
today_time = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
end subroutine get_today
!===============================================================================
! PRINT_PARTICLE displays the attributes of a particle
!===============================================================================
subroutine print_particle()
integer :: i
integer :: i ! index for coordinate levels
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
@ -263,6 +230,7 @@ contains
write(ou,*) ' Surface = ' // to_str(sign(s % id, p % surface))
end if
! Display weight, energy, grid index, and interpolation factor
write(ou,*) ' Weight = ' // to_str(p % wgt)
write(ou,*) ' Energy = ' // to_str(p % E)
write(ou,*) ' IE = ' // to_str(p % IE)
@ -298,27 +266,35 @@ contains
subroutine print_cell(c, unit)
type(Cell), pointer :: c
integer, optional :: unit
integer, optional :: unit ! specified unit to write to
integer :: temp
integer :: i
integer :: unit_
integer :: index_cell ! index in cells array
integer :: i ! loop index for surfaces
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
! Set unit to stdout if not already set
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write user-specified id for cell
write(unit_,*) 'Cell ' // to_str(c % id)
temp = dict_get_key(cell_dict, c % id)
write(unit_,*) ' Array Index = ' // to_str(temp)
! Find index in cells array and write
index_cell = dict_get_key(cell_dict, c % id)
write(unit_,*) ' Array Index = ' // to_str(index_cell)
! Write what universe this cell is in
u => universes(c % universe)
write(unit_,*) ' Universe = ' // to_str(u % id)
! Write information on fill for cell
select case (c % type)
case (CELL_NORMAL)
write(unit_,*) ' Fill = NONE'
@ -329,12 +305,16 @@ contains
l => lattices(c % fill)
write(unit_,*) ' Fill = Lattice ' // to_str(l % id)
end select
! Write information on material
if (c % material == 0) then
write(unit_,*) ' Material = NONE'
else
m => materials(c % material)
write(unit_,*) ' Material = ' // to_str(m % id)
end if
! Write surface specification
string = ""
do i = 1, c % n_surfaces
select case (c % surfaces(i))
@ -364,24 +344,31 @@ contains
type(Universe), pointer :: univ
integer, optional :: unit
integer :: i
integer :: unit_
integer :: i ! loop index for cells in this universe
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Universe), pointer :: base_u
type(Universe), pointer :: base_u => null()
! Set default unit to stdout if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Get a pointer to the base universe
base_u => universes(BASE_UNIVERSE)
! Write user-specified id for this universe
write(unit_,*) 'Universe ' // to_str(univ % id)
! If this is the base universe, indicate so
if (associated(univ, base_u)) then
write(unit_,*) ' Base Universe'
end if
! Write list of cells in this universe
string = ""
do i = 1, univ % n_cells
c => cells(univ % cells(i))
@ -401,14 +388,16 @@ contains
type(Lattice), pointer :: lat
integer, optional :: unit
integer :: unit_
integer :: unit_ ! unit to write to
! set default unit if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write information about lattice
write(unit_,*) 'Lattice ' // to_str(lat % id)
write(unit_,*) ' n_x = ' // to_str(lat % n_x)
write(unit_,*) ' n_y = ' // to_str(lat % n_y)
@ -427,19 +416,23 @@ contains
subroutine print_surface(surf, unit)
type(Surface), pointer :: surf
integer, optional :: unit
integer, optional :: unit ! specified unit to write to
integer :: i
integer :: unit_
integer :: i ! loop index for coefficients
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
! set default unit if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write user-specified id of surface
write(unit_,*) 'Surface ' // to_str(surf % id)
! Write type of surface
select case (surf % type)
case (SURF_PX)
string = "X Plane"
@ -466,12 +459,14 @@ contains
end select
write(unit_,*) ' Type = ' // trim(string)
! Write coefficients for this surface
string = ""
do i = 1, size(surf % coeffs)
string = trim(string) // ' ' // to_str(surf % coeffs(i), 4)
end do
write(unit_,*) ' Coefficients = ' // trim(string)
! Write neighboring cells on positive side of this surface
string = ""
if (allocated(surf % neighbor_pos)) then
do i = 1, size(surf % neighbor_pos)
@ -480,6 +475,7 @@ contains
end if
write(unit_,*) ' Positive Neighbors = ' // trim(string)
! Write neighboring cells on negative side of this surface
string = ""
if (allocated(surf % neighbor_neg)) then
do i = 1, size(surf % neighbor_neg)
@ -487,6 +483,8 @@ contains
end do
end if
write(unit_,*) ' Negative Neighbors =' // trim(string)
! Write boundary condition for this surface
select case (surf % bc)
case (BC_TRANSMIT)
write(unit_,*) ' Boundary Condition = Transmission'
@ -510,12 +508,13 @@ contains
type(Material), pointer :: mat
integer, optional :: unit
integer :: i
integer :: unit_
real(8) :: density
integer :: i ! loop index for nuclides
integer :: unit_ ! unit to write to
real(8) :: density ! density in atom/b-cm
character(MAX_LINE_LEN) :: string
type(Nuclide), pointer :: nuc => null()
! set default unit to stdout if not specified
if (present(unit)) then
unit_ = unit
else
@ -556,9 +555,9 @@ contains
type(TallyObject), pointer :: t
integer, optional :: unit
integer :: i
integer :: id
integer :: unit_
integer :: i ! index for filter or score bins
integer :: id ! user-specified id
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
@ -566,14 +565,17 @@ contains
type(Material), pointer :: m => null()
type(StructuredMesh), pointer :: sm => null()
! set default unit to stdout if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write user-specified id of tally
write(unit_,*) 'Tally ' // to_str(t % id)
! Write any cells bins if present
if (t % n_filter_bins(FILTER_CELL) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_CELL)
@ -584,6 +586,7 @@ contains
write(unit_, *) ' Cell Bins:' // trim(string)
end if
! Write any surface bins if present
if (t % n_filter_bins(FILTER_SURFACE) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_SURFACE)
@ -594,6 +597,7 @@ contains
write(unit_, *) ' Surface Bins:' // trim(string)
end if
! Write any universe bins if present
if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_UNIVERSE)
@ -604,6 +608,7 @@ contains
write(unit_, *) ' Material Bins:' // trim(string)
end if
! Write any material bins if present
if (t % n_filter_bins(FILTER_MATERIAL) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_MATERIAL)
@ -614,6 +619,7 @@ contains
write(unit_, *) ' Material Bins:' // trim(string)
end if
! Write any mesh bins if present
if (t % n_filter_bins(FILTER_MESH) > 0) then
string = ""
id = t % mesh
@ -625,6 +631,7 @@ contains
write(unit_, *) ' Mesh Bins:' // trim(string)
end if
! Write any birth region bins if present
if (t % n_filter_bins(FILTER_CELLBORN) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_CELLBORN)
@ -635,6 +642,7 @@ contains
write(unit_, *) ' Birth Region Bins:' // trim(string)
end if
! Write any incoming energy bins if present
if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_ENERGYIN) + 1
@ -644,6 +652,7 @@ contains
write(unit_,*) ' Incoming Energy Bins:' // trim(string)
end if
! Write any outgoing energy bins if present
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_ENERGYOUT) + 1
@ -653,6 +662,7 @@ contains
write(unit_,*) ' Outgoing Energy Bins:' // trim(string)
end if
! Write any score bins if present
if (t % n_score_bins > 0) then
string = ""
do i = 1, t % n_score_bins
@ -671,7 +681,7 @@ contains
string = trim(string) // ' nu-fission'
end select
end do
write(unit_,*) ' Macro Reactions:' // trim(string)
write(unit_,*) ' Scores:' // trim(string)
end if
write(unit_,*)
@ -684,7 +694,7 @@ contains
subroutine print_geometry()
integer :: i
integer :: i ! loop index for various arrays
type(Surface), pointer :: s => null()
type(Cell), pointer :: c => null()
type(Universe), pointer :: u => null()
@ -732,12 +742,12 @@ contains
type(Nuclide), pointer :: nuc
integer, optional :: unit
integer :: i
integer :: unit_
integer :: size_total
integer :: size_xs
integer :: size_angle
integer :: size_energy
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_total ! memory used by nuclide (bytes)
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_angle ! memory used for angle distributions (bytes)
integer :: size_energy ! memory used for energy distributions (bytes)
type(Reaction), pointer :: rxn => null()
type(UrrData), pointer :: urr => null()
@ -821,7 +831,7 @@ contains
subroutine print_summary()
integer :: i
integer :: i ! loop index
character(15) :: string
type(Material), pointer :: m => null()
type(TallyObject), pointer :: t => null()
@ -830,8 +840,9 @@ contains
call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY)
if (problem_type == PROB_CRITICALITY) then
write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality'
write(UNIT_SUMMARY,101) 'Number of Cycles:', n_cycles
write(UNIT_SUMMARY,101) 'Number of Inactive Cycles:', n_inactive
write(UNIT_SUMMARY,101) 'Number of Batches:', n_batches
write(UNIT_SUMMARY,101) 'Number of Inactive Batches:', n_inactive
write(UNIT_SUMMARY,101) 'Generations per Batch:', gen_per_batch
elseif (problem_type == PROB_SOURCE) then
write(UNIT_SUMMARY,100) 'Problem type:', 'External Source'
end if
@ -892,41 +903,32 @@ contains
subroutine print_plot()
integer i
type(Plot), pointer :: pl => null()
integer :: i ! loop index for plots
type(Plot), pointer :: pl => null()
! Display header for plotting
call header("PLOTTING SUMMARY")
do i=1,n_plots
do i = 1, n_plots
pl => plots(i)
! Print plot id
! Write plot id
write(ou,100) "Plot ID:", trim(to_str(pl % id))
! Print plotting origin
! Write plotting origin
write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
! Print plotting width
! Write plotting width
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2)))
write(ou,100) "Coloring:", trim(to_str(pl % color))
write(ou,100) "Coloring:", trim(to_str(pl % color_by))
write(ou,100) "Basis:", trim(to_str(pl % basis))
write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
trim(to_str(pl % pixels(2)))
else if (pl % type == PLOT_TYPE_POINTS) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2))) // " " &
// trim(to_str(pl % width(3)))
write(ou,100) "Coloring:", trim(to_str(pl % color))
write(ou,100) "Ray Spacing:", trim(to_str(pl % aspect))
end if
write(ou,*)
@ -945,8 +947,8 @@ contains
subroutine print_runtime()
integer(8) :: total_particles
real(8) :: speed
integer(8) :: total_particles ! total # of particles simulated
real(8) :: speed ! # of neutrons/second
character(15) :: string
! display header block
@ -959,10 +961,16 @@ contains
write(ou,100) "Total time in simulation", time_inactive % elapsed + &
time_active % elapsed
write(ou,100) " Time in transport only", time_transport % elapsed
<<<<<<< HEAD:src/output.F90
write(ou,100) " Time in inactive cycles", time_inactive % elapsed
write(ou,100) " Time in active cycles", time_active % elapsed
write(ou,100) " Time between cycles", time_intercycle % elapsed
if(cmfd_on) write(ou,100) "Total CMFD time", time_cmfd % elapsed
=======
write(ou,100) " Time in inactive batches", time_inactive % elapsed
write(ou,100) " Time in active batches", time_active % elapsed
write(ou,100) " Time between generations", time_intercycle % elapsed
>>>>>>> master:src/output.F90
write(ou,100) " Accumulating tallies", time_ic_tallies % elapsed
write(ou,100) " Sampling source sites", time_ic_sample % elapsed
write(ou,100) " SEND/RECV source sites", time_ic_sendrecv % elapsed
@ -970,7 +978,7 @@ contains
write(ou,100) "Total time elapsed", time_total % elapsed
! display calculate rate and final keff
total_particles = n_particles * n_cycles
total_particles = n_particles * n_batches * gen_per_batch
speed = real(total_particles) / (time_inactive % elapsed + &
time_active % elapsed)
string = to_str(speed)
@ -1023,7 +1031,8 @@ contains
end subroutine create_summary_file
!===============================================================================
! CREATE_XS_SUMMARY_FILE
! CREATE_XS_SUMMARY_FILE creates an output file to write information about the
! cross section tables used in the simulation.
!===============================================================================
subroutine create_xs_summary_file()

3
src/particle_header.F90
View file

@ -52,7 +52,8 @@ module particle_header
real(8) :: last_E ! last energy
! What event last took place
integer :: event
integer :: event ! scatter, absorption, fission
integer :: event_nuclide ! index in nuclides array
! Post-collision physical data
integer :: n_bank ! number of fission sites banked

112
src/physics.F90
View file

@ -3,7 +3,7 @@ module physics
use ace_header, only: Nuclide, Reaction, DistEnergy
use constants
use cross_section, only: calculate_xs
use endf, only: reaction_name, is_fission, is_scatter
use endf, only: reaction_name
use error, only: fatal_error, warning
use fission, only: nu_total, nu_prompt, nu_delayed
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
@ -11,6 +11,7 @@ module physics
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use output, only: write_message
use particle_header, only: LocalCoord
use random_lcg, only: prn
@ -39,6 +40,11 @@ contains
logical :: found_cell ! found cell which particle is in?
type(LocalCoord), pointer :: coord => null()
if (verbosity >= 9 .or. trace) then
message = "Simulating Particle " // trim(to_str(p % id))
call write_message()
end if
if (p % coord % cell == NONE) then
call find_cell(found_cell)
! Particle couldn't be located
@ -51,17 +57,6 @@ contains
p % cell_born = p % coord % cell
end if
if (verbosity >= 9 .or. trace) then
message = "Simulating Particle " // trim(to_str(p % id))
call write_message()
end if
if (verbosity >= 10 .or. trace) then
message = " Born in cell " // trim(to_str(&
cells(p % coord % cell) % id))
call write_message()
end if
! Initialize number of events to zero
n_event = 0
@ -133,7 +128,7 @@ contains
p % surface = NONE
call collision()
! Save coordinates at collision for tallying purposes
! Save coordinates for tallying purposes
p % last_xyz = p % coord0 % xyz
! Set last material to none since cross sections will need to be
@ -188,9 +183,18 @@ contains
! Display information about collision
if (verbosity >= 10 .or. trace) then
!!$ message = " " // trim(reaction_name(MT)) // ". Energy = " // &
!!$ trim(to_str(p % E * 1e6_8)) // " eV."
!!$ call write_message()
select case (p % event)
case (EVENT_SCATTER)
message = " Scattered off of"
case (EVENT_FISSION)
message = " Fissioned in"
case (EVENT_ABSORB)
message = " Absorbed in"
end select
message = trim(message) // " " // trim(nuclides(&
p % event_nuclide) % name) // ". Energy = " // &
trim(to_str(p % E * 1e6_8)) // " eV."
call write_message()
end if
! check for very low energy
@ -268,6 +272,9 @@ contains
IE = micro_xs(index_nuclide) % index_grid
f = micro_xs(index_nuclide) % interp_factor
! Save which nuclide particle had collision with
p % event_nuclide = index_nuclide
! ==========================================================================
! DISAPPEARANCE REACTIONS (ANALOG) OR IMPLICIT CAPTURE (SURVIVAL BIASING)
@ -2055,77 +2062,4 @@ contains
end function watt_spectrum
!===============================================================================
! WIGNER samples a Wigner distribution of energy level spacings. Note that this
! scheme is C50 in the Monte Carlo Sampler from Los Alamos (LA-9721-MS).
!===============================================================================
function wigner(D_avg) result (D)
real(8), intent(in) :: D_avg ! average level spacing
real(8) :: D ! sampled level spacing
real(8) :: c
c = -4.*D_avg*D_avg/PI * log(prn())
D = sqrt(c)
end function wigner
!===============================================================================
! CHI_SQUARED samples a chi-squared distribution with n degrees of freedom. The
! distribution of resonance widths in the unresolved region is given by a
! chi-squared distribution. For the special case of n=1, this is a Porter-Thomas
! distribution. For cases with n odd, rule C64 is used whereas for cases with n
! even, rule C45 is used.
!===============================================================================
function chi_squared(n, G_avg) result(G)
integer, intent(in) :: n ! number of degrees of freedom
real(8), intent(in), optional :: G_avg ! average resonance width
integer :: i ! loop index
real(8) :: G ! sampled random variable (or resonance width)
real(8) :: x, y, c ! dummy variables
real(8) :: r1, r2 ! psuedorandom numbers
select case (mod(n,2))
case (0)
! Even number of degrees of freedom can be sampled via rule C45. We can
! sample x as -2/n*log(product(r_i, i = 1 to n/2))
x = ONE
do i = 1, n/2
x = x * prn()
end do
x = -2./n * log(x)
case (1)
! Odd number of degrees of freedom can be sampled via rule C64. We can
! sample x as -2/n*(log(r)*cos^2(pi/2*r) + log(product(r_i, i = 1 to
! floor(n/2)))
! Note that we take advantage of integer division on n/2
y = ONE
do i = 1, n/2
y = y * prn()
end do
r1 = prn()
r2 = prn()
c = cos(PI/2.*r2)
x = -2./n * (log(y) + log(r1)*c*c)
end select
! If sampling a chi-squared distribution for a resonance width and the
! average resonance width has been given, return the sampled resonance
! width.
if (present(G_avg)) then
G = x * G_avg
else
G = x
end if
end function chi_squared
end module physics

123
src/plot.F90
View file

@ -2,12 +2,11 @@ module plotter
use constants
use error, only: fatal_error
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, cell_contains
use geometry, only: find_cell
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord
use particle_header, only: deallocate_coord
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
@ -18,17 +17,16 @@ module plotter
contains
!===============================================================================
! RUN_PLOT
! RUN_PLOT controls the logic for making one or many plots
!===============================================================================
subroutine run_plot()
integer :: i
type(Plot), pointer :: pl => null()
integer :: i ! loop index for plots
type(Plot), pointer :: pl => null()
do i=1,n_plots
do i = 1, n_plots
pl => plots(i)
! Display output message
@ -36,31 +34,34 @@ contains
call write_message(5)
if (pl % type == PLOT_TYPE_SLICE) then
! create 2d image
call create_ppm(pl)
end if
end do
end subroutine run_plot
!===============================================================================
! create_ppm
! CREATE_PPM creates an image based on user input from a plots.xml <plot>
! specification in the portable pixmap format (PPM)
!===============================================================================
subroutine create_ppm(pl)
type(Plot), pointer :: pl
type(Plot), pointer :: pl
type(Image) :: img
integer :: in_i, out_i
integer :: x, y
integer :: r, g, b
real(8) :: in_pixel, out_pixel
real(8) :: xyz(3)
logical :: found_cell
type(Cell), pointer :: c => null()
integer :: in_i
integer :: out_i
integer :: x, y ! pixel location
integer :: r, g, b ! colors (red, green, blue) from 0-255
real(8) :: in_pixel
real(8) :: out_pixel
real(8) :: xyz(3)
logical :: found_cell
type(Image) :: img
type(Cell), pointer :: c => null()
! Initialize and allocate space for image
call init_image(img)
call allocate_image(img, pl % pixels(1), pl % pixels(2))
@ -70,7 +71,7 @@ contains
in_pixel = pl % width(1)/dble(pl % pixels(1))
out_pixel = pl % width(2)/dble(pl % pixels(2))
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
xyz(2) = pl % origin(2) - pl % width(2) / 2.0
xyz(2) = pl % origin(2) + pl % width(2) / 2.0
xyz(3) = pl % origin(3)
else if (pl % basis == PLOT_BASIS_XZ) then
in_i = 1
@ -79,7 +80,7 @@ contains
out_pixel = pl % width(2)/dble(pl % pixels(2))
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
xyz(2) = pl % origin(2)
xyz(3) = pl % origin(3) - pl % width(2) / 2.0
xyz(3) = pl % origin(3) + pl % width(2) / 2.0
else if (pl % basis == PLOT_BASIS_YZ) then
in_i = 2
out_i = 3
@ -87,7 +88,7 @@ contains
out_pixel = pl % width(2)/dble(pl % pixels(2))
xyz(1) = pl % origin(1)
xyz(2) = pl % origin(2) - pl % width(1) / 2.0
xyz(3) = pl % origin(3) - pl % width(2) / 2.0
xyz(3) = pl % origin(3) + pl % width(2) / 2.0
end if
! allocate and initialize particle
@ -97,8 +98,8 @@ contains
p % coord % uvw = (/ 1, 0, 0 /)
p % coord % universe = BASE_UNIVERSE
do y=1, img%height
do x=1, img%width
do y = 1, img % height
do x = 1, img % width
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
@ -106,65 +107,74 @@ contains
call find_cell(found_cell)
if (.not. found_cell) then
r = 255
g = 255
b = 255
! If no cell, revert to default color
r = pl % not_found % rgb(1)
g = pl % not_found % rgb(2)
b = pl % not_found % rgb(3)
else
c => cells(p % coord % cell)
if (pl % color == PLOT_COLOR_MATS) then
r = materials(c % material) % rgb(1)
g = materials(c % material) % rgb(2)
b = materials(c % material) % rgb(3)
else if (pl % color == PLOT_COLOR_CELLS) then
r = c % rgb(1)
g = c % rgb(2)
b = c % rgb(3)
else
r = 0
g = 0
b = 0
end if
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
c => cells(p % coord % cell)
r = pl % colors(c % material) % rgb(1)
g = pl % colors(c % material) % rgb(2)
b = pl % colors(c % material) % rgb(3)
else if (pl % color_by == PLOT_COLOR_CELLS) then
! Assign color based on cell
r = pl % colors(p % coord % cell) % rgb(1)
g = pl % colors(p % coord % cell) % rgb(2)
b = pl % colors(p % coord % cell) % rgb(3)
else
r = 0
g = 0
b = 0
end if
end if
! Create a pixel at (x,y) with color (r,g,b)
call set_pixel(img, x, y, r, g, b)
! Advance pixel in first direction
p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel
end do
! Advance pixel in second direction
p % coord0 % xyz(in_i) = xyz(in_i)
p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) + out_pixel
p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel
end do
! Write out the ppm to a file
call output_ppm(pl,img)
! Free up space
call deallocate_image(img)
end subroutine create_ppm
!===============================================================================
! output_ppm
! OUTPUT_PPM writes out a previously generated image to a PPM file
!===============================================================================
subroutine output_ppm(pl,img)
type(Plot), pointer :: pl
type(Image), intent(in) :: img
subroutine output_ppm(pl, img)
integer :: i, j
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
type(Plot), pointer :: pl
type(Image), intent(in) :: img
path_plot = trim(path_input) // "slice" // trim(to_str(pl % id)) // ".ppm"
open(UNIT=UNIT_PLOT, FILE=path_plot)
integer :: i ! loop index for height
integer :: j ! loop index for width
! Open PPM file for writing
open(UNIT=UNIT_PLOT, FILE=pl % path_plot)
! Write header
write(UNIT_PLOT, '(A2)') 'P6'
write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height
write(UNIT_PLOT, '(A)') '255'
do j=1, img%height
do i=1, img%width
! Write color for each pixel
do j = 1, img % height
do i = 1, img % width
write(UNIT_PLOT, '(3A1)', advance='no') achar(img%red(i,j)), &
achar(img%green(i,j)), &
achar(img%blue(i,j))
achar(img%green(i,j)), achar(img%blue(i,j))
end do
end do
@ -173,5 +183,4 @@ contains
end subroutine output_ppm
end module plotter

18
src/plot_header.F90
View file

@ -1,20 +1,32 @@
module plot_header
use constants
implicit none
!===============================================================================
! PLOT hold plot information
! ObjectColor holds color information for plotted objects
!===============================================================================
type ObjectColor
integer :: rgb(3)
end type ObjectColor
!===============================================================================
! PLOT holds plot information
!===============================================================================
type Plot
integer :: id ! Unique ID
character(MAX_LINE_LEN) :: path_plot ! path for plot file
integer :: type ! Type
integer :: color ! quantity to color regions by
integer :: color_by ! quantity to color regions by
real(8) :: origin(3) ! xyz center of plot location
real(8) :: aspect ! spacing between rays in raytracer
real(8) :: width(3) ! xyz widths of plot
integer :: basis ! direction of plot slice
integer :: pixels(2) ! pixel width/height of plot slice
type(ObjectColor) :: not_found ! color for positions where no cell found
type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats
end type Plot
integer :: PLOT_TYPE_SLICE = 1

66
src/ppmlib.F90
View file

@ -13,6 +13,10 @@ module ppmlib
contains
!===============================================================================
! INIT_IMAGE initializes the Image derived type
!===============================================================================
subroutine init_image(img)
type(Image), intent(out) :: img
@ -26,29 +30,46 @@ contains
end subroutine init_image
!===============================================================================
! ALLOCATE_IMAGE sets the width and height of an image and allocates color
! arrays
!===============================================================================
subroutine allocate_image(img, w, h)
type(Image), intent(inout) :: img
integer, intent(in) :: w, h
integer, intent(in) :: w ! width of image
integer, intent(in) :: h ! height of image
! allocate red, green, and blue array
allocate(img % red(w, h))
allocate(img % green(w, h))
allocate(img % blue(w, h))
! set width and height
img % width = w
img % height = h
end subroutine allocate_image
!===============================================================================
! DEALLOCATE_IMAGE
!===============================================================================
subroutine deallocate_image(img)
type(Image) :: img
if ( associated(img % red) ) deallocate(img % red)
if ( associated(img % green) ) deallocate(img % green)
if ( associated(img % blue) ) deallocate(img % blue)
if (associated(img % red)) deallocate(img % red)
if (associated(img % green)) deallocate(img % green)
if (associated(img % blue)) deallocate(img % blue)
end subroutine deallocate_image
!===============================================================================
! INSIDE_IMAGE determines whether a point (x,y) is inside the image
!===============================================================================
function inside_image(img, x, y) result(inside)
@ -58,13 +79,15 @@ contains
inside = .false.
if ( ( x < img % width) .and. &
( y < img % height ) .and. &
( x >= 0 ) .and. &
( y >= 0 ) ) inside = .true.
if ((x < img % width) .and. (y < img % height) .and. &
(x >= 0) .and. (y >= 0)) inside = .true.
end function inside_image
!===============================================================================
! VALID_IMAGE checks whether the image has a width and height and if its color
! arrays are allocated
!===============================================================================
function valid_image(img) result(valid)
@ -73,29 +96,32 @@ contains
valid = .false.
if ( img % width == 0 ) return
if ( img % height == 0 ) return
if ( .not. associated(img % red) .or. &
if (img % width == 0) return
if (img % height == 0) return
if (.not. associated(img % red) .or. &
.not. associated(img % green) .or. &
.not. associated(img % blue) ) return
.not. associated(img % blue)) return
valid = .true.
end function valid_image
!===============================================================================
! SET_PIXEL sets the
!===============================================================================
subroutine set_pixel(img, x, y, r, g, b)
type(Image), intent(inout) :: img
integer, intent(in) :: x, y
integer, intent(in) :: r, g, b
integer, intent(in) :: x, y ! coordinates
integer, intent(in) :: r, g, b ! red, green, and blue
if ( inside_image(img, x, y) .and. valid_image(img)) then
img % red(x+1,y+1) = mod(abs(r), 256)
img % green(x+1, y+1) = mod(abs(g), 256)
img % blue(x+1, y+1) = mod(abs(b), 256)
if (inside_image(img, x, y) .and. valid_image(img)) then
img % red(x+1,y+1) = mod(abs(r), 256)
img % green(x+1, y+1) = mod(abs(g), 256)
img % blue(x+1, y+1) = mod(abs(b), 256)
end if
end subroutine set_pixel
end subroutine set_pixel
end module ppmlib

1
src/search.F90
View file

@ -1,6 +1,5 @@
module search
use constants, only: ONE, MAX_LINE_LEN
use error, only: fatal_error
use global, only: message

158
src/source.F90
View file

@ -1,7 +1,7 @@
module source
use bank_header, only: Bank
use constants, only: ONE
use constants
use error, only: fatal_error
use geometry_header, only: BASE_UNIVERSE
use global
@ -11,6 +11,10 @@ module source
use random_lcg, only: prn, set_particle_seed
use string, only: to_str
#ifdef MPI
use mpi
#endif
implicit none
contains
@ -43,6 +47,13 @@ contains
! Determine maximum amount of particles to simulate on each processor
maxwork = ceiling(real(n_particles)/n_procs,8)
! ID's of first and last source particles
bank_first = rank*maxwork + 1
bank_last = min((rank+1)*maxwork, n_particles)
! number of particles for this processor
work = bank_last - bank_first + 1
! Allocate source bank
allocate(source_bank(maxwork), STAT=alloc_err)
if (alloc_err /= 0) then
@ -69,16 +80,16 @@ contains
call fatal_error()
end if
! Read the source from a binary file instead of sampling from some assumed
! source distribution
if (external_source % type == SRC_FILE) then
call read_source_binary()
return
end if
! Initialize first cycle source bank
do i = 0, n_procs - 1
if (rank == i) then
! ID's of first and last source particles
bank_first = rank*maxwork + 1
bank_last = min((rank+1)*maxwork, n_particles)
! number of particles for this processor
work = bank_last - bank_first + 1
do j = 1, work
id = bank_first + j - 1
source_bank(j) % id = id
@ -93,6 +104,8 @@ contains
p_max = external_source % values(4:6)
r = (/ (prn(), k = 1,3) /)
source_bank(j) % xyz = p_min + r*(p_max - p_min)
case (SRC_POINT)
source_bank(j) % xyz = external_source % values
end select
! sample angle
@ -145,12 +158,13 @@ contains
p % id = bank_first + index_source - 1
! set random number seed
particle_seed = (current_cycle - 1)*n_particles + p % id
particle_seed = ((current_batch - 1)*gen_per_batch + &
current_gen - 1)*n_particles + p % id
call set_particle_seed(particle_seed)
! set particle trace
trace = .false.
if (current_cycle == trace_cycle .and. &
if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
p % id == trace_particle) trace = .true.
end subroutine get_source_particle
@ -186,4 +200,128 @@ contains
end subroutine initialize_particle
!===============================================================================
! WRITE_SOURCE writes out the final source distribution to a binary file that
! can be used as a starting source in a new simulation
!===============================================================================
subroutine write_source_binary()
#ifdef MPI
integer :: fh ! file handle
integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file)
! ==========================================================================
! PARALLEL I/O USING MPI-2 ROUTINES
! Open binary source file for reading
call MPI_FILE_OPEN(MPI_COMM_WORLD, 'source.binary', MPI_MODE_CREATE + &
MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpi_err)
if (master) then
offset = 0
call MPI_FILE_WRITE_AT(fh, offset, n_particles, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpi_err)
end if
! Set proper offset for source data on this processor
offset = 8*(1 + rank*maxwork*8)
! Write all source sites
call MPI_FILE_WRITE_AT(fh, offset, source_bank(1), work, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Close binary source file
call MPI_FILE_CLOSE(fh, mpi_err)
#else
! ==========================================================================
! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES
! Open binary source file for writing
open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='replace', &
ACCESS='stream')
! Write the number of particles
write(UNIT=UNIT_SOURCE) n_particles
! Write information from the source bank
write(UNIT=UNIT_SOURCE) source_bank(1:work)
! Close binary source file
close(UNIT=UNIT_SOURCE)
#endif
end subroutine write_source_binary
!===============================================================================
! READ_SOURCE reads a source distribution from a source.binary file and
! initializes the source bank
!===============================================================================
subroutine read_source_binary()
integer(8) :: n_sites ! number of sites in binary file
#ifdef MPI
integer :: fh ! file handle
integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file)
#endif
#ifdef MPI
! ==========================================================================
! PARALLEL I/O USING MPI-2 ROUTINES
! Open binary source file for reading
call MPI_FILE_OPEN(MPI_COMM_WORLD, 'source.binary', MPI_MODE_RDONLY, &
MPI_INFO_NULL, fh, mpi_err)
! Read number of source sites in file
offset = 0
call MPI_FILE_READ_AT(fh, offset, n_sites, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpi_err)
! Check that number of source sites matches
if (n_sites /= n_particles) then
message = "No support yet for source files of different size than &
&specified number of particles per generation."
call fatal_error()
end if
! Set proper offset for source data on this processor
offset = 8*(1 + rank*maxwork*8)
! Read all source sites
call MPI_FILE_READ_AT(fh, offset, source_bank(1), work, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Close binary source file
call MPI_FILE_CLOSE(fh, mpi_err)
#else
! ==========================================================================
! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES
! Open binary source file for reading
open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='old', &
ACCESS='stream')
! Read number of source sites in file
read(UNIT=UNIT_SOURCE) n_sites
! Check that number of source sites matches
if (n_sites /= n_particles) then
message = "No support yet for source files of different size than &
&specified number of particles per generation."
call fatal_error()
end if
! Read position, angle, and energy
read(UNIT=UNIT_SOURCE) source_bank(1:work)
! Close binary source file
close(UNIT=UNIT_SOURCE)
#endif
end subroutine read_source_binary
end module source

2
src/source_header.F90
View file

@ -8,7 +8,7 @@ module source_header
!===============================================================================
type ExtSource
integer :: type ! type of source, e.g. 'box' or 'cell'
integer :: type ! type of source, e.g. 'box' or 'point'
real(8), allocatable :: values(:) ! values for particular source type
end type ExtSource

304
src/tally.F90
View file

@ -3,7 +3,8 @@ module tally
use constants
use error, only: fatal_error
use global
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, &
mesh_indices_to_bin, mesh_intersects
use mesh_header, only: StructuredMesh
use output, only: header
use search, only: binary_search
@ -345,6 +346,20 @@ contains
score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu)
case (SCORE_DIFFUSION)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle
! Temporarily store the scattering cross section
score = material_xs % total - material_xs % absorption
! Since this only gets tallied at every scattering event, the flux
! estimator is 1/Sigma_s. Therefore, the diffusion coefficient
! times flux is 1/(3*Sigma_s*(Sigma_t - mu*Sigma_s)).
score = last_wgt / (3.0_8 * score * (material_xs % total - &
mu * score))
case (SCORE_N_1N)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle
@ -533,6 +548,14 @@ contains
do i = 1, n_tracklength_tallies
t => tallies(tracklength_tallies(i))
! Check if this tally has a mesh filter -- if so, we treat it separately
! since multiple bins can be scored to with a single track
if (t % n_filter_bins(FILTER_MESH) > 0) then
call score_tl_on_mesh(tracklength_tallies(i), distance)
cycle
end if
! =======================================================================
! DETERMINE SCORING BIN COMBINATION
@ -592,6 +615,239 @@ contains
end subroutine score_tracklength_tally
!===============================================================================
! SCORE_TL_ON_MESH calculate fluxes and reaction rates based on the track-length
! estimate of the flux specifically for tallies that have mesh filters. For
! these tallies, it is possible to score to multiple mesh cells for each track.
!===============================================================================
subroutine score_tl_on_mesh(index_tally, d_track)
integer, intent(in) :: index_tally
real(8), intent(in) :: d_track
integer :: i ! loop index for filter/score bins
integer :: j ! loop index for direction
integer :: k ! loop index for mesh cell crossings
integer :: ijk0(3) ! indices of starting coordinates
integer :: ijk1(3) ! indices of ending coordinates
integer :: ijk_cross(3) ! indices of mesh cell crossed
integer :: n_cross ! number of surface crossings
integer :: bins(N_FILTER_TYPES) ! scoring bin combination
integer :: score_index ! single index for single bin
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
real(8) :: flux ! tracklength estimate of flux
real(8) :: score ! actual score (e.g., flux*xs)
real(8) :: uvw(3) ! cosine of angle of particle
real(8) :: xyz0(3) ! starting/intermediate coordinates
real(8) :: xyz1(3) ! ending coordinates of particle
real(8) :: xyz_cross(3) ! coordinates of next boundary
real(8) :: d(3) ! distance to each bounding surface
real(8) :: distance ! distance traveled in mesh cell
logical :: found_bin ! was a scoring bin found?
logical :: start_in_mesh ! starting coordinates inside mesh?
logical :: end_in_mesh ! ending coordinates inside mesh?
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
found_bin = .true.
bins = 1
t => tallies(index_tally)
! ==========================================================================
! CHECK IF THIS TRACK INTERSECTS THE MESH
! Copy starting and ending location of particle
xyz0 = p % coord0 % xyz - (d_track - TINY_BIT) * p % coord0 % uvw
xyz1 = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw
! Determine indices for starting and ending location
m => meshes(t % mesh)
call get_mesh_indices(m, xyz0, ijk0, start_in_mesh)
call get_mesh_indices(m, xyz1, ijk1, end_in_mesh)
! Check if start or end is in mesh -- if not, check if track still
! intersects with mesh
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
if (.not. mesh_intersects(m, xyz0, xyz1)) return
end if
! Reset starting and ending location
xyz0 = p % coord0 % xyz - d_track * p % coord0 % uvw
xyz1 = p % coord0 % xyz
! =========================================================================
! CHECK FOR SCORING COMBINATION FOR FILTERS OTHER THAN MESH
FILTER_LOOP: do i = 1, t % n_filters
select case (t % filters(i))
case (FILTER_UNIVERSE)
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
bins(FILTER_UNIVERSE) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, index_tally)
if (bins(FILTER_UNIVERSE) == NO_BIN_FOUND) then
found_bin = .false.
return
end if
case (FILTER_MATERIAL)
bins(FILTER_MATERIAL) = get_next_bin(FILTER_MATERIAL, &
p % material, index_tally)
if (bins(FILTER_MATERIAL) == NO_BIN_FOUND) then
found_bin = .false.
return
end if
case (FILTER_CELL)
! determine next cell bin
! TODO: Account for cells in multiple levels when performing this filter
bins(FILTER_CELL) = get_next_bin(FILTER_CELL, &
p % coord % cell, index_tally)
if (bins(FILTER_CELL) == NO_BIN_FOUND) then
found_bin = .false.
return
end if
case (FILTER_CELLBORN)
! determine next cellborn bin
bins(FILTER_CELLBORN) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, index_tally)
if (bins(FILTER_CELLBORN) == NO_BIN_FOUND) then
found_bin = .false.
return
end if
case (FILTER_SURFACE)
! determine next surface bin
bins(FILTER_SURFACE) = get_next_bin(FILTER_SURFACE, &
p % surface, index_tally)
if (bins(FILTER_SURFACE) == NO_BIN_FOUND) then
found_bin = .false.
return
end if
case (FILTER_ENERGYIN)
! determine incoming energy bin
k = t % n_filter_bins(FILTER_ENERGYIN)
! check if energy of the particle is within energy bins
if (p % E < t % energy_in(1) .or. &
p % E > t % energy_in(k + 1)) then
found_bin = .false.
return
end if
! search to find incoming energy bin
bins(FILTER_ENERGYIN) = binary_search(t % energy_in, k + 1, p % E)
end select
end do FILTER_LOOP
! ==========================================================================
! DETERMINE WHICH MESH CELLS TO SCORE TO
! Calculate number of surface crossings
n_cross = sum(abs(ijk1 - ijk0)) + 1
! Copy particle's direction
uvw = p % coord0 % uvw
! Bounding coordinates
do j = 1, 3
if (uvw(j) > 0) then
xyz_cross(j) = m % lower_left(j) + ijk0(j) * m % width(j)
else
xyz_cross(j) = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
end if
end do
MESH_LOOP: do k = 1, n_cross
found_bin = .false.
! Calculate distance to each bounding surface. We need to treat special
! case where the cosine of the angle is zero since this would result in a
! divide-by-zero.
if (k == n_cross) xyz_cross = xyz1
do j = 1, 3
if (uvw(j) == 0) then
d(j) = INFINITY
else
d(j) = (xyz_cross(j) - xyz0(j))/uvw(j)
end if
end do
! Determine the closest bounding surface of the mesh cell by calculating
! the minimum distance
j = minloc(d, 1)
distance = d(j)
! Now use the minimum distance and diretion of the particle to determine
! which surface was crossed
if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then
ijk_cross = ijk0
found_bin = .true.
end if
! Increment indices and determine new crossing point
if (uvw(j) > 0) then
ijk0(j) = ijk0(j) + 1
xyz_cross(j) = xyz_cross(j) + m % width(j)
else
ijk0(j) = ijk0(j) - 1
xyz_cross(j) = xyz_cross(j) - m % width(j)
end if
! =======================================================================
! SCORE TO THIS MESH CELL
if (found_bin) then
! Calculate track-length estimate of flux
flux = p % wgt * distance
! Determine mesh bin
bins(FILTER_MESH) = mesh_indices_to_bin(m, ijk_cross)
! Determining scoring index
score_index = sum((bins - 1) * t % stride) + 1
! Determine score for each bin
SCORE_LOOP: do i = 1, t % n_score_bins
! determine what type of score bin
score_bin = t % score_bins(i) % scalar
! Determine cross section
select case(score_bin)
case (SCORE_FLUX)
score = flux
case (SCORE_TOTAL)
score = material_xs % total * flux
case (SCORE_SCATTER)
score = (material_xs % total - material_xs % absorption) * flux
case (SCORE_ABSORPTION)
score = material_xs % absorption * flux
case (SCORE_FISSION)
score = material_xs % fission * flux
case (SCORE_NU_FISSION)
score = material_xs % nu_fission * flux
end select
! Add score to tally
call add_to_score(t % scores(score_index, i), score)
end do SCORE_LOOP
end if
! Calculate new coordinates
xyz0 = xyz0 + distance * uvw
end do MESH_LOOP
end subroutine score_tl_on_mesh
!===============================================================================
! GET_SCORING_BINS determines a combination of filter bins that should be scored
! for a tally based on the particle's current attributes.
@ -606,6 +862,7 @@ contains
integer :: i ! loop index for filters
integer :: n ! number of bins for single filter
integer :: mesh_bin ! index for mesh bin
real(8) :: E ! particle energy
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
@ -677,15 +934,22 @@ contains
case (FILTER_ENERGYIN)
! determine incoming energy bin
n = t % n_filter_bins(FILTER_ENERGYIN)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
else
E = p % last_E
end if
! check if energy of the particle is within energy bins
if (p % last_E < t % energy_in(1) .or. &
p % last_E > t % energy_in(n + 1)) then
if (E < t % energy_in(1) .or. E > t % energy_in(n + 1)) then
found_bin = .false.
return
end if
! search to find incoming energy bin
bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E)
bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, E)
case (FILTER_ENERGYOUT)
! determine outgoing energy bin
@ -1048,28 +1312,32 @@ contains
end subroutine add_to_score
!===============================================================================
! ACCUMULATE_CYCLE_ESTIMATE
! ACCUMULATE_BATCH_ESTIMATE accumulates scores from many histories (or many
! generations) into a single realization of a random variable.
!===============================================================================
elemental subroutine accumulate_cycle_estimate(score)
elemental subroutine accumulate_batch_estimate(score)
type(TallyScore), intent(inout) :: score
real(8) :: val
! Add the sum and square of the sum of contributions from each cycle
! within a cycle to the variables sum and sum_sq. This will later allow us
! to calculate a variance on the tallies
score % sum = score % sum + score % value/n_particles
score % sum_sq = score % sum_sq + (score % value/n_particles)**2
val = score % value/(n_particles*gen_per_batch)
score % sum = score % sum + val
score % sum_sq = score % sum_sq + val*val
! Reset the single cycle estimate
! Reset the single batch estimate
score % value = ZERO
end subroutine accumulate_cycle_estimate
end subroutine accumulate_batch_estimate
!===============================================================================
! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history
! within the cycle to a new random variable
! within the batch to a new random variable
!===============================================================================
subroutine synchronize_tallies()
@ -1086,7 +1354,7 @@ contains
t => tallies(i)
! Loop over all filter and scoring bins
call accumulate_cycle_estimate(t % scores)
call accumulate_batch_estimate(t % scores)
end do
end subroutine synchronize_tallies
@ -1183,6 +1451,7 @@ contains
score_name(abs(SCORE_SCATTER_1)) = "First Scattering Moment"
score_name(abs(SCORE_SCATTER_2)) = "Second Scattering Moment"
score_name(abs(SCORE_SCATTER_3)) = "Third Scattering Moment"
score_name(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient"
score_name(abs(SCORE_N_1N)) = "(n,1n) Rate"
score_name(abs(SCORE_N_2N)) = "(n,2n) Rate"
score_name(abs(SCORE_N_3N)) = "(n,3n) Rate"
@ -1522,16 +1791,21 @@ contains
end function get_label
!===============================================================================
! CALCULATE_STATISTICS
! CALCULATE_STATISTICS determines the sample mean and the standard deviation of
! the mean for a TallyScore.
!===============================================================================
elemental subroutine calculate_statistics(score)
type(TallyScore), intent(inout) :: score
! Calculate mean and standard deviation of the mean
! Calculate sample mean and standard deviation of the mean -- note that we
! have used Bessel's correction so that the estimator of the variance of the
! sample mean is unbiased.
score % sum = score % sum/n_active
score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2)/n_active)
score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / &
(n_active - 1))
end subroutine calculate_statistics

10
src/timing.F90
View file

@ -13,7 +13,7 @@ module timing
type Timer
logical :: running = .false. ! is timer running?
integer :: start_counts = 0 ! counts when started
real(8) :: elapsed = 0. ! total time elapsed in seconds
real(8) :: elapsed = ZERO ! total time elapsed in seconds
!!$ contains
!!$ procedure :: start => timer_start
!!$ procedure :: get_value => timer_get_value
@ -24,7 +24,7 @@ module timing
contains
!===============================================================================
! TIMER_START
! TIMER_START starts running a timer and measures the current time
!===============================================================================
subroutine timer_start(self)
@ -38,7 +38,7 @@ contains
end subroutine timer_start
!===============================================================================
! TIMER_GET_VALUE
! TIMER_GET_VALUE returns the current value of the timer
!===============================================================================
function timer_get_value(self) result(elapsed)
@ -61,7 +61,7 @@ contains
end function timer_get_value
!===============================================================================
! TIMER_STOP
! TIMER_STOP stops the timer and sets the elapsed time
!===============================================================================
subroutine timer_stop(self)
@ -78,7 +78,7 @@ contains
end subroutine timer_stop
!===============================================================================
! TIMER_RESET
! TIMER_RESET resets a timer to have a zero value
!===============================================================================
subroutine timer_reset(self)

33
src/utils/read_source.py Executable file
View file

@ -0,0 +1,33 @@
#!/usr/bin/env python
import struct
class SourceFile(object):
def __init__(self, filename):
# Open source file for reading
self.f = open(filename, 'r')
# Read number of source sites
self.n_sites = struct.unpack('q', self.f.read(8))[0]
# Create list to store source sites
self.sites = []
for i in range(self.n_sites):
# Read position, angle, and energy
uid = struct.unpack('q', self.f.read(8))
xyz = struct.unpack('3d', self.f.read(24))
uvw = struct.unpack('3d', self.f.read(24))
E = struct.unpack('d', self.f.read(8))
# Create source site and append to list
self.sites.append(SourceSite(xyz, uvw, E))
class SourceSite(object):
def __init__(self, xyz, uvw, E):
self.xyz = xyz
self.uvw = uvw
self.E = E
def __repr__(self):
return "<SourceSite: xyz={0} at E={1}>".format(self.xyz, self.E)

15
src/xml-fortran/templates/plots_t.xml
View file

@ -3,15 +3,28 @@
<options rootname="plots" />
<typedef name="col_spec_xml">
<component name="id" type="integer" />
<component name="rgb" type="integer-array" />
</typedef>
<typedef name="mask_xml">
<component name="components" type="integer-array" />
<component name="background" type="integer-array" />
</typedef>
<typedef name="plot_xml">
<component name="id" type="integer" />
<component name="filename" type="word" length="50" default="'plot'" />
<component name="type" type="word" length="10" default="'slice'"/>
<component name="color" type="word" length="10" default="'cell'"/>
<component name="origin" type="double-array" />
<component name="width" type="double-array" />
<component name="basis" type="word" length="3" default="'xy'" />
<component name="pixels" type="integer-array" />
<component name="aspect" type="double" default="0.1" />
<component name="background" type="integer-array"/>
<component name="col_spec_" tag="col_spec" type="col_spec_xml" dimension="1" />
<component name="mask_" tag="mask" type="mask_xml" dimension="1" />
</typedef>
<variable name="plot_" tag="plot" type="plot_xml" dimension="1" />

16
src/xml-fortran/templates/settings_t.xml
View file

@ -4,10 +4,11 @@
<options rootname="settings" />
<typedef name="criticality_xml">
<component name="cycles" type="integer" />
<component name="batches" type="integer" default="0" />
<component name="inactive" type="integer" />
<component name="particles" type="word" length="25" default="''" />
<component name="seed" type="integer" />
<component name="generations_per_batch" type="integer" default="1" />
</typedef>
<typedef name="source_xml">
@ -27,13 +28,16 @@
</typedef>
<variable name="criticality" type="criticality_xml" />
<variable name="verbosity_" tag="verbosity" type="integer" />
<variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
<variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
<variable name="energy_grid_" tag="energy_grid" type="word" length="7" />
<variable name="entropy_" tag="entropy" type="entropy_xml" dimension="1" />
<variable name="ptables_" tag="ptables" type="word" length="3" />
<variable name="seed_" tag="seed" type="integer" />
<variable name="source_" tag="source" type="source_xml" />
<variable name="survival_" tag="survival_biasing" type="word" length="3" />
<variable name="ptables_" tag="ptables" type="word" length="3" />
<variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
<variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
<variable name="verbosity_" tag="verbosity" type="integer" />
<variable name="trace_" tag="trace" type="integer-array" />
<variable name="entropy_" tag="entropy" type="entropy_xml" dimension="1" />
<variable name="write_source_" tag="write_source" type="word" length="3" />
</template>