crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merged master branch to cmfd

Browse source
This commit is contained in:
Bryan Herman 2011-12-10 14:09:54 -05:00
commit bf9cbb0c16
63 changed files with 5270 additions and 2924 deletions
Download patch file Download diff file Expand all files Collapse all files

2076
cross_sections.xml Normal file
View file

File diff suppressed because it is too large Load diff

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

@ -64,14 +64,14 @@ number of cells, surfaces, and lattices. Let us look at the following example::
<!-- This is a comment -->
<surface>
<uid>1</uid>
<id>1</id>
<type>sphere</type>
<coeffs>0.0 0.0 0.0 5.0</coeffs>
<boundary>vacuum</boundary>
<surface>
<cell>
<uid>1</uid>
<id>1</id>
<universe>0</universe>
<material>1</material>
<surfaces>-1</surfaces>
@ -89,14 +89,17 @@ could be written as::
<geometry>
<!-- This is a comment -->
<surface uid="1" type="sphere" coeffs="0.0 0.0 0.0 5.0" boundary="vacuum" />
<cell uid="1" universe="0" material="1" surfaces="-1" />
<surface id="1" type="sphere" coeffs="0.0 0.0 0.0 5.0" boundary="vacuum" />
<cell id="1" universe="0" material="1" surfaces="-1" />
</geometry>
``surface`` Element
-------------------
Each ``surface`` element can have the following attributes or sub-elements:
:uid:
:id:
A unique integer that can be used to identify the surface.
*Default*: None
@ -118,75 +121,116 @@ Each ``surface`` element can have the following attributes or sub-elements:
*Default*: ``reflective``
The following quadratic surfaces can be modeled:
:x-plane:
A plane perpendicular to the x axis, i.e. a surface of the form :math:`x -
x_0 = 0`. The coefficients specified are ":math:`x_0`".
:y-plane:
A plane perpendicular to the y axis, i.e. a surface of the form :math:`y -
y_0 = 0`. The coefficients specified are ":math:`y_0`".
:z-plane:
A plane perpendicular to the z axis, i.e. a surface of the form :math:`z -
z_0 = 0`. The coefficients specified are ":math:`z_0`".
:plane:
An arbitrary plane of the form :math:`Ax + By + Cz = D`. The coefficients
specified are ":math:`A \: B \: C \: D`".
:x-cylinder:
An infinite cylinder whose length is paralle to the x-axis. This is a
quadratic surface of the form :math:`(y - y_0)^2 + (z - z_0)^2 = R^2`. The
coefficients specified are ":math:`y_0 \: z_0 \: R`".
:y-cylinder:
An infinite cylinder whose length is paralle to the y-axis. This is a
quadratic surface of the form :math:`(x - x_0)^2 + (z - z_0)^2 = R^2`. The
coefficients specified are ":math:`x_0 \: z_0 \: R`".
:z-cylinder:
An infinite cylinder whose length is paralle to the z-axis. This is a
quadratic surface of the form :math:`(x - x_0)^2 + (y - y_0)^2 = R^2`. The
coefficients specified are ":math:`x_0 \: y_0 \: R`".
:sphere:
A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 =
R^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0 \: R`".
``cell`` Element
----------------
Each ``cell`` element can have the following attributes or sub-elements:
:uid:
:id:
A unique integer that can be used to identify the surface.
*Default*: None
:universe:
The ``uid`` of the universe that this cell is contained in.
The ``id`` of the universe that this cell is contained in.
*Default*: 0
:fill:
The ``uid`` of the universe that fills this cell.
The ``id`` of the universe that fills this cell.
.. note:: If a fill is specified, no material should be given.
*Default*: None
:material:
The ``uid`` of the material that this cell contains.
The ``id`` of the material that this cell contains.
.. note:: If a material is specified, no fill should be given.
*Default*: None
:surfaces:
A list of the ``uids`` for surfaces that bound this cell, e.g. if the cell
A list of the ``ids`` for surfaces that bound this cell, e.g. if the cell
is on the negative side of surface 3 and the positive side of surface 5, the
bounding surfaces would be given as "-3 5".
*Default*: None
The following quadratic surfaces can be modeled:
``lattice`` Element
-------------------
:x-plane:
A plane perpendicular to the x axis, i.e. a surface of the form :math:`x - x_0
= 0`. The coefficients specified are ":math:`x_0`".
The ``lattice`` can be used to represent repeating structures (e.g. fuel pins in
an assembly) or other geometry which naturally fits into a two-dimensional
structured mesh. Each cell within the lattice is filled with a specified
universe. A ``lattice`` accepts the following attributes or sub-elements:
:y-plane:
A plane perpendicular to the y axis, i.e. a surface of the form :math:`y - y_0
= 0`. The coefficients specified are ":math:`y_0`".
:id:
A unique integer that can be used to identify the surface.
:z-plane:
A plane perpendicular to the z axis, i.e. a surface of the form :math:`z - z_0
= 0`. The coefficients specified are ":math:`z_0`".
:type:
A string indicating the arrangement of lattice cells. Accepted options are
"rectangular" and "hexagonal".
:plane:
An arbitrary plane of the form :math:`Ax + By + Cz = D`. The coefficients
specified are ":math:`A \: B \: C \: D`".
*Default*: rectangular
:x-cylinder:
An infinite cylinder whose length is paralle to the x-axis. This is a
quadratic surface of the form :math:`(y - y_0)^2 + (z - z_0)^2 = R^2`. The
coefficients specified are ":math:`y_0 \: z_0 \: R`".
:dimension:
Two integers representing the number of lattice cells in the x- and y-
directions, respectively.
:y-cylinder:
An infinite cylinder whose length is paralle to the y-axis. This is a
quadratic surface of the form :math:`(x - x_0)^2 + (z - z_0)^2 = R^2`. The
coefficients specified are ":math:`x_0 \: z_0 \: R`".
*Default*: None
:z-cylinder:
An infinite cylinder whose length is paralle to the z-axis. This is a
quadratic surface of the form :math:`(x - x_0)^2 + (y - y_0)^2 = R^2`. The
coefficients specified are ":math:`x_0 \: y_0 \: R`".
:origin:
The coordinates of the lower-left corner of the lattice.
:sphere:
A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 =
R^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0 \: R`".
*Default*: None
:width:
The width of the lattice cell in the x- and y- directions.
*Default*: None
:universes:
A list of the universe numbers that fill each cell of the lattice.
*Default*: None
.. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
@ -196,8 +240,14 @@ The following quadratic surfaces can be modeled:
Materials Specification -- materials.xml
----------------------------------------
``material`` Element
--------------------
Each ``material`` element can have the following attributes or sub-elements:
:id:
A unique integer that can be used to identify the material.
:density:
An element with attributes/sub-elements called ``value`` and ``units``. The
``value`` attribute is the numeric value of the density while the ``units``
@ -209,7 +259,6 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: None
:nuclide:
An element with attributes/sub-elements called ``name``, ``xs``, and ``ao``
or ``wo``. The ``name`` attribute is the name of the cross-section for a
desired nuclide while the ``xs`` attribute is the cross-section
@ -225,25 +274,35 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: None
:sab:
Associates an S(a,b) table with the material. This element has
attributes/sub-elements called ``name`` and ``xs``. The ``name`` attribute
is the name of the S(a,b) table that should be associated with the material,
and ``xs`` is the cross-section identifier for the table.
*Default*: None
``default_xs`` Element
----------------------
In some circumstances, the cross-section identifier may be the same for many or
all nuclides in a given problem. In this case, rather than specifying the
``xs=...`` attribute on every nuclide, a ``default_xs`` element can be used to
set the default cross-section identifier for any nuclide without an identifier
explicitly listed. This element has no attributes and accepts a 3-letter string
that indicates the default cross-section identifier, e.g. "70c".
*Default*: None
--------------------------------------
Settings Specification -- settings.xml
--------------------------------------
All simulation parameters and miscellaneous options are specified in the
settings.xml file. The following elements can be specified:
settings.xml file.
- ``xslibrary``
- ``criticality``
- ``verbosity``
- ``source``
The ``xslibrary`` element has the following attributes:
:path:
The absolute or relative path of the xsdata file which lists cross sections
to be used in the simulation.
*Default*: None
``criticality`` Element
-----------------------
The ``criticality`` element indicates that a criticality calculation should be
performed. It has the following attributes/sub-elements:
@ -266,14 +325,25 @@ performed. It has the following attributes/sub-elements:
*Default*: None
The ``verbosity`` element tells the code how much information to display to the
standard output. A higher verbosity corresponds to more information being
displayed. This element takes the following attributes:
``cross_sections`` Element
--------------------------
:value:
The specified verbosity between 1 and 10.
The ``cross_sections`` element has no attributes and simply indicates the path
to an XML cross section listing file (usually named cross_sections.xml). If this
element is absent from the settings.xml file, the environment variable
``CROSS_SECTIONS`` will be used to find the path to the XML cross section
listing.
*Default*: 5
``cutoff`` Element
------------------
The ``cutoff`` element has no attributes and indicates the weight cutoff used
below which particles undergo Russian roulette.
*Default*: 0.25
``source`` Element
------------------
The ``source`` element gives information on an initial source guess for
criticality calculations. It takes the following attributes:
@ -288,6 +358,27 @@ criticality calculations. It takes the following attributes:
and the last three of which specify the upper-right corner. Source sites are
sampled uniformly through that parallelepiped.
``survival_biasing`` Element
----------------------------
The ``survival_biasing`` element has no attributes and assumes wither the
value ``on`` or ``off``. If turned on, this option will enable the use of
survival biasing, otherwise known as implicit capture or absorption.
*Default*: off
``verbosity`` Element
---------------------
The ``verbosity`` element tells the code how much information to display to the
standard output. A higher verbosity corresponds to more information being
displayed. This element takes the following attributes:
:value:
The specified verbosity between 1 and 10.
*Default*: 5
------------------------------------
Tallies Specification -- tallies.xml
------------------------------------
@ -306,8 +397,12 @@ filters can be used for a tally. The following types of filter are available:
cell, universe, material, surface, birth region, pre-collision energy,
post-collision energy, and an arbitrary structured mesh.
The two valid elements in the tallies.xml file are ``tally`` and ``mesh``. The
``tally`` element accepts the following sub-elements:
The two valid elements in the tallies.xml file are ``tally`` and ``mesh``.
``tally`` Element
-----------------
The ``tally`` element accepts the following sub-elements:
:filters:
A list of filters to specify what region of phase space should contribute to
@ -387,6 +482,9 @@ The following responses can be tallied.
:nu-fission:
Total production of neutrons due to fission
``mesh`` Element
----------------
If a structured mesh is desired as a filter for a tally, it must be specified in
a separate element with the tag name ``mesh``. This element has the following
attributes/sub-elements:
@ -404,3 +502,47 @@ attributes/sub-elements:
:width:
The width of mesh cells in each direction.
-------------------------------------------
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.
``origin`` Element
------------------
The ``origin`` element has no attributes/sub-elements and indicates the
Cartesian coordinates of the center of the plot.
*Default*: None
``width`` Element
-----------------
The ``width`` element has no attributes/sub-elements and indicates the width of
the plot in each of the basis directions.
*Default*: None
``basis`` Element
-----------------
The ``basis`` element has no attributes/sub-elements and indicates the specified
basis for plotting.
.. note:: The only accepted option currently is "xy"
*Default*: xy
``pixel`` Element
-----------------
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.
*Default*: 0.01

58
examples/lattice/geometry.xml
View file

@ -1,58 +0,0 @@
<?xml version="1.0"?>
<geometry>
<cell id="1">
<universe>0</universe>
<fill>111</fill>
<surfaces>1 -2 3 -4</surfaces>
</cell>
<cell id="100">
<universe>3</universe>
<material>40</material>
<surfaces>-5</surfaces>
</cell>
<cell id="101">
<universe>3</universe>
<material>41</material>
<surfaces>5</surfaces>
</cell>
<lattice id="111">
<type>rectangular</type>
<dimension>20 20</dimension>
<origin>-10.0 -10.0</origin>
<width>1.0 1.0</width>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
</universes>
</lattice>
<surface id="1" type="x-plane" coeffs="-9.9" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="9.9" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-9.9" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="9.9" boundary="vacuum" />
<surface id="5">
<type>z-cylinder</type>
<coeffs>0.0 0.0 0.4</coeffs>
</surface>
</geometry>

17
examples/lattice/materials.xml
View file

@ -1,17 +0,0 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="03c" 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" />
</material>
</materials>

45
examples/lattice/nested/geometry.xml Normal file
View file

@ -0,0 +1,45 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="6" surfaces="1 -2 3 -4" />
<cell id="2" universe="5" fill="4" surfaces="1 -2 3 -4" />
<cell id="101" universe="1" material="1" surfaces="-5" />
<cell id="102" universe="1" material="2" surfaces="5" />
<cell id="201" universe="2" material="1" surfaces="-6" />
<cell id="202" universe="2" material="2" surfaces="6" />
<cell id="301" universe="3" material="1" surfaces="-7" />
<cell id="302" universe="3" material="2" surfaces="7" />
<!-- 4 x 4 assembly -->
<lattice id="4">
<type>rectangular</type>
<dimension>2 2</dimension>
<origin>-1.0 -1.0</origin>
<width>1.0 1.0</width>
<universes>
1 2
2 3
</universes>
</lattice>
<!-- 4 x 4 core -->
<lattice id="6">
<type>rectangular</type>
<dimension>2 2</dimension>
<origin>-2.0 -2.0</origin>
<width>2.0 2.0</width>
<universes>
5 5
5 5
</universes>
</lattice>
<surface id="1" type="x-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="2.0" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="2.0" boundary="vacuum" />
<surface id="5" type="z-cylinder" coeffs="0.0 0.0 0.4" />
<surface id="6" type="z-cylinder" coeffs="0.0 0.0 0.3" />
<surface id="7" type="z-cylinder" coeffs="0.0 0.0 0.2" />
</geometry>

19
examples/lattice/nested/materials.xml Normal file
View file

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

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

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

9
examples/lattice/settings.xml → examples/lattice/nested/settings.xml
View file

@ -1,11 +1,6 @@
<?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>20</cycles>
@ -13,10 +8,12 @@
<particles>10000</particles>
</criticality>
<verbosity>7</verbosity>
<!-- Starting source -->
<source>
<type>box</type>
<coeffs>-4 -4 -4 4 4 4</coeffs>
<coeffs>-1 -1 -1 1 1 1</coeffs>
</source>
</settings>

19
examples/lattice/nested/tallies.xml Normal file
View file

@ -0,0 +1,19 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>rectangular</type>
<dimension>4 4</dimension>
<origin>-2.0 -2.0</origin>
<width>1.0 1.0</width>
</mesh>
<tally id="1">
<filters>
<mesh>1</mesh>
</filters>
<macros>total</macros>
</tally>
</tallies>

33
examples/lattice/simple/geometry.xml Normal file
View file

@ -0,0 +1,33 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="5" surfaces="1 -2 3 -4" />
<cell id="101" universe="1" material="1" surfaces="-5" />
<cell id="102" universe="1" material="2" surfaces="5" />
<cell id="201" universe="2" material="1" surfaces="-6" />
<cell id="202" universe="2" material="2" surfaces="6" />
<cell id="301" universe="3" material="1" surfaces="-7" />
<cell id="302" universe="3" material="2" surfaces="7" />
<lattice id="5">
<type>rectangular</type>
<dimension>4 4</dimension>
<origin>-2.0 -2.0</origin>
<width>1.0 1.0</width>
<universes>
1 2 1 2
2 3 2 3
1 2 1 2
2 3 2 3
</universes>
</lattice>
<surface id="1" type="x-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="2.0" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="2.0" boundary="vacuum" />
<surface id="5" type="z-cylinder" coeffs="0.0 0.0 0.4" />
<surface id="6" type="z-cylinder" coeffs="0.0 0.0 0.3" />
<surface id="7" type="z-cylinder" coeffs="0.0 0.0 0.2" />
</geometry>

19
examples/lattice/simple/materials.xml Normal file
View file

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

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

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

19
examples/lattice/simple/settings.xml Normal file
View file

@ -0,0 +1,19 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Parameters for criticality calculation -->
<criticality>
<cycles>20</cycles>
<inactive>10</inactive>
<particles>10000</particles>
</criticality>
<verbosity>7</verbosity>
<!-- Starting source -->
<source>
<type>box</type>
<coeffs>-1 -1 -1 1 1 1</coeffs>
</source>
</settings>

19
examples/lattice/simple/tallies.xml Normal file
View file

@ -0,0 +1,19 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>rectangular</type>
<dimension>4 4</dimension>
<origin>-2.0 -2.0</origin>
<width>1.0 1.0</width>
</mesh>
<tally id="1">
<filters>
<mesh>1</mesh>
</filters>
<macros>total</macros>
</tally>
</tallies>

17
src/DEPENDENCIES
View file

@ -19,6 +19,7 @@ cross_section.o: datatypes_header.o
cross_section.o: endf.o
cross_section.o: error.o
cross_section.o: fileio.o
cross_section.o: fission.o
cross_section.o: global.o
cross_section.o: material_header.o
cross_section.o: output.o
@ -87,9 +88,9 @@ initialize.o: geometry_header.o
initialize.o: global.o
initialize.o: input_xml.o
initialize.o: logging.o
initialize.o: mcnp_random.o
initialize.o: mpi_routines.o
initialize.o: output.o
initialize.o: random_lcg.o
initialize.o: source.o
initialize.o: string.o
initialize.o: tally.o
@ -116,6 +117,7 @@ interpolation.o: endf_header.o
interpolation.o: error.o
interpolation.o: global.o
interpolation.o: search.o
interpolation.o: string.o
logging.o: constants.o
logging.o: global.o
@ -124,12 +126,12 @@ main.o: cmfd_execute.o
main.o: constants.o
main.o: global.o
main.o: initialize.o
main.o: mcnp_random.o
main.o: mpi_routines.o
main.o: output.o
main.o: particle_header.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
@ -140,10 +142,11 @@ mesh.o: mesh_header.o
mpi_routines.o: constants.o
mpi_routines.o: error.o
mpi_routines.o: global.o
mpi_routines.o: mcnp_random.o
mpi_routines.o: output.o
mpi_routines.o: particle_header.o
mpi_routines.o: random_lcg.o
mpi_routines.o: tally_header.o
mpi_routines.o: timing.o
output.o: constants.o
output.o: datatypes.o
@ -151,10 +154,12 @@ output.o: endf.o
output.o: geometry_header.o
output.o: global.o
output.o: mesh_header.o
output.o: particle_header.o
output.o: string.o
output.o: tally_header.o
particle_header.o: constants.o
particle_header.o: geometry_header.o
physics.o: constants.o
physics.o: cross_section_header.o
@ -165,9 +170,9 @@ physics.o: geometry.o
physics.o: geometry_header.o
physics.o: global.o
physics.o: interpolation.o
physics.o: mcnp_random.o
physics.o: output.o
physics.o: particle_header.o
physics.o: random_lcg.o
physics.o: search.o
physics.o: string.o
physics.o: tally.o
@ -186,11 +191,13 @@ search.o: global.o
source.o: bank_header.o
source.o: constants.o
source.o: cross_section_header.o
source.o: error.o
source.o: global.o
source.o: mcnp_random.o
source.o: output.o
source.o: particle_header.o
source.o: physics.o
source.o: random_lcg.o
source.o: string.o
string.o: constants.o
string.o: error.o

146
src/Makefile
View file

@ -20,27 +20,60 @@ DEBUG = no
PROFILE = no
OPTIMIZE = no
USE_MPI = no
USE_OPENMP = no
USE_COARRAY = no
#===============================================================================
# Compiler Options
#===============================================================================
# Intel Fortran compiler options
#===============================================================================
ifeq ($(COMPILER),intel)
F90 = ifort
F90FLAGS := -fpp -warn
F90FLAGS := -fpp -warn -assume byterecl
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -traceback -ftrapuv -fp-stack-check -check all
LDFLAGS += -g
endif
# Profiling
ifeq ($(PROFILE),yes)
F90FLAGS += -pg
LDFLAGS += -pg
endif
# Optimization
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -O3 -ipo
endif
endif
#===============================================================================
# GNU Fortran compiler options
#===============================================================================
ifeq ($(COMPILER),gfortran)
F90 = gfortran
F90FLAGS := -cpp -Wall
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -pedantic -std=f2008 -fbacktrace -fbounds-check \
-ffpe-trap=invalid,zero,overflow,underflow
LDFLAGS += -g
endif
# Profiling
ifeq ($(PROFILE),yes)
F90FLAGS += -pg
LDFLAGS += -pg
endif
# Optimization
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -O3
endif
endif
# use petsc compiler
@ -53,36 +86,78 @@ ifeq ($(COMPILER),petsc)
LDFLAGS = ${PETSC_SYS_LIB}
endif
# Set compiler flags for debugging
#===============================================================================
# PGI compiler options
#===============================================================================
ifeq ($(DEBUG),yes)
F90FLAGS += -g
LDFLAGS += -g
ifeq ($(COMPILER),intel)
F90FLAGS += -traceback -ftrapuv -fp-stack-check -check all
ifeq ($(COMPILER),pgi)
F90 = pgf90
F90FLAGS := -Mpreprocess -DNO_F2008 -Minform=inform
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -Mbounds -Mchkptr -Mchkstk -traceback
LDFLAGS += -g
endif
ifeq ($(COMPILER),gfortran)
F90FLAGS += -pedantic -std=f2008 -fbacktrace -fbounds-check \
-ffpe-trap=invalid,zero,overflow,underflow
# Profiling
ifeq ($(PROFILE),yes)
F90FLAGS += -pg
LDFLAGS += -pg
endif
# Optimization
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -fast -Mipa
endif
endif
# Set compiler flags for profiling
#===============================================================================
# IBM XL compiler options
#===============================================================================
ifeq ($(PROFILE),yes)
F90FLAGS += -pg
LDFLAGS += -pg
endif
ifeq ($(COMPILER),ibm)
F90 = xlf2003
F90FLAGS := -WF,-DNO_F2008
# Set compiler flags for high optimization
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -C -qflag=i:i -u
LDFLAGS += -g
endif
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -O3
ifeq ($(COMPILER),intel)
F90FLAGS += -ipo
# Profiling
ifeq ($(PROFILE),yes)
F90FLAGS += -p
LDFLAGS += -p
endif
# Optimization
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -O3
endif
endif
#===============================================================================
# Cray compiler options
#===============================================================================
ifeq ($(COMPILER),cray)
F90 = ftn
F90FLAGS := -e Z -m 0
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -R abcnsp -O0
LDFLAGS += -g
endif
endif
#===============================================================================
# Miscellaneous compiler options
#===============================================================================
# Use MPI for parallelism
ifeq ($(USE_MPI),yes)
@ -91,17 +166,14 @@ ifeq ($(USE_MPI),yes)
F90FLAGS += -DMPI
endif
# Use OpenMP for shared-memory parallelism
#===============================================================================
# Special options for ORNL Jaguar supercomputer
#===============================================================================
ifeq ($(USE_OPENMP),yes)
F90FLAGS += -openmp
LDFLAGS += -openmp
endif
# Use Fortran 2008 Coarrays for parallelism
ifeq ($(USE_COARRAY),yes)
F90FLAGS += -coarray
HOSTNAME = $(shell hostname)
ifneq (,$(findstring jaguar,$(HOSTNAME)))
F90 = ftn
F90FLAGS += -DMPI
endif
#===============================================================================
@ -126,10 +198,10 @@ neat:
# Rules
#===============================================================================
.SUFFIXES: .f90 .o
.SUFFIXES: .F90 .o
.PHONY: all xml-fortran clean neat distclean
%.o: %.f90
%.o: %.F90
$(F90) -Ixml-fortran -Ixml-fortran/templates $(F90FLAGS) -c $<
#===============================================================================

