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

Merge branch 'develop' into mpif08

Browse source 
Conflicts:
	CMakeLists.txt
This commit is contained in:
Paul Romano 2015-07-16 07:03:34 +02:00
commit f8a6203d7f
23 changed files with 556 additions and 447 deletions
Download patch file Download diff file Expand all files Collapse all files

21
.travis.yml
View file

@ -1,11 +1,21 @@
sudo: false
language: python
python:
- "2.7"
- "3.4"
addons:
apt:
packages:
- gfortran
- g++
cache:
directories:
- $HOME/mpich_install
- $HOME/hdf5_install
- $HOME/phdf5_install
before_install:
# ============== Handle Python third-party packages ==============
- sudo apt-get update
- if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then
wget https://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh -O miniconda.sh;
else
@ -20,13 +30,12 @@ before_install:
- conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy h5py
- source activate test-environment
# ============== Install GCC, MPICH, HDF5, PHDF5 ==============
- sudo apt-get install -qq -y gfortran g++
# Install GCC, MPICH, HDF5, PHDF5
- ./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 MPI_DIR=$HOME/mpich_install
- export PHDF5_DIR=$HOME/phdf5_install
- export HDF5_DIR=$HOME/hdf5_install
install: true

6
CMakeLists.txt
View file

@ -24,13 +24,17 @@ option(profile "Compile with profiling flags" OFF)
option(debug "Compile with debug flags" OFF)
option(optimize "Turn on all compiler optimization flags" OFF)
option(verbose "Create verbose Makefiles" OFF)
option(coverage "Compile with flags" OFF)
option(coverage "Compile with coverage analysis flags" OFF)
option(mpif08 "Use Fortran 2008 MPI interface" OFF)
if (verbose)
set(CMAKE_VERBOSE_MAKEFILE on)
endif()
# Maximum number of nested coordinates levels
set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels")
add_definitions(-DMAX_COORD=${maxcoord})
#===============================================================================
# MPI for distributed-memory parallelism / HDF5 for binary output
#===============================================================================

12
docs/source/devguide/styleguide.rst
View file

@ -125,8 +125,16 @@ program, subroutine, function, if, associate, etc. Emacs users should set the
variables f90-if-indent, f90-do-indent, f90-continuation-indent,
f90-type-indent, f90-associate-indent, and f90-program indent to 2.
Continuation lines should be indented by an extra 5 spaces. This is the default
value of f90-continuation-indent in Emacs.
Continuation lines should be indented by at least 5 spaces. They may be indented
more in order to make the content match the context. For example, either of
these are valid continuation indentations:
.. code-block:: fortran
local_xyz(1) = xyz(1) - (this % lower_left(1) + &
(i_xyz(1) - HALF)*this % pitch(1))
call which_data(scatt_type, get_scatt, get_nuscatt, get_chi_t, get_chi_p, &
get_chi_d, scatt_order)
Whitespace in Expressions
-------------------------

9
docs/source/usersguide/install.rst
View file

@ -160,6 +160,13 @@ openmp
Enables shared-memory parallelism using the OpenMP API. The Fortran compiler
being used must support OpenMP.
coverage
Compile and link code instrumented for coverage analysis. This is typically
used in conjunction with gcov_.
maxcoord
Maximum number of nested coordinate levels in geometry. Defaults to 10.
To set any of these options (e.g. turning on debug mode), the following form
should be used:
@ -167,6 +174,8 @@ should be used:
cmake -Ddebug=on /path/to/openmc
.. _gcov: https://gcc.gnu.org/onlinedocs/gcc/Gcov.html
Compiling with MPI
++++++++++++++++++

2
openmc/geometry.py
View file

@ -205,7 +205,7 @@ class GeometryFile(object):
sort_xml_elements(self._geometry_file)
clean_xml_indentation(self._geometry_file)
# Write the XML Tree to the materials.xml file
# Write the XML Tree to the geometry.xml file
tree = ET.ElementTree(self._geometry_file)
tree.write("geometry.xml", xml_declaration=True,
encoding='utf-8', method="xml")

8
openmc/universe.py
View file

@ -935,6 +935,7 @@ class RectLattice(Lattice):
if self._outer is not None:
outer = ET.SubElement(lattice_subelement, "outer")
outer.text = '{0}'.format(self._outer._id)
self._outer.create_xml_subelement(xml_element)
# Export Lattice cell dimensions
dimension = ET.SubElement(lattice_subelement, "dimension")
@ -1059,7 +1060,7 @@ class HexLattice(Lattice):
@Lattice.pitch.setter
def pitch(self, pitch):
check_type('lattice pitch', pitch, Iterable, Real)
check_length('lattice pitch', pitch, 2, 3)
check_length('lattice pitch', pitch, 1, 2)
for dim in pitch:
check_greater_than('lattice pitch', dim, 0)
self._pitch = pitch
@ -1197,6 +1198,7 @@ class HexLattice(Lattice):
if self._outer is not None:
outer = ET.SubElement(lattice_subelement, "outer")
outer.text = '{0}'.format(self._outer._id)
self._outer.create_xml_subelement(xml_element)
lattice_subelement.set("n_rings", str(self._num_rings))
@ -1236,8 +1238,8 @@ class HexLattice(Lattice):
universe.create_xml_subelement(xml_element)
# Initialize the remaining universes.
for r in range(self._num_rings-1):
for theta in range(2*(self._num_rings - r)):
for r in range(self._num_rings - 1):
for theta in range(6*(self._num_rings - 1 - r)):
universe = self._universes[r][theta]
universe.create_xml_subelement(xml_element)

20
src/eigenvalue.F90
View file

@ -171,6 +171,26 @@ contains
subroutine finalize_generation()
! Update global tallies with the omp private accumulation variables
!$omp parallel
!$omp critical
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + global_tally_collision
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + global_tally_leakage
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + global_tally_absorption
!$omp end critical
! reset private tallies
global_tally_tracklength = 0
global_tally_collision = 0
global_tally_leakage = 0
global_tally_absorption = 0
!$omp end parallel
#ifdef _OPENMP
! Join the fission bank from each thread into one global fission bank
call join_bank_from_threads()

373
src/geometry.F90
View file

