crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #256 from walshjon/res_scat

Browse source 
Improved resonance scattering methods
This commit is contained in:
Paul Romano 2014-09-10 20:03:28 -04:00
commit 55fb5a582b
15 changed files with 791 additions and 72 deletions
Download patch file Download diff file Expand all files Collapse all files

4
docs/source/publications.rst
View file

@ -4,6 +4,10 @@
Publications
============
- Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Accelerated sampling
of the free gas resonance elastic scattering kernel," *Ann. Nucl. Energy*,
**69**, 116--124 (2014). `<http://dx.doi.org/10.1016/j.anucene.2014.01.017>`_
- Benoit Forget, Sheng Xu, and Kord Smith, "Direct Doppler broadening in Monte
Carlo simulations using the multipole representation," *Ann. Nucl. Energy*,
**64**, 78--85 (2014). `<http://dx.doi.org/10.1016/j.anucene.2013.09.043>`_

59
docs/source/usersguide/input.rst
View file

@ -260,6 +260,65 @@ or sub-elements and can be set to either "false" or "true".
*Default*: true
``<resonance_scattering>`` Element
----------------------
The ``resonance_scattering`` element can contain one or more of the following
attributes or sub-elements:
:scatterer:
An element with attributes/sub-elements called ``nuclide``, ``method``,
``xs_label``, ``xs_label_0K``, ``E_min``, and ``E_max``. The ``nuclide``
attribute is the name, as given by the ``name`` attribute within the
``nuclide`` sub-element of the ``material`` element in ``materials.xml``,
of the nuclide to which a resonance scattering treatment is to be applied.
The ``method`` attribute gives the type of resonance scattering treatment
that is to be applied to the ``nuclide``. Acceptable inputs - none of
which are case-sensitive - for the ``method`` attribute are ``ARES``,
``CXS``, ``WCM``, and ``DBRC``. Descriptions of each of these methods
are documented here_. The ``xs_label`` attribute gives the label for the
cross section data of the ``nuclide`` at a given temperature. The
``xs_label_0K`` gives the label for the 0 K cross section data for the
``nuclide``. The ``E_min`` attribute gives the minimum energy above
which the ``method`` is applied. The ``E_max`` attribute gives the
maximum energy below which the ``method`` is applied. One example would
be as follows:
.. _here: http://dx.doi.org/10.1016/j.anucene.2014.01.017
.. code-block:: xml
<resonance_scattering>
<scatterer>
<nuclide>U-238</nuclide>
<method>ARES</method>
<xs_label>92238.72c</xs_label>
<xs_label_0K>92238.00c</xs_label_0K>
<E_min>5.0e-6</E_min>
<E_max>40.0e-6</E_max>
</scatterer>
<scatterer>
<nuclide>Pu-239</nuclide>
<method>dbrc</method>
<xs_label>94239.72c</xs_label>
<xs_label_0K>94239.00c</xs_label_0K>
<E_min>0.01e-6</E_min>
<E_max>210.0e-6</E_max>
</scatterer>
</resonance_scattering>
.. note:: If the ``resonance_scattering`` element is not given, the free gas,
constant cross section (``cxs``) scattering model, which has
historically been used by Monte Carlo codes to sample target
velocities, is used to treat the target motion of all nuclides. If
``resonance_scattering`` is present, the ``cxs`` method is applied
below ``E_min`` and the target-at-rest (asymptotic) kernel is used
above ``E_max``. An arbitrary number of ``scatterer`` elements may
be specified, each corresponding to a single nuclide at a single
temperature.
*Defaults*: None (scatterer), ARES (method), 0.01 eV (E_min), 1.0 keV (E_max)
``<run_cmfd>`` Element
----------------------

151
src/ace.F90
View file