65
src/Makefile.jaguar
View file

@ -1,65 +0,0 @@
program = openmc
#===============================================================================
# Object Files
#===============================================================================
include OBJECTS
#===============================================================================
# User Options
#===============================================================================
DEBUG = no
PROFILE = no
OPTIMIZE = no
#===============================================================================
# Compiler Options
#===============================================================================
F90 = ftn
F90FLAGS = -cpp -DMPI
LDFLAGS =
ifeq ($(DEBUG),yes)
F90FLAGS += -g -traceback -ftrapuv -fp-stack-check -check all
LDFLAGS += -g
endif
ifeq ($(PROFILE),yes)
F90FLAGS += -pg
LDFLAGS += -pg
endif
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -ipo -O3
endif
#===============================================================================
# Targets
#===============================================================================
all: $(program)
$(program): $(objects)
$(F90) $(objects) -o $@ $(LDFLAGS)
clean:
@rm -f *.o *.mod $(program)
neat:
@rm -f *.o *.mod
#===============================================================================
# Rules
#===============================================================================
.SUFFIXES: .f90 .o
.PHONY: all clean neat
%.o: %.f90
$(F90) $(F90FLAGS) -c $<
#===============================================================================
# Dependencies
#===============================================================================
include DEPENDENCIES

2
src/OBJECTS
View file

@ -23,7 +23,6 @@ input_xml.o \
logging.o \
main.o \
material_header.o \
mcnp_random.o \
mesh_header.o \
mesh.o \
mpi_routines.o \
@ -31,6 +30,7 @@ output.o \
particle_header.o \
physics.o \
plot.o \
random_lcg.o \
search.o \
source.o \
source_header.o \

0
src/bank_header.f90 → src/bank_header.F90
View file

118
src/constants.f90 → src/constants.F90
View file

@ -2,12 +2,42 @@ module constants
implicit none
! Versioning numbers
! ============================================================================
! VERSIONING NUMBERS
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 3
integer, parameter :: VERSION_RELEASE = 3
! Physical constants
! ============================================================================
! ADJUSTABLE PARAMETERS
! NOTE: This is the only section of the constants module that should ever be
! adjusted. Modifying constants in other sections may cause the code to fail.
! Monoatomic ideal-gas scattering treatment threshold
real(8), parameter :: FREE_GAS_THRESHOLD = 400.0
! Used for surface current tallies
real(8), parameter :: TINY_BIT = 1e-8_8
! User for precision in geometry
real(8), parameter :: FP_PRECISION = 1e-7_8
! Maximum number of collisions/crossings
integer, parameter :: MAX_EVENTS = 10000
integer, parameter :: MAX_SAMPLE = 100000
! Maximum number of words in a single line, length of line, and length of
! single word
integer, parameter :: MAX_WORDS = 500
integer, parameter :: MAX_LINE_LEN = 250
integer, parameter :: MAX_WORD_LEN = 150
integer, parameter :: MAX_FILE_LEN = 255
! ============================================================================
! PHYSICAL CONSTANTS
real(8), parameter :: &
PI = 3.1415926535898_8, & ! pi
MASS_NEUTRON = 1.0086649156, & ! mass of a neutron
@ -20,8 +50,8 @@ module constants
ONE = 1.0_8, &
TWO = 2.0_8
! Monoatomic ideal-gas scattering treatment threshold
real(8), parameter :: FREE_GAS_THRESHOLD = 400.0
! ============================================================================
! GEOMETRY-RELATED CONSTANTS
! Boundary conditions
integer, parameter :: &
@ -70,23 +100,8 @@ module constants
& SENSE_POSITIVE = 1, &
& SENSE_NEGATIVE = -1
! Used for surface current tallies
real(8), parameter :: TINY_BIT = 1e-8
! Codes for read errors -- better hope these numbers are never used in an
! input file!
integer, parameter :: ERROR_INT = -huge(0)
real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8
! 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
integer, parameter :: &
PROB_SOURCE = 1, & ! External source problem
PROB_CRITICALITY = 2 ! Criticality problem
! ============================================================================
! CROSS SECTION RELATED CONSTANTS
! Interpolation flag
integer, parameter :: &
@ -175,6 +190,23 @@ module constants
ACE_THERMAL = 2, & ! thermal S(a,b) scattering data
ACE_DOSIMETRY = 3 ! dosimetry cross sections
! Fission neutron emission (nu) type
integer, parameter :: &
NU_NONE = 0, & ! No nu values (non-fissionable)
NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
NU_TABULAR = 2 ! Nu values given by tabular distribution
! Cross section filetypes
integer, parameter :: &
ASCII = 1, & ! ASCII cross section file
BINARY = 2 ! Binary cross section file
! Maximum number of partial fission reactions
integer, parameter :: PARTIAL_FISSION_MAX = 4
! ============================================================================
! TALLY-RELATED CONSTANTS
! Tally macro reactions
integer, parameter :: N_MACRO_TYPES = 15
integer, parameter :: &
@ -209,6 +241,14 @@ module constants
T_ENERGYIN = 7, &
T_ENERGYOUT = 8
! Filter types for surface current tallies
integer, parameter :: &
TS_MESH_X = 1, &
TS_MESH_Y = 2, &
TS_MESH_Z = 3, &
TS_ENERGYIN = 4, &
TS_SURFACE = 5
! Tally surface current directions
integer, parameter :: &
IN_RIGHT = 1, &
@ -218,26 +258,26 @@ module constants
IN_TOP = 5, &
OUT_TOP = 6
! Fission neutron emission (nu) type
! ============================================================================
! MISCELLANEOUS CONSTANTS
! indicates that an array index hasn't been set
integer, parameter :: NONE = 0
! Codes for read errors -- better hope these numbers are never used in an
! input file!
integer, parameter :: ERROR_INT = -huge(0)
real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8
! Source types
integer, parameter :: &
NU_NONE = 0, & ! No nu values (non-fissionable)
NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
NU_TABULAR = 2 ! Nu values given by tabular distribution
SRC_BOX = 1, & ! Source in a rectangular prism
SRC_CELL = 2, & ! Source in a cell
SRC_SURFACE = 3 ! Source on a surface
! Cross section filetypes
integer, parameter :: &
ASCII = 1, & ! ASCII cross section file
BINARY = 2 ! Binary cross section file
! Maximum number of partial fission reactions
integer, parameter :: PARTIAL_FISSION_MAX = 4
! Maximum number of words in a single line, length of line, and length of
! single word
integer, parameter :: MAX_WORDS = 500
integer, parameter :: MAX_LINE_LEN = 250
integer, parameter :: MAX_WORD_LEN = 150
integer, parameter :: MAX_FILE_LEN = 255
integer, parameter :: &
PROB_SOURCE = 1, & ! External source problem
PROB_CRITICALITY = 2 ! Criticality problem
! Unit numbers
integer, parameter :: UNIT_LOG = 11 ! unit # for writing log file

252
src/cross_section.f90 → src/cross_section.F90
View file

@ -9,6 +9,7 @@ module cross_section
use endf, only: reaction_name
use error, only: fatal_error
use fileio, only: read_line, skip_lines
use fission, only: nu_total
use global
use material_header, only: Material
use output, only: write_message
@ -39,11 +40,11 @@ contains
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: index ! index in xs_listings array
integer :: index_list ! index in xs_listings array
integer :: index_nuclides ! index in nuclides
integer :: index_sab ! index in sab_tables
character(10) :: name ! name of isotope, e.g. 92235.03c
character(10) :: alias ! alias of nuclide, e.g. U-235.03c
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
type(SAB_Table), pointer :: sab => null()
@ -76,9 +77,9 @@ contains
! Find index in xs_listing and set the name and alias according to the
! listing
index = dict_get_key(xs_listing_dict, name)
name = xs_listings(index) % name
alias = xs_listings(index) % alias
index_list = dict_get_key(xs_listing_dict, name)
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
! If this nuclide hasn't been encountered yet, we need to add its name
! and alias to the nuclide_dict
@ -107,9 +108,9 @@ contains
! Find index in xs_listing and set the name and alias according to the
! listing
index = dict_get_key(xs_listing_dict, name)
name = xs_listings(index) % name
alias = xs_listings(index) % alias
index_list = dict_get_key(xs_listing_dict, name)
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
! If this S(a,b) table hasn't been encountered yet, we need to add its
! name and alias to the sab_dict
@ -146,12 +147,12 @@ contains
name = mat % names(j)
if (.not. dict_has_key(already_read, name)) then
index = dict_get_key(xs_listing_dict, name)
index_list = dict_get_key(xs_listing_dict, name)
index_nuclides = dict_get_key(nuclide_dict, name)
name = xs_listings(index) % name
alias = xs_listings(index) % alias
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
call read_ace_table(index_nuclides, index)
call read_ace_table(index_nuclides, index_list)
call dict_add_key(already_read, name, 0)
call dict_add_key(already_read, alias, 0)
@ -162,10 +163,10 @@ contains
name = mat % sab_name
if (.not. dict_has_key(already_read, name)) then
index = dict_get_key(xs_listing_dict, name)
index_sab = dict_get_key(sab_dict, name)
index_list = dict_get_key(xs_listing_dict, name)
index_sab = dict_get_key(sab_dict, name)
call read_ace_table(index_sab, index)
call read_ace_table(index_sab, index_list)
call dict_add_key(already_read, name, 0)
end if
@ -210,10 +211,10 @@ contains
! subroutines to parse the actual data.
!===============================================================================
subroutine read_ace_table(index_table, index)
subroutine read_ace_table(index_table, index_list)
integer, intent(in) :: index_table ! index in nuclides/sab_tables
integer, intent(in) :: index ! index in xs_listings
integer, intent(in) :: index_list ! index in xs_listings
integer :: i ! loop index for XSS records
integer :: j, j1, j2 ! indices in XSS
@ -230,7 +231,7 @@ contains
logical :: file_exists ! does ACE library exist?
character(7) :: readable ! is ACE library readable?
character(10) :: name ! name of ACE table
character(10) :: date ! date ACE library was processed
character(10) :: date_ ! date ACE library was processed
character(10) :: mat ! material identifier
character(70) :: comment ! comment for ACE table
character(MAX_FILE_LEN) :: filename ! path to ACE cross section library
@ -239,7 +240,7 @@ contains
type(XsListing), pointer :: listing => null()
! determine path, record length, and location of table
listing => xs_listings(index)
listing => xs_listings(index_list)
filename = listing % path
record_length = listing % recl
location = listing % location
@ -273,7 +274,7 @@ contains
end do
! Read first line of header
read(UNIT=in, FMT='(A10,2E12.0,1X,A10)') name, awr, kT, date
read(UNIT=in, FMT='(A10,2E12.0,1X,A10)') name, awr, kT, date_
! Read more header and NXS and JXS
read(UNIT=in, FMT=100) comment, mat, &
@ -296,7 +297,7 @@ contains
ACCESS='direct', RECL=record_length)
! Read all header information
read(UNIT=in, REC=location) name, awr, kT, date, &
read(UNIT=in, REC=location) name, awr, kT, date_, &
comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS
! determine table length
@ -322,12 +323,20 @@ contains
nuc % kT = kT
nuc % zaid = NXS(2)
! read all blocks
call read_esz(nuc)
call read_nu_data(nuc)
call read_reactions(nuc)
call read_angular_dist(nuc)
call read_energy_dist(nuc)
call read_unr_res(nuc)
! for fissionable nuclides, precalculate microscopic nu-fission cross
! sections so that we don't need to call the nu_total function during
! cross section lookups
if (nuc % fissionable) call generate_nu_fission(nuc)
case (ACE_THERMAL)
sab => sab_tables(index_table)
sab % name = name
@ -366,6 +375,7 @@ contains
allocate(nuc % total(NE))
allocate(nuc % elastic(NE))
allocate(nuc % fission(NE))
allocate(nuc % nu_fission(NE))
allocate(nuc % absorption(NE))
allocate(nuc % heating(NE))
@ -373,6 +383,7 @@ contains
nuc % total = ZERO
nuc % elastic = ZERO
nuc % fission = ZERO
nuc % nu_fission = ZERO
nuc % absorption = ZERO
nuc % heating = ZERO
@ -414,7 +425,7 @@ contains
integer :: LED ! location of energy distribution locators
integer :: LDIS ! location of all energy distributions
integer :: LOCC ! location of energy distributions for given MT
integer :: loc ! locator
integer :: lc ! locator
integer :: length ! length of data to allocate
type(DistEnergy), pointer :: edist => null()
@ -464,21 +475,21 @@ contains
! =======================================================================
! PROMPT AND TOTAL NU DATA -- read prompt data first
KNU = JXS2 + 1
LNU = XSS(KNU)
LNU = int(XSS(KNU))
if (LNU == 1) then
! Polynomial data
nuc % nu_p_type = NU_POLYNOMIAL
! allocate determine how many coefficients for polynomial
NC = XSS(KNU+1)
NC = int(XSS(KNU+1))
length = NC + 1
elseif (LNU == 2) then
! Tabular data
nuc % nu_p_type = NU_TABULAR
! determine number of interpolation regions and number of energies
NR = XSS(KNU+1)
NE = XSS(KNU+2+2*NR)
NR = int(XSS(KNU+1))
NE = int(XSS(KNU+2+2*NR))
length = 2 + 2*NR + 2*NE
end if
@ -490,8 +501,8 @@ contains
nuc % nu_p_data = get_real(length)
! Now read total nu data
KNU = JXS2 + abs(XSS(JXS2)) + 1
LNU = XSS(KNU)
KNU = JXS2 + int(abs(XSS(JXS2))) + 1
LNU = int(XSS(KNU))
if (LNU == 1) then
! Polynomial data
nuc % nu_t_type = NU_POLYNOMIAL
@ -550,7 +561,7 @@ contains
! Loop over all delayed neutron precursor groups
do i = 1, NPCR
! find location of energy distribution data
LOCC = XSS(LED + i - 1)
LOCC = int(XSS(LED + i - 1))
! read energy distribution data
edist => nuc % nu_d_edist(i)
@ -562,10 +573,10 @@ contains
! determine length of all precursor constants/yields/interp data
length = 0
loc = JXS(25)
lc = JXS(25)
do i = 1, NPCR
NR = XSS(loc + length + 1)
NE = XSS(loc + length + 2 + 2*NR)
NR = int(XSS(lc + length + 1))
NE = int(XSS(lc + length + 2 + 2*NR))
length = length + 3 + 2*NR + 2*NE
end do
@ -573,7 +584,7 @@ contains
allocate(nuc % nu_d_precursor_data(length))
! read delayed neutron precursor data
XSS_index = loc
XSS_index = lc
nuc % nu_d_precursor_data = get_real(length)
else
@ -640,17 +651,17 @@ contains
rxn => nuc % reactions(i+1)
! read MT number, Q-value, and neutrons produced
rxn % MT = XSS(LMT + i - 1)
rxn % MT = int(XSS(LMT + i - 1))
rxn % Q_value = XSS(JXS4 + i - 1)
rxn % TY = XSS(JXS5 + i - 1)
rxn % TY = int(XSS(JXS5 + i - 1))
! read starting energy index
LOCA = XSS(LXS + i - 1)
IE = XSS(JXS7 + LOCA - 1)
LOCA = int(XSS(LXS + i - 1))
IE = int(XSS(JXS7 + LOCA - 1))
rxn % IE = IE
! read number of energies cross section values
NE = XSS(JXS7 + LOCA)
NE = int(XSS(JXS7 + LOCA))
allocate(rxn % sigma(NE))
XSS_index = JXS7 + LOCA + 1
rxn % sigma = get_real(NE)
@ -727,7 +738,7 @@ contains
rxn => nuc%reactions(i)
! find location of angular distribution
LOCB = XSS(JXS8 + i - 1)
LOCB = int(XSS(JXS8 + i - 1))
if (LOCB == -1) then
! Angular distribution data are specified through LAWi = 44 in the DLW
! block
@ -740,7 +751,7 @@ contains
rxn % has_angle_dist = .true.
! allocate space for incoming energies and locations
NE = XSS(JXS9 + LOCB - 1)
NE = int(XSS(JXS9 + LOCB - 1))
rxn % adist % n_energy = NE
allocate(rxn % adist % energy(NE))
allocate(rxn % adist % type(NE))
@ -765,7 +776,7 @@ contains
elseif (LC < 0) then
! tabular distribution
rxn % adist % type(j) = ANGLE_TABULAR
NP = XSS(JXS9 + abs(LC))
NP = int(XSS(JXS9 + abs(LC)))
length = length + 2 + 3*NP
end if
end do
@ -800,7 +811,7 @@ contains
integer :: LED ! location of energy distribution locators
integer :: LOCC ! location of energy distributions for given MT
integer :: i ! loop index
type(Reaction), pointer :: rxn => null()
type(Reaction), pointer :: rxn => null()
LED = JXS(10)
@ -810,7 +821,7 @@ contains
rxn % has_energy_dist = .true.
! find location of energy distribution data
LOCC = XSS(LED + i - 1)
LOCC = int(XSS(LED + i - 1))
! allocate energy distribution
allocate(rxn % edist)
@ -839,7 +850,7 @@ contains
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: IDAT ! location of first energy distribution for given MT
integer :: loc ! locator
integer :: lc ! locator
integer :: length ! length of data to allocate
integer :: length_interp_data ! length of interpolation data
@ -851,10 +862,10 @@ contains
end if
! locator for next law and information on this law
LNW = XSS(LDIS + loc_law - 1)
LAW = XSS(LDIS + loc_law)
IDAT = XSS(LDIS + loc_law + 1)
NR = XSS(LDIS + loc_law + 2)
LNW = int(XSS(LDIS + loc_law - 1))
LAW = int(XSS(LDIS + loc_law))
IDAT = int(XSS(LDIS + loc_law + 1))
NR = int(XSS(LDIS + loc_law + 2))
edist % law = LAW
edist % p_valid % n_regions = NR
@ -867,12 +878,12 @@ contains
! read ENDF interpolation parameters
XSS_index = LDIS + loc_law + 3
if (NR > 0) then
edist % p_valid % nbt = get_real(NR)
edist % p_valid % int = get_real(NR)
edist % p_valid % nbt = int(get_real(NR))
edist % p_valid % int = int(get_real(NR))
end if
! allocate space for law validity data
NE = XSS(LDIS + loc_law + 3 + 2*NR)
NE = int(XSS(LDIS + loc_law + 3 + 2*NR))
edist % p_valid % n_pairs = NE
allocate(edist % p_valid % x(NE))
allocate(edist % p_valid % y(NE))
@ -885,16 +896,16 @@ contains
edist % p_valid % y = get_real(NE)
! Set index to beginning of IDAT array
loc = LDIS + IDAT - 2
lc = LDIS + IDAT - 2
! determine length of energy distribution
length = length_energy_dist(loc, LAW, loc_law, length_interp_data)
length = length_energy_dist(lc, LAW, loc_law, length_interp_data)
! allocate secondary energy distribution array
allocate(edist % data(length))
! read secondary energy distribution
XSS_index = loc + 1
XSS_index = lc + 1
edist % data = get_real(length)
! read next energy distribution if present
@ -910,9 +921,9 @@ contains
! distribution array based on the secondary energy law and location in XSS
!===============================================================================
function length_energy_dist(loc, law, LOCC, lid) result(length)
function length_energy_dist(lc, law, LOCC, lid) result(length)
integer, intent(in) :: loc ! location in XSS array
integer, intent(in) :: lc ! location in XSS array
integer, intent(in) :: law ! energy distribution law
integer, intent(in) :: LOCC ! location of energy distribution
integer, intent(in) :: lid ! length of interpolation data
@ -934,9 +945,9 @@ contains
select case (law)
case (1)
! Tabular equiprobable energy bins
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NP = XSS(loc + 3 + 2*NR + NE)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + NE))
length = 3 + 2*NR + NE + 3*NP*NE
case (2)
@ -949,72 +960,72 @@ contains
case (4)
! Continuous tabular distribution
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! determine length
NP = XSS(loc + length + 2)
NP = int(XSS(lc + length + 2))
length = length + 2 + 3*NP
! adjust location for this block
j = loc + 2 + 2*NR + NE + i
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
case (5)
! General evaporation spectrum
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NP = XSS(loc + 3 + 2*NR + 2*NE)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + 2*NE))
length = 3 + 2*NR + 2*NE + NP
case (7)
! Maxwell fission spectrum
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
case (9)
! Evaporation spectrum
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
case (11)
! Watt spectrum
NRa = XSS(loc + 1)
NEa = XSS(loc + 2 + 2*NRa)
NRb = XSS(loc + 3 + 2*(NRa+NEa))
NEb = XSS(loc + 4 + 2*(NRa+NEa+NRb))
NRa = int(XSS(lc + 1))
NEa = int(XSS(lc + 2 + 2*NRa))
NRb = int(XSS(lc + 3 + 2*(NRa+NEa)))
NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb)))
length = 5 + 2*(NRa + NEa + NRb + NEb)
case (44)
! Kalbach-Mann correlated scattering
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
NP = XSS(loc + length + 2)
NP = int(XSS(lc + length + 2))
length = length + 2 + 5*NP
! adjust location for this block
j = loc + 2 + 2*NR + NE + i
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
case (61)
! Correlated energy and angle distribution
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! outgoing energy distribution
NP = XSS(loc + length + 2)
NP = int(XSS(lc + length + 2))
! adjust locators for angular distribution
do j = 1, NP
k = loc + length + 2 + 3*NP + j
k = lc + length + 2 + 3*NP + j
if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
end do
@ -1022,12 +1033,12 @@ contains
do j = 1, NP
! outgoing angle distribution -- NMU here is actually
! referred to as NP in the MCNP documentation
NMU = XSS(loc + length + 2)
NMU = int(XSS(lc + length + 2))
length = length + 2 + 3*NMU
end do
! adjust locators for energy distribution
j = loc + 2 + 2*NR + NE + i
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
@ -1037,9 +1048,9 @@ contains
case (67)
! Laboratory energy-angle law
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NMU = XSS(loc + 4 + 2*NR + 2*NE)
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NMU = int(XSS(lc + 4 + 2*NR + 2*NE))
length = 4 + 2*(NR + NE + NMU)
end select
@ -1055,7 +1066,7 @@ contains
type(Nuclide), pointer :: nuc
integer :: JXS23 ! location of URR data
integer :: loc ! locator
integer :: lc ! locator
integer :: N ! # of incident energies
integer :: M ! # of probabilities
integer :: i ! index over incoming energies
@ -1070,19 +1081,19 @@ contains
nuc % urr_present = .true.
allocate(nuc % urr_data)
allocate(nuc % urr_data % params(6))
loc = JXS23
lc = JXS23
else
nuc % urr_present = .false.
return
end if
! read parameters
nuc % urr_data % params(1) = XSS(loc) ! # of incident energies
nuc % urr_data % params(2) = XSS(loc + 1) ! # of probabilities
nuc % urr_data % params(3) = XSS(loc + 2) ! interpolation parameter
nuc % urr_data % params(4) = XSS(loc + 3) ! inelastic competition flag
nuc % urr_data % params(5) = XSS(loc + 4) ! other absorption flag
nuc % urr_data % params(6) = XSS(loc + 5) ! factors flag
nuc % urr_data % params(1) = int(XSS(lc)) ! # of incident energies
nuc % urr_data % params(2) = int(XSS(lc + 1)) ! # of probabilities
nuc % urr_data % params(3) = int(XSS(lc + 2)) ! interpolation parameter
nuc % urr_data % params(4) = int(XSS(lc + 3)) ! inelastic competition flag
nuc % urr_data % params(5) = int(XSS(lc + 4)) ! other absorption flag
nuc % urr_data % params(6) = int(XSS(lc + 5)) ! factors flag
! allocate incident energies and probability tables
N = nuc % urr_data % params(1)
@ -1091,7 +1102,7 @@ contains
allocate(nuc % urr_data % prob(N,6,M))
! read incident energies
XSS_index = loc + 6
XSS_index = lc + 6
nuc % urr_data % energy = get_real(N)
! read probability tables
@ -1106,6 +1117,33 @@ contains
end subroutine read_unr_res
!===============================================================================
! GENERATE_NU_FISSION precalculates the microscopic nu-fission cross section for
! a given nuclide. This is done so that the nu_total function does not need to
! be called during cross section lookups.
!===============================================================================
subroutine generate_nu_fission(nuc)
type(Nuclide), pointer :: nuc
integer :: i ! index on nuclide energy grid
real(8) :: E ! energy
real(8) :: nu ! # of neutrons per fission
do i = 1, nuc % n_grid
! determine energy
E = nuc % energy(i)
! determine total nu at given energy
nu = nu_total(nuc, E)
! determine nu-fission microscopic cross section
nuc % nu_fission(i) = nu * nuc % fission(i)
end do
end subroutine generate_nu_fission
!===============================================================================
! READ_THERMAL_DATA reads elastic and inelastic cross sections and corresponding
! secondary energy/angle distributions derived from experimental S(a,b)
@ -1120,7 +1158,7 @@ contains
integer :: i ! index for incoming energies
integer :: j ! index for outgoing energies
integer :: k ! index for outoging angles
integer :: loc ! location in XSS array
integer :: lc ! location in XSS array
integer :: NE_in ! number of incoming energies
integer :: NE_out ! number of outgoing energies
integer :: NMU ! number of outgoing angles
@ -1130,7 +1168,7 @@ contains
table % secondary_mode = NXS(7)
! read number of inelastic energies and allocate arrays
NE_in = XSS(JXS(1))
NE_in = int(XSS(JXS(1)))
table % n_inelastic_e_in = NE_in
allocate(table % inelastic_e_in(NE_in))
allocate(table % inelastic_sigma(NE_in))
@ -1153,26 +1191,26 @@ contains
allocate(table % inelastic_mu(NMU, NE_out, NE_in))
! read outgoing energy/angle distribution for inelastic scattering
loc = JXS(3) - 1
lc = JXS(3) - 1
do i = 1, NE_in
do j = 1, NE_out
! read outgoing energy
table % inelastic_e_out(j,i) = XSS(loc + 1)
table % inelastic_e_out(j,i) = XSS(lc + 1)
! read outgoing angles for this outgoing energy
do k = 1, NMU
table % inelastic_mu(k,j,i) = XSS(loc + 1 + k)
table % inelastic_mu(k,j,i) = XSS(lc + 1 + k)
end do
! advance pointer
loc = loc + 1 + NMU
lc = lc + 1 + NMU
end do
end do
! read number of elastic energies and allocate arrays
JXS4 = JXS(4)
if (JXS4 /= 0) then
NE_in = XSS(JXS4)
NE_in = int(XSS(JXS4))
table % n_elastic_e_in = NE_in
allocate(table % elastic_e_in(NE_in))
allocate(table % elastic_P(NE_in))
@ -1202,12 +1240,12 @@ contains
! read equiprobable outgoing cosines for elastic scattering each
! incoming energy
if (JXS4 /= 0 .and. NMU /= 0) then
loc = JXS(6) - 1
lc = JXS(6) - 1
do i = 1, NE_in
do j = 1, NMU
table % elastic_mu(j,i) = XSS(loc + j)
table % elastic_mu(j,i) = XSS(lc + j)
end do
loc = loc + NMU
lc = lc + NMU
end do
end if