@ -6,7 +6,7 @@ module geometry
&RectLattice, HexLattice
use global
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord, Particle
use particle_header, only: LocalCoord, Particle
use particle_restart_write, only: write_particle_restart
use string, only: to_str
use tally, only: score_surface_current
@ -76,20 +76,18 @@ contains
type(Particle), intent(inout) :: p
integer :: i ! cell loop index on a level
integer :: j ! coordinate level index
integer :: n_coord ! saved number of coordinate levels
integer :: n ! number of cells to search on a level
integer :: index_cell ! index in cells array
type(Cell), pointer :: c ! pointer to cell
type(Universe), pointer :: univ ! universe to search in
type(LocalCoord), pointer :: coord ! particle coordinate to search on
coord => p % coord0
! loop through each coordinate level
do while (associated(coord))
p % coord => coord
univ => universes(coord % universe)
n_coord = p % n_coord
do j = 1, n_coord
p % n_coord = j
univ => universes(p % coord(j) % universe)
n = univ % n_cells
! loop through each cell on this level
@ -99,10 +97,10 @@ contains
if (simple_cell_contains(c, p)) then
! the particle should only be contained in one cell per level
if (index_cell /= coord % cell) then
if (index_cell /= p % coord(j) % cell) then
call fatal_error("Overlapping cells detected: " &
&// trim(to_str(cells(index_cell) % id)) // ", " &
&// trim(to_str(cells(coord % cell) % id)) &
&// trim(to_str(cells(p % coord(j) % cell) % id)) &
&// " on universe " // trim(to_str(univ % id)))
end if
@ -111,9 +109,6 @@ contains
end if
end do
coord => coord % next
end do
end subroutine check_cell_overlap
@ -130,6 +125,7 @@ contains
logical, intent(inout) :: found
integer, optional :: search_cells(:)
integer :: i ! index over cells
integer :: j ! coordinate level index
integer :: i_xyz(3) ! indices in lattice
integer :: n ! number of cells to search
integer :: index_cell ! index in cells array
@ -138,8 +134,10 @@ contains
class(Lattice), pointer :: lat ! pointer to lattice
type(Universe), pointer :: univ ! universe to search in
! Remove coordinates for any lower levels
call deallocate_coord(p % coord % next)
do j = p % n_coord + 1, MAX_COORD
call p % coord(j) % reset()
end do
j = p % n_coord
! set size of list to search
if (present(search_cells)) then
@ -147,7 +145,7 @@ contains
n = size(search_cells)
else
use_search_cells = .false.
univ => universes(p % coord % universe)
univ => universes(p % coord(j) % universe)
n = univ % n_cells
end if
@ -157,7 +155,7 @@ contains
if (use_search_cells) then
index_cell = search_cells(i)
! check to make sure search cell is in same universe
if (cells(index_cell) % universe /= p % coord % universe) cycle
if (cells(index_cell) % universe /= p % coord(j) % universe) cycle
else
index_cell = univ % cells(i)
end if
@ -169,7 +167,7 @@ contains
if (.not. simple_cell_contains(c, p)) cycle
! Set cell on this level
p % coord % cell = index_cell
p % coord(j) % cell = index_cell
! Show cell information on trace
if (verbosity >= 10 .or. trace) then
@ -188,28 +186,29 @@ contains
! ======================================================================
! 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
! Store lower level coordinates
p % coord(j + 1) % xyz = p % coord(j) % xyz
p % coord(j + 1) % uvw = p % coord(j) % uvw
! Move particle to next level and set universe
p % coord => p % coord % next
p % coord % universe = c % fill
j = j + 1
p % n_coord = j
p % coord(j) % universe = c % fill
! Apply translation
if (allocated(c % translation)) then
p % coord % xyz = p % coord % xyz - c % translation
p % coord(j) % xyz = p % coord(j) % 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.
p % coord(j) % xyz = matmul(c % rotation_matrix, p % coord(j) % xyz)
p % coord(j) % uvw = matmul(c % rotation_matrix, p % coord(j) % uvw)
p % coord(j) % rotated = .true.
end if
call find_cell(p, found)
j = p % n_coord
if (.not. found) exit
elseif (c % type == CELL_LATTICE) then CELL_TYPE
@ -220,40 +219,41 @@ contains
lat => lattices(c % fill) % obj
! Determine lattice indices
i_xyz = lat % get_indices(p % coord % xyz + TINY_BIT * p % coord % uvw)
i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % 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
! Store lower level coordinates
p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz)
p % coord(j + 1) % uvw = p % coord(j) % 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)
p % coord(j + 1) % lattice = c % fill
p % coord(j + 1) % lattice_x = i_xyz(1)
p % coord(j + 1) % lattice_y = i_xyz(2)
p % coord(j + 1) % 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))
p % coord(j + 1) % 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 &
// trim(to_str(lat % id)) // " but the lattice has no &
&defined outer universe.")
else
p % coord % next % universe = lat % outer
p % coord(j + 1) % universe = lat % outer
end if
end if
! Move particle to next level and search for the lower cells.
p % coord => p % coord % next
j = j + 1
p % n_coord = j
call find_cell(p, found)
j = p % n_coord
if (.not. found) exit
end if CELL_TYPE
@ -314,15 +314,13 @@ contains
! TODO: Find a better solution to score surface currents than
! physically moving the particle forward slightly
p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
call score_surface_current(p)
end if
! Score to global leakage tally
if (tallies_on) then
!$omp atomic
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + p % wgt
global_tally_leakage = global_tally_leakage + p % wgt
end if
! Display message
@ -337,7 +335,7 @@ contains
! PARTICLE REFLECTS FROM SURFACE
! Do not handle reflective boundary conditions on lower universes
if (.not. associated(p % coord, p % coord0)) then
if (p % n_coord /= 1) then
call handle_lost_particle(p, "Cannot reflect particle " &
&// trim(to_str(p % id)) // " off surface in a lower universe.")
return
@ -348,15 +346,15 @@ contains
! case the surface crossing in coincident with a mesh boundary
if (active_current_tallies % size() > 0) then
p % coord0 % xyz = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
call score_surface_current(p)
p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
end if
! Copy particle's direction cosines
u = p % coord0 % uvw(1)
v = p % coord0 % uvw(2)
w = p % coord0 % uvw(3)
u = p % coord(1) % uvw(1)
v = p % coord(1) % uvw(2)
w = p % coord(1) % uvw(3)
select case (surf%type)
case (SURF_PX)
@ -383,8 +381,8 @@ contains
case (SURF_CYL_X)
! Find y-y0, z-z0 and dot product of direction and surface normal
y = p % coord0 % xyz(2) - surf % coeffs(1)
z = p % coord0 % xyz(3) - surf % coeffs(2)
y = p % coord(1) % xyz(2) - surf % coeffs(1)
z = p % coord(1) % xyz(3) - surf % coeffs(2)
R = surf % coeffs(3)
dot_prod = v*y + w*z
@ -394,8 +392,8 @@ contains
case (SURF_CYL_Y)
! Find x-x0, z-z0 and dot product of direction and surface normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
z = p % coord0 % xyz(3) - surf % coeffs(2)
x = p % coord(1) % xyz(1) - surf % coeffs(1)
z = p % coord(1) % xyz(3) - surf % coeffs(2)
R = surf % coeffs(3)
dot_prod = u*x + w*z
@ -405,8 +403,8 @@ contains
case (SURF_CYL_Z)
! Find x-x0, y-y0 and dot product of direction and surface normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
x = p % coord(1) % xyz(1) - surf % coeffs(1)
y = p % coord(1) % xyz(2) - surf % coeffs(2)
R = surf % coeffs(3)
dot_prod = u*x + v*y
@ -417,9 +415,9 @@ contains
case (SURF_SPHERE)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
z = p % coord0 % xyz(3) - surf % coeffs(3)
x = p % coord(1) % xyz(1) - surf % coeffs(1)
y = p % coord(1) % xyz(2) - surf % coeffs(2)
z = p % coord(1) % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
dot_prod = u*x + v*y + w*z
@ -431,9 +429,9 @@ contains
case (SURF_CONE_X)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
z = p % coord0 % xyz(3) - surf % coeffs(3)
x = p % coord(1) % xyz(1) - surf % coeffs(1)
y = p % coord(1) % xyz(2) - surf % coeffs(2)
z = p % coord(1) % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
dot_prod = (v*y + w*z - R*u*x)/((R + ONE)*R*x*x)
@ -445,9 +443,9 @@ contains
case (SURF_CONE_Y)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
z = p % coord0 % xyz(3) - surf % coeffs(3)
x = p % coord(1) % xyz(1) - surf % coeffs(1)
y = p % coord(1) % xyz(2) - surf % coeffs(2)
z = p % coord(1) % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
dot_prod = (u*x + w*z - R*v*y)/((R + ONE)*R*y*y)
@ -459,9 +457,9 @@ contains
case (SURF_CONE_Z)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
y = p % coord0 % xyz(2) - surf % coeffs(2)
z = p % coord0 % xyz(3) - surf % coeffs(3)
x = p % coord(1) % xyz(1) - surf % coeffs(1)
y = p % coord(1) % xyz(2) - surf % coeffs(2)
z = p % coord(1) % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
dot_prod = (u*x + v*y - R*w*z)/((R + ONE)*R*z*z)
@ -477,28 +475,26 @@ contains
! Set new particle direction
norm = sqrt(u*u + v*v + w*w)
p % coord0 % uvw = [u, v, w] / norm
p % coord(1) % uvw = [u, v, w] / norm
! Reassign particle's cell and surface
p % coord0 % cell = last_cell
p % coord(1) % cell = last_cell
p % surface = -p % surface
! If a reflective surface is coincident with a lattice or universe
! boundary, it is necessary to redetermine the particle's coordinates in
! the lower universes.
if (associated(p % coord0 % next)) then
call deallocate_coord(p % coord0 % next)
call find_cell(p, found)
if (.not. found) then
call handle_lost_particle(p, "Couldn't find particle after reflecting&
& from surface.")
return
end if
p % n_coord = 1
call find_cell(p, found)
if (.not. found) then
call handle_lost_particle(p, "Couldn't find particle after reflecting&
& from surface.")
return
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then
@ -532,8 +528,7 @@ contains
! Remove lower coordinate levels and assignment of surface
p % surface = NONE
p % coord => p % coord0
call deallocate_coord(p % coord % next)
p % n_coord = 1
call find_cell(p, found)
if (run_mode /= MODE_PLOTTING .and. (.not. found)) then
@ -542,9 +537,8 @@ contains
! the particle is really traveling tangent to a surface, if we move it
! forward a tiny bit it should fix the problem.
p % coord => p % coord0
call deallocate_coord(p % coord % next)
p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
p % n_coord = 1
p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
call find_cell(p, found)
! Couldn't find next cell anywhere! This probably means there is an actual
@ -568,52 +562,46 @@ contains
type(Particle), intent(inout) :: p
integer, intent(in) :: lattice_translation(3)
integer :: j
integer :: i_xyz(3) ! indices in lattice
logical :: found ! particle found in cell?
class(Lattice), pointer :: lat
type(LocalCoord), pointer :: parent_coord
lat => lattices(p % coord % lattice) % obj
j = p % n_coord
lat => lattices(p % coord(j) % lattice) % obj
if (verbosity >= 10 .or. trace) then
call write_message(" Crossing lattice " // trim(to_str(lat % id)) &
&// ". Current position (" // trim(to_str(p % coord % lattice_x)) &
&// "," // trim(to_str(p % coord % lattice_y)) // "," &
&// trim(to_str(p % coord % lattice_z)) // ")")
&// ". Current position (" // trim(to_str(p % coord(j) % lattice_x)) &
&// "," // trim(to_str(p % coord(j) % lattice_y)) // "," &
&// trim(to_str(p % coord(j) % lattice_z)) // ")")
end if
! 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
! 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
p % coord(j) % lattice_x = p % coord(j) % lattice_x + lattice_translation(1)
p % coord(j) % lattice_y = p % coord(j) % lattice_y + lattice_translation(2)
p % coord(j) % lattice_z = p % coord(j) % lattice_z + lattice_translation(3)
i_xyz(1) = p % coord(j) % lattice_x
i_xyz(2) = p % coord(j) % lattice_y
i_xyz(3) = p % coord(j) % lattice_z
! Set the new coordinate position.
p % coord % xyz = lat % get_local_xyz(parent_coord % xyz, i_xyz)
p % coord(j) % xyz = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz)
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
! The particle is outside the lattice. Search for it from base coord
p % n_coord = 1
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.")
// trim(to_str(p % id)) // " after crossing a lattice boundary.")
return
end if
else OUTSIDE_LAT
! Find cell in next lattice element
p % coord % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
p % coord(j) % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
call find_cell(p, found)
if (.not. found) then
@ -622,15 +610,14 @@ contains
! off all lower-level coordinates and search from universe zero
! Remove lower coordinates
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
p % n_coord = 1
! 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.")
// trim(to_str(p % id)) &
// " after crossing a lattice boundary.")
return
end if
end if
@ -644,14 +631,17 @@ 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_translation)
subroutine distance_to_boundary(p, dist, surface_crossed, lattice_translation, &
next_level)
type(Particle), intent(inout) :: p
real(8), intent(out) :: dist
integer, intent(out) :: surface_crossed
integer, intent(out) :: lattice_translation(3)
integer, intent(out) :: next_level
integer :: i ! index for surface in cell
integer :: j
integer :: index_surf ! index in surfaces array (with sign)
integer :: i_xyz(3) ! lattice indices
integer :: level_surf_cross ! surface crossed on current level
@ -675,30 +665,25 @@ contains
type(Cell), pointer :: cl
type(Surface), pointer :: surf
class(Lattice), pointer :: lat
type(LocalCoord), pointer :: coord
type(LocalCoord), pointer :: final_coord
type(LocalCoord), pointer :: parent_coord
! inialize distance to infinity (huge)
dist = INFINITY
d_lat = INFINITY
d_surf = INFINITY
lattice_translation(:) = [0, 0, 0]
nullify(final_coord)
! Get pointer to top-level coordinates
coord => p % coord0
next_level = 0
! Loop over each universe level
LEVEL_LOOP: do while(associated(coord))
LEVEL_LOOP: do j = 1, p % n_coord
! get pointer to cell on this level
cl => cells(coord % cell)
cl => cells(p % coord(j) % cell)
! copy directional cosines
u = coord % uvw(1)
v = coord % uvw(2)
w = coord % uvw(3)
u = p % coord(j) % uvw(1)
v = p % coord(j) % uvw(2)
w = p % coord(j) % uvw(3)
! =======================================================================
! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL
@ -706,9 +691,9 @@ contains
SURFACE_LOOP: do i = 1, cl % n_surfaces
! copy local coordinates of particle
x = coord % xyz(1)
y = coord % xyz(2)
z = coord % xyz(3)
x = p % coord(j) % xyz(1)
y = p % coord(j) % xyz(2)
z = p % coord(j) % xyz(3)
! check for coincident surface -- note that we can't skip the
! calculation in general because a particle could be on one side of a
@ -1129,16 +1114,16 @@ contains
! =======================================================================
! FIND MINIMUM DISTANCE TO LATTICE SURFACES
LAT_COORD: if (coord % lattice /= NONE) then
lat => lattices(coord % lattice) % obj
LAT_COORD: if (p % coord(j) % lattice /= NONE) then
lat => lattices(p % coord(j) % 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)
x = p % coord(j) % xyz(1)
y = p % coord(j) % xyz(2)
z = p % coord(j) % xyz(3)
! determine oncoming edge
x0 = sign(lat % pitch(1) * HALF, u)
@ -1202,14 +1187,10 @@ contains
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
z = p % coord(j) % xyz(3)
i_xyz(1) = p % coord(j) % lattice_x
i_xyz(2) = p % coord(j) % lattice_y
i_xyz(3) = p % coord(j) % lattice_z
! Compute velocities along the hexagonal axes.
beta_dir = u*sqrt(THREE)/TWO + v/TWO
@ -1225,9 +1206,9 @@ contains
! Upper right and lower left sides.
edge = -sign(lat % pitch(1)/TWO, beta_dir) ! Oncoming edge
if (beta_dir > ZERO) then
xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, 0, 0])
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, 0, 0])
else
xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 0, 0])
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 0, 0])
end if
beta = xyz_t(1)*sqrt(THREE)/TWO + xyz_t(2)/TWO
if (abs(beta - edge) < FP_PRECISION) then
@ -1248,9 +1229,9 @@ contains
! Lower right and upper left sides.
edge = -sign(lat % pitch(1)/TWO, gama_dir) ! Oncoming edge
if (gama_dir > ZERO) then
xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, -1, 0])
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, -1, 0])
else
xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 1, 0])
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 1, 0])
end if
gama = xyz_t(1)*sqrt(THREE)/TWO - xyz_t(2)/TWO
if (abs(gama - edge) < FP_PRECISION) then
@ -1273,9 +1254,9 @@ contains
! Upper and lower sides.
edge = -sign(lat % pitch(1)/TWO, v) ! Oncoming edge
if (v > ZERO) then
xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, 1, 0])
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, 1, 0])
else
xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, -1, 0])
xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, -1, 0])
end if
if (abs(xyz_t(2) - edge) < FP_PRECISION) then
d = INFINITY
@ -1333,24 +1314,19 @@ contains
dist = d_surf
surface_crossed = level_surf_cross
lattice_translation(:) = [0, 0, 0]
final_coord => coord
next_level = j
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
next_level = j
end if
end if
coord => coord % next
end do LEVEL_LOOP
! Move particle to appropriate coordinate level
if (associated(final_coord)) p % coord => final_coord
end subroutine distance_to_boundary
!===============================================================================
@ -1365,6 +1341,7 @@ contains
type(Surface), pointer :: surf ! surface
logical :: s ! sense of particle
integer :: j
real(8) :: x,y,z ! coordinates of particle
real(8) :: func ! surface function evaluated at point
real(8) :: A ! coefficient on x for plane
@ -1374,9 +1351,10 @@ contains
real(8) :: x0,y0,z0 ! coefficients for quadratic surfaces / box
real(8) :: r ! radius for quadratic surfaces
x = p % coord % xyz(1)
y = p % coord % xyz(2)
z = p % coord % xyz(3)
j = p % n_coord
x = p % coord(j) % xyz(1)
y = p % coord(j) % xyz(2)
z = p % coord(j) % xyz(3)
select case (surf % type)
case (SURF_PX)
@ -1468,7 +1446,7 @@ contains
if (abs(func) < FP_COINCIDENT) then
! Particle may be coincident with this surface. Artifically move the
! particle forward a tiny bit.
p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
p % coord(j) % xyz = p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw
s = sense(p, surf)
elseif (func > 0) then
s = .true.
@ -1907,5 +1885,82 @@ contains
end subroutine count_instance
!===============================================================================
! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in
! the geometry
!===============================================================================
recursive function maximum_levels(univ) result(levels)
type(Universe), intent(in) :: univ ! universe to search through
integer :: levels ! maximum number of levels for this universe
integer :: i ! index over cells
integer :: j, k, m ! indices in lattice
integer :: levels_below ! max levels below this universe
type(Cell), pointer :: c ! pointer to current cell
type(Universe), pointer :: next_univ ! next universe to loop through
class(Lattice), pointer :: lat ! pointer to current lattice
levels_below = 0
do i = 1, univ % n_cells
c => cells(univ % cells(i))
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == CELL_FILL) then
next_univ => universes(c % fill)
levels_below = max(levels_below, maximum_levels(next_univ))
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
type is (RectLattice)
! Loop over lattice coordinates
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
next_univ => universes(lat % universes(j, k, m))
levels_below = max(levels_below, maximum_levels(next_univ))
end do
end do
end do
type is (HexLattice)
! Loop over lattice coordinates
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
! This array location is never used
if (j + k < lat % n_rings + 1) then
cycle
! This array location is never used
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
next_univ => universes(lat % universes(j, k, m))
levels_below = max(levels_below, maximum_levels(next_univ))
end if
end do
end do
end do
end select
end if
end do
levels = 1 + levels_below
end function maximum_levels
end module geometry

