diff --git a/.travis.yml b/.travis.yml index 21d110b40..7af617f85 100644 --- a/.travis.yml +++ b/.travis.yml @@ -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 diff --git a/CMakeLists.txt b/CMakeLists.txt index e732d9f79..3bbc353f6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 #=============================================================================== diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index 19e44ec7f..25c58ae77 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -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 ------------------------- diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 848c6ecaa..7d3cda173 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -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 ++++++++++++++++++ diff --git a/openmc/geometry.py b/openmc/geometry.py index f1c038403..2ea9705fa 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -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") diff --git a/openmc/universe.py b/openmc/universe.py index 62055049d..84ce3696e 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -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) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 07b6e472b..aefb597d0 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -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() diff --git a/src/geometry.F90 b/src/geometry.F90 index d88265c53..39d95817b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 4a74bbeee..f7d9f184e 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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(:) diff --git a/src/initialize.F90 b/src/initialize.F90 index 578253203..99a66e2ce 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index 009e5f7e2..7b5e98564 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/particle_header.F90 b/src/particle_header.F90 index bf03c5635..7c2586c3e 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -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 diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 6ff5d63c2..0ed25fdc8 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index 04be42fd9..6f9457ff6 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/plot.F90 b/src/plot.F90 index 01f4a6021..c6341e3a9 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -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 diff --git a/src/source.F90 b/src/source.F90 index 6dd73422e..7cfdcf655 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -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 diff --git a/src/tally.F90 b/src/tally.F90 index dc9169fd5..e78e58f28 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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) diff --git a/src/track_output.F90 b/src/track_output.F90 index 2159003ce..fba699e8e 100644 --- a/src/track_output.F90 +++ b/src/track_output.F90 @@ -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 !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index c14d85b5d..e5faeb209 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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 diff --git a/tests/test_filter_distribcell/case-1/geometry.xml b/tests/test_filter_distribcell/case-1/geometry.xml index ddbd23afe..507559615 100644 --- a/tests/test_filter_distribcell/case-1/geometry.xml +++ b/tests/test_filter_distribcell/case-1/geometry.xml @@ -4,21 +4,23 @@ + 2 2 -2.0 -2.0 2.0 2.0 + 3 2 2 2 2 - - - - + + + + diff --git a/tests/test_filter_distribcell/case-1/results_true.dat b/tests/test_filter_distribcell/case-1/results_true.dat index fb68a38f9..49bf3ed4e 100644 --- a/tests/test_filter_distribcell/case-1/results_true.dat +++ b/tests/test_filter_distribcell/case-1/results_true.dat @@ -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 diff --git a/tests/test_filter_universe/results_true.dat b/tests/test_filter_universe/results_true.dat index 1442d52b7..c63da0083 100644 --- a/tests/test_filter_universe/results_true.dat +++ b/tests/test_filter_universe/results_true.dat @@ -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 diff --git a/tests/travis_install.sh b/tests/travis_install.sh index 4976ad092..c01a980c0 100755 --- a/tests/travis_install.sh +++ b/tests/travis_install.sh @@ -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