99
src/cross_section_header.f90 → src/cross_section_header.F90
View file

@ -67,28 +67,29 @@ module cross_section_header
!===============================================================================
type Nuclide
character(20) :: name
integer :: zaid
real(8) :: awr
real(8) :: kT
character(10) :: name ! name of nuclide, e.g. 92235.03c
integer :: zaid ! Z and A identifier, e.g. 92235
real(8) :: awr ! weight of nucleus in neutron masses
real(8) :: kT ! temperature in MeV (k*T)
! Energy grid information
integer :: n_grid
integer, allocatable :: grid_index(:)
real(8), allocatable :: energy(:)
integer :: n_grid ! # of nuclide grid points
integer, allocatable :: grid_index(:) ! pointers to union grid
real(8), allocatable :: energy(:) ! energy values corresponding to xs
! Cross sections
real(8), allocatable :: total(:)
real(8), allocatable :: elastic(:)
real(8), allocatable :: fission(:)
real(8), allocatable :: absorption(:)
real(8), allocatable :: heating(:)
! Microscopic cross sections
real(8), allocatable :: total(:) ! total cross section
real(8), allocatable :: elastic(:) ! elastic scattering
real(8), allocatable :: fission(:) ! fission
real(8), allocatable :: nu_fission(:) ! neutron production
real(8), allocatable :: absorption(:) ! absorption (MT > 100)
real(8), allocatable :: heating(:) ! heating
! Fission information
logical :: fissionable
logical :: has_partial_fission
integer :: n_fission
integer, allocatable :: index_fission(:)
logical :: fissionable ! nuclide is fissionable?
logical :: has_partial_fission ! nuclide has partial fission reactions?
integer :: n_fission ! # of fission reactions
integer, allocatable :: index_fission(:) ! indices in reactions
! Total fission neutron emission
integer :: nu_t_type
@ -100,7 +101,7 @@ module cross_section_header
! Delayed fission neutron emission
integer :: nu_d_type
integer :: n_precursor
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
type(DistEnergy), pointer :: nu_d_edist(:) => null()
@ -110,7 +111,7 @@ module cross_section_header
type(UrrData), pointer :: urr_data => null()
! Reactions
integer :: n_reaction
integer :: n_reaction ! # of reactions
type(Reaction), pointer :: reactions(:) => null()
end type Nuclide
@ -121,29 +122,29 @@ module cross_section_header
!===============================================================================
type SAB_Table
character(20) :: name
integer :: zaid
real(8) :: awr
real(8) :: kT
character(10) :: name ! name of table, e.g. lwtr.10t
integer :: zaid ! Z and A identifier, e.g. 6012 for Carbon-12
real(8) :: awr ! weight of nucleus in neutron masses
real(8) :: kT ! temperature in MeV (k*T)
! threshold for S(a,b) treatment (usually ~4 eV)
real(8) :: threshold_inelastic
real(8) :: threshold_elastic = 0.0
! Inelastic scattering data
integer :: n_inelastic_e_in
integer :: n_inelastic_e_out
integer :: n_inelastic_mu
integer :: secondary_mode
integer :: n_inelastic_e_in ! # of incoming E for inelastic
integer :: n_inelastic_e_out ! # of outgoing E for inelastic
integer :: n_inelastic_mu ! # of outgoing angles for inelastic
integer :: secondary_mode ! secondary mode (equal/skewed)
real(8), allocatable :: inelastic_e_in(:)
real(8), allocatable :: inelastic_sigma(:)
real(8), allocatable :: inelastic_e_out(:,:)
real(8), allocatable :: inelastic_mu(:,:,:)
! Elastic scattering data
integer :: elastic_mode
integer :: n_elastic_e_in
integer :: n_elastic_mu
integer :: elastic_mode ! elastic mode (discrete/exact)
integer :: n_elastic_e_in ! # of incoming E for elastic
integer :: n_elastic_mu ! # of outgoing angles for elastic
real(8), allocatable :: elastic_e_in(:)
real(8), allocatable :: elastic_P(:)
real(8), allocatable :: elastic_mu(:,:)
@ -154,8 +155,8 @@ module cross_section_header
!===============================================================================
type XsListing
character(10) :: name ! table name, e.g. 92235.70c
character(10) :: alias ! table alias, e.g. U-235.70c
character(12) :: name ! table name, e.g. 92235.70c
character(12) :: alias ! table alias, e.g. U-235.70c
integer :: type ! type of table (cont-E neutron, S(A,b), etc)
integer :: zaid ! ZAID identifier = 1000*Z + A
integer :: filetype ! ASCII or BINARY
@ -174,20 +175,19 @@ module cross_section_header
!===============================================================================
type NuclideMicroXS
integer :: index_grid
integer :: index_temp
integer :: last_index_grid
integer :: last_index_temp
real(8) :: interp_factor
real(8) :: total
real(8) :: elastic
real(8) :: absorption
real(8) :: fission
real(8) :: nu_fission
integer :: index_grid ! index on nuclide energy grid
integer :: index_temp ! temperature index for nuclide
real(8) :: last_E = 0.0 ! last evaluated energy
real(8) :: interp_factor ! interpolation factor on nuc. energy grid
real(8) :: total ! microscropic total xs
real(8) :: elastic ! microscopic elastic scattering xs
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
! Information for S(a,b) use
logical :: use_sab
real(8) :: elastic_sab
logical :: use_sab ! in S(a,b) energy range?
real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
end type NuclideMicroXS
!===============================================================================
@ -196,12 +196,11 @@ module cross_section_header
!===============================================================================
type MaterialMacroXS
real(8) :: total
real(8) :: scatter
real(8) :: elastic
real(8) :: absorption
real(8) :: fission
real(8) :: nu_fission
real(8) :: total ! macroscopic total xs
real(8) :: elastic ! macroscopic elastic scattering xs
real(8) :: absorption ! macroscopic absorption xs
real(8) :: fission ! macroscopic fission xs
real(8) :: nu_fission ! macroscopic production xs
end type MaterialMacroXS
end module cross_section_header

0
src/datatypes.f90 → src/datatypes.F90
View file

0
src/datatypes_header.f90 → src/datatypes_header.F90
View file

28
src/doppler.f90 → src/doppler.F90
View file

@ -210,4 +210,32 @@ contains
end subroutine calculate_F
#ifdef NO_F2008
!===============================================================================
! ERFC computes the complementary error function of x
!===============================================================================
function erfc(x) result(y)
real(8), intent(in) :: x
real(8) :: y
real(8) :: a1 = 0.254829592_8
real(8) :: a2 = -0.284496736_8
real(8) :: a3 = 1.421413741_8
real(8) :: a4 = -1.453152027_8
real(8) :: a5 = 1.061405429_8
real(8) :: p = 0.3275911_8
real(8) :: t
! Abramowitz and Stegun formula 7.1.26
t = 1.0_8/(1.0_8 + p*abs(x))
y = (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x)
! Account for negative values of x
y = sign(y,x)
end function erfc
#endif
end module doppler

0
src/endf.f90 → src/endf.F90
View file

0
src/endf_header.f90 → src/endf_header.F90
View file

98
src/energy_grid.f90 → src/energy_grid.F90
View file

@ -1,8 +1,10 @@
module energy_grid
use constants, only: MAX_LINE_LEN
use datatypes, only: list_insert, list_size, list_delete
use datatypes_header, only: ListReal
use datatypes, only: list_insert, list_size, list_delete, &
dict_create, dict_get_key, dict_has_key, &
dict_add_key, dict_delete
use datatypes_header, only: ListReal, DictionaryCI
use global
use output, only: write_message
@ -18,25 +20,44 @@ contains
subroutine unionized_grid()
type(ListReal), pointer :: list => null()
type(ListReal), pointer :: current => null()
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
integer :: i, j
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: index_list ! index in xs_listings array
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
type(ListReal), pointer :: list => null()
type(ListReal), pointer :: current => null()
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
type(DictionaryCI), pointer :: already_added => null()
message = "Creating unionized energy grid..."
call write_message(5)
! loop over all materials
! Create dictionary for keeping track of cross sections already added
call dict_create(already_added)
! Loop over all files
do i = 1, n_materials
mat => materials(i)
! loop over all isotopes
do j = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(j))
! loop over energy points
call add_grid_points(list, nuc % energy)
do j = 1, mat % n_nuclides
name = mat % names(j)
if (.not. dict_has_key(already_added, name)) then
! loop over energy points
nuc => nuclides(mat % nuclide(j))
call add_grid_points(list, nuc % energy)
! determine name and alias from xs_listings
index_list = dict_get_key(xs_listing_dict, name)
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
! add name and alias to dictionary
call dict_add_key(already_added, name, 0)
call dict_add_key(already_added, alias, 0)
end if
end do
end do
@ -49,8 +70,9 @@ contains
current => current % next
end do
! delete linked list
! delete linked list and dictionary
call list_delete(list)
call dict_delete(already_added)
end subroutine unionized_grid
@ -64,7 +86,7 @@ contains
type(ListReal), pointer :: list
real(8), intent(in) :: energy(:)
integer :: index
integer :: i
integer :: n
real(8) :: E
type(ListReal), pointer :: current => null()
@ -72,7 +94,7 @@ contains
type(ListReal), pointer :: head => null()
type(ListReal), pointer :: tmp => null()
index = 1
i = 1
n = size(energy)
! if the original list is empty, we need to allocate the first element and
@ -80,9 +102,9 @@ contains
if (list_size(list) == 0) then
allocate(list)
current => list
do index = 1, n
current % data = energy(index)
if (index == n) then
do i = 1, n
current % data = energy(i)
if (i == n) then
current % next => null()
return
end if
@ -94,19 +116,19 @@ contains
current => list
head => list
do while (index <= n)
E = energy(index)
do while (i <= n)
E = energy(i)
! If we've reached the end of the grid energy list, add the remaining
! energy points to the end
if (.not. associated(current)) then
! finish remaining energies
do while (index <= n)
do while (i <= n)
allocate(previous % next)
current => previous % next
current % data = energy(index)
current % data = energy(i)
previous => current
index = index + 1
i = i + 1
end do
current%next => null()
exit
@ -127,12 +149,12 @@ contains
nullify(tmp)
! advance index
index = index + 1
i = i + 1
elseif (E == current % data) then
! found the exact same energy, no need to store duplicates so just
! skip and move to next index
index = index + 1
i = i + 1
else
previous => current
current => current % next
@ -152,32 +174,30 @@ contains
subroutine original_indices()
integer :: i, j
integer :: index
integer :: i
integer :: j
integer :: index_e
integer :: n_grid_nuclide
type(Nuclide), pointer :: nuc
real(8) :: union_energy
real(8) :: energy
type(Nuclide), pointer :: nuc
do i = 1, n_nuclides_total
nuc => nuclides(i)
n_grid_nuclide = size(nuc % energy)
allocate(nuc % grid_index(n_grid))
index = 1
energy = nuc % energy(index)
index_e = 1
energy = nuc % energy(index_e)
do j = 1, n_grid
union_energy = e_grid(j)
if (union_energy >= energy .and. index < n_grid_nuclide) then
index = index + 1
energy = nuc % energy(index)
if (union_energy >= energy .and. index_e < n_grid_nuclide) then
index_e = index_e + 1
energy = nuc % energy(index_e)
end if
nuc % grid_index(j) = index-1
nuc % grid_index(j) = index_e - 1
end do
end do
end subroutine original_indices

5
src/error.f90 → src/error.F90
View file

@ -64,9 +64,12 @@ contains
write(eu,*)
end if
! All processors abort
! Release memory from all allocatable arrays
call free_memory()
! Abort program
stop
end subroutine fatal_error
end module error

16
src/fileio.f90 → src/fileio.F90
View file

@ -37,9 +37,9 @@ contains
character(MAX_LINE_LEN) :: line ! single line
character(MAX_WORD_LEN) :: local_words(MAX_WORDS) ! words on one line
integer :: index ! index of words
integer :: index_word ! index of words
index = 0
index_word = 0
do
! read line from file
read(UNIT=unit, FMT='(A100)', IOSTAT=ioError) line
@ -55,17 +55,17 @@ contains
! Check whether there is a continuation line
if (local_words(n) == '&') then
words(index+1:index+n-1) = local_words(1:n-1)
index = index + n - 1
words(index_word+1:index_word+n-1) = local_words(1:n-1)
index_word = index_word + n - 1
else
words(index+1:index+n) = local_words(1:n)
index = index + n
words(index_word+1:index_word+n) = local_words(1:n)
index_word = index_word + n
exit
end if
end do
! set total number of words
n = index
n = index_word
end subroutine get_next_line
@ -79,7 +79,7 @@ contains
integer, intent(in) :: n_lines ! number of lines to skip
integer, intent(out) :: ioError ! error status
integer :: i ! index for number of lines
integer :: i ! index for number of lines
do i = 1, n_lines
read(UNIT=unit, FMT=*, IOSTAT=ioError)

0
src/fission.f90 → src/fission.F90
View file

1144
src/geometry.f90 → src/geometry.F90
View file

File diff suppressed because it is too large Load diff

0
src/geometry_header.f90 → src/geometry_header.F90
View file

26
src/global.f90 → src/global.F90
View file

@ -119,14 +119,22 @@ module global
integer :: rank ! rank of process
logical :: master ! master process?
logical :: mpi_enabled ! is MPI in use and initialized?
integer :: mpi_err ! MPI error code
! ============================================================================
! TIMING VARIABLES
type(Timer) :: time_total ! timer for total run
type(Timer) :: time_init ! timer for initialization
type(Timer) :: time_initialize ! timer for initialization
type(Timer) :: time_read_xs ! timer for reading cross sections
type(Timer) :: time_unionize ! timer for unionizing energy grid
type(Timer) :: time_intercycle ! timer for intercycle synchronization
type(Timer) :: time_ic_tallies ! timer for intercycle accumulate tallies
type(Timer) :: time_ic_sample ! timer for intercycle sampling
type(Timer) :: time_ic_sendrecv ! timer for intercycle SEND/RECV
type(Timer) :: time_ic_rebuild ! timer for intercycle source bank rebuild
type(Timer) :: time_inactive ! timer for inactive cycles
type(Timer) :: time_active ! timer for active cycles
type(Timer) :: time_compute ! timer for computation
! ===========================================================================
@ -148,7 +156,7 @@ module global
! ============================================================================
! MISCELLANEOUS VARIABLES
character(MAX_WORD_LEN) :: path_input ! Path to input file
character(MAX_FILE_LEN) :: path_input ! Path to input file
character(MAX_FILE_LEN) :: path_cross_sections ! Path to cross_sections.xml
! Message used in message/warning/fatal_error
@ -172,6 +180,11 @@ module global
type(Timer) :: time_mat ! timer for mat building
type(Timer) :: time_power ! timer for power iteration
! Trace for single particle
logical :: trace
integer :: trace_cycle
integer(8) :: trace_particle
contains
!===============================================================================
@ -181,10 +194,6 @@ contains
subroutine free_memory()
#ifdef MPI
integer :: ierr
#endif
! Deallocate cells, surfaces, materials
if (allocated(cells)) deallocate(cells)
if (allocated(surfaces)) deallocate(surfaces)
@ -205,12 +214,9 @@ contains
#ifdef MPI
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(ierr)
call MPI_FINALIZE(mpi_err)
#endif
! End program
stop
end subroutine free_memory
end module global

107
src/initialize.f90 → src/initialize.F90
View file