@ -36,6 +36,7 @@ contains
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: k ! index over S(a,b) tables in material
integer :: n ! index over resonant scatterers
integer :: i_listing ! index in xs_listings array
integer :: i_nuclide ! index in nuclides
integer :: i_sab ! index in sab_tables
@ -80,6 +81,30 @@ contains
! array
call read_ace_table(i_nuclide, i_listing)
! 0K resonant scatterer information, if treating resonance scattering
if (treat_res_scat) then
do n = 1, n_res_scatterers_total
if (name == nuclides_0K(n) % name) then
nuclides(i_nuclide) % resonant = .true.
nuclides(i_nuclide) % name_0K = nuclides_0K(n) % name_0K
nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % &
& name_0K)
nuclides(i_nuclide) % scheme = nuclides_0K(n) % scheme
nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % &
& scheme)
nuclides(i_nuclide) % E_min = nuclides_0K(n) % E_min
nuclides(i_nuclide) % E_max = nuclides_0K(n) % E_max
if (.not. already_read % contains(nuclides(i_nuclide) % &
& name_0K)) then
i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % &
& name_0K)
call read_ace_table(i_nuclide, i_listing)
end if
exit
end if
end do
end if
! Add name and alias to dictionary
call already_read % add(name)
call already_read % add(alias)
@ -217,6 +242,7 @@ contains
real(8) :: awrs(16) ! list of atomic weight ratios (not used)
real(8) :: awr ! atomic weight ratio for table
logical :: file_exists ! does ACE library exist?
logical :: data_0K ! are we reading 0K data?
character(7) :: readable ! is ACE library readable?
character(10) :: name ! name of ACE table
character(10) :: date_ ! date ACE library was processed
@ -320,19 +346,32 @@ contains
select case(listing % type)
case (ACE_NEUTRON)
! only read in a resonant scatterers info once
nuc => nuclides(i_table)
nuc % name = name
nuc % awr = awr
nuc % kT = kT
nuc % zaid = NXS(2)
data_0K = .false.
if (trim(adjustl(name)) == nuc % name_0K) then
data_0K = .true.
else
nuc % name = name
nuc % awr = awr
nuc % kT = kT
nuc % zaid = NXS(2)
end if
! read all blocks
call read_esz(nuc)
call read_nu_data(nuc)
call read_reactions(nuc)
call read_angular_dist(nuc)
call read_energy_dist(nuc)
call read_unr_res(nuc)
call read_esz(nuc, data_0K)
! don't read unnecessary 0K data for resonant scatterers
if (data_0K) then
continue
else
call read_nu_data(nuc)
call read_reactions(nuc)
call read_angular_dist(nuc)
call read_energy_dist(nuc)
call read_unr_res(nuc)
end if
! Currently subcritical fixed source calculations are not allowed. Thus,
! if any fissionable material is found in a fixed source calculation,
@ -344,9 +383,12 @@ contains
! for fissionable nuclides, precalculate microscopic nu-fission cross
! sections so that we don't need to call the nu_total function during
! cross section lookups
! cross section lookups (except if we're dealing w/ 0K data for resonant
! scatterers)
if (nuc % fissionable) call generate_nu_fission(nuc)
if (nuc % fissionable .and. .not. data_0K) then
call generate_nu_fission(nuc)
end if
case (ACE_THERMAL)
sab => sab_tables(i_table)
@ -377,43 +419,82 @@ contains
! total xs, absorption xs, elastic scattering xs, and heating numbers.
!===============================================================================
subroutine read_esz(nuc)
subroutine read_esz(nuc, data_0K)
type(Nuclide), pointer :: nuc
logical :: data_0K ! are we reading 0K data?
integer :: NE ! number of energy points for total and elastic cross sections
integer :: i ! index in 0K elastic xs array for this nuclide
real(8) :: xs_cdf_sum = ZERO ! xs cdf value
! determine number of energy points
NE = NXS(3)
nuc % n_grid = NE
! allocate storage for energy grid and cross section arrays
allocate(nuc % energy(NE))
allocate(nuc % total(NE))
allocate(nuc % elastic(NE))
allocate(nuc % fission(NE))
allocate(nuc % nu_fission(NE))
allocate(nuc % absorption(NE))
! initialize cross sections
nuc % total = ZERO
nuc % elastic = ZERO
nuc % fission = ZERO
nuc % nu_fission = ZERO
nuc % absorption = ZERO
! read in 0K data if we've already read in non-0K data
if (data_0K) then
nuc % n_grid_0K = NE
allocate(nuc % energy_0K(NE))
allocate(nuc % elastic_0K(NE))
allocate(nuc % xs_cdf(NE))
nuc % elastic_0K = ZERO
nuc % xs_cdf = ZERO
XSS_index = 1
nuc % energy_0K = get_real(NE)
! Read data from XSS -- only the energy grid, elastic scattering and heating
! cross section values are actually read from here. The total and absorption
! cross sections are reconstructed from the partial reaction data.
! Skip total and absorption
XSS_index = XSS_index + 2*NE
! Continue reading elastic scattering and heating
nuc % elastic_0K = get_real(NE)
XSS_index = 1
nuc % energy = get_real(NE)
do i = 1, nuc % n_grid_0K - 1
! Skip total and absorption
XSS_index = XSS_index + 2*NE
! Negative cross sections result in a CDF that is not monotonically
! increasing. Set all negative xs values to ZERO.
if (nuc % elastic_0K(i) < ZERO) nuc % elastic_0K(i) = ZERO
! Continue reading elastic scattering and heating
nuc % elastic = get_real(NE)
! build xs cdf
xs_cdf_sum = xs_cdf_sum + (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) &
& + sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO &
& * (nuc % energy_0K(i+1) - nuc % energy_0K(i))
nuc % xs_cdf(i) = xs_cdf_sum
end do
else ! read in non-0K data
nuc % n_grid = NE
allocate(nuc % energy(NE))
allocate(nuc % total(NE))
allocate(nuc % elastic(NE))
allocate(nuc % fission(NE))
allocate(nuc % nu_fission(NE))
allocate(nuc % absorption(NE))
! initialize cross sections
nuc % total = ZERO
nuc % elastic = ZERO
nuc % fission = ZERO
nuc % nu_fission = ZERO
nuc % absorption = ZERO
! Read data from XSS -- only the energy grid, elastic scattering and heating
! cross section values are actually read from here. The total and absorption
! cross sections are reconstructed from the partial reaction data.
XSS_index = 1
nuc % energy = get_real(NE)
! Skip total and absorption
XSS_index = XSS_index + 2*NE
! Continue reading elastic scattering and heating
nuc % elastic = get_real(NE)
end if
end subroutine read_esz