17
src/global.F90
View file

@ -115,11 +115,24 @@ module global
! Global tallies
! 1) collision estimate of k-eff
! 2) track-length estimate of k-eff
! 3) leakage fraction
! 2) absorption estimate of k-eff
! 3) track-length estimate of k-eff
! 4) leakage fraction
type(TallyResult), allocatable, target :: global_tallies(:)
! It is possible to protect accumulate operations on global tallies by using
! an atomic update. However, when multiple threads accumulate to the same
! global tally, it can cause a higher cache miss rate due to
! invalidation. Thus, we use threadprivate variables to accumulate global
! tallies and then reduce at the end of a generation.
real(8) :: global_tally_collision = ZERO
real(8) :: global_tally_absorption = ZERO
real(8) :: global_tally_tracklength = ZERO
real(8) :: global_tally_leakage = ZERO
!$omp threadprivate(global_tally_collision, global_tally_absorption, &
!$omp& global_tally_tracklength, global_tally_leakage)
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)

66
src/initialize.F90
View file

@ -7,7 +7,8 @@ module initialize
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method, unionized_grid
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
use global
@ -95,11 +96,20 @@ contains
! Initialize distribcell_filters
call prepare_distribcell()
! After reading input and basic geometry setup is complete, build lists of
! neighboring cells for efficient tracking
call neighbor_lists()
! Check to make sure there are not too many nested coordinate levels in the
! geometry since the coordinate list is statically allocated for performance
! reasons
if (maximum_levels(universes(BASE_UNIVERSE)) > MAX_COORD) then
call fatal_error("Too many nested coordinate levels in the geometry. &
&Try increasing the maximum number of coordinate levels by &
&providing the CMake -Dmaxcoord= option.")
end if
if (run_mode /= MODE_PLOTTING) then
! With the AWRs from the xs_listings, change all material specifications
! so that they contain atom percents summing to 1
@ -942,7 +952,7 @@ contains
count_all = .false.
! Loop over tallies
! Loop over tallies
do i = 1, n_tallies
! Get pointer to tally
@ -960,25 +970,25 @@ contains
if (size(tally % filters(j) % int_bins) > 1) then
call fatal_error("A distribcell filter was specified with &
&multiple bins. This feature is not supported.")
end if
end if
end if
end do
end do
if (count_all) then
univ => universes(BASE_UNIVERSE)
! sum the number of occurrences of all cells
call count_instance(univ)
! Loop over tallies
do i = 1, n_tallies
! Loop over tallies
do i = 1, n_tallies
! Get pointer to tally
tally => tallies(i)
tally => tallies(i)
! Initialize the filters
do j = 1, tally % n_filters
@ -999,7 +1009,7 @@ contains
! Calculate offsets for each target distribcell
do i = 1, n_maps
do j = 1, n_universes
do j = 1, n_universes
univ => universes(j)
call calc_offsets(univ_list(i), i, univ, counts, found)
end do
@ -1009,7 +1019,7 @@ contains
deallocate(counts)
deallocate(found)
deallocate(univ_list)
end subroutine prepare_distribcell
!===============================================================================
@ -1024,31 +1034,31 @@ contains
logical, intent(out), allocatable :: found(:,:) ! Target found
integer :: i, j, k, l, m ! Loop counters
type(SetInt) :: cell_list ! distribells to track
type(SetInt) :: cell_list ! distribells to track
type(Universe), pointer :: univ ! pointer to universe
class(Lattice), pointer :: lat ! pointer to lattice
type(TallyObject), pointer :: tally ! pointer to tally
type(TallyFilter), pointer :: filter ! pointer to filter
! Begin gathering list of cells in distribcell tallies
n_maps = 0
! Populate list of distribcells to track
do i = 1, n_tallies
tally => tallies(i)
do j = 1, tally % n_filters
filter => tally % filters(j)
filter => tally % filters(j)
if (filter % type == FILTER_DISTRIBCELL) then
if (.not. cell_list % contains(filter % int_bins(1))) then
call cell_list % add(filter % int_bins(1))
end if
end if
end if
end do
end do
! Compute the number of unique universes containing these distribcells
! to determine the number of offset tables to allocate
do i = 1, n_universes
@ -1059,7 +1069,7 @@ contains
end if
end do
end do
! Allocate the list of offset tables for each unique universe
allocate(univ_list(n_maps))
@ -1077,34 +1087,34 @@ contains
univ => universes(i)
do j = 1, univ % n_cells
if (cell_list % contains(univ % cells(j))) then
! Loop over all tallies
! Loop over all tallies
do l = 1, n_tallies
tally => tallies(l)
do m = 1, tally % n_filters
filter => tally % filters(m)
! Loop over only distribcell filters
! If filter points to cell we just found, set offset index
if (filter % type == FILTER_DISTRIBCELL) then
if (filter % type == FILTER_DISTRIBCELL) then
if (filter % int_bins(1) == univ % cells(j)) then
filter % offset = k
end if
end if
end do
end do
end do
univ_list(k) = univ % id
k = k + 1
end if
end do
end do
! Allocate the offset tables for lattices
! Allocate the offset tables for lattices
do i = 1, n_lattices
lat => lattices(i) % obj

70
src/output.F90
View file

@ -258,7 +258,6 @@ contains
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
class(Lattice), pointer :: l => null()
type(LocalCoord), pointer :: coord => null()
! display type of particle
select case (p % type)
@ -273,39 +272,34 @@ contains
end select
! loop through each level of universes
coord => p % coord0
i = 0
do while(associated(coord))
do i = 1, p % n_coord
! Print level
write(ou,*) ' Level ' // trim(to_str(i))
write(ou,*) ' Level ' // trim(to_str(i - 1))
! Print cell for this level
if (coord % cell /= NONE) then
c => cells(coord % cell)
if (p % coord(i) % cell /= NONE) then
c => cells(p % coord(i) % cell)
write(ou,*) ' Cell = ' // trim(to_str(c % id))
end if
! Print universe for this level
if (coord % universe /= NONE) then
u => universes(coord % universe)
if (p % coord(i) % universe /= NONE) then
u => universes(p % coord(i) % universe)
write(ou,*) ' Universe = ' // trim(to_str(u % id))
end if
! Print information on lattice
if (coord % lattice /= NONE) then
l => lattices(coord % lattice) % obj
if (p % coord(i) % lattice /= NONE) then
l => lattices(p % coord(i) % lattice) % obj
write(ou,*) ' Lattice = ' // trim(to_str(l % id))
write(ou,*) ' Lattice position = (' // trim(to_str(&
p % coord % lattice_x)) // ',' // trim(to_str(&
p % coord % lattice_y)) // ')'
p % coord(i) % lattice_x)) // ',' // trim(to_str(&
p % coord(i) % lattice_y)) // ')'
end if
! Print local coordinates
write(ou,'(1X,A,3ES12.4)') ' xyz = ', coord % xyz
write(ou,'(1X,A,3ES12.4)') ' uvw = ', coord % uvw
coord => coord % next
i = i + 1
write(ou,'(1X,A,3ES12.4)') ' xyz = ', p % coord(i) % xyz
write(ou,'(1X,A,3ES12.4)') ' uvw = ', p % coord(i) % uvw
end do
! Print surface
@ -2181,9 +2175,9 @@ contains
end select
end function get_label
!===============================================================================
! FIND_OFFSET uses a given map number, a target cell ID, and a target offset
! FIND_OFFSET uses a given map number, a target cell ID, and a target offset
! to build a string which is the path from the base universe to the target cell
! with the given offset
!===============================================================================
@ -2196,7 +2190,7 @@ contains
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(100) :: path ! Path to offset
integer :: i, j ! Index over cells
integer :: k, l, m ! Indices in lattice
integer :: old_k, old_l, old_m ! Previous indices in lattice
@ -2212,7 +2206,7 @@ contains
class(Lattice), pointer :: lat ! Pointer to current lattice
n = univ % n_cells
! Write to the geometry stack
if (univ%id == 0) then
path = trim(path) // to_str(univ%id)
@ -2223,31 +2217,31 @@ contains
! Look through all cells in this universe
do i = 1, n
cell_index = univ % cells(i)
cell_index = univ % cells(i)
c => cells(cell_index)
! If the cell ID matches the goal and the offset matches final,
! write to the geometry stack
if (cell_dict % get_key(c % id) == goal .AND. offset == final) then
path = trim(path) // "->" // to_str(c%id)
return
end if
end do
! Find the fill cell or lattice cell that we need to enter
do i = 1, n
later_cell = .false.
cell_index = univ % cells(i)
cell_index = univ % cells(i)
c => cells(cell_index)
this_cell = .false.
this_cell = .false.
! If we got here, we still think the target is in this universe
! or further down, but it's not this exact cell.
! Compare offset to next cell to see if we should enter this cell
! or further down, but it's not this exact cell.
! Compare offset to next cell to see if we should enter this cell
if (i /= n) then
do j = i+1, n
@ -2260,8 +2254,8 @@ contains
cycle
end if
! Break loop once we've found the next cell with an offset
exit
! Break loop once we've found the next cell with an offset
exit
end do
! Ensure we didn't just end the loop by iteration
@ -2278,13 +2272,13 @@ contains
else
lat => lattices(c % fill) % obj
temp_offset = lat % offset(map, 1, 1, 1)
end if
end if
! If the final offset is in the range of offset - temp_offset+offset
! then the goal is in this cell
if (final < temp_offset + offset) then
this_cell = .true.
end if
end if
end if
end if
@ -2341,7 +2335,7 @@ contains
! Loop over lattice coordinates
do k = 1, n_x
do l = 1, n_y
do m = 1, n_z
do m = 1, n_z
if (final >= lat % offset(map, k, l, m) + offset) then
if (k == n_x .and. l == n_y .and. m == n_z) then
@ -2364,14 +2358,14 @@ contains
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
next_univ => universes(lat % universes(old_k, old_l, old_m))
path = trim(path) // "(" // trim(to_str(old_k)) // &
"," // trim(to_str(old_l)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
return
end if
end do
end do
end do
@ -2429,7 +2423,7 @@ contains
end if
end if
end do
end do
end subroutine find_offset
end module output

