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 @@
-
-
-
-
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 7132a67cdb..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
@@ -1022,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",
@@ -1572,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 177956529f..1672405e8d 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/global.F90 b/src/global.F90
index 874b8d0d1a..ec11216608 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -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.
@@ -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
@@ -462,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)
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 69677d1757..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
@@ -2330,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)
@@ -2821,7 +2833,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)))
@@ -2927,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
@@ -2935,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)
@@ -2943,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", &
@@ -2955,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
@@ -2978,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')
@@ -2994,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')
@@ -3039,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 2a7d6890e1..26a4a317d7 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
@@ -1500,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 38cbeb5ef3..9b32c64d92 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -1832,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
@@ -1849,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 ddc4641b2a..d9f6b55e2b 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 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 7c7747ff51..c96f1ecd6e 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -1837,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