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

fixes per @smharpers comments

Browse source
This commit is contained in:
Adam Nelson 2016-02-05 05:37:50 -05:00
parent 989569ac5e
commit 86d5600eb2
13 changed files with 124 additions and 377 deletions
Download patch file Download diff file Expand all files Collapse all files

25
docs/source/usersguide/mgxs_library.rst
View file

@ -114,20 +114,25 @@ attributes/sub-elements required to describe the meta-data:
*Default*: "isotropic"
:num_azimuthal:
This element provides the number of equi-width bins that the azimuthal
angular domain is subdivided in the case of angle-dependent cross sections
(i.e., "angle" is passed to the ``representation`` element).
This element provides the number of equal width angular bins that the
azimuthal angular domain is subdivided in the case of angle-dependent
cross sections (i.e., "angle" is passed to the ``representation`` element).
Note that these bins are equal in azimuthal angle widths, not equal in the
cosine of the azimuthal angle widths.
*Default*: If ``representation`` is "angle", this must be provided. If
not, this parameter is not used.
*Default*: If ``representation`` is "angle", this must be provided. This
parameter is not used for other ``representation`` types.
:num_polar:
This element provides the number of equi-width bins that the polar angular
domain is subdivided in the case of angle-dependent cross sections
(i.e., "angle" is passed to the ``representation`` element).
This element provides the number of equal width angular bins that the
polar angular domain is subdivided in the case of angle-dependent
cross sections (i.e., "angle" is passed to the ``representation`` element).
Note that these bins are equal in polar angle widths, not equal in the
cosine of the polar angle widths.
*Default*: If ``representation`` is "angle", this must be provided. If
not, this parameter is not used.
*Default*: If ``representation`` is "angle", this must be provided. This
parameter is not used for other ``representation`` types.
:scatt_type:
This element provides the representation of the angular distribution

2
examples/python/pincell_multigroup/build-xml.py
View file

@ -72,7 +72,7 @@ mg_cross_sections_file.export_to_xml()
uo2_data = openmc.Macroscopic('UO2', '300K')
h2o_data = openmc.Macroscopic('LWTR', '300K')
# Instantiate some Materials and register the appropriate Nuclides
# Instantiate some Materials and register the appropriate Macroscopic objects
uo2 = openmc.Material(material_id=1, name='UO2 fuel')
uo2.set_density('macro', 1.0)
uo2.add_macroscopic(uo2_data)

2
openmc/macroscopic.py
View file

@ -8,7 +8,7 @@ if sys.version_info[0] >= 3:
class Macroscopic(object):
"""A nuclide that can be used in a material.
"""A Macroscopic object that can be used in a material.
Parameters
----------

4
src/global.F90
View file

@ -490,10 +490,6 @@ contains
end if
if (allocated(macro_xs)) then
! First call the clear routines
do i = 1, size(macro_xs)
call macro_xs(i) % obj % clear()
end do
deallocate(macro_xs)
end if

2
src/input_xml.F90
View file