59
src/particle_header.F90
View file

@ -26,9 +26,8 @@ module particle_header
! Is this level rotated?
logical :: rotated = .false.
! Pointer to next (more local) set of coordinates
type(LocalCoord), pointer :: next => null()
contains
procedure :: reset => reset_coord
end type LocalCoord
!===============================================================================
@ -42,8 +41,8 @@ module particle_header
integer :: type ! Particle type (n, p, e, etc)
! Particle coordinates
type(LocalCoord), pointer :: coord0 => null() ! coordinates on universe 0
type(LocalCoord), pointer :: coord => null() ! coordinates on lowest universe
integer :: n_coord ! number of current coordinates
type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels
! Other physical data
real(8) :: wgt ! particle weight
@ -87,26 +86,6 @@ module particle_header
contains
!===============================================================================
! 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 this coord
deallocate(coord)
end if
end subroutine deallocate_coord
!===============================================================================
! INITIALIZE_PARTICLE sets default attributes for a particle from the source
! bank
@ -137,9 +116,8 @@ contains
this % fission = .false.
! Set up base level coordinates
allocate(this % coord0)
this % coord0 % universe = BASE_UNIVERSE
this % coord => this % coord0
this % coord(1) % universe = BASE_UNIVERSE
this % n_coord = 1
end subroutine initialize_particle
@ -150,13 +128,30 @@ contains
subroutine clear_particle(this)
class(Particle) :: this
integer :: i
! remove any coordinate levels
call deallocate_coord(this % coord0)
! Make sure coord pointer is nullified
nullify(this % coord)
do i = 1, MAX_COORD
call this % coord(i) % reset()
end do
end subroutine clear_particle
!===============================================================================
! RESET_COORD
!===============================================================================
elemental subroutine reset_coord(this)
class(LocalCoord), intent(inout) :: this
this % cell = NONE
this % universe = NONE
this % lattice = NONE
this % lattice_x = NONE
this % lattice_y = NONE
this % lattice_z = NONE
this % rotated = .false.
end subroutine reset_coord
end module particle_header

