mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Revised the mgxs_combine procedure to work with the sparse format of scattering data
This commit is contained in:
Adam Nelson
parent
0cc5f06114
commit
3b6f836ec4
2 changed files with 543 additions and 313 deletions
|
|
@ -380,18 +380,16 @@ module mgxs_header
|
|||
integer, intent(in) :: legendre_to_tabular_points ! Number of points to use
|
||||
! in that conversion
|
||||
|
||||
character(MAX_LINE_LEN) :: temp_str
|
||||
integer(HID_T) :: xsdata_grp, scatt_grp
|
||||
real(8), allocatable :: temp_arr(:), temp_2d(:, :)
|
||||
real(8), allocatable :: scatt_coeffs(:, :, :)
|
||||
real(8), allocatable :: temp_scatt(:, :, :)
|
||||
real(8) :: dmu, mu, norm
|
||||
integer :: order, order_dim, gin, gout, l, imu, length
|
||||
type(VectorInt) :: temps_to_read
|
||||
integer :: t
|
||||
type(Jagged2D) :: input_scatt(:)
|
||||
type(Jagged1D) :: temp_mult(:)
|
||||
integer, allocatable :: gmin(:), gmax(:)
|
||||
character(MAX_LINE_LEN) :: temp_str
|
||||
integer(HID_T) :: xsdata_grp, scatt_grp
|
||||
real(8), allocatable :: temp_arr(:), temp_2d(:, :)
|
||||
real(8) :: dmu, mu, norm
|
||||
integer :: order, order_dim, gin, gout, l, imu, length
|
||||
type(VectorInt) :: temps_to_read
|
||||
integer :: t
|
||||
type(Jagged2D), allocatable :: input_scatt(:), scatt_coeffs(:)
|
||||
type(Jagged1D), allocatable :: temp_mult(:)
|
||||
integer, allocatable :: gmin(:), gmax(:)
|
||||
|
||||
! Call generic data gathering routine (will populate the metadata)
|
||||
call mgxs_from_hdf5(this, xs_id, temperature, method, tolerance, &
|
||||
|
|
@ -480,28 +478,20 @@ module mgxs_header
|
|||
call fatal_error("Must provide absorption!")
|
||||
end if
|
||||
|
||||
! Get multiplication data if present
|
||||
allocate(temp_mult(groups, groups))
|
||||
if (check_dataset(xsdata_grp, "multiplicity matrix")) then
|
||||
call read_dataset(temp_mult, xsdata_grp, "multiplicity matrix")
|
||||
else
|
||||
temp_mult(:, :) = ONE
|
||||
end if
|
||||
|
||||
! Get scattering data
|
||||
if (.not. check_group(xsdata_grp, "scatter data")) &
|
||||
call fatal_error("Must provide 'scatter data'")
|
||||
scatt_grp = open_group(xsdata_grp, 'scatter data')
|
||||
! First get the outgoing group boundary indices
|
||||
if (check_dataset(xsdata_grp, "g_min")) then
|
||||
allocate(g_min(groups))
|
||||
call read_dataset(g_min, scatt_grp, "g_min")
|
||||
allocate(gmin(groups))
|
||||
call read_dataset(gmin, scatt_grp, "g_min")
|
||||
else
|
||||
call fatal_error("'g_min' for the scatter matrix must be provided")
|
||||
end if
|
||||
if (check_dataset(xsdata_grp, "g_max")) then
|
||||
allocate(g_max(groups))
|
||||
call read_dataset(g_max, scatt_grp, "g_max")
|
||||
allocate(gmax(groups))
|
||||
call read_dataset(gmax, scatt_grp, "g_max")
|
||||
else
|
||||
call fatal_error("'g_max' for the scatter matrix must be provided")
|
||||
end if
|
||||
|
|
@ -510,11 +500,11 @@ module mgxs_header
|
|||
! to hold the flattened data
|
||||
length = 0
|
||||
do gin = 1, groups
|
||||
length = length + order_dim * (g_max(gin) - gmin(gin) + 1)
|
||||
length = length + order_dim * (gmax(gin) - gmin(gin) + 1)
|
||||
end do
|
||||
! Allocate flattened array
|
||||
allocate(temp_arr(length))
|
||||
jf (.not. check_dataset(scatt_grp, 'scatter matrix') &
|
||||
if (.not. check_dataset(scatt_grp, 'scatter matrix')) &
|
||||
call fatal_error("'scatter matrix' must be provided")
|
||||
call read_dataset(temp_arr, scatt_grp, "scatter matrix")
|
||||
|
||||
|
|
@ -542,7 +532,6 @@ module mgxs_header
|
|||
end if
|
||||
|
||||
allocate(scatt_coeffs(gin))
|
||||
|
||||
if (this % scatter_type == ANGLE_LEGENDRE .and. &
|
||||
legendre_to_tabular) then
|
||||
this % scatter_type = ANGLE_TABULAR
|
||||
|
|
@ -576,7 +565,7 @@ module mgxs_header
|
|||
! Now that we have the integral, lets ensure that the distribution
|
||||
! is normalized such that it preserves the original scattering xs
|
||||
if (norm > ZERO) then
|
||||
scatt_coeffs(gin) % data(i:, gout) = &
|
||||
scatt_coeffs(gin) % data(:, gout) = &
|
||||
scatt_coeffs(gin) % data(:, gout) * &
|
||||
input_scatt(gin) % data(1, gout)
|
||||
end if
|
||||
|
|
@ -597,11 +586,11 @@ module mgxs_header
|
|||
! to hold the flattened data
|
||||
length = 0
|
||||
do gin = 1, groups
|
||||
length = length + (g_max(gin) - gmin(gin) + 1)
|
||||
length = length + (gmax(gin) - gmin(gin) + 1)
|
||||
end do
|
||||
! Allocate flattened array
|
||||
allocate(temp_arr(length))
|
||||
jf (.not. check_dataset(scatt_grp, 'multiplicity matrix') &
|
||||
if (.not. check_dataset(scatt_grp, 'multiplicity matrix')) &
|
||||
call fatal_error("'multiplicity matrix' must be provided")
|
||||
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
|
||||
|
||||
|
|
@ -674,17 +663,16 @@ module mgxs_header
|
|||
integer, intent(in) :: legendre_to_tabular_points ! Number of points to use
|
||||
! in that conversion
|
||||
|
||||
character(MAX_LINE_LEN) :: temp_str
|
||||
integer(HID_T) :: xsdata_grp
|
||||
real(8), allocatable :: temp_arr(:), temp_4d(:, :, :, :)
|
||||
real(8), allocatable :: temp_mult(:, :, :, :)
|
||||
real(8), allocatable :: scatt_coeffs(:, :, :, :, :)
|
||||
real(8), allocatable :: input_scatt(:, :, :, :, :)
|
||||
real(8), allocatable :: temp_scatt(:, :, :, :, :)
|
||||
real(8) :: dmu, mu, norm
|
||||
integer :: order, order_dim, gin, gout, l, imu, ipol, iazi
|
||||
type(VectorInt) :: temps_to_read
|
||||
integer :: t
|
||||
character(MAX_LINE_LEN) :: temp_str
|
||||
integer(HID_T) :: xsdata_grp, scatt_grp
|
||||
real(8), allocatable :: temp_arr(:), temp_4d(:, :, :, :)
|
||||
real(8) :: dmu, mu, norm
|
||||
integer :: order, order_dim, gin, gout, l, imu
|
||||
type(VectorInt) :: temps_to_read
|
||||
integer :: t, length, ipol, iazi
|
||||
type(Jagged2D), allocatable :: input_scatt(:, :, :), scatt_coeffs(:, :, :)
|
||||
type(Jagged1D), allocatable :: temp_mult(:, :, :)
|
||||
integer, allocatable :: gmin(:, :, :), gmax(:, :, :)
|
||||
|
||||
! Call generic data gathering routine (will populate the metadata)
|
||||
call mgxs_from_hdf5(this, xs_id, temperature, method, tolerance, &
|
||||
|
|
@ -807,115 +795,167 @@ module mgxs_header
|
|||
call fatal_error("Must provide absorption!")
|
||||
end if
|
||||
|
||||
! Get multiplication data if present
|
||||
allocate(temp_mult(groups,groups, this % n_azi, this % n_pol))
|
||||
if (check_dataset(xsdata_grp, "multiplicity matrix")) then
|
||||
call read_dataset(temp_mult, xsdata_grp, "multiplicity matrix")
|
||||
else
|
||||
temp_mult(:, :, :, :) = ONE
|
||||
end if
|
||||
|
||||
! Get scattering data
|
||||
! The input is gathered in the more user-friendly facing format of
|
||||
! Gout x Gin x Order x Azi x Pol. We will get it in that format in
|
||||
! input_scatt, but then need to convert it to a more useful ordering
|
||||
! for processing (Order x Gout x Gin x Azi x Pol).
|
||||
allocate(input_scatt(groups, groups, order_dim, this % n_azi, &
|
||||
this % n_pol))
|
||||
if (check_dataset(xsdata_grp, "scatter matrix")) then
|
||||
allocate(temp_arr(groups * groups * order_dim * this % n_azi * &
|
||||
this % n_pol))
|
||||
call read_dataset(temp_arr, xsdata_grp, "scatter matrix")
|
||||
input_scatt(:, :, :, :, :) = reshape(temp_arr, (/groups, groups, &
|
||||
order_dim, this % n_azi, this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
|
||||
! Compare the number of orders given with the maximum order of the
|
||||
! problem. Strip off the supefluous orders if needed.
|
||||
if (this % scatter_type == ANGLE_LEGENDRE) then
|
||||
order = min(order_dim - 1, max_order)
|
||||
order_dim = order + 1
|
||||
end if
|
||||
|
||||
allocate(temp_scatt(groups, groups, order_dim, this % n_azi, &
|
||||
this % n_pol))
|
||||
temp_scatt(:, :, :, :, :) = input_scatt(:, :, 1:order_dim, :, :)
|
||||
|
||||
! Take input format (groups, groups, order) and convert to
|
||||
! the more useful format needed for scattdata: (order, groups, groups)
|
||||
! However, if scatt_type was ANGLE_LEGENDRE (i.e., the data was
|
||||
! provided as Legendre coefficients), and the user requested that
|
||||
! these legendres be converted to tabular form (note xs is also
|
||||
! the default behavior), convert that now.
|
||||
if (this % scatter_type == ANGLE_LEGENDRE .and. legendre_to_tabular) then
|
||||
|
||||
! Convert input parameters to what we need for the rest.
|
||||
this % scatter_type = ANGLE_TABULAR
|
||||
order_dim = legendre_to_tabular_points
|
||||
order = order_dim
|
||||
dmu = TWO / real(order - 1, 8)
|
||||
|
||||
allocate(scatt_coeffs(order_dim, groups, groups, this % n_azi, &
|
||||
this % n_pol))
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
norm = ZERO
|
||||
do imu = 1, order_dim
|
||||
if (imu == 1) then
|
||||
mu = -ONE
|
||||
else if (imu == order_dim) then
|
||||
mu = ONE
|
||||
else
|
||||
mu = -ONE + real(imu - 1, 8) * dmu
|
||||
end if
|
||||
scatt_coeffs(imu, gout, gin, iazi, ipol) = &
|
||||
evaluate_legendre(temp_scatt(gout, gin, :, iazi, ipol), mu)
|
||||
! Ensure positivity of distribution
|
||||
if (scatt_coeffs(imu, gout, gin, iazi, ipol) < ZERO) &
|
||||
scatt_coeffs(imu, gout, gin, iazi, ipol) = ZERO
|
||||
! And accrue the integral
|
||||
if (imu > 1) then
|
||||
norm = norm + HALF * dmu * &
|
||||
(scatt_coeffs(imu - 1, gout, gin, iazi, ipol) + &
|
||||
scatt_coeffs(imu, gout, gin, iazi, ipol))
|
||||
end if
|
||||
end do
|
||||
! Now that we have the integral, lets ensure that the distribution
|
||||
! is normalized such that it preserves the original scattering xs
|
||||
if (norm > ZERO) then
|
||||
scatt_coeffs(:, gout, gin, iazi, ipol) = &
|
||||
scatt_coeffs(:, gout, gin, iazi, ipol) * &
|
||||
temp_scatt(gout, gin, 1, iazi, ipol) / norm
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
else
|
||||
! Sticking with current representation, carry forward but change
|
||||
! the array ordering
|
||||
allocate(scatt_coeffs(order_dim, groups, groups, this % n_azi, &
|
||||
this % n_pol))
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
do l = 1, order_dim
|
||||
scatt_coeffs(l, gout, gin, iazi, ipol) = &
|
||||
temp_scatt(gout, gin, l, iazi, ipol)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
deallocate(temp_scatt)
|
||||
if (.not. check_group(xsdata_grp, "scatter data")) &
|
||||
call fatal_error("Must provide 'scatter data'")
|
||||
scatt_grp = open_group(xsdata_grp, 'scatter data')
|
||||
! First get the outgoing group boundary indices
|
||||
if (check_dataset(xsdata_grp, "g_min")) then
|
||||
allocate(gmin(groups, this % n_azi, this % n_pol))
|
||||
call read_dataset(gmin, scatt_grp, "g_min")
|
||||
else
|
||||
call fatal_error("Must provide scatter matrix!")
|
||||
call fatal_error("'g_min' for the scatter matrix must be provided")
|
||||
end if
|
||||
if (check_dataset(xsdata_grp, "g_max")) then
|
||||
allocate(gmax(groups, this % n_azi, this % n_pol))
|
||||
call read_dataset(gmax, scatt_grp, "g_max")
|
||||
else
|
||||
call fatal_error("'g_max' for the scatter matrix must be provided")
|
||||
end if
|
||||
|
||||
! Now use this information to find the length of a container array
|
||||
! to hold the flattened data
|
||||
length = 0
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
length = length + order_dim * (gmax(gin, iazi, ipol) - &
|
||||
gmin(gin, iazi, ipol) + 1)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
! Allocate flattened array
|
||||
allocate(temp_arr(length))
|
||||
if (.not. check_dataset(scatt_grp, 'scatter matrix')) &
|
||||
call fatal_error("'scatter matrix' must be provided")
|
||||
call read_dataset(temp_arr, scatt_grp, "scatter matrix")
|
||||
|
||||
! Convert temp_arr to a jagged array ((gin) % data(l, gout)) for passing
|
||||
! to ScattData
|
||||
allocate(input_scatt(groups, this % n_azi, this % n_pol))
|
||||
index = 1
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
allocate(input_scatt(gin, iazi, ipol) % data(order_dim, &
|
||||
gmin(gin, iazi, ipol):gmax(gin, iazi, ipol)))
|
||||
do l = 1, order_dim
|
||||
do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol)
|
||||
input_scatt(gin, iazi, ipol) % data(l, gout) = &
|
||||
temp_arr(index)
|
||||
index = index + 1
|
||||
end do ! gout
|
||||
end do ! order
|
||||
end do ! gin
|
||||
end do ! iazi
|
||||
end do ! ipol
|
||||
deallocate(temp_arr)
|
||||
|
||||
! Finally convert the legendre to tabular if needed
|
||||
! Compare the number of orders given with the maximum order of the
|
||||
! problem. Strip off the supefluous orders if needed.
|
||||
if (this % scatter_type == ANGLE_LEGENDRE) then
|
||||
order = min(order_dim - 1, max_order)
|
||||
order_dim = order + 1
|
||||
end if
|
||||
|
||||
allocate(scatt_coeffs(gin, this % n_azi, this % n_pol))
|
||||
if (this % scatter_type == ANGLE_LEGENDRE .and. &
|
||||
legendre_to_tabular) then
|
||||
this % scatter_type = ANGLE_TABULAR
|
||||
order_dim = legendre_to_tabular_points
|
||||
order = order_dim
|
||||
dmu = TWO / real(order - 1, 8)
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
allocate(scatt_coeffs(gin, iazi, ipol) % data(&
|
||||
order_dim, groups))
|
||||
do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol)
|
||||
norm = ZERO
|
||||
do imu = 1, order_dim
|
||||
if (imu == 1) then
|
||||
mu = -ONE
|
||||
else if (imu == order_dim) then
|
||||
mu = ONE
|
||||
else
|
||||
mu = -ONE + real(imu - 1, 8) * dmu
|
||||
end if
|
||||
scatt_coeffs(gin, iazi, ipol) % data(imu, gout) = &
|
||||
evaluate_legendre(&
|
||||
input_scatt(gin, iazi, ipol) % data(:, gout), mu)
|
||||
! Ensure positivity of distribution
|
||||
if (scatt_coeffs(gin, iazi, ipol) % data(imu, gout) < ZERO) &
|
||||
scatt_coeffs(gin, iazi, ipol) % data(imu, gout) = ZERO
|
||||
! And accrue the integral
|
||||
if (imu > 1) then
|
||||
norm = norm + HALF * dmu * &
|
||||
(scatt_coeffs(gin, iazi, ipol) % data(imu - 1, gout) + &
|
||||
scatt_coeffs(gin, iazi, ipol) % data(imu, gout))
|
||||
end if
|
||||
end do ! mu
|
||||
! Now that we have the integral, lets ensure that the distribution
|
||||
! is normalized such that it preserves the original scattering xs
|
||||
if (norm > ZERO) then
|
||||
scatt_coeffs(gin, iazi, ipol) % data(:, gout) = &
|
||||
scatt_coeffs(gin, iazi, ipol) % data(:, gout) * &
|
||||
input_scatt(gin, iazi, ipol) % data(1, gout)
|
||||
end if
|
||||
end do ! gout
|
||||
end do ! gin
|
||||
end do ! iazi
|
||||
end do ! ipol
|
||||
else
|
||||
! Sticking with current representation, carry forward but change
|
||||
! the array ordering
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
allocate(scatt_coeffs(gin, iazi, ipol) % data(order_dim, groups))
|
||||
scatt_coeffs(gin, iazi, ipol) % data(:, :) = &
|
||||
input_scatt(gin, iazi, ipol) % data(:, :)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
deallocate(input_scatt)
|
||||
|
||||
! Now get the multiplication matrix
|
||||
! Now use this information to find the length of a container array
|
||||
! to hold the flattened data
|
||||
length = 0
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
length = length + (gmax(gin, iazi, ipol) - gmin(gin, iazi, ipol) + 1)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
! Allocate flattened array
|
||||
allocate(temp_arr(length))
|
||||
if (.not. check_dataset(scatt_grp, 'multiplicity matrix')) &
|
||||
call fatal_error("'multiplicity matrix' must be provided")
|
||||
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
|
||||
|
||||
! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing
|
||||
! to ScattData
|
||||
allocate(temp_mult(groups, this % n_azi, this % n_pol))
|
||||
index = 1
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
allocate(temp_mult(gin, iazi, ipol) % data( &
|
||||
gmin(gin, iazi, ipol):gmax(gin, iazi, ipol)))
|
||||
do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol)
|
||||
temp_mult(gin, iazi, ipol) % data(gout) = temp_arr(index)
|
||||
index = index + 1
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
deallocate(temp_arr)
|
||||
|
||||
! Allocate and initialize our ScattData Object.
|
||||
allocate(xs % scatter(this % n_azi, this % n_pol))
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
|
|
@ -929,9 +969,22 @@ module mgxs_header
|
|||
end if
|
||||
|
||||
! Initialize the ScattData Object
|
||||
call xs % scatter(iazi, ipol) % obj % init(&
|
||||
temp_mult(:, :, iazi, ipol), &
|
||||
scatt_coeffs(:, :, :, iazi, ipol))
|
||||
call xs % scatter(iazi, ipol) % obj % init(gmin(:, iazi, ipol), &
|
||||
gmax(:, iazi, ipol), temp_mult(:, iazi, ipol), &
|
||||
scatt_coeffs(:, iazi, ipol))
|
||||
end do
|
||||
end do
|
||||
|
||||
! Check sigA to ensure it is not 0 since it is
|
||||
! often divided by in the tally routines
|
||||
! (This may happen with Helium data)
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
if (xs % absorption(gin, iazi, ipol) == ZERO) then
|
||||
xs % absorption(gin, iazi, ipol) = 1E-10_8
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
|
@ -947,6 +1000,18 @@ module mgxs_header
|
|||
end do
|
||||
end if
|
||||
|
||||
! Finally, check sigT to ensure it is not 0 since it is
|
||||
! often divided by in the tally routines
|
||||
do ipol = 1, this % n_pol
|
||||
do iazi = 1, this % n_azi
|
||||
do gin = 1, groups
|
||||
if (xs % total(gin, iazi, ipol) == ZERO) then
|
||||
xs % total(gin, iazi, ipol) = 1E-10_8
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Close the groups we have opened and deallocate
|
||||
call close_group(xsdata_grp)
|
||||
deallocate(input_scatt, scatt_coeffs, temp_mult)
|
||||
|
|
@ -1087,6 +1152,10 @@ module mgxs_header
|
|||
integer :: nuc_t
|
||||
real(8) :: temp_actual, temp_desired
|
||||
integer :: scatter_type
|
||||
type(Jagged2D), allocatable :: nuc_matrix(:)
|
||||
integer, allocatable :: gmin(:), gmax(:)
|
||||
type(Jagged2D), allocatable :: jagged_scatt(:)
|
||||
type(Jagged1D), allocatable :: jagged_mult(:)
|
||||
|
||||
! Set the meta-data
|
||||
call mgxs_combine(this, temps, mat, nuclides, max_order, scatter_type, &
|
||||
|
|
@ -1170,7 +1239,13 @@ module mgxs_header
|
|||
end if
|
||||
end if
|
||||
|
||||
! Get the multiplication matrix
|
||||
! We will next gather the multiplicaity and scattering matrices.
|
||||
! To avoid multiple re-allocations as we resize the storage
|
||||
! matrix (and/or to avoidlots of duplicate code), we will use a
|
||||
! dense matrix for this storage, with a reduction to the sparse
|
||||
! format at the end.
|
||||
|
||||
! Get the multiplicity matrix
|
||||
! To combine from nuclidic data we need to use the final relationship
|
||||
! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
|
||||
! sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
|
||||
|
|
@ -1199,10 +1274,18 @@ module mgxs_header
|
|||
|
||||
! Get the complete scattering matrix
|
||||
nuc_order_dim = size(nuc % xs(nuc_t) % scatter % dist(1) % data, dim=1)
|
||||
scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :) = &
|
||||
scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :) + &
|
||||
atom_density * &
|
||||
nuc % xs(nuc_t) % scatter % get_matrix(min(nuc_order_dim, order_dim))
|
||||
nuc_order_dim = min(nuc_order_dim, order_dim)
|
||||
call nuc % xs(nuc_t) % scatter % get_matrix(nuc_order_dim, &
|
||||
nuc_matrix)
|
||||
do gin = 1, groups
|
||||
do gout = nuc % xs(nuc_t) % scatter % gmin(gin), &
|
||||
nuc % xs(nuc_t) % scatter % gmax(gin)
|
||||
scatt_coeffs(1:nuc_order_dim, gout, gin) = &
|
||||
scatt_coeffs(1: nuc_order_dim, gout, gin) + &
|
||||
atom_density * nuc_matrix(gin) % data(1:nuc_order_dim, gout)
|
||||
end do
|
||||
end do
|
||||
|
||||
type is (MgxsAngle)
|
||||
call fatal_error("Invalid passing of MgxsAngle to MgxsIso object")
|
||||
end select
|
||||
|
|
@ -1218,8 +1301,13 @@ module mgxs_header
|
|||
end do
|
||||
end do
|
||||
|
||||
! Now create our jagged data from the dense data
|
||||
call jagged_from_dense_2D(scatt_coeffs, jagged_scatt)
|
||||
call jagged_from_dense_1D(temp_mult, jagged_mult, gmin, gmax)
|
||||
|
||||
! Initialize the ScattData Object
|
||||
call this % xs(t) % scatter % init_from_dense(temp_mult, scatt_coeffs)
|
||||
call this % xs(t) % scatter % init(gmin, gmax, jagged_mult, &
|
||||
jagged_scatt)
|
||||
|
||||
! Now normalize chi
|
||||
if (mat % fissionable) then
|
||||
|
|
@ -1232,7 +1320,8 @@ module mgxs_header
|
|||
end if
|
||||
|
||||
! Deallocate temporaries
|
||||
deallocate(scatt_coeffs, temp_mult, mult_num, mult_denom)
|
||||
deallocate(jagged_mult, jagged_scatt, gmin, gmax, scatt_coeffs, &
|
||||
temp_mult, mult_num, mult_denom)
|
||||
end associate ! nuc
|
||||
end do NUC_LOOP
|
||||
end do TEMP_LOOP
|
||||
|
|
@ -1262,6 +1351,10 @@ module mgxs_header
|
|||
integer :: nuc_t
|
||||
real(8) :: temp_actual, temp_desired
|
||||
integer :: scatter_type
|
||||
type(Jagged2D), allocatable :: nuc_matrix(:)
|
||||
integer, allocatable :: gmin(:), gmax(:)
|
||||
type(Jagged2D), allocatable :: jagged_scatt(:)
|
||||
type(Jagged1D), allocatable :: jagged_mult(:)
|
||||
|
||||
! Set the meta-data
|
||||
call mgxs_combine(this, temps, mat, nuclides, max_order, scatter_type, &
|
||||
|
|
@ -1375,7 +1468,13 @@ module mgxs_header
|
|||
end if
|
||||
end if
|
||||
|
||||
! Get the multiplication matrix
|
||||
! We will next gather the multiplicaity and scattering matrices.
|
||||
! To avoid multiple re-allocations as we resize the storage
|
||||
! matrix (and/or to avoidlots of duplicate code), we will use a
|
||||
! dense matrix for this storage, with a reduction to the sparse
|
||||
! format at the end.
|
||||
|
||||
! Get the multiplicity matrix
|
||||
! To combine from nuclidic data we need to use the final relationship
|
||||
! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
|
||||
! sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
|
||||
|
|
@ -1411,20 +1510,31 @@ module mgxs_header
|
|||
|
||||
! Get the complete scattering matrix
|
||||
nuc_order_dim = &
|
||||
size(nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % dist(1) % data, dim=1)
|
||||
scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :, iazi, ipol) = &
|
||||
scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :, iazi, ipol) + &
|
||||
atom_density * &
|
||||
nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % get_matrix(min(nuc_order_dim, order_dim))
|
||||
end do
|
||||
end do
|
||||
size(nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % &
|
||||
dist(1) % data, dim=1)
|
||||
nuc_order_dim = min(nuc_order_dim, order_dim)
|
||||
call nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % &
|
||||
get_matrix(nuc_order_dim, nuc_matrix)
|
||||
do gin = 1, groups
|
||||
do gout = &
|
||||
nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % gmin(gin), &
|
||||
nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % gmax(gin)
|
||||
scatt_coeffs(1:nuc_order_dim, gout, gin, iazi, ipol) = &
|
||||
scatt_coeffs(1: nuc_order_dim, gout, gin, iazi, ipol) + &
|
||||
atom_density * nuc_matrix(gin) % data(1:nuc_order_dim, gout)
|
||||
end do ! gout
|
||||
end do ! gin
|
||||
end do ! iazi
|
||||
end do ! ipol
|
||||
type is (MgxsIso)
|
||||
call fatal_error("Invalid passing of MgxsIso to MgxsAngle object")
|
||||
end select
|
||||
|
||||
! Obtain temp_mult
|
||||
! Obtain temp_mult, create jaged arrays and initialize the
|
||||
! ScattData object.
|
||||
do ipol = 1, n_pol
|
||||
do iazi = 1, n_azi
|
||||
! Obtain temp_mult
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
if (mult_denom(gout, gin, iazi, ipol) > ZERO) then
|
||||
|
|
@ -1436,15 +1546,17 @@ module mgxs_header
|
|||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Initialize the ScattData Object
|
||||
do ipol = 1, n_pol
|
||||
do iazi = 1, n_azi
|
||||
call this % xs(t) % scatter(iazi, ipol) % obj % init_from_dense(&
|
||||
temp_mult(:, :, iazi, ipol), &
|
||||
scatt_coeffs(:, :, :, iazi, ipol))
|
||||
! Create the Jagged arrays
|
||||
call jagged_from_dense_2D(scatt_coeffs(:, :, :, iazi, ipol), &
|
||||
jagged_scatt)
|
||||
call jagged_from_dense_1D(temp_mult(:, :, iazi, ipol), &
|
||||
jagged_mult, gmin, gmax)
|
||||
|
||||
! Initialize the ScattData Object
|
||||
! Initialize the ScattData Object
|
||||
call this % xs(t) % scatter(iazi, ipol) % obj % init(gmin, &
|
||||
gmax, jagged_mult, jagged_scatt)
|
||||
end do
|
||||
end do
|
||||
|
||||
|
|
@ -1464,7 +1576,8 @@ module mgxs_header
|
|||
end if
|
||||
|
||||
! Deallocate temporaries
|
||||
deallocate(scatt_coeffs, temp_mult, mult_num, mult_denom)
|
||||
deallocate(jagged_mult, jagged_scatt, gmin, gmax, scatt_coeffs, &
|
||||
temp_mult, mult_num, mult_denom)
|
||||
end associate ! nuc
|
||||
end do NUC_LOOP
|
||||
end do TEMP_LOOP
|
||||
|
|
|
|||
|
|
@ -41,21 +41,21 @@ module scattdata_header
|
|||
real(8), allocatable :: scattxs(:) ! Isotropic Sigma_{s,g_{in}}
|
||||
|
||||
contains
|
||||
procedure(scattdata_init_), deferred :: init ! Initializes ScattData
|
||||
! Initializes ScattData from a dense matrix
|
||||
procedure(scattdata_init_dense_), deferred :: init_from_dense
|
||||
procedure(scattdata_init_), deferred :: init ! Initializes ScattData
|
||||
procedure(scattdata_calc_f_), deferred :: calc_f ! Calculates f, given mu
|
||||
procedure(scattdata_sample_), deferred :: sample ! sample the scatter event
|
||||
procedure :: get_matrix => scattdata_get_matrix ! Rebuild scattering matrix
|
||||
end type ScattData
|
||||
|
||||
abstract interface
|
||||
subroutine scattdata_init_dense_(this, mult, coeffs)
|
||||
import ScattData
|
||||
class(ScattData), intent(inout) :: this ! Object to work with
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
end subroutine scattdata_init_dense_
|
||||
subroutine scattdata_init_(this, gmin, gmax, mult, coeffs)
|
||||
import ScattData, Jagged1D, Jagged2D
|
||||
class(ScattData), intent(inout) :: this ! Object to work with
|
||||
integer, intent(in) :: gmin(:) ! Min Gout
|
||||
integer, intent(in) :: gmax(:) ! Max Gout
|
||||
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
|
||||
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
|
||||
end subroutine scattdata_init_
|
||||
|
||||
pure function scattdata_calc_f_(this, gin, gout, mu) result(f)
|
||||
import ScattData
|
||||
|
|
@ -81,7 +81,7 @@ module scattdata_header
|
|||
! Maximal value for rejection sampling from rectangle
|
||||
type(Jagged1D), allocatable :: max_val(:) ! (Gin % data(Gout))
|
||||
contains
|
||||
procedure :: init_from_dense => scattdatalegendre_init_from_dense
|
||||
procedure :: init => scattdatalegendre_init
|
||||
procedure :: calc_f => scattdatalegendre_calc_f
|
||||
procedure :: sample => scattdatalegendre_sample
|
||||
end type ScattDataLegendre
|
||||
|
|
@ -92,7 +92,7 @@ module scattdata_header
|
|||
! Histogram of f(mu) (dist has CDF)
|
||||
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
|
||||
contains
|
||||
procedure :: init_from_dense => scattdatahistogram_init_from_dense
|
||||
procedure :: init => scattdatahistogram_init
|
||||
procedure :: calc_f => scattdatahistogram_calc_f
|
||||
procedure :: sample => scattdatahistogram_sample
|
||||
procedure :: get_matrix => scattdatahistogram_get_matrix
|
||||
|
|
@ -104,7 +104,7 @@ module scattdata_header
|
|||
! PDF of f(mu) (dist has CDF)
|
||||
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
|
||||
contains
|
||||
procedure :: init_from_dense => scattdatatabular_init_from_dense
|
||||
procedure :: init => scattdatatabular_init
|
||||
procedure :: calc_f => scattdatatabular_calc_f
|
||||
procedure :: sample => scattdatatabular_sample
|
||||
procedure :: get_matrix => scattdatatabular_get_matrix
|
||||
|
|
@ -124,13 +124,15 @@ contains
|
|||
! SCATTDATA*_INIT builds the scattdata object
|
||||
!===============================================================================
|
||||
|
||||
subroutine scattdata_init_from_dense(this, order, energy, mult)
|
||||
class(ScattData), intent(inout) :: this ! Object to work on
|
||||
integer, intent(in) :: order ! Data Order
|
||||
real(8), intent(inout) :: energy(:, :) ! Energy Transfer Matrix
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
subroutine scattdata_init(this, order, gmin, gmax, energy, mult)
|
||||
class(ScattData), intent(inout) :: this ! Object to work on
|
||||
integer, intent(in) :: order ! Data Order
|
||||
integer, intent(in) :: gmin(:) ! Min Gout
|
||||
integer, intent(in) :: gmax(:) ! Max Gout
|
||||
type(Jagged1D), intent(inout) :: energy(:) ! Energy Transfer Matrix
|
||||
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
|
||||
|
||||
integer :: groups, gmin, gmax, gin
|
||||
integer :: groups, gin
|
||||
real(8) :: norm
|
||||
|
||||
groups = size(energy, dim=1)
|
||||
|
|
@ -140,83 +142,79 @@ contains
|
|||
allocate(this % energy(groups))
|
||||
allocate(this % mult(groups))
|
||||
allocate(this % dist(groups))
|
||||
! Use energy to find the gmin and gmax values
|
||||
! Also set energy values when doing it
|
||||
|
||||
this % gmin = gmin
|
||||
this % gmax = gmax
|
||||
|
||||
! Set the outgoing energy PDF values
|
||||
do gin = 1, groups
|
||||
! Make sure energy is normalized (i.e., CDF is 1)
|
||||
norm = sum(energy(:, gin))
|
||||
if (norm /= ZERO) energy(:, gin) = energy(:, gin) / norm
|
||||
! Find gmin by checking the P0 moment
|
||||
do gmin = 1, groups
|
||||
if (energy(gmin, gin) > ZERO) exit
|
||||
end do
|
||||
! Find gmax by checking the P0 moment
|
||||
do gmax = groups, 1, -1
|
||||
if (energy(gmax, gin) > ZERO) exit
|
||||
end do
|
||||
! Treat the case of all zeros
|
||||
if (gmin > gmax) then
|
||||
gmin = gin
|
||||
gmax = gin
|
||||
! By not changing energy(gin) here we are leaving it as zero
|
||||
end if
|
||||
allocate(this % energy(gin) % data(gmin:gmax))
|
||||
this % energy(gin) % data(gmin:gmax) = energy(gmin:gmax, gin)
|
||||
allocate(this % mult(gin) % data(gmin:gmax))
|
||||
this % mult(gin) % data(gmin:gmax) = mult(gmin:gmax, gin)
|
||||
allocate(this % dist(gin) % data(order, gmin:gmax))
|
||||
norm = sum(energy(gin) % data(:))
|
||||
if (norm /= ZERO) energy(gin) % data(:) = energy(gin) % data(:) / norm
|
||||
! Set the values
|
||||
allocate(this % energy(gin) % data(gmin(gin):gmax(gin)))
|
||||
this % energy(gin) % data(:) = energy(gin) % data(:)
|
||||
allocate(this % mult(gin) % data(gmin(gin):gmax(gin)))
|
||||
this % mult(gin) % data(:) = mult(gin) % data(:)
|
||||
allocate(this % dist(gin) % data(order, gmin(gin):gmax(gin)))
|
||||
this % dist(gin) % data = ZERO
|
||||
this % gmin(gin) = gmin
|
||||
this % gmax(gin) = gmax
|
||||
end do
|
||||
end subroutine scattdata_init_from_dense
|
||||
end subroutine scattdata_init
|
||||
|
||||
subroutine scattdatalegendre_init_from_dense(this, mult, coeffs)
|
||||
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
subroutine scattdatalegendre_init(this, gmin, gmax, mult, coeffs)
|
||||
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
|
||||
integer, intent(in) :: gmin(:) ! Min Gout
|
||||
integer, intent(in) :: gmax(:) ! Max Gout
|
||||
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
|
||||
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
|
||||
|
||||
real(8) :: dmu, mu, f, norm
|
||||
integer :: imu, Nmu, gout, gin, groups, order
|
||||
real(8), allocatable :: energy(:, :)
|
||||
real(8), allocatable :: matrix(:, :, :)
|
||||
type(Jagged1D), allocatable :: energy(:)
|
||||
type(Jagged2D), allocatable :: matrix(:)
|
||||
|
||||
groups = size(coeffs, dim=3)
|
||||
order = size(coeffs, dim=1)
|
||||
groups = size(coeffs)
|
||||
order = size(coeffs(1) % data, dim=1)
|
||||
|
||||
! make a copy of coeffs that we can use to extract data and normalize
|
||||
allocate(matrix(order, groups, groups))
|
||||
matrix (:, :, :)= coeffs
|
||||
allocate(matrix(groups))
|
||||
do gin = 1, groups
|
||||
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
|
||||
matrix(gin) % data = coeffs(gin) % data
|
||||
end do
|
||||
|
||||
! Get scattxs value
|
||||
allocate(this % scattxs(groups))
|
||||
! Get this by summing the un-normalized P0 coefficient in matrix
|
||||
! over all outgoing groups
|
||||
this % scattxs(:) = sum(matrix(1, :, :), dim=1)
|
||||
do gin = 1, groups
|
||||
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
|
||||
end do
|
||||
|
||||
allocate(energy(groups, groups))
|
||||
energy(:, :) = ZERO
|
||||
allocate(energy(groups))
|
||||
! Build energy transfer probability matrix from data in matrix
|
||||
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
norm = matrix(1, gout, gin)
|
||||
energy(gout, gin) = norm
|
||||
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
|
||||
energy(gin) % data = ZERO
|
||||
do gout = gmin(gin), gmax(gin)
|
||||
norm = matrix(gin) % data(1, gout)
|
||||
energy(gin) % data(gout) = norm
|
||||
if (norm /= ZERO) then
|
||||
matrix(:, gout, gin) = matrix(:, gout, gin) / norm
|
||||
matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
call scattdata_init_from_dense(this, order, energy, mult)
|
||||
call scattdata_init(this, order, gmin, gmax, energy, mult)
|
||||
|
||||
allocate(this % max_val(groups))
|
||||
! Set dist values from matrix and initialize max_val
|
||||
do gin = 1, groups
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
this % dist(gin) % data(:, gout) = matrix(:, gout, gin)
|
||||
do gout = gmin(gin), gmax(gin)
|
||||
this % dist(gin) % data(:, gout) = matrix(gin) % data(:, gout)
|
||||
end do
|
||||
allocate(this % max_val(gin) % data(this % gmin(gin):this % gmax(gin)))
|
||||
allocate(this % max_val(gin) % data(gmin(gin):gmax(gin)))
|
||||
this % max_val(gin) % data(:) = ZERO
|
||||
end do
|
||||
|
||||
|
|
@ -225,7 +223,7 @@ contains
|
|||
Nmu = 1001
|
||||
dmu = TWO / real(Nmu - 1, 8)
|
||||
do gin = 1, groups
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
do gout = gmin(gin), gmax(gin)
|
||||
do imu = 1, Nmu
|
||||
! Update mu. Do first and last seperate to avoid float errors
|
||||
if (imu == 1) then
|
||||
|
|
@ -246,46 +244,52 @@ contains
|
|||
this % max_val(gin) % data(gout) * 1.1_8
|
||||
end do
|
||||
end do
|
||||
end subroutine scattdatalegendre_init_from_dense
|
||||
end subroutine scattdatalegendre_init
|
||||
|
||||
subroutine scattdatahistogram_init_from_dense(this, mult, coeffs)
|
||||
class(ScattDataHistogram), intent(inout) :: this ! Object to work on
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
subroutine scattdatahistogram_init(this, gmin, gmax, mult, coeffs)
|
||||
class(ScattDataHistogram), intent(inout) :: this ! Object to work on
|
||||
integer, intent(in) :: gmin(:) ! Min Gout
|
||||
integer, intent(in) :: gmax(:) ! Max Gout
|
||||
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
|
||||
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
|
||||
|
||||
integer :: imu, gin, gout, groups, order
|
||||
real(8) :: norm
|
||||
real(8), allocatable :: energy(:, :)
|
||||
real(8), allocatable :: matrix(:, :, :)
|
||||
type(Jagged1D), allocatable :: energy(:)
|
||||
type(Jagged2D), allocatable :: matrix(:)
|
||||
|
||||
groups = size(coeffs, dim=3)
|
||||
order = size(coeffs, dim=1)
|
||||
groups = size(coeffs)
|
||||
order = size(coeffs(1) % data, dim=1)
|
||||
|
||||
! make a copy of coeffs that we can use to extract data and normalize
|
||||
allocate(matrix(order, groups, groups))
|
||||
matrix(:, :, :) = coeffs
|
||||
allocate(matrix(groups))
|
||||
do gin = 1, groups
|
||||
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
|
||||
matrix(gin) % data = coeffs(gin) % data
|
||||
end do
|
||||
|
||||
! Get scattxs value
|
||||
allocate(this % scattxs(groups))
|
||||
! Get this by summing the un-normalized P0 coefficient in matrix
|
||||
! over all outgoing groups
|
||||
this % scattxs(:) = sum(sum(matrix(:, :, :), dim=1), dim=1)
|
||||
do gin = 1, groups
|
||||
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
|
||||
end do
|
||||
|
||||
allocate(energy(groups, groups))
|
||||
energy(:, :) = ZERO
|
||||
allocate(energy(groups))
|
||||
! Build energy transfer probability matrix from data in matrix
|
||||
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
norm = sum(matrix(:, gout, gin))
|
||||
energy(gout, gin) = norm
|
||||
norm = sum(matrix(gin) % data(:, gout))
|
||||
energy(gin) % data(gout) = norm
|
||||
if (norm /= ZERO) then
|
||||
matrix(:, gout, gin) = matrix(:, gout, gin) / norm
|
||||
matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
call scattdata_init_from_dense(this, order, energy, mult)
|
||||
call scattdata_init(this, order, gmin, gmax, energy, mult)
|
||||
|
||||
allocate(this % mu(order))
|
||||
this % dmu = TWO / real(order, 8)
|
||||
|
|
@ -298,17 +302,16 @@ contains
|
|||
! also saving the original histogram in fmu
|
||||
allocate(this % fmu(groups))
|
||||
do gin = 1, groups
|
||||
allocate(this % fmu(gin) % data(order, &
|
||||
this % gmin(gin):this % gmax(gin)))
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
|
||||
do gout = gmin(gin), gmax(gin)
|
||||
! Store the histogram
|
||||
this % fmu(gin) % data(:, gout) = matrix(:, gout, gin)
|
||||
this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
|
||||
! Integrate the histogram
|
||||
this % dist(gin) % data(1, gout) = &
|
||||
this % dmu * matrix(1, gout, gin)
|
||||
this % dmu * matrix(gin) % data(1, gout)
|
||||
do imu = 2, order
|
||||
this % dist(gin) % data(imu, gout) = &
|
||||
this % dmu * matrix(imu, gout, gin) + &
|
||||
this % dmu * matrix(gin) % data(imu, gout) + &
|
||||
this % dist(gin) % data(imu - 1, gout)
|
||||
end do
|
||||
|
||||
|
|
@ -323,24 +326,29 @@ contains
|
|||
end do
|
||||
end do
|
||||
|
||||
end subroutine scattdatahistogram_init_from_dense
|
||||
end subroutine scattdatahistogram_init
|
||||
|
||||
subroutine scattdatatabular_init_from_dense(this, mult, coeffs)
|
||||
class(ScattDataTabular), intent(inout) :: this ! Object to work on
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
subroutine scattdatatabular_init(this, gmin, gmax, mult, coeffs)
|
||||
class(ScattDataTabular), intent(inout) :: this ! Object to work on
|
||||
integer, intent(in) :: gmin(:) ! Min Gout
|
||||
integer, intent(in) :: gmax(:) ! Max Gout
|
||||
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
|
||||
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
|
||||
|
||||
integer :: imu, gin, gout, groups, order
|
||||
real(8) :: norm
|
||||
real(8), allocatable :: energy(:, :)
|
||||
real(8), allocatable :: matrix(:, :, :)
|
||||
type(Jagged1D), allocatable :: energy(:)
|
||||
type(Jagged2D), allocatable :: matrix(:)
|
||||
|
||||
groups = size(coeffs, dim=3)
|
||||
order = size(coeffs, dim=1)
|
||||
groups = size(coeffs)
|
||||
order = size(coeffs(1) % data, dim=1)
|
||||
|
||||
! make a copy of coeffs that we can use to extract data and normalize
|
||||
allocate(matrix(order, groups, groups))
|
||||
matrix(:, :, :) = coeffs
|
||||
allocate(matrix(groups))
|
||||
do gin = 1, groups
|
||||
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
|
||||
matrix(gin) % data = coeffs(gin) % data
|
||||
end do
|
||||
|
||||
! Build the angular distribution mu values
|
||||
allocate(this % mu(order))
|
||||
|
|
@ -358,39 +366,39 @@ contains
|
|||
! over all outgoing groups
|
||||
do gin = 1, groups
|
||||
norm = ZERO
|
||||
do gout = 1, groups
|
||||
do gout = gmin(gin), gmax(gin)
|
||||
do imu = 2, order
|
||||
norm = norm + HALF * this % dmu * (matrix(imu - 1, gout, gin) + &
|
||||
matrix(imu, gout, gin))
|
||||
norm = norm + HALF * this % dmu * &
|
||||
(matrix(gin) % data(imu - 1, gout) + &
|
||||
matrix(gin) % data(imu, gout))
|
||||
end do
|
||||
end do
|
||||
this % scattxs(gin) = norm
|
||||
end do
|
||||
|
||||
allocate(energy(groups, groups))
|
||||
energy(:, :) = ZERO
|
||||
allocate(energy(groups))
|
||||
! Build energy transfer probability matrix from data in matrix
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
do gout = gmin(gin), gmax(gin)
|
||||
norm = ZERO
|
||||
do imu = 2, order
|
||||
norm = norm + HALF * this % dmu * &
|
||||
(matrix(imu - 1, gout, gin) + matrix(imu, gout, gin))
|
||||
(matrix(gin) % data(imu - 1, gout) + &
|
||||
matrix(gin) % data(imu, gout))
|
||||
end do
|
||||
energy(gout, gin) = norm
|
||||
energy(gin) % data(gout) = norm
|
||||
end do
|
||||
end do
|
||||
call scattdata_init_from_dense(this, order, energy, mult)
|
||||
call scattdata_init(this, order, gmin, gmax, energy, mult)
|
||||
|
||||
! Calculate f(mu) and integrate it so we can avoid rejection sampling
|
||||
allocate(this % fmu(groups))
|
||||
do gin = 1, groups
|
||||
allocate(this % fmu(gin) % data(order, &
|
||||
this % gmin(gin):this % gmax(gin)))
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
|
||||
do gout = gmin(gin), gmax(gin)
|
||||
! Coeffs contain f(mu), put in f(mu) as that is where the
|
||||
! PDF lives
|
||||
this % fmu(gin) % data(:, gout) = matrix(:, gout, gin)
|
||||
this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
|
||||
|
||||
! Force positivity
|
||||
do imu = 1, order
|
||||
|
|
@ -417,7 +425,7 @@ contains
|
|||
end if
|
||||
end do
|
||||
end do
|
||||
end subroutine scattdatatabular_init_from_dense
|
||||
end subroutine scattdatatabular_init
|
||||
|
||||
!===============================================================================
|
||||
! SCATTDATA_*_CALC_F Calculates the value of f given mu (and gin,gout pair)
|
||||
|
|
@ -642,10 +650,10 @@ contains
|
|||
! using ScattData's information of fmu/dist, energy, and scattxs
|
||||
!===============================================================================
|
||||
|
||||
pure function scattdata_get_matrix(this, req_order) result(matrix)
|
||||
class(ScattData), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
|
||||
subroutine scattdata_get_matrix(this, req_order, matrix)
|
||||
class(ScattData), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
|
||||
|
||||
integer :: order, groups, gin, gout
|
||||
|
||||
|
|
@ -653,66 +661,175 @@ contains
|
|||
! Set gin and gout for getting the order
|
||||
order = min(req_order, size(this % dist(1) % data, dim=1))
|
||||
|
||||
allocate(matrix(order, groups, groups))
|
||||
if (allocated(matrix)) deallocate(matrix)
|
||||
allocate(matrix(groups))
|
||||
! Initialize to 0; this way the zero entries in the dense matrix dont
|
||||
! need to be explicitly set, requiring a significant increase in the
|
||||
! lines of code.
|
||||
matrix(:, :, :) = ZERO
|
||||
do gin = 1, groups
|
||||
allocate(matrix(gin) % data(order, groups))
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
matrix(:, gout, gin) = this % scattxs(gin) * &
|
||||
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
|
||||
this % energy(gin) % data(gout) * &
|
||||
this % dist(gin) % data(1:order, gout)
|
||||
end do
|
||||
end do
|
||||
end function scattdata_get_matrix
|
||||
end subroutine scattdata_get_matrix
|
||||
|
||||
pure function scattdatahistogram_get_matrix(this, req_order) result(matrix)
|
||||
class(ScattDataHistogram), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
|
||||
subroutine scattdatahistogram_get_matrix(this, req_order, matrix)
|
||||
class(ScattDataHistogram), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
|
||||
|
||||
integer :: order, groups, gin, gout
|
||||
|
||||
groups = size(this % energy)
|
||||
order = min(req_order, size(this % dist(1) % data, dim=1))
|
||||
|
||||
allocate(matrix(order, groups, groups))
|
||||
if (allocated(matrix)) deallocate(matrix)
|
||||
allocate(matrix(groups))
|
||||
! Initialize to 0; this way the zero entries in the dense matrix dont
|
||||
! need to be explicitly set, requiring a significant increase in the
|
||||
! lines of code.
|
||||
matrix(:, :, :) = ZERO
|
||||
do gin = 1, groups
|
||||
allocate(matrix(gin) % data(order, groups))
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
matrix(:, gout, gin) = this % scattxs(gin) * &
|
||||
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
|
||||
this % energy(gin) % data(gout) * &
|
||||
this % fmu(gin) % data(1:order, gout)
|
||||
end do
|
||||
end do
|
||||
end function scattdatahistogram_get_matrix
|
||||
end subroutine scattdatahistogram_get_matrix
|
||||
|
||||
pure function scattdatatabular_get_matrix(this, req_order) result(matrix)
|
||||
class(ScattDataTabular), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
|
||||
subroutine scattdatatabular_get_matrix(this, req_order, matrix)
|
||||
class(ScattDataTabular), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
|
||||
|
||||
integer :: order, groups, gin, gout
|
||||
|
||||
groups = size(this % energy)
|
||||
order = min(req_order, size(this % dist(1) % data, dim=1))
|
||||
|
||||
allocate(matrix(order, groups, groups))
|
||||
if (allocated(matrix)) deallocate(matrix)
|
||||
allocate(matrix(groups))
|
||||
! Initialize to 0; this way the zero entries in the dense matrix dont
|
||||
! need to be explicitly set, requiring a significant increase in the
|
||||
! lines of code.
|
||||
matrix(:, :, :) = ZERO
|
||||
do gin = 1, groups
|
||||
allocate(matrix(gin) % data(order, groups))
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
matrix(:, gout, gin) = this % scattxs(gin) * &
|
||||
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
|
||||
this % energy(gin) % data(gout) * &
|
||||
this % fmu(gin) % data(1:order, gout)
|
||||
end do
|
||||
end do
|
||||
end function scattdatatabular_get_matrix
|
||||
end subroutine scattdatatabular_get_matrix
|
||||
|
||||
!===============================================================================
|
||||
! JAGGED_FROM_DENSE_*D Creates a jagged array from a sparse dense matrix.
|
||||
! The user can supply a key which indicates the values to remove, but the
|
||||
! default is ZERO
|
||||
!===============================================================================
|
||||
|
||||
subroutine jagged_from_dense_1D(dense, jagged, lo_bounds, hi_bounds, key_)
|
||||
real(8), intent(in) :: dense(:, :)
|
||||
type(Jagged1D), allocatable, intent(inout) :: jagged(:)
|
||||
real(8), intent(in), optional :: key_
|
||||
integer, intent(inout), allocatable, optional :: lo_bounds(:)
|
||||
integer, intent(inout), allocatable, optional :: hi_bounds(:)
|
||||
|
||||
real(8) :: key
|
||||
integer :: i, jmin, jmax
|
||||
|
||||
if (present(key_)) then
|
||||
key = key_
|
||||
else
|
||||
key = ZERO
|
||||
end if
|
||||
|
||||
if (present(lo_bounds)) then
|
||||
if (allocated(lo_bounds)) deallocate(lo_bounds)
|
||||
allocate(lo_bounds(size(dense, dim=2)))
|
||||
end if
|
||||
if (present(hi_bounds)) then
|
||||
if (allocated(hi_bounds)) deallocate(hi_bounds)
|
||||
allocate(hi_bounds(size(dense, dim=2)))
|
||||
end if
|
||||
|
||||
allocate(jagged(size(dense, dim=2)))
|
||||
do i = 1, size(dense, dim=2)
|
||||
! Find the min and max j values
|
||||
do jmin = 1, size(dense, dim=1)
|
||||
if (dense(jmin, i) > key) exit
|
||||
end do
|
||||
do jmax = size(dense, dim=1), 1, -1
|
||||
if (dense(jmax, i) > key) exit
|
||||
end do
|
||||
! Treat the case of all values matching the key
|
||||
if (jmin > jmax) then
|
||||
jmin = i
|
||||
jmax = i
|
||||
end if
|
||||
|
||||
! Now store the jagged row
|
||||
allocate(jagged(i) % data(jmin:jmax))
|
||||
jagged(i) % data(jmin:jmax) = dense(jmin:jmax, i)
|
||||
|
||||
if (present(lo_bounds)) lo_bounds(i) = jmin
|
||||
if (present(hi_bounds)) hi_bounds(i) = jmax
|
||||
end do
|
||||
|
||||
end subroutine jagged_from_dense_1D
|
||||
|
||||
subroutine jagged_from_dense_2D(dense, jagged, lo_bounds, hi_bounds, key_)
|
||||
real(8), intent(in) :: dense(:, :, :)
|
||||
type(Jagged2D), allocatable, intent(inout) :: jagged(:)
|
||||
real(8), intent(in), optional :: key_
|
||||
integer, intent(inout), allocatable, optional :: lo_bounds(:)
|
||||
integer, intent(inout), allocatable, optional :: hi_bounds(:)
|
||||
|
||||
real(8) :: key
|
||||
integer :: i, jmin, jmax
|
||||
|
||||
if (present(key_)) then
|
||||
key = key_
|
||||
else
|
||||
key = ZERO
|
||||
end if
|
||||
|
||||
if (present(lo_bounds)) then
|
||||
if (allocated(lo_bounds)) deallocate(lo_bounds)
|
||||
allocate(lo_bounds(size(dense, dim=3)))
|
||||
end if
|
||||
if (present(hi_bounds)) then
|
||||
if (allocated(hi_bounds)) deallocate(hi_bounds)
|
||||
allocate(hi_bounds(size(dense, dim=3)))
|
||||
end if
|
||||
|
||||
allocate(jagged(size(dense, dim=3)))
|
||||
do i = 1, size(dense, dim=3)
|
||||
! Find the min and max j values
|
||||
do jmin = 1, size(dense, dim=1)
|
||||
if (sum(dense(:, jmin, i)) > key) exit
|
||||
end do
|
||||
do jmax = size(dense, dim=1), 1, -1
|
||||
if (sum(dense(:, jmax, i)) > key) exit
|
||||
end do
|
||||
! Treat the case of all values matching the key
|
||||
if (jmin > jmax) then
|
||||
jmin = i
|
||||
jmax = i
|
||||
end if
|
||||
|
||||
! Now store the jagged row
|
||||
allocate(jagged(i) % data(size(dense, dim=1), jmin:jmax))
|
||||
jagged(i) % data(:, jmin:jmax) = dense(:, jmin:jmax, i)
|
||||
|
||||
if (present(lo_bounds)) lo_bounds(i) = jmin
|
||||
if (present(hi_bounds)) hi_bounds(i) = jmax
|
||||
end do
|
||||
|
||||
end subroutine jagged_from_dense_2D
|
||||
|
||||
end module scattdata_header
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue