diff --git a/docs/source/_images/uniongrid.svg b/docs/source/_images/uniongrid.svg deleted file mode 100644 index 27c3922fdf..0000000000 --- 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 a193dde96d..db2e5156ee 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 1515ffa891..1ba29b79ef 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -421,7 +421,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 54dc913455..db6cfd89e9 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 9dd6be6d1d..3fb9036d39 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 f7e70ef209..c834bbc92e 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 @@ -757,7 +769,10 @@ Each ```` element can have the following attributes or sub-elements: is on the negative side of surface 3 and the positive side of surface 5, the bounding surfaces would be given as "-3 5". - *Default*: None + .. note:: The surface attribute/element can be omitted to make a cell fill + its entire universe. + + *Default*: No surfaces :rotation: If the cell is filled with a universe, this element specifies the angles in @@ -1019,6 +1034,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", @@ -1170,7 +1195,7 @@ implemented in openMC: ```` Element ------------------ -Each plot must contain a combination of the following attributes or +Each plot is specified by a combination of the following attributes or sub-elements: :id: @@ -1191,6 +1216,18 @@ sub-elements: *Default*: ``cell`` + :level: + Universe depth to plot at (optional). This parameter controls how many + universe levels deep to pull cell and material ids from when setting plot + colors. If a given location does not have as many levels as specified, + colors will be taken from the lowest level at that location. For example, if + ``level`` is set to zero colors will be taken from top-level (universe zero) + cells only. However, if ``level`` is set to 1 colors will be taken from + cells in universes that fill top-level fill-cells, and from top-level cells + that contain materials. + + *Default*: Whatever the deepest universe is in the model + :origin: Specifies the (x,y,z) coordinate of the center of the plot. Should be three floats separated by spaces. @@ -1557,3 +1594,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 76492dfe80..afed796e29 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 af12b97be2..7e30bd07ea 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 8e4c19cccd..7c95dfb37f 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 @@ -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 4c4061a201..234a856c1e 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 043ded05b5..c79f0e6119 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -11,7 +11,7 @@ module geometry use tally, only: score_surface_current implicit none - + contains !=============================================================================== @@ -59,7 +59,8 @@ contains end if end do SURFACE_LOOP - ! If we've reached here, then the sense matched on every surface + ! If we've reached here, then the sense matched on every surface or there + ! are no surfaces. in_cell = .true. end function simple_cell_contains @@ -264,13 +265,13 @@ contains lattice_edge = .true. end if end if - + if (lattice_edge) then - + ! In this case the neutron is leaving the lattice, so we move it ! out, remove all lower coordinate levels and then search from ! universe 0. - + p % coord => p % coord0 call deallocate_coord(p % coord % next) @@ -287,7 +288,7 @@ contains p % last_material = p % material p % material = c % material - ! We'll still make a new coordinate for the particle, as + ! We'll still make a new coordinate for the particle, as ! distance_to_boundary will still need to track through lattice ! widths even though there's nothing in them but this material @@ -406,10 +407,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 end if ! Display message @@ -644,7 +644,7 @@ contains return end if end if - + end subroutine cross_surface !=============================================================================== @@ -906,7 +906,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 @@ -955,7 +955,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 @@ -1004,7 +1004,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 @@ -1051,7 +1051,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 @@ -1098,7 +1098,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 @@ -1117,7 +1117,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 @@ -1147,7 +1147,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 @@ -1166,7 +1166,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 @@ -1196,7 +1196,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 @@ -1215,7 +1215,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 @@ -1273,7 +1273,7 @@ contains ! logic here checks whether the relative difference is within floating ! point precision. - if (d < dist) then + if (d < dist) then if (abs(d - dist)/dist >= FP_REL_PRECISION) then dist = d if (u > 0) then @@ -1565,9 +1565,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 f4c50b8330..ec11216608 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 c9e367200c..fc4af5bc7b 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -67,7 +67,7 @@ contains end if ! Terminate access to the file. - call su % file_close() + call su % file_close() end subroutine hdf5_write_summary @@ -80,7 +80,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") @@ -88,7 +88,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 @@ -144,12 +144,12 @@ 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))) call su % write_data(lattices(c % fill) % id, "lattice", & - group="geometry/cells/cell " // trim(to_str(c % id))) + group="geometry/cells/cell " // trim(to_str(c % id))) end select ! Write list of bounding surfaces @@ -303,7 +303,7 @@ contains else n_z = 1 end if - + ! Write lattice universes allocate(lattice_universes(n_x, n_y, n_z)) do j = 1, n_x @@ -371,7 +371,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 @@ -659,8 +659,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", & @@ -685,8 +683,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 a418b3ec32..1ab34ca0bf 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, BASE_UNIVERSE @@ -108,11 +108,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 e992a0eae0..e5c21bdc1c 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 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 @@ -1011,17 +1023,18 @@ contains call fatal_error("Cannot specify material and fill simultaneously") end if - ! Check to make sure that surfaces were specified - if (.not. check_for_node(node_cell, "surfaces")) then - call fatal_error("No surfaces specified for cell " & - &// trim(to_str(c % id))) - end if - ! Allocate array for surfaces and copy - n = get_arraysize_integer(node_cell, "surfaces") + if (check_for_node(node_cell, "surfaces")) then + n = get_arraysize_integer(node_cell, "surfaces") + else + n = 0 + end if c % n_surfaces = n - allocate(c % surfaces(n)) - call get_node_array(node_cell, "surfaces", c % surfaces) + + if (n > 0) then + allocate(c % surfaces(n)) + call get_node_array(node_cell, "surfaces", c % surfaces) + end if ! Rotation matrix if (check_for_node(node_cell, "rotation")) then @@ -2329,7 +2342,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) @@ -2730,16 +2743,14 @@ contains end select ! Set output file path - filename = "plot" + filename = trim(to_str(pl % id)) // "_plot" if (check_for_node(node_plot, "filename")) & call get_node_value(node_plot, "filename", filename) select case (pl % type) case (PLOT_TYPE_SLICE) - pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & - "_" // trim(filename) // ".ppm" + pl % path_plot = trim(path_input) // trim(filename) // ".ppm" case (PLOT_TYPE_VOXEL) - pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & - "_" // trim(filename) // ".voxel" + pl % path_plot = trim(path_input) // trim(filename) // ".voxel" end select ! Copy plot pixel size @@ -2819,6 +2830,18 @@ contains end if end if + ! 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))) + end if + else + pl % level = PLOT_LEVEL_LOWEST + end if + ! Copy plot color type and initialize all colors randomly temp_str = "cell" if (check_for_node(node_plot, "color")) & @@ -2916,7 +2939,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 @@ -2924,7 +2947,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) @@ -2932,7 +2955,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", & @@ -2944,19 +2967,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 @@ -2967,7 +2990,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') @@ -2983,17 +3006,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') @@ -3028,9 +3051,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 4d0ef7edd4..8deee0e5c7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1226,12 +1226,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 @@ -1428,6 +1422,12 @@ contains ! Plot id write(ou,100) "Plot ID:", trim(to_str(pl % id)) + ! Plot filename + write(ou,100) "Plot file:", trim(pl % path_plot) + + ! Plot level + write(ou,100) "Universe depth:", trim(to_str(pl % level)) + ! Plot type if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Plot Type:", "Slice" @@ -1494,7 +1494,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 8fa0575cdc..3439a83905 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 @@ -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 1297b746c4..d9f6b55e2b 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -7,7 +7,7 @@ module plot use global use mesh, only: get_mesh_indices use output, only: write_message - use particle_header, only: deallocate_coord, Particle + use particle_header, only: deallocate_coord, Particle, LocalCoord use plot_header use ppmlib, only: Image, init_image, allocate_image, & deallocate_image, set_pixel @@ -32,7 +32,7 @@ contains ! Display output message call write_message("Processing plot " // trim(to_str(pl % id)) & - &// "...", 5) + &// ": " // trim(pl % path_plot) // " ...", 5) if (pl % type == PLOT_TYPE_SLICE) then ! create 2d image @@ -58,7 +58,9 @@ contains integer, intent(out) :: id logical :: found_cell - type(Cell), pointer :: c => null() + integer :: level + type(Cell), pointer :: c => null() + type(LocalCoord), pointer :: coord => null() call deallocate_coord(p % coord0 % next) p % coord => p % coord0 @@ -66,6 +68,16 @@ contains call find_cell(p, found_cell) if (check_overlaps) call check_cell_overlap(p) + ! Loop through universes and stop on any specified level + level = 0 + coord => p % coord0 + do + if (level == pl % level) exit + if (.not. associated(coord % next)) exit + coord => coord % next + level = level + 1 + end do + if (.not. found_cell) then ! If no cell, revert to default color rgb = pl % not_found % rgb @@ -73,19 +85,23 @@ contains else if (pl % color_by == PLOT_COLOR_MATS) then ! Assign color based on material - c => cells(p % coord % cell) + c => cells(coord % cell) if (c % material == MATERIAL_VOID) then ! By default, color void cells white rgb = 255 id = -1 + else if (c % type == CELL_FILL) then + ! If we stopped on a middle universe level, treat as if not found + rgb = pl % not_found % rgb + id = -1 else rgb = pl % colors(c % material) % rgb id = materials(c % material) % id end if else if (pl % color_by == PLOT_COLOR_CELLS) then ! Assign color based on cell - rgb = pl % colors(p % coord % cell) % rgb - id = cells(p % coord % cell) % id + rgb = pl % colors(coord % cell) % rgb + id = cells(coord % cell) % id else rgb = 0 id = -1 @@ -256,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/plot_header.F90 b/src/plot_header.F90 index 68eb89a11a..7ce0aa4500 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -27,6 +27,7 @@ module plot_header integer :: basis ! direction of plot slice integer :: pixels(3) ! pixel width/height of plot slice integer :: meshlines_width ! pixel width of meshlines + integer :: level ! universe depth to plot the cells of type(StructuredMesh), pointer :: meshlines_mesh => null() ! mesh to plot type(ObjectColor) :: meshlines_color ! Color for meshlines type(ObjectColor) :: not_found ! color for positions where no cell found @@ -36,6 +37,9 @@ module plot_header ! Plot type integer, parameter :: PLOT_TYPE_SLICE = 1 integer, parameter :: PLOT_TYPE_VOXEL = 2 + + ! Plot level + integer, parameter :: PLOT_LEVEL_LOWEST = -1 ! Plot basis plane integer, parameter :: PLOT_BASIS_XY = 1 diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index b7975d05f3..62827f04c1 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -7,7 +7,7 @@ element geometry { (element material { ( xsd:int | "void" ) } | attribute material { ( xsd:int | "void" ) }) ) & - (element surfaces { list { xsd:int+ } } | attribute surfaces { list { xsd:int+ } }) & + (element surfaces { list { xsd:int* } } | attribute surfaces { list { xsd:int* } })? & (element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? & (element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })? }* @@ -23,8 +23,8 @@ element geometry { & element lattice { (element id { xsd:int } | attribute id { xsd:int }) & - (element type { ( "rectangular" | "hexagonal" ) } | - attribute type { ( "rectangular" | "hexagonal" ) })? & + (element type { ( "rect" | "rectangle" | "rectangular" | "hexagonal" ) } | + attribute type { ( "rect" | "rectangle" | "rectangular" | "hexagonal" ) })? & (element dimension { list { xsd:positiveInteger+ } } | attribute dimension { list { xsd:positiveInteger+ } }) & (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) & diff --git a/src/relaxng/plots.rnc b/src/relaxng/plots.rnc index 27b2ae7f72..5ef14f2297 100644 --- a/src/relaxng/plots.rnc +++ b/src/relaxng/plots.rnc @@ -7,6 +7,7 @@ element plots { attribute type { "slice" | "voxel" })? & (element color { ( "cell" | "mat" | "material" ) } | attribute color { ( "cell" | "mat" | "material" ) })? & + (element level { xsd:int } | attribute level { xsd:int })? & (element origin { list { xsd:double+ } } | attribute origin { list { xsd:double+ } })? & (element width { list { xsd:double+ } } | diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 707ecbbc80..a7b92e1738 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 07c1bfb411..3a12a3f718 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/tally.F90 b/src/tally.F90 index 87f5a2cab7..7c7747ff51 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 @@ -2021,10 +2013,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 +2026,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 +2043,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 +2056,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 +2073,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 +2086,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 +2210,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 23ab5b5634..b7dc5ed98a 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 d108325fa0..8eca3871c7 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 @@ -112,11 +111,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 ! ==================================================================== @@ -140,11 +138,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 diff --git a/tests/test_infinite_cell/geometry.xml b/tests/test_infinite_cell/geometry.xml new file mode 100644 index 0000000000..77ef6110a9 --- /dev/null +++ b/tests/test_infinite_cell/geometry.xml @@ -0,0 +1,17 @@ + + + + + + + + 11 12 + 12 11 + + + + + + + diff --git a/tests/test_infinite_cell/materials.xml b/tests/test_infinite_cell/materials.xml new file mode 100644 index 0000000000..2e5b48381a --- /dev/null +++ b/tests/test_infinite_cell/materials.xml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + diff --git a/tests/test_infinite_cell/results.py b/tests/test_infinite_cell/results.py new file mode 100644 index 0000000000..be13ee66f1 --- /dev/null +++ b/tests/test_infinite_cell/results.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python + +import sys + +# import statepoint +sys.path.insert(0, '../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_infinite_cell/results_true.dat b/tests/test_infinite_cell/results_true.dat new file mode 100644 index 0000000000..45eaa4f917 --- /dev/null +++ b/tests/test_infinite_cell/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +9.998895E-02 2.846817E-04 diff --git a/tests/test_infinite_cell/settings.xml b/tests/test_infinite_cell/settings.xml new file mode 100644 index 0000000000..a6fd5da19e --- /dev/null +++ b/tests/test_infinite_cell/settings.xml @@ -0,0 +1,16 @@ + + + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/test_infinite_cell/test_infinite_cell.py b/tests/test_infinite_cell/test_infinite_cell.py new file mode 100644 index 0000000000..6fdbf87459 --- /dev/null +++ b/tests/test_infinite_cell/test_infinite_cell.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + try: + test_run() + test_created_statepoint() + test_results() + finally: + teardown()