45
src/ace_header.F90
View file

@ -112,6 +112,17 @@ module ace_header
real(8), allocatable :: absorption(:) ! absorption (MT > 100)
real(8), allocatable :: heating(:) ! heating
! Resonance scattering info
logical :: resonant = .false. ! resonant scatterer?
character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c
character(16) :: scheme ! target velocity sampling scheme
integer :: n_grid_0K ! number of 0K energy grid points
real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs
real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section
real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section
real(8) :: E_min ! lower cutoff energy for res scattering
real(8) :: E_max ! upper cutoff energy for res scattering
! Fission information
logical :: fissionable ! nuclide is fissionable?
logical :: has_partial_fission ! nuclide has partial fission reactions?
@ -147,6 +158,22 @@ module ace_header
procedure :: clear => nuclide_clear ! Deallocates Nuclide
end type Nuclide
!===============================================================================
! NUCLIDE0K temporarily contains all 0K cross section data and other parameters
! needed to treat resonance scattering before transferring them to NUCLIDE
!===============================================================================
type Nuclide0K
character(10) :: nuclide ! name of nuclide, e.g. U-238
character(16) :: scheme = 'ares' ! target velocity sampling scheme
character(10) :: name ! name of nuclide, e.g. 92235.03c
character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c
real(8) :: E_min = 0.01e-6 ! lower cutoff energy for res scattering
real(8) :: E_max = 1000.0e-6 ! upper cutoff energy for res scattering
end type Nuclide0K
!===============================================================================
! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
! thermal scattering collisions which utilize a continuous secondary energy
@ -342,12 +369,24 @@ module ace_header
integer :: i ! Loop counter
if (allocated(this % grid_index)) deallocate(this % grid_index)
if (allocated(this % grid_index)) &
deallocate(this % grid_index)
if (allocated(this % energy)) &
deallocate(this % energy, this % total, this % elastic, &
this % fission, this % nu_fission, this % absorption)
if (allocated(this % heating)) deallocate(this % heating)
& this % fission, this % nu_fission, this % absorption)
if (allocated(this % energy_0K)) &
deallocate(this % energy_0K)
if (allocated(this % elastic_0K)) &
deallocate(this % elastic_0K)
if (allocated(this % xs_cdf)) &
deallocate(this % xs_cdf)
if (allocated(this % heating)) &
deallocate(this % heating)
if (allocated(this % index_fission)) deallocate(this % index_fission)

