diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000000..dd137b950e
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,27 @@
+before_install:
+ - sudo apt-get update -qq
+ - sudo apt-get install -qq -y python-numpy
+ - sudo apt-get install -qq -y python-scipy
+ - sudo apt-get install -qq -y python-h5py
+ - sudo apt-get install -qq -y gfortran
+ - sudo apt-get install -qq -y g++
+ - ./tests/travis_install.sh
+ - export FC=gfortran
+ - export MPI_DIR=$PWD/mpich_install
+ - export PHDF5_DIR=$PWD/phdf5_install
+ - export HDF5_DIR=$PWD/hdf5_install
+ - export PETSC_DIR=$PWD/petsc_install
+
+before_script:
+ - cd data
+ - git clone --branch=master git://github.com/bhermanmit/nndc_xs nndc_xs
+ - cat nndc_xs/nndc.tar.gza* | tar xzvf -
+ - rm -rf nndc_xs
+ - export CROSS_SECTIONS=$PWD/nndc/cross_sections.xml
+ - cd ..
+
+script:
+ - cd tests
+ - export OMP_NUM_THREADS=3
+ - ./travis.sh
+ - cd ..
diff --git a/data/get_nndc_data.py b/data/get_nndc_data.py
index 12a6de6e0b..5df69f3ac2 100755
--- a/data/get_nndc_data.py
+++ b/data/get_nndc_data.py
@@ -8,6 +8,12 @@ import sys
import tarfile
import glob
import hashlib
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument('-b', '--batch', action = 'store_true',
+ help = 'supresses standard in')
+args = parser.parse_args()
try:
from urllib.request import urlopen
@@ -73,7 +79,7 @@ for f in files:
print('Verifying MD5 checksums...')
for f, checksum in zip(files, checksums):
- downloadsum = hashlib.md5(open(f, 'r').read()).hexdigest()
+ downloadsum = hashlib.md5(open(f, 'rb').read()).hexdigest()
if downloadsum != checksum:
raise IOError("MD5 checksum for {} does not match. If this is your first "
"time receiving this message, please re-run the script. "
@@ -118,10 +124,13 @@ shutil.copyfile('cross_sections_nndc.xml', 'nndc/cross_sections.xml')
# PROMPT USER TO DELETE .TAR.GZ FILES
# Ask user to delete
-if sys.version_info[0] < 3:
- response = raw_input('Delete *.tar.gz files? ([y]/n) ')
+if not args.batch:
+ if sys.version_info[0] < 3:
+ response = raw_input('Delete *.tar.gz files? ([y]/n) ')
+ else:
+ response = input('Delete *.tar.gz files? ([y]/n) ')
else:
- response = input('Delete *.tar.gz files? ([y]/n) ')
+ response = 'y'
# Delete files if requested
if not response or response.lower().startswith('y'):
@@ -134,10 +143,13 @@ if not response or response.lower().startswith('y'):
# PROMPT USER TO CONVERT ASCII TO BINARY
# Ask user to convert
-if sys.version_info[0] < 3:
- response = raw_input('Convert ACE files to binary? ([y]/n) ')
+if not args.batch:
+ if sys.version_info[0] < 3:
+ response = raw_input('Convert ACE files to binary? ([y]/n) ')
+ else:
+ response = input('Convert ACE files to binary? ([y]/n) ')
else:
- response = input('Convert ACE files to binary? ([y]/n) ')
+ response = 'y'
# Convert files if requested
if not response or response.lower().startswith('y'):
diff --git a/docs/source/_images/hex_lat.svg b/docs/source/_images/hex_lat.svg
new file mode 100644
index 0000000000..cc150e711e
--- /dev/null
+++ b/docs/source/_images/hex_lat.svg
@@ -0,0 +1,780 @@
+
+
+
+
diff --git a/docs/source/_images/rect_lat.svg b/docs/source/_images/rect_lat.svg
new file mode 100644
index 0000000000..cec6260c63
--- /dev/null
+++ b/docs/source/_images/rect_lat.svg
@@ -0,0 +1,711 @@
+
+
+
+
diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst
index 1ba29b79ef..8dda10e013 100644
--- a/docs/source/methods/geometry.rst
+++ b/docs/source/methods/geometry.rst
@@ -385,6 +385,101 @@ is found that contains the specified point.
.. _cell-contains:
+----------------------
+Finding a Lattice Tile
+----------------------
+
+If a particle is inside a lattice, its position inside the lattice must be
+determined before assigning it to a cell. Throughout this section, the
+volumetric units of the lattice will be referred to as "tiles". Tiles are
+identified by thier indices, and the process of discovering which tile contains
+the particle is referred to as "indexing".
+
+Rectilinear Lattice Indexing
+----------------------------
+
+Indices are assigned to tiles in a rectilinear lattice based on the tile's
+position along the :math:`x`, :math:`y`, and :math:`z` axes. The figure below
+maps the indices for a 2D lattice. The indices, (1, 1), map to the
+lower-left tile. (5, 1) and (5, 5) map to the lower-right and upper-right
+tiles, respectively.
+
+.. figure:: ../_images/rect_lat.*
+ :align: center
+ :figclass: align-center
+ :width: 400px
+
+ Rectilinear lattice tile indices.
+
+In general, a lattice tile is specified by the three indices,
+:math:`(i_x, i_y, i_z)`. If a particle's current coordinates are
+:math:`(x, y, z)` then the indices can be determined from these formulas:
+
+.. math::
+ :label: rect_indexing
+
+ i_x = \left \lceil \frac{x - x_0}{p_0} \right \rceil
+
+ i_y = \left \lceil \frac{y - y_0}{p_1} \right \rceil
+
+ i_z = \left \lceil \frac{z - z_0}{p_2} \right \rceil
+
+where :math:`(x_0, y_0, z_0)` are the coordinates to the lower-left-bottom
+corner of the lattice, and :math:`p_0, p_1, p_2` are the pitches along the
+:math:`x`, :math:`y`, and :math:`z` axes, respectively.
+
+Hexagonal Lattice Indexing
+--------------------------
+
+A skewed coordinate system is used for indexing hexagonal lattice tiles. Rather
+than a :math:`y`-axis, another axis is used that is rotated 30 degrees
+counter-clockwise from the :math:`y`-axis. This axis is referred to as the
+:math:`\alpha`-axis. The figure below shows how 2D hexagonal tiles are mapped
+with the :math:`(x, \alpha)` basis. In this system, (0, 0) maps to the center
+tile, (0, 2) to the top tile, and (2, -1) to the middle tile on the right side.
+
+.. figure:: ../_images/hex_lat.*
+ :align: center
+ :figclass: align-center
+ :width: 400px
+
+ Hexagonal lattice tile indices.
+
+Unfortunately, the indices cannot be determined with one simple formula as
+before. Indexing requires a two-step process, a coarse step which determines a
+set of four tiles that contains the particle and a fine step that determines
+which of those four tiles actually contains the particle.
+
+In the first step, indices are found using these formulas:
+
+.. math::
+ :label: hex_indexing
+
+ \alpha = -\frac{x}{\sqrt{3}} + y
+
+ i_x^* = \left \lfloor \frac{x}{p_0 \sqrt{3} / 2} \right \rfloor
+
+ i_\alpha^* = \left \lfloor \frac{\alpha}{p_0} \right \rfloor
+
+where :math:`p_0` is the lattice pitch (in the :math:`x`-:math:`y` plane). The
+true index of the particle could be :math:`(i_x^*, i_\alpha^*)`,
+:math:`(i_x^* + 1, i_\alpha^*)`, :math:`(i_x^*, i_\alpha^* + 1)`, or
+:math:`(i_x^* + 1, i_\alpha^* + 1)`.
+
+The second step selects the correct tile from that neighborhood of 4. OpenMC
+does this by calculating the distance between the particle and the centers of
+each of the 4 tiles, and then picking the closest tile. This works because
+regular hexagonal tiles form a Voronoi tessellation which means that all of the
+points within a tile are closest to the center of that same tile.
+
+Indexing along the :math:`z`-axis uses the same method from rectilinear
+lattices, i.e.
+
+.. math::
+ :label: hex_indexing_z
+
+ i_z = \left \lceil \frac{z - z_0}{p_2} \right \rceil
+
----------------------------------------
Determining if a Coordinate is in a Cell
----------------------------------------
diff --git a/docs/source/usersguide/beginners.rst b/docs/source/usersguide/beginners.rst
index 69f44c2ca1..d65ee1f831 100644
--- a/docs/source/usersguide/beginners.rst
+++ b/docs/source/usersguide/beginners.rst
@@ -146,13 +146,12 @@ and `Volume II`_. You may also find it helpful to review the following terms:
.. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
.. _git: http://git-scm.com/
.. _git tutorials: http://git-scm.com/documentation
-.. _Reactor Concepts Manual: http://web.mit.edu/romano7/www/reactor_concepts.pdf
-.. _Volume I: http://www.hss.doe.gov/nuclearsafety/techstds/docs/handbook/h1019v1.pdf
-.. _Volume II: http://www.hss.doe.gov/nuclearsafety/techstds/docs/handbook/h1019v2.pdf
+.. _Reactor Concepts Manual: http://www.tayloredge.com/periodic/trivia/ReactorConcepts.pdf
+.. _Volume I: http://energy.gov/sites/prod/files/2013/06/f2/h1019v1.pdf
+.. _Volume II: http://energy.gov/sites/prod/files/2013/06/f2/h1019v2.pdf
.. _OpenMC source code: https://github.com/mit-crpg/openmc
.. _GitHub: https://github.com/
.. _bug reports: https://github.com/mit-crpg/openmc/issues
.. _Neutron cross section: http://en.wikipedia.org/wiki/Neutron_cross_section
.. _Effective multiplication factor: http://en.wikipedia.org/wiki/Effective_multiplication_factor
.. _Flux: http://en.wikipedia.org/wiki/Neutron_flux
-
diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst
index 20fd66e2f7..cb1a577039 100644
--- a/docs/source/usersguide/input.rst
+++ b/docs/source/usersguide/input.rst
@@ -51,6 +51,55 @@ files are called:
* ``plots.xml``
* ``cmfd.xml``
+--------------------
+Validating XML Files
+--------------------
+
+Input files can be checked before executing OpenMC using the ``xml_validate``
+script. It is located in ``src/utils/xml_validate.py`` in the source code or in
+``bin/xml_validate`` in the install directory. Before use, the third party
+package TRANG_ must be installed and in your ``PATH`` to convert the compact
+RelaxNG schema to standard RelaxNG format. For Ubuntu, you can install with:
+
+.. code-block:: bash
+
+ sudo apt-get install trang
+
+Two command line arguments can be set when running ``xml_validate``:
+
+* ``-i``, ``--input-path`` - Location of OpenMC input files.
+ *Default*: current working directory
+* ``-r``, ``--relaxng-path`` - Location of OpenMC RelaxNG files.
+ *Default*: None
+
+If the RelaxNG path is not set, ``xml_validate`` will search for these files
+because it expects that the user is either running the script located in the
+install directory ``bin`` folder or in ``src/utils``. Once executed, it will
+match OpenMC XML files with their RelaxNG schema and check if they are valid.
+Below is a table of the messages that will be printed after each file is
+checked.
+
+======================== ===================================
+Message Description
+======================== ===================================
+[XML ERROR] Cannot parse XML file.
+[NO RELAXNG FOUND] No RelaxNG file found for XML file.
+[TRANG FAILED] TRANG not installed properly.
+[NOT VALID] XML file does not match RelaxNG.
+[VALID] XML file matches RelaxNG.
+======================== ===================================
+
+As an example, if OpenMC is installed in the directory
+``/opt/openmc/0.6.2`` and the current working directory is where
+OpenMC XML input files are located, they can be validated using
+the following command:
+
+.. code-block:: bash
+
+ /opt/openmc/0.6.2/bin/xml_validate
+
+.. _TRANG: http://www.thaiopensource.com/relaxng/trang.html
+
--------------------------------------
Settings Specification -- settings.xml
--------------------------------------
@@ -802,19 +851,12 @@ Each ```` element can have the following attributes or sub-elements:
---------------------
The ```` can be used to represent repeating structures (e.g. fuel pins
-in an assembly) or other geometry which naturally fits into a two- or
-three-dimensional structured mesh. Each cell within the lattice is filled with a
-specified universe. A ```` accepts the following attributes or
-sub-elements:
+in an assembly) or other geometry which fits onto a rectilinear grid. Each cell
+within the lattice is filled with a specified universe. A ```` accepts
+the following attributes or sub-elements:
:id:
- A unique integer that can be used to identify the surface.
-
- :type:
- A string indicating the arrangement of lattice cells. Currently, the only
- accepted option is "rectangular".
-
- *Default*: rectangular
+ A unique integer that can be used to identify the lattice.
:dimension:
Two or three integers representing the number of lattice cells in the x- and
@@ -828,8 +870,10 @@ sub-elements:
*Default*: None
- :width:
- The width of the lattice cell in the x- and y- (and z-) directions.
+ :pitch:
+ If the lattice is 3D, then three real numbers that express the distance
+ between the centers of lattice cells in the x-, y-, and z- directions. If
+ the lattice is 2D, then omit the third value.
*Default*: None
@@ -846,6 +890,92 @@ sub-elements:
*Default*: None
+Here is an example of a properly defined 2d rectangular lattice:
+
+.. code-block:: xml
+
+
+ -1.5 -1.5
+ 1.0 1.0
+
+ 2 2 2
+ 2 1 2
+ 2 2 2
+
+
+
+```` Element
+-------------------------
+
+The ```` can be used to represent repeating structures (e.g. fuel
+pins in an assembly) or other geometry which naturally fits onto a hexagonal
+grid or hexagonal prism grid. Each cell within the lattice is filled with a
+specified universe. This lattice uses the "flat-topped hexagon" scheme where two
+of the six edges are perpendicular to the y-axis. A ```` accepts
+the following attributes or sub-elements:
+
+ :id:
+ A unique integer that can be used to identify the lattice.
+
+ :n_rings:
+ An integer representing the number of radial ring positions in the xy-plane.
+ Note that this number includes the degenerate center ring which only has one
+ element.
+
+ *Default*: None
+
+ :n_axial:
+ An integer representing the number of positions along the z-axis. This
+ element is optional.
+
+ *Default*: None
+
+ :center:
+ The coordinates of the center of the lattice. If the lattice does not have
+ axial sections then only the x- and y-coordinates are specified.
+
+ *Default*: None
+
+ :pitch:
+ If the lattice is 3D, then two real numbers that express the distance
+ between the centers of lattice cells in the xy-plane and along the z-axis,
+ respectively. If the lattice is 2D, then omit the second value.
+
+ *Default*: None
+
+ :outer:
+ The unique integer identifier of a universe that will be used to fill all
+ space outside of the lattice. The universe will be tiled repeatedly as if
+ it were placed in a lattice of infinite size. This element is optional.
+
+ *Default*: An error will be raised if a particle leaves a lattice with no
+ outer universe.
+
+ :universes:
+ A list of the universe numbers that fill each cell of the lattice.
+
+ *Default*: None
+
+Here is an example of a properly defined 2d hexagonal lattice:
+
+.. code-block:: xml
+
+
+
0.0 0.0
+ 1.0
+
+ 202
+ 202 202
+ 202 202 202
+ 202 202
+ 202 101 202
+ 202 202
+ 202 202 202
+ 202 202
+ 202
+
+
+
.. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
.. _quadratic surfaces: http://en.wikipedia.org/wiki/Quadric
diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst
index 7cb7ac7a45..6fc002b1b8 100644
--- a/docs/source/usersguide/install.rst
+++ b/docs/source/usersguide/install.rst
@@ -88,11 +88,11 @@ Prerequisites
* PETSc_ for CMFD acceleration
- To enable CMFD acceleration, you will need to have PETSc_ (3.4.2 or higher)
- installed on your computer. The installed version will need to have been
- compiled with the same compiler you intend to compile OpenMC with. OpenMC
- requires PETSc_ to be configured with Fortran datatypes. An example of
- configuring PETSc_ is listed below::
+ To enable some features of CMFD acceleration, you will need to have
+ PETSc_ (3.4.2 or higher) installed on your computer. The installed version
+ will need to have been compiled with the same compiler you intend to
+ compile OpenMC with. OpenMC requires PETSc_ to be configured with Fortran
+ datatypes. An example of configuring PETSc_ is listed below::
./configure --prefix=/opt/petsc/3.4.4 --download-f-blas-lapack \
--with-mpi-dir=/opt/mpich/3.1 --with-shared-libraries \
diff --git a/examples/xml/lattice/nested/geometry.xml b/examples/xml/lattice/nested/geometry.xml
index f223cf3bb4..9df9b9e931 100644
--- a/examples/xml/lattice/nested/geometry.xml
+++ b/examples/xml/lattice/nested/geometry.xml
@@ -12,10 +12,9 @@
- rectangular2 2-1.0 -1.0
- 1.0 1.0
+ 1.0 1.0
1 2
2 3
@@ -24,10 +23,9 @@
- rectangular2 2-2.0 -2.0
- 2.0 2.0
+ 2.0 2.0
5 5
5 5
diff --git a/examples/xml/lattice/simple/geometry.xml b/examples/xml/lattice/simple/geometry.xml
index b5bad90b9c..c1f8d78644 100644
--- a/examples/xml/lattice/simple/geometry.xml
+++ b/examples/xml/lattice/simple/geometry.xml
@@ -10,10 +10,9 @@
- rectangular4 4-2.0 -2.0
- 1.0 1.0
+ 1.0 1.0
1 2 1 2
2 3 2 3
diff --git a/readme.rst b/readme.rst
index a010974375..ab13e36f5c 100644
--- a/readme.rst
+++ b/readme.rst
@@ -2,6 +2,9 @@
OpenMC Monte Carlo Particle Transport Code
==========================================
+.. image:: https://travis-ci.org/mit-crpg/openmc.svg?branch=develop
+ :target: https://travis-ci.org/mit-crpg/openmc
+
The OpenMC project aims to provide a fully-featured Monte Carlo particle
transport code based on modern methods. It is a constructive solid geometry,
continuous-energy transport code that uses ACE format cross sections. The
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index dd920bfcc0..b3e1a88779 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -272,6 +272,10 @@ install(PROGRAMS utils/statepoint_meshplot.py
install(PROGRAMS utils/update_inputs.py
DESTINATION bin
RENAME update_inputs)
+install(PROGRAMS utils/xml_validate.py
+ DESTINATION bin
+ RENAME xml_validate)
+install(DIRECTORY relaxng DESTINATION share)
install(FILES ../man/man1/openmc.1 DESTINATION share/man/man1)
install(FILES ../LICENSE DESTINATION "share/doc/${program}/copyright")
diff --git a/src/geometry.F90 b/src/geometry.F90
index 3aeaf3afc4..c4837e1057 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -2,7 +2,8 @@ module geometry
use constants
use error, only: fatal_error, warning
- use geometry_header, only: Cell, Surface, Universe, Lattice
+ use geometry_header, only: Cell, Surface, Universe, Lattice, &
+ &RectLattice, HexLattice
use global
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord, Particle
@@ -131,18 +132,13 @@ contains
logical, intent(inout) :: found
integer, optional :: search_cells(:)
- integer :: i ! index over cells
- integer :: i_x, i_y, i_z ! indices in lattice
- integer :: n_x, n_y, n_z ! size of lattice
- integer :: n ! number of cells to search
- integer :: index_cell ! index in cells array
- real(8) :: xyz(3) ! temporary location
- real(8) :: upper_right(3) ! lattice upper_right
- logical :: use_search_cells ! use cells provided as argument
- logical :: outside_lattice ! if particle is not inside lattice bounds
- logical :: lattice_edge ! if particle is on a lattice edge
+ integer :: i ! index over cells
+ integer :: i_xyz(3) ! indices in lattice
+ integer :: n ! number of cells to search
+ integer :: index_cell ! index in cells array
+ logical :: use_search_cells ! use cells provided as argument
type(Cell), pointer, save :: c => null() ! pointer to cell
- type(Lattice), pointer, save :: lat => null() ! pointer to lattice
+ class(Lattice), pointer, save :: lat => null() ! pointer to lattice
type(Universe), pointer, save :: univ => null() ! universe to search in
!$omp threadprivate(c, lat, univ)
@@ -159,7 +155,7 @@ contains
n = univ % n_cells
end if
- do i = 1, n
+ CELL_LOOP: do i = 1, n
! select cells based on whether we are searching a universe or a provided
! list of cells (this would be for lists of neighbor cells)
if (use_search_cells) then
@@ -173,170 +169,102 @@ contains
! get pointer to cell
c => cells(index_cell)
- if (simple_cell_contains(c, p)) then
- ! Set cell on this level
- p % coord % cell = index_cell
+ ! Move on to the next cell if the particle is not inside this cell
+ if (.not. simple_cell_contains(c, p)) cycle
- ! Show cell information on trace
- if (verbosity >= 10 .or. trace) then
- call write_message(" Entering cell " // trim(to_str(c % id)))
- end if
+ ! Set cell on this level
+ p % coord % cell = index_cell
- if (c % type == CELL_NORMAL) then
- ! ====================================================================
- ! AT LOWEST UNIVERSE, TERMINATE SEARCH
-
- ! set material
- p % last_material = p % material
- p % material = c % material
-
- elseif (c % type == CELL_FILL) then
- ! ====================================================================
- ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-
- ! Create new level of coordinates
- allocate(p % coord % next)
- p % coord % next % xyz = p % coord % xyz
- p % coord % next % uvw = p % coord % uvw
-
- ! Move particle to next level and set universe
- p % coord => p % coord % next
- p % coord % universe = c % fill
-
- ! Apply translation
- if (allocated(c % translation)) then
- p % coord % xyz = p % coord % xyz - c % translation
- end if
-
- ! Apply rotation
- if (allocated(c % rotation_matrix)) then
- p % coord % xyz = matmul(c % rotation_matrix, p % coord % xyz)
- p % coord % uvw = matmul(c % rotation_matrix, p % coord % uvw)
- p % coord % rotated = .true.
- end if
-
- call find_cell(p, found)
- if (.not. found) exit
-
- elseif (c % type == CELL_LATTICE) then
- ! ====================================================================
- ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-
- ! Set current lattice
- lat => lattices(c % fill)
-
- outside_lattice = .false.
- lattice_edge = .false.
-
- ! determine lattice index based on position
- xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
- i_x = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
- i_y = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
- n_x = lat % dimension(1)
- n_y = lat % dimension(2)
- if (lat % n_dimension == 3) then
- i_z = ceiling((xyz(3) - lat % lower_left(3))/lat % width(3))
- n_z = lat % dimension(3)
- else
- i_z = 1
- n_z = 1
- end if
-
- ! Check if lattice coordinates are within bounds
- if (i_x < 1 .or. i_x > n_x .or. i_y < 1 .or. i_y > n_y .or. &
- i_z < 1 .or. i_z > n_z) then
-
- ! Check for when particle is on lattice edge
- upper_right(1) = lat % lower_left(1) + &
- lat % width(1) * dble(lat % dimension(1))
- upper_right(2) = lat % lower_left(2) + &
- lat % width(2) * dble(lat % dimension(2))
- if ( abs(xyz(1) - lat % lower_left(1)) < FP_COINCIDENT .or. &
- abs(xyz(2) - lat % lower_left(2)) < FP_COINCIDENT .or. &
- abs(upper_right(1) - xyz(1)) < FP_COINCIDENT .or. &
- abs(upper_right(2) - xyz(2)) < FP_COINCIDENT) then
- lattice_edge = .true.
- end if
- if (lat % n_dimension == 3) then
- upper_right(3) = lat % lower_left(3) + &
- lat % width(3) * dble(lat % dimension(3))
- if (abs(xyz(3) - lat % lower_left(3)) < FP_COINCIDENT .or. &
- abs(upper_right(3) - xyz(3)) < FP_COINCIDENT) then
- lattice_edge = .true.
- end if
- end if
-
- if (lattice_edge) then
-
- ! In this case the neutron is leaving the lattice, so we move it
- ! out, remove all lower coordinate levels and then search from
- ! universe 0.
-
- p % coord => p % coord0
- call deallocate_coord(p % coord % next)
-
- ! Reset surface and advance particle a tiny bit
- p % surface = NONE
- p % coord % xyz = xyz
-
- else
- outside_lattice = .true.
-
- end if
-
- end if
-
- if (.not. lattice_edge) then
-
- ! Create new level of coordinates
- allocate(p % coord % next)
-
- ! adjust local position of particle
- p % coord % next % xyz(1) = p % coord % xyz(1) - &
- (lat % lower_left(1) + (i_x - 0.5_8)*lat % width(1))
- p % coord % next % xyz(2) = p % coord % xyz(2) - &
- (lat % lower_left(2) + (i_y - 0.5_8)*lat % width(2))
- if (lat % n_dimension == 3) then
- p % coord % next % xyz(3) = p % coord % xyz(3) - &
- (lat % lower_left(3) + (i_z - 0.5_8)*lat % width(3))
- else
- p % coord % next % xyz(3) = p % coord % xyz(3)
- end if
- p % coord % next % uvw = p % coord % uvw
-
- ! set particle lattice indices
- p % coord % next% lattice = c % fill
- p % coord % next% lattice_x = i_x
- p % coord % next% lattice_y = i_y
- p % coord % next% lattice_z = i_z
- if (.not. outside_lattice) then
- p % coord % next % universe = lat % universes(i_x,i_y,i_z)
- else
- if (lat % outer == NO_OUTER_UNIVERSE) then
- call fatal_error("A particle is outside latttice " &
- &// trim(to_str(lat % id)) // " but the lattice has no &
- &defined outer universe.")
- else
- p % coord % next % universe = lat % outer
- end if
- end if
-
- ! Move particle to next level
- p % coord => p % coord % next
-
- end if
-
- ! Find in the next lowest coordinate level.
- call find_cell(p, found)
- if (.not. found) exit
-
- end if
-
- ! Found cell so we can return
- found = .true.
- return
+ ! Show cell information on trace
+ if (verbosity >= 10 .or. trace) then
+ call write_message(" Entering cell " // trim(to_str(c % id)))
end if
- end do
+
+ CELL_TYPE: if (c % type == CELL_NORMAL) then
+ ! ======================================================================
+ ! AT LOWEST UNIVERSE, TERMINATE SEARCH
+
+ ! set material
+ p % last_material = p % material
+ p % material = c % material
+
+ elseif (c % type == CELL_FILL) then CELL_TYPE
+ ! ======================================================================
+ ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
+
+ ! Create new level of coordinates
+ allocate(p % coord % next)
+ p % coord % next % xyz = p % coord % xyz
+ p % coord % next % uvw = p % coord % uvw
+
+ ! Move particle to next level and set universe
+ p % coord => p % coord % next
+ p % coord % universe = c % fill
+
+ ! Apply translation
+ if (allocated(c % translation)) then
+ p % coord % xyz = p % coord % xyz - c % translation
+ end if
+
+ ! Apply rotation
+ if (allocated(c % rotation_matrix)) then
+ p % coord % xyz = matmul(c % rotation_matrix, p % coord % xyz)
+ p % coord % uvw = matmul(c % rotation_matrix, p % coord % uvw)
+ p % coord % rotated = .true.
+ end if
+
+ call find_cell(p, found)
+ if (.not. found) exit
+
+ elseif (c % type == CELL_LATTICE) then CELL_TYPE
+ ! ======================================================================
+ ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
+
+ ! Set current lattice
+ lat => lattices(c % fill) % obj
+
+ ! Determine lattice indices
+ i_xyz = lat % get_indices(p % coord % xyz + TINY_BIT * p % coord % uvw)
+
+ ! Create new level of coordinates
+ allocate(p % coord % next)
+ p % coord % next % xyz = lat % get_local_xyz(p % coord % xyz, i_xyz)
+ p % coord % next % uvw = p % coord % uvw
+
+ ! set particle lattice indices
+ p % coord % next% lattice = c % fill
+ p % coord % next% lattice_x = i_xyz(1)
+ p % coord % next% lattice_y = i_xyz(2)
+ p % coord % next% lattice_z = i_xyz(3)
+
+ ! Set the next lowest coordinate level.
+ if (lat % are_valid_indices(i_xyz)) then
+ ! Particle is inside the lattice.
+ p % coord % next % universe = &
+ &lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
+
+ else
+ ! Particle is outside the lattice.
+ if (lat % outer == NO_OUTER_UNIVERSE) then
+ call fatal_error("A particle is outside latttice " &
+ &// trim(to_str(lat % id)) // " but the lattice has no &
+ &defined outer universe.")
+ else
+ p % coord % next % universe = lat % outer
+ end if
+ end if
+
+ ! Move particle to next level and search for the lower cells.
+ p % coord => p % coord % next
+ call find_cell(p, found)
+ if (.not. found) exit
+
+ end if CELL_TYPE
+
+ ! Found cell so we can return
+ found = .true.
+ return
+ end do CELL_LOOP
found = .false.
@@ -640,19 +568,18 @@ contains
! CROSS_LATTICE moves a particle into a new lattice element
!===============================================================================
- subroutine cross_lattice(p, lattice_crossed)
+ subroutine cross_lattice(p, lattice_translation)
type(Particle), intent(inout) :: p
- integer, intent(in) :: lattice_crossed
+ integer, intent(in) :: lattice_translation(3)
- integer :: i_x, i_y, i_z ! indices in lattice
- integer :: n_x, n_y, n_z ! size of lattice
- real(8) :: x0, y0, z0 ! half width of lattice element
- logical :: found ! particle found in cell?
- type(Lattice), pointer, save :: lat => null()
-!$omp threadprivate(lat)
+ integer :: i_xyz(3) ! indices in lattice
+ logical :: found ! particle found in cell?
+ class(Lattice), pointer, save :: lat => null()
+ type(LocalCoord), pointer, save :: parent_coord => null()
+!$omp threadprivate(lat, parent_coord)
- lat => lattices(p % coord % lattice)
+ lat => lattices(p % coord % lattice) % obj
if (verbosity >= 10 .or. trace) then
call write_message(" Crossing lattice " // trim(to_str(lat % id)) &
@@ -661,80 +588,42 @@ contains
&// trim(to_str(p % coord % lattice_z)) // ")")
end if
- if (lat % type == LATTICE_RECT) then
- x0 = lat % width(1) * 0.5_8
- y0 = lat % width(2) * 0.5_8
- if (lat % n_dimension == 3) z0 = lat % width(3) * 0.5_8
+ ! Find the coordiante level just above the current one.
+ parent_coord => p % coord0
+ do while(.not. associated(parent_coord % next, p % coord))
+ parent_coord => parent_coord % next
+ end do
- select case (lattice_crossed)
- case (LATTICE_LEFT)
- ! Move particle to left element
- p % coord % lattice_x = p % coord % lattice_x - 1
- p % coord % xyz(1) = x0
+ ! Set the lattice indices.
+ p % coord % lattice_x = p % coord % lattice_x + lattice_translation(1)
+ p % coord % lattice_y = p % coord % lattice_y + lattice_translation(2)
+ p % coord % lattice_z = p % coord % lattice_z + lattice_translation(3)
+ i_xyz(1) = p % coord % lattice_x
+ i_xyz(2) = p % coord % lattice_y
+ i_xyz(3) = p % coord % lattice_z
- case (LATTICE_RIGHT)
- ! Move particle to right element
- p % coord % lattice_x = p % coord % lattice_x + 1
- p % coord % xyz(1) = -x0
+ ! Set the new coordinate position.
+ p % coord % xyz = lat % get_local_xyz(parent_coord % xyz, i_xyz)
- case (LATTICE_BACK)
- ! Move particle to bottom element
- p % coord % lattice_y = p % coord % lattice_y - 1
- p % coord % xyz(2) = y0
-
- case (LATTICE_FRONT)
- ! Move particle to top element
- p % coord % lattice_y = p % coord % lattice_y + 1
- p % coord % xyz(2) = -y0
-
- case (LATTICE_BOTTOM)
- ! Move particle to bottom element
- p % coord % lattice_z = p % coord % lattice_z - 1
- p % coord % xyz(3) = z0
-
- case (LATTICE_TOP)
- ! Move particle to top element
- p % coord % lattice_z = p % coord % lattice_z + 1
- p % coord % xyz(3) = -z0
-
- end select
- elseif (lat % type == LATTICE_HEX) then
- ! TODO: Add hex lattice support
- end if
-
- ! Check to make sure still in lattice
- i_x = p % coord % lattice_x
- i_y = p % coord % lattice_y
- i_z = p % coord % lattice_z
- n_x = lat % dimension(1)
- n_y = lat % dimension(2)
- if (lat % n_dimension == 3) then
- n_z = lat % dimension(3)
- else
- n_z = 1
- end if
- if (i_x < 1 .or. i_x > n_x .or. i_y < 1 .or. i_y > n_y .or. &
- i_z < 1 .or. i_z > n_z) then
+ OUTSIDE_LAT: if (.not. lat % are_valid_indices(i_xyz)) then
+ ! The particle is outside the lattice. Search for it from coord0.
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
-
- ! Search for particle
call find_cell(p, found)
if (.not. found) then
call handle_lost_particle(p, "Could not locate particle " &
&// trim(to_str(p % id)) // " after crossing a lattice boundary.")
return
end if
- else
- ! Find universe for next lattice element
- p % coord % universe = lat % universes(i_x, i_y, i_z)
+ else OUTSIDE_LAT
! Find cell in next lattice element
+ p % coord % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
call find_cell(p, found)
if (.not. found) then
- ! In some circumstances, a particle crossing the corner of a cell may not
- ! be able to be found in the next universe. In this scenario we cut off
- ! all lower-level coordinates and search from universe zero
+ ! In some circumstances, a particle crossing the corner of a cell may
+ ! not be able to be found in the next universe. In this scenario we cut
+ ! off all lower-level coordinates and search from universe zero
! Remove lower coordinates
call deallocate_coord(p % coord0 % next)
@@ -749,7 +638,7 @@ contains
return
end if
end if
- end if
+ end if OUTSIDE_LAT
end subroutine cross_lattice
@@ -759,34 +648,47 @@ contains
! that has a parent cell, also include the surfaces of the edge of the universe.
!===============================================================================
- subroutine distance_to_boundary(p, dist, surface_crossed, lattice_crossed)
+ subroutine distance_to_boundary(p, dist, surface_crossed, lattice_translation)
type(Particle), intent(inout) :: p
real(8), intent(out) :: dist
integer, intent(out) :: surface_crossed
- integer, intent(out) :: lattice_crossed
+ integer, intent(out) :: lattice_translation(3)
- integer :: i ! index for surface in cell
- integer :: index_surf ! index in surfaces array (with sign)
- real(8) :: x,y,z ! particle coordinates
- real(8) :: u,v,w ! particle directions
- real(8) :: d ! evaluated distance
- real(8) :: x0,y0,z0 ! coefficients for surface
- real(8) :: r ! radius for quadratic surfaces
- real(8) :: tmp ! dot product of surface normal with direction
- real(8) :: a,b,c,k ! quadratic equation coefficients
- real(8) :: quad ! discriminant of quadratic equation
- logical :: on_surface ! is particle on surface?
+ integer :: i ! index for surface in cell
+ integer :: index_surf ! index in surfaces array (with sign)
+ integer :: i_xyz(3) ! lattice indices
+ integer :: level_surf_cross ! surface crossed on current level
+ integer :: level_lat_trans(3) ! lattice translation on current level
+ real(8) :: x,y,z ! particle coordinates
+ real(8) :: xyz_t(3) ! local particle coordinates
+ real(8) :: beta, gama ! skewed particle coordiantes
+ real(8) :: u,v,w ! particle directions
+ real(8) :: beta_dir ! skewed particle direction
+ real(8) :: gama_dir ! skewed particle direction
+ real(8) :: edge ! distance to oncoming edge
+ real(8) :: d ! evaluated distance
+ real(8) :: d_lat ! distance to lattice boundary
+ real(8) :: d_surf ! distance to surface
+ real(8) :: x0,y0,z0 ! coefficients for surface
+ real(8) :: r ! radius for quadratic surfaces
+ real(8) :: tmp ! dot product of surface normal with direction
+ real(8) :: a,b,c,k ! quadratic equation coefficients
+ real(8) :: quad ! discriminant of quadratic equation
+ logical :: on_surface ! is particle on surface?
type(Cell), pointer, save :: cl => null()
type(Surface), pointer, save :: surf => null()
- type(Lattice), pointer, save :: lat => null()
+ class(Lattice), pointer, save :: lat => null()
type(LocalCoord), pointer, save :: coord => null()
type(LocalCoord), pointer, save :: final_coord => null()
-!$omp threadprivate(cl, surf, lat, coord, final_coord)
+ type(LocalCoord), pointer, save :: parent_coord => null()
+!$omp threadprivate(cl, surf, lat, coord, final_coord, parent_coord)
! inialize distance to infinity (huge)
dist = INFINITY
- lattice_crossed = NONE
+ d_lat = INFINITY
+ d_surf = INFINITY
+ lattice_translation(:) = [0, 0, 0]
nullify(final_coord)
! Get pointer to top-level coordinates
@@ -1220,12 +1122,10 @@ contains
end select
! Check is calculated distance is new minimum
- if (d < dist) then
- if (abs(d - dist)/dist >= FP_PRECISION) then
- dist = d
- surface_crossed = -cl % surfaces(i)
- lattice_crossed = NONE
- final_coord => coord
+ if (d < d_surf) then
+ if (abs(d - d_surf)/d_surf >= FP_PRECISION) then
+ d_surf = d
+ level_surf_cross = -cl % surfaces(i)
end if
end if
@@ -1234,17 +1134,20 @@ contains
! =======================================================================
! FIND MINIMUM DISTANCE TO LATTICE SURFACES
- if (coord % lattice /= NONE) then
- lat => lattices(coord % lattice)
- if (lat % type == LATTICE_RECT) then
+ LAT_COORD: if (coord % lattice /= NONE) then
+ lat => lattices(coord % lattice) % obj
+
+ LAT_TYPE: select type(lat)
+
+ type is (RectLattice)
! copy local coordinates
x = coord % xyz(1)
y = coord % xyz(2)
z = coord % xyz(3)
! determine oncoming edge
- x0 = sign(lat % width(1) * 0.5_8, u)
- y0 = sign(lat % width(2) * 0.5_8, v)
+ x0 = sign(lat % pitch(1) * 0.5_8, u)
+ y0 = sign(lat % pitch(2) * 0.5_8, v)
! left and right sides
if (abs(x - x0) < FP_PRECISION) then
@@ -1255,23 +1158,11 @@ contains
d = (x0 - x)/u
end if
- ! If the lattice boundary is coincident with the parent cell boundary,
- ! we need to make sure that the lattice is not selected. This is
- ! complicated by the fact that floating point may determine that one
- ! is closer than the other (can't check direct equality). Thus, the
- ! logic here checks whether the relative difference is within floating
- ! point precision.
-
- if (d < dist) then
- if (abs(d - dist)/dist >= FP_REL_PRECISION) then
- dist = d
- if (u > 0) then
- lattice_crossed = LATTICE_RIGHT
- else
- lattice_crossed = LATTICE_LEFT
- end if
- final_coord => coord
- end if
+ d_lat = d
+ if (u > 0) then
+ level_lat_trans(:) = [1, 0, 0]
+ else
+ level_lat_trans(:) = [-1, 0, 0]
end if
! front and back sides
@@ -1283,20 +1174,17 @@ contains
d = (y0 - y)/v
end if
- if (d < dist) then
- if (abs(d - dist)/dist >= FP_REL_PRECISION) then
- dist = d
- if (v > 0) then
- lattice_crossed = LATTICE_FRONT
- else
- lattice_crossed = LATTICE_BACK
- end if
- final_coord => coord
+ if (d < d_lat) then
+ d_lat = d
+ if (v > 0) then
+ level_lat_trans(:) = [0, 1, 0]
+ else
+ level_lat_trans(:) = [0, -1, 0]
end if
end if
- if (lat % n_dimension == 3) then
- z0 = sign(lat % width(3) * 0.5_8, w)
+ if (lat % is_3d) then
+ z0 = sign(lat % pitch(3) * 0.5_8, w)
! top and bottom sides
if (abs(z - z0) < FP_PRECISION) then
@@ -1307,21 +1195,157 @@ contains
d = (z0 - z)/w
end if
- if (d < dist) then
- if (abs(d - dist)/dist >= FP_REL_PRECISION) then
- dist = d
- if (w > 0) then
- lattice_crossed = LATTICE_TOP
- else
- lattice_crossed = LATTICE_BOTTOM
- end if
- final_coord => coord
+ if (d < d_lat) then
+ d_lat = d
+ if (w > 0) then
+ level_lat_trans(:) = [0, 0, 1]
+ else
+ level_lat_trans(:) = [0, 0, -1]
end if
end if
end if
- elseif (lat % type == LATTICE_HEX) then
- ! TODO: Add hex lattice support
+ type is (HexLattice) LAT_TYPE
+ ! Copy local coordinates.
+ z = coord % xyz(3)
+ i_xyz(1) = coord % lattice_x
+ i_xyz(2) = coord % lattice_y
+ i_xyz(3) = coord % lattice_z
+ parent_coord => p % coord0
+ do while(.not. associated(parent_coord % next, coord))
+ parent_coord => parent_coord % next
+ end do
+
+ ! Compute velocities along the hexagonal axes.
+ beta_dir = u*sqrt(3.0_8)/2.0_8 + v/2.0_8
+ gama_dir = u*sqrt(3.0_8)/2.0_8 - v/2.0_8
+
+ ! Note that hexagonal lattice distance calculations are performed
+ ! using the particle's coordinates relative to the neighbor lattice
+ ! cells, not relative to the particle's current cell. This is done
+ ! because there is significant disagreement between neighboring cells
+ ! on where the lattice boundary is due to the worse finite precision
+ ! of hex lattices.
+
+ ! Upper right and lower left sides.
+ edge = -sign(lat % pitch(1)/2.0_8, beta_dir) ! Oncoming edge
+ if (beta_dir > 0.0) then
+ xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, 0, 0])
+ else
+ xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 0, 0])
+ end if
+ beta = xyz_t(1)*sqrt(3.0_8)/2.0_8 + xyz_t(2)/2.0_8
+ if (abs(beta - edge) < FP_PRECISION) then
+ d = INFINITY
+ else if (beta_dir == ZERO) then
+ d = INFINITY
+ else
+ d = (edge - beta)/beta_dir
+ end if
+
+ d_lat = d
+ if (beta_dir > 0) then
+ level_lat_trans(:) = [1, 0, 0]
+ else
+ level_lat_trans(:) = [-1, 0, 0]
+ end if
+
+ ! Lower right and upper left sides.
+ edge = -sign(lat % pitch(1)/2.0_8, gama_dir) ! Oncoming edge
+ if (gama_dir > 0.0) then
+ xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, -1, 0])
+ else
+ xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 1, 0])
+ end if
+ gama = xyz_t(1)*sqrt(3.0_8)/2.0_8 - xyz_t(2)/2.0_8
+ if (abs(gama - edge) < FP_PRECISION) then
+ d = INFINITY
+ else if (gama_dir == ZERO) then
+ d = INFINITY
+ else
+ d = (edge - gama)/gama_dir
+ end if
+
+ if (d < d_lat) then
+ d_lat = d
+ if (gama_dir > 0) then
+ level_lat_trans(:) = [1, -1, 0]
+ else
+ level_lat_trans(:) = [-1, 1, 0]
+ end if
+ end if
+
+ ! Upper and lower sides.
+ edge = -sign(lat % pitch(1)/2.0_8, v) ! Oncoming edge
+ if (v > 0.0) then
+ xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, 1, 0])
+ else
+ xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, -1, 0])
+ end if
+ if (abs(xyz_t(2) - edge) < FP_PRECISION) then
+ d = INFINITY
+ else if (v == ZERO) then
+ d = INFINITY
+ else
+ d = (edge - xyz_t(2))/v
+ end if
+
+ if (d < d_lat) then
+ d_lat = d
+ if (v > 0) then
+ level_lat_trans(:) = [0, 1, 0]
+ else
+ level_lat_trans(:) = [0, -1, 0]
+ end if
+ end if
+
+ ! Top and bottom sides.
+ if (lat % is_3d) then
+ z0 = sign(lat % pitch(2) * 0.5_8, w)
+
+ if (abs(z - z0) < FP_PRECISION) then
+ d = INFINITY
+ elseif (w == ZERO) then
+ d = INFINITY
+ else
+ d = (z0 - z)/w
+ end if
+
+ if (d < d_lat) then
+ d_lat = d
+ if (w > 0) then
+ level_lat_trans(:) = [0, 0, 1]
+ else
+ level_lat_trans(:) = [0, 0, -1]
+ end if
+ end if
+ end if
+ end select LAT_TYPE
+
+ if (d_lat < 0.0) then
+ call handle_lost_particle(p, "Particle " // trim(to_str(p % id)) &
+ &//" had a negative distance to a lattice boundary. d = " &
+ &//trim(to_str(d_lat)))
+ end if
+ end if LAT_COORD
+
+ ! If the boundary on this lattice level is coincident with a boundary on
+ ! a higher level then we need to make sure that the higher level boundary
+ ! is selected. This logic must include consideration of floating point
+ ! precision.
+ if (d_surf < d_lat) then
+ if ((dist - d_surf)/dist >= FP_REL_PRECISION) then
+ dist = d_surf
+ surface_crossed = level_surf_cross
+ lattice_translation(:) = [0, 0, 0]
+ final_coord => coord
+ end if
+ else
+ if ((dist - d_lat)/dist >= FP_REL_PRECISION) then
+ dist = d_lat
+ surface_crossed = None
+ lattice_translation(:) = level_lat_trans
+ final_coord => coord
end if
end if
diff --git a/src/geometry_header.F90 b/src/geometry_header.F90
index 57cb5f5906..6f6cfaecec 100644
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@@ -17,22 +17,98 @@ module geometry_header
end type Universe
!===============================================================================
-! LATTICE is an ordered array of elements (either rectangular, hexagonal, or
-! triangular)
+! LATTICE abstract type for ordered array of universes.
!===============================================================================
- type Lattice
- integer :: id ! Universe number for lattice
- integer :: type ! Type of lattice (rectangular, hex, etc)
- integer :: level ! Level of lattice
- integer :: n_dimension ! Number of dimensions
- integer, allocatable :: dimension(:) ! number of cells in each direction
- real(8), allocatable :: lower_left(:) ! lower-left corner of lattice
- real(8), allocatable :: width(:) ! width of each lattice cell
- integer, allocatable :: universes(:,:,:) ! specified universes
- integer :: outer ! universe to tile outside the lat
+ type, abstract :: Lattice
+ integer :: id ! Universe number for lattice
+ real(8), allocatable :: pitch(:) ! Pitch along each axis
+ integer, allocatable :: universes(:,:,:) ! Specified universes
+ integer :: outside ! Material to fill area outside
+ integer :: outer ! universe to tile outside the lat
+ logical :: is_3d ! Lattice has cells on z axis
+
+ contains
+
+ procedure(are_valid_indices_), deferred :: are_valid_indices
+ procedure(get_indices_), deferred :: get_indices
+ procedure(get_local_xyz_), deferred :: get_local_xyz
end type Lattice
+ abstract interface
+
+!===============================================================================
+! ARE_VALID_INDICES returns .true. if the given lattice indices fit within the
+! bounds of the lattice. Returns false otherwise.
+
+ function are_valid_indices_(this, i_xyz) result(is_valid)
+ import Lattice
+ class(Lattice), intent(in) :: this
+ integer, intent(in) :: i_xyz(3)
+ logical :: is_valid
+ end function are_valid_indices_
+
+!===============================================================================
+! GET_INDICES returns the indices in a lattice for the given global xyz.
+
+ function get_indices_(this, global_xyz) result(i_xyz)
+ import Lattice
+ class(Lattice), intent(in) :: this
+ real(8), intent(in) :: global_xyz(3)
+ integer :: i_xyz(3)
+ end function get_indices_
+
+!===============================================================================
+! GET_LOCAL_XYZ returns the translated local version of the given global xyz.
+
+ function get_local_xyz_(this, global_xyz, i_xyz) result(local_xyz)
+ import Lattice
+ class(Lattice), intent(in) :: this
+ real(8), intent(in) :: global_xyz(3)
+ integer, intent(in) :: i_xyz(3)
+ real(8) :: local_xyz(3)
+ end function get_local_xyz_
+ end interface
+
+!===============================================================================
+! RECTLATTICE extends LATTICE for rectilinear arrays.
+!===============================================================================
+
+ type, extends(Lattice) :: RectLattice
+ integer :: n_cells(3) ! Number of cells along each axis
+ real(8), allocatable :: lower_left(:) ! Global lower-left corner of lat
+
+ contains
+
+ procedure :: are_valid_indices => valid_inds_rect
+ procedure :: get_indices => get_inds_rect
+ procedure :: get_local_xyz => get_local_rect
+ end type RectLattice
+
+!===============================================================================
+! HEXLATTICE extends LATTICE for hexagonal (sometimes called triangular) arrays.
+!===============================================================================
+
+ type, extends(Lattice) :: HexLattice
+ integer :: n_rings ! Number of radial ring cell positoins
+ integer :: n_axial ! Number of axial cell positions
+ real(8), allocatable :: center(:) ! Global center of lattice
+
+ contains
+
+ procedure :: are_valid_indices => valid_inds_hex
+ procedure :: get_indices => get_inds_hex
+ procedure :: get_local_xyz => get_local_hex
+ end type HexLattice
+
+!===============================================================================
+! LATTICECONTAINER pointer array for storing lattices
+!===============================================================================
+
+ type LatticeContainer
+ class(Lattice), allocatable :: obj
+ end type LatticeContainer
+
!===============================================================================
! SURFACE type defines a first- or second-order surface that can be used to
! construct closed volumes (cells)
@@ -73,4 +149,166 @@ module geometry_header
! array index of universe 0
integer :: BASE_UNIVERSE
+contains
+
+!===============================================================================
+
+ function valid_inds_rect(this, i_xyz) result(is_valid)
+ class(RectLattice), intent(in) :: this
+ integer, intent(in) :: i_xyz(3)
+ logical :: is_valid
+
+ is_valid = all(i_xyz > 0 .and. i_xyz <= this % n_cells)
+ end function valid_inds_rect
+
+!===============================================================================
+
+ function valid_inds_hex(this, i_xyz) result(is_valid)
+ class(HexLattice), intent(in) :: this
+ integer, intent(in) :: i_xyz(3)
+ logical :: is_valid
+
+ is_valid = (all(i_xyz > 0) .and. &
+ &i_xyz(1) < 2*this % n_rings .and. &
+ &i_xyz(2) < 2*this % n_rings .and. &
+ &i_xyz(1) + i_xyz(2) > this % n_rings .and. &
+ &i_xyz(1) + i_xyz(2) < 3*this % n_rings .and. &
+ &i_xyz(3) <= this % n_axial)
+ end function valid_inds_hex
+
+!===============================================================================
+
+ function get_inds_rect(this, global_xyz) result(i_xyz)
+ class(RectLattice), intent(in) :: this
+ real(8), intent(in) :: global_xyz(3)
+ integer :: i_xyz(3)
+
+ real(8) :: xyz(3) ! global_xyz alias
+
+ xyz = global_xyz
+
+ i_xyz(1) = ceiling((xyz(1) - this % lower_left(1))/this % pitch(1))
+ i_xyz(2) = ceiling((xyz(2) - this % lower_left(2))/this % pitch(2))
+ if (this % is_3d) then
+ i_xyz(3) = ceiling((xyz(3) - this % lower_left(3))/this % pitch(3))
+ else
+ i_xyz(3) = 1
+ end if
+ end function get_inds_rect
+
+!===============================================================================
+
+ function get_inds_hex(this, global_xyz) result(i_xyz)
+ class(HexLattice), intent(in) :: this
+ real(8), intent(in) :: global_xyz(3)
+ integer :: i_xyz(3)
+
+ real(8) :: xyz(3) ! global_xyz alias
+ real(8) :: alpha ! Skewed coord axis
+ real(8) :: xyz_t(3) ! Local xyz
+ real(8) :: dists(4) ! Squared distances from cell centers
+ integer :: i, j, k ! Iterators
+ integer :: loc(1) ! Minimum distance index
+
+ xyz = global_xyz
+
+ ! Index z direction.
+ if (this % is_3d) then
+ i_xyz(3) = ceiling((xyz(3) - this % center(3))/this % pitch(2) + 0.5_8)&
+ &+ this % n_axial/2
+ else
+ i_xyz(3) = 1
+ end if
+
+ ! Convert coordinates into skewed bases. The (x, alpha) basis is used to
+ ! find the index of the global coordinates to within 4 cells.
+ alpha = xyz(2) - xyz(1) / sqrt(3.0_8)
+ i_xyz(1) = floor(xyz(1) / (sqrt(3.0_8) / 2.0_8 * this % pitch(1)))
+ i_xyz(2) = floor(alpha / this % pitch(1))
+
+ ! Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
+ ! the array is offset so that the indices never go below 1).
+ i_xyz(1) = i_xyz(1) + this % n_rings
+ i_xyz(2) = i_xyz(2) + this % n_rings
+
+ ! Calculate the (squared) distance between the particle and the centers of
+ ! the four possible cells. Regular hexagonal tiles form a centroidal
+ ! Voronoi tessellation so the global xyz should be in the hexagonal cell
+ ! that it is closest to the center of. This method is used over a
+ ! method that uses the remainders of the floor divisions above becasue it
+ ! provides better finite precision performance. Squared distances are
+ ! used becasue they are more computationally efficient than normal
+ ! distances.
+ k = 1
+ do i=0,1
+ do j=0,1
+ xyz_t = this % get_local_xyz(xyz, i_xyz + (/j, i, 0/))
+ dists(k) = xyz_t(1)**2 + xyz_t(2)**2
+ k = k + 1
+ end do
+ end do
+
+ ! Select the minimum squared distance which corresponds to the cell the
+ ! coordinates are in.
+ loc = minloc(dists)
+ if (loc(1) == 2) then
+ i_xyz = i_xyz + (/1, 0, 0/)
+ else if (loc(1) == 3) then
+ i_xyz = i_xyz + (/0, 1, 0/)
+ else if (loc(1) == 4) then
+ i_xyz = i_xyz + (/1, 1, 0/)
+ end if
+ end function get_inds_hex
+
+!===============================================================================
+
+ function get_local_rect(this, global_xyz, i_xyz) result(local_xyz)
+ class(RectLattice), intent(in) :: this
+ real(8), intent(in) :: global_xyz(3)
+ integer, intent(in) :: i_xyz(3)
+ real(8) :: local_xyz(3)
+
+ real(8) :: xyz(3) ! global_xyz alias
+
+ xyz = global_xyz
+
+ local_xyz(1) = xyz(1) - (this % lower_left(1) + &
+ &(i_xyz(1) - 0.5_8)*this % pitch(1))
+ local_xyz(2) = xyz(2) - (this % lower_left(2) + &
+ &(i_xyz(2) - 0.5_8)*this % pitch(2))
+ if (this % is_3d) then
+ local_xyz(3) = xyz(3) - (this % lower_left(3) + &
+ &(i_xyz(3) - 0.5_8)*this % pitch(3))
+ else
+ local_xyz(3) = xyz(3)
+ end if
+ end function get_local_rect
+
+!===============================================================================
+
+ function get_local_hex(this, global_xyz, i_xyz) result(local_xyz)
+ class(HexLattice), intent(in) :: this
+ real(8), intent(in) :: global_xyz(3)
+ integer, intent(in) :: i_xyz(3)
+ real(8) :: local_xyz(3)
+
+ real(8) :: xyz(3) ! global_xyz alias
+
+ xyz = global_xyz
+
+ ! x_l = x_g - (center + pitch_x*cos(30)*index_x)
+ local_xyz(1) = xyz(1) - (this % center(1) + &
+ &sqrt(3.0_8) / 2.0_8 * (i_xyz(1) - this % n_rings) * this % pitch(1))
+ ! y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
+ local_xyz(2) = xyz(2) - (this % center(2) + &
+ &(i_xyz(2) - this % n_rings) * this % pitch(1) + &
+ &(i_xyz(1) - this % n_rings) * this % pitch(1) / 2.0_8)
+ if (this % is_3d) then
+ local_xyz(3) = xyz(3) - this % center(3) &
+ &+ (this % n_axial/2 - i_xyz(3) + 1) * this % pitch(2)
+ else
+ local_xyz(3) = xyz(3)
+ end if
+ end function get_local_hex
+
end module geometry_header
diff --git a/src/global.F90 b/src/global.F90
index ec11216608..004c43746e 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -6,7 +6,7 @@ module global
use cmfd_header
use constants
use dict_header, only: DictCharInt, DictIntInt
- use geometry_header, only: Cell, Universe, Lattice, Surface
+ use geometry_header, only: Cell, Universe, Lattice, LatticeContainer, Surface
use material_header, only: Material
use mesh_header, only: StructuredMesh
use plot_header, only: ObjectPlot
@@ -26,12 +26,12 @@ module global
! GEOMETRY-RELATED VARIABLES
! Main arrays
- type(Cell), allocatable, target :: cells(:)
- type(Universe), allocatable, target :: universes(:)
- type(Lattice), allocatable, target :: lattices(:)
- type(Surface), allocatable, target :: surfaces(:)
- type(Material), allocatable, target :: materials(:)
- type(ObjectPlot),allocatable, target :: plots(:)
+ type(Cell), allocatable, target :: cells(:)
+ type(Universe), allocatable, target :: universes(:)
+ type(LatticeContainer), allocatable, target :: lattices(:)
+ type(Surface), allocatable, target :: surfaces(:)
+ type(Material), allocatable, target :: materials(:)
+ type(ObjectPlot), allocatable, target :: plots(:)
! Size of main arrays
integer :: n_cells ! # of cells
diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90
index 9d6a603f91..f318a47eaf 100644
--- a/src/hdf5_summary.F90
+++ b/src/hdf5_summary.F90
@@ -5,7 +5,8 @@ module hdf5_summary
use ace_header, only: Reaction, UrrData, Nuclide
use constants
use endf, only: reaction_name
- use geometry_header, only: Cell, Surface, Universe, Lattice
+ use geometry_header, only: Cell, Surface, Universe, Lattice, RectLattice, &
+ &HexLattice
use global
use material_header, only: Material
use mesh_header, only: StructuredMesh
@@ -99,13 +100,11 @@ contains
subroutine hdf5_write_geometry()
integer :: i, j, k, m
- integer :: n_x, n_y, n_z
- integer :: length(3)
integer, allocatable :: lattice_universes(:,:,:)
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
- type(Lattice), pointer :: lat => null()
+ class(Lattice), pointer :: lat => null()
! Use H5LT interface to write number of geometry objects
call su % write_data(n_cells, "n_cells", group="geometry")
@@ -174,8 +173,8 @@ contains
case (CELL_LATTICE)
call su % write_data("lattice", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
- call su % write_data(lattices(c % fill) % id, "lattice", &
- group="geometry/cells/cell " // trim(to_str(c % id)))
+ call su % write_data(lattices(c % fill) % obj % id, "lattice", &
+ group="geometry/cells/cell " // trim(to_str(c % id)))
end select
! Write list of bounding surfaces
@@ -306,7 +305,7 @@ contains
! Write information on each lattice
LATTICE_LOOP: do i = 1, n_lattices
- lat => lattices(i)
+ lat => lattices(i) % obj
! Write internal OpenMC index for this lattice
call su % write_data(i, "index", &
@@ -315,50 +314,125 @@ contains
! Write lattice type
select case(lat % type)
case (LATTICE_RECT)
+ select type (lat)
+ type is (RectLattice)
+ ! Write lattice type.
call su % write_data("rectangular", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
- case (LATTICE_HEX)
- call su % write_data("hexagonal", "type", &
+
+ ! Write number of lattice cells.
+ call su % write_data(lat % n_cells, "n_cells", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
- end select
- ! Write lattice dimensions, lower left corner, and width of element
- call su % write_data(lat % dimension, "dimension", &
- length=lat % n_dimension, &
- group="geometry/lattices/lattice " // trim(to_str(lat % id)))
- call su % write_data(lat % lower_left, "lower_left", &
- length=lat % n_dimension, &
- group="geometry/lattices/lattice " // trim(to_str(lat % id)))
- call su % write_data(lat % width, "width", &
- length=lat % n_dimension, &
- group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ ! Write lattice dimensions, lower left corner, and width of element
+ call su % write_data(lat % dimension, "dimension", &
+ length=lat % n_dimension, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ call su % write_data(lat % lower_left, "lower_left", &
+ length=lat % n_dimension, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ call su % write_data(lat % width, "width", &
+ length=lat % n_dimension, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
- call su % write_data(lat % outer, "outer", &
- group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ call su % write_data(lat % outer, "outer", &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
- ! Determine dimensions of lattice
- n_x = lat % dimension(1)
- n_y = lat % dimension(2)
- if (lat % n_dimension == 3) then
- n_z = lat % dimension(3)
- else
- n_z = 1
- end if
+ ! Determine dimensions of lattice
+ n_x = lat % dimension(1)
+ n_y = lat % dimension(2)
+ if (lat % n_dimension == 3) then
+ n_z = lat % dimension(3)
+ else
+ n_z = 1
+ end if
- ! Write lattice universes
- allocate(lattice_universes(n_x, n_y, n_z))
- do j = 1, n_x
- do k = 1, n_y
- do m = 1, n_z
- lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id
+ ! Write lattice lower-left.
+ if (lat % is_3d) then
+ call su % write_data(lat % lower_left, "lower_left", length=3, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ else
+ call su % write_data(lat % lower_left, "lower_left", length=2, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ end if
+
+ ! Write lattice pitch.
+ if (lat % is_3d) then
+ call su % write_data(lat % pitch, "pitch", length=3, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ else
+ call su % write_data(lat % pitch, "pitch", length=2, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ end if
+
+ ! Write lattice universes.
+ allocate(lattice_universes(lat % n_cells(1), lat % n_cells(2), &
+ &lat % n_cells(3)))
+ do j = 1, lat % n_cells(1)
+ do k = 1, lat % n_cells(2)
+ do m = 1, lat % n_cells(3)
+ lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id
+ end do
end do
end do
- end do
- length = [n_x, n_y, n_z]
- call su % write_data(lattice_universes, "universes", length=length, &
- group="geometry/lattices/lattice " // trim(to_str(lat % id)))
- deallocate(lattice_universes)
+ call su % write_data(lattice_universes, "universes", &
+ length=lat % n_cells, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ deallocate(lattice_universes)
+ type is (HexLattice)
+ ! Write lattice type.
+ call su % write_data("hexagonal", "type", &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+
+ ! Write number of lattice cells.
+ call su % write_data(lat % n_rings, "n_rings", &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ call su % write_data(lat % n_rings, "n_axial", &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+
+ ! Write lattice center.
+ if (lat % is_3d) then
+ call su % write_data(lat % center, "center", length=3, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ else
+ call su % write_data(lat % center, "center", length=2, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ end if
+
+ ! Write lattice pitch.
+ if (lat % is_3d) then
+ call su % write_data(lat % pitch, "pitch", length=2, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ else
+ call su % write_data(lat % pitch, "pitch", length=1, &
+ group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ end if
+
+ ! Write lattice universes.
+ allocate(lattice_universes(2*lat % n_rings - 1, 2*lat % n_rings - 1, &
+ &lat % n_axial))
+ do j = 1, lat % n_axial
+ do k = 1, 2*lat % n_rings - 1
+ do m = 1, 2*lat % n_rings - 1
+ if (j + k < lat % n_rings + 1) then
+ ! This array position is never used; put a -1 to indicate this
+ lattice_universes(j,k,m) = -1
+ cycle
+ else if (j + k > 3*lat % n_rings - 1) then
+ ! This array position is never used; put a -1 to indicate this
+ lattice_universes(j,k,m) = -1
+ cycle
+ end if
+ lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id
+ end do
+ end do
+ end do
+ call su % write_data(lattice_universes, "universes", &
+ &length=(/lat % n_axial, 2*lat % n_rings-1, 2*lat % n_rings-1/), &
+ &group="geometry/lattices/lattice " // trim(to_str(lat % id)))
+ deallocate(lattice_universes)
+ end select
end do LATTICE_LOOP
end subroutine hdf5_write_geometry
diff --git a/src/initialize.F90 b/src/initialize.F90
index ec64a1f19b..e39331675c 100644
--- a/src/initialize.F90
+++ b/src/initialize.F90
@@ -7,7 +7,8 @@ module initialize
use energy_grid, only: logarithmic_grid, grid_method
use error, only: fatal_error, warning
use geometry, only: neighbor_lists
- use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE
+ use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
+ &BASE_UNIVERSE
use global
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict, read_plots_xml
@@ -545,16 +546,15 @@ contains
subroutine adjust_indices()
- integer :: i ! index for various purposes
- integer :: j ! index for various purposes
- integer :: k ! loop index for lattices
- integer :: m ! loop index for lattices
- integer :: lid ! lattice IDs
- integer :: n_x, n_y, n_z ! size of lattice
- integer :: i_array ! index in surfaces/materials array
- integer :: id ! user-specified id
+ integer :: i ! index for various purposes
+ integer :: j ! index for various purposes
+ integer :: k ! loop index for lattices
+ integer :: m ! loop index for lattices
+ integer :: lid ! lattice IDs
+ integer :: i_array ! index in surfaces/materials array
+ integer :: id ! user-specified id
type(Cell), pointer :: c => null()
- type(Lattice), pointer :: lat => null()
+ class(Lattice), pointer :: lat => null()
type(TallyObject), pointer :: t => null()
do i = 1, n_cells
@@ -621,28 +621,47 @@ contains
! ADJUST UNIVERSE INDICES FOR EACH LATTICE
do i = 1, n_lattices
- lat => lattices(i)
- n_x = lat % dimension(1)
- n_y = lat % dimension(2)
- if (lat % n_dimension == 3) then
- n_z = lat % dimension(3)
- else
- n_z = 1
- end if
+ lat => lattices(i) % obj
+ select type (lat)
- do m = 1, n_z
- do k = 1, n_y
- do j = 1, n_x
- id = lat % universes(j,k,m)
- if (universe_dict % has_key(id)) then
- lat % universes(j,k,m) = universe_dict % get_key(id)
- else
- call fatal_error("Invalid universe number " // trim(to_str(id)) &
- &// " specified on lattice " // trim(to_str(lat % id)))
- end if
+ type is (RectLattice)
+ do m = 1, lat % n_cells(3)
+ do k = 1, lat % n_cells(2)
+ do j = 1, lat % n_cells(1)
+ id = lat % universes(j,k,m)
+ if (universe_dict % has_key(id)) then
+ lat % universes(j,k,m) = universe_dict % get_key(id)
+ else
+ call fatal_error("Invalid universe number " &
+ &// trim(to_str(id)) // " specified on lattice " &
+ &// trim(to_str(lat % id)))
+ end if
+ end do
end do
end do
- end do
+
+ type is (HexLattice)
+ do m = 1, lat % n_axial
+ do k = 1, 2*lat % n_rings - 1
+ do j = 1, 2*lat % n_rings - 1
+ if (j + k < lat % n_rings + 1) then
+ cycle
+ else if (j + k > 3*lat % n_rings - 1) then
+ cycle
+ end if
+ id = lat % universes(j, k, m)
+ if (universe_dict % has_key(id)) then
+ lat % universes(j, k, m) = universe_dict % get_key(id)
+ else
+ call fatal_error("Invalid universe number " &
+ &// trim(to_str(id)) // " specified on lattice " &
+ &// trim(to_str(lat % id)))
+ end if
+ end do
+ end do
+ end do
+
+ end select
if (lat % outer /= NO_OUTER_UNIVERSE) then
if (universe_dict % has_key(lat % outer)) then
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 45c3298781..17a5226063 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -5,7 +5,7 @@ module input_xml
use dict_header, only: DictIntInt, ElemKeyValueCI
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
- use geometry_header, only: Cell, Surface, Lattice
+ use geometry_header, only: Cell, Surface, Lattice, RectLattice, HexLattice
use global
use list_header, only: ListChar, ListReal
use mesh_header, only: StructuredMesh
@@ -898,29 +898,30 @@ contains
subroutine read_geometry_xml()
- integer :: i, j, k, m
+ integer :: i, j, k, m, i_x, i_a, input_index
integer :: n
- integer :: n_x, n_y, n_z
+ integer :: n_x, n_y, n_z, n_rings, n_rlats, n_hlats
integer :: universe_num
integer :: n_cells_in_univ
integer :: coeffs_reqd
- integer :: temp_int_array3(3)
+ integer :: temp_double_array3(3)
integer, allocatable :: temp_int_array(:)
real(8) :: phi, theta, psi
logical :: file_exists
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
- type(Cell), pointer :: c => null()
- type(Surface), pointer :: s => null()
- type(Lattice), pointer :: lat => null()
+ type(Cell), pointer :: c => null()
+ type(Surface), pointer :: s => null()
+ class(Lattice), pointer :: lat => null()
type(Node), pointer :: doc => null()
type(Node), pointer :: node_cell => null()
type(Node), pointer :: node_surf => null()
type(Node), pointer :: node_lat => null()
type(NodeList), pointer :: node_cell_list => null()
type(NodeList), pointer :: node_surf_list => null()
- type(NodeList), pointer :: node_lat_list => null()
+ type(NodeList), pointer :: node_rlat_list => null()
+ type(NodeList), pointer :: node_hlat_list => null()
! Display output message
call write_message("Reading geometry XML file...", 5)
@@ -1052,11 +1053,11 @@ contains
end if
! Copy rotation angles in x,y,z directions
- call get_node_array(node_cell, "rotation", temp_int_array3)
- c % rotation = temp_int_array3
- phi = -temp_int_array3(1) * PI/180.0_8
- theta = -temp_int_array3(2) * PI/180.0_8
- psi = -temp_int_array3(3) * PI/180.0_8
+ call get_node_array(node_cell, "rotation", temp_double_array3)
+ c % rotation = temp_double_array3
+ phi = -temp_double_array3(1) * PI/180.0_8
+ theta = -temp_double_array3(2) * PI/180.0_8
+ psi = -temp_double_array3(3) * PI/180.0_8
! Calculate rotation matrix based on angles given
allocate(c % rotation_matrix(3,3))
@@ -1241,17 +1242,23 @@ contains
! READ LATTICES FROM GEOMETRY.XML
! Get pointer to list of XML
- call get_node_list(doc, "lattice", node_lat_list)
+ call get_node_list(doc, "lattice", node_rlat_list)
+ call get_node_list(doc, "hex_lattice", node_hlat_list)
! Allocate lattices array
- n_lattices = get_list_size(node_lat_list)
+ n_rlats = get_list_size(node_rlat_list)
+ n_hlats = get_list_size(node_hlat_list)
+ n_lattices = n_rlats + n_hlats
allocate(lattices(n_lattices))
- do i = 1, n_lattices
- lat => lattices(i)
+ RECT_LATTICES: do i = 1, n_rlats
+ allocate(RectLattice::lattices(i) % obj)
+ lat => lattices(i) % obj
+ select type(lat)
+ type is (RectLattice)
! Get pointer to i-th lattice
- call get_list_item(node_lat_list, i, node_lat)
+ call get_list_item(node_rlat_list, i, node_lat)
! ID of lattice
if (check_for_node(node_lat, "id")) then
@@ -1266,32 +1273,21 @@ contains
&// to_str(lat % id))
end if
- ! Read lattice type
- word = ''
- if (check_for_node(node_lat, "type")) &
- call get_node_value(node_lat, "type", word)
- select case (to_lower(word))
- case ('rect', 'rectangle', 'rectangular')
- lat % type = LATTICE_RECT
- case ('hex', 'hexagon', 'hexagonal')
- lat % type = LATTICE_HEX
- case default
- call fatal_error("Invalid lattice type: " // trim(word))
- end select
-
! Read number of lattice cells in each dimension
n = get_arraysize_integer(node_lat, "dimension")
- if (n /= 2 .and. n /= 3) then
- call fatal_error("Lattice must be two or three dimensions.")
+ if (n == 2) then
+ call get_node_array(node_lat, "dimension", lat % n_cells(1:2))
+ lat % n_cells(3) = 1
+ lat % is_3d = .false.
+ else if (n == 3) then
+ call get_node_array(node_lat, "dimension", lat % n_cells)
+ lat % is_3d = .true.
+ else
+ call fatal_error("Rectangular lattice must be two or three dimensions.")
end if
- lat % n_dimension = n
- allocate(lat % dimension(n))
- call get_node_array(node_lat, "dimension", lat % dimension)
-
! Read lattice lower-left location
- if (size(lat % dimension) /= &
- get_arraysize_double(node_lat, "lower_left")) then
+ if (get_arraysize_double(node_lat, "lower_left") /= n) then
call fatal_error("Number of entries on must be the same &
&as the number of entries on .")
end if
@@ -1299,24 +1295,42 @@ contains
allocate(lat % lower_left(n))
call get_node_array(node_lat, "lower_left", lat % lower_left)
- ! Read lattice widths
- if (size(lat % dimension) /= &
- get_arraysize_double(node_lat, "width")) then
- call fatal_error("Number of entries on must be the same as &
- &the number of entries on .")
+ ! Read lattice pitches.
+ ! TODO: Remove this deprecation warning in a future release.
+ if (check_for_node(node_lat, "width")) then
+ call warning("The use of 'width' is deprecated and will be disallowed &
+ &in a future release. Use 'pitch' instead. The utility openmc/&
+ &src/utils/update_inputs.py can be used to automatically update &
+ &geometry.xml files.")
+ if (get_arraysize_double(node_lat, "width") /= n) then
+ call fatal_error("Number of entries on must be the same as &
+ &the number of entries on .")
+ end if
+
+ else if (get_arraysize_double(node_lat, "pitch") /= n) then
+ call fatal_error("Number of entries on must be the same as &
+ &the number of entries on .")
end if
- allocate(lat % width(n))
- call get_node_array(node_lat, "width", lat % width)
+ allocate(lat % pitch(n))
+ ! TODO: Remove the 'width' code in a future release.
+ if (check_for_node(node_lat, "width")) then
+ call get_node_array(node_lat, "width", lat % pitch)
+ else
+ call get_node_array(node_lat, "pitch", lat % pitch)
+ end if
+
+ ! TODO: Remove deprecation warning in a future release.
+ if (check_for_node(node_lat, "type")) then
+ call warning("The use of 'type' is no longer needed. The utility &
+ &openmc/src/utils/update_inputs.py can be used to automatically &
+ &update geometry.xml files.")
+ end if
! Copy number of dimensions
- n_x = lat % dimension(1)
- n_y = lat % dimension(2)
- if (lat % n_dimension == 3) then
- n_z = lat % dimension(3)
- else
- n_z = 1
- end if
+ n_x = lat % n_cells(1)
+ n_y = lat % n_cells(2)
+ n_z = lat % n_cells(3)
allocate(lat % universes(n_x, n_y, n_z))
! Check that number of universes matches size
@@ -1334,7 +1348,7 @@ contains
do k = 0, n_y - 1
do j = 1, n_x
lat % universes(j, n_y - k, m) = &
- temp_int_array(j + n_x*k + n_x*n_y*(m-1))
+ &temp_int_array(j + n_x*k + n_x*n_y*(m-1))
end do
end do
end do
@@ -1357,7 +1371,182 @@ contains
! Add lattice to dictionary
call lattice_dict % add_key(lat % id, i)
- end do
+ end select
+ end do RECT_LATTICES
+
+ HEX_LATTICES: do i = 1, n_hlats
+ allocate(HexLattice::lattices(n_rlats + i) % obj)
+ lat => lattices(n_rlats + i) % obj
+ select type (lat)
+ type is (HexLattice)
+
+ ! Get pointer to i-th lattice
+ call get_list_item(node_hlat_list, i, node_lat)
+
+ ! ID of lattice
+ if (check_for_node(node_lat, "id")) then
+ call get_node_value(node_lat, "id", lat % id)
+ else
+ call fatal_error("Must specify id of lattice in geometry XML file.")
+ end if
+
+ ! Check to make sure 'id' hasn't been used
+ if (lattice_dict % has_key(lat % id)) then
+ call fatal_error("Two or more lattices use the same unique ID: " &
+ &// to_str(lat % id))
+ end if
+
+ ! Read number of lattice cells in each dimension
+ call get_node_value(node_lat, "n_rings", lat % n_rings)
+ if (check_for_node(node_lat, "n_axial")) then
+ call get_node_value(node_lat, "n_axial", lat % n_axial)
+ lat % is_3d = .true.
+ else
+ lat % n_axial = 1
+ lat % is_3d = .false.
+ end if
+
+ ! Read lattice lower-left location
+ n = get_arraysize_double(node_lat, "center")
+ if (lat % is_3d .and. n /= 3) then
+ call fatal_error("A hexagonal lattice with must have &
+ &
specified by 3 numbers.")
+ else if ((.not. lat % is_3d) .and. n /= 2) then
+ call fatal_error("A hexagonal lattice without must have &
+ &