@ -13,15 +13,19 @@ module initialize
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict
use logging, only: create_log
use mcnp_random, only: RN_init_problem
use mpi_routines, only: setup_mpi
use output, only: title, echo_input, message, print_summary, &
print_particle, header, print_plot
use output, only: title, header, print_summary, print_geometry, &
print_plot
use random_lcg, only: initialize_prng
use source, only: initialize_source
use string, only: int_to_str, starts_with, ends_with, lower_case
use tally, only: create_tally_map, TallyObject
use timing, only: timer_start, timer_stop
implicit none
type(DictionaryII), pointer :: build_dict => null()
contains
!===============================================================================
@ -36,7 +40,7 @@ contains
type(Universe), pointer :: univ
! Start initialization timer
call timer_start(time_init)
call timer_start(time_initialize)
! Setup MPI
call setup_mpi()
@ -51,10 +55,8 @@ contains
! Print initialization header block
if (master) call header("INITIALIZATION", 1)
! Initialize random number generator. The first argument corresponds to
! which random number generator to use- in this case one of the L'Ecuyer
! 63-bit RNGs.
call RN_init_problem(3, 0_8, 0_8, 0_8, 0)
! Initialize random number generator
call initialize_prng()
! set up dictionaries
call create_dictionaries()
@ -86,11 +88,15 @@ contains
call normalize_ao()
! Read ACE-format cross sections
call timer_start(time_read_xs)
call read_xs()
call timer_stop(time_read_xs)
! Construct unionized energy grid from cross-sections
call timer_start(time_unionize)
call unionized_grid()
call original_indices()
call timer_stop(time_unionize)
! Create tally map
call create_tally_map()
@ -102,15 +108,15 @@ contains
! stop timer for initialization
if (master) then
if (plotting) then
call print_geometry()
call print_plot()
else
call echo_input()
call print_summary()
end if
end if
! Stop initialization timer
call timer_stop(time_init)
call timer_stop(time_initialize)
end subroutine initialize_run
@ -123,8 +129,8 @@ contains
integer :: i ! loop index
integer :: argc ! number of command line arguments
integer :: last_flag ! index of last flag
character(MAX_LINE_LEN) :: pwd ! present working directory
character(MAX_WORD_LEN) :: argv(10) ! command line arguments
character(MAX_FILE_LEN) :: pwd ! present working directory
character(MAX_WORD_LEN) :: argv(10) ! command line arguments
! Get working directory
call GET_ENVIRONMENT_VARIABLE("PWD", pwd)
@ -148,7 +154,7 @@ contains
end do
! Determine directory where XML input files are
if (argc > 0) then
if (argc > 0 .and. last_flag < argc) then
path_input = argv(last_flag + 1)
! Need to add working directory if the given path is a relative path
if (.not. starts_with(path_input, "/")) then
@ -189,6 +195,9 @@ 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
@ -199,9 +208,9 @@ contains
subroutine prepare_universes()
integer :: i ! index in cells array
integer :: index ! index in universes array
integer :: count ! number of cells in a universe
integer :: i ! index in cells array
integer :: i_univ ! index in universes array
integer :: n_cells_in_univ ! number of cells in a universe
integer, allocatable :: index_cell_in_univ(:) ! the index in the univ%cells
! array for each universe
type(ListKeyValueII), pointer :: key_list => null()
@ -218,19 +227,19 @@ contains
key_list => dict_keys(universe_dict)
do while (associated(key_list))
! find index of universe in universes array
index = key_list%data%value
univ => universes(index)
i_univ = key_list%data%value
univ => universes(i_univ)
univ % id = key_list%data%key
! check for lowest level universe
if (univ % id == 0) BASE_UNIVERSE = index
if (univ % id == 0) BASE_UNIVERSE = i_univ
! find cell count for this universe
count = dict_get_key(cells_in_univ_dict, key_list%data%key)
n_cells_in_univ = dict_get_key(cells_in_univ_dict, key_list%data%key)
! allocate cell list for universe
allocate(univ % cells(count))
univ % n_cells = count
allocate(univ % cells(n_cells_in_univ))
univ % n_cells = n_cells_in_univ
! move to next universe
key_list => key_list%next
@ -246,13 +255,13 @@ contains
c => cells(i)
! get pointer to corresponding universe
index = dict_get_key(universe_dict, c % universe)
univ => universes(index)
i_univ = dict_get_key(universe_dict, c % universe)
univ => universes(i_univ)
! increment the index for the cells array within the Universe object and
! then store the index of the Cell object in that array
index_cell_in_univ(index) = index_cell_in_univ(index) + 1
univ % cells(index_cell_in_univ(index)) = i
index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1
univ % cells(index_cell_in_univ(i_univ)) = i
end do
end subroutine prepare_universes
@ -270,7 +279,7 @@ contains
integer :: i ! index in cells array
integer :: j ! index over surface list
integer :: k
integer :: index ! index in surfaces/materials array
integer :: i_array ! index in surfaces/materials array
integer :: id ! user-specified id
type(Cell), pointer :: c => null()
type(Lattice), pointer :: l => null()
@ -285,8 +294,8 @@ contains
id = c % surfaces(j)
if (id < OP_DIFFERENCE) then
if (dict_has_key(surface_dict, abs(id))) then
index = dict_get_key(surface_dict, abs(id))
c % surfaces(j) = sign(index, id)
i_array = dict_get_key(surface_dict, abs(id))
c % surfaces(j) = sign(i_array, id)
else
message = "Could not find surface " // trim(int_to_str(abs(id))) // &
& " specified on cell " // trim(int_to_str(c % id))
@ -466,46 +475,42 @@ contains
integer, intent(in) :: parent ! cell containing universe
integer, intent(in) :: level ! level of universe
integer :: i ! index for cells in universe
integer :: x,y ! indices for lattice positions
integer :: index ! index in cells array
integer :: i ! index for cells in universe
integer :: x,y ! indices for lattice positions
integer :: i_cell ! index in cells array
integer :: universe_num
type(Cell), pointer :: c => null()
type(Universe), pointer :: subuniverse => null()
type(Lattice), pointer :: lat => null()
type(DictionaryII), pointer :: dict => null()
! set level of the universe
univ % level = level
! loop over all cells in the universe
do i = 1, univ % n_cells
index = univ % cells(i)
c => cells(index)
c%parent = parent
i_cell = univ % cells(i)
c => cells(i_cell)
c % parent = parent
! if this cell is filled with another universe, recursively
! call this subroutine
if (c % type == CELL_FILL) then
subuniverse => universes(c % fill)
call build_universe(subuniverse, index, level + 1)
call build_universe(subuniverse, i_cell, level + 1)
end if
! if this cell is filled by a lattice, need to build the
! universe for each unique lattice element
if (c % type == CELL_LATTICE) then
lat => lattices(c % fill)
call dict_create(dict)
do x = 1, lat % n_x
do y = 1, lat % n_y
lat => lattices(cells(i_cell) % fill)
universe_num = lat % element(x,y)
if (dict_has_key(dict, universe_num)) then
cycle
else
call dict_add_key(dict, universe_num, 0)
if (.not. dict_has_key(build_dict, universe_num)) then
call dict_add_key(build_dict, universe_num, 0)
subuniverse => universes(universe_num)
call build_universe(subuniverse, index, level + 1)
call build_universe(subuniverse, i_cell, level + 1)
end if
end do
end do
@ -521,7 +526,7 @@ contains
subroutine normalize_ao()
integer :: index ! index used for several purposes
integer :: index_list ! index in xs_listings array
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: n ! length of string
@ -530,7 +535,7 @@ contains
real(8) :: x ! atom percent
logical :: percent_in_atom ! nuclides specified in atom percent?
logical :: density_in_atom ! density specified in atom/b-cm?
character(10) :: key ! name of nuclide, e.g. 92235.03c
character(12) :: key ! name of nuclide, e.g. 92235.03c
type(Material), pointer :: mat => null()
! first find the index in the xs_listings array for each nuclide in each
@ -564,8 +569,8 @@ contains
end if
if (dict_has_key(xs_listing_dict, key)) then
index = dict_get_key(xs_listing_dict, key)
mat % xs_listing(j) = index
index_list = dict_get_key(xs_listing_dict, key)
mat % xs_listing(j) = index_list
else
message = "Cannot find cross-section " // trim(key) // &
" in specified cross_sections.xml file."
@ -573,7 +578,7 @@ contains
end if
! determine atomic weight ratio
awr = xs_listings(index) % awr
awr = xs_listings(index_list) % awr
! if given weight percent, convert all values so that they are divided
! by awr. thus, when a sum is done over the values, it's actually
@ -594,8 +599,8 @@ contains
if (.not. density_in_atom) then
sum_percent = ZERO
do j = 1, mat % n_nuclides
index = mat % xs_listing(j)
awr = xs_listings(index) % awr
index_list = mat % xs_listing(j)
awr = xs_listings(index_list) % awr
x = mat % atom_percent(j)
sum_percent = sum_percent + x*awr
end do

50
src/input_xml.f90 → src/input_xml.F90
View file

@ -138,6 +138,12 @@ contains
weight_survive = cutoff_(1) % weight_avg
end if
! Particle trace
if (associated(trace_)) then
trace_cycle = trace_(1)
trace_particle = trace_(2)
end if
end subroutine read_settings_xml
!===============================================================================
@ -153,7 +159,7 @@ contains
integer :: n
integer :: n_x, n_y
integer :: universe_num
integer :: count
integer :: n_cells_in_univ
integer :: coeffs_reqd
logical :: file_exists
character(MAX_LINE_LEN) :: filename
@ -230,12 +236,12 @@ contains
universe_num = cell_(i) % universe
if (.not. dict_has_key(cells_in_univ_dict, universe_num)) then
n_universes = n_universes + 1
count = 1
n_cells_in_univ = 1
call dict_add_key(universe_dict, universe_num, n_universes)
else
count = 1 + dict_get_key(cells_in_univ_dict, universe_num)
n_cells_in_univ = 1 + dict_get_key(cells_in_univ_dict, universe_num)
end if
call dict_add_key(cells_in_univ_dict, universe_num, count)
call dict_add_key(cells_in_univ_dict, universe_num, n_cells_in_univ)
end do
@ -427,7 +433,7 @@ contains
logical :: file_exists
character(3) :: default_xs
character(MAX_WORD_LEN) :: units
character(MAX_WORD_LEN) :: name
character(12) :: name
character(MAX_LINE_LEN) :: filename
type(Material), pointer :: m => null()
type(nuclide_xml), pointer :: nuc => null()
@ -577,7 +583,7 @@ contains
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
integer :: id ! user-specified identifier
integer :: index ! index in meshes array
integer :: i_mesh ! index in meshes array
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
logical :: file_exists ! does tallies.xml file exist?
@ -658,6 +664,7 @@ contains
allocate(m % dimension(n))
allocate(m % origin(n))
allocate(m % width(n))
allocate(m % upper_right(n))
! Read dimensions in each direction
m % dimension = mesh_(i) % dimension
@ -678,6 +685,9 @@ contains
end if
m % width = mesh_(i) % width
! Set upper right coordinate
m % upper_right = m % origin + m % dimension * m % width
! Add mesh to dictionary
call dict_add_key(mesh_dict, m % id, i)
end do
@ -714,7 +724,7 @@ contains
call split_string(tally_(i) % filters % cell, words, n_words)
allocate(t % cell_bins(n_words))
do j = 1, n_words
t % cell_bins(j) % scalar = str_to_int(words(j))
t % cell_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_CELL) = n_words
end if
@ -724,7 +734,7 @@ contains
call split_string(tally_(i) % filters % surface, words, n_words)
allocate(t % surface_bins(n_words))
do j = 1, n_words
t % surface_bins(j) % scalar = str_to_int(words(j))
t % surface_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_SURFACE) = n_words
end if
@ -734,7 +744,7 @@ contains
call split_string(tally_(i) % filters % universe, words, n_words)
allocate(t % universe_bins(n_words))
do j = 1, n_words
t % universe_bins(j) % scalar = str_to_int(words(j))
t % universe_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_UNIVERSE) = n_words
end if
@ -744,7 +754,7 @@ contains
call split_string(tally_(i) % filters % material, words, n_words)
allocate(t % material_bins(n_words))
do j = 1, n_words
t % material_bins(j) % scalar = str_to_int(words(j))
t % material_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_MATERIAL) = n_words
end if
@ -755,8 +765,8 @@ contains
! Determine index in mesh array for this bin
id = t % mesh
if (dict_has_key(mesh_dict, id)) then
index = dict_get_key(mesh_dict, id)
m => meshes(index)
i_mesh = dict_get_key(mesh_dict, id)
m => meshes(i_mesh)
else
message = "Could not find mesh " // trim(int_to_str(id)) // &
" specified on tally " // trim(int_to_str(t % id))
@ -771,7 +781,7 @@ contains
call split_string(tally_(i) % filters % cellborn, words, n_words)
allocate(t % cellborn_bins(n_words))
do j = 1, n_words
t % cellborn_bins(j) % scalar = str_to_int(words(j))
t % cellborn_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_CELLBORN) = n_words
end if
@ -865,11 +875,11 @@ contains
end if
! Check to make sure that only the mesh filter was specified
if (t % mesh == 0 .or. t % n_bins(T_MESH) /= &
product(t % n_bins, t % n_bins > 0)) then
message = "Surface currents must be used with a mesh filter only."
call fatal_error()
end if
!!$ if (t % mesh == 0 .or. t % n_bins(T_MESH) /= &
!!$ product(t % n_bins, t % n_bins > 0)) then
!!$ message = "Surface currents must be used with a mesh filter only."
!!$ call fatal_error()
!!$ end if
! Since the number of bins for the mesh filter was already set
! assuming it was a flux tally, we need to adjust the number of
@ -878,8 +888,8 @@ contains
! Get pointer to mesh
id = t % mesh
index = dict_get_key(mesh_dict, id)
m => meshes(index)
i_mesh = dict_get_key(mesh_dict, id)
m => meshes(i_mesh)
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.

31
src/interpolation.f90 → src/interpolation.F90
View file

