diff --git a/docs/source/publications.rst b/docs/source/publications.rst index 9d3491b8a..99a00d13f 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -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). ``_ + - 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). ``_ diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 4a92305bc..de3cbc7f2 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -260,6 +260,65 @@ or sub-elements and can be set to either "false" or "true". *Default*: true +```` 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 + + + + U-238 + ARES + 92238.72c + 92238.00c + 5.0e-6 + 40.0e-6 + + + Pu-239 + dbrc + 94239.72c + 94239.00c + 0.01e-6 + 210.0e-6 + + + + .. 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) + ```` Element ---------------------- diff --git a/src/ace.F90 b/src/ace.F90 index 3e9ae7394..fb42180da 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -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 diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 2067030c1..6d339891c 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -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) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 18a23a6b0..0620cf11f 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 20a26dce3..b136c2c16 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c42943e5c..12bd01b5e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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) diff --git a/src/physics.F90 b/src/physics.F90 index 47da60f2b..04cc66afa 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index ea3ecabf5..dd4b91dc7 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -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 })? + }* }? } diff --git a/tests/test_resonance_scattering/geometry.xml b/tests/test_resonance_scattering/geometry.xml new file mode 100644 index 000000000..612e46132 --- /dev/null +++ b/tests/test_resonance_scattering/geometry.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/tests/test_resonance_scattering/materials.xml b/tests/test_resonance_scattering/materials.xml new file mode 100644 index 000000000..52a8c04be --- /dev/null +++ b/tests/test_resonance_scattering/materials.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/test_resonance_scattering/results.py b/tests/test_resonance_scattering/results.py new file mode 100644 index 000000000..8ff10971c --- /dev/null +++ b/tests/test_resonance_scattering/results.py @@ -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) diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat new file mode 100644 index 000000000..62b0ac223 --- /dev/null +++ b/tests/test_resonance_scattering/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +6.452021E-02 1.738736E-03 diff --git a/tests/test_resonance_scattering/settings.xml b/tests/test_resonance_scattering/settings.xml new file mode 100644 index 000000000..7ce4f23ac --- /dev/null +++ b/tests/test_resonance_scattering/settings.xml @@ -0,0 +1,27 @@ + + + + + + U-238 + cxs + 92238.71c + 92238.71c + 5.0e-6 + 40.0e-6 + + + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py new file mode 100644 index 000000000..6fdbf8745 --- /dev/null +++ b/tests/test_resonance_scattering/test_resonance_scattering.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()