diff --git a/docs/source/_images/uniongrid.svg b/docs/source/_images/uniongrid.svg deleted file mode 100644 index 27c3922fd..000000000 --- a/docs/source/_images/uniongrid.svg +++ /dev/null @@ -1,792 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - image/svg+xml - - - - - - - - E1 - - E2 - - E3 - - E4 - - E5 - - E6 - - E7 - - E8 - UnionEnergy Grid - - 0 - NuclidePointers - - 0 - - 1 - - 1 - - 1 - - 2 - - 3 - - 3 - - - - - - - - - - E1 - - E2 - - E3 - NuclideEnergy Grid - - σ1 - - σ2 - - σ3 - NuclideCross Sections - - - - diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index a193dde96..db2e5156e 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -47,45 +47,28 @@ there would be for burnup calculations. Thus, there is a strong motive to implement a method of reducing the number of energy grid searches in order to speed up the calculation. -Unionized Energy Grid ---------------------- +Logarithmic Mapping +------------------- -The most naïve method to reduce the number of energy grid searches is to -construct a new energy grid that consists of the union of the energy points of -each nuclide and use this energy grid for all nuclides. This method is -computationally very efficient as it only requires one energy grid search at -each collision as well as one interpolation between cross section values since -the interpolation factor can be used for all nuclides. However, it requires -redundant storage of cross section values at points which were added to each -nuclide grid. This additional burden on memory storage can become quite -prohibitive. To lessen that burden, the unionized energy grid can be thinned -with cross sections reconstructed on the thinned energy grid. This method is -currently used by default in the Serpent Monte Carlo code. +To speed up energy grid searches, OpenMC uses logarithmic mapping technique +[Brown]_ to limit the range of energies that must be searched for each +nuclide. The entire energy range is divided up into equal-lethargy segments, and +the bounding energies of each segment are mapped to bounding indices on each of +the nuclide energy grids. By default, OpenMC uses 8000 equal-lethargy segments +as recommended by Brown. -Unionized Energy Grid with Nuclide Pointers -------------------------------------------- +Other Methods +------------- -While having a unionized grid that is used for all nuclides allows for very fast -lookup of cross sections, the burden on memory is in many circumstances -unacceptable. The OpenMC Monte Carlo code utilizes a method that allows for a -single energy grid search to be performed at every collision while avoiding the -redundant storage of cross section values. Instead of using the unionized grid -for every nuclide, the original energy grid of each nuclide is kept and a list -of pointers (of the same length as the unionized energy grid) is constructed for -each nuclide that gives the corresponding grid index on the nuclide grid for a -given grid index on the unionized grid. One must still interpolate on cross -section values for each nuclide since the interpolation factors will generally -be different. The figure below illustrates this method. All values within the -dashed box would need to be stored on a per-nuclide basis, and the union grid -would need to be stored once. This method is also referred to as *double -indexing* and is available as an option in Serpent (see paper by Leppanen_). +A good survey of other energy grid techniques, including unionized energy grids, +can be found in a paper by Leppanen_. -.. figure:: ../_images/uniongrid.* - :width: 600px - :align: center - :figclass: align-center +---------- +References +---------- - Mapping of union energy grid to nuclide energy grid through pointers. +.. [Brown] Forrest B. Brown, "New Hash-based Energy Lookup Algorithm for Monte + Carlo codes," LA-UR-14-24530, Los Alamos National Laboratory (2014). .. _MCNP: http://mcnp.lanl.gov .. _Serpent: http://montecarlo.vtt.fi diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index 3439687b3..8de7e9bd5 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -516,7 +516,7 @@ satisfy the following equations x^2 + y^2 + z^2 - 10^2 < 0 \\ x - (-3) > 0 \\ - x - 2 < 0 + y - 2 < 0 In order to determine if a point is inside the cell, we would substitute its coordinates into equation :eq:`cell-contains-example`. If the inequalities are diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index 54dc91345..db6cfd89e 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -682,17 +682,20 @@ nuclear temperature, which is a function of the incoming energy of the neutron. The ACE format contains a list of nuclear temperatures versus incoming energies. The nuclear temperature is interpolated between neighboring incoming energies using a specified interpolation law. Once the temperature :math:`T` is -determined, we then calculate a candidate outgoing energy based on rule C45 in -the `Monte Carlo Sampler`_: +determined, we then calculate a candidate outgoing energy based on the algorithm +given in LA-UR-14-27694_: .. math:: :label: evaporation-E - E' = -T \log (\xi_1 \xi_2) + E' = -T \log ((1 - g\xi_1)(1 - g\xi_2)) -where :math:`\xi_1, \xi_2` are random numbers sampled on the unit -interval. The outgoing energy is only accepted according to a specified -restriction energy as in equation :eq:`maxwell-restriction`. +where :math:`g = 1 - e^{-w}`, :math:`w = (E - U)/T`, :math:`U` is the +restriction energy, and :math:`\xi_1, \xi_2` are random numbers sampled on the +unit interval. The outgoing energy is only accepted according to the restriction +energy as in equation :eq:`maxwell-restriction`. This algorithm has a much +higher rejection efficiency than the standard technique, i.e. rule C45 in the +`Monte Carlo Sampler`_. ACE Law 11 - Energy-Dependent Watt Spectrum +++++++++++++++++++++++++++++++++++++++++++ @@ -1591,6 +1594,8 @@ References .. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721_3rdmcsampler.pdf +.. _LA-UR-14-27694: http://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-14-27694 + .. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf .. _Sutton and Brown: http://www.osti.gov/bridge/product.biblio.jsp?osti_id=307911 diff --git a/docs/source/quickinstall.rst b/docs/source/quickinstall.rst index 9dd6be6d1..3fb9036d3 100644 --- a/docs/source/quickinstall.rst +++ b/docs/source/quickinstall.rst @@ -35,9 +35,9 @@ OpenMC from source as is described in :ref:`usersguide_install`. Installing from Source on Linux or Mac OS X ------------------------------------------- -All OpenMC source code is hosted on GitHub_. If you have git_ and the gfortran_ -compiler installed, you can download and install OpenMC be entering the -following commands in a terminal: +All OpenMC source code is hosted on GitHub_. If you have git_, the gfortran_ +compiler, and CMake_ installed, you can download and install OpenMC be entering +the following commands in a terminal: .. code-block:: sh @@ -58,3 +58,4 @@ can be replaced with a local install, e.g. .. _GitHub: https://github.com/mit-crpg/openmc .. _git: http://git-scm.com .. _gfortran: http://gcc.gnu.org/wiki/GFortran +.. _CMake: http://www.cmake.org diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 231cfac75..393004dc4 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -134,13 +134,15 @@ should be performed. It has the following attributes/sub-elements: ------------------------- The ```` element determines the treatment of the energy grid during -a simulation. Setting this element to "nuclide" will cause OpenMC to use a -nuclide's energy grid when determining what points to interpolate between for -determining cross sections (i.e. non-unionized energy grid). To use a unionized -energy grid, set this element to "union". Note that the unionized energy grid -treatment is slightly different than that employed in Serpent. +a simulation. The valid options are "nuclide" and "logarithm". Setting this +element to "nuclide" will cause OpenMC to use a nuclide's energy grid when +determining what points to interpolate between for determining cross sections +(i.e. non-unionized energy grid). Setting this element to "logarithm" causes +OpenMC to use a logarithmic mapping technique described in LA-UR-14-24530_. - *Default*: union + *Default*: logarithm + +.. _LA-UR-14-24530: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf ```` Element --------------------- @@ -182,6 +184,16 @@ performed. It has the following attributes/sub-elements: *Default*: None +```` Element +--------------------------- + +The ```` element indicates the number of bins to use for the +logarithmic-mapped energy grid. Using more bins will result in energy grid +searches over a smaller range at the expense of more memory. The default is +based on the recommended value in LA-UR-14-24530_. + + *Default*: 8000 + .. _natural_elements: ```` Element @@ -1097,6 +1109,16 @@ The ```` element accepts the following sub-elements: *Default*: total + :estimator: + The estimator element is used to force the use of either ``analog`` or + ``tracklength`` tally estimation. ''analog'' is generally less efficient + though it can be used with every score type. ''tracklength'' is generally + the most efficient, though its usage is restricted to tallies that do not + score particle information which requires a collision to have occured, such + as a scattering tally which utilizes outgoing energy filters. + + *Default*: ``tracklength`` but will revert to analog if necessary. + :scores: A space-separated list of the desired responses to be accumulated. Accepted options are "flux", "total", "scatter", "absorption", "fission", @@ -1647,3 +1669,15 @@ into MATLAB using PETSc-MATLAB utilities. This option can be turned on with "true" and off with "false". *Default*: false + +------------------------------------ +ERSN-OpenMC Graphical User Interface +------------------------------------ + +A third-party Java-based user-friendly graphical user interface for creating XML +input files called ERSN-OpenMC_ is developed and maintained by members of the +Radiation and Nuclear Systems Group at the Faculty of Sciences Tetouan, Morocco. +The GUI also allows one to automatically download prerequisites for installing and +running OpenMC. + +.. _ERSN-OpenMC: https://github.com/EL-Bakkali-Jaafar/ERSN-OpenMC diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 76492dfe8..afed796e2 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -101,7 +101,7 @@ module ace_header ! Energy grid information integer :: n_grid ! # of nuclide grid points - integer, allocatable :: grid_index(:) ! pointers to union grid + integer, allocatable :: grid_index(:) ! union grid pointers / log grid mapping real(8), allocatable :: energy(:) ! energy values corresponding to xs ! Microscopic cross sections diff --git a/src/constants.F90 b/src/constants.F90 index af12b97be..7e30bd07e 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -365,9 +365,8 @@ module constants ! Energy grid methods integer, parameter :: & - GRID_NUCLIDE = 1, & ! non-unionized energy grid - GRID_UNION = 2, & ! union grid with pointers - GRID_LETHARGY = 3 ! lethargy mapping + GRID_NUCLIDE = 1, & ! non-unionized energy grid + GRID_LOGARITHM = 2 ! logarithmic mapping ! Running modes integer, parameter :: & diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 8e4c19ccc..1672405e8 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -2,6 +2,7 @@ module cross_section use ace_header, only: Nuclide, SAlphaBeta, Reaction, UrrData use constants + use energy_grid, only: grid_method, log_spacing use error, only: fatal_error use fission, only: nu_total use global @@ -14,9 +15,6 @@ module cross_section implicit none save - integer :: union_grid_index -!$omp threadprivate(union_grid_index) - contains !=============================================================================== @@ -50,9 +48,6 @@ contains mat => materials(p % material) - ! Find energy index on unionized grid - if (grid_method == GRID_UNION) call find_energy_index(p % E) - ! Determine if this material has S(a,b) tables check_sab = (mat % n_sab > 0) @@ -142,8 +137,10 @@ contains integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy - integer :: i_grid ! index on nuclide energy grid - real(8) :: f ! interp factor on nuclide energy grid + integer :: i_grid ! index on nuclide energy grid + integer :: i_low, i_high ! bounding indices from logarithmic mapping + integer :: u ! index into logarithmic mapping array + real(8) :: f ! interp factor on nuclide energy grid type(Nuclide), pointer, save :: nuc => null() !$omp threadprivate(nuc) @@ -152,16 +149,29 @@ contains ! Determine index on nuclide energy grid select case (grid_method) - case (GRID_UNION) - ! If we're using the unionized grid with pointers, finding the index on - ! the nuclide energy grid is as simple as looking up the pointer + case (GRID_LOGARITHM) + ! Determine the energy grid index using a logarithmic mapping to reduce + ! the energy range over which a binary search needs to be performed - i_grid = nuc % grid_index(union_grid_index) + if (E < nuc % energy(1)) then + i_grid = 1 + elseif (E > nuc % energy(nuc % n_grid)) then + i_grid = nuc % n_grid - 1 + else + ! Determine bounding indices based on which equal log-spaced interval + ! the energy is in + u = int(log(E/1.0e-11_8)/log_spacing) + i_low = nuc % grid_index(u) + i_high = nuc % grid_index(u + 1) + 1 + + ! Perform binary search over reduced range + i_grid = binary_search(nuc % energy(i_low:i_high), & + i_high - i_low + 1, E) + i_low - 1 + end if case (GRID_NUCLIDE) - ! If we're not using the unionized grid, we have to do a binary search on - ! the nuclide energy grid in order to determine which points to - ! interpolate between + ! Perform binary search on the nuclide energy grid in order to determine + ! which points to interpolate between if (E < nuc % energy(1)) then i_grid = 1 @@ -198,7 +208,7 @@ contains micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) & + f * nuc % total(i_grid+1) - ! Calculate microscopic nuclide total cross section + ! Calculate microscopic nuclide elastic cross section micro_xs(i_nuclide) % elastic = (ONE - f) * nuc % elastic(i_grid) & + f * nuc % elastic(i_grid+1) @@ -496,27 +506,6 @@ contains end subroutine calculate_urr_xs -!=============================================================================== -! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain -! energy -!=============================================================================== - - subroutine find_energy_index(E) - - real(8), intent(in) :: E ! energy of particle - - ! if particle's energy is outside of energy grid range, set to first or last - ! index. Otherwise, do a binary search through the union energy grid. - if (E < e_grid(1)) then - union_grid_index = 1 - elseif (E > e_grid(n_grid)) then - union_grid_index = n_grid - 1 - else - union_grid_index = binary_search(e_grid, n_grid, E) - end if - - end subroutine find_energy_index - !=============================================================================== ! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section ! for a given nuclide at the trial relative energy used in resonance scattering diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 4c4061a20..234a856c1 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -1,152 +1,62 @@ module energy_grid - use constants, only: MAX_LINE_LEN use global - use list_header, only: ListReal - use output, only: write_message implicit none + integer :: grid_method ! how to treat the energy grid + integer :: n_log_bins ! number of bins for logarithmic grid + real(8) :: log_spacing ! spacing on logarithmic grid + contains !=============================================================================== -! UNIONIZED_GRID creates a single unionized energy grid combined from each -! nuclide of each material. Right now, the grid for each nuclide is added into a -! linked list one at a time with an effective insertion sort. Could be done with -! a hash for all energy points and then a quicksort at the end (what hash -! function to use?) +! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding +! indices on a nuclide energy grid !=============================================================================== - subroutine unionized_grid() + subroutine logarithmic_grid() - integer :: i ! index in nuclides array - type(ListReal), pointer :: list => null() - type(Nuclide), pointer :: nuc => null() - - call write_message("Creating unionized energy grid...", 5) - - ! Add grid points for each nuclide in the problem - do i = 1, n_nuclides_total - nuc => nuclides(i) - call add_grid_points(list, nuc % energy) - end do - - ! Set size of unionized energy grid - n_grid = list % size() - - ! create allocated array from linked list - allocate(e_grid(n_grid)) - do i = 1, n_grid - e_grid(i) = list % get_item(i) - end do - - ! delete linked list and dictionary - call list % clear() - deallocate(list) - - ! Set pointers to unionized energy grid for each nuclide - call grid_pointers() - - end subroutine unionized_grid - -!=============================================================================== -! ADD_GRID_POINTS adds energy points from the 'energy' array into a linked list -! of points already stored from previous arrays. -!=============================================================================== - - subroutine add_grid_points(list, energy) - - type(ListReal), pointer :: list - real(8), intent(in) :: energy(:) - - integer :: i ! index in energy array - integer :: n ! size of energy array - integer :: current ! current index - real(8) :: E ! actual energy value - - i = 1 - n = size(energy) - - ! If the original list is empty, we need to allocate the first element and - ! store first energy point - if (.not. associated(list)) then - allocate(list) - do i = 1, n - call list % append(energy(i)) - end do - return - end if - - ! Set current index to beginning of the list - current = 1 - - do while (i <= n) - E = energy(i) - - ! If we've reached the end of the grid energy list, add the remaining - ! energy points to the end - if (current > list % size()) then - ! Finish remaining energies - do while (i <= n) - call list % append(energy(i)) - i = i + 1 - end do - exit - end if - - if (E < list % get_item(current)) then - - ! Insert new energy in this position - call list % insert(current, E) - - ! Advance index in linked list and in new energy grid - i = i + 1 - current = current + 1 - - elseif (E == list % get_item(current)) then - ! Found the exact same energy, no need to store duplicates so just - ! skip and move to next index - i = i + 1 - current = current + 1 - else - current = current + 1 - end if - - end do - - end subroutine add_grid_points - -!=============================================================================== -! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link -! each point on the nuclide energy grid to one on the unionized energy grid -!=============================================================================== - - subroutine grid_pointers() - - integer :: i ! loop index for nuclides - integer :: j ! loop index for nuclide energy grid - integer :: index_e ! index on union energy grid - real(8) :: union_energy ! energy on union grid - real(8) :: energy ! energy on nuclide grid + integer :: i, j, k ! Loop indices + integer :: M ! Number of equally log-spaced bins + real(8) :: E_max ! Maximum energy in MeV + real(8) :: E_min ! Minimum energy in MeV + real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid type(Nuclide), pointer :: nuc => null() + ! Set minimum/maximum energies + E_max = 20.0_8 + E_min = 1.0e-11_8 + + ! Determine equal-logarithmic energy spacing + M = n_log_bins + log_spacing = log(E_max/E_min)/M + + ! Create equally log-spaced energy grid + allocate(umesh(0:M)) + umesh(:) = [(i*log_spacing, i=0, M)] + do i = 1, n_nuclides_total + ! Allocate logarithmic mapping for nuclide nuc => nuclides(i) - allocate(nuc % grid_index(n_grid)) + allocate(nuc % grid_index(0:M)) - index_e = 1 - energy = nuc % energy(index_e) - - do j = 1, n_grid - union_energy = e_grid(j) - if (union_energy >= energy .and. index_e < nuc % n_grid) then - index_e = index_e + 1 - energy = nuc % energy(index_e) - end if - nuc % grid_index(j) = index_e - 1 + ! Determine corresponding indices in nuclide grid to energies on + ! equal-logarithmic grid + j = 1 + do k = 0, M - 1 + do while (log(nuc%energy(j + 1)/E_min) <= umesh(k)) + j = j + 1 + end do + nuc % grid_index(k) = j end do + + ! Set the last point explicitly so that we don't have out-of-bounds issues + nuc % grid_index(M) = size(nuc % energy) - 1 end do - end subroutine grid_pointers + deallocate(umesh) + + end subroutine logarithmic_grid end module energy_grid diff --git a/src/geometry.F90 b/src/geometry.F90 index a103f4011..0b6aea38d 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -348,10 +348,9 @@ contains ! Score to global leakage tally if (tallies_on) then -!$omp critical +!$omp atomic global_tallies(LEAKAGE) % value = & - global_tallies(LEAKAGE) % value + p % wgt -!$omp end critical + global_tallies(LEAKAGE) % value + p % wgt end if ! Display message @@ -586,7 +585,7 @@ contains return end if end if - + end subroutine cross_surface !=============================================================================== @@ -823,7 +822,7 @@ contains if (quad < ZERO) then ! no intersection with cylinder - d = INFINITY + d = INFINITY elseif (on_surface) then ! particle is on the cylinder, thus one distance is @@ -872,7 +871,7 @@ contains if (quad < ZERO) then ! no intersection with cylinder - d = INFINITY + d = INFINITY elseif (on_surface) then ! particle is on the cylinder, thus one distance is @@ -921,7 +920,7 @@ contains if (quad < ZERO) then ! no intersection with cylinder - d = INFINITY + d = INFINITY elseif (on_surface) then ! particle is on the cylinder, thus one distance is @@ -968,7 +967,7 @@ contains if (quad < ZERO) then ! no intersection with sphere - d = INFINITY + d = INFINITY elseif (on_surface) then ! particle is on the sphere, thus one distance is @@ -1015,7 +1014,7 @@ contains if (quad < ZERO) then ! no intersection with cone - d = INFINITY + d = INFINITY elseif (on_surface) then ! particle is on the cone, thus one distance is positive/negative @@ -1034,7 +1033,7 @@ contains d = (-k - quad)/a b = (-k + quad)/a - ! determine the smallest positive solution + ! determine the smallest positive solution if (d < ZERO) then if (b > ZERO) then d = b @@ -1064,7 +1063,7 @@ contains if (quad < ZERO) then ! no intersection with cone - d = INFINITY + d = INFINITY elseif (on_surface) then ! particle is on the cone, thus one distance is positive/negative @@ -1083,7 +1082,7 @@ contains d = (-k - quad)/a b = (-k + quad)/a - ! determine the smallest positive solution + ! determine the smallest positive solution if (d < ZERO) then if (b > ZERO) then d = b @@ -1113,7 +1112,7 @@ contains if (quad < ZERO) then ! no intersection with cone - d = INFINITY + d = INFINITY elseif (on_surface) then ! particle is on the cone, thus one distance is positive/negative @@ -1132,7 +1131,7 @@ contains d = (-k - quad)/a b = (-k + quad)/a - ! determine the smallest positive solution + ! determine the smallest positive solution if (d < ZERO) then if (b > ZERO) then d = b @@ -1604,9 +1603,8 @@ contains ! Increment number of lost particles p % alive = .false. -!$omp critical +!$omp atomic n_lost_particles = n_lost_particles + 1 -!$omp end critical ! Abort the simulation if the maximum number of lost particles has been ! reached diff --git a/src/global.F90 b/src/global.F90 index 1953e9291..004c43746 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -62,7 +62,7 @@ module global ! Cross section arrays type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables - type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings + type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings ! Cross section caches type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide @@ -77,11 +77,6 @@ module global type(DictCharInt) :: sab_dict type(DictCharInt) :: xs_listing_dict - ! Unionized energy grid - integer :: grid_method ! how to treat the energy grid - integer :: n_grid ! number of points on unionized grid - real(8), allocatable :: e_grid(:) ! energies on unionized grid - ! Unreoslved resonance probablity tables logical :: urr_ptables_on = .true. @@ -119,7 +114,7 @@ module global ! 2) track-length estimate of k-eff ! 3) leakage fraction - type(TallyResult), target :: global_tallies(N_GLOBAL_TALLIES) + type(TallyResult), allocatable, target :: global_tallies(:) ! Tally map structure type(TallyMap), allocatable :: tally_maps(:) @@ -223,7 +218,6 @@ module global type(Timer) :: time_total ! timer for total run type(Timer) :: time_initialize ! timer for initialization type(Timer) :: time_read_xs ! timer for reading cross sections - type(Timer) :: time_unionize ! timer for unionizing energy grid type(Timer) :: time_bank ! timer for fission bank synchronization type(Timer) :: time_bank_sample ! timer for fission bank sampling type(Timer) :: time_bank_sendrecv ! timer for fission bank SEND/RECV @@ -300,7 +294,7 @@ module global logical :: write_initial_source = .false. ! ============================================================================ - ! CMFD VARIABLES + ! CMFD VARIABLES ! Main object type(cmfd_type) :: cmfd @@ -310,11 +304,11 @@ module global ! CMFD communicator integer :: cmfd_comm - + ! Timing objects type(Timer) :: time_cmfd ! timer for whole cmfd calculation type(Timer) :: time_cmfdbuild ! timer for matrix build - type(Timer) :: time_cmfdsolve ! timer for solver + type(Timer) :: time_cmfdsolve ! timer for solver ! Flag for active core map logical :: cmfd_coremap = .false. @@ -390,7 +384,7 @@ module global ! RESONANCE SCATTERING VARIABLES logical :: treat_res_scat = .false. ! is resonance scattering treated? - integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers + integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, & @@ -399,14 +393,14 @@ module global contains !=============================================================================== -! FREE_MEMORY deallocates and clears all global allocatable arrays in the +! FREE_MEMORY deallocates and clears all global allocatable arrays in the ! program !=============================================================================== subroutine free_memory() - + integer :: i ! Loop Index - + ! Deallocate cells, surfaces, materials if (allocated(cells)) deallocate(cells) if (allocated(universes)) deallocate(universes) @@ -449,6 +443,7 @@ contains if (allocated(entropy_p)) deallocate(entropy_p) ! Deallocate tally-related arrays + if (allocated(global_tallies)) deallocate(global_tallies) if (allocated(meshes)) deallocate(meshes) if (allocated(tallies)) then ! First call the clear routines @@ -461,9 +456,6 @@ contains 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) @@ -488,7 +480,7 @@ contains ! Deallocate track_identifiers if (allocated(track_identifiers)) deallocate(track_identifiers) - + ! Deallocate dictionaries call cell_dict % clear() call universe_dict % clear() @@ -525,7 +517,7 @@ contains if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width) deallocate(ufs_mesh) end if - + end subroutine free_memory end module global diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index faf6df6e2..e10ddc982 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -68,7 +68,7 @@ contains end if ! Terminate access to the file. - call su % file_close() + call su % file_close() end subroutine hdf5_write_summary @@ -81,7 +81,7 @@ contains ! Write version information call su % write_data(VERSION_MAJOR, "version_major") call su % write_data(VERSION_MINOR, "version_minor") - call su % write_data(VERSION_RELEASE, "version_release") + call su % write_data(VERSION_RELEASE, "version_release") ! Write current date and time call su % write_data(time_stamp(), "date_and_time") @@ -89,7 +89,7 @@ contains ! Write MPI information call su % write_data(n_procs, "n_procs") call su % write_attribute_string("n_procs", "description", & - "Number of MPI processes") + "Number of MPI processes") end subroutine hdf5_write_header @@ -143,7 +143,7 @@ contains call su % write_data("universe", "fill_type", & group="geometry/cells/cell " // trim(to_str(c % id))) call su % write_data(universes(c % fill) % id, "material", & - group="geometry/cells/cell " // trim(to_str(c % id))) + group="geometry/cells/cell " // trim(to_str(c % id))) case (CELL_LATTICE) call su % write_data("lattice", "fill_type", & group="geometry/cells/cell " // trim(to_str(c % id))) @@ -422,7 +422,7 @@ contains group="materials/material " // trim(to_str(m % id))) call su % write_data(m % i_sab_tables, "i_sab_tables", & length=m % n_sab, & - group="materials/material " // trim(to_str(m % id))) + group="materials/material " // trim(to_str(m % id))) end if end do @@ -710,8 +710,6 @@ contains group="timing") call su % write_data(time_read_xs % elapsed, "time_read_xs", & group="timing") - call su % write_data(time_unionize % elapsed, "time_unionize", & - group="timing") call su % write_data(time_transport % elapsed, "time_transport", & group="timing") call su % write_data(time_bank % elapsed, "time_bank", & @@ -736,8 +734,6 @@ contains "Total time elapsed for initialization (s)", group="timing") call su % write_attribute_string("time_read_xs", "description", & "Time reading cross-section libraries (s)", group="timing") - call su % write_attribute_string("time_unionize", "description", & - "Time unionizing energy grid (s)", group="timing") call su % write_attribute_string("time_transport", "description", & "Time in transport only (s)", group="timing") call su % write_attribute_string("time_bank", "description", & diff --git a/src/initialize.F90 b/src/initialize.F90 index 35d7c3ea5..9909ae620 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -4,7 +4,7 @@ module initialize use bank_header, only: Bank use constants use dict_header, only: DictIntInt, ElemKeyValueII - use energy_grid, only: unionized_grid + use energy_grid, only: logarithmic_grid, grid_method use error, only: fatal_error, warning use geometry, only: neighbor_lists use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& @@ -109,11 +109,9 @@ contains ! Create linked lists for multiple instances of the same nuclide call same_nuclide_list() - ! Construct unionized energy grid from cross-sections - if (grid_method == GRID_UNION) then - call time_unionize % start() - call unionized_grid() - call time_unionize % stop() + ! Construct logarithmic energy grid for cross-sections + if (grid_method == GRID_LOGARITHM) then + call logarithmic_grid() end if ! Allocate and setup tally stride, matching_bins, and tally maps diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 025ab80e5..9bf1b7119 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3,6 +3,7 @@ module input_xml use cmfd_input, only: configure_cmfd use constants use dict_header, only: DictIntInt, ElemKeyValueCI + use energy_grid, only: grid_method, n_log_bins use error, only: fatal_error, warning use geometry_header, only: Cell, Surface, Lattice, RectLattice, HexLattice use global @@ -207,19 +208,30 @@ contains if (check_for_node(doc, "energy_grid")) then call get_node_value(doc, "energy_grid", temp_str) else - temp_str = 'union' + temp_str = 'logarithm' end if select case (trim(temp_str)) case ('nuclide') grid_method = GRID_NUCLIDE case ('union') - grid_method = GRID_UNION - case ('lethargy') - call fatal_error("Lethargy mapped energy grid not yet supported.") + call fatal_error("Union energy grid is no longer supported.") + case ('logarithm', 'logarithmic', 'log') + grid_method = GRID_LOGARITHM case default call fatal_error("Unknown energy grid method: " // trim(temp_str)) end select + ! Number of bins for logarithmic grid + if (check_for_node(doc, "log_grid_bins")) then + call get_node_value(doc, "log_grid_bins", n_log_bins) + if (n_log_bins < 1) then + call fatal_error("Number of bins for logarithmic grid must be & + &greater than zero.") + end if + else + n_log_bins = 8000 + end if + ! Verbosity if (check_for_node(doc, "verbosity")) then call get_node_ptr(doc, "verbosity", node_verb) @@ -524,11 +536,11 @@ contains ! Copy dimensions call get_node_array(node_entropy, "dimension", entropy_mesh % dimension) - + ! Calculate width entropy_mesh % width = (entropy_mesh % upper_right - & entropy_mesh % lower_left) / entropy_mesh % dimension - + end if ! Turn on Shannon entropy calculation @@ -2516,7 +2528,7 @@ contains j = j + 1 ! Get the input string in scores(l) but if score is one of the moment ! scores then strip off the n and store it as an integer to be used - ! later. Then perform the select case on this modified (number + ! later. Then perform the select case on this modified (number ! removed) string score_name = sarray(l) do imomstr = 1, size(MOMENT_STRS) @@ -3007,7 +3019,7 @@ contains ! Copy plot cell universe level if (check_for_node(node_plot, "level")) then call get_node_value(node_plot, "level", pl % level) - + if (pl % level < 0) then call fatal_error("Bad universe level in plot " & &// trim(to_str(pl % id))) @@ -3113,7 +3125,7 @@ contains call warning("Meshlines ignored in voxel plot " & &// trim(to_str(pl % id))) end if - + select case(n_meshlines) case (0) ! Skip if no meshlines are specified @@ -3121,7 +3133,7 @@ contains ! Get pointer to meshlines call get_list_item(node_meshline_list, 1, node_meshlines) - + ! Check mesh type if (check_for_node(node_meshlines, "meshtype")) then call get_node_value(node_meshlines, "meshtype", meshtype) @@ -3129,7 +3141,7 @@ contains call fatal_error("Must specify a meshtype for meshlines & &specification in plot " // trim(to_str(pl % id))) end if - + ! Ensure that there is a linewidth for this meshlines specification if (check_for_node(node_meshlines, "linewidth")) then call get_node_value(node_meshlines, "linewidth", & @@ -3141,19 +3153,19 @@ contains ! Check for color if (check_for_node(node_meshlines, "color")) then - + ! Check and make sure 3 values are specified for RGB if (get_arraysize_double(node_meshlines, "color") /= 3) then call fatal_error("Bad RGB for meshlines color in plot " & &// trim(to_str(pl % id))) end if - + call get_node_array(node_meshlines, "color", & pl % meshlines_color % rgb) else - + pl % meshlines_color % rgb = (/ 0, 0, 0 /) - + end if ! Set mesh based on type @@ -3164,7 +3176,7 @@ contains call fatal_error("No UFS mesh for meshlines on plot " & &// trim(to_str(pl % id))) end if - + pl % meshlines_mesh => ufs_mesh case ('cmfd') @@ -3180,17 +3192,17 @@ contains pl % meshlines_mesh => meshes(i_mesh) case ('entropy') - + if (.not. associated(entropy_mesh)) then call fatal_error("No entropy mesh for meshlines on plot " & &// trim(to_str(pl % id))) end if - + if (.not. allocated(entropy_mesh % dimension)) then call fatal_error("No dimension specified on entropy mesh & &for meshlines on plot " // trim(to_str(pl % id))) end if - + pl % meshlines_mesh => entropy_mesh case ('tally') @@ -3225,9 +3237,9 @@ contains call fatal_error("Mutliple meshlines specified in plot " & &// trim(to_str(pl % id))) end select - + end if - + ! Deal with masks call get_node_list(node_plot, "mask", node_mask_list) n_masks = get_list_size(node_mask_list) diff --git a/src/output.F90 b/src/output.F90 index d742baab1..ff669af08 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -51,7 +51,7 @@ contains ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Copyright: 2011-2014 Massachusetts Institute of Technology' + ' Copyright: 2011-2015 Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT=*) & ' License: http://mit-crpg.github.io/openmc/license.html' write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') & @@ -158,7 +158,7 @@ contains if (master) then write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') & "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE - write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2013 & + write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2015 & &Massachusetts Institute of Technology" write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at & &" @@ -1208,7 +1208,7 @@ contains call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2011-2013 Massachusetts Institute of Technology" + "Copyright: 2011-2015 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 @@ -1261,12 +1261,6 @@ contains end do end if - ! print summary of unionized energy grid - call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY) - write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid)) - write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(& - n_grid*n_nuclides_total*4)) // " bytes" - ! print summary of variance reduction call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY) if (survival_biasing) then @@ -1535,7 +1529,6 @@ contains ! display time elapsed for various sections write(ou,100) "Total time for initialization", time_initialize % elapsed write(ou,100) " Reading cross sections", time_read_xs % elapsed - write(ou,100) " Unionizing energy grid", time_unionize % elapsed write(ou,100) "Total time in simulation", time_inactive % elapsed + & time_active % elapsed write(ou,100) " Time in transport only", time_transport % elapsed diff --git a/src/physics.F90 b/src/physics.F90 index 8fa0575cd..9b32c64d9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -231,7 +231,7 @@ contains prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + f*(rxn%sigma(i_grid - rxn%threshold + 2))) - ! Create fission bank sites if fission occus + ! Create fission bank sites if fission occurs if (prob > cutoff) exit FISSION_REACTION_LOOP end do FISSION_REACTION_LOOP @@ -256,22 +256,19 @@ contains p % last_wgt = p % wgt ! Score implicit absorption estimate of keff -!$omp critical +!$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 -!$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 critical +!$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 -!$omp end critical p % alive = .false. p % event = EVENT_ABSORB @@ -388,7 +385,7 @@ contains + f*(rxn%sigma(i_grid - rxn%threshold + 2))) end do - ! Perform collision physics for inelastics scattering + ! Perform collision physics for inelastic scattering call inelastic_scatter(nuc, rxn, p % E, p % coord0 % uvw, & p % mu, p % wgt) p % event_MT = rxn % MT @@ -797,7 +794,7 @@ contains sampling_scheme = 'cxs' end if - ! otherwise, use free gas model + ! otherwise, use free gas model else if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then v_target = ZERO @@ -859,7 +856,7 @@ contains m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) - + ! get max 0K xs value over range of practical relative energies xs_max = max(xs_low, & & maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up) @@ -972,7 +969,7 @@ contains case default call fatal_error("Not a recognized resonance scattering treatment!") end select - + end subroutine sample_target_velocity !=============================================================================== @@ -1835,14 +1832,15 @@ contains lc = 2 + 2*NR + 2*NE U = edist % data(lc + 1) + y = (E_in - U)/T + v = 1 - exp(-y) + ! sample outgoing energy based on evaporation spectrum probability ! density function n_sample = 0 do - r1 = prn() - r2 = prn() - E_out = -T * log(r1*r2) - if (E_out <= E_in - U) exit + x = -log((1 - v*prn())*(1 - v*prn())) + if (x <= y) exit ! check for large number of rejections n_sample = n_sample + 1 @@ -1852,6 +1850,8 @@ contains end if end do + E_out = x*T + case (11) ! ======================================================================= ! ENERGY-DEPENDENT WATT SPECTRUM diff --git a/src/plot.F90 b/src/plot.F90 index ddc4641b2..d9f6b55e2 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -272,11 +272,11 @@ contains outrange(1) = int(frac * real(img % width, 8)) frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer) outrange(2) = int(frac * real(img % width, 8)) - - frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner) - inrange(1) = int(frac * real(img % height, 8)) + frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner) - inrange(2) = int(frac * real(img % height, 8)) + inrange(1) = int((1. - frac) * real(img % height, 8)) + frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner) + inrange(2) = int((1. - frac) * real(img % height, 8)) ! draw lines do out_ = outrange(1), outrange(2) diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 707ecbbc8..a7b92e173 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -3,19 +3,19 @@ element settings { ( element eigenvalue { - (element batches { xsd:positiveInteger } | + (element batches { xsd:positiveInteger } | attribute batches { xsd:positiveInteger }) & - (element inactive { xsd:nonNegativeInteger } | + (element inactive { xsd:nonNegativeInteger } | attribute inactive { xsd:nonNegativeInteger }) & - (element particles { xsd:positiveInteger } | + (element particles { xsd:positiveInteger } | attribute particles { xsd:positiveInteger }) & - (element generations_per_batch { xsd:positiveInteger } | + (element generations_per_batch { xsd:positiveInteger } | attribute generations_per_batch { xsd:positiveInteger })? } | element fixed_source { - (element batches { xsd:positiveInteger } | + (element batches { xsd:positiveInteger } | attribute batches { xsd:positiveInteger }) & - (element particles { xsd:positiveInteger } | + (element particles { xsd:positiveInteger } | attribute particles { xsd:positiveInteger }) } ) & @@ -27,17 +27,19 @@ element settings { (element weight_avg { xsd:double } | attribute weight_avg { xsd:double })? }? & - element energy_grid { ( "nuclide" | "union" | "lethargy" ) }? & + element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" ) }? & element entropy { - (element dimension { list { xsd:int+ } } | + (element dimension { list { xsd:int+ } } | attribute dimension { list { xsd:int+ } })? & - (element lower_left { list { xsd:double+ } } | + (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) & - (element upper_right { list { xsd:double+ } } | + (element upper_right { list { xsd:double+ } } | attribute upper_right { list { xsd:double+ } }) }? & + element log_grid_bins { xsd:positiveInteger }? & + element natural_elements { xsd:string { maxLength = "20" } }? & element no_reduce { xsd:boolean }? & @@ -93,7 +95,7 @@ element settings { ( (element batches { list { xsd:positiveInteger+ } } | attribute batches { list { xsd:positiveInteger+ } }) | - (element interval { xsd:positiveInteger } | + (element interval { xsd:positiveInteger } | attribute interval { xsd:positiveInteger }) ) }? & @@ -102,12 +104,12 @@ element settings { ( (element batches { list { xsd:positiveInteger+ } } | attribute batches { list { xsd:positiveInteger+ } }) | - (element interval { xsd:positiveInteger } | + (element interval { xsd:positiveInteger } | attribute interval { xsd:positiveInteger }) )? & - (element separate { xsd:boolean } | + (element separate { xsd:boolean } | attribute separate { xsd:boolean })? & - (element write { xsd:boolean } | + (element write { xsd:boolean } | attribute write { xsd:boolean })? & (element overwrite_latest { xsd:boolean} | attribute overwrite_latest {xsd:boolean})? @@ -124,11 +126,11 @@ element settings { element verbosity { xsd:positiveInteger }? & element uniform_fs{ - (element dimension { list { xsd:positiveInteger+ } } | + (element dimension { list { xsd:positiveInteger+ } } | attribute dimension { list { xsd:positiveInteger+ } }) & - (element lower_left { list { xsd:double+ } } | + (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) & - (element upper_right { list { xsd:double+ } } | + (element upper_right { list { xsd:double+ } } | attribute upper_right { list { xsd:double+ } }) }? & diff --git a/src/source.F90 b/src/source.F90 index 07c1bfb41..3a12a3f71 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -211,6 +211,9 @@ contains case (SRC_ENERGY_MONO) ! Monoenergtic source site % E = external_source % params_energy(1) + if (site % E >= 20) then + call fatal_error("Source energies above 20 MeV not allowed.") + end if case (SRC_ENERGY_MAXWELL) a = external_source % params_energy(1) diff --git a/src/string.F90 b/src/string.F90 index cff69bb16..5071493f8 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -148,7 +148,7 @@ contains end function concatenate !=============================================================================== -! LOWER_CASE converts a string to all lower case characters +! TO_LOWER converts a string to all lower case characters !=============================================================================== elemental function to_lower(word) result(word_lower) @@ -171,7 +171,7 @@ contains end function to_lower !=============================================================================== -! UPPER_CASE converts a string to all upper case characters +! TO_UPPER converts a string to all upper case characters !=============================================================================== elemental function to_upper(word) result(word_upper) diff --git a/src/tally.F90 b/src/tally.F90 index 87f5a2cab..c96f1ecd6 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -281,10 +281,9 @@ contains ! get the score and tally it score = last_wgt * calc_pn(n, mu) -!$omp critical +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score -!$omp end critical end do j = j + t % moment_order(j) cycle SCORE_LOOP @@ -347,10 +346,9 @@ contains ! get the score and tally it score = wgt * calc_pn(n, mu) -!$omp critical +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score -!$omp end critical end do j = j + t % moment_order(j) cycle SCORE_LOOP @@ -542,10 +540,9 @@ contains end select ! Add score to tally -!$omp critical +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score -!$omp end critical end do SCORE_LOOP @@ -617,10 +614,9 @@ contains i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! Add score to tally -!$omp critical +!$omp atomic 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 @@ -1015,10 +1011,9 @@ contains end if ! Add score to tally -!$omp critical +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score -!$omp end critical end do SCORE_LOOP @@ -1214,10 +1209,9 @@ contains end select ! Add score to tally -!$omp critical +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score -!$omp end critical end do SCORE_LOOP @@ -1366,10 +1360,9 @@ contains end select ! Add score to tally -!$omp critical +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score -!$omp end critical end do MATERIAL_SCORE_LOOP @@ -1787,10 +1780,9 @@ contains end if ! Add score to tally -!$omp critical +!$omp atomic t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score -!$omp end critical end do SCORE_LOOP @@ -1845,8 +1837,10 @@ contains p % coord % universe, i_tally) case (FILTER_MATERIAL) - matching_bins(i) = get_next_bin(FILTER_MATERIAL, & - p % material, i_tally) + if (p % material /= MATERIAL_VOID) then + matching_bins(i) = get_next_bin(FILTER_MATERIAL, & + p % material, i_tally) + endif case (FILTER_CELL) ! determine next cell bin @@ -2021,10 +2015,9 @@ contains matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 -!$omp critical +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt -!$omp end critical end if end do else @@ -2035,10 +2028,9 @@ contains matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 -!$omp critical +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt -!$omp end critical end if end do end if @@ -2053,10 +2045,9 @@ contains matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 -!$omp critical +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt -!$omp end critical end if end do else @@ -2067,10 +2058,9 @@ contains matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 -!$omp critical +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt -!$omp end critical end if end do end if @@ -2085,10 +2075,9 @@ contains matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 -!$omp critical +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt -!$omp end critical end if end do else @@ -2099,10 +2088,9 @@ contains matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 -!$omp critical +!$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt -!$omp end critical end if end do end if @@ -2224,10 +2212,9 @@ contains end if ! Add to surface current tally -!$omp critical +!$omp atomic 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/tally_initialize.F90 b/src/tally_initialize.F90 index 23ab5b563..b7dc5ed98 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -19,6 +19,9 @@ contains subroutine configure_tallies() + ! Allocate global tallies + allocate(global_tallies(N_GLOBAL_TALLIES)) + call setup_tally_arrays() call setup_tally_maps() diff --git a/src/tracking.F90 b/src/tracking.F90 index 87e6e0f69..c939f7991 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -62,9 +62,8 @@ contains n_event = 0 ! Add paricle's starting weight to count for normalizing tallies later -!$omp critical +!$omp atomic 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 @@ -113,11 +112,10 @@ contains call score_tracklength_tally(p, distance) ! Score track-length estimate of k-eff -!$omp critical +!$omp atomic 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 ! ==================================================================== @@ -141,11 +139,10 @@ contains ! PARTICLE HAS COLLISION ! Score collision estimate of keff -!$omp critical +!$omp atomic 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