8
src/particle_restart.F90
View file

@ -89,13 +89,13 @@ contains
call pr % read_data(p % id, 'id')
call pr % read_data(p % wgt, 'weight')
call pr % read_data(p % E, 'energy')
call pr % read_data(p % coord % xyz, 'xyz', length=3)
call pr % read_data(p % coord % uvw, 'uvw', length=3)
call pr % read_data(p % coord(1) % xyz, 'xyz', length=3)
call pr % read_data(p % coord(1) % uvw, 'uvw', length=3)
! Set particle last attributes
p % last_wgt = p % wgt
p % last_xyz = p % coord % xyz
p % last_uvw = p % coord % uvw
p % last_xyz = p % coord(1) % xyz
p % last_uvw = p % coord(1) % uvw
p % last_E = p % E
! Close hdf5 file

28
src/physics.F90
View file

@ -37,7 +37,7 @@ contains
! Store pre-collision particle properties
p % last_wgt = p % wgt
p % last_E = p % E
p % last_uvw = p % coord0 % uvw
p % last_uvw = p % coord(1) % uvw
! Add to collision counter for particle
p % n_collision = p % n_collision + 1
@ -253,19 +253,19 @@ contains
p % last_wgt = p % wgt
! Score implicit absorption estimate of keff
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % absorb_wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
if (run_mode == MODE_EIGENVALUE) then
global_tally_absorption = global_tally_absorption + p % absorb_wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
end if
else
! See if disappearance reaction happens
if (micro_xs(i_nuclide) % absorption > &
prn() * micro_xs(i_nuclide) % total) then
! Score absorption estimate of keff
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
if (run_mode == MODE_EIGENVALUE) then
global_tally_absorption = global_tally_absorption + p % wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
end if
p % alive = .false.
p % event = EVENT_ABSORB
@ -332,7 +332,7 @@ contains
if (micro_xs(i_nuclide) % index_sab /= NONE) then
! S(a,b) scattering
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, &
p % E, p % coord0 % uvw, p % mu)
p % E, p % coord(1) % uvw, p % mu)
else
! get pointer to elastic scattering reaction
@ -340,7 +340,7 @@ contains
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, rxn, &
p % E, p % coord0 % uvw, p % mu, p % wgt)
p % E, p % coord(1) % uvw, p % mu, p % wgt)
end if
p % event_MT = ELASTIC
@ -382,7 +382,7 @@ contains
end do
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, rxn, p % E, p % coord0 % uvw, &
call inelastic_scatter(nuc, rxn, p % E, p % coord(1) % uvw, &
p % mu, p % wgt)
p % event_MT = rxn % MT
@ -1074,7 +1074,7 @@ contains
if (ufs) then
! Determine indices on ufs mesh for current location
call get_mesh_indices(ufs_mesh, p % coord0 % xyz, ijk, in_mesh)
call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")
@ -1112,7 +1112,7 @@ contains
p % fission = .true. ! Fission neutrons will be banked
do i = int(n_bank,4) + 1, int(min(n_bank + nu, int(size(fission_bank),8)),4)
! Bank source neutrons by copying particle data
fission_bank(i) % xyz = p % coord0 % xyz
fission_bank(i) % xyz = p % coord(1) % xyz
! Set weight of fission bank site
fission_bank(i) % wgt = ONE/weight

