crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

merged develop into cmfd-additions

Browse source
This commit is contained in:
Bryan Herman 2014-09-16 16:20:20 -04:00
commit d692c05145
24 changed files with 958 additions and 166 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>`_

65
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
----------------------
@ -311,6 +370,12 @@ attributes/sub-elements:
parallelepiped and the last three of which specify the upper-right
corner. Source sites are sampled uniformly through that parallelepiped.
To filter a "box" spatial distribution by fissionable material, specify
"fission" tag instead of "box". The ``parameters`` should be given as six
real numbers, the first three of which specify the lower-left corner of a
parallelepiped and the last three of which specify the upper-right
corner. Source sites are sampled uniformly through that parallelepiped.
For a "point" spatial distribution, ``parameters`` should be given as
three real numbers which specify the (x,y,z) location of an isotropic
point source

2
src/CMakeLists.txt
View file

@ -287,7 +287,7 @@ file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_*.py)
# Check to see if PETSC is compiled for CMFD tests
if (NOT ${PETSC_ENABLED})
file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd*.py)
file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd_jfnk.py)
foreach(cmfd_test in ${CMFD_TESTS})
list(REMOVE_ITEM TESTS ${cmfd_test})
endforeach(cmfd_test)

179
src/ace.F90
View file

@ -7,7 +7,7 @@ module ace
use error, only: fatal_error, warning
use fission, only: nu_total
use global
use list_header, only: ListElemInt, ListInt
use list_header, only: ListInt
use material_header, only: Material
use output, only: write_message
use set_header, only: SetChar
@ -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)
@ -172,6 +197,25 @@ contains
! Avoid some valgrind leak errors
call already_read % clear()
! Loop around material
MATERIAL_LOOP3: do i = 1, n_materials
! Get material
mat => materials(i)
! Loop around nuclides in material
NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
! Check for fission in nuclide
if (nuclides(mat % nuclide(j)) % fissionable) then
mat % fissionable = .true.
exit NUCLIDE_LOOP2
end if
end do NUCLIDE_LOOP2
end do MATERIAL_LOOP3
end subroutine read_xs
!===============================================================================
@ -198,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
@ -301,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,
@ -325,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)
@ -358,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
@ -1463,16 +1563,11 @@ contains
integer :: i ! index in nuclides array
integer :: j ! index in nuclides array
type(ListElemInt), pointer :: nuc_list => null() ! pointer to nuclide list
do i = 1, n_nuclides_total
allocate(nuclides(i) % nuc_list)
nuc_list => nuclides(i) % nuc_list
do j = 1, n_nuclides_total
if (nuclides(i) % zaid == nuclides(j) % zaid) then
nuc_list % data = j
allocate(nuc_list % next)
nuc_list => nuc_list % next
call nuclides(i) % nuc_list % append(j)
end if
end do
end do

51
src/ace_header.F90
View file

@ -2,7 +2,7 @@ module ace_header
use constants, only: MAX_FILE_LEN
use endf_header, only: Tab1
use list_header, only: ListElemInt
use list_header, only: ListInt
implicit none
@ -97,7 +97,7 @@ module ace_header
real(8) :: kT ! temperature in MeV (k*T)
! Linked list of indices in nuclides array of instances of this same nuclide
type(ListElemInt), pointer :: nuc_list => null()
type(ListInt) :: nuc_list
! Energy grid information
integer :: n_grid ! # of nuclide grid points
@ -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)
@ -377,6 +416,8 @@ module ace_header
deallocate(this % reactions)
end if
call this % nuc_list % clear()
end subroutine nuclide_clear
end module ace_header

34
src/cmfd_input.F90
View file

@ -155,7 +155,7 @@ contains
call get_node_value(doc, "feedback", temp_str)
call lower_case(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_feedback = .true.
cmfd_feedback = .true.
end if
! Set downscatter logical
@ -163,7 +163,7 @@ contains
call get_node_value(doc, "downscatter", temp_str)
call lower_case(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_downscatter = .true.
cmfd_downscatter = .true.
end if
! Reset dhat parameters
@ -184,26 +184,26 @@ contains
! Set the solver type
if (check_for_node(doc, "solver")) &
call get_node_value(doc, "solver", cmfd_solver_type)
call get_node_value(doc, "solver", cmfd_solver_type)
! Set monitoring
if (check_for_node(doc, "snes_monitor")) then
call get_node_value(doc, "snes_monitor", temp_str)
call lower_case(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_snes_monitor = .true.
cmfd_snes_monitor = .true.
end if
if (check_for_node(doc, "ksp_monitor")) then
call get_node_value(doc, "ksp_monitor", temp_str)
call lower_case(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_ksp_monitor = .true.
cmfd_ksp_monitor = .true.
end if
if (check_for_node(doc, "power_monitor")) then
call get_node_value(doc, "power_monitor", temp_str)
call lower_case(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_power_monitor = .true.
cmfd_power_monitor = .true.
end if
! Output logicals
@ -211,7 +211,7 @@ contains
call get_node_value(doc, "write_matrices", temp_str)
call lower_case(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_write_matrices = .true.
cmfd_write_matrices = .true.
end if
! Run an adjoint calc
@ -228,7 +228,7 @@ contains
! Batch to begin cmfd
if (check_for_node(doc, "begin")) &
call get_node_value(doc, "begin", cmfd_begin)
call get_node_value(doc, "begin", cmfd_begin)
! Check for cmfd tally resets
if (check_for_node(doc, "tally_reset")) then
@ -247,9 +247,9 @@ contains
! Get display
if (check_for_node(doc, "display")) &
call get_node_value(doc, "display", cmfd_display)
call get_node_value(doc, "display", cmfd_display)
if (trim(cmfd_display) == 'dominance' .and. &
trim(cmfd_solver_type) /= 'power') then
trim(cmfd_solver_type) /= 'power') then
message = 'Dominance Ratio only aviable with power iteration solver'
call warning()
cmfd_display = ''
@ -257,17 +257,17 @@ contains
! Read in spectral radius estimate and tolerances
if (check_for_node(doc, "spectral")) &
call get_node_value(doc, "spectral", cmfd_spectral)
call get_node_value(doc, "spectral", cmfd_spectral)
if (check_for_node(doc, "shift")) &
call get_node_value(doc, "shift", cmfd_shift)
call get_node_value(doc, "shift", cmfd_shift)
if (check_for_node(doc, "ktol")) &
call get_node_value(doc, "ktol", cmfd_ktol)
call get_node_value(doc, "ktol", cmfd_ktol)
if (check_for_node(doc, "stol")) &
call get_node_value(doc, "stol", cmfd_stol)
call get_node_value(doc, "stol", cmfd_stol)
if (check_for_node(doc, "atoli")) &
call get_node_value(doc, "atoli", cmfd_atoli)
call get_node_value(doc, "atoli", cmfd_atoli)
if (check_for_node(doc, "rtoli")) &
call get_node_value(doc, "rtoli", cmfd_rtoli)
call get_node_value(doc, "rtoli", cmfd_rtoli)
! Create tally objects
call create_cmfd_tally(doc)
@ -440,7 +440,7 @@ contains
call get_node_value(doc, "reset", temp_str)
call lower_case(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
t % reset = .true.
t % reset = .true.
end if
! Set up mesh filter

77
src/cmfd_solver.F90
View file

@ -11,7 +11,6 @@ module cmfd_solver
private
public :: cmfd_solver_execute
logical :: iconv ! did the problem converged
real(8) :: k_n ! new k-eigenvalue
real(8) :: k_o ! old k-eigenvalue
real(8) :: k_s ! shift of eigenvalue
@ -38,11 +37,11 @@ module cmfd_solver
subroutine linsolve(A, b, x, tol, i)
import :: Matrix
import :: Vector
type(Matrix) :: A
type(Vector) :: b
type(Vector) :: x
real(8) :: tol
integer :: i
type(Matrix), intent(inout) :: A
type(Vector), intent(inout) :: b
type(Vector), intent(inout) :: x
real(8), intent(in) :: tol
integer, intent(out) :: i
end subroutine linsolve
end interface
@ -56,7 +55,7 @@ contains
use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve
logical, optional :: adjoint ! adjoint calc
logical, optional, intent(in) :: adjoint ! adjoint calc
logical :: physical_adjoint = .false.
@ -105,7 +104,7 @@ contains
use global, only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, &
cmfd_write_matrices
logical :: adjoint
logical, intent(in) :: adjoint
integer :: n ! problem size
real(8) :: guess ! initial guess
@ -144,7 +143,7 @@ contains
! Fill in production matrix
call build_prod_matrix(prod, adjoint=adjoint)
! Setup petsc for everything
! Finalize setup of CSR matrices
call loss % assemble()
call prod % assemble()
if (cmfd_write_matrices) then
@ -210,11 +209,13 @@ contains
subroutine execute_power_iter()
use constants, only: ONE
use global, only: cmfd_atoli, cmfd_rtoli
use error, only: fatal_error
use global, only: cmfd_atoli, cmfd_rtoli, message
integer :: i ! iteration counter
integer :: innerits ! # of inner iterations
integer :: totalits ! total number of inners
logical :: iconv ! did the problem converged
real(8) :: atoli ! absolute minimum tolerance
real(8) :: rtoli ! relative tolerance based on source conv
real(8) :: toli ! the current tolerance of inners
@ -234,6 +235,12 @@ contains
! Begin power iteration
do i = 1, 10000
! Check if reached iteration 10000
if (i == 10000) then
message = 'Reached maximum iterations in CMFD power iteration solver.'
call fatal_error()
end if
! Compute source vector
call prod % vector_multiply(phi_o, s_o)
@ -256,7 +263,7 @@ contains
s_o % val = s_o % val * k_lo
! Check convergence
call convergence(i, innerits)
call convergence(i, innerits, iconv)
totalits = totalits + innerits
! Break loop if converged
@ -305,14 +312,15 @@ contains
! CONVERGENCE checks the convergence of the CMFD problem
!===============================================================================
subroutine convergence(iter, innerits)
subroutine convergence(iter, innerits, iconv)
use constants, only: ONE, ZERO
use global, only: cmfd_power_monitor, master
use, intrinsic :: ISO_FORTRAN_ENV
integer :: iter ! outer iteration number
integer :: innerits ! inner iteration nubmer
integer, intent(in) :: iter ! outer iteration number
integer, intent(in) :: innerits ! inner iteration nubmer
logical, intent(out) :: iconv ! convergence logical
! Reset convergence flag
iconv = .false.
@ -350,11 +358,11 @@ contains
use constants, only: ONE, ZERO
use global, only: cmfd, cmfd_spectral
type(Matrix) :: A ! coefficient matrix
type(Vector) :: b ! right hand side vector
type(Vector) :: x ! unknown vector
real(8) :: tol ! tolerance on final error
integer :: its ! number of inner iterations
type(Matrix), intent(inout) :: A ! coefficient matrix
type(Vector), intent(inout) :: b ! right hand side vector
type(Vector), intent(inout) :: x ! unknown vector
real(8), intent(in) :: tol ! tolerance on final error
integer, intent(out) :: its ! number of inner iterations
integer :: g ! group index
integer :: i ! loop counter for x
@ -430,7 +438,6 @@ contains
! Check convergence
err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n)
its = igs
!print *, igs, err, w
if (err < tol) exit
! Calculation new overrelaxation parameter
@ -451,11 +458,11 @@ contains
use constants, only: ONE, ZERO
use global, only: cmfd, cmfd_spectral
type(Matrix) :: A ! coefficient matrix
type(Vector) :: b ! right hand side vector
type(Vector) :: x ! unknown vector
real(8) :: tol ! tolerance on final error
integer :: its ! number of inner iterations
type(Matrix), intent(inout) :: A ! coefficient matrix
type(Vector), intent(inout) :: b ! right hand side vector
type(Vector), intent(inout) :: x ! unknown vector
real(8), intent(in) :: tol ! tolerance on final error
integer, intent(out) :: its ! number of inner iterations
integer :: g ! group index
integer :: i ! loop counter for x
@ -570,8 +577,6 @@ contains
! Check convergence
err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n)
its = igs
!print *, igs, err, w
!read *
if (err < tol) exit
! Calculation new overrelaxation parameter
@ -744,15 +749,15 @@ contains
use global, only: cmfd, cmfd_coremap
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: irow ! iteration counter over row (0 reference)
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer, intent(out) :: i ! iteration counter for x
integer, intent(out) :: j ! iteration counter for y
integer, intent(out) :: k ! iteration counter for z
integer, intent(out) :: g ! iteration counter for groups
integer, intent(in) :: irow ! iteration counter over row (0 reference)
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
! Check for core map
if (cmfd_coremap) then

5
src/constants.F90
View file

@ -335,8 +335,9 @@ module constants
! Source spatial distribution types
integer, parameter :: &
SRC_SPACE_BOX = 1, & ! Source in a rectangular prism
SRC_SPACE_POINT = 2 ! Source at a single point
SRC_SPACE_BOX = 1, & ! Source in a rectangular prism
SRC_SPACE_POINT = 2, & ! Source at a single point
SRC_SPACE_FISSION = 3 ! Source in prism filtered by fissionable mats
! Source angular distribution types
integer, parameter :: &

72
src/cross_section.F90
View file

@ -342,21 +342,22 @@ contains
integer, intent(in) :: i_nuclide ! index into nuclides array
real(8), intent(in) :: E ! energy
integer :: i_energy ! index for energy
integer :: i_low ! band index at lower bounding energy
integer :: i_up ! band index at upper bounding energy
real(8) :: f ! interpolation factor
real(8) :: r ! pseudo-random number
real(8) :: elastic ! elastic cross section
real(8) :: capture ! (n,gamma) cross section
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
logical :: same_nuc ! do we know the xs for this nuclide at this energy?
integer :: i ! loop index
integer :: i_energy ! index for energy
integer :: i_low ! band index at lower bounding energy
integer :: i_up ! band index at upper bounding energy
integer :: same_nuc_idx ! index of same nuclide
real(8) :: f ! interpolation factor
real(8) :: r ! pseudo-random number
real(8) :: elastic ! elastic cross section
real(8) :: capture ! (n,gamma) cross section
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
logical :: same_nuc ! do we know the xs for this nuclide at this energy?
type(UrrData), pointer, save :: urr => null()
type(Nuclide), pointer, save :: nuc => null()
type(Reaction), pointer, save :: rxn => null()
type(ListElemInt), pointer :: nuc_list => null()
!$omp threadprivate(urr, nuc, rxn, nuc_list)
!$omp threadprivate(urr, nuc, rxn)
micro_xs(i_nuclide) % use_ptable = .true.
@ -381,18 +382,16 @@ contains
! this energy but a different temperature, use the original random number to
! preserve correlation of temperature in probability tables
same_nuc = .false.
nuc_list => nuc % nuc_list
do
if (E /= ZERO .and. E == micro_xs(nuc_list % data) % last_E) then
do i = 1, nuc % nuc_list % size()
if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % get_item(i)) % last_E) then
same_nuc = .true.
same_nuc_idx = i
exit
end if
nuc_list => nuc_list % next
if (.not. associated(nuc_list % next)) exit
end do
if (same_nuc) then
r = micro_xs(nuc_list % data) % last_prn
r = micro_xs(nuc % nuc_list % get_item(same_nuc_idx)) % last_prn
else
r = prn()
micro_xs(i_nuclide) % last_prn = r
@ -513,4 +512,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
@ -392,6 +392,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)
@ -425,6 +432,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)

1
src/initialize.F90
View file

@ -11,6 +11,7 @@ module initialize
use global
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict, read_plots_xml
use material_header, only: Material
use output, only: title, header, write_summary, print_version, &
print_usage, write_xs_summary, print_plot, &
write_message

117
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..."
@ -300,6 +303,9 @@ contains
case ('box')
external_source % type_space = SRC_SPACE_BOX
coeffs_reqd = 6
case ('fission')
external_source % type_space = SRC_SPACE_FISSION
coeffs_reqd = 6
case ('point')
external_source % type_space = SRC_SPACE_POINT
coeffs_reqd = 3
@ -788,6 +794,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)
@ -3189,6 +3277,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)

4
src/material_header.F90
View file

@ -21,6 +21,10 @@ module material_header
! Temporary names read during initialization
character(12), allocatable :: names(:) ! isotope names
character(12), allocatable :: sab_names(:) ! name of S(a,b) table
! Does this material contain fissionable nuclides?
logical :: fissionable = .false.
end type Material
end module material_header

18
src/matrix_header.F90
View file

@ -368,11 +368,11 @@ contains
subroutine matrix_search_indices(self, row, col, idx, found)
class(Matrix) :: self
integer :: row
integer :: col
integer :: idx
logical :: found
class(Matrix), intent(inout) :: self
integer, intent(in) :: row
integer, intent(in) :: col
integer, intent(out) :: idx
logical, intent(out) :: found
integer :: j
@ -396,8 +396,8 @@ contains
subroutine matrix_write(self, filename)
character(*) :: filename
class(Matrix) :: self
character(*), intent(in) :: filename
class(Matrix), intent(inout) :: self
integer :: unit_
integer :: i
@ -421,8 +421,8 @@ contains
subroutine matrix_copy(self, mattocopy)
class(Matrix) :: self
type(Matrix) :: mattocopy
class(Matrix), intent(inout) :: self
type(Matrix), intent(in) :: mattocopy
! Set n and nnz
self % n_count = mattocopy % n_count

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 })?
}*
}?
}

36
src/source.F90
View file

@ -146,6 +146,42 @@ contains
end do
call p % clear()
case (SRC_SPACE_FISSION)
! Repeat sampling source location until a good site has been found
found = .false.
do while (.not.found)
! Set particle defaults
call p % initialize()
! Coordinates sampled uniformly over a box
p_min = external_source % params_space(1:3)
p_max = external_source % params_space(4:6)
r = (/ (prn(), i = 1,3) /)
site % xyz = p_min + r*(p_max - p_min)
! Fill p with needed data
p % coord0 % xyz = site % xyz
p % coord0 % uvw = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
if (.not. found) then
num_resamples = num_resamples + 1
if (num_resamples == MAX_EXTSRC_RESAMPLES) then
message = "Maximum number of external source spatial resamples &
&reached!"
call fatal_error()
end if
cycle
end if
if (p % material == MATERIAL_VOID) then
found = .false.
cycle
end if
if (.not. materials(p % material) % fissionable) found = .false.
end do
call p % clear()
case (SRC_SPACE_POINT)
! Point source
site % xyz = external_source % params_space

6
src/vector_header.F90
View file

@ -28,7 +28,9 @@ module vector_header
# endif
end type Vector
#ifdef PETSC
integer :: petsc_err ! petsc error code
#endif
contains
@ -132,8 +134,8 @@ contains
subroutine vector_copy(self, vectocopy)
class(Vector), target :: self
type(Vector) :: vectocopy
class(Vector), target, intent(inout) :: self
type(Vector), intent(in) :: vectocopy
! Preallocate vector
if (.not.allocated(self % data)) allocate(self % data(vectocopy % n))

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