@ -5,6 +5,7 @@ module interpolation
use error, only: fatal_error
use global, only: message
use search, only: binary_search
use string, only: int_to_str
implicit none
@ -29,6 +30,7 @@ contains
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: loc_0 ! starting location
integer :: n_regions ! number of interpolation regions
integer :: n_points ! number of tabulated values
@ -49,14 +51,14 @@ contains
end if
! determine number of interpolation regions
n_regions = data(loc_0 + 1)
n_regions = int(data(loc_0 + 1))
! set locations for breakpoints and interpolation schemes
loc_breakpoints = loc_0 + 1
loc_interp = loc_breakpoints + n_regions
! determine number of tabulated values
n_points = data(loc_interp + n_regions + 1)
n_points = int(data(loc_interp + n_regions + 1))
! set locations for x's and y's
loc_x = loc_interp + n_regions + 1
@ -79,10 +81,14 @@ contains
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = data(loc_interp + 1)
interp = int(data(loc_interp + 1))
elseif (n_regions > 1) then
message = "Multiple interpolation regions not yet supported."
call fatal_error()
do j = 1, n_regions
if (i < data(loc_breakpoints + j)) then
interp = int(data(loc_interp + j))
exit
end if
end do
end if
! handle special case of histogram interpolation
@ -111,6 +117,9 @@ contains
case (LOG_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
message = "Unsupported interpolation scheme: " // int_to_str(interp)
call fatal_error()
end select
end function interpolate_tab1_array
@ -129,6 +138,7 @@ contains
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: n_regions ! number of interpolation regions
integer :: n_pairs ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
@ -159,8 +169,12 @@ contains
elseif (n_regions == 1) then
interp = obj % int(1)
elseif (n_regions > 1) then
message = "Multiple interpolation regions not yet supported."
call fatal_error()
do j = 1, n_regions
if (i < obj % nbt(j)) then
interp = obj % int(j)
exit
end if
end do
end if
! handle special case of histogram interpolation
@ -189,6 +203,9 @@ contains
case (LOG_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
message = "Unsupported interpolation scheme: " // int_to_str(interp)
call fatal_error()
end select
end function interpolate_tab1_object

0
src/logging.f90 → src/logging.F90
View file

29
src/main.f90 → src/main.F90
View file

@ -4,17 +4,16 @@ program main
use constants
use global
use initialize, only: initialize_run
use mcnp_random, only: RN_init_particle
use mpi_routines, only: synchronize_bank
use output, only: write_message, header, print_runtime
use particle_header, only: Particle
use plot, only: run_plot
use physics, only: transport
use tally, only: calculate_keff
use random_lcg, only: set_particle_seed
use source, only: get_source_particle
use string, only: int_to_str
use tally, only: synchronize_tallies, write_tallies, &
tally_statistics
tally_statistics, calculate_keff
use timing, only: timer_start, timer_stop
#ifdef MPI
@ -35,10 +34,13 @@ program main
else
call run_problem()
! Calculate statistics for tallies and write to tallies.out
call tally_statistics()
if (master) call write_tallies()
! show timing statistics
call timer_stop(time_total)
if (master) call print_runtime()
if (master) call write_tallies()
end if
! call cmfd calculation
@ -91,7 +93,12 @@ contains
! set random number seed
i_particle = (i_cycle-1)*n_particles + p % id
call RN_init_particle(i_particle)
call set_particle_seed(i_particle)
! set particle trace
trace = .false.
if (i_cycle == trace_cycle .and. &
p % id == trace_particle) trace = .true.
! transport particle
call transport(p)
@ -108,7 +115,11 @@ contains
call timer_start(time_intercycle)
! Collect tallies
if (tallies_on) call synchronize_tallies()
if (tallies_on) then
call timer_start(time_ic_tallies)
call synchronize_tallies()
call timer_stop(time_ic_tallies)
end if
! Distribute fission bank across processors evenly
call synchronize_bank(i_cycle)
@ -122,6 +133,7 @@ contains
if (i_cycle == n_inactive) then
tallies_on = .true.
call timer_stop(time_inactive)
call timer_start(time_active)
end if
! Stop timer for inter-cycle synchronization
@ -129,12 +141,11 @@ contains
end do CYCLE_LOOP
call timer_stop(time_active)
! ==========================================================================
! END OF RUN WRAPUP
! Calculate statistics for tallies
call tally_statistics()
if (master) call header("SIMULATION FINISHED", 1)
end subroutine run_problem

4
src/material_header.f90 → src/material_header.F90
View file

@ -9,7 +9,7 @@ module material_header
type Material
integer :: id ! unique identifier
integer :: n_nuclides ! number of nuclides
character(10), allocatable :: names(:) ! isotope names
character(12), allocatable :: names(:) ! isotope names
integer, allocatable :: xs_listing(:) ! index in xs_listings list
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
@ -18,7 +18,7 @@ module material_header
! S(a,b) data references
logical :: has_sab_table = .false.
character(10) :: sab_name ! name of S(a,b) table
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
end type Material

627
src/mcnp_random.f90
View file

@ -1,627 +0,0 @@
module mcnp_random
!=======================================================================
! Description:
! mcnp_random.F90 -- random number generation routines
!=======================================================================
! This module contains:
!
! * Constants for the RN generator, including initial RN seed for the
! problem & the current RN seed
!
! * MCNP interface routines:
! - random number function: rang()
! - RN initialization for problem: RN_init_problem
! - RN initialization for particle: RN_init_particle
! - get info on RN parameters: RN_query
! - get RN seed for n-th history: RN_query_first
! - set new RN parameters: RN_set
! - skip-ahead in the RN sequence: RN_skip_ahead
! - Unit tests: RN_test_basic, RN_test_skip, RN_test_mixed
!
! * For interfacing with the rest of MCNP, arguments to/from these
! routines will have types of I8 or I4.
! Any args which are to hold random seeds, multipliers,
! skip-distance will be type I8, so that 63 bits can be held without
! truncation.
!
! Revisions:
! * 10-04-2001 - F Brown, initial mcnp version
! * 06-06-2002 - F Brown, mods for extended generators
! * 12-21-2004 - F Brown, added 3 of LeCuyer's 63-bit mult. RNGs
! * 01-29-2005 - J Sweezy, Modify to use mcnp modules prior to automatic
! io unit numbers.
! * 12-02-2005 - F Brown, mods for consistency with C version
!=======================================================================
!-------------------
! MCNP output units
!-------------------
integer, parameter :: iuo = 6
integer, parameter :: jtty = 6
PRIVATE
!---------------------------------------------------
! Kinds for LONG INTEGERS (64-bit) & REAL*8 (64-bit)
!---------------------------------------------------
integer, parameter :: R8 = selected_real_kind(15,307)
integer, parameter :: I8 = selected_int_kind(18)
!-----------------------------------
! Public functions and subroutines for this module
!-----------------------------------
PUBLIC :: rang
PUBLIC :: RN_init_problem
PUBLIC :: RN_init_particle
PUBLIC :: RN_set
PUBLIC :: RN_query
PUBLIC :: RN_query_first
PUBLIC :: RN_update_stats
PUBLIC :: RN_test_basic
PUBLIC :: RN_test_skip
PUBLIC :: RN_test_mixed
PUBLIC :: RN_skip
!-------------------------------------
! Constants for standard RN generators
!-------------------------------------
type :: RN_GEN
integer :: index
integer(I8) :: mult ! generator (multiplier)
integer(I8) :: add ! additive constant
integer :: log2mod ! log2 of modulus, must be <64
integer(I8) :: stride ! stride for particle skip-ahead
integer(I8) :: initseed ! default seed for problem
character(len=8) :: name
end type RN_GEN
! parameters for standard generators
integer, parameter :: n_RN_GEN = 7
type(RN_GEN), SAVE :: standard_generator(n_RN_GEN)
data standard_generator / &
& RN_GEN( 1, 19073486328125_I8, 0_I8, 48, 152917_I8, 19073486328125_I8 , 'mcnp std' ), &
& RN_GEN( 2, 9219741426499971445_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer1' ), &
& RN_GEN( 3, 2806196910506780709_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer2' ), &
& RN_GEN( 4, 3249286849523012805_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer3' ), &
& RN_GEN( 5, 3512401965023503517_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer4' ), &
& RN_GEN( 6, 2444805353187672469_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer5' ), &
& RN_GEN( 7, 1987591058829310733_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer6' ) &
& /
!-----------------------------------------------------------------
! * Linear multiplicative congruential RN algorithm:
!
! RN_SEED = RN_SEED*RN_MULT + RN_ADD mod RN_MOD
!
! * Default values listed below will be used, unless overridden
!-----------------------------------------------------------------
integer, SAVE :: RN_INDEX = 1
integer(I8), SAVE :: RN_MULT = 19073486328125_I8
integer(I8), SAVE :: RN_ADD = 0_I8
integer, SAVE :: RN_BITS = 48
integer(I8), SAVE :: RN_STRIDE = 152917_I8
integer(I8), SAVE :: RN_SEED0 = 19073486328125_I8
integer(I8), SAVE :: RN_MOD = 281474976710656_I8
integer(I8), SAVE :: RN_MASK = 281474976710655_I8
integer(I8), SAVE :: RN_PERIOD = 70368744177664_I8
real(R8), SAVE :: RN_NORM = 1._R8 / 281474976710656._R8
!------------------------------------
! Private data for a single particle
!------------------------------------
integer(I8), save :: RN_SEED = 19073486328125_I8 ! current seed
integer(I8), save :: RN_COUNT = 0_I8 ! current counter
integer(I8), save :: RN_NPS = 0_I8 ! current particle number
!$OMP THREADPRIVATE (RN_SEED, RN_COUNT, RN_NPS)
!------------------------------------------
! Shared data, to collect info on RN usage
!------------------------------------------
integer(I8), SAVE :: RN_COUNT_TOTAL = 0 ! total RN count all particles
integer(I8), SAVE :: RN_COUNT_STRIDE = 0 ! count for stride exceeded
integer(I8), SAVE :: RN_COUNT_MAX = 0 ! max RN count all particles
integer(I8), SAVE :: RN_COUNT_MAX_NPS = 0 ! part index for max count
!---------------------------------------------------------------------
! Reference data: Seeds for case of init.seed = 1,
! Seed numbers for index 1-5, 123456-123460
!---------------------------------------------------------------------
integer(I8), dimension(10,n_RN_GEN) :: RN_CHECK
data RN_CHECK / &
! ***** 1 ***** mcnp standard gen *****
& 19073486328125_I8, 29763723208841_I8, 187205367447973_I8, &
& 131230026111313_I8, 264374031214925_I8, 260251000190209_I8, &
& 106001385730621_I8, 232883458246025_I8, 97934850615973_I8, &
& 163056893025873_I8, &
! ***** 2 *****
& 9219741426499971446_I8, 666764808255707375_I8, 4935109208453540924_I8, &
& 7076815037777023853_I8, 5594070487082964434_I8, 7069484152921594561_I8, &
& 8424485724631982902_I8, 19322398608391599_I8, 8639759691969673212_I8, &
& 8181315819375227437_I8, &
! ***** 3 *****
& 2806196910506780710_I8, 6924308458965941631_I8, 7093833571386932060_I8, &
& 4133560638274335821_I8, 678653069250352930_I8, 6431942287813238977_I8, &
& 4489310252323546086_I8, 2001863356968247359_I8, 966581798125502748_I8, &
& 1984113134431471885_I8, &
! ***** 4 *****
& 3249286849523012806_I8, 4366192626284999775_I8, 4334967208229239068_I8, &
& 6386614828577350285_I8, 6651454004113087106_I8, 2732760390316414145_I8, &
& 2067727651689204870_I8, 2707840203503213343_I8, 6009142246302485212_I8, &
& 6678916955629521741_I8, &
! ***** 5 *****
& 3512401965023503517_I8, 5461769869401032777_I8, 1468184805722937541_I8, &
& 5160872062372652241_I8, 6637647758174943277_I8, 794206257475890433_I8, &
& 4662153896835267997_I8, 6075201270501039433_I8, 889694366662031813_I8, &
& 7299299962545529297_I8, &
! ***** 6 *****
& 2444805353187672469_I8, 316616515307798713_I8, 4805819485453690029_I8, &
& 7073529708596135345_I8, 3727902566206144773_I8, 1142015043749161729_I8, &
& 8632479219692570773_I8, 2795453530630165433_I8, 5678973088636679085_I8, &
& 3491041423396061361_I8, &
! ***** 7 *****
& 1987591058829310733_I8, 5032889449041854121_I8, 4423612208294109589_I8, &
& 3020985922691845009_I8, 5159892747138367837_I8, 8387642107983542529_I8, &
& 8488178996095934477_I8, 708540881389133737_I8, 3643160883363532437_I8, &
& 4752976516470772881_I8 /
!---------------------------------------------------------------------
CONTAINS
!-------------------------------------------------------------------
function rang()
! MCNP random number generator
!
! ***************************************
! ***** modifies RN_SEED & RN_COUNT *****
! ***************************************
implicit none
real(R8) :: rang
RN_SEED = iand( iand( RN_MULT*RN_SEED, RN_MASK) + RN_ADD, RN_MASK)
rang = RN_SEED * RN_NORM
RN_COUNT = RN_COUNT + 1
return
end function rang
!-------------------------------------------------------------------
subroutine RN_skip( skip )
! initialize MCNP random number parameters for particle "nps"
!
! * generate a new particle seed from the base seed
! & particle index
! * set the RN count to zero
implicit none
integer(I8), intent(in) :: skip
RN_SEED = RN_skip_ahead( RN_SEED, skip )
end subroutine RN_skip
!-------------------------------------------------------------------
function RN_skip_ahead( seed, skip )
! advance the seed "skip" RNs: seed*RN_MULT^n mod RN_MOD
implicit none
integer(I8) :: RN_skip_ahead
integer(I8), intent(in) :: seed, skip
integer(I8) :: nskip, gen, g, inc, c, gp, rn, seed_old
seed_old = seed
! add period till nskip>0
nskip = skip
do while( nskip<0_I8 )
if( RN_PERIOD>0_I8 ) then
nskip = nskip + RN_PERIOD
else
nskip = nskip + RN_MASK
nskip = nskip + 1_I8
endif
enddo
! get gen=RN_MULT^n, in log2(n) ops, not n ops !
nskip = iand( nskip, RN_MASK )
gen = 1
g = RN_MULT
inc = 0
c = RN_ADD
do while( nskip>0_I8 )
if( btest(nskip,0) ) then
gen = iand( gen*g, RN_MASK )
inc = iand( inc*g, RN_MASK )
inc = iand( inc+c, RN_MASK )
endif
gp = iand( g+1, RN_MASK )
g = iand( g*g, RN_MASK )
c = iand( gp*c, RN_MASK )
nskip = ishft( nskip, -1 )
enddo
rn = iand( gen*seed_old, RN_MASK )
rn = iand( rn + inc, RN_MASK )
RN_skip_ahead = rn
return
end function RN_skip_ahead
!-------------------------------------------------------------------
subroutine RN_init_problem( new_standard_gen, new_seed, &
& new_stride, new_part1, print_info )
! * initialize MCNP random number parameters for problem,
! based on user input. This routine should be called
! only from the main thread, if OMP threading is being used.
!
! * for initial & continue runs, these args should be set:
! new_standard_gen - index of built-in standard RN generator,
! from RAND gen= (or dbcn(14)
! new_seed - from RAND seed= (or dbcn(1))
! output - logical, print RN seed & mult if true
!
! new_stride - from RAND stride= (or dbcn(13))
! new_part1 - from RAND hist= (or dbcn(8))
!
! * for continue runs only, these should also be set:
! new_count_total - from "rnr" at end of previous run
! new_count_stride - from nrnh(1) at end of previous run
! new_count_max - from nrnh(2) at end of previous run
! new_count_max_nps - from nrnh(3) at end of previous run
!
! * check on size of long-ints & long-int arithmetic
! * check the multiplier
! * advance the base seed for the problem
! * set the initial particle seed
! * initialize the counters for RN stats
implicit none
integer, intent(in) :: new_standard_gen
integer(I8), intent(in) :: new_seed
integer(I8), intent(in) :: new_stride
integer(I8), intent(in) :: new_part1
integer, intent(in) :: print_info
character(len=20) :: printseed
integer(I8) :: itemp1, itemp2, itemp3, itemp4
if( new_standard_gen<1 .or. new_standard_gen>n_RN_GEN ) then
call expire( 'RN_init_problem', &
& ' ***** ERROR: illegal index for built-in RN generator')
endif
! set defaults, override if input supplied: seed, mult, stride
RN_INDEX = new_standard_gen
RN_MULT = standard_generator(RN_INDEX)%mult
RN_ADD = standard_generator(RN_INDEX)%add
RN_STRIDE = standard_generator(RN_INDEX)%stride
RN_SEED0 = standard_generator(RN_INDEX)%initseed
RN_BITS = standard_generator(RN_INDEX)%log2mod
RN_MOD = ishft( 1_I8, RN_BITS )
RN_MASK = ishft( not(0_I8), RN_BITS-64 )
RN_NORM = 2._R8**(-RN_BITS)
if( RN_ADD==0_I8) then
RN_PERIOD = ishft( 1_I8, RN_BITS-2 )
else
RN_PERIOD = ishft( 1_I8, RN_BITS )
endif
if( new_seed>0_I8 ) then
RN_SEED0 = new_seed
endif
if( new_stride>0_I8 ) then
RN_STRIDE = new_stride
endif
RN_COUNT_TOTAL = 0
RN_COUNT_STRIDE = 0
RN_COUNT_MAX = 0
RN_COUNT_MAX_NPS = 0
if( print_info /= 0 ) then
write(printseed,'(i20)') RN_SEED0
write( iuo,1) RN_INDEX, RN_SEED0, RN_MULT, RN_ADD, RN_BITS, RN_STRIDE
write(jtty,2) RN_INDEX, adjustl(printseed)
1 format( &
& /,' ***************************************************', &
& /,' * Random Number Generator = ',i20, ' *', &
& /,' * Random Number Seed = ',i20, ' *', &
& /,' * Random Number Multiplier = ',i20, ' *', &
& /,' * Random Number Adder = ',i20, ' *', &
& /,' * Random Number Bits Used = ',i20, ' *', &
& /,' * Random Number Stride = ',i20, ' *', &
& /,' ***************************************************',/)
2 format(' comment. using random number generator ',i2, &
& ', initial seed = ',a20)
endif
! double-check on number of bits in a long int
if( bit_size(RN_SEED)<64 ) then
call expire( 'RN_init_problem', &
& ' ***** ERROR: <64 bits in long-int, can-t generate RN-s')
endif
itemp1 = 5_I8**25
itemp2 = 5_I8**19
itemp3 = ishft(2_I8**62-1_I8,1) + 1_I8
itemp4 = itemp1*itemp2
if( iand(itemp4,itemp3)/=8443747864978395601_I8 ) then
call expire( 'RN_init_problem', &
& ' ***** ERROR: can-t do 64-bit integer ops for RN-s')
endif
if( new_part1>1_I8 ) then
! advance the problem seed to that for part1
RN_SEED0 = RN_skip_ahead( RN_SEED0, (new_part1-1_I8)*RN_STRIDE )
itemp1 = RN_skip_ahead( RN_SEED0, RN_STRIDE )
if( print_info /= 0 ) then
write(printseed,'(i20)') itemp1
write( iuo,3) new_part1, RN_SEED0, itemp1
write(jtty,4) new_part1, adjustl(printseed)
3 format( &
& /,' ***************************************************', &
& /,' * Random Number Seed will be advanced to that for *', &
& /,' * previous particle number = ',i20, ' *', &
& /,' * New RN Seed for problem = ',i20, ' *', &
& /,' * Next Random Number Seed = ',i20, ' *', &
& /,' ***************************************************',/)
4 format(' comment. advancing random number to particle ',i12, &
& ', initial seed = ',a20)
endif
endif
! set the initial particle seed
RN_SEED = RN_SEED0
RN_COUNT = 0
RN_NPS = 0
return
end subroutine RN_init_problem
!-------------------------------------------------------------------
subroutine RN_init_particle( nps )
! initialize MCNP random number parameters for particle "nps"
!
! * generate a new particle seed from the base seed
! & particle index
! * set the RN count to zero
implicit none
integer(I8), intent(in) :: nps
RN_SEED = RN_skip_ahead( RN_SEED0, nps*RN_STRIDE )
RN_COUNT = 0
RN_NPS = nps
return
end subroutine RN_init_particle
!-------------------------------------------------------------------
subroutine RN_set( key, value )
implicit none
character(len=*), intent(in) :: key
integer(I8), intent(in) :: value
character(len=20) :: printseed
integer(I8) :: itemp1
if( key == "stride" ) then
if( value>0_I8 ) then
RN_STRIDE = value
endif
endif
if( key == "count_total" ) RN_COUNT_TOTAL = value
if( key == "count_stride" ) RN_COUNT_STRIDE = value
if( key == "count_max" ) RN_COUNT_MAX = value
if( key == "count_max_nps" ) RN_COUNT_MAX_NPS = value
if( key == "seed" ) then
if( value>0_I8 ) then
RN_SEED0 = value
RN_SEED = RN_SEED0
RN_COUNT = 0
RN_NPS = 0
endif
endif
if( key == "part1" ) then
if( value>1_I8 ) then
! advance the problem seed to that for part1
RN_SEED0 = RN_skip_ahead( RN_SEED0, (value-1_I8)*RN_STRIDE )
itemp1 = RN_skip_ahead( RN_SEED0, RN_STRIDE )
write(printseed,'(i20)') itemp1
write( iuo,3) value, RN_SEED0, itemp1
write(jtty,4) value, adjustl(printseed)
3 format( &
& /,' ***************************************************', &
& /,' * Random Number Seed will be advanced to that for *', &
& /,' * previous particle number = ',i20, ' *', &
& /,' * New RN Seed for problem = ',i20, ' *', &
& /,' * Next Random Number Seed = ',i20, ' *', &
& /,' ***************************************************',/)
4 format(' comment. advancing random number to particle ',i12, &
& ', initial seed = ',a20)
RN_SEED = RN_SEED0
RN_COUNT = 0
RN_NPS = 0
endif
endif
return
end subroutine RN_set
!-------------------------------------------------------------------
function RN_query( key )
implicit none
integer(I8) :: RN_query
character(len=*), intent(in) :: key
RN_query = 0_I8
if( key == "seed" ) RN_query = RN_SEED
if( key == "stride" ) RN_query = RN_STRIDE
if( key == "mult" ) RN_query = RN_MULT
if( key == "add" ) RN_query = RN_ADD
if( key == "count" ) RN_query = RN_COUNT
if( key == "period" ) RN_query = RN_PERIOD
if( key == "count_total" ) RN_query = RN_COUNT_TOTAL
if( key == "count_stride" ) RN_query = RN_COUNT_STRIDE
if( key == "count_max" ) RN_query = RN_COUNT_MAX
if( key == "count_max_nps" ) RN_query = RN_COUNT_MAX_NPS
if( key == "first" ) RN_query = RN_SEED0
return
end function RN_query
!-------------------------------------------------------------------
function RN_query_first( nps )
implicit none
integer(I8) :: RN_query_first
integer(I8), intent(in) :: nps
RN_query_first = RN_skip_ahead( RN_SEED0, nps*RN_STRIDE )
return
end function RN_query_first
!-------------------------------------------------------------------
subroutine RN_update_stats()
! update overall RN count info
implicit none
!$OMP CRITICAL (RN_STATS)
RN_COUNT_TOTAL = RN_COUNT_TOTAL + RN_COUNT
if( RN_COUNT>RN_COUNT_MAX ) then
RN_COUNT_MAX = RN_COUNT
RN_COUNT_MAX_NPS = RN_NPS
endif
if( RN_COUNT>RN_STRIDE ) then
RN_COUNT_STRIDE = RN_COUNT_STRIDE + 1
endif
!$OMP END CRITICAL (RN_STATS)
RN_COUNT = 0
RN_NPS = 0
return
end subroutine RN_update_stats
!-------------------------------------------------------------------
subroutine expire( c1, c2 )
character(len=*), intent(in) :: c1, c2
write(*,*) ' ********** error: ',c1
write(*,*) ' ********** error: ',c2
stop '**error**'
end subroutine expire
!-------------------------------------------------------------------
!###################################################################
!#
!# Unit tests
!#
!###################################################################
subroutine RN_test_basic( new_gen )
! test routine for basic random number generator
implicit none
integer, intent(in) :: new_gen
real(R8) :: s
integer(I8) :: seeds(10)
integer :: i, j
write(jtty,"(/,a)") " ***** random number - basic test *****"
! set the seed
call RN_init_problem( new_gen, 1_I8, 0_I8, 0_I8, 1 )
! get the first 5 seeds, then skip a few, get 5 more - directly
s = 0.0_R8
do i = 1,5
s = s + rang()
seeds(i) = RN_query( "seed" )
enddo
do i = 6,123455
s = s + rang()
enddo
do i = 6,10
s = s + rang()
seeds(i) = RN_query( "seed" )
enddo
! compare
do i = 1,10
j = i
if( i>5 ) j = i + 123450
write(jtty,"(1x,i6,a,i20,a,i20)") &
& j, " reference: ", RN_CHECK(i,new_gen), " computed: ", seeds(i)
if( seeds(i)/=RN_CHECK(i,new_gen) ) then
write(jtty,"(a)") " ***** basic_test of RN generator failed:"
endif
enddo
return
end subroutine RN_test_basic
!-------------------------------------------------------------------
subroutine RN_test_skip( new_gen )
! test routine for basic random number generation & skip-ahead
implicit none
integer, intent(in) :: new_gen
integer(I8) :: seeds(10)
integer :: i, j
! set the seed
call RN_init_problem( new_gen, 1_I8, 0_I8, 0_I8, 0 )
! use the skip-ahead function to get first 5 seeds, then 5 more
do i = 1,10
j = i
if( i>5 ) j = i + 123450
seeds(i) = RN_skip_ahead( 1_I8, int(j,I8) )
enddo
! compare
write(jtty,"(/,a)") " ***** random number - skip test *****"
do i = 1,10
j = i
if( i>5 ) j = i + 123450
write(jtty,"(1x,i6,a,i20,a,i20)") &
& j, " reference: ", RN_CHECK(i,new_gen), " computed: ", seeds(i)
if( seeds(i)/=RN_CHECK(i,new_gen) ) then
write(jtty,"(a)") " ***** skip_test of RN generator failed:"
endif
enddo
return
end subroutine RN_test_skip
!-------------------------------------------------------------------
subroutine RN_test_mixed( new_gen )
! test routine -- print RN's 1-5 & 123456-123460,
! with reference vals
implicit none
integer, intent(in) :: new_gen
integer(I8) :: r
integer :: i, j
write(jtty,"(/,a)") " ***** random number - mixed test *****"
! set the seed & set the stride to 1
call RN_init_problem( new_gen, 1_I8, 1_I8, 0_I8, 0 )
write(jtty,"(a,i20,z20)") " RN_MULT = ", RN_MULT, RN_MULT
write(jtty,"(a,i20,z20)") " RN_ADD = ", RN_ADD, RN_ADD
write(jtty,"(a,i20,z20)") " RN_MOD = ", RN_MOD, RN_MOD
write(jtty,"(a,i20,z20)") " RN_MASK = ", RN_MASK, RN_MASK
write(jtty,"(a,i20)") " RN_BITS = ", RN_BITS
write(jtty,"(a,i20)") " RN_PERIOD = ", RN_PERIOD
write(jtty,"(a,es20.13)") " RN_NORM = ", RN_NORM
write(jtty,"(a)") " "
do i = 1,10
j = i
if( i>5 ) j = i + 123450
call RN_init_particle( int(j,I8) )
r = RN_query( "seed" )
write(jtty,"(1x,i6,a,i20,a,i20)") &
& j, " reference: ", RN_CHECK(i,new_gen)," computed: ", r
if( r/=RN_CHECK(i,new_gen) ) then
write(jtty,"(a)") " ***** mixed test of RN generator failed:"
endif
enddo
return
end subroutine RN_test_mixed
!-------------------------------------------------------------------
end module mcnp_random

46
src/mesh.f90 → src/mesh.F90
View file

@ -12,30 +12,52 @@ contains
! GET_MESH_BIN determines the tally bin for a particle in a structured mesh
!===============================================================================
subroutine get_mesh_bin(m, xyz, bin, in_mesh)
subroutine get_mesh_bin(m, xyz, bin)
type(StructuredMesh), pointer :: m
real(8), intent(in) :: xyz(:)
integer, intent(out) :: bin
logical, intent(out) :: in_mesh
integer :: n
integer, allocatable :: ijk(:)
integer :: n
integer :: ijk(3)
logical :: in_mesh
! Get number of dimensions
n = m % n_dimension
! Create indices array same size as xyz
allocate(ijk(n))
! Check for cases where particle is outside of mesh
if (xyz(1) < m % origin(1)) then
bin = NO_BIN_FOUND
return
elseif (xyz(1) > m % upper_right(1)) then
bin = NO_BIN_FOUND
return
elseif (xyz(2) < m % origin(2)) then
bin = NO_BIN_FOUND
return
elseif (xyz(2) > m % upper_right(2)) then
bin = NO_BIN_FOUND
return
end if
if (n > 2) then
if (xyz(3) < m % origin(3)) then
bin = NO_BIN_FOUND
return
elseif (xyz(3) > m % upper_right(3)) then
bin = NO_BIN_FOUND
return
end if
end if
! Determine indices
call get_mesh_indices(m, xyz(1:n), ijk, in_mesh)
call get_mesh_indices(m, xyz(1:n), ijk(1:n), in_mesh)
! Convert indices to bin
bin = mesh_indices_to_bin(m, ijk)
! Release memory for ijk
deallocate(ijk)
if (in_mesh) then
bin = mesh_indices_to_bin(m, ijk)
else
bin = NO_BIN_FOUND
end if
end subroutine get_mesh_bin

5
src/mesh_header.f90 → src/mesh_header.F90
View file

@ -12,8 +12,9 @@ module mesh_header
integer :: type
integer :: n_dimension
integer, allocatable :: dimension(:)
integer, allocatable :: origin(:)
integer, allocatable :: width(:)
real(8), allocatable :: origin(:)
real(8), allocatable :: upper_right(:)
real(8), allocatable :: width(:)
end type StructuredMesh
end module mesh_header

183
src/mpi_routines.f90 → src/mpi_routines.F90
View file

@ -3,10 +3,11 @@ module mpi_routines
use constants, only: MAX_LINE_LEN
use error, only: fatal_error
use global
use mcnp_random, only: rang, RN_init_particle, RN_skip
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use random_lcg, only: prn, set_particle_seed, prn_skip
use tally_header, only: TallyObject
use timing, only: timer_start, timer_stop
#ifdef MPI
use mpi
@ -16,7 +17,6 @@ module mpi_routines
integer :: MPI_BANK ! MPI datatype for fission bank
integer(8) :: bank_index ! Fission bank site unique identifier
real(8) :: t_sync(4) ! synchronization time
contains
@ -30,7 +30,6 @@ contains
#ifdef MPI
integer :: i
integer :: ierr ! Error status
integer :: bank_blocks(4) ! Count for each datatype
integer :: bank_types(4) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements
@ -40,22 +39,22 @@ contains
mpi_enabled = .true.
! Initialize MPI
call MPI_INIT(ierr)
if (ierr /= MPI_SUCCESS) then
call MPI_INIT(mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Failed to initialize MPI."
call fatal_error()
end if
! Determine number of processors
call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, ierr)
if (ierr /= MPI_SUCCESS) then
call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Could not determine number of processors."
call fatal_error()
end if
! Determine rank of each processor
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
if (ierr /= MPI_SUCCESS) then
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Could not determine MPI rank."
call fatal_error()
end if
@ -68,10 +67,10 @@ contains
end if
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % id, bank_disp(1), ierr)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), ierr)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), ierr)
call MPI_GET_ADDRESS(b % E, bank_disp(4), ierr)
call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err)
! Adjust displacements
base = bank_disp(1)
@ -83,10 +82,9 @@ contains
bank_blocks = (/ 1, 3, 3, 1 /)
bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /)
call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, &
& bank_types, MPI_BANK, ierr)
call MPI_TYPE_COMMIT(MPI_BANK, ierr)
& bank_types, MPI_BANK, mpi_err)
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
t_sync = ZERO
#else
! if no MPI, set processor to master
mpi_enabled = .false.
@ -111,8 +109,7 @@ contains
integer(8) :: start ! starting index in local fission bank
integer(8) :: finish ! ending index in local fission bank
integer(8) :: total ! total sites in global fission bank
integer(8) :: count ! index for source bank
integer(8) :: index ! index for id -- accounts for all nodes
integer(8) :: index_local ! index for source bank
integer :: send_to_left ! # of bank sites to send/recv to or from left
integer :: send_to_right ! # of bank sites to send/recv to or from right
integer(8) :: sites_needed ! # of sites to be sampled
@ -123,32 +120,25 @@ contains
& right_bank(:) ! bank sites to send/recv to or fram right node
#ifdef MPI
integer :: ierr
integer :: status(MPI_STATUS_SIZE) ! message status
integer :: request ! communication request for sending sites
integer :: request_left ! communication request for recv sites from left
integer :: request_right ! communication request for recv sites from right
real(8) :: t0, t1, t2, t3, t4
#endif
message = "Collecting number of fission sites..."
call write_message(8)
#ifdef MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
t0 = MPI_WTIME()
! Determine starting index for fission bank and total sites in fission bank
start = 0_8
call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, &
& MPI_COMM_WORLD, ierr)
& MPI_COMM_WORLD, mpi_err)
finish = start + n_bank
total = finish
call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, &
& MPI_COMM_WORLD, ierr)
& MPI_COMM_WORLD, mpi_err)
t1 = MPI_WTIME()
t_sync(1) = t_sync(1) + (t1 - t0)
#else
start = 0_8
finish = n_bank
@ -162,14 +152,13 @@ contains
end if
! Make sure all processors start at the same point for random sampling
call RN_init_particle(int(i_cycle,8))
call set_particle_seed(int(i_cycle,8))
! Skip ahead however many random numbers are needed
call RN_skip(start)
call prn_skip(start)
allocate(temp_sites(2*work))
count = 0_8 ! Index for local source_bank
index = 0_8 ! Index for global source id -- must account for all nodes
index_local = 0_8 ! Index for local source_bank
if (total < n_particles) then
sites_needed = mod(n_particles,total)
@ -181,9 +170,11 @@ contains
message = "Sampling fission sites..."
call write_message(8)
call timer_start(time_ic_sample)
! ==========================================================================
! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
do i = 1, n_bank
do i = 1, int(n_bank,4)
! If there are less than n_particles particles banked, automatically add
! int(n_particles/total) sites to temp_sites. For example, if you need
@ -192,15 +183,15 @@ contains
if (total < n_particles) then
do j = 1,int(n_particles/total)
! If index is within this node's range, add site to source
count = count + 1
temp_sites(count) = fission_bank(i)
index_local = index_local + 1
temp_sites(index_local) = fission_bank(i)
end do
end if
! Randomly sample sites needed
if (rang() < p_sample) then
count = count + 1
temp_sites(count) = fission_bank(i)
if (prn() < p_sample) then
index_local = index_local + 1
temp_sites(index_local) = fission_bank(i)
end if
end do
@ -208,37 +199,44 @@ contains
! the source bank
#ifdef MPI
start = 0_8
call MPI_EXSCAN(count, start, 1, MPI_INTEGER8, MPI_SUM, &
& MPI_COMM_WORLD, ierr)
finish = start + count
call MPI_EXSCAN(index_local, start, 1, MPI_INTEGER8, MPI_SUM, &
& MPI_COMM_WORLD, mpi_err)
finish = start + index_local
total = finish
call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, &
& MPI_COMM_WORLD, ierr)
& MPI_COMM_WORLD, mpi_err)
#else
start = 0_8
finish = count
total = count
finish = index_local
total = index_local
#endif
! Determine how many sites to send to adjacent nodes
send_to_left = bank_first - 1 - start
send_to_right = finish - bank_last
send_to_left = int(bank_first - 1_8 - start, 4)
send_to_right = int(finish - bank_last, 4)
! Check to make sure number of sites is not more than size of bank
if (abs(send_to_left) > work .or. abs(send_to_right) > work) then
message = "Tried sending sites to neighboring process greater than " &
// "the size of the source bank."
call fatal_error()
end if
if (rank == n_procs - 1) then
if (total > n_particles) then
! If we have extra sites sampled, we will simply discard the extra
! ones on the last processor
if (rank == n_procs - 1) then
count = count - send_to_right
index_local = index_local - send_to_right
end if
elseif (total < n_particles) then
! If we have too few sites, grab sites from the very end of the
! fission bank
sites_needed = n_particles - total
do i = 1, sites_needed
count = count + 1
temp_sites(count) = fission_bank(n_bank - sites_needed + i)
do i = 1, int(sites_needed,4)
index_local = index_local + 1
temp_sites(index_local) = fission_bank(n_bank - sites_needed + i)
end do
end if
@ -246,10 +244,10 @@ contains
send_to_right = 0
end if
#ifdef MPI
t2 = MPI_WTIME()
t_sync(2) = t_sync(2) + (t2 - t1)
call timer_stop(time_ic_sample)
call timer_start(time_ic_sendrecv)
#ifdef MPI
message = "Sending fission sites..."
call write_message(8)
@ -259,73 +257,72 @@ contains
allocate(right_bank(abs(send_to_right)))
if (send_to_right > 0) then
i = count - send_to_right + 1
i = index_local - send_to_right + 1
call MPI_ISEND(temp_sites(i), send_to_right, MPI_BANK, rank+1, 0, &
& MPI_COMM_WORLD, request, ierr)
& MPI_COMM_WORLD, request, mpi_err)
else if (send_to_right < 0) then
call MPI_IRECV(right_bank, -send_to_right, MPI_BANK, rank+1, 1, &
& MPI_COMM_WORLD, request_right, ierr)
& MPI_COMM_WORLD, request_right, mpi_err)
end if
if (send_to_left < 0) then
call MPI_IRECV(left_bank, -send_to_left, MPI_BANK, rank-1, 0, &
& MPI_COMM_WORLD, request_left, ierr)
& MPI_COMM_WORLD, request_left, mpi_err)
else if (send_to_left > 0) then
call MPI_ISEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, &
& MPI_COMM_WORLD, request, ierr)
& MPI_COMM_WORLD, request, mpi_err)
end if
t3 = MPI_WTIME()
t_sync(3) = t_sync(3) + (t3 - t2)
#endif
call timer_stop(time_ic_sendrecv)
call timer_start(time_ic_rebuild)
message = "Constructing source bank..."
call write_message(8)
! ==========================================================================
! RECONSTRUCT SOURCE BANK
if (send_to_left < 0 .and. send_to_right >= 0) then
i = -send_to_left ! size of first block
j = count - send_to_right ! size of second block
i = -send_to_left ! size of first block
j = int(index_local,4) - send_to_right ! size of second block
call copy_from_bank(temp_sites, i+1, j)
#ifdef MPI
call MPI_WAIT(request_left, status, ierr)
call MPI_WAIT(request_left, status, mpi_err)
#endif
call copy_from_bank(left_bank, 1, i)
else if (send_to_left >= 0 .and. send_to_right < 0) then
i = count - send_to_left ! size of first block
j = -send_to_right ! size of second block
i = int(index_local,4) - send_to_left ! size of first block
j = -send_to_right ! size of second block
call copy_from_bank(temp_sites(1+send_to_left), 1, i)
#ifdef MPI
call MPI_WAIT(request_right, status, ierr)
call MPI_WAIT(request_right, status, mpi_err)
#endif
call copy_from_bank(right_bank, i+1, j)
else if (send_to_left >= 0 .and. send_to_right >= 0) then
i = count - send_to_left - send_to_right
i = int(index_local,4) - send_to_left - send_to_right
call copy_from_bank(temp_sites(1+send_to_left), 1, i)
else if (send_to_left < 0 .and. send_to_right < 0) then
i = -send_to_left
j = count
j = int(index_local,4)
k = -send_to_right
call copy_from_bank(temp_sites, i+1, j)
#ifdef MPI
call MPI_WAIT(request_left, status, ierr)
call MPI_WAIT(request_left, status, mpi_err)
#endif
call copy_from_bank(left_bank, 1, i)
#ifdef MPI
call MPI_WAIT(request_right, status, ierr)
call MPI_WAIT(request_right, status, mpi_err)
#endif
call copy_from_bank(right_bank, i+j+1, k)
end if
! Reset source index
source_index = 0_8
call timer_stop(time_ic_rebuild)
#ifdef MPI
t4 = MPI_WTIME()
t_sync(4) = t_sync(4) + (t4 - t3)
deallocate(left_bank )
deallocate(left_bank)
deallocate(right_bank)
#endif
deallocate(temp_sites)
@ -336,30 +333,29 @@ contains
! COPY_FROM_BANK
!===============================================================================
subroutine copy_from_bank(temp_bank, index, n_sites)
subroutine copy_from_bank(temp_bank, i_start, n_sites)
integer, intent(in) :: n_sites ! # of bank sites to copy
type(Bank), intent(in) :: temp_bank(n_sites)
integer, intent(in) :: index ! starting index in source_bank
integer, intent(in) :: i_start ! starting index in source_bank
integer :: i ! index in temp_bank
integer :: index_source ! index in source_bank
integer :: i ! index in temp_bank
integer :: i_source ! index in source_bank
type(Particle), pointer :: p
do i = 1, n_sites
index_source = index + i - 1
p => source_bank(index_source)
p % xyz = temp_bank(i) % xyz
p % xyz_local = temp_bank(i) % xyz
p % last_xyz = temp_bank(i) % xyz
p % uvw = temp_bank(i) % uvw
p % E = temp_bank(i) % E
p % last_E = p % E
i_source = i_start + i - 1
p => source_bank(i_source)
! set defaults
call initialize_particle(p)
p % coord % xyz = temp_bank(i) % xyz
p % coord % uvw = temp_bank(i) % uvw
p % last_xyz = temp_bank(i) % xyz
p % E = temp_bank(i) % E
p % last_E = temp_bank(i) % E
end do
end subroutine copy_from_bank
@ -374,8 +370,7 @@ contains
integer :: i
integer :: n
integer :: m
integer :: count
integer :: ierr
integer :: n_bins
real(8), allocatable :: tally_temp(:,:)
type(TallyObject), pointer :: t
@ -384,7 +379,7 @@ contains
n = t % n_total_bins
m = t % n_macro_bins
count = n*m
n_bins = n*m
allocate(tally_temp(n,m))
@ -392,15 +387,15 @@ contains
if (master) then
! Description of MPI_IN_PLANE
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, count, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, ierr)
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Transfer values to val_history on master
t % scores(:,:) % val_history = tally_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, tally_temp, count, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, ierr)
call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Reset val_history on other processors
t % scores(:,:) % val_history = 0