58
src/plot.F90
View file

@ -7,7 +7,7 @@ module plot
use global
use mesh, only: get_mesh_indices
use output, only: write_message
use particle_header, only: deallocate_coord, Particle, LocalCoord
use particle_header, only: Particle, LocalCoord
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
@ -57,26 +57,18 @@ contains
integer, intent(out) :: rgb(3)
integer, intent(out) :: id
integer :: j
logical :: found_cell
integer :: level
type(Cell), pointer :: c => null()
type(LocalCoord), pointer :: coord => null()
type(Cell), pointer :: c
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
p % n_coord = 1
call find_cell(p, found_cell)
j = p % n_coord
if (check_overlaps) call check_cell_overlap(p)
! Loop through universes and stop on any specified level
level = 0
coord => p % coord0
do
if (level == pl % level) exit
if (.not. associated(coord % next)) exit
coord => coord % next
level = level + 1
end do
! Set coordinate level if specified
if (pl % level >= 0) j = pl % level + 1
if (.not. found_cell) then
! If no cell, revert to default color
@ -85,7 +77,7 @@ contains
else
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
c => cells(coord % cell)
c => cells(p % coord(j) % cell)
if (c % material == MATERIAL_VOID) then
! By default, color void cells white
rgb = 255
@ -100,8 +92,8 @@ contains
end if
else if (pl % color_by == PLOT_COLOR_CELLS) then
! Assign color based on cell
rgb = pl % colors(coord % cell) % rgb
id = cells(coord % cell) % id
rgb = pl % colors(p % coord(j) % cell) % rgb
id = cells(p % coord(j) % cell) % id
else
rgb = 0
id = -1
@ -160,9 +152,9 @@ contains
! allocate and initialize particle
call p % initialize()
p % coord % xyz = xyz
p % coord % uvw = [ HALF, HALF, HALF ]
p % coord % universe = BASE_UNIVERSE
p % coord(1) % xyz = xyz
p % coord(1) % uvw = [ HALF, HALF, HALF ]
p % coord(1) % universe = BASE_UNIVERSE
do y = 1, img % height
call progress % set_value(dble(y)/dble(img % height)*100)
@ -175,12 +167,12 @@ contains
call set_pixel(img, x-1, y-1, rgb(1), rgb(2), rgb(3))
! Advance pixel in first direction
p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel
p % coord(1) % xyz(in_i) = p % coord(1) % xyz(in_i) + in_pixel
end do
! Advance pixel in second direction
p % coord0 % xyz(in_i) = xyz(in_i)
p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel
p % coord(1) % xyz(in_i) = xyz(in_i)
p % coord(1) % xyz(out_i) = p % coord(1) % xyz(out_i) - out_pixel
end do
! Draw tally mesh boundaries on the image if requested
@ -367,9 +359,9 @@ contains
! allocate and initialize particle
call p % initialize()
p % coord0 % xyz = ll
p % coord0 % uvw = [ HALF, HALF, HALF ]
p % coord0 % universe = BASE_UNIVERSE
p % coord(1) % xyz = ll
p % coord(1) % uvw = [ HALF, HALF, HALF ]
p % coord(1) % universe = BASE_UNIVERSE
! Open binary plot file for writing
open(UNIT=UNIT_PLOT, FILE=pl % path_plot, STATUS='replace', &
@ -393,20 +385,20 @@ contains
write(UNIT_PLOT) id
! advance particle in z direction
p % coord0 % xyz(3) = p % coord0 % xyz(3) + vox(3)
p % coord(1) % xyz(3) = p % coord(1) % xyz(3) + vox(3)
end do
! advance particle in y direction
p % coord0 % xyz(2) = p % coord0 % xyz(2) + vox(2)
p % coord0 % xyz(3) = ll(3)
p % coord(1) % xyz(2) = p % coord(1) % xyz(2) + vox(2)
p % coord(1) % xyz(3) = ll(3)
end do
! advance particle in y direction
p % coord0 % xyz(1) = p % coord0 % xyz(1) + vox(1)
p % coord0 % xyz(2) = ll(2)
p % coord0 % xyz(3) = ll(3)
p % coord(1) % xyz(1) = p % coord(1) % xyz(1) + vox(1)
p % coord(1) % xyz(2) = ll(2)
p % coord(1) % xyz(3) = ll(3)
end do

12
src/source.F90
View file

@ -132,8 +132,8 @@ contains
site % xyz = p_min + r*(p_max - p_min)
! Fill p with needed data
p % coord0 % xyz = site % xyz
p % coord0 % uvw = [ ONE, ZERO, ZERO ]
p % coord(1) % xyz = site % xyz
p % coord(1) % uvw = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
@ -161,8 +161,8 @@ contains
site % xyz = p_min + r*(p_max - p_min)
! Fill p with needed data
p % coord0 % xyz = site % xyz
p % coord0 % uvw = [ ONE, ZERO, ZERO ]
p % coord(1) % xyz = site % xyz
p % coord(1) % uvw = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
@ -306,8 +306,8 @@ contains
! copy attributes from source bank site
p % wgt = src % wgt
p % last_wgt = src % wgt
p % coord % xyz = src % xyz
p % coord % uvw = src % uvw
p % coord(1) % xyz = src % xyz
p % coord(1) % uvw = src % uvw
p % last_xyz = src % xyz
p % last_uvw = src % uvw
p % E = src % E

68
src/tally.F90
View file

@ -11,7 +11,7 @@ module tally
mesh_intersects_2d, mesh_intersects_3d
use mesh_header, only: StructuredMesh
use output, only: header
use particle_header, only: LocalCoord, Particle, deallocate_coord
use particle_header, only: LocalCoord, Particle
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyResult, TallyMapItem, TallyMapElement
@ -487,7 +487,7 @@ contains
if (t % estimator == ESTIMATOR_ANALOG) then
uvw = p % last_uvw
else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
uvw = p % coord0 % uvw
uvw = p % coord(1) % uvw
end if
! Find the order for a collection of requested moments
! and store the moment contribution of each
@ -909,7 +909,6 @@ contains
type(TallyObject), pointer :: t
type(StructuredMesh), pointer :: m
type(Material), pointer :: mat
type(LocalCoord), pointer :: coord
t => tallies(i_tally)
matching_bins(1:t%n_filters) = 1
@ -918,8 +917,8 @@ contains
! CHECK IF THIS TRACK INTERSECTS THE MESH
! Copy starting and ending location of particle
xyz0 = p % coord0 % xyz - (d_track - TINY_BIT) * p % coord0 % uvw
xyz1 = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw
xyz0 = p % coord(1) % xyz - (d_track - TINY_BIT) * p % coord(1) % uvw
xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
! Get index for mesh filter
i_filter_mesh = t % find_filter(FILTER_MESH)
@ -940,8 +939,8 @@ contains
end if
! Reset starting and ending location
xyz0 = p % coord0 % xyz - d_track * p % coord0 % uvw
xyz1 = p % coord0 % xyz
xyz0 = p % coord(1) % xyz - d_track * p % coord(1) % uvw
xyz1 = p % coord(1) % xyz
! =========================================================================
! CHECK FOR SCORING COMBINATION FOR FILTERS OTHER THAN MESH
@ -953,7 +952,7 @@ contains
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
p % coord(p % n_coord) % universe, i_tally)
case (FILTER_MATERIAL)
matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
@ -961,15 +960,12 @@ contains
case (FILTER_CELL)
! determine next cell bin
coord => p % coord0
do while(associated(coord))
do j = 1, p % n_coord
position(FILTER_CELL) = 0
matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
p % coord(j) % cell, i_tally)
if (matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
@ -1009,7 +1005,7 @@ contains
n_cross = sum(abs(ijk1(:m % n_dimension) - ijk0(:m % n_dimension))) + 1
! Copy particle's direction
uvw = p % coord0 % uvw
uvw = p % coord(1) % uvw
! Bounding coordinates
do j = 1, m % n_dimension
@ -1135,12 +1131,12 @@ contains
logical, intent(out) :: found_bin
integer :: i ! loop index for filters
integer :: j
integer :: n ! number of bins for single filter
integer :: offset ! offset for distribcell
real(8) :: E ! particle energy
type(TallyObject), pointer :: t
type(StructuredMesh), pointer :: m
type(LocalCoord), pointer :: coord
found_bin = .true.
t => tallies(i_tally)
@ -1154,13 +1150,13 @@ contains
m => meshes(t % filters(i) % int_bins(1))
! Determine if we're in the mesh first
call get_mesh_bin(m, p % coord0 % xyz, matching_bins(i))
call get_mesh_bin(m, p % coord(1) % xyz, matching_bins(i))
case (FILTER_UNIVERSE)
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
p % coord(p % n_coord) % universe, i_tally)
case (FILTER_MATERIAL)
if (p % material /= MATERIAL_VOID) then
@ -1170,37 +1166,39 @@ contains
case (FILTER_CELL)
! determine next cell bin
coord => p % coord0
do while(associated(coord))
do j = 1, p % n_coord
position(FILTER_CELL) = 0
matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
p % coord(j) % cell, i_tally)
if (matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_DISTRIBCELL)
! determine next distribcell bin
matching_bins(i) = NO_BIN_FOUND
coord => p % coord0
offset = 0
do while(associated(coord))
if (cells(coord % cell) % type == CELL_FILL) then
offset = offset + cells(coord % cell) % &
do j = 1, p % n_coord
if (cells(p % coord(j) % cell) % type == CELL_FILL) then
offset = offset + cells(p % coord(j) % cell) % &
offset(t % filters(i) % offset)
elseif(cells(coord % cell) % type == CELL_LATTICE) then
offset = offset + lattices(coord % next % lattice) % obj % &
offset(t % filters(i) % offset, coord % next % lattice_x, &
coord % next % lattice_y, coord % next % lattice_z)
elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then
if (lattices(p % coord(j + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(j + 1) % lattice_x, &
p % coord(j + 1) % lattice_y, &
p % coord(j + 1) % lattice_z])) then
offset = offset + lattices(p % coord(j + 1) % lattice) % obj % &
offset(t % filters(i) % offset, &
p % coord(j + 1) % lattice_x, &
p % coord(j + 1) % lattice_y, &
p % coord(j + 1) % lattice_z)
end if
end if
if (t % filters(i) % int_bins(1) == coord % cell) then
if (t % filters(i) % int_bins(1) == p % coord(j) % cell) then
matching_bins(i) = offset + 1
exit
end if
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
@ -1296,7 +1294,7 @@ contains
TALLY_LOOP: do i = 1, active_current_tallies % size()
! Copy starting and ending location of particle
xyz0 = p % last_xyz
xyz1 = p % coord0 % xyz
xyz1 = p % coord(1) % xyz
! Get pointer to tally
i_tally = active_current_tallies % get_item(i)
@ -1328,7 +1326,7 @@ contains
end if
! Copy particle's direction
uvw = p % coord0 % uvw
uvw = p % coord(1) % uvw
! determine incoming energy bin
j = t % find_filter(FILTER_ENERGYIN)