37
src/cross_section.F90
View file

@ -513,4 +513,41 @@ contains
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
!===============================================================================
function elastic_xs_0K(E, nuc) result(xs_out)
type(Nuclide), pointer :: nuc ! target nuclide at temperature
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
real(8), intent(inout) :: E ! trial energy
real(8) :: xs_out ! 0K xs at trial energy
! Determine index on nuclide energy grid
if (E < nuc % energy_0K(1)) then
i_grid = 1
elseif (E > nuc % energy_0K(nuc % n_grid_0K)) then
i_grid = nuc % n_grid_0K - 1
else
i_grid = binary_search(nuc % energy_0K, nuc % n_grid_0K, E)
end if
! check for rare case where two energy points are the same
if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then
i_grid = i_grid + 1
end if
! calculate interpolation factor
f = (E - nuc % energy_0K(i_grid)) &
& / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid))
! Calculate microscopic nuclide elastic cross section
xs_out = (ONE - f) * nuc % elastic_0K(i_grid) &
& + f * nuc % elastic_0K(i_grid + 1)
end function elastic_xs_0K
end module cross_section

14
src/global.F90
View file

@ -1,7 +1,7 @@
module global
use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, &
MaterialMacroXS
MaterialMacroXS, Nuclide0K
use bank_header, only: Bank
use cmfd_header
use constants
@ -384,6 +384,13 @@ module global
logical :: output_xs = .false.
logical :: output_tallies = .true.
! ============================================================================
! RESONANCE SCATTERING VARIABLES
logical :: treat_res_scat = .false. ! is resonance scattering treated?
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, message, &
!$omp& trace, thread_id, current_work, matching_bins)
@ -417,6 +424,11 @@ contains
end do
deallocate(nuclides)
end if
if (allocated(nuclides_0K)) then
deallocate(nuclides_0K)
end if
if (allocated(sab_tables)) deallocate(sab_tables)
if (allocated(xs_listings)) deallocate(xs_listings)
if (allocated(micro_xs)) deallocate(micro_xs)

114
src/input_xml.F90
View file