@ -2071,7 +2071,7 @@ contains
call get_list_item(node_macro_list, 1, node_nuc)
! Check for empty name on nuclide
if (.not.check_for_node(node_nuc, "name")) then
if (.not. check_for_node(node_nuc, "name")) then
call fatal_error("No name specified on macroscopic data in material " &
// trim(to_str(mat % id)))
end if

2
src/macroxs.F90
View file

@ -140,7 +140,7 @@ contains
!===============================================================================
subroutine macroxs_sample_scatter(scatt, gin, gout, mu, wgt)
Class(ScattData_Base), intent(in) :: scatt ! Scattering Object to Use
class(ScattData_Base), intent(in) :: scatt ! Scattering Object to Use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle

213
src/macroxs_header.F90
View file

@ -19,10 +19,9 @@ module macroxs_header
integer :: order
! Type-Bound procedures
contains
procedure(macroxs_init_), deferred, pass :: init ! initializes object
procedure(macroxs_clear_), deferred, pass :: clear ! Deallocates object
procedure(macroxs_get_xs_), deferred, pass :: get_xs ! Return xs
contains
procedure(macroxs_init_), deferred, pass :: init ! initializes object
procedure(macroxs_get_xs_), deferred, pass :: get_xs ! Return xs
end type MacroXS_Base
abstract interface
@ -79,11 +78,10 @@ module macroxs_header
real(8), allocatable :: scattxs(:) ! scattering xs
real(8), allocatable :: chi(:,:) ! fission spectra
! Type-Bound procedures
contains
procedure, pass :: init => macroxs_iso_init ! inits object
procedure, pass :: clear => macroxs_iso_clear ! Deallocates object
procedure, pass :: get_xs => macroxs_iso_get_xs ! Returns xs
! Type-Bound procedures
contains
procedure, pass :: init => macroxs_iso_init ! inits object
procedure, pass :: get_xs => macroxs_iso_get_xs ! Returns xs
end type MacroXS_Iso
type, extends(MacroXS_Base) :: MacroXS_Angle
@ -99,11 +97,10 @@ module macroxs_header
real(8), allocatable :: polar(:) ! polar angles
real(8), allocatable :: azimuthal(:) ! azimuthal angles
! Type-Bound procedures
contains
procedure, pass :: init => macroxs_angle_init ! inits object
procedure, pass :: clear => macroxs_angle_clear ! Deallocates object
procedure, pass :: get_xs => macroxs_angle_get_xs ! Returns xs
! Type-Bound procedures
contains
procedure, pass :: init => macroxs_angle_init ! inits object
procedure, pass :: get_xs => macroxs_angle_get_xs ! Returns xs
end type MacroXS_Angle
!===============================================================================
@ -657,68 +654,6 @@ contains
end subroutine macroxs_angle_init
!===============================================================================
! MACROXS*_CLEAR resets and deallocates data in MacroXS.
!===============================================================================
subroutine macroxs_iso_clear(this)
class(MacroXS_Iso), intent(inout) :: this ! The MacroXS to clear
if (allocated(this % total)) then
deallocate(this % total, this % absorption, &
this % nu_fission)
end if
if (allocated(this % fission)) then
deallocate(this % fission)
end if
if (allocated(this % k_fission)) then
deallocate(this % k_fission)
end if
call this % scatter % clear()
if (allocated(this % chi)) then
deallocate(this % chi)
end if
end subroutine macroxs_iso_clear
subroutine macroxs_angle_clear(this)
class(MacroXS_Angle), intent(inout) :: this ! The MacroXS to clear
integer :: i, j
if (allocated(this % total)) then
deallocate(this % total, this % absorption, &
this % nu_fission)
end if
if (allocated(this % fission)) then
deallocate(this % fission)
end if
if (allocated(this % k_fission)) then
deallocate(this % k_fission)
end if
do i = 1, size(this % scatter,dim=2)
do j = 1, size(this % scatter,dim=1)
call this % scatter(j,i) % obj % clear()
end do
end do
if (allocated(this % scatter)) then
deallocate(this % scatter)
end if
if (allocated(this % chi)) then
deallocate(this % chi)
end if
end subroutine macroxs_angle_clear
!===============================================================================
! MACROXS_*_GET_XS returns the requested data type
!===============================================================================
@ -790,130 +725,4 @@ contains
end function macroxs_angle_get_xs
!===============================================================================
! THIN_GRID thins an (x,y) set while also thinning an associated y2
!===============================================================================
subroutine thin_grid(xout, yout, yout2, tol, compression, maxerr)
real(8), allocatable, intent(inout) :: xout(:) ! Resultant x grid
real(8), allocatable, intent(inout) :: yout(:) ! Resultant y values
real(8), allocatable, intent(inout) :: yout2(:) ! Secondary y values
real(8), intent(in) :: tol ! Desired fractional error to maintain
real(8), intent(out) :: compression ! Data reduction fraction
real(8), intent(inout) :: maxerr ! Maximum error due to compression
real(8), allocatable :: xin(:) ! Incoming x grid
real(8), allocatable :: yin(:) ! Incoming y values
real(8), allocatable :: yin2(:) ! Secondary Incoming y values
integer :: k, klo, khi
integer :: all_ok
real(8) :: x1, y1, x2, y2, x, y, testval
integer :: num_keep, remove_it
real(8) :: initial_size
real(8) :: error
real(8) :: x_frac
initial_size = real(size(xout), 8)
allocate(xin(size(xout)))
xin = xout
allocate(yin(size(yout)))
yin = yout
allocate(yin2(size(yout2)))
yin2 = yout2
all_ok = size(yin)
maxerr = 0.0_8
! This loop will step through each entry in dim==3 and check to see if
! all of the values in other 2 dims can be replaced with linear interp.
! If not, the value will be saved to a new array, if so, it will be
! skipped.
xout = 0.0_8
yout = 0.0_8
! Keep first point's data
xout(1) = xin(1)
yout(1) = yin(1)
yout2(1) = yin2(1)
! Initialize data
num_keep = 1
klo = 1
khi = 3
k = 2
do while (khi <= size(xin))
remove_it = 0
x1 = xin(klo)
x2 = xin(khi)
x = xin(k)
x_frac = 1.0_8 / (x2 - x1) * (x - x1) ! Linear interp.
! Check for removal. Otherwise, it stays. This is accomplished by leaving
! remove_it as 0, entering the else portion of if(remove_it==all_ok)
y1 = yin(klo)
y2 = yin(khi)
y = yin(k)
testval = y1 + (y2 - y1) * x_frac
error = abs(testval - y)
if (y /= 0.0_8) then
error = error / y
end if
if (error <= tol) then
remove_it = remove_it + 1
if (error > maxerr) then
maxerr = abs(testval - y)
end if
end if
! Now place the point in to the proper bin and advance iterators.
if (remove_it /= 0) then
! Then don't put it in the new grid but advance iterators
k = k + 1
khi = khi + 1
else
! Put it in new grid and advance iterators accordingly
num_keep = num_keep + 1
xout(num_keep) = xin(k)
yout(num_keep) = yin(k)
yout2(num_keep) = yin2(k)
klo = k
k = k + 1
khi = khi + 1
end if
end do
! Save the last point's data
num_keep = num_keep + 1
xout(num_keep) = xin(size(xin))
yout(num_keep) = yin(size(xin))
yout2(num_keep) = yin2(size(xin))
! Finally, xout and yout were sized to match xin and yin since we knew
! they would be no larger than those. Now we must resize these arrays
! and copy only the useful data in. Will use xin/yin for temp arrays.
xin = xout(1:num_keep)
yin = yout(1:num_keep)
yin2 = yout2(1:num_keep)
deallocate(xout)
deallocate(yout)
deallocate(yout2)
allocate(xout(num_keep))
allocate(yout(size(yin)))
allocate(yout2(size(yin2)))
xout = xin(1:num_keep)
yout = yin(1:num_keep)
yout2 = yin2(1:num_keep)
! Clean up
deallocate(xin)
deallocate(yin)
deallocate(yin2)
compression = (initial_size - real(size(xout),8)) / initial_size
end subroutine thin_grid
end module macroxs_header

2
src/math.F90
View file

@ -558,7 +558,7 @@ contains
end function calc_rn
!===============================================================================
! EXPAND_HARMONIC expands a given series of harmonics
! EXPAND_HARMONIC expands a given series of real spherical harmonics
!===============================================================================
pure function expand_harmonic(data, order, uvw) result(val)
real(8), intent(in) :: data(:)

3
src/mgxs_data.F90
View file

@ -76,8 +76,7 @@ contains
get_fiss = .true.
end if
end do
if (get_kfiss .and. get_fiss) &
exit
if (get_kfiss .and. get_fiss) exit
end do
! ==========================================================================

10
src/nuclide_header.F90
View file

@ -155,7 +155,7 @@ module nuclide_header
end function nuclide_calc_f_
end interface
!===============================================================================
!===============================================================================
! NUCLIDE_ISO contains the base MGXS data for a nuclide specifically for
! isotropically weighted MGXS
!===============================================================================
@ -688,11 +688,12 @@ module nuclide_header
end function nuclide_angle_get_xs
!===============================================================================
! NUCLIDE_*_CALC_F Finds the value of f(mu), the scattering probability, given mu
! NUCLIDE_*_CALC_F Finds the value of f(mu), the scattering angle probability,
! given mu
!===============================================================================
pure function nuclide_mg_iso_calc_f(this, gin, gout, mu, uvw, i_azi, i_pol) &
result(f)
pure function nuclide_mg_iso_calc_f(this, gin, gout, mu, uvw, i_azi, i_pol) &
result(f)
class(Nuclide_Iso), intent(in) :: this
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
@ -791,6 +792,7 @@ module nuclide_header
end if
end function nuclide_mg_angle_calc_f
!===============================================================================
! find_angle finds the closest angle on the data grid and returns that index
!===============================================================================

1
src/particle_header.F90
View file

@ -132,7 +132,6 @@ contains
this % fission = .false.
this % delayed_group = 0
this % n_delayed_bank(:) = 0
! Initialize this % g so there is always at least some initialized value
this % g = 1
! Set up base level coordinates

87
src/scattdata_header.F90
View file

@ -16,11 +16,10 @@ module scattdata_header
real(8), allocatable :: mult(:,:) ! (Gout x Gin)
real(8), allocatable :: data(:,:,:) ! (Order/Nmu x Gout x Gin)
! Type-Bound procedures
contains
procedure(init_), deferred, pass :: init ! Initializes ScattData
procedure(calc_f_), deferred, pass :: calc_f ! Calculates f, given mu
procedure(clear_), deferred, pass :: clear ! Deallocates ScattData
! Type-Bound procedures
contains
procedure(init_), deferred, pass :: init ! Initializes ScattData
procedure(calc_f_), deferred, pass :: calc_f ! Calculates f, given mu
end type ScattData_Base
abstract interface
@ -42,37 +41,29 @@ module scattdata_header
real(8) :: f ! Return value of f(mu)
end function calc_f_
subroutine clear_(this)
import ScattData_Base
class(ScattData_Base), intent(inout) :: this ! The ScattData to clear
end subroutine clear_
end interface
type, extends(ScattData_Base) :: ScattData_Legendre
contains
procedure, pass :: init => scattdata_legendre_init
procedure, pass :: calc_f => scattdata_legendre_calc_f
procedure, pass :: clear => scattdata_legendre_clear
contains
procedure, pass :: init => scattdata_legendre_init
procedure, pass :: calc_f => scattdata_legendre_calc_f
end type ScattData_Legendre
type, extends(ScattData_Base) :: ScattData_Histogram
real(8), allocatable :: mu(:) ! Mu bins
real(8) :: dmu ! Mu spacing
contains
procedure, pass :: init => scattdata_histogram_init
procedure, pass :: calc_f => scattdata_histogram_calc_f
procedure, pass :: clear => scattdata_histogram_clear
contains
procedure, pass :: init => scattdata_histogram_init
procedure, pass :: calc_f => scattdata_histogram_calc_f
end type ScattData_Histogram
type, extends(ScattData_Base) :: ScattData_Tabular
real(8), allocatable :: mu(:) ! Mu bins
real(8) :: dmu ! Mu spacing
real(8), allocatable :: fmu(:,:,:) ! PDF of f(mu)
contains
procedure, pass :: init => scattdata_tabular_init
procedure, pass :: calc_f => scattdata_tabular_calc_f
procedure, pass :: clear => scattdata_tabular_clear
contains
procedure, pass :: init => scattdata_tabular_init
procedure, pass :: calc_f => scattdata_tabular_calc_f
end type ScattData_Tabular
!===============================================================================
@ -239,56 +230,6 @@ contains
end subroutine scattdata_tabular_init
!===============================================================================
! SCATTDATA_CLEAR resets and deallocates data in ScattData.
!===============================================================================
subroutine scattdata_base_clear(this)
class(ScattData_Base), intent(inout) :: this
if (allocated(this % energy)) then
deallocate(this % energy)
end if
if (allocated(this % mult)) then
deallocate(this % mult)
end if
if (allocated(this % data)) then
deallocate(this % data)
end if
end subroutine scattdata_base_clear
subroutine scattdata_legendre_clear(this)
class(ScattData_Legendre), intent(inout) :: this
call scattdata_base_clear(this)
end subroutine scattdata_legendre_clear
subroutine scattdata_histogram_clear(this)
class(ScattData_Histogram), intent(inout) :: this
call scattdata_base_clear(this)
if (allocated(this % mu)) then
deallocate(this % mu)
end if
end subroutine scattdata_histogram_clear
subroutine scattdata_tabular_clear(this)
class(ScattData_Tabular), intent(inout) :: this
call scattdata_base_clear(this)
if (allocated(this % mu)) then
deallocate(this % mu)
end if
end subroutine scattdata_tabular_clear
!===============================================================================
! SCATTDATA_*_CALC_F Calculates the value of f given mu (and gin,gout pair)
!===============================================================================
@ -350,6 +291,4 @@ contains
end function scattdata_tabular_calc_f
end module scattdata_header

148
src/tally.F90
View file

@ -26,6 +26,8 @@ module tally
integer :: position(N_FILTER_TYPES - 3) = 0 ! Tally map positioning array
!$omp threadprivate(position)
procedure(score_general_intfc), pointer :: score_general => null()
procedure(get_scoring_bins_intfc), pointer :: get_scoring_bins => null()
@ -52,8 +54,6 @@ module tally
end interface
!$omp threadprivate(position)
contains
!===============================================================================
@ -1255,97 +1255,95 @@ contains
real(8) :: uvw(3)
select case(score_bin)
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
! Find the scattering order for a singly requested moment, and
! store its moment contribution.
if (t % moment_order(i) == 1) then
score = score * p % mu ! avoid function call overhead
else
score = score * calc_pn(t % moment_order(i), p % mu)
endif
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
! Find the scattering order for a singly requested moment, and
! store its moment contribution.
if (t % moment_order(i) == 1) then
score = score * p % mu ! avoid function call overhead
else
score = score * calc_pn(t % moment_order(i), p % mu)
endif
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN)
score_index = score_index - 1
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + num_nm
! Update number of total n,m bins for this n (m = [-n: n])
num_nm = 2 * n + 1
case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN)
score_index = score_index - 1
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + num_nm
! Update number of total n,m bins for this n (m = [-n: n])
num_nm = 2 * n + 1
! multiply score by the angular flux moments and store
! multiply score by the angular flux moments and store
!$omp critical (score_general_scatt_yn)
t % results(score_index: score_index + num_nm - 1, filter_index) &
% value = t &
% results(score_index: score_index + num_nm - 1, filter_index)&
% value &
+ score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw)
t % results(score_index: score_index + num_nm - 1, filter_index) &
% value = t &
% results(score_index: score_index + num_nm - 1, filter_index)&
% value &
+ score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw)
!$omp end critical (score_general_scatt_yn)
end do
i = i + (t % moment_order(i) + 1)**2 - 1
end do
i = i + (t % moment_order(i) + 1)**2 - 1
case(SCORE_FLUX_YN, SCORE_TOTAL_YN)
score_index = score_index - 1
num_nm = 1
if (t % estimator == ESTIMATOR_ANALOG .or. &
t % estimator == ESTIMATOR_COLLISION) then
uvw = p % last_uvw
else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
uvw = p % coord(1) % uvw
end if
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + num_nm
! Update number of total n,m bins for this n (m = [-n: n])
num_nm = 2 * n + 1
case(SCORE_FLUX_YN, SCORE_TOTAL_YN)
score_index = score_index - 1
num_nm = 1
if (t % estimator == ESTIMATOR_ANALOG .or. &
t % estimator == ESTIMATOR_COLLISION) then
uvw = p % last_uvw
else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
uvw = p % coord(1) % uvw
end if
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + num_nm
! Update number of total n,m bins for this n (m = [-n: n])
num_nm = 2 * n + 1
! multiply score by the angular flux moments and store
! multiply score by the angular flux moments and store
!$omp critical (score_general_flux_tot_yn)
t % results(score_index: score_index + num_nm - 1, filter_index) &
% value = t &
% results(score_index: score_index + num_nm - 1, filter_index)&
% value &
+ score * calc_rn(n, uvw)
t % results(score_index: score_index + num_nm - 1, filter_index) &
% value = t &
% results(score_index: score_index + num_nm - 1, filter_index)&
% value &
+ score * calc_rn(n, uvw)
!$omp end critical (score_general_flux_tot_yn)
end do
i = i + (t % moment_order(i) + 1)**2 - 1
end do
i = i + (t % moment_order(i) + 1)**2 - 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
! Find the scattering order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
! Find the scattering order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + 1
! get the score and tally it
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value &
+ score * calc_pn(n, p % mu)
end do
i = i + t % moment_order(i)
case default
! get the score and tally it
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
t % results(score_index, filter_index) % value &
+ score * calc_pn(n, p % mu)
end do
i = i + t % moment_order(i)
end select
case default
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
end select
end subroutine expand_and_score