2
src/track_output.F90
View file

@ -46,7 +46,7 @@ contains
end if
! Write current coordinates into the newest column.
coords(:, n_tracks) = p % coord0 % xyz
coords(:, n_tracks) = p % coord(1) % xyz
end subroutine write_particle_track
!===============================================================================

53
src/tracking.F90
View file

@ -1,5 +1,6 @@
module tracking
use constants, only: MODE_EIGENVALUE
use cross_section, only: calculate_xs
use error, only: fatal_error, warning
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
@ -28,6 +29,8 @@ contains
type(Particle), intent(inout) :: p
integer :: j ! coordinate level
integer :: next_level ! next coordinate level to check
integer :: surface_crossed ! surface which particle is on
integer :: lattice_translation(3) ! in-lattice translation vector
integer :: last_cell ! most recent cell particle was in
@ -36,7 +39,6 @@ contains
real(8) :: d_collision ! sampled distance to collision
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
type(LocalCoord), pointer :: coord
! Display message if high verbosity or trace is on
if (verbosity >= 9 .or. trace) then
@ -45,7 +47,7 @@ contains
! If the cell hasn't been determined based on the particle's location,
! initiate a search for the current cell
if (p % coord % cell == NONE) then
if (p % coord(p % n_coord) % cell == NONE) then
call find_cell(p, found_cell)
! Particle couldn't be located
@ -54,7 +56,7 @@ contains
end if
! set birth cell attribute
p % cell_born = p % coord % cell
p % cell_born = p % coord(p % n_coord) % cell
end if
! Initialize number of events to zero
@ -87,7 +89,7 @@ contains
! Find the distance to the nearest boundary
call distance_to_boundary(p, d_boundary, surface_crossed, &
&lattice_translation)
lattice_translation, next_level)
! Sample a distance to collision
if (material_xs % total == ZERO) then
@ -100,10 +102,8 @@ contains
distance = min(d_boundary, d_collision)
! Advance particle
coord => p % coord0
do while (associated(coord))
coord % xyz = coord % xyz + distance * coord % uvw
coord => coord % next
do j = 1, p % n_coord
p % coord(j) % xyz = p % coord(j) % xyz + distance * p % coord(j) % uvw
end do
! Score track-length tallies
@ -111,17 +111,18 @@ contains
call score_tracklength_tally(p, distance)
! Score track-length estimate of k-eff
!$omp atomic
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * &
material_xs % nu_fission
if (run_mode == MODE_EIGENVALUE) then
global_tally_tracklength = global_tally_tracklength + p % wgt * &
distance * material_xs % nu_fission
end if
if (d_collision > d_boundary) then
! ====================================================================
! PARTICLE CROSSES SURFACE
last_cell = p % coord % cell
p % coord % cell = NONE
if (next_level > 0) p % n_coord = next_level
last_cell = p % coord(p % n_coord) % cell
p % coord(p % n_coord) % cell = NONE
if (any(lattice_translation /= 0)) then
! Particle crosses lattice boundary
p % surface = NONE
@ -138,10 +139,10 @@ contains
! PARTICLE HAS COLLISION
! Score collision estimate of keff
!$omp atomic
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + p % wgt * &
material_xs % nu_fission / material_xs % total
if (run_mode == MODE_EIGENVALUE) then
global_tally_collision = global_tally_collision + p % wgt * &
material_xs % nu_fission / material_xs % total
end if
! score surface current tallies -- this has to be done before the collision
! since the direction of the particle will change and we need to use the
@ -168,7 +169,7 @@ contains
p % fission = .false.
! Save coordinates for tallying purposes
p % last_xyz = p % coord0 % xyz
p % last_xyz = p % coord(1) % xyz
! Set last material to none since cross sections will need to be
! re-evaluated
@ -176,19 +177,15 @@ contains
! Set all uvws to base level -- right now, after a collision, only the
! base level uvws are changed
coord => p % coord0
do while(associated(coord % next))
if (coord % next % rotated) then
do j = 1, p % n_coord - 1
if (p % coord(j + 1) % rotated) then
! If next level is rotated, apply rotation matrix
coord % next % uvw = matmul(cells(coord % cell) % &
rotation_matrix, coord % uvw)
p % coord(j + 1) % uvw = matmul(cells(p % coord(j) % cell) % &
rotation_matrix, p % coord(j) % uvw)
else
! Otherwise, copy this level's direction
coord % next % uvw = coord % uvw
p % coord(j + 1) % uvw = p % coord(j) % uvw
end if
! Advance coordinate level
coord => coord % next
end do
end if