@ -61,16 +61,19 @@ contains
character(MAX_FILE_LEN) :: env_variable
character(MAX_WORD_LEN) :: type
character(MAX_LINE_LEN) :: filename
type(Node), pointer :: doc => null()
type(Node), pointer :: node_mode => null()
type(Node), pointer :: node_source => null()
type(Node), pointer :: node_dist => null()
type(Node), pointer :: node_cutoff => null()
type(Node), pointer :: node_entropy => null()
type(Node), pointer :: node_ufs => null()
type(Node), pointer :: node_sp => null()
type(Node), pointer :: node_output => null()
type(Node), pointer :: node_verb => null()
type(Node), pointer :: doc => null()
type(Node), pointer :: node_mode => null()
type(Node), pointer :: node_source => null()
type(Node), pointer :: node_dist => null()
type(Node), pointer :: node_cutoff => null()
type(Node), pointer :: node_entropy => null()
type(Node), pointer :: node_ufs => null()
type(Node), pointer :: node_sp => null()
type(Node), pointer :: node_output => null()
type(Node), pointer :: node_verb => null()
type(Node), pointer :: node_res_scat => null()
type(Node), pointer :: node_scatterer => null()
type(NodeList), pointer :: node_scat_list => null()
! Display output message
message = "Reading settings XML file..."
@ -789,6 +792,88 @@ contains
end if
end if
! Resonance scattering parameters
if (check_for_node(doc, "resonance_scattering")) then
call get_node_ptr(doc, "resonance_scattering", node_res_scat)
call get_node_list(node_res_scat, "scatterer", node_scat_list)
! check that a nuclide is specified
if (get_list_size(node_scat_list) >= 1) then
treat_res_scat = .true.
n_res_scatterers_total = get_list_size(node_scat_list)
! store 0K info for resonant scatterers
allocate(nuclides_0K(n_res_scatterers_total))
do i = 1, n_res_scatterers_total
call get_list_item(node_scat_list, i, node_scatterer)
! check to make sure a nuclide is specified
if (.not. check_for_node(node_scatterer, "nuclide")) then
message = "No nuclide specified for scatterer " // trim(to_str(i)) &
// " in settings.xml file!"
call fatal_error()
end if
call get_node_value(node_scatterer, "nuclide", &
nuclides_0K(i) % nuclide)
if (check_for_node(node_scatterer, "method")) then
call get_node_value(node_scatterer, "method", &
nuclides_0K(i) % scheme)
end if
! check to make sure xs name for which method is applied is given
if (.not. check_for_node(node_scatterer, "xs_label")) then
message = "Must specify the temperature dependent name of " // '' &
//"scatterer " // trim(to_str(i)) // " given in cross_sections.xml"
call fatal_error()
end if
call get_node_value(node_scatterer, "xs_label", &
nuclides_0K(i) % name)
! check to make sure 0K xs name for which method is applied is given
if (.not. check_for_node(node_scatterer, "xs_label_0K")) then
message = "Must specify the 0K name of " // '' &
//"scatterer "// trim(to_str(i)) // " given in cross_sections.xml"
call fatal_error()
end if
call get_node_value(node_scatterer, "xs_label_0K", &
nuclides_0K(i) % name_0K)
if (check_for_node(node_scatterer, "E_min")) then
call get_node_value(node_scatterer, "E_min", &
nuclides_0K(i) % E_min)
end if
! check that E_min is non-negative
if (nuclides_0K(i) % E_min < ZERO) then
message = "Lower resonance scattering energy bound is negative"
call fatal_error()
end if
if (check_for_node(node_scatterer, "E_max")) then
call get_node_value(node_scatterer, "E_max", &
nuclides_0K(i) % E_max)
end if
! check that E_max is not less than E_min
if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then
message = "Lower resonance scattering energy bound exceeds upper"
call fatal_error()
end if
nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide)
nuclides_0K(i) % scheme = trim(nuclides_0K(i) % scheme)
call lower_case(nuclides_0K(i) % scheme)
nuclides_0K(i) % name = trim(nuclides_0K(i) % name)
nuclides_0K(i) % name_0K = trim(nuclides_0K(i) % name_0K)
end do
else
message = "No resonant scatterers are specified within the " // "" &
// "resonance_scattering element in settings.xml"
call fatal_error()
end if
end if
! Natural element expansion option
if (check_for_node(doc, "natural_elements")) then
call get_node_value(doc, "natural_elements", temp_str)
@ -3190,6 +3275,15 @@ contains
end if
end do
! Check that 0K nuclides are listed in the cross_sections.xml file
do i = 1, n_res_scatterers_total
if (.not. xs_listing_dict % has_key(trim(nuclides_0K(i) % name_0K))) then
message = "Could not find nuclide " // trim(nuclides_0K(i) % name_0K) // &
" in cross_sections.xml file!"
call fatal_error()
end if
end do
! Close cross sections XML file
call close_xmldoc(doc)