241
src/output.f90 → src/output.F90
View file

@ -9,6 +9,7 @@ module output
use geometry_header, only: Cell, Universe, Surface
use global
use mesh_header, only: StructuredMesh
use particle_header, only: Particle, LocalCoord
use string, only: upper_case, int_to_str, real_to_str
use tally_header, only: TallyObject
@ -27,8 +28,8 @@ contains
subroutine title()
character(10) :: date
character(8) :: time
character(10) :: today_date
character(8) :: today_time
write(ou,*)
write(ou,*) ' .d88888b. 888b d888 .d8888b.'
@ -51,42 +52,13 @@ contains
100 format (6X,"Version:",9X,I1,".",I1,".",I1)
! Write the date and time
call get_today(date, time)
write(ou,101) trim(date), trim(time)
call get_today(today_date, today_time)
write(ou,101) trim(today_date), trim(today_time)
101 format (6X,"Date/Time:",7X,A,1X,A)
write(ou,*)
end subroutine title
!===============================================================================
! ECHO_INPUT displays summary information about the problem about to be run
! after reading all the input.
!===============================================================================
subroutine echo_input()
! Display problem summary
call header("PROBLEM SUMMARY")
if (problem_type == PROB_CRITICALITY) then
write(ou,100) 'Problem type:', 'Criticality'
write(ou,100) 'Number of Cycles:', int_to_str(n_cycles)
write(ou,100) 'Number of Inactive Cycles:', int_to_str(n_inactive)
elseif (problem_type == PROB_SOURCE) then
write(ou,100) 'Problem type:', 'External Source'
end if
write(ou,100) 'Number of Particles:', int_to_str(n_particles)
! Display geometry summary
call header("GEOMETRY SUMMARY")
write(ou,100) 'Number of Cells:', int_to_str(n_cells)
write(ou,100) 'Number of Surfaces:', int_to_str(n_surfaces)
write(ou,100) 'Number of Materials:', int_to_str(n_materials)
! Format descriptor for columns
100 format (1X,A,T35,A)
end subroutine echo_input
!===============================================================================
! HEADER displays a header block according to a specified level. If no level is
! specified, it is assumed to be a minor header block (H3).
@ -186,11 +158,11 @@ contains
character(8), intent(out) :: today_time
integer :: val(8)
character(8) :: date
character(10) :: time
character(8) :: date_
character(10) :: time_
character(5) :: zone
call date_and_time(date, time, zone, val)
call date_and_time(date_, time_, zone, val)
! val(1) = year (YYYY)
! val(2) = month (MM)
! val(3) = day (DD)
@ -202,18 +174,18 @@ contains
if (val(2) < 10) then
if (val(3) < 10) then
today_date = date(6:6) // "/" // date(8:8) // "/" // date(1:4)
today_date = date_(6:6) // "/" // date_(8:8) // "/" // date_(1:4)
else
today_date = date(6:6) // "/" // date(7:8) // "/" // date(1:4)
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)
today_date = date_(5:6) // "/" // date_(8:8) // "/" // date_(1:4)
else
today_date = date(5:6) // "/" // date(7:8) // "/" // date(1:4)
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)
today_time = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
end subroutine get_today
@ -223,12 +195,16 @@ contains
subroutine print_particle(p)
type(Particle), pointer :: p
type(Particle), pointer :: p
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
integer :: i
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(LocalCoord), pointer :: coord => null()
! display type of particle
select case (p % type)
case (NEUTRON)
write(ou,*) 'Neutron ' // int_to_str(p % id)
@ -239,44 +215,55 @@ contains
case default
write(ou,*) 'Unknown Particle ' // int_to_str(p % id)
end select
write(ou,*) ' x = ' // real_to_str(p % xyz(1))
write(ou,*) ' y = ' // real_to_str(p % xyz(2))
write(ou,*) ' z = ' // real_to_str(p % xyz(3))
write(ou,*) ' x local = ' // real_to_str(p % xyz_local(1))
write(ou,*) ' y local = ' // real_to_str(p % xyz_local(2))
write(ou,*) ' z local = ' // real_to_str(p % xyz_local(3))
write(ou,*) ' u = ' // real_to_str(p % uvw(1))
write(ou,*) ' v = ' // real_to_str(p % uvw(2))
write(ou,*) ' w = ' // real_to_str(p % uvw(3))
write(ou,*) ' Weight = ' // real_to_str(p % wgt)
write(ou,*) ' Energy = ' // real_to_str(p % E)
write(ou,*) ' x index = ' // int_to_str(p % index_x)
write(ou,*) ' y index = ' // int_to_str(p % index_y)
write(ou,*) ' IE = ' // int_to_str(p % IE)
write(ou,*) ' Interpolation factor = ' // real_to_str(p % interp)
if (p % cell > 0) then
c => cells(p % cell)
write(ou,*) ' Cell = ' // int_to_str(c % id)
else
write(ou,*) ' Cell not determined'
end if
! loop through each level of universes
coord => p % coord0
i = 0
do while(associated(coord))
! Print level
write(ou,*) ' Level ' // trim(int_to_str(i))
if (p % surface > 0) then
! Print cell for this level
if (coord % cell /= NONE) then
c => cells(coord % cell)
write(ou,*) ' Cell = ' // trim(int_to_str(c % id))
end if
! Print universe for this level
if (coord % universe /= NONE) then
u => universes(coord % universe)
write(ou,*) ' Universe = ' // trim(int_to_str(u % id))
end if
! Print information on lattice
if (coord % lattice /= NONE) then
l => lattices(coord % lattice)
write(ou,*) ' Lattice = ' // trim(int_to_str(l % id))
write(ou,*) ' Lattice position = (' // trim(int_to_str(&
p % coord % lattice_x)) // ',' // trim(int_to_str(&
p % coord % lattice_y)) // ')'
end if
! Print local coordinates
write(ou,'(1X,A,3ES11.4)') ' xyz = ', coord % xyz
write(ou,'(1X,A,3ES11.4)') ' uvw = ', coord % uvw
coord => coord % next
i = i + 1
end do
! Print surface
if (p % surface /= NONE) then
s => surfaces(p % surface)
write(ou,*) ' Surface = ' // int_to_str(s % id)
else
write(ou,*) ' Surface = None'
end if
u => universes(p % universe)
write(ou,*) ' Universe = ' // int_to_str(u % id)
write(ou,*) ' Weight = ' // real_to_str(p % wgt)
write(ou,*) ' Energy = ' // real_to_str(p % E)
write(ou,*) ' IE = ' // int_to_str(p % IE)
write(ou,*) ' Interpolation factor = ' // real_to_str(p % interp)
write(ou,*)
nullify(c)
nullify(s)
nullify(u)
end subroutine print_particle
!===============================================================================
@ -354,10 +341,6 @@ contains
write(ou,*) ' Surface Specification:' // trim(string)
write(ou,*)
! nullify associated pointers
nullify(u)
nullify(m)
end subroutine print_cell
!===============================================================================
@ -382,8 +365,6 @@ contains
write(ou,*) ' Cells =' // trim(string)
write(ou,*)
nullify(c)
end subroutine print_universe
!===============================================================================
@ -514,8 +495,6 @@ contains
end if
write(ou,*)
nullify(nuc)
end subroutine print_material
!===============================================================================
@ -641,20 +620,17 @@ contains
end subroutine print_tally
!===============================================================================
! PRINT_SUMMARY displays the attributes of all cells, universes,
! surfaces and materials read in the input file. Very useful for
! debugging!
! PRINT_GEOMETRY displays the attributes of all cells, surfaces, universes,
! surfaces, and lattices read in the input files.
!===============================================================================
subroutine print_summary()
subroutine print_geometry()
integer :: i
type(Surface), pointer :: s => null()
type(Cell), pointer :: c => null()
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
type(TallyObject), pointer :: t => null()
integer :: i
! print summary of cells
call header("CELL SUMMARY")
@ -686,6 +662,40 @@ contains
call print_surface(s)
end do
end subroutine print_geometry
!===============================================================================
! PRINT_SUMMARY displays summary information about the problem about to be run
! after reading all input files
!===============================================================================
subroutine print_summary()
integer :: i
character(15) :: string
type(Material), pointer :: m => null()
type(TallyObject), pointer :: t => null()
! Display problem summary
call header("PROBLEM SUMMARY")
if (problem_type == PROB_CRITICALITY) then
write(ou,100) 'Problem type:', 'Criticality'
write(ou,101) 'Number of Cycles:', n_cycles
write(ou,101) 'Number of Inactive Cycles:', n_inactive
elseif (problem_type == PROB_SOURCE) then
write(ou,100) 'Problem type:', 'External Source'
end if
write(ou,101) 'Number of Particles:', n_particles
! Display geometry summary
call header("GEOMETRY SUMMARY")
write(ou,101) 'Number of Cells:', n_cells
write(ou,101) 'Number of Surfaces:', n_surfaces
write(ou,101) 'Number of Materials:', n_materials
! print summary of all geometry
call print_geometry()
! print summary of materials
call header("MATERIAL SUMMARY")
do i = 1, n_materials
@ -702,6 +712,12 @@ contains
end do
end if
! print summary of unionized energy grid
call header("UNIONIZED ENERGY GRID")
write(ou,*) "Points on energy grid: " // trim(int_to_str(n_grid))
write(ou,*) "Extra storage required: " // trim(int_to_str(&
n_grid*n_nuclides_total*4)) // " bytes"
! print summary of variance reduction
call header("VARIANCE REDUCTION")
if (survival_biasing) then
@ -709,19 +725,16 @@ contains
else
write(ou,100) "Survival Biasing:", "off"
end if
write(ou,100) "Weight Cutoff:", trim(real_to_str(weight_cutoff))
write(ou,100) "Survival weight:", trim(real_to_str(weight_survive))
string = real_to_str(weight_cutoff)
write(ou,100) "Weight Cutoff:", trim(string)
string = real_to_str(weight_survive)
write(ou,100) "Survival weight:", trim(string)
write(ou,*)
! Format descriptor for columns
100 format (1X,A,T35,A)
101 format (1X,A,T35,I11)
nullify(s)
nullify(c)
nullify(u)
nullify(l)
nullify(m)
nullify(t)
end subroutine print_summary
@ -759,17 +772,26 @@ contains
subroutine print_runtime()
integer :: total_particles
real(8) :: speed
integer(8) :: total_particles
real(8) :: speed
character(15) :: string
! display header block
call header("Time Elapsed")
call header("Timing Statistics")
! display time elapsed for various sections
write(ou,100) "Total time elapsed", trim(real_to_str(time_total % elapsed))
write(ou,100) "Total time for initialization", trim(real_to_str(time_init % elapsed))
write(ou,100) "Total time in computation", trim(real_to_str(time_compute % elapsed))
write(ou,100) "Total time between cycles", trim(real_to_str(time_intercycle % elapsed))
write(ou,100) "Total time for initialization", time_initialize % elapsed
write(ou,100) " Reading cross sections", time_read_xs % elapsed
write(ou,100) " Unionizing energy grid", time_unionize % elapsed
write(ou,100) "Total time in computation", time_compute % elapsed
write(ou,100) "Total time between cycles", time_intercycle % elapsed
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
write(ou,100) " Reconstruct source bank", time_ic_rebuild % elapsed
write(ou,100) "Total time in inactive cycles", time_inactive % elapsed
write(ou,100) "Total time in active cycles", time_active % elapsed
write(ou,100) "Total time elapsed", time_total % elapsed
! display header block
call header("Run Statistics")
@ -777,14 +799,15 @@ contains
! display calculate rate and final keff
total_particles = n_particles * n_cycles
speed = real(total_particles) / time_compute % elapsed
write(ou,101) "Calculation Rate", trim(real_to_str(speed))
write(ou,102) "Final Keff", trim(real_to_str(keff)), trim(real_to_str(keff_std))
string = real_to_str(speed)
write(ou,101) "Calculation Rate", trim(string)
write(ou,102) "Final Keff", keff, keff_std
write(ou,*)
! format for write statements
100 format (1X,A,T33,"= ",A," seconds")
100 format (1X,A,T35,"= ",ES11.4," seconds")
101 format (1X,A,T20,"= ",A," neutrons/second")
102 format (1X,A,T20,"= ",A," +/- ",A)
102 format (1X,A,T20,"= ",F8.5," +/- ",F8.5)
end subroutine print_runtime

134
src/particle_header.F90 Normal file
View file

@ -0,0 +1,134 @@
module particle_header
use constants, only: NEUTRON, ONE, NONE
use geometry_header, only: BASE_UNIVERSE
implicit none
!===============================================================================
! LOCALCOORD describes the location of a particle local to a single
! universe. When the geometry consists of nested universes, a particle will have
! a list of coordinates in each level
!===============================================================================
type LocalCoord
! Indices in various arrays for this level
integer :: cell = NONE
integer :: universe = NONE
integer :: lattice = NONE
integer :: lattice_x = NONE
integer :: lattice_y = NONE
! Particle position and direction for this level
real(8) :: xyz(3)
real(8) :: uvw(3)
! Pointers to next (lower) and previous (higher) universe
type(LocalCoord), pointer :: next => null()
type(LocalCoord), pointer :: prev => null()
end type LocalCoord
!===============================================================================
! PARTICLE describes the state of a particle being transported through the
! geometry
!===============================================================================
type Particle
! Basic data
integer(8) :: id ! Unique ID
integer :: type ! Particle type (n, p, e, etc)
! Particle coordinates
logical :: in_lower_universe ! is particle in lower universe?
type(LocalCoord), pointer :: coord0 => null() ! coordinates on universe 0
type(LocalCoord), pointer :: coord => null() ! coordinates on lowest universe
! Other physical data
real(8) :: wgt ! particle weight
real(8) :: E ! energy
real(8) :: mu ! angle of scatter
logical :: alive ! is particle alive?
! Pre-collision physical data
real(8) :: last_xyz(3) ! previous coordinates
real(8) :: last_wgt ! last particle weight
real(8) :: last_E ! last energy
! Post-collision physical data
integer :: n_bank ! number of fission sites banked
! Energy grid data
integer :: IE ! index on energy grid
real(8) :: interp ! interpolation factor for energy grid
! Indices for various arrays
integer :: surface ! index for surface particle is on
integer :: cell_born ! index for cell particle was born in
integer :: material ! index for current material
integer :: last_material ! index for last material
! Statistical data
integer :: n_collision ! # of collisions
end type Particle
contains
!===============================================================================
! INITIALIZE_PARTICLE sets default attributes for a particle from the source
! bank
!===============================================================================
subroutine initialize_particle(p)
type(Particle), pointer :: p
! TODO: if information on the cell, lattice, universe, and material is
! passed through the fission bank to the source bank, no lookup would be
! needed at the beginning of a cycle
p % type = NEUTRON
p % alive = .true.
! clear attributes
p % surface = NONE
p % cell_born = NONE
p % material = NONE
p % last_material = NONE
p % wgt = ONE
p % last_wgt = ONE
p % n_bank = 0
p % n_collision = 0
! remove any original coordinates
call deallocate_coord(p % coord0)
! Set up base level coordinates
allocate(p % coord0)
p % coord0 % universe = BASE_UNIVERSE
p % coord => p % coord0
p % in_lower_universe = .false.
end subroutine initialize_particle
!===============================================================================
! DEALLOCATE_COORD removes all levels of coordinates below a given level. This
! is used in distance_to_boundary when the particle moves from a lower universe
! to a higher universe since the data for the lower one is not needed anymore.
!===============================================================================
recursive subroutine deallocate_coord(coord)
type(LocalCoord), pointer :: coord
if (associated(coord)) then
! recursively deallocate lower coordinates
if (associated(coord % next)) call deallocate_coord(coord%next)
! deallocate original coordinate
deallocate(coord)
end if
end subroutine deallocate_coord
end module particle_header

87
src/particle_header.f90
View file

@ -1,87 +0,0 @@
module particle_header
use constants, only: NEUTRON, ONE
implicit none
!===============================================================================
! PARTICLE describes the state of a particle being transported through the
! geometry
!===============================================================================
type Particle
! Basic data
integer(8) :: id ! Unique ID
integer :: type ! Particle type (n, p, e, etc)
! Physical data
real(8) :: xyz(3) ! location
real(8) :: xyz_local(3) ! local location (after transformations)
real(8) :: uvw(3) ! directional cosines
real(8) :: wgt ! particle weight
real(8) :: E ! energy
real(8) :: mu ! angle of scatter
logical :: alive ! is particle alive?
! Pre-collision physical data
real(8) :: last_xyz(3) ! previous coordinates
real(8) :: last_wgt ! last particle weight
real(8) :: last_E ! last energy
! Post-collision physical data
integer :: n_bank ! number of fission sites banked
! Energy grid data
integer :: IE ! index on energy grid
real(8) :: interp ! interpolation factor for energy grid
! Indices for various arrays
integer :: cell ! index for current cell
integer :: cell_born ! index for cell particle was born in
integer :: universe ! index for current universe
integer :: lattice ! index for current lattice
integer :: surface ! index for current surface
integer :: material ! index for current material
integer :: last_material ! index for last material
integer :: index_x ! lattice index for x direction
integer :: index_y ! lattice index for y direction
! Statistical data
integer :: n_collision ! # of collisions
end type Particle
contains
!===============================================================================
! INITIALIZE_PARTICLE sets default attributes for a particle from the source
! bank
!===============================================================================
subroutine initialize_particle(p)
type(Particle), pointer :: p
! TODO: if information on the cell, lattice, universe, and material is
! passed through the fission bank to the source bank, no lookup would be
! needed at the beginning of a cycle
p % type = NEUTRON
p % wgt = ONE
p % last_wgt = ONE
p % alive = .true.
p % n_bank = 0
p % cell = 0
p % cell_born = 0
p % universe = 0
p % lattice = 0
p % surface = 0
p % material = 0
p % last_material = 0
p % index_x = 0
p % index_y = 0
p % n_collision = 0
end subroutine initialize_particle
end module particle_header