10
tests/test_filter_distribcell/case-1/geometry.xml
View file

@ -4,21 +4,23 @@
<cell id="1" fill="5" surfaces="1 -2 3 -4" />
<cell id="201" universe="2" material="1" surfaces="-5" />
<cell id="202" universe="2" material="2" surfaces="5" />
<cell id="301" universe="3" material="2"/>
<lattice id="5">
<dimension>2 2</dimension>
<lower_left>-2.0 -2.0</lower_left>
<pitch>2.0 2.0</pitch>
<outer>3</outer>
<universes>
2 2
2 2
</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="1" type="x-plane" coeffs="-3.0" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="3.0" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-3.0" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="3.0" boundary="vacuum" />
<surface id="5" type="z-cylinder" coeffs="0.0 0.0 0.3" />
</geometry>

20
tests/test_filter_distribcell/case-1/results_true.dat
View file

@ -1,14 +1,14 @@
k-combined:
0.000000E+00 0.000000E+00
tally 1:
8.565980E-03
7.337601E-05
8.045879E-03
6.473617E-05
9.027116E-03
8.148882E-05
7.326285E-03
5.367445E-05
1.440759E-02
2.075788E-04
1.222930E-02
1.495558E-04
1.407292E-02
1.980471E-04
1.034365E-02
1.069911E-04
tally 2:
3.296526E-02
1.086708E-03
5.105347E-02
2.606457E-03

12
tests/test_filter_universe/results_true.dat
View file

@ -1,11 +1,11 @@
k-combined:
1.093844E+00 1.626801E-02
tally 1:
6.144371E+01
7.795909E+02
7.330533E+00
1.137816E+01
4.850651E+01
4.933509E+02
6.369043E+01
8.385422E+02
7.330762E+00
1.137874E+01
5.011385E+01
5.251276E+02
5.132453E+00
5.801019E+00

69
tests/travis_install.sh
View file

@ -1,38 +1,39 @@
#!/bin/sh
#!/bin/bash
set -ev
# Build MPICH and HDF5 for rest of debug tests
if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
# Build MPICH
wget -q http://www.mpich.org/static/downloads/3.1.3/mpich-3.1.3.tar.gz
tar -xzvf mpich-3.1.3.tar.gz >/dev/null 2>&1
cd mpich-3.1.3
./configure --prefix=$PWD/../mpich_install -q
make -j >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
# Build PHDF5
wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
mv hdf5-1.8.15 phdf5-1.8.15; cd phdf5-1.8.15
CC=$PWD/../mpich_install/bin/mpicc FC=$PWD/../mpich_install/bin/mpif90 \
./configure \
--prefix=$PWD/../phdf5_install -q --enable-fortran \
--enable-fortran2003 --enable-parallel
make -j >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
# Build HDF5
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
cd hdf5-1.8.15
CC=gcc FC=gfortran ./configure --prefix=$PWD/../hdf5_install -q \
--enable-fortran --enable-fortran2003
make -j >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
# Build MPICH
if [[ ! -e $HOME/mpich_install/bin/mpiexec ]]; then
wget -q http://www.mpich.org/static/downloads/3.1.3/mpich-3.1.3.tar.gz
tar -xzvf mpich-3.1.3.tar.gz >/dev/null 2>&1
cd mpich-3.1.3
./configure --prefix=$HOME/mpich_install -q
make >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
fi
# Build PHDF5
if [[ ! -e $HOME/phdf5_install/bin/h5pfc ]]; then
wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
mv hdf5-1.8.15 phdf5-1.8.15; cd phdf5-1.8.15
CC=$HOME/mpich_install/bin/mpicc FC=$HOME/mpich_install/bin/mpif90 \
./configure \
--prefix=$HOME/phdf5_install -q --enable-fortran \
--enable-fortran2003 --enable-parallel
make >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
fi
# Build HDF5
if [[ ! -e $HOME/hdf5_install/bin/h5fc ]]; then
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
cd hdf5-1.8.15
CC=gcc FC=gfortran ./configure --prefix=$HOME/hdf5_install -q \
--enable-fortran --enable-fortran2003
make -j >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
fi