292
src/physics.F90
View file

@ -2,6 +2,7 @@ module physics
use ace_header, only: Nuclide, Reaction, DistEnergy
use constants
use cross_section, only: elastic_xs_0K
use endf, only: reaction_name
use error, only: fatal_error, warning
use fission, only: nu_total, nu_delayed
@ -349,7 +350,7 @@ contains
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, rxn, &
p % E, p % coord0 % uvw, p % mu)
p % E, p % coord0 % uvw, p % mu, p % wgt)
end if
p % event_MT = ELASTIC
@ -408,13 +409,14 @@ contains
! target.
!===============================================================================
subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab)
subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt)
integer, intent(in) :: i_nuclide
type(Reaction), pointer :: rxn
real(8), intent(inout) :: E
real(8), intent(inout) :: uvw(3)
real(8), intent(out) :: mu_lab
real(8), intent(inout) :: wgt
real(8) :: awr ! atomic weight ratio of target
real(8) :: mu_cm ! cosine of polar angle in center-of-mass
@ -437,7 +439,8 @@ contains
! Sample velocity of target nucleus
if (.not. micro_xs(i_nuclide) % use_ptable) then
call sample_target_velocity(nuc, v_t, E, uvw)
call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, &
& micro_xs(i_nuclide) % elastic)
else
v_t = ZERO
end if
@ -729,19 +732,270 @@ contains
end subroutine sab_scatter
!===============================================================================
! SAMPLE_TARGET_VELOCITY samples the target velocity based on the free gas
! scattering formulation used by most Monte Carlo codes. Excellent documentation
! for this method can be found in FRA-TM-123.
! SAMPLE_TARGET_VELOCITY samples the target velocity. The constant cross section
! free gas model is the default method. Methods for correctly accounting
! for the energy dependence of cross sections in treating resonance elastic
! scattering such as the DBRC, WCM, and a new, accelerated scheme are also
! implemented here.
!===============================================================================
subroutine sample_target_velocity(nuc, v_target, E, uvw)
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff)
type(Nuclide), pointer :: nuc
type(Nuclide), pointer :: nuc ! target nuclide at temperature T
real(8), intent(out) :: v_target(3) ! target velocity
real(8), intent(in) :: v_neut(3) ! neutron velocity
real(8), intent(in) :: E ! particle energy
real(8), intent(in) :: uvw(3) ! direction cosines
real(8), intent(inout) :: wgt ! particle weight
real(8) :: awr ! target/neutron mass ratio
real(8) :: kT ! equilibrium temperature of target in MeV
real(8) :: E_rel ! trial relative energy
real(8) :: xs_0K ! 0K xs at E_rel
real(8) :: xs_eff ! effective elastic xs at temperature T
real(8) :: wcf ! weight correction factor
real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1)
real(8) :: E_low ! lowest practical relative energy
real(8) :: E_up ! highest practical relative energy
real(8) :: E_mode ! most probable Maxwellian energy
real(8) :: E_t_max ! highest practical target energy
real(8) :: E_t ! trial target energy
real(8) :: xs_max ! max 0K xs over practical relative energies
real(8) :: xs_low ! 0K xs at lowest practical relative energy
real(8) :: xs_up ! 0K xs at highest practical relative energy
real(8) :: m ! slope for interpolation
real(8) :: R_dbrc ! DBRC rejection criterion
real(8) :: R_speed ! target speed rejection criterion
real(8) :: cdf_low ! xs cdf at lowest practical relative energy
real(8) :: cdf_up ! xs cdf at highest practical relative energy
real(8) :: cdf_rel ! trial xs cdf value
real(8) :: p_mode ! probability at most probable energy
real(8) :: p_t ! probability at trial target energy
real(8) :: mu ! cosine between neutron and target velocities
integer :: i_E_low ! 0K index to lowest practical relative energy
integer :: i_E_up ! 0K index to highest practical relative energy
integer :: i_E_rel ! index to trial relative energy
logical :: reject ! resample if true
character(80) :: sampling_scheme ! method of target velocity sampling
kT = nuc % kT
awr = nuc % awr
! check if nuclide is a resonant scatterer
if (nuc % resonant) then
! sampling scheme to use
sampling_scheme = nuc % scheme
! upper resonance scattering energy bound (target is at rest above this E)
if (E > nuc % E_max) then
v_target = ZERO
return
! lower resonance scattering energy bound (should be no resonances below)
else if (E < nuc % E_min) then
sampling_scheme = 'cxs'
end if
! otherwise, use free gas model
else
if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then
v_target = ZERO
return
else
sampling_scheme = 'cxs'
end if
end if
! use appropriate target velocity sampling method
select case (sampling_scheme)
case ('cxs')
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw)
case ('wcm')
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw)
! adjust weight as prescribed by the weight correction method (wcm)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
xs_0K = elastic_xs_0K(E_rel, nuc)
wcf = xs_0K / xs_eff
wgt = wcf * wgt
case ('dbrc')
E_red = sqrt((awr * E) / kT)
E_low = (((E_red - 4.0_8)**2) * kT) / awr
E_up = (((E_red + 4.0_8)**2) * kT) / awr
! find lower and upper energy bound indices
! lower index
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_low = nuc % n_grid_0K - 1
else
i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low)
end if
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_up = nuc % n_grid_0K - 1
else
i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up)
end if
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
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)
reject = .true.
! sample target velocities until one is accepted by the DBRC
do
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw)
! perform Doppler broadening rejection correction (dbrc)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
xs_0K = elastic_xs_0K(E_rel, nuc)
R_dbrc = xs_0K / xs_max
if (prn() < R_dbrc) reject = .false.
if (.not. reject) exit
end do
case ('ares')
E_red = sqrt((awr * E) / kT)
E_low = (((E_red - 4.0_8)**2) * kT) / awr
E_up = (((E_red + 4.0_8)**2) * kT) / awr
! find lower and upper energy bound indices
! lower index
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_low = nuc % n_grid_0K - 1
else
i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low)
end if
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then
i_E_up = nuc % n_grid_0K - 1
else
i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up)
end if
! interpolate xs CDF since we're not exactly at the energy indices
! cdf value at lower bound attainable energy
if (i_E_low > 1) then
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
& + m * (E_low - nuc % energy_0K(i_E_low))
else
m = nuc % xs_cdf(i_E_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
end if
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
& + m * (E_up - nuc % energy_0K(i_E_up))
! values used to sample the Maxwellian
E_mode = kT
p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_mode / kT)
E_t_max = 16.0_8 * E_mode
reject = .true.
do
! perform Maxwellian rejection sampling
E_t = E_t_max * prn()**2
p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_t / kT)
R_speed = p_t / p_mode
if (prn() < R_speed) then
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
& i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
& - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
& / (nuc % energy_0K(i_E_low + i_E_rel) &
& - nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
! perform rejection sampling on cosine between
! neutron and target velocities
mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t))
if (abs(mu) < ONE) then
! set and accept target velocity
E_t = E_t / awr
v_target = sqrt(E_t) * rotate_angle(uvw, mu)
reject = .false.
end if
end if
if (.not. reject) exit
end do
case default
message = "Not a recognized resonance scattering treatment!"
call fatal_error()
end select
end subroutine sample_target_velocity
!===============================================================================
! SAMPLE_CXS_TARGET_VELOCITY samples a target velocity based on the free gas
! scattering formulation, used by most Monte Carlo codes, in which cross section
! is assumed to be constant in energy. Excellent documentation for this method
! can be found in FRA-TM-123.
!===============================================================================
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw)
type(Nuclide), pointer :: nuc ! target nuclide at temperature
real(8), intent(out) :: v_target(3)
real(8), intent(in) :: E
real(8), intent(in) :: uvw(3)
real(8) :: kT ! equilibrium temperature of target in MeV
real(8) :: awr ! target/neutron mass ratio
real(8) :: alpha ! probability of sampling f2 over f1
real(8) :: mu ! cosine of angle between neutron and target vel
real(8) :: r1, r2 ! pseudo-random numbers
@ -752,18 +1006,10 @@ contains
real(8) :: beta_vt_sq ! (beta * speed of target)^2
real(8) :: vt ! speed of target
! Determine equilibrium temperature in MeV
kT = nuc % kT
awr = nuc % awr
! Check if energy is above threshold
if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then
v_target = ZERO
return
end if
! calculate beta
beta_vn = sqrt(nuc%awr * E / kT)
beta_vn = sqrt(awr * E / kT)
alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO)
do
@ -795,20 +1041,20 @@ contains
! Determine rejection probability
accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) &
/(beta_vn + beta_vt)
/(beta_vn + beta_vt)
! Perform rejection sampling on vt and mu
if (prn() < accept_prob) exit
end do
! determine speed of target nucleus
vt = sqrt(beta_vt_sq*kT/nuc % awr)
! Determine speed of target nucleus
vt = sqrt(beta_vt_sq*kT/awr)
! determine velocity vector of target nucleus based on neutron's velocity
! Determine velocity vector of target nucleus based on neutron's velocity
! and the sampled angle between them
v_target = vt * rotate_angle(uvw, mu)
end subroutine sample_target_velocity
end subroutine sample_cxs_target_velocity
!===============================================================================
! CREATE_FISSION_SITES determines the average total, prompt, and delayed