782
src/physics.f90 → src/physics.F90
View file

File diff suppressed because it is too large Load diff

191
src/plot.F90 Normal file
View file

@ -0,0 +1,191 @@
module plot
use constants
use error, only: fatal_error
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, cell_contains
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use particle_header, only: Particle, initialize_particle, LocalCoord, &
deallocate_coord
implicit none
contains
!===============================================================================
! RUN_PLOT generates a binary stream file containing a list of surface/lattice
! crossings and what cell was traveled through. A Python script can then be used
! to generate a plot based on the recorded crossings and cells
!===============================================================================
subroutine run_plot()
integer :: i ! loop index
integer :: surface_crossed ! surface which particle is on
integer :: last_cell ! most recent cell particle was in
integer :: enter_surface ! entrance surface
real(8) :: xyz(3) ! starting coordinates
real(8) :: last_x_coord ! bounding x coordinate
real(8) :: last_y_coord ! bounding y coordinate
real(8) :: d ! distance to boundary
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
logical :: lattice_crossed ! is surface crossing in lattice?
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
type(Cell), pointer :: c => null()
type(Universe), pointer :: univ => null()
type(Particle), pointer :: p => null()
type(LocalCoord), pointer :: coord => null()
! Open plot file for binary writing
path_plot = trim(path_input) // "plot.out"
open(UNIT=UNIT_PLOT, FILE=path_plot, STATUS="replace", ACCESS="stream")
! Write origin, width, basis, and pixel width to file
write(UNIT=UNIT_PLOT) plot_origin
write(UNIT=UNIT_PLOT) plot_width
write(UNIT=UNIT_PLOT) plot_basis
write(UNIT=UNIT_PLOT) pixel
! Determine coordinates of the upper-left corner of the plot
xyz(1) = plot_origin(1) - plot_width(1) / 2.0
xyz(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
xyz(3) = plot_origin(3)
! Determine bounding x and y coordinates for plot
last_x_coord = plot_origin(1) + plot_width(1) / 2.0
last_y_coord = plot_origin(2) - plot_width(2) / 2.0
! allocate and initialize particle
allocate(p)
! loop over horizontal rays
do while(xyz(2) > last_y_coord)
! initialize the particle and set starting coordinate and direction
call initialize_particle(p)
p % coord % xyz = xyz
p % coord % uvw = (/ 1, 0, 0 /)
! write starting coordinate to file
write(UNIT=UNIT_PLOT) p % coord % xyz
! Find cell that particle is currently in
call find_cell(p, found_cell)
! =======================================================================
! MOVE PARTICLE FORWARD TO NEXT CELL
if (.not. found_cell) then
! Clear any coordinates beyond first level
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
distance = INFINITY
univ => universes(BASE_UNIVERSE)
do i = 1, univ % n_cells
p % coord0 % xyz = xyz
p % coord0 % cell = univ % cells(i)
call distance_to_boundary(p, d, surface_crossed, lattice_crossed)
if (d < distance) then
! Check to make sure particle is actually going into this cell
! by moving it slightly forward and seeing if the cell contains
! that coordinate
p % coord0 % xyz = p % coord0 % xyz + (d + TINY_BIT) * p % coord0 % uvw
c => cells(p % coord0 % cell)
if (.not. cell_contains(c, p)) cycle
! Set new distance and retain pointer to this cell
distance = d
enter_surface = surface_crossed
end if
end do
! No cell was found on this horizontal ray
if (distance == INFINITY) then
p % coord0 % xyz(1) = last_x_coord
write(UNIT_PLOT) p % coord0 % xyz, 0
! Move to next horizontal ray
xyz(2) = xyz(2) - pixel
cycle
end if
! Write coordinate where next cell begins
p % coord0 % xyz = xyz + distance * p % coord0 % uvw
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
! Process surface crossing for next cell
p % coord0 % cell = NONE
p % surface = -enter_surface
call cross_surface(p, enter_surface)
end if
! =======================================================================
! MOVE PARTICLE ACROSS HORIZONTAL TRACK
do while (p % alive)
! save particle's current cell
last_cell = p % coord % cell
! Calculate distance to next boundary
call distance_to_boundary(p, distance, surface_crossed, lattice_crossed)
! Advance particle
coord => p % coord0
do while (associated(coord))
coord % xyz = coord % xyz + distance * coord % uvw
coord => coord % next
end do
! If next boundary crossing is out of range of the plot, only include
! the visible portion and move to next horizontal ray
if (p % coord0 % xyz(1) >= last_x_coord) then
p % alive = .false.
p % coord0 % xyz(1) = last_x_coord
! If there is no cell beyond this boundary, mark it as cell 0
if (distance == INFINITY) p % coord % cell = 0
! Write ending coordinates to file
write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
cycle
end if
! Write boundary crossing coordinates to file
write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
p % coord % cell = 0
if (lattice_crossed) then
p % surface = NONE
call cross_lattice(p)
else
p % surface = surface_crossed
call cross_surface(p, last_cell)
! Since boundary conditions are disabled in plotting mode, we need
! to manually add the last segment
if (surfaces(abs(surface_crossed)) % bc == BC_VACUUM) then
p % coord0 % xyz(1) = last_x_coord
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
exit
end if
end if
end do
! Move y-coordinate to next position
xyz(2) = xyz(2) - pixel
end do
! Close plot file
close(UNIT=UNIT_PLOT)
end subroutine run_plot
end module plot

188
src/plot.f90
View file

@ -1,188 +0,0 @@
module plot
use constants
use error, only: fatal_error
use geometry, only: find_cell, dist_to_boundary, cross_surface, &
cross_lattice, cell_contains
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use particle_header, only: Particle, initialize_particle
implicit none
contains
!===============================================================================
! RUN_PLOT
!===============================================================================
subroutine run_plot()
integer :: i ! loop index
integer :: surf ! surface which particle is on
integer :: last_cell ! most recent cell particle was in
real(8) :: coord(3) ! starting coordinates
real(8) :: last_x_coord ! bounding x coordinate
real(8) :: last_y_coord ! bounding y coordinate
real(8) :: d ! distance to boundary
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
logical :: in_lattice ! is surface crossing in lattice?
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
type(Cell), pointer :: c => null()
type(Universe), pointer :: univ => null()
type(Particle), pointer :: p => null()
! Open plot file for binary writing
path_plot = trim(path_input) // "plot.out"
open(UNIT=UNIT_PLOT, FILE=path_plot, STATUS="replace", ACCESS="stream")
! Write origin, width, basis, and pixel width to file
write(UNIT=UNIT_PLOT) plot_origin
write(UNIT=UNIT_PLOT) plot_width
write(UNIT=UNIT_PLOT) plot_basis
write(UNIT=UNIT_PLOT) pixel
! Determine coordinates of the upper-left corner of the plot
coord(1) = plot_origin(1) - plot_width(1) / 2.0
coord(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
coord(3) = plot_origin(3)
! Determine bounding x and y coordinates for plot
last_x_coord = plot_origin(1) + plot_width(1) / 2.0
last_y_coord = plot_origin(2) - plot_width(2) / 2.0
! allocate and initialize particle
allocate(p)
! loop over horizontal rays
do while(coord(2) > last_y_coord)
! initialize the particle and set starting coordinate and direction
call initialize_particle(p)
p % xyz = coord
p % xyz_local = coord
p % uvw = (/ 1, 0, 0 /)
! write starting coordinate to file
write(UNIT=UNIT_PLOT) p % xyz
! Find cell that particle is currently in
univ => universes(BASE_UNIVERSE)
call find_cell(univ, p, found_cell)
! =======================================================================
! MOVE PARTICLE FORWARD TO NEXT CELL
if (.not. found_cell) then
univ => universes(BASE_UNIVERSE)
do i = 1, univ % n_cells
p % xyz = coord
p % xyz_local = coord
p % cell = univ % cells(i)
distance = INFINITY
call dist_to_boundary(p, d, surf, in_lattice)
if (d < distance) then
! Move particle forward to next surface
p % xyz = p % xyz + d * p % uvw
! Check to make sure particle is actually going into this cell
! by moving it slightly forward and seeing if the cell contains
! that coordinate
p % xyz = p % xyz + 1e-4 * p % uvw
p % xyz_local = p % xyz
c => cells(p % cell)
if (.not. cell_contains(c, p)) cycle
! Reset coordinate to surface crossing
p % xyz = p % xyz - 1e-4 * p % uvw
p % xyz_local = p % xyz
! Set new distance and retain pointer to this cell
distance = d
last_cell = p % cell
end if
end do
! No cell was found on this horizontal ray
if (distance == INFINITY) then
p % xyz(1) = last_x_coord
p % cell = 0
write(UNIT_PLOT) p % xyz, p % cell
! Move to next horizontal ray
coord(2) = coord(2) - pixel
cycle
end if
! Write coordinate where next cell begins
write(UNIT=UNIT_PLOT) p % xyz, 0
! Process surface crossing for next cell
p % cell = 0
p % surface = -surf
call cross_surface(p, last_cell)
end if
! =======================================================================
! MOVE PARTICLE ACROSS HORIZONTAL TRACK
do while (p % alive)
! Calculate distance to next boundary
call dist_to_boundary(p, distance, surf, in_lattice)
! Advance particle
p%xyz = p%xyz + distance * p%uvw
p%xyz_local = p%xyz_local + distance * p%uvw
! If next boundary crossing is out of range of the plot, only include
! the visible portion and move to next horizontal ray
if (p % xyz(1) >= last_x_coord) then
p % alive = .false.
p % xyz(1) = last_x_coord
! If there is no cell beyond this boundary, mark it as cell 0
if (distance == INFINITY) p % cell = 0
! Write ending coordinates to file
write(UNIT=UNIT_PLOT) p % xyz, p % cell
cycle
end if
! Write boundary crossing coordinates to file
write(UNIT=UNIT_PLOT) p % xyz, p % cell
last_cell = p % cell
p % cell = 0
if (in_lattice) then
p % surface = 0
call cross_lattice(p)
else
p % surface = surf
call cross_surface(p, last_cell)
! Since boundary conditions are disabled in plotting mode, we need
! to manually add the last segment
if (surfaces(surf) % bc == BC_VACUUM) then
p % xyz(1) = last_x_coord
write(UNIT=UNIT_PLOT) p % xyz, 0
exit
end if
end if
end do
! Move y-coordinate to next position
coord(2) = coord(2) - pixel
end do
! Close plot file
close(UNIT=UNIT_PLOT)
end subroutine run_plot
end module plot

147
src/random_lcg.F90 Normal file
View file

@ -0,0 +1,147 @@
module random_lcg
implicit none
private
save
integer(8) :: prn_seed0 ! original seed
integer(8) :: prn_seed ! current seed
integer(8) :: prn_mult ! multiplication factor, g
integer(8) :: prn_add ! additive factor, c
integer :: prn_bits ! number of bits, M
integer(8) :: prn_mod ! 2^M
integer(8) :: prn_mask ! 2^M - 1
integer(8) :: prn_stride ! stride between particles
real(8) :: prn_norm ! 2^(-M)
public :: prn
public :: initialize_prng
public :: set_particle_seed
public :: prn_skip
contains
!===============================================================================
! PRN generates a pseudo-random number using a linean congruential generator
!===============================================================================
function prn() result(pseudo_rn)
real(8) :: pseudo_rn
! This algorithm uses bit-masking to find the next integer(8) value to be
! used to calculate the random number
prn_seed = iand(prn_mult*prn_seed + prn_add, prn_mask)
! Once the integer is calculated, we just need to divide by 2**m,
! represented here as multiplying by a pre-calculated factor
pseudo_rn = prn_seed * prn_norm
end function prn
!===============================================================================
! INITIALIZE_PRNG sets up the random number generator, determining the seed and
! values for g, c, and m.
!===============================================================================
subroutine initialize_prng()
prn_seed0 = 1_8
prn_seed = prn_seed
prn_mult = 2806196910506780709_8
prn_add = 1_8
prn_bits = 63
prn_mod = ibset(0_8, prn_bits) ! clever way of calculating 2**bits
prn_mask = prn_mod - 1_8
prn_stride = 152917_8
prn_norm = 2._8**(-prn_bits)
end subroutine initialize_prng
!===============================================================================
! SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the
! particle
!===============================================================================
subroutine set_particle_seed(id)
integer(8), intent(in) :: id
prn_seed = prn_skip_ahead(id*prn_stride, prn_seed0)
end subroutine set_particle_seed
!===============================================================================
! PRN_SKIP advances the random number seed 'n' times from the current seed
!===============================================================================
subroutine prn_skip(n)
integer(8), intent(in) :: n ! number of seeds to skip
prn_seed = prn_skip_ahead(n, prn_seed)
end subroutine prn_skip
!===============================================================================
! PRN_SKIP_AHEAD advances the random number seed 'skip' times. This is usually
! used to skip a fixed number of random numbers (the stride) so that a given
! particle always has the same starting seed regardless of how many processors
! are used
!===============================================================================
function prn_skip_ahead(n, seed) result(new_seed)
integer(8), intent(in) :: n ! number of seeds to skip
integer(8), intent(in) :: seed ! original seed
integer(8) :: new_seed ! new seed
integer(8) :: nskip
integer(8) :: gen
integer(8) :: g
integer(8) :: inc
integer(8) :: c
integer(8) :: gp
! In cases where we want to skip backwards, we add the period of the random
! number generator until the number of PRNs to skip is positive since
! skipping ahead that much is the same as skipping backwards by the original
! amount
nskip = n
do while (nskip < 0_8)
nskip = nskip + prn_mod
enddo
! The algorithm here to determine the parameters used to skip ahead is
! described in F. Brown, "Random Number Generation with Arbitrary Stride,"
! Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
! O(log2(N)) operations instead of O(N). Basically, it computes parameters G
! and C which can then be used to find x_N = G*x_0 + C mod 2^M.
nskip = iand(nskip, prn_mask)
gen = 1
g = prn_mult
inc = 0
c = prn_add
do while (nskip > 0_8)
if (btest(nskip,0)) then
gen = iand(gen*g, prn_mask)
inc = iand(inc*g, prn_mask)
inc = iand(inc+c, prn_mask)
endif
gp = iand(g+1, prn_mask)
g = iand(g*g, prn_mask)
c = iand(gp*c, prn_mask)
nskip = ishft(nskip, -1)
enddo
! With G and C, we can now find the new seed
new_seed = iand(gen*seed + inc, prn_mask)
end function prn_skip_ahead
end module random_lcg

47
src/search.f90 → src/search.F90
View file

@ -4,6 +4,8 @@ module search
use error, only: fatal_error
use global, only: message
integer, parameter :: MAX_ITERATION = 64
contains
!===============================================================================
@ -11,15 +13,16 @@ contains
! value lies in the array. This is used extensively for energy grid searching
!===============================================================================
function binary_search(array, n, val) result(index)
function binary_search(array, n, val) result(array_index)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
real(8), intent(in) :: val
integer :: index
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
real(8) :: testval
L = 1
@ -30,45 +33,37 @@ contains
call fatal_error()
end if
n_iteration = 0
do while (R - L > 1)
! Check boundaries
if (val > array(L) .and. val < array(L+1)) then
index = L
array_index = L
return
elseif (val > array(R-1) .and. val < array(R)) then
index = R-1
array_index = R - 1
return
end if
! Find values at midpoint
index = L + (R - L)/2
testval = array(index)
if (val > testval) then
L = index
array_index = L + (R - L)/2
testval = array(array_index)
if (val >= testval) then
L = array_index
elseif (val < testval) then
R = index
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
message = "Reached maximum number of iterations on binary search."
call fatal_error()
end if
end do
index = L
array_index = L
end function binary_search
!===============================================================================
! INTERPOLATE
!===============================================================================
function interpolate(array, n, index, f) result(val)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
integer, intent(in) :: index
real(8), intent(in) :: f
real(8) :: val
val = (ONE-f) * array(index) + f * array(index+1)
end function interpolate
end module search

68
src/source.f90 → src/source.F90
View file

@ -3,11 +3,13 @@ module source
use bank_header, only: Bank
use constants, only: ONE, MAX_LINE_LEN
use cross_section_header, only: Nuclide
use error, only: fatal_error
use global
use mcnp_random, only: rang, RN_init_particle
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use physics, only: watt_spectrum
use random_lcg, only: prn, set_particle_seed
use string, only: int_to_str
implicit none
@ -24,12 +26,15 @@ contains
integer(8) :: j ! loop index over bank sites
integer :: k ! dummy loop index
integer(8) :: maxwork ! maxinum # of particles per processor
integer :: alloc_err ! allocation error code
integer(8) :: bytes ! size of fission/source bank
real(8) :: r(3) ! sampled coordinates
real(8) :: phi ! azimuthal angle
real(8) :: mu ! cosine of polar angle
real(8) :: E ! outgoing energy
real(8) :: p_min(3) ! minimum coordinates of source
real(8) :: p_max(3) ! maximum coordinates of source
type(Bank) :: bank_obj
message = "Initializing source particles..."
call write_message(6)
@ -37,9 +42,31 @@ contains
! Determine maximum amount of particles to simulate on each processor
maxwork = ceiling(real(n_particles)/n_procs,8)
! Allocate fission and source banks
allocate(source_bank(maxwork))
allocate(fission_bank(3*maxwork))
! Allocate source bank
allocate(source_bank(maxwork), STAT=alloc_err)
if (alloc_err /= 0) then
#ifndef NO_F2008
bytes = maxwork * storage_size(bank_obj) / 8
#else
bytes = maxwork * 64 / 8
#endif
message = "Could not allocate source bank. Attempted to allocate " &
// trim(int_to_str(bytes)) // " bytes."
call fatal_error()
end if
! Allocate fission bank
allocate(fission_bank(3*maxwork), STAT=alloc_err)
if (alloc_err /= 0) then
#ifndef NO_F2008
bytes = 3 * maxwork * storage_size(bank_obj) / 8
#else
bytes = 3 * maxwork * 64 / 8
#endif
message = "Could not allocate fission bank. Attempted to allocate " &
// trim(int_to_str(bytes)) // " bytes."
call fatal_error()
end if
! Check external source type
if (external_source%type == SRC_BOX) then
@ -60,26 +87,25 @@ contains
do j = bank_first, bank_last
p => source_bank(j - bank_first + 1)
! initialize random number seed
call RN_init_particle(int(j,8))
! sample position
r = (/ (rang(), k = 1,3) /)
p % id = j
p % xyz = p_min + r*(p_max - p_min)
p % xyz_local = p % xyz
p % last_xyz = p % xyz
! sample angle
phi = TWO*PI*rang()
mu = TWO*rang() - ONE
p % uvw(1) = mu
p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! set defaults
call initialize_particle(p)
! initialize random number seed
call set_particle_seed(int(j,8))
! sample position
r = (/ (prn(), k = 1,3) /)
p % id = j
p % coord0 % xyz = p_min + r*(p_max - p_min)
p % last_xyz = p % coord0 % xyz
! sample angle
phi = TWO*PI*prn()
mu = TWO*prn() - ONE
p % coord0 % uvw(1) = mu
p % coord0 % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
p % coord0 % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! sample energy from Watt fission energy spectrum for U-235
do
E = watt_spectrum(0.988_8, 2.249_8)

0
src/source_header.f90 → src/source_header.F90
View file

20
src/string.f90 → src/string.F90
View file

@ -28,7 +28,7 @@ contains
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
character(1) :: char ! current character
character(1) :: chr ! current character
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
@ -37,14 +37,14 @@ contains
i_end = 0
n = 0
do i = 1, len_trim(string)
char = string(i:i)
chr = string(i:i)
! Note that ACHAR(9) is a horizontal tab
if ((i_start == 0) .and. (char /= ' ') .and. (char /= achar(9))) then
if ((i_start == 0) .and. (chr /= ' ') .and. (chr /= achar(9))) then
i_start = i
end if
if (i_start > 0) then
if ((char == ' ') .or. (char == achar(9))) i_end = i - 1
if ((chr == ' ') .or. (chr == achar(9))) i_end = i - 1
if (i == len_trim(string)) i_end = i
if (i_end > 0) then
n = n + 1
@ -80,7 +80,7 @@ contains
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
character(1) :: char ! current character
character(1) :: chr ! current character
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
@ -89,27 +89,27 @@ contains
i_end = 0
n = 0
do i = 1, len_trim(string)
char = string(i:i)
chr = string(i:i)
! Check for special characters
if (index('():#', char) > 0) then
if (index('():#', chr) > 0) then
if (i_start > 0) then
i_end = i - 1
n = n + 1
words(n) = string(i_start:i_end)
end if
n = n + 1
words(n) = char
words(n) = chr
i_start = 0
i_end = 0
cycle
end if
if ((i_start == 0) .and. (char /= ' ')) then
if ((i_start == 0) .and. (chr /= ' ')) then
i_start = i
end if
if (i_start > 0) then
if (char == ' ') i_end = i - 1
if (chr == ' ') i_end = i - 1
if (i == len_trim(string)) i_end = i
if (i_end > 0) then
n = n + 1

533
src/tally.f90 → src/tally.F90
View file

@ -1,5 +1,7 @@
module tally
use ISO_FORTRAN_ENV
use constants
use error, only: fatal_error
use global
@ -22,67 +24,77 @@ module tally
contains
!===============================================================================
! CALCULATE_KEFF
! 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
!===============================================================================
subroutine calculate_keff(i_cycle)
integer, intent(in) :: i_cycle ! index of current cycle
integer(8) :: total_bank ! total number of source sites
integer :: n ! active cycle number
real(8) :: kcoll ! keff collision estimator
real(8), save :: k1 = 0. ! accumulated keff
real(8), save :: k2 = 0. ! accumulated keff**2
real(8) :: std ! stdev of keff over active cycles
#ifdef MPI
integer :: ierr
#endif
integer(8) :: total_bank ! total number of source sites
integer :: n ! active cycle number
real(8) :: k_cycle ! single cycle estimate of keff
real(8), save :: k_sum ! accumulated keff
real(8), save :: k_sum_sq ! accumulated keff**2
message = "Calculate cycle keff..."
call write_message(8)
! set k1 and k2 at beginning of run
! initialize sum and square of sum at beginning of run
if (i_cycle == 1) then
k1 = ZERO
k2 = ZERO
k_sum = ZERO
k_sum_sq = ZERO
end if
#ifdef MPI
! Collect number bank sites onto master process
call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, &
& MPI_COMM_WORLD, ierr)
& MPI_COMM_WORLD, mpi_err)
#else
total_bank = n_bank
#endif
! Collect statistics and print output
if (master) then
kcoll = real(total_bank)/real(n_particles)*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
k_cycle = real(total_bank)/real(n_particles)*keff
if (i_cycle > n_inactive) then
! Active cycle number
n = i_cycle - n_inactive
k1 = k1 + kcoll
k2 = k2 + kcoll**2
keff = k1/n
std = sqrt((k2/n-keff**2)/n)
keff_std = std
! Accumulate cycle estimate of k
k_sum = k_sum + k_cycle
k_sum_sq = k_sum_sq + k_cycle*k_cycle
! Determine mean and standard deviation of mean
keff = k_sum/n
keff_std = sqrt((k_sum_sq/n - keff*keff)/n)
! Display output for this cycle
if (i_cycle > n_inactive+1) then
write(6,101) i_cycle, kcoll, keff, std
write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std
else
write(6,100) i_cycle, kcoll
write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
end if
else
write(6,100) i_cycle, kcoll
keff = kcoll
! Display output for inactive cycle
write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
keff = k_cycle
end if
end if
#ifdef MPI
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
! Broadcast new keff value to all processors
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
#endif
100 format (2X,I4,2X,F8.5)
101 format (2X,I4,2X,F8.5,9X,F8.5,1X,F8.5)
100 format (2X,I5,2X,F8.5)
101 format (2X,I5,2X,F8.5,9X,F8.5,1X,F8.5)
end subroutine calculate_keff
@ -95,12 +107,12 @@ contains
subroutine create_tally_map()
integer :: i ! loop index for tallies
integer :: j ! loop index for filter arrays
integer :: index ! filter bin entries
integer :: n ! number of bins
integer :: filter_bins ! running total of number of filter bins
integer :: score_bins ! number of scoring bins
integer :: i ! loop index for tallies
integer :: j ! loop index for filter arrays
integer :: i_item ! filter bin entries
integer :: n ! number of bins
integer :: filter_bins ! running total of number of filter bins
integer :: score_bins ! number of scoring bins
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
@ -129,23 +141,31 @@ contains
if (t % surface_current) then
m => meshes(t % mesh)
t % stride(TS_SURFACE) = filter_bins
! Set stride for surface/direction
if (m % n_dimension == 2) then
filter_bins = filter_bins * 4
elseif (m % n_dimension == 3) then
filter_bins = filter_bins * 6
end if
! Add filter for incoming energy
n = t % n_bins(T_ENERGYIN)
t % stride(TS_ENERGYIN) = filter_bins
if (n > 0) then
filter_bins = filter_bins * n
end if
! account for z direction
t % stride(3) = filter_bins
t % stride(TS_MESH_Z) = filter_bins
filter_bins = filter_bins * (m % dimension(3) + 1)
! account for y direction
t % stride(2) = filter_bins
t % stride(TS_MESH_Y) = filter_bins
filter_bins = filter_bins * (m % dimension(2) + 1)
! account for z direction
t % stride(1) = filter_bins
t % stride(TS_MESH_X) = filter_bins
filter_bins = filter_bins * (m % dimension(1) + 1)
! Finally add scoring bins for the macro tallies and allocate scores
@ -185,8 +205,8 @@ contains
t % stride(T_SURFACE) = filter_bins
if (n > 0) then
do j = 1, n
index = t % surface_bins(j) % scalar
call add_map_element(tally_maps(T_SURFACE) % items(index), i, j)
i_item = t % surface_bins(j) % scalar
call add_map_element(tally_maps(T_SURFACE) % items(i_item), i, j)
end do
filter_bins = filter_bins * n
end if
@ -196,8 +216,8 @@ contains
t % stride(T_CELLBORN) = filter_bins
if (n > 0) then
do j = 1, n
index = t % cellborn_bins(j) % scalar
call add_map_element(tally_maps(T_CELLBORN) % items(index), i, j)
i_item = t % cellborn_bins(j) % scalar
call add_map_element(tally_maps(T_CELLBORN) % items(i_item), i, j)
end do
filter_bins = filter_bins * n
end if
@ -207,8 +227,8 @@ contains
t % stride(T_CELL) = filter_bins
if (n > 0) then
do j = 1, n
index = t % cell_bins(j) % scalar
call add_map_element(tally_maps(T_CELL) % items(index), i, j)
i_item = t % cell_bins(j) % scalar
call add_map_element(tally_maps(T_CELL) % items(i_item), i, j)
end do
filter_bins = filter_bins * n
end if
@ -218,8 +238,8 @@ contains
t % stride(T_MATERIAL) = filter_bins
if (n > 0) then
do j = 1, n
index = t % material_bins(j) % scalar
call add_map_element(tally_maps(T_MATERIAL) % items(index), i, j)
i_item = t % material_bins(j) % scalar
call add_map_element(tally_maps(T_MATERIAL) % items(i_item), i, j)
end do
filter_bins = filter_bins * n
end if
@ -229,8 +249,8 @@ contains
t % stride(T_UNIVERSE) = filter_bins
if (n > 0) then
do j = 1, n
index = t % universe_bins(j) % scalar
call add_map_element(tally_maps(T_UNIVERSE) % items(index), i, j)
i_item = t % universe_bins(j) % scalar
call add_map_element(tally_maps(T_UNIVERSE) % items(i_item), i, j)
end do
filter_bins = filter_bins * n
end if
@ -313,7 +333,6 @@ contains
real(8) :: wgt
real(8) :: mu
real(8) :: E_out
logical :: in_mesh
logical :: has_energyout_bin
logical :: analog
type(TallyObject), pointer :: t
@ -336,9 +355,21 @@ contains
! =======================================================================
! DETERMINE SCORING BIN COMBINATION
! determine mesh bin
if (t % n_bins(T_MESH) > 0) then
m => meshes(t % mesh)
! Determine if we're in the mesh first
call get_mesh_bin(m, p % coord0 % xyz, mesh_bin)
if (mesh_bin == NO_BIN_FOUND) cycle
bins(T_MESH) = mesh_bin
else
bins(T_MESH) = 1
end if
! determine next universe bin
if (t % n_bins(T_UNIVERSE) > 0) then
bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % universe, i)
bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % coord % universe, i)
if (bins(T_UNIVERSE) == NO_BIN_FOUND) cycle
else
bins(T_UNIVERSE) = 1
@ -354,7 +385,7 @@ contains
! determine next cell bin
if (t % n_bins(T_CELL) > 0) then
bins(T_CELL) = get_next_bin(T_CELL, p % cell, i)
bins(T_CELL) = get_next_bin(T_CELL, p % coord % cell, i)
if (bins(T_CELL) == NO_BIN_FOUND) cycle
else
bins(T_CELL) = 1
@ -376,19 +407,6 @@ contains
bins(T_SURFACE) = 1
end if
! determine mesh bin
if (t % n_bins(T_MESH) > 0) then
m => meshes(t % mesh)
! Determine if we're in the mesh first
call get_mesh_bin(m, p % xyz, mesh_bin, in_mesh)
if (.not. in_mesh) cycle
bins(T_MESH) = mesh_bin
else
bins(T_MESH) = 1
end if
! determine incoming energy bin
n = t % n_bins(T_ENERGYIN)
if (n > 0) then
@ -507,6 +525,13 @@ contains
bin_energyout = bins(T_ENERGYOUT)
score_index0 = score_index
! Since the creation of fission sites is weighted such that it
! is expected to create n_particles sites, we need to multiply
! the score by keff to get the true nu-fission rate. Otherwise,
! the sum of all nu-fission rates would be ~1.0.
score = keff
! loop over number of particles banked
do k = 1, p % n_bank
! determine outgoing energy from fission bank
@ -518,14 +543,6 @@ contains
! determine scoring index
score_index = sum((bins - 1) * t % stride) + 1
! Since the creation of fission sites is weighted such that
! it is expected to create n_particles sites, we need to
! multiply the score by keff to get the true nu-fission
! rate. Otherwise, the sum of all nu-fission rates would be
! ~1.0.
score = keff
! Add score to tally
call add_to_score(t % scores(score_index, j), score)
end do
@ -578,31 +595,35 @@ contains
type(Particle), pointer :: p
integer :: i ! loop indices
integer :: j ! loop indices
integer :: k ! loop indices
integer :: ijk0(3) ! indices of starting coordinates
integer :: ijk1(3) ! indices of ending coordinates
integer :: n_cross ! number of surface crossings
integer :: score_index ! index of scoring bin
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 bounding surfaces
real(8) :: d(3) ! distance to each bounding surface
real(8) :: distance ! actual distance traveled
logical :: start_in_mesh ! particle's starting xyz in mesh?
logical :: end_in_mesh ! particle's ending xyz in mesh?
logical :: x_same ! same starting/ending x index (i)
logical :: y_same ! same starting/ending y index (j)
logical :: z_same ! same starting/ending z index (k)
integer :: i ! loop indices
integer :: j ! loop indices
integer :: k ! loop indices
integer :: ijk0(3) ! indices of starting coordinates
integer :: ijk1(3) ! indices of ending coordinates
integer :: n_cross ! number of surface crossings
integer :: n ! number of incoming energy bins
integer :: bins(TALLY_TYPES) ! scoring bin combination
integer :: score_index ! index of scoring bin
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 bounding surfaces
real(8) :: d(3) ! distance to each bounding surface
real(8) :: distance ! actual distance traveled
logical :: start_in_mesh ! particle's starting xyz in mesh?
logical :: end_in_mesh ! particle's ending xyz in mesh?
logical :: x_same ! same starting/ending x index (i)
logical :: y_same ! same starting/ending y index (j)
logical :: z_same ! same starting/ending z index (k)
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
bins = 1
do i = 1, n_tallies
! Copy starting and ending location of particle
xyz0 = p % last_xyz
xyz1 = p % xyz
xyz1 = p % coord0 % xyz
! Get pointer to tally
t => tallies(i)
@ -623,7 +644,20 @@ contains
if (n_cross == 0) cycle
! Copy particle's direction
uvw = p % uvw
uvw = p % coord0 % uvw
! determine incoming energy bin
n = t % n_bins(T_ENERGYIN)
if (n > 0) then
! 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)) cycle
! search to find incoming energy bin
bins(TS_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E)
else
bins(TS_ENERGYIN) = 1
end if
! =======================================================================
! SPECIAL CASES WHERE TWO INDICES ARE THE SAME
@ -638,16 +672,20 @@ contains
do j = ijk0(3), ijk1(3) - 1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP
call add_to_score(t % scores(score_index, 1), p % last_wgt)
bins(TS_SURFACE) = OUT_TOP
bins(1:3) = ijk0 + 1
score_index = sum((bins - 1) * t % stride) + 1
call add_to_score(t % scores(score_index, 1), p % wgt)
end if
end do
else
do j = ijk0(3) - 1, ijk1(3), -1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_TOP
call add_to_score(t % scores(score_index, 1), p % last_wgt)
bins(TS_SURFACE) = IN_TOP
bins(1:3) = ijk0 + 1
score_index = sum((bins - 1) * t % stride) + 1
call add_to_score(t % scores(score_index, 1), p % wgt)
end if
end do
end if
@ -658,16 +696,20 @@ contains
do j = ijk0(2), ijk1(2) - 1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
bins(TS_SURFACE) = OUT_FRONT
bins(1:3) = ijk0 + 1
score_index = sum((bins - 1) * t % stride) + 1
call add_to_score(t % scores(score_index, 1), p % wgt)
end if
end do
else
do j = ijk0(2) - 1, ijk1(2), -1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
bins(TS_SURFACE) = IN_FRONT
bins(1:3) = ijk0 + 1
score_index = sum((bins - 1) * t % stride) + 1
call add_to_score(t % scores(score_index, 1), p % wgt)
end if
end do
end if
@ -678,16 +720,20 @@ contains
do j = ijk0(1), ijk1(1) - 1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
bins(TS_SURFACE) = OUT_RIGHT
bins(1:3) = ijk0 + 1
score_index = sum((bins - 1) * t % stride) + 1
call add_to_score(t % scores(score_index, 1), p % wgt)
end if
end do
else
do j = ijk0(1) - 1, ijk1(1), -1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
bins(TS_SURFACE) = IN_RIGHT
bins(1:3) = ijk0 + 1
score_index = sum((bins - 1) * t % stride) + 1
call add_to_score(t % scores(score_index, 1), p % wgt)
end if
end do
end if
@ -708,7 +754,7 @@ contains
do k = 1, n_cross
! Reset scoring bin index
score_index = 0
bins(TS_SURFACE) = 0
! Calculate distance to each bounding surface. We need to treat
! special case where the cosine of the angle is zero since this would
@ -735,7 +781,8 @@ contains
! Crossing into right mesh cell -- this is treated as outgoing
! current from (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT
bins(TS_SURFACE) = OUT_RIGHT
bins(1:3) = ijk0 + 1
end if
ijk0(1) = ijk0(1) + 1
xyz_cross(1) = xyz_cross(1) + m % width(1)
@ -745,7 +792,8 @@ contains
ijk0(1) = ijk0(1) - 1
xyz_cross(1) = xyz_cross(1) - m % width(1)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT
bins(TS_SURFACE) = IN_RIGHT
bins(1:3) = ijk0 + 1
end if
end if
elseif (distance == d(2)) then
@ -753,7 +801,8 @@ contains
! Crossing into front mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT
bins(TS_SURFACE) = OUT_FRONT
bins(1:3) = ijk0 + 1
end if
ijk0(2) = ijk0(2) + 1
xyz_cross(2) = xyz_cross(2) + m % width(2)
@ -763,7 +812,8 @@ contains
ijk0(2) = ijk0(2) - 1
xyz_cross(2) = xyz_cross(2) - m % width(2)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT
bins(TS_SURFACE) = IN_FRONT
bins(1:3) = ijk0 + 1
end if
end if
else if (distance == d(3)) then
@ -771,7 +821,8 @@ contains
! Crossing into top mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP
bins(TS_SURFACE) = OUT_TOP
bins(1:3) = ijk0 + 1
end if
ijk0(3) = ijk0(3) + 1
xyz_cross(3) = xyz_cross(3) + m % width(3)
@ -781,20 +832,24 @@ contains
ijk0(3) = ijk0(3) - 1
xyz_cross(3) = xyz_cross(3) - m % width(3)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_TOP
bins(TS_SURFACE) = IN_TOP
bins(1:3) = ijk0 + 1
end if
end if
end if
! Check for errors
if (score_index < 0 .or. score_index > t % n_total_bins) then
message = "Score index outside range."
call fatal_error()
end if
! Determine scoring index
if (bins(TS_SURFACE) > 0) then
score_index = sum((bins - 1) * t % stride) + 1
! Add to surface current tally
if (score_index > 0) then
call add_to_score(t % scores(score_index, 1), p % last_wgt)
! Check for errors
if (score_index <= 0 .or. score_index > t % n_total_bins) then
message = "Score index outside range."
call fatal_error()
end if
! Add to surface current tally
call add_to_score(t % scores(score_index, 1), p % wgt)
end if
! Calculate new coordinates
@ -937,10 +992,9 @@ contains
integer :: score_index ! index in scores array for filters
logical :: file_exists ! does tallies.out file already exists?
logical :: has_filter(TALLY_TYPES) ! does tally have this filter?
character(MAX_LINE_LEN) :: filename ! name of output file
character(MAX_FILE_LEN) :: filename ! name of output file
character(15) :: filter_name(TALLY_TYPES) ! names of tally filters
character(27) :: macro_name(N_MACRO_TYPES) ! names of macro scores
character(80) :: space = " " ! spaces
type(TallyObject), pointer :: t
! Skip if there are no tallies
@ -1046,7 +1100,7 @@ contains
else
if (has_filter(j)) then
! Print current filter information
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') space(1:indent), &
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(j)), trim(get_label(t, j, bins(j)))
indent = indent + 2
end if
@ -1057,7 +1111,7 @@ contains
end do find_bin
! Print filter information
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') space(1:indent), &
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(j)), trim(get_label(t, j, bins(j)))
! Determine scoring index for this bin combination -- note that unlike
@ -1070,7 +1124,7 @@ contains
indent = indent + 2
do k = 1, t % n_macro_bins
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
space(1:indent), macro_name(abs(t % macro_bins(k) % scalar)), &
repeat(" ", indent), macro_name(abs(t % macro_bins(k) % scalar)), &
real_to_str(t % scores(score_index,k) % val), &
trim(real_to_str(t % scores(score_index,k) % val_sq))
end do
@ -1093,19 +1147,34 @@ contains
type(TallyObject), pointer :: t
integer :: i ! mesh index for x
integer :: j ! mesh index for y
integer :: k ! mesh index for z
integer :: ijk(3) ! indices in mesh
integer :: len1 ! length of string
integer :: len2 ! length of string
integer :: score_index ! index in scores array for filters
integer :: i ! mesh index for x
integer :: j ! mesh index for y
integer :: k ! mesh index for z
integer :: l ! mesh index for energy
integer :: bins(TALLY_TYPES) ! bin combination
integer :: n ! number of incoming energy bins
integer :: len1 ! length of string
integer :: len2 ! length of string
integer :: score_index ! index in scores array for filters
logical :: print_ebin ! should incoming energy bin be displayed?
character(MAX_LINE_LEN) :: string
type(StructuredMesh), pointer :: m => null()
! Get pointer to mesh
m => meshes(t % mesh)
! initialize bins array
bins = 1
! determine how many energy in bins there are
n = t % n_bins(T_ENERGYIN)
if (n > 0) then
print_ebin = .true.
else
print_ebin = .false.
n = 1
end if
do i = 1, m % dimension(1)
string = "Mesh Index (" // trim(int_to_str(i)) // ", "
len1 = len_trim(string)
@ -1117,83 +1186,113 @@ contains
string = string(1:len2+1) // trim(int_to_str(k)) // ")"
write(UNIT=UNIT_TALLY, FMT='(1X,A)') trim(string)
! Left Surface
ijk = (/ i-1, j, k /)
score_index = sum(t % stride(1:3) * ijk) + IN_RIGHT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
score_index = sum(t % stride(1:3) * ijk) + OUT_RIGHT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Right Surface
ijk = (/ i, j, k /)
score_index = sum(t % stride(1:3) * ijk) + IN_RIGHT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
score_index = sum(t % stride(1:3) * ijk) + OUT_RIGHT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Back Surface
ijk = (/ i, j-1, k /)
score_index = sum(t % stride(1:3) * ijk) + IN_FRONT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
score_index = sum(t % stride(1:3) * ijk) + OUT_FRONT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Front Surface
ijk = (/ i, j, k /)
score_index = sum(t % stride(1:3) * ijk) + IN_FRONT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
score_index = sum(t % stride(1:3) * ijk) + OUT_FRONT
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Bottom Surface
ijk = (/ i, j, k-1 /)
score_index = sum(t % stride(1:3) * ijk) + IN_TOP
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
score_index = sum(t % stride(1:3) * ijk) + OUT_TOP
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Top Surface
ijk = (/ i, j, k /)
score_index = sum(t % stride(1:3) * ijk) + IN_TOP
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
score_index = sum(t % stride(1:3) * ijk) + OUT_TOP
write(UNIT=UNIT_TALLY, FMT='(3X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
do l = 1, n
! Write incoming energy bin
if (print_ebin) then
write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') &
"Incoming Energy", trim(get_label(t, T_ENERGYIN, l))
end if
! Set incoming energy bin
bins(TS_ENERGYIN) = l
! Left Surface
bins(1:3) = (/ i-1, j, k /) + 1
bins(TS_SURFACE) = IN_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Right Surface
bins(1:3) = (/ i, j, k /) + 1
bins(TS_SURFACE) = IN_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Back Surface
bins(1:3) = (/ i, j-1, k /) + 1
bins(TS_SURFACE) = IN_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Front Surface
bins(1:3) = (/ i, j, k /) + 1
bins(TS_SURFACE) = IN_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Bottom Surface
bins(1:3) = (/ i, j, k-1 /) + 1
bins(TS_SURFACE) = IN_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
! Top Surface
bins(1:3) = (/ i, j, k /) + 1
bins(TS_SURFACE) = IN_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
end do
end do
end do
end do
@ -1212,7 +1311,7 @@ contains
integer, intent(in) :: bin ! bin in filter array
character(30) :: label ! user-specified identifier
integer :: index ! index in cells/surfaces/etc array
integer :: i ! index in cells/surfaces/etc array
integer, allocatable :: ijk(:) ! indices in mesh
real(8) :: E0 ! lower bound for energy bin
real(8) :: E1 ! upper bound for energy bin
@ -1220,20 +1319,20 @@ contains
select case(filter_type)
case (T_UNIVERSE)
index = t % universe_bins(bin) % scalar
label = int_to_str(universes(index) % id)
i = t % universe_bins(bin) % scalar
label = int_to_str(universes(i) % id)
case (T_MATERIAL)
index = t % material_bins(bin) % scalar
label = int_to_str(materials(index) % id)
i = t % material_bins(bin) % scalar
label = int_to_str(materials(i) % id)
case (T_CELL)
index = t % cell_bins(bin) % scalar
label = int_to_str(cells(index) % id)
i = t % cell_bins(bin) % scalar
label = int_to_str(cells(i) % id)
case (T_CELLBORN)
index = t % cellborn_bins(bin) % scalar
label = int_to_str(cells(index) % id)
i = t % cellborn_bins(bin) % scalar
label = int_to_str(cells(i) % id)
case (T_SURFACE)
index = t % surface_bins(bin) % scalar
label = int_to_str(surfaces(index) % id)
i = t % surface_bins(bin) % scalar
label = int_to_str(surfaces(i) % id)
case (T_MESH)
m => meshes(t % mesh)
allocate(ijk(m % n_dimension))

0
src/tally_header.f90 → src/tally_header.F90
View file

0
src/timing.f90 → src/timing.F90
View file

45
src/utils/convert_xsdir.py
View file

@ -4,6 +4,18 @@ import os
import sys
from xml.dom.minidom import getDOMImplementation
elements = [None, "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na",
"Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V",
"Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se",
"Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh",
"Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba",
"La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho",
"Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt",
"Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac",
"Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
"Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg",
"Cn"]
class Xsdir(object):
def __init__(self, filename):
@ -192,11 +204,30 @@ class XsdirTable(object):
# All other cases
A = int(self.zaid) % 1000
if A > 600:
if A > 300:
return 1
else:
return 0
@property
def alias(self):
zaid = self.zaid
Z = int(zaid[:-3])
A = zaid[-3:]
if A == '000':
s = 'Nat'
elif zaid == '95242':
s = '242m'
elif zaid == '95642':
s = '242'
elif int(A) > 300:
s = str(int(A) - 400) + "m"
else:
s = str(int(A))
return "{0}-{1}.{2}".format(elements[Z], s, self.xs)
@property
def zaid(self):
if self.name.endswith('c'):
@ -204,19 +235,17 @@ class XsdirTable(object):
else:
return 0
@property
def xs(self):
return self.name[self.name.find('.')+1:]
def to_xml_node(self, doc):
node = doc.createElement("ace_table")
node.setAttribute("name", self.name)
for attribute in ["alias", "zaid", "type", "metastable", "awr",
"temperature", "path", "location"]:
if hasattr(self, attribute):
# Join string for alias attribute
if attribute == "alias":
if not self.alias:
continue
string = " ".join(self.alias)
else:
string = str(getattr(self,attribute))
string = str(getattr(self,attribute))
# Skip metastable and binary if 0
if attribute == "metastable" and self.metastable == 0:

2
src/xml-fortran/templates/geometry_t.xml
View file

@ -22,7 +22,7 @@
<typedef name="lattice_xml">
<component name="id" type="integer" />
<component name="type" type="word" length="12" />
<component name="type" type="word" length="12" default="'rectangular'" />
<component name="dimension" type="integer-array" />
<component name="origin" type="double-array" />
<component name="width" type="double-array" />

2
src/xml-fortran/templates/materials_t.xml
View file

@ -15,7 +15,7 @@
<!-- Type for specifying a nuclide -->
<typedef name="nuclide_xml">
<component name="name" type="word" length="10" default="''" />
<component name="name" type="word" length="7" default="''" />
<component name="xs" type="word" length="3" default="''" />
<component name="ao" type="double" default="0.0" />
<component name="wo" type="double" default="0.0" />

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

@ -25,5 +25,6 @@
<variable name="survival_" tag="survival_biasing" 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="trace_" tag="trace" type="integer-array" />
</template>

1
src/xml-fortran/xmlreader.f90
View file

@ -552,6 +552,7 @@ subroutine write_prolog
& ' use WRITE_XML_PRIMITIVES', &
& ' use XMLPARSE', &
& ' implicit none', &
& ' save', &
& ' integer, private :: lurep_', &
& ' logical, private :: strict_'