From 077431f1e5a1a6211dde36b1e83307b47ee204c7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Aug 2013 09:18:02 -0500 Subject: [PATCH 01/16] Added OpenMP shared memory parallelism based on changes in threading branch. Works for ifort, but not for gfortran. --- src/Makefile | 24 ++++++++++++++++-- src/ace.F90 | 2 ++ src/cross_section.F90 | 17 ++++++++----- src/eigenvalue.F90 | 52 ++++++++++++++++++++++++++++++++++++++- src/geometry.F90 | 35 +++++++++++++++----------- src/global.F90 | 11 +++++++++ src/initialize.F90 | 26 ++++++++++++++++++-- src/physics.F90 | 32 +++++++++++++++--------- src/random_lcg.F90 | 2 ++ src/source.F90 | 3 ++- src/tally.F90 | 53 ++++++++++++++++++++++++++++------------ src/tracking.F90 | 6 ++++- src/xml-fortran/Makefile | 2 +- 13 files changed, 210 insertions(+), 55 deletions(-) diff --git a/src/Makefile b/src/Makefile index 6107c37adc..2369e29fba 100644 --- a/src/Makefile +++ b/src/Makefile @@ -17,6 +17,7 @@ DEBUG = no PROFILE = no OPTIMIZE = no MPI = no +OPENMP = no HDF5 = no PETSC = no @@ -179,6 +180,25 @@ else endif endif +# OpenMP for shared-memory parallelism + +ifeq ($(OPENMP),yes) + ifeq ($(COMPILER),intel) + F90FLAGS += -openmp -DOPENMP + LDFLAGS += -openmp + endif + + ifeq ($(COMPILER),gnu) + F90FLAGS += -fopenmp -DOPENMP + LDFLAGS += -fopenmp + endif + + ifeq ($(COMPILER),ibm) + F90FLAGS += -qsmp=omp -WF,-DOPENMP + LDFLAGS += -qsmp=omp + endif +endif + # PETSC for CMFD functionality ifeq ($(PETSC),yes) @@ -225,8 +245,8 @@ endif all: xml-fortran $(program) xml-fortran: - cd xml-fortran; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)" - cd templates; make F90=$(F90) F90FLAGS="$(F90FLAGS)" + cd xml-fortran; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)" LDFLAGS="$(LDFLAGS)" + cd templates; make F90=$(F90) F90FLAGS="$(F90FLAGS)" LDFLAGS="$(LDFLAGS)" $(program): $(objects) $(F90) $(objects) $(templates) $(xml_fort) $(LDFLAGS) -o $@ install: diff --git a/src/ace.F90 b/src/ace.F90 index 2617bcbed4..ab67b3e09d 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -51,7 +51,9 @@ contains ! allocate arrays for ACE table storage and cross section cache allocate(nuclides(n_nuclides_total)) allocate(sab_tables(n_sab_tables)) +!$omp parallel allocate(micro_xs(n_nuclides_total)) +!$omp end parallel ! ========================================================================== ! READ ALL ACE CROSS SECTION TABLES diff --git a/src/cross_section.F90 b/src/cross_section.F90 index be9f0740a1..4c9ffb8c88 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -14,6 +14,7 @@ module cross_section save integer :: union_grid_index +!$omp threadprivate(union_grid_index) contains @@ -32,7 +33,8 @@ contains integer :: j ! index in mat % i_sab_nuclides real(8) :: atom_density ! atom density of a nuclide logical :: check_sab ! should we check for S(a,b) table? - type(Material), pointer :: mat => null() ! current material + type(Material), pointer, save :: mat => null() ! current material +!$omp threadprivate(mat) ! Set all material macroscopic cross sections to zero material_xs % total = ZERO @@ -139,7 +141,8 @@ contains integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid - type(Nuclide), pointer :: nuc => null() + type(Nuclide), pointer, save :: nuc => null() +!$omp threadprivate(nuc) ! Set pointer to nuclide nuc => nuclides(i_nuclide) @@ -258,7 +261,8 @@ contains real(8) :: f ! interp factor on S(a,b) energy grid real(8) :: inelastic ! S(a,b) inelastic cross section real(8) :: elastic ! S(a,b) elastic cross section - type(SAlphaBeta), pointer :: sab => null() + type(SAlphaBeta), pointer, save :: sab => null() +!$omp threadprivate(sab) ! Set flag that S(a,b) treatment should be used for scattering micro_xs(i_nuclide) % index_sab = i_sab @@ -346,9 +350,10 @@ contains real(8) :: capture ! (n,gamma) cross section real(8) :: fission ! fission cross section real(8) :: inelastic ! inelastic cross section - type(UrrData), pointer :: urr => null() - type(Nuclide), pointer :: nuc => null() - type(Reaction), pointer :: rxn => null() + type(UrrData), pointer, save :: urr => null() + type(Nuclide), pointer, save :: nuc => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(urr, nuc, rxn) micro_xs(i_nuclide) % use_ptable = .true. diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 99e758bf38..f52ccef810 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -39,6 +39,7 @@ contains subroutine run_eigenvalue() type(Particle) :: p + integer :: i_work if (master) call header("K EIGENVALUE SIMULATION", level=1) @@ -71,7 +72,9 @@ contains ! ==================================================================== ! LOOP OVER PARTICLES - PARTICLE_LOOP: do current_work = 1, work +!$omp parallel do schedule(static) firstprivate(p) + PARTICLE_LOOP: do i_work = 1, work + current_work = i_work ! grab source particle from bank call get_source_particle(p, current_work) @@ -80,6 +83,7 @@ contains call transport(p) end do PARTICLE_LOOP +!$omp end parallel do ! Accumulate time for transport call time_transport % stop() @@ -157,6 +161,11 @@ contains subroutine finalize_generation() +#ifdef OPENMP + ! Join the fission bank from each thread into one global fission bank + call join_bank_from_threads() +#endif + ! Distribute fission bank across processors evenly call time_bank % start() call synchronize_bank() @@ -808,4 +817,45 @@ contains end subroutine replay_batch_history +#ifdef OPENMP +!=============================================================================== +! JOIN_BANK_FROM_THREADS +!=============================================================================== + + subroutine join_bank_from_threads() + + integer :: total ! total number of fission bank sites + integer :: i ! loop index for threads + + ! Initialize the total number of fission bank sites + total = 0 + +!$omp parallel + + ! Copy thread fission bank sites to one shared copy +!$omp do ordered schedule(static) + do i = 1, n_threads +!$omp ordered + master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank) + total = total + n_bank +!$omp end ordered + end do +!$omp end do + + ! Make sure all threads have made it to this point +!$omp barrier + + ! Now copy the shared fission bank sites back to the master thread's copy. + if (thread_id == 0) then + n_bank = total + fission_bank(1:n_bank) = master_fission_bank(1:n_bank) + else + n_bank = 0 + end if + +!$omp end parallel + + end subroutine join_bank_from_threads +#endif + end module eigenvalue diff --git a/src/geometry.F90 b/src/geometry.F90 index a50e55252a..547ad21d36 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -29,7 +29,8 @@ contains integer :: i_surface ! index in surfaces array (with sign) logical :: specified_sense ! specified sense of surface in list logical :: actual_sense ! sense of particle wrt surface - type(Surface), pointer :: s => null() + type(Surface), pointer, save :: s => null() +!$omp threadprivate(s) SURFACE_LOOP: do i = 1, c % n_surfaces ! Lookup surface @@ -76,9 +77,10 @@ contains integer :: i ! cell loop index on a level 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 + type(Cell), pointer, save :: c => null() ! pointer to cell + type(Universe), pointer, save :: univ => null() ! universe to search in + type(LocalCoord), pointer, save :: coord => null() ! particle coordinate to search on +!$omp threadprivate(c, univ, coord) coord => p % coord0 @@ -139,9 +141,10 @@ contains logical :: use_search_cells ! use cells provided as argument logical :: outside_lattice ! if particle is not inside lattice bounds logical :: lattice_edge ! if particle is on a lattice edge - type(Cell), pointer :: c ! pointer to cell - type(Lattice), pointer :: lat ! pointer to lattice - type(Universe), pointer :: univ ! universe to search in + type(Cell), pointer, save :: c => null() ! pointer to cell + type(Lattice), pointer, save :: lat => null() ! pointer to lattice + type(Universe), pointer, save :: univ => null() ! universe to search in +!$omp threadprivate(c, lat, univ) ! Remove coordinates for any lower levels call deallocate_coord(p % coord % next) @@ -375,7 +378,8 @@ contains real(8) :: norm ! "norm" of surface normal integer :: i_surface ! index in surfaces logical :: found ! particle found in universe? - type(Surface), pointer :: surf => null() + type(Surface), pointer, save :: surf => null() +!$omp threadprivate(surf) i_surface = abs(p % surface) surf => surfaces(i_surface) @@ -658,7 +662,8 @@ contains integer :: n_x, n_y, n_z ! size of lattice real(8) :: x0, y0, z0 ! half width of lattice element logical :: found ! particle found in cell? - type(Lattice), pointer :: lat => null() + type(Lattice), pointer, save :: lat => null() +!$omp threadprivate(lat) lat => lattices(p % coord % lattice) @@ -787,11 +792,12 @@ contains real(8) :: a,b,c,k ! quadratic equation coefficients real(8) :: quad ! discriminant of quadratic equation logical :: on_surface ! is particle on surface? - type(Cell), pointer :: cl => null() - type(Surface), pointer :: surf => null() - type(Lattice), pointer :: lat => null() - type(LocalCoord), pointer :: coord => null() - type(LocalCoord), pointer :: final_coord => null() + type(Cell), pointer, save :: cl => null() + type(Surface), pointer, save :: surf => null() + type(Lattice), pointer, save :: lat => null() + type(LocalCoord), pointer, save :: coord => null() + type(LocalCoord), pointer, save :: final_coord => null() +!$omp threadprivate(cl, surf, lat, coord, final_coord) ! inialize distance to infinity (huge) dist = INFINITY @@ -1562,6 +1568,7 @@ contains ! Increment number of lost particles p % alive = .false. +!$omp atomic n_lost_particles = n_lost_particles + 1 ! Abort the simulation if the maximum number of lost particles has been diff --git a/src/global.F90 b/src/global.F90 index a7a7bfb972..a84e9b00b8 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -159,6 +159,9 @@ module global ! Source and fission bank type(Bank), allocatable, target :: source_bank(:) type(Bank), allocatable, target :: fission_bank(:) +#ifdef OPENMP + type(Bank), allocatable, target :: master_fission_bank(:) +#endif integer(8) :: n_bank ! # of sites in fission bank integer(8) :: bank_first ! index of first particle in bank integer(8) :: bank_last ! index of last particle in bank @@ -205,6 +208,11 @@ module global integer :: MPI_BANK ! MPI datatype for fission bank integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult +#ifdef OPENMP + integer :: n_threads ! number of OpenMP threads + integer :: thread_id ! ID of a given thread +#endif + ! No reduction at end of batch logical :: reduce_tallies = .true. @@ -375,6 +383,9 @@ module global logical :: output_xs = .false. logical :: output_tallies = .true. +!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, & +!$omp& trace, thread_id, current_work) + contains !=============================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index da9194915c..03ee0045a1 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -26,6 +26,10 @@ module initialize use mpi #endif +#ifdef OPENMP + use omp_lib +#endif + #ifdef HDF5 use hdf5_interface use hdf5_summary, only: hdf5_write_summary @@ -761,7 +765,7 @@ contains subroutine allocate_banks() - integer :: alloc_err ! allocation error code + integer :: alloc_err ! allocation error code ! Allocate source bank allocate(source_bank(maxwork), STAT=alloc_err) @@ -772,8 +776,26 @@ contains call fatal_error() end if - ! Allocate fission bank +#ifdef OPENMP + ! If OpenMP is being used, each thread needs its own private fission + ! bank. Since the private fission banks need to be combined at the end of a + ! generation, there is also a 'master_fission_bank' that is used to collect + ! the sites from each thread. + +!$omp parallel + n_threads = omp_get_num_threads() + thread_id = omp_get_thread_num() + + if (thread_id == 0) then + allocate(fission_bank(3*maxwork)) + else + allocate(fission_bank(3*maxwork/n_threads)) + end if +!$omp end parallel + allocate(master_fission_bank(3*maxwork), STAT=alloc_err) +#else allocate(fission_bank(3*maxwork), STAT=alloc_err) +#endif ! Check for allocation errors if (alloc_err /= 0) then diff --git a/src/physics.F90 b/src/physics.F90 index a8325b2d87..500d32477a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -74,7 +74,8 @@ contains integer :: i_nuclide ! index in nuclides array integer :: i_reaction ! index in nuc % reactions array - type(Nuclide), pointer :: nuc => null() + type(Nuclide), pointer, save :: nuc => null() +!$omp threadprivate(nuc) i_nuclide = sample_nuclide(p, 'total ') @@ -129,7 +130,8 @@ contains real(8) :: cutoff real(8) :: atom_density ! atom density of nuclide in atom/b-cm real(8) :: sigma ! microscopic total xs for nuclide - type(Material), pointer :: mat => null() + type(Material), pointer, save :: mat => null() +!$omp threadprivate(mat) ! Get pointer to current material mat => materials(p % material) @@ -191,8 +193,9 @@ contains real(8) :: f real(8) :: prob real(8) :: cutoff - type(Nuclide), pointer :: nuc => null() - type(Reaction), pointer :: rxn => null() + type(Nuclide), pointer, save :: nuc => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(nuc, rxn) ! Get pointer to nuclide nuc => nuclides(i_nuclide) @@ -251,6 +254,7 @@ 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 @@ -260,6 +264,7 @@ contains 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 @@ -306,8 +311,9 @@ contains real(8) :: f real(8) :: prob real(8) :: cutoff - type(Nuclide), pointer :: nuc => null() - type(Reaction), pointer :: rxn => null() + type(Nuclide), pointer, save :: nuc => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(nuc, rxn) ! Get pointer to nuclide and grid index/interpolation factor nuc => nuclides(i_nuclide) @@ -410,7 +416,8 @@ contains real(8) :: v_cm(3) ! velocity of center-of-mass real(8) :: v_t(3) ! velocity of target nucleus real(8) :: uvw_cm(3) ! directional cosines in center-of-mass - type(Nuclide), pointer :: nuc => null() + type(Nuclide), pointer, save :: nuc => null() +!$omp threadprivate(nuc) ! get pointer to nuclide nuc => nuclides(i_nuclide) @@ -487,7 +494,8 @@ contains real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j) real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j) real(8) :: prob ! probability for sampling Bragg edge - type(SAlphaBeta), pointer :: sab => null() + type(SAlphaBeta), pointer, save :: sab => null() +!$omp threadprivate(sab) ! Get pointer to S(a,b) table sab => sab_tables(i_sab) @@ -715,8 +723,9 @@ contains real(8) :: phi ! fission neutron azimuthal angle real(8) :: weight ! weight adjustment for ufs method logical :: in_mesh ! source site in ufs mesh? - type(Nuclide), pointer :: nuc => null() - type(Reaction), pointer :: rxn => null() + type(Nuclide), pointer, save :: nuc => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(nuc, rxn) ! Get pointers nuc => nuclides(i_nuclide) @@ -815,7 +824,8 @@ contains real(8) :: xi ! random number real(8) :: yield ! delayed neutron precursor yield real(8) :: prob ! cumulative probability - type(DistEnergy), pointer :: edist => null() + type(DistEnergy), pointer, save :: edist => null() +!$omp threadprivate(edist) ! Determine total nu nu_t = nu_total(nuc, E) diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 3130d42a53..83ff6409f9 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -15,6 +15,8 @@ module random_lcg integer(8) :: prn_stride ! stride between particles real(8) :: prn_norm ! 2^(-M) +!$omp threadprivate(prn_seed) + public :: prn public :: initialize_prng public :: set_particle_seed diff --git a/src/source.F90 b/src/source.F90 index c951c11ffa..6f2bad2477 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -155,7 +155,8 @@ contains integer(8), intent(in) :: index_source integer(8) :: particle_seed ! unique index for particle - type(Bank), pointer :: src => null() + type(Bank), pointer, save :: src => null() +!$omp threadprivate(src) ! set defaults call p % initialize() diff --git a/src/tally.F90 b/src/tally.F90 index 8550473e4a..49a06140bb 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -23,6 +23,7 @@ module tally ! Tally map positioning array integer :: position(N_FILTER_TYPES - 3) = 0 +!$omp threadprivate(position) contains @@ -54,7 +55,8 @@ contains real(8) :: macro_total ! material macro total xs real(8) :: macro_scatt ! material macro scatt xs logical :: found_bin ! scoring bin found? - type(TallyObject), pointer :: t => null() + type(TallyObject), pointer, save :: t => null() +!$omp threadprivate(t) ! Copy particle's pre- and post-collision weight and angle last_wgt = p % last_wgt @@ -200,6 +202,7 @@ contains ! get the score and tally it score = last_wgt * calc_pn(n, mu) +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score end do @@ -338,6 +341,7 @@ contains end select ! Add score to tally +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score @@ -411,6 +415,7 @@ contains i_filter = sum((t % matching_bins - 1) * t % stride) + 1 ! Add score to tally +!$omp atomic t % results(i_score, i_filter) % value = & t % results(i_score, i_filter) % value + score end do @@ -449,9 +454,10 @@ contains real(8) :: score ! actual score (e.g., flux*xs) real(8) :: atom_density ! atom density of single nuclide in atom/b-cm logical :: found_bin ! scoring bin found? - type(TallyObject), pointer :: t => null() - type(Material), pointer :: mat => null() - type(Reaction), pointer :: rxn => null() + type(TallyObject), pointer, save :: t => null() + type(Material), pointer, save :: mat => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(t, mat, rxn) ! Determine track-length estimate of flux flux = p % wgt * distance @@ -698,6 +704,7 @@ contains score_index = (k - 1)*t % n_score_bins + j ! Add score to tally +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score @@ -742,9 +749,10 @@ contains real(8) :: f ! interpolation factor real(8) :: score ! actual scoring tally value real(8) :: atom_density ! atom density of single nuclide in atom/b-cm - type(TallyObject), pointer :: t => null() - type(Material), pointer :: mat => null() - type(Reaction), pointer :: rxn => null() + type(TallyObject), pointer, save :: t => null() + type(Material), pointer, save :: mat => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(t, mat, rxn) ! Get pointer to tally t => tallies(i_tally) @@ -840,6 +848,7 @@ contains score_index = (i_nuclide - 1)*t % n_score_bins + j ! Add score to tally +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score @@ -938,6 +947,7 @@ contains score_index = n_nuclides_total*t % n_score_bins + j ! Add score to tally +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score @@ -982,10 +992,11 @@ contains logical :: found_bin ! was a scoring bin found? logical :: start_in_mesh ! starting coordinates inside mesh? logical :: end_in_mesh ! ending coordinates inside mesh? - type(TallyObject), pointer :: t => null() - type(StructuredMesh), pointer :: m => null() - type(Material), pointer :: mat => null() - type(LocalCoord), pointer :: coord => null() + type(TallyObject), pointer, save :: t => null() + type(StructuredMesh), pointer, save :: m => null() + type(Material), pointer, save :: mat => null() + type(LocalCoord), pointer, save :: coord => null() +!$omp threadprivate(t, m, mat, coord) t => tallies(i_tally) t % matching_bins = 1 @@ -1245,6 +1256,7 @@ contains score_index = (b - 1)*t % n_score_bins + j ! Add score to tally +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score @@ -1275,9 +1287,10 @@ contains integer :: i ! loop index for filters integer :: n ! number of bins for single filter real(8) :: E ! particle energy - type(TallyObject), pointer :: t => null() - type(StructuredMesh), pointer :: m => null() - type(LocalCoord), pointer :: coord => null() + type(TallyObject), pointer, save :: t => null() + type(StructuredMesh), pointer, save :: m => null() + type(LocalCoord), pointer, save :: coord => null() +!$omp threadprivate(t, m, coord) found_bin = .true. t => tallies(i_tally) @@ -1403,8 +1416,9 @@ contains logical :: x_same ! same starting/ending x index (i) logical :: y_same ! same starting/ending y index (j) logical :: z_same ! same starting/ending z index (k) - type(TallyObject), pointer :: t => null() - type(StructuredMesh), pointer :: m => null() + type(TallyObject), pointer, save :: t => null() + type(StructuredMesh), pointer, save :: m => null() +!$omp threadprivate(t, m) TALLY_LOOP: do i = 1, active_current_tallies % size() ! Copy starting and ending location of particle @@ -1475,6 +1489,7 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if @@ -1487,6 +1502,7 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if @@ -1503,6 +1519,7 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if @@ -1515,6 +1532,7 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if @@ -1531,6 +1549,7 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if @@ -1543,6 +1562,7 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if @@ -1667,6 +1687,7 @@ contains end if ! Add to surface current tally +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if diff --git a/src/tracking.F90 b/src/tracking.F90 index aee90ea52c..6080b2bae1 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -32,7 +32,8 @@ 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 => null() + type(LocalCoord), pointer, save :: coord => null() +!$omp threadprivate(coord) ! Display message if high verbosity or trace is on if (verbosity >= 9 .or. trace) then @@ -59,6 +60,7 @@ contains n_event = 0 ! Add paricle's starting weight to count for normalizing tallies later +!$omp atomic total_weight = total_weight + p % wgt ! Force calculation of cross-sections by setting last energy to zero @@ -99,6 +101,7 @@ 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 @@ -125,6 +128,7 @@ 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 diff --git a/src/xml-fortran/Makefile b/src/xml-fortran/Makefile index c0514533a3..6bba7a1486 100644 --- a/src/xml-fortran/Makefile +++ b/src/xml-fortran/Makefile @@ -22,7 +22,7 @@ endif all: $(reader) $(reader): $(objects) - $(F90) $(objects) -o $@ + $(F90) $(objects) $(LDFLAGS) -o $@ clean: @rm -f *.o *.mod $(reader) neat: From 31db2828022caeabd60ec8a04561ba8ceab4edf7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Aug 2013 17:31:56 -0500 Subject: [PATCH 02/16] Changed all atomics to criticals. gfortran now produce correct k-effective. --- src/physics.F90 | 6 ++++-- src/tally.F90 | 42 ++++++++++++++++++++++++++++-------------- src/tracking.F90 | 9 ++++++--- 3 files changed, 38 insertions(+), 19 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 500d32477a..e475655e3d 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -254,20 +254,22 @@ contains p % last_wgt = p % wgt ! Score implicit absorption estimate of keff -!$omp atomic +!$omp critical global_tallies(K_ABSORPTION) % value = & global_tallies(K_ABSORPTION) % value + p % absorb_wgt * & micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption +!$omp end critical 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 +!$omp critical global_tallies(K_ABSORPTION) % value = & global_tallies(K_ABSORPTION) % value + p % wgt * & micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption +!$omp end critical p % alive = .false. p % event = EVENT_ABSORB diff --git a/src/tally.F90 b/src/tally.F90 index 49a06140bb..2fafacf320 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -202,9 +202,10 @@ contains ! get the score and tally it score = last_wgt * calc_pn(n, mu) -!$omp atomic +!$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score +!$omp end critical end do j = j + t % scatt_order(j) cycle SCORE_LOOP @@ -341,9 +342,10 @@ contains end select ! Add score to tally -!$omp atomic +!$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score +!$omp end critical end do SCORE_LOOP @@ -415,9 +417,10 @@ contains i_filter = sum((t % matching_bins - 1) * t % stride) + 1 ! Add score to tally -!$omp atomic +!$omp critical t % results(i_score, i_filter) % value = & t % results(i_score, i_filter) % value + score +!$omp end critical end do ! reset outgoing energy bin and score index @@ -704,9 +707,10 @@ contains score_index = (k - 1)*t % n_score_bins + j ! Add score to tally -!$omp atomic +!$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score +!$omp end critical end do SCORE_LOOP @@ -848,9 +852,10 @@ contains score_index = (i_nuclide - 1)*t % n_score_bins + j ! Add score to tally -!$omp atomic +!$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score +!$omp end critical end do SCORE_LOOP @@ -947,9 +952,10 @@ contains score_index = n_nuclides_total*t % n_score_bins + j ! Add score to tally -!$omp atomic +!$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score +!$omp end critical end do MATERIAL_SCORE_LOOP @@ -1256,9 +1262,10 @@ contains score_index = (b - 1)*t % n_score_bins + j ! Add score to tally -!$omp atomic +!$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score +!$omp end critical end do SCORE_LOOP @@ -1489,9 +1496,10 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 -!$omp atomic +!$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt +!$omp end critical end if end do else @@ -1502,9 +1510,10 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 -!$omp atomic +!$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt +!$omp end critical end if end do end if @@ -1519,9 +1528,10 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 -!$omp atomic +!$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt +!$omp end critical end if end do else @@ -1532,9 +1542,10 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 -!$omp atomic +!$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt +!$omp end critical end if end do end if @@ -1549,9 +1560,10 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 -!$omp atomic +!$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt +!$omp end critical end if end do else @@ -1562,9 +1574,10 @@ contains t % matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 -!$omp atomic +!$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt +!$omp end critical end if end do end if @@ -1687,9 +1700,10 @@ contains end if ! Add to surface current tally -!$omp atomic +!$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt +!$omp end critical end if ! Calculate new coordinates diff --git a/src/tracking.F90 b/src/tracking.F90 index 6080b2bae1..0ceb5d7e8f 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -60,8 +60,9 @@ contains n_event = 0 ! Add paricle's starting weight to count for normalizing tallies later -!$omp atomic +!$omp critical total_weight = total_weight + p % wgt +!$omp end critical ! Force calculation of cross-sections by setting last energy to zero micro_xs % last_E = ZERO @@ -101,10 +102,11 @@ contains call score_tracklength_tally(p, distance) ! Score track-length estimate of k-eff -!$omp atomic +!$omp critical global_tallies(K_TRACKLENGTH) % value = & global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * & material_xs % nu_fission +!$omp end critical if (d_collision > d_boundary) then ! ==================================================================== @@ -128,10 +130,11 @@ contains ! PARTICLE HAS COLLISION ! Score collision estimate of keff -!$omp atomic +!$omp critical global_tallies(K_COLLISION) % value = & global_tallies(K_COLLISION) % value + p % wgt * & material_xs % nu_fission / material_xs % total +!$omp end critical ! 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 From be05b1ce55a62614a604ec10628eac7435e450e6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Aug 2013 19:34:51 -0500 Subject: [PATCH 03/16] Removed matching_bins from TallyObject type and made it threadprivate. Still have race condition with Lists. --- src/cmfd_data.F90 | 85 ++++++++++++------------- src/global.F90 | 3 +- src/output.F90 | 76 +++++++++++------------ src/tally.F90 | 130 +++++++++++++++++++-------------------- src/tally_header.F90 | 3 - src/tally_initialize.F90 | 14 +++-- 6 files changed, 158 insertions(+), 153 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 9943379784..a552b9d000 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -57,7 +57,8 @@ contains FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT, & OUT_FRONT, IN_TOP, OUT_TOP, CMFD_NOACCEL, ZERO, ONE use error, only: fatal_error - use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes + use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes,& + matching_bins use mesh, only: mesh_indices_to_bin use mesh_header, only: StructuredMesh use tally_header, only: TallyObject @@ -137,21 +138,21 @@ contains TALLY: if (ital == 1) then ! reset all bins to 1 - t % matching_bins = 1 + matching_bins(1:t%n_filters) = 1 ! set ijk as mesh indices ijk = (/ i, j, k /) ! get bin number for mesh indices - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk) + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk) ! apply energy in filter if (i_filter_ein > 0) then - t % matching_bins(i_filter_ein) = ng - h + 1 + matching_bins(i_filter_ein) = ng - h + 1 end if ! calculate score index from bins - score_index = sum((t % matching_bins - 1) * t%stride) + 1 + score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1 ! get flux flux = t % results(1,score_index) % sum @@ -180,24 +181,24 @@ contains INGROUP: do g = 1, ng ! reset all bins to 1 - t % matching_bins = 1 + matching_bins(1:t%n_filters) = 1 ! set ijk as mesh indices ijk = (/ i, j, k /) ! get bin number for mesh indices - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk) + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk) if (i_filter_ein > 0) then ! apply energy in filter - t % matching_bins(i_filter_ein) = ng - h + 1 + matching_bins(i_filter_ein) = ng - h + 1 ! set energy out bin - t % matching_bins(i_filter_eout) = ng - g + 1 + matching_bins(i_filter_eout) = ng - g + 1 end if ! calculate score index from bins - score_index = sum((t % matching_bins - 1) * t%stride) + 1 + score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1 ! get scattering cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /& @@ -216,69 +217,69 @@ contains else if (ital == 3) then ! initialize and filter for energy - t % matching_bins = 1 + matching_bins(1:t%n_filters) = 1 if (i_filter_ein > 0) then - t % matching_bins(i_filter_ein) = ng - h + 1 + matching_bins(i_filter_ein) = ng - h + 1 end if ! left surface - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i-1, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_RIGHT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing + matching_bins(i_filter_surf) = IN_RIGHT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum - t % matching_bins(i_filter_surf) = OUT_RIGHT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming + matching_bins(i_filter_surf) = OUT_RIGHT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum ! right surface - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_RIGHT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming + matching_bins(i_filter_surf) = IN_RIGHT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum - t % matching_bins(i_filter_surf) = OUT_RIGHT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing + matching_bins(i_filter_surf) = OUT_RIGHT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum ! back surface - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j-1, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_FRONT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing + matching_bins(i_filter_surf) = IN_FRONT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum - t % matching_bins(i_filter_surf) = OUT_FRONT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming + matching_bins(i_filter_surf) = OUT_FRONT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum ! front surface - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_FRONT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming + matching_bins(i_filter_surf) = IN_FRONT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum - t % matching_bins(i_filter_surf) = OUT_FRONT - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing + matching_bins(i_filter_surf) = OUT_FRONT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum ! bottom surface - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k-1 /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_TOP - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing + matching_bins(i_filter_surf) = IN_TOP + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum - t % matching_bins(i_filter_surf) = OUT_TOP - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming + matching_bins(i_filter_surf) = OUT_TOP + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum ! top surface - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_TOP - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming + matching_bins(i_filter_surf) = IN_TOP + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum - t % matching_bins(i_filter_surf) = OUT_TOP - score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing + matching_bins(i_filter_surf) = OUT_TOP + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum end if TALLY diff --git a/src/global.F90 b/src/global.F90 index a84e9b00b8..dd436fefc2 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -97,6 +97,7 @@ module global type(StructuredMesh), allocatable, target :: meshes(:) type(TallyObject), allocatable, target :: tallies(:) + integer, allocatable :: matching_bins(:) ! Pointers for different tallies type(TallyObject), pointer :: user_tallies(:) => null() @@ -384,7 +385,7 @@ module global logical :: output_tallies = .true. !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, & -!$omp& trace, thread_id, current_work) +!$omp& trace, thread_id, current_work, matching_bins) contains diff --git a/src/output.F90 b/src/output.F90 index 4e1d903b8e..472cb5574b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1644,7 +1644,7 @@ contains ! to be used for a given tally. ! Initialize bins, filter level, and indentation - t % matching_bins = 0 + matching_bins(1:t%n_filters) = 0 j = 1 indent = 0 @@ -1654,16 +1654,16 @@ contains if (t % n_filters == 0) exit find_bin ! Increment bin combination - t % matching_bins(j) = t % matching_bins(j) + 1 + matching_bins(j) = matching_bins(j) + 1 ! ================================================================= ! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER - if (t % matching_bins(j) > t % filters(j) % n_bins) then + if (matching_bins(j) > t % filters(j) % n_bins) then ! If this is the first filter, then exit if (j == 1) exit print_bin - t % matching_bins(j) = 0 + matching_bins(j) = 0 j = j - 1 indent = indent - 2 @@ -1696,7 +1696,7 @@ contains ! bins below the lowest filter level will be zeros if (t % n_filters > 0) then - filter_index = sum((max(t % matching_bins,1) - 1) * t % stride) + 1 + filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * t % stride) + 1 else filter_index = 1 end if @@ -1805,7 +1805,7 @@ contains m => meshes(t % filters(i_filter_mesh) % int_bins(1)) ! initialize bins array - t % matching_bins = 1 + matching_bins(1:t%n_filters) = 1 ! determine how many energy in bins there are i_filter_ein = t % find_filter(FILTER_ENERGYIN) @@ -1831,7 +1831,7 @@ contains do l = 1, n if (print_ebin) then ! Set incoming energy bin - t % matching_bins(i_filter_ein) = l + matching_bins(i_filter_ein) = l ! Write incoming energy bin write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') & @@ -1839,102 +1839,102 @@ contains end if ! Left Surface - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_RIGHT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = IN_RIGHT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - t % matching_bins(i_filter_surf) = OUT_RIGHT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = OUT_RIGHT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Right Surface - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_RIGHT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = IN_RIGHT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - t % matching_bins(i_filter_surf) = OUT_RIGHT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = OUT_RIGHT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Back Surface - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_FRONT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = IN_FRONT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - t % matching_bins(i_filter_surf) = OUT_FRONT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = OUT_FRONT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Front Surface - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_FRONT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = IN_FRONT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Front", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - t % matching_bins(i_filter_surf) = OUT_FRONT - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = OUT_FRONT + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Front", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Bottom Surface - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_TOP - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = IN_TOP + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - t % matching_bins(i_filter_surf) = OUT_TOP - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = OUT_TOP + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Top Surface - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) - t % matching_bins(i_filter_surf) = IN_TOP - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = IN_TOP + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Top", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - t % matching_bins(i_filter_surf) = OUT_TOP - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + matching_bins(i_filter_surf) = OUT_TOP + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Top", & to_str(t % results(1,filter_index) % sum), & @@ -1965,7 +1965,7 @@ contains real(8) :: E1 ! upper bound for energy bin type(StructuredMesh), pointer :: m => null() - bin = t % matching_bins(i_filter) + bin = matching_bins(i_filter) select case(t % filters(i_filter) % type) case (FILTER_UNIVERSE) diff --git a/src/tally.F90 b/src/tally.F90 index 2fafacf320..e1fcc2ff06 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -85,7 +85,7 @@ contains ! be accumulating the tally values ! Determine scoring index for this filter combination - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! Check for nuclide bins k = 0 @@ -388,7 +388,7 @@ contains ! save original outgoing energy bin and score index i = t % find_filter(FILTER_ENERGYOUT) - bin_energyout = t % matching_bins(i) + bin_energyout = matching_bins(i) ! Get number of energies on filter n = size(t % filters(i) % real_bins) @@ -411,10 +411,10 @@ contains E_out > t % filters(i) % real_bins(n)) cycle ! change outgoing energy bin - t % matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out) + matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out) ! determine scoring index - i_filter = sum((t % matching_bins - 1) * t % stride) + 1 + i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! Add score to tally !$omp critical @@ -424,7 +424,7 @@ contains end do ! reset outgoing energy bin and score index - t % matching_bins(i) = bin_energyout + matching_bins(i) = bin_energyout end subroutine score_fission_eout @@ -495,7 +495,7 @@ contains ! be accumulating the tally values ! Determine scoring index for this filter combination - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 if (t % all_nuclides) then call score_all_nuclides(p, i_tally, flux, filter_index) @@ -1005,7 +1005,7 @@ contains !$omp threadprivate(t, m, mat, coord) t => tallies(i_tally) - t % matching_bins = 1 + matching_bins(1:t%n_filters) = 1 ! ========================================================================== ! CHECK IF THIS TRACK INTERSECTS THE MESH @@ -1045,11 +1045,11 @@ contains case (FILTER_UNIVERSE) ! determine next universe bin ! TODO: Account for multiple universes when performing this filter - t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & + matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & p % coord % universe, i_tally) case (FILTER_MATERIAL) - t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & + matching_bins(i) = get_next_bin(FILTER_MATERIAL, & p % material, i_tally) case (FILTER_CELL) @@ -1057,21 +1057,21 @@ contains coord => p % coord0 do while(associated(coord)) position(FILTER_CELL) = 0 - t % matching_bins(i) = get_next_bin(FILTER_CELL, & + matching_bins(i) = get_next_bin(FILTER_CELL, & coord % cell, i_tally) - if (t % matching_bins(i) /= NO_BIN_FOUND) exit + if (matching_bins(i) /= NO_BIN_FOUND) exit coord => coord % next end do nullify(coord) case (FILTER_CELLBORN) ! determine next cellborn bin - t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & + matching_bins(i) = get_next_bin(FILTER_CELLBORN, & p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin - t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & + matching_bins(i) = get_next_bin(FILTER_SURFACE, & p % surface, i_tally) case (FILTER_ENERGYIN) @@ -1081,17 +1081,17 @@ contains ! check if energy of the particle is within energy bins if (p % E < t % filters(i) % real_bins(1) .or. & p % E > t % filters(i) % real_bins(k + 1)) then - t % matching_bins(i) = NO_BIN_FOUND + matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin - t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & + matching_bins(i) = binary_search(t % filters(i) % real_bins, & k + 1, p % E) end if end select ! Check if no matching bin was found - if (t % matching_bins(i) == NO_BIN_FOUND) return + if (matching_bins(i) == NO_BIN_FOUND) return end do FILTER_LOOP @@ -1161,10 +1161,10 @@ contains flux = p % wgt * distance ! Determine mesh bin - t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross) + matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross) ! Determining scoring index - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 if (t % all_nuclides) then ! Score reaction rates for each nuclide in material @@ -1301,7 +1301,7 @@ contains found_bin = .true. t => tallies(i_tally) - t % matching_bins = 1 + matching_bins(1:t%n_filters) = 1 FILTER_LOOP: do i = 1, t % n_filters @@ -1311,16 +1311,16 @@ 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, t % matching_bins(i)) + call get_mesh_bin(m, p % coord0 % xyz, matching_bins(i)) case (FILTER_UNIVERSE) ! determine next universe bin ! TODO: Account for multiple universes when performing this filter - t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & + matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & p % coord % universe, i_tally) case (FILTER_MATERIAL) - t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & + matching_bins(i) = get_next_bin(FILTER_MATERIAL, & p % material, i_tally) case (FILTER_CELL) @@ -1328,21 +1328,21 @@ contains coord => p % coord0 do while(associated(coord)) position(FILTER_CELL) = 0 - t % matching_bins(i) = get_next_bin(FILTER_CELL, & + matching_bins(i) = get_next_bin(FILTER_CELL, & coord % cell, i_tally) - if (t % matching_bins(i) /= NO_BIN_FOUND) exit + if (matching_bins(i) /= NO_BIN_FOUND) exit coord => coord % next end do nullify(coord) case (FILTER_CELLBORN) ! determine next cellborn bin - t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & + matching_bins(i) = get_next_bin(FILTER_CELLBORN, & p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin - t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & + matching_bins(i) = get_next_bin(FILTER_SURFACE, & p % surface, i_tally) case (FILTER_ENERGYIN) @@ -1359,10 +1359,10 @@ contains ! check if energy of the particle is within energy bins if (E < t % filters(i) % real_bins(1) .or. & E > t % filters(i) % real_bins(n + 1)) then - t % matching_bins(i) = NO_BIN_FOUND + matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin - t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & + matching_bins(i) = binary_search(t % filters(i) % real_bins, & n + 1, E) end if @@ -1373,17 +1373,17 @@ contains ! check if energy of the particle is within energy bins if (p % E < t % filters(i) % real_bins(1) .or. & p % E > t % filters(i) % real_bins(n + 1)) then - t % matching_bins(i) = NO_BIN_FOUND + matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin - t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & + matching_bins(i) = binary_search(t % filters(i) % real_bins, & n + 1, p % E) end if end select ! If the current filter didn't match, exit this subroutine - if (t % matching_bins(i) == NO_BIN_FOUND) then + if (matching_bins(i) == NO_BIN_FOUND) then found_bin = .false. return end if @@ -1475,7 +1475,7 @@ contains end if ! search to find incoming energy bin - t % matching_bins(j) = binary_search(t % filters(j) % real_bins, & + matching_bins(j) = binary_search(t % filters(j) % real_bins, & n + 1, p % E) end if @@ -1492,10 +1492,10 @@ contains do j = ijk0(3), ijk1(3) - 1 ijk0(3) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = OUT_TOP - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt @@ -1506,10 +1506,10 @@ contains do j = ijk0(3) - 1, ijk1(3), -1 ijk0(3) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = IN_TOP - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = IN_TOP + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt @@ -1524,10 +1524,10 @@ contains do j = ijk0(2), ijk1(2) - 1 ijk0(2) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = OUT_FRONT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt @@ -1538,10 +1538,10 @@ contains do j = ijk0(2) - 1, ijk1(2), -1 ijk0(2) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = IN_FRONT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = IN_FRONT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt @@ -1556,10 +1556,10 @@ contains do j = ijk0(1), ijk1(1) - 1 ijk0(1) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = OUT_RIGHT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt @@ -1570,10 +1570,10 @@ contains do j = ijk0(1) - 1, ijk1(1), -1 ijk0(1) = j if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = IN_RIGHT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = IN_RIGHT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp critical t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt @@ -1598,7 +1598,7 @@ contains do k = 1, n_cross ! Reset scoring bin index - t % matching_bins(i_filter_surf) = 0 + matching_bins(i_filter_surf) = 0 ! Calculate distance to each bounding surface. We need to treat ! special case where the cosine of the angle is zero since this would @@ -1625,8 +1625,8 @@ contains ! Crossing into right mesh cell -- this is treated as outgoing ! current from (i,j,k) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = OUT_RIGHT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(1) = ijk0(1) + 1 @@ -1637,8 +1637,8 @@ contains ijk0(1) = ijk0(1) - 1 xyz_cross(1) = xyz_cross(1) - m % width(1) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = IN_RIGHT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = IN_RIGHT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if @@ -1647,8 +1647,8 @@ contains ! Crossing into front mesh cell -- this is treated as outgoing ! current in (i,j,k) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = OUT_FRONT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(2) = ijk0(2) + 1 @@ -1659,8 +1659,8 @@ contains ijk0(2) = ijk0(2) - 1 xyz_cross(2) = xyz_cross(2) - m % width(2) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = IN_FRONT - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = IN_FRONT + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if @@ -1669,8 +1669,8 @@ contains ! Crossing into top mesh cell -- this is treated as outgoing ! current in (i,j,k) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = OUT_TOP - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(3) = ijk0(3) + 1 @@ -1681,16 +1681,16 @@ contains ijk0(3) = ijk0(3) - 1 xyz_cross(3) = xyz_cross(3) - m % width(3) if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - t % matching_bins(i_filter_surf) = IN_TOP - t % matching_bins(i_filter_mesh) = & + matching_bins(i_filter_surf) = IN_TOP + matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if end if ! Determine scoring index - if (t % matching_bins(i_filter_surf) > 0) then - filter_index = sum((t % matching_bins - 1) * t % stride) + 1 + if (matching_bins(i_filter_surf) > 0) then + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! Check for errors if (filter_index <= 0 .or. filter_index > & diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 41c5d79715..17cb1cb0de 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -86,7 +86,6 @@ module tally_header ! mapped onto one dimension in the results array, the stride attribute gives ! the stride for a given filter type within the results array - integer, allocatable :: matching_bins(:) integer, allocatable :: stride(:) ! This array provides a way to lookup what index in the filters array a @@ -169,8 +168,6 @@ module tally_header deallocate(this % filters) end if - if (allocated(this % matching_bins)) & - deallocate(this % matching_bins) if (allocated(this % stride)) & deallocate(this % stride) diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 04928499e7..23ab5b5634 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -31,9 +31,10 @@ contains subroutine setup_tally_arrays() - integer :: i ! loop index for tallies - integer :: j ! loop index for filters - integer :: n ! temporary stride + integer :: i ! loop index for tallies + integer :: j ! loop index for filters + integer :: n ! temporary stride + integer :: max_n_filters = 0 ! maximum number of filters type(TallyObject), pointer :: t => null() TALLY_LOOP: do i = 1, n_tallies @@ -42,7 +43,7 @@ contains ! Allocate stride and matching_bins arrays allocate(t % stride(t % n_filters)) - allocate(t % matching_bins(t % n_filters)) + max_n_filters = max(max_n_filters, t % n_filters) ! The filters are traversed in opposite order so that the last filter has ! the shortest stride in memory and the first filter has the largest @@ -63,6 +64,11 @@ contains end do TALLY_LOOP + ! Allocate array for matching filter bins +!$omp parallel + allocate(matching_bins(max_n_filters)) +!$omp end parallel + end subroutine setup_tally_arrays !=============================================================================== From 68d5ad00a2ed6617ba97e1895b52c01fb2a7c758 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Aug 2013 10:50:25 -0400 Subject: [PATCH 04/16] Make active tally lists threadprivate. --- src/eigenvalue.F90 | 2 ++ src/global.F90 | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index f52ccef810..e46d5c996b 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -127,7 +127,9 @@ contains tallies_on = .true. ! Add user tallies to active tallies list +!$omp parallel call setup_active_usertallies() +!$omp end parallel end if ! check CMFD initialize batch diff --git a/src/global.F90 b/src/global.F90 index dd436fefc2..9466589b1b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -112,6 +112,8 @@ module global type(SetInt) :: active_tracklength_tallies type(SetInt) :: active_current_tallies type(SetInt) :: active_tallies +!$omp threadprivate(active_analog_tallies, active_tracklength_tallies, & +!$omp& active_current_tallies, active_tallies) ! Global tallies ! 1) collision estimate of k-eff From 8f0f884eec7ea422fccf25a869fabdb42179dd15 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Aug 2013 11:13:47 -0400 Subject: [PATCH 05/16] Added OpenMP threading for fixed source runs. --- src/fixed_source.F90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/fixed_source.F90 b/src/fixed_source.F90 index b057008728..cac08a1559 100644 --- a/src/fixed_source.F90 +++ b/src/fixed_source.F90 @@ -12,6 +12,7 @@ module fixed_source use tracking, only: transport type(Bank), pointer :: source_site => null() +!$omp threadprivate(source_site) contains @@ -23,11 +24,15 @@ contains if (master) call header("FIXED SOURCE TRANSPORT SIMULATION", level=1) ! Allocate particle and dummy source site +!$omp parallel allocate(source_site) +!$omp end parallel ! Turn timer and tallies on tallies_on = .true. +!$omp parallel call setup_active_usertallies() +!$omp end parallel call time_active % start() ! ========================================================================== @@ -47,6 +52,7 @@ contains ! ======================================================================= ! LOOP OVER PARTICLES +!$omp parallel do schedule(static) firstprivate(p) PARTICLE_LOOP: do i = 1, work ! Set unique particle ID @@ -67,6 +73,7 @@ contains call transport(p) end do PARTICLE_LOOP +!$omp end parallel do ! Accumulate time for transport call time_transport % stop() From a0ddce284e6fe8c6cbcab56c7114b5242156fdda Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Aug 2013 12:50:34 -0400 Subject: [PATCH 06/16] Enclose global leakage tally in omp critical. --- src/geometry.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 547ad21d36..7e1bc283c1 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -408,8 +408,12 @@ contains end if ! Score to global leakage tally - if (tallies_on) global_tallies(LEAKAGE) % value = & + if (tallies_on) then +!$omp critical + global_tallies(LEAKAGE) % value = & global_tallies(LEAKAGE) % value + p % wgt +!$omp end critical + end if ! Display message if (verbosity >= 10 .or. trace) then From ad20baa03f2a4b2dbe5aed57cfa69054076fe8b3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Aug 2013 16:02:48 -0400 Subject: [PATCH 07/16] Add OPENMP Makefile variable to installation instructions. --- docs/source/usersguide/install.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index ef01aeff7f..ee7697de3e 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -135,6 +135,10 @@ MPI Enables parallel runs using the Message Passing Interface. The MPI_DIR variable should be set to the base directory of the MPI implementation. +OPENMP + Enables shared-memory parallelism using the OpenMP API. The Fortran compiler + being used must support OpenMP. + HDF5 Enables HDF5 output in addition to normal screen and text file output. The HDF5_DIR variable should be set to the base directory of the HDF5 From dc096bda0f0aaf4d36aae50c7dcc0152b1c9b8cf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 23 Aug 2013 12:31:37 -0400 Subject: [PATCH 08/16] Make sure thread fission banks are deallocated. --- src/global.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/global.F90 b/src/global.F90 index b8a0b4bcfd..1530f13a09 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -446,13 +446,19 @@ contains ! Now deallocate the tally array deallocate(tallies) end if + if (allocated(matching_bins)) deallocate(matching_bins) if (allocated(tally_maps)) deallocate(tally_maps) ! Deallocate energy grid if (allocated(e_grid)) deallocate(e_grid) ! Deallocate fission and source bank and entropy +!$omp parallel if (allocated(fission_bank)) deallocate(fission_bank) +!$omp end parallel +#ifdef OPENMP + if (allocated(master_fission_bank)) deallocate(master_fission_bank) +#endif if (allocated(source_bank)) deallocate(source_bank) if (allocated(entropy_p)) deallocate(entropy_p) From b7e4fd6d8416184ff398974983ce96ff6320e221 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 28 Aug 2013 01:05:43 -0400 Subject: [PATCH 09/16] Change last instance of omp atomic to omp critical. --- src/geometry.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 7e1bc283c1..0e6700d851 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -1572,8 +1572,9 @@ contains ! Increment number of lost particles p % alive = .false. -!$omp atomic +!$omp critical n_lost_particles = n_lost_particles + 1 +!$omp end critical ! Abort the simulation if the maximum number of lost particles has been ! reached From f6766c169d1a7e133bd244ee6ac8ea44680b8d54 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 1 Sep 2013 23:45:38 -0400 Subject: [PATCH 10/16] Add option to set number of threads in settings.xml and command line. --- src/global.F90 | 2 +- src/initialize.F90 | 17 +++++++++++++++++ src/input_xml.F90 | 18 ++++++++++++++++++ src/output.F90 | 1 + src/templates/settings.rnc | 2 ++ src/templates/settings_t.xml | 1 + 6 files changed, 40 insertions(+), 1 deletion(-) diff --git a/src/global.F90 b/src/global.F90 index 1530f13a09..f96f59cf08 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -212,7 +212,7 @@ module global integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult #ifdef OPENMP - integer :: n_threads ! number of OpenMP threads + integer :: n_threads = NONE ! number of OpenMP threads integer :: thread_id ! ID of a given thread #endif diff --git a/src/initialize.F90 b/src/initialize.F90 index d82b1a7276..da2ace4d09 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -367,6 +367,23 @@ contains case ('-g', '-geometry-debug', '--geometry-debug') check_overlaps = .true. + case ('-t', '--threads') + ! Read number of threads + i = i + 1 + +#ifdef OPENMP + ! Read and set number of OpenMP threads + n_threads = str_to_int(argv(i)) + if (n_threads < 1) then + message = "Invalid number of threads specified on command line." + call fatal_error() + end if + call omp_set_num_threads(n_threads) +#else + message = "Ignoring number of threads specified on command line." + call warning() +#endif + case ('-?', '-help', '--help') call print_usage() stop diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 617788e472..c2d1a86d16 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -73,6 +73,7 @@ contains cross_sections_ = '' output_path_ = '' verbosity_ = 0 + threads_ = NONE energy_grid_ = 'union' seed_ = 0_8 source_ % file = '' @@ -207,6 +208,23 @@ contains ! Verbosity if (verbosity_ > 0) verbosity = verbosity_ + ! Number of OpenMP threads + if (threads_ /= NONE) then +#ifdef OPENMP + if (n_threads == NONE) then + n_threads = threads_ + if (n_threads < 1) then + message = "Invalid number of threads: " // to_str(n_threads) + call fatal_error() + end if + call omp_set_num_threads(n_threads) + end if +#else + message = "Ignoring number of threads." + call warning() +#endif + end if + ! ========================================================================== ! EXTERNAL SOURCE diff --git a/src/output.F90 b/src/output.F90 index 472cb5574b..ffffb745c5 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -170,6 +170,7 @@ contains write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode' write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point' write(OUTPUT_UNIT,*) ' or a particle restart file' + write(OUTPUT_UNIT,*) ' -t, --threads Number of OpenMP threads' write(OUTPUT_UNIT,*) ' -v, --version Show version information' write(OUTPUT_UNIT,*) ' -?, --help Show this message' end if diff --git a/src/templates/settings.rnc b/src/templates/settings.rnc index 47dd690bcc..cd4ed4beb9 100644 --- a/src/templates/settings.rnc +++ b/src/templates/settings.rnc @@ -97,6 +97,8 @@ element settings { element survival_biasing { xsd:boolean }? & + element threads { xsd:positiveInteger }? & + element trace { list { xsd:positiveInteger+ } }? & element verbosity { xsd:positiveInteger }? & diff --git a/src/templates/settings_t.xml b/src/templates/settings_t.xml index 3afa61745e..8f1cc7b00c 100644 --- a/src/templates/settings_t.xml +++ b/src/templates/settings_t.xml @@ -58,6 +58,7 @@ + From 7e5e146f4ebd0cb5c58b05f2d895ef83633725a1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Sep 2013 19:10:52 -0400 Subject: [PATCH 11/16] Change command-line options. --- src/initialize.F90 | 4 ++-- src/output.F90 | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 8307c1e3d1..dea901fb66 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -335,7 +335,7 @@ contains run_mode = MODE_PLOTTING check_overlaps = .true. - case ('-n', '-n_particles', '--n_particles') + case ('-n', '-particles', '--particles') ! Read number of particles per cycle i = i + 1 n_particles = str_to_int(argv(i)) @@ -368,7 +368,7 @@ contains case ('-g', '-geometry-debug', '--geometry-debug') check_overlaps = .true. - case ('-t', '--threads') + case ('-s', '--threads') ! Read number of threads i = i + 1 diff --git a/src/output.F90 b/src/output.F90 index ffffb745c5..8b54dc54a9 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -167,10 +167,11 @@ contains write(OUTPUT_UNIT,*) write(OUTPUT_UNIT,*) 'Options:' write(OUTPUT_UNIT,*) ' -g, --geometry-debug Run in geometry debugging mode' + write(OUTPUT_UNIT,*) ' -n, --particles Number of particles per generation' write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode' write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point' write(OUTPUT_UNIT,*) ' or a particle restart file' - write(OUTPUT_UNIT,*) ' -t, --threads Number of OpenMP threads' + write(OUTPUT_UNIT,*) ' -s, --threads Number of OpenMP threads' write(OUTPUT_UNIT,*) ' -v, --version Show version information' write(OUTPUT_UNIT,*) ' -?, --help Show this message' end if From d44dfd20d9c77be545e3560fa3f05e949d860db6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Sep 2013 20:29:18 -0400 Subject: [PATCH 12/16] Add OpenMP tests (also fix PEP8 compliance in test_compile.py). --- tests/run_tests.py | 13 +++--- tests/test_compile/test_compile.py | 75 +++++++++++++++++++++++------- 2 files changed, 66 insertions(+), 22 deletions(-) diff --git a/tests/run_tests.py b/tests/run_tests.py index 5e261fba73..1a1a05e839 100755 --- a/tests/run_tests.py +++ b/tests/run_tests.py @@ -55,7 +55,7 @@ def run_suite(name=None, mpi=False): results.append((name, result)) # set mpiexec path -if os.environ.has_key('COMPILER'): +if 'COMPILER' in os.environ: compiler = os.environ['COMPILER'] else: compiler = 'gnu' @@ -67,8 +67,9 @@ sys.path.append(pwd) # Set list of tests, either default or from command line flags = [] -tests = ['compile', 'normal', 'debug', 'optimize', 'hdf5', 'mpi', 'phdf5', - 'petsc', 'phdf5-petsc', 'phdf5-petsc-optimize'] +tests = ['compile', 'normal', 'debug', 'optimize', 'mpi', 'omp', 'hdf5', + 'petsc', 'mpi-omp', 'omp-hdf5', 'phdf5', 'phdf5-petsc', + 'phdf5-petsc-optimize', 'omp-phdf5-petsc-optimize'] if len(sys.argv) > 1: flags = [i for i in sys.argv[1:] if i.startswith('-')] tests_ = [i for i in sys.argv[1:] if not i.startswith('-')] @@ -79,10 +80,10 @@ results = [] for name in tests: if name == 'compile': run_compile() - elif name in ['normal', 'debug', 'optimize', 'hdf5']: + elif name in ['normal', 'debug', 'optimize', 'omp', 'hdf5', 'omp-hdf5']: run_suite(name=name) - elif name in ['mpi', 'phdf5', 'petsc', 'phdf5-petsc', - 'phdf5-petsc-optimize']: + elif name in ['mpi', 'mpi-omp' 'phdf5', 'petsc', 'phdf5-petsc', + 'phdf5-petsc-optimize', 'omp-phdf5-petsc-optimize']: run_suite(name=name, mpi=True) # print out summary of results diff --git a/tests/test_compile/test_compile.py b/tests/test_compile/test_compile.py index 72e4dd6535..66368c08ad 100644 --- a/tests/test_compile/test_compile.py +++ b/tests/test_compile/test_compile.py @@ -3,8 +3,8 @@ """Compilation tests This set of tests makes sure that OpenMC can compile for many different -combinations of compilation flags and options. By default, only the gfortran and -ifort compilers are tested. This requires that the MPI, HDF5, and PETSC +combinations of compilation flags and options. By default, only the gfortran +and ifort compilers are tested. This requires that the MPI, HDF5, and PETSC libraries are already set up appropriately in the Makefile. """ @@ -15,80 +15,123 @@ import shutil pwd = os.path.dirname(__file__) -if os.environ.has_key('COMPILER'): +if 'COMPILER' in os.environ: compiler = 'COMPILER=' + os.environ['COMPILER'] else: compiler = 'COMPILER=gnu' + def setup(): # Change to source directory os.chdir(pwd + '/../../src') + def test_normal(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler]) assert returncode == 0 shutil.move('openmc', 'openmc-normal') + def test_debug(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'DEBUG=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-debug') + def test_profile(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'PROFILE=yes']) assert returncode == 0 + def test_optimize(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'OPTIMIZE=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-optimize') + def test_mpi(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'MPI=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-mpi') + +def test_omp(): + returncode = run(['make', 'distclean']) + returncode = run(['make', compiler, 'OPENMP=yes']) + assert returncode == 0 + shutil.move('openmc', 'openmc-omp') + + def test_hdf5(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'HDF5=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-hdf5') + def test_petsc(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'MPI=yes', 'PETSC=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-petsc') + +def test_mpi_omp(): + returncode = run(['make', 'distclean']) + returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes']) + assert returncode == 0 + shutil.move('openmc', 'openmc-mpi-omp') + + +def test_omp_hdf5(): + returncode = run(['make', 'distclean']) + returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'HDF5=yes']) + assert returncode == 0 + shutil.move('openmc', 'openmc-omp-hdf5') + + def test_mpi_hdf5(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-phdf5') + def test_mpi_hdf5_petsc(): - returncode = run(['make','distclean']) + returncode = run(['make', 'distclean']) returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-phdf5-petsc') + def test_mpi_hdf5_petsc_optimize(): - returncode = run(['make','distclean']) - returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes', 'OPTIMIZE=yes']) + returncode = run(['make', 'distclean']) + returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes', + 'OPTIMIZE=yes']) assert returncode == 0 shutil.move('openmc', 'openmc-phdf5-petsc-optimize') + +def test_mpi_omp_hdf5_petsc_optimize(): + returncode = run(['make', 'distclean']) + returncode = run(['make', compiler, 'MPI=yes', 'OMP=yes', 'HDF5=yes', + 'PETSC=yes', 'OPTIMIZE=yes']) + assert returncode == 0 + shutil.move('openmc', 'openmc-omp-phdf5-petsc-optimize') + + def run(commands): proc = Popen(commands, stderr=STDOUT, stdout=PIPE) returncode = proc.wait() print(proc.communicate()[0]) return returncode + def teardown(commands): - returncode = run(['make','distclean']) - shutil.copy('openmc-normal','openmc') + returncode = run(['make', 'distclean']) + shutil.copy('openmc-normal', 'openmc') From 71bfed38fb12fe9f5b8fab7ced55a20bf0441666 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Sep 2013 22:34:22 -0400 Subject: [PATCH 13/16] Add description of element in user's guide. --- docs/source/usersguide/input.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index d9b43c8214..3ce8934108 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -381,6 +381,14 @@ survival biasing, otherwise known as implicit capture or absorption. .. _trace: +```` Element +--------------------- + +The ```` element indicates the number of OpenMP threads to be used for +a simulation. It has no attributes and accept a positive integer value. + + *Default*: None (Determined by environment variable :envvar:`OMP_NUM_THREADS`) + ```` Element ------------------- From 45fb774d6d6ff1dec3a7683fb86814e862400b52 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Sep 2013 22:56:46 -0400 Subject: [PATCH 14/16] Mention command-line flags in user's guide. --- docs/source/usersguide/install.rst | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index d855aca499..159c689e98 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -330,6 +330,20 @@ Alternatively, you could run from any directory: Note that in the latter case, any output files will be placed in the present working directory which may be different from ``/home/username/somemodel``. +Command-Line Flags +------------------ + +OpenMC accepts the following command line flags: + +-g, --geometry-debug Run in geometry debugging mode, where cell overlaps are + checked for after each move of a particle +-n, --particles N Use *N* particles per generation or batch +-p, --plot Run in plotting mode +-r, --restart file Restart a previous run from a state point or a particle + restart file +-s, --threads N Run with *N* OpenMP threads +-v, --version Show version information + ----------------------------------------------------- Configuring Input Validation with GNU Emacs nXML mode ----------------------------------------------------- From 68f636294e98402d41bdcbcea0a4850867f67147 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Sep 2013 22:58:28 -0400 Subject: [PATCH 15/16] Update man page with -s, --threads option. --- man/man1/openmc.1 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index 1614e48751..fb0fdd7177 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -24,6 +24,9 @@ Run in plotting mode Restart a previous run from a state point or a particle restart file named \fIbinaryFile\fP. .TP +.BI \-s " N" "\fR,\fP \-\-threads" " N" +Use \fIN\fP OpenMP threads. +.TP .B "\-v\fR, \fP\-\-version" Show version information. .TP From 8aa03c1cbd0fdb6fbbc5230035175e269189cd11 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 11 Sep 2013 10:35:07 -0400 Subject: [PATCH 16/16] Fixed a typo in addition to docs entry. --- docs/source/usersguide/input.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 3ce8934108..425df5029a 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -385,7 +385,7 @@ survival biasing, otherwise known as implicit capture or absorption. --------------------- The ```` element indicates the number of OpenMP threads to be used for -a simulation. It has no attributes and accept a positive integer value. +a simulation. It has no attributes and accepts a positive integer value. *Default*: None (Determined by environment variable :envvar:`OMP_NUM_THREADS`)