17
src/relaxng/settings.rnc
View file

@ -129,5 +129,22 @@ element settings {
attribute lower_left { list { xsd:double+ } }) &
(element upper_right { list { xsd:double+ } } |
attribute upper_right { list { xsd:double+ } })
}? &
element resonance_scattering {
element scatterer {
(element nuclide { xsd:string { maxLength = "12" } } |
attribute nuclide { xsd:string { maxLength = "12" } }) &
(element method { xsd:string { maxLength = "16" } } |
attribute method { xsd:string { maxLength = "16" } }) &
(element xs_label { xsd:string { maxLength = "12" } } |
attribute xs_label { xsd:string { maxLength = "12" } }) &
(element xs_label_0K { xsd:string { maxLength = "12" } } |
attribute xs_label_0K { xsd:string { maxLength = "12" } }) &
(element E_min { xsd:double } |
attribute E_min { xsd:double }) &
(element E_max { xsd:double } |
attribute E_max { xsd:double })?
}*
}?
}

8
tests/test_resonance_scattering/geometry.xml Normal file
View file

@ -0,0 +1,8 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
</geometry>

9
tests/test_resonance_scattering/materials.xml Normal file
View file

@ -0,0 +1,9 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-238" xs="71c" ao="1.0" />
</material>
</materials>

25
tests/test_resonance_scattering/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../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)

2
tests/test_resonance_scattering/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
6.452021E-02 1.738736E-03

27
tests/test_resonance_scattering/settings.xml Normal file
View file

@ -0,0 +1,27 @@
<?xml version="1.0"?>
<settings>
<resonance_scattering>
<scatterer>
<nuclide>U-238</nuclide>
<method>cxs</method>
<xs_label>92238.71c</xs_label>
<xs_label_0K>92238.71c</xs_label_0K>
<E_min>5.0e-6</E_min>
<E_max>40.0e-6</E_max>
</scatterer>
</resonance_scattering>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

59
tests/test_resonance_scattering/test_resonance_scattering.py Normal file
